program cbis
c*********************************************************************
c  Version 6-26-2009   M. J. Hinich
c  The data can be on two columns on one data file or on one column of one
c  data file and another on a column on a second data file.
c  File 3 is always opened. The output in 3 are statistics for the cross
c   bispectrum for the data used and summary statistics for this data
c  File 1 data file and if there are two, 2 - data file,
c  File 7 - xxy cross bispec file, File 4 - graph file
c======================================================================
c  Parameters in control file RUNCBIS.CNL. Use spaces between entries!
c  A delimiter symbol * must be placed on each line.
c  Place comments on the lines of runcbis.cnl file after the symbol *.
c  Runcbis.cnl allows runs a number of control files with different file
c  names containing different data files and parameters.
c  The parameters for each run are echoed in filename.out where filename
c  is the name of the data file used.
c===================================================================

c  Lines for the control file RUNCBIS.CNL
c  The parameters & commands are found using key words whose first letter
c  may be either upper or lower case. The rest of the key word must be in
c  lower case. The symbol ! is a wildcard for the value to be set.
c  Examples: !nrun below is the wildcard for the number of runs of the program
c  using different cnl files called by runcbis.cnl. The integer to the right
c  of the = sign in runcbis.cnl is read by the program and sets the number
c  of runs.
c  In any *.cnl called by runcbis.cnl, !r is the wildcard for the resolution
c  bandwidth. The real number to the right of !r= is the resolution bandwidth
c  that is used.
c
c  IMPORTANT: Leave at least ONE SPACE after the paramter
c
c  ALWAYS use 1.0 rather than .0 since the program searches for the number
c   to the left of the period in the parameter.
c
c  A delimiter symbol * must be placed on each line.
c
c  Use the keywords followed by an = sign and the setting. Some settings are
c   character strings while others are real numbers.
c  Comments on the lines of the cnl file are placed to the right of the
c   symbol *. These are not read by the program.
c
c===================================================================
c
c  Name & path of the fid.cnl control file=path/name *
c  Example:
c  name & path of the control file=\prog\cbis-electric.cnl  *
c
c============================
c
c  Number of data files=nfs!  *
c  Example:
c  number of data files=1 or 2  *
c
c=======================================================================
c
c  These entries are followed by the line
c  # * End of Runcbis.cnl file reads
c
c  The program seeks the delimiter symbol # to end the reading of the
c  control file names.
c
c============================
c
c  The run parameters are now read for each control file
c
c  The parameter lines are followed by the line
c  # */ End of run parameter reads from the control file
c
c  The program seaks the delimiter symbol # to end the reading of the
c  parameters for each control file set in Runcbis.cnl
c
c=========================
c
c  Data file name1=data file name with path  *
c
c  If there are two data files then enter below this line:
c
c  Data file name2=data file name with path  *
c
c=========================
c
c  Column name of x=idx  *
c  Column name of y=idy  *
c
c=========================
c
c  Column number for data x=!ncox   *
c  Column number for data y=!ncoy   *
c  !ncox & !ncoy are the column numbers of the nc columns in the data file
c  to be used for the cross bispectrum.
c
c=========================
c
c  Resolution bandwidth=!rb  *  ( !r - real wildcard)
c   This bandwidth is in Hz or the unit used below if frequency is not typed
c  Note: rb determines the frame length lb used to compute the spectrum
c
c=========================
c
c  Sampling unit=!unit *
c  kHz or Hz to interpret the sampling rate as a multiple of kHz or Hz
c  You can also use the lower case letters k or h.
c
c  If khz is used the spectrum output unit is kHz, or if hz is used the
c  spectrum output unit is Hz   OR
c
c  a 10 character unit in lower case to interpret the sampling rate
c  as a multipled of this unit. The output will be in periods using that
c  unit as for example  unit='sec'.
c
c=========================
c
c  Lower frequency=!rl  upper frequency=!ru or  *
c  Longest period=!rl   shortest period=!ru     *
c   The passband is (!rl,!ru)
c   !rl - lower frequency, !ru - upper frequency of band in units of kHz
c   if frequency+kHz  is typed or if OR Hz if frequency+Hz is typed on a line
c   !rl - longest period - the longest period
c   !ru - shortest period - the shortest period if a unit is used
c
c   If !rl=0.0 then the longest period will be the frame length, which is the
c   inverse of the resolution bandwidth
c   If !ru=0.0 then the normalized upper frequencey is 0.5 or the
c    normalized shortest period is lb/2
c
c=========================
c
c  File read start=!ns  file read end=!ne   *
c   !ns - index of datum first read, !ne - index of datum last read
c    If !ne = 0, all the data is read
c
c=========================
c
c  % taper=!rpt                 *
c   !rpt = % taper of frame for sidelobe reduction (0. < rpt < 25.)
c   If !rpt = 0, frames are NOT tapered
c
c=======================================================================
c
c  Normalization=divide/no      *
c   Type 'divide' after = if you want to divide the bispectrum at (f1,f2) by
c   sqrt[S(f1)S(f2)S(f1+f2)] where S(f) is the estimate of spectral
c   density at frequency f. This is a whitening operation.
c   If the character string 'divide' is not present or =no then the
c   bispectrum is then divided by the cube of the sample standard deviation.
c
c=======================================================================
c
c  Data file format for the first two lines (records)
c
c  Line 1: A character string that identifies the data in thie file .
c  Line 2: If the unit is in frequency then line 2 using an EXAMPLE format is
c  rows=!nfl, columns=!nc, sampling rate=!sr, format=(f10.5,2x,a12) unit=kHz
c  Line 2: If the unit is a time unit then line 2 is
c  rows=!nfl, columns=!nc, sampling interval=!sr, format=(f10.5,2x,a12) unit=day
c
c  The data file contains nc time series in column form where the data
c  format is set in the second line of the data file.
c  Up to 50 columns can be read.
c
c   nfl - no. of observations, sr - sampling rate in kHz if freq. is
c   set, or in multiples of the time unit set in cbis-fn.cnl
c
c   Data format is a character*50 string including ( )
c   unit=kHz or khz if that is the sr unit
c   unit=Hz or hz if that is the sr unit  OR unit=unit in time for sr
c   This character string is used to check the read from %.cnl
c
c   If there is an 'a' in the format then a time & data character string
c   is read. The time & data string must be less than 21 characters long.
c   An example format for a data file with two columns and a date & time
c   column 7 spaces wide is: (2(f7.2,1x),2x,a7).
c
c*****************************************************
      parameter(mt=12,ndy=30,ny=2008)
      real gs,rt
      integer run,nr
      character*700 name
      character*80 buf,cnl,id
      character*50 charead,fmt,fname
      character*24 dt
      character*20 par
      character*10 su
      character*2 delim
      logical exs
c**********************
c  Start time
      call cpu_time(gs)
      open(9,file='Cbis-error.out')
c  Inquire if file exists
      inquire(file='Runcbis.cnl',exist=exs)
      if(.not.exs) then
       write(6,'(/)')
       write(6,*)'Runcbis.cnl file does not exist'
       write(9,*)'Runcbis.cnl file does not exist'
       stop
      endif
c  Open summary output file
      open(3,file='Cbis.out',err=1,iostat=io1)
      write(3,'(11x,''XXY Crosbipectrum Statistics by M. J. Hinich'')')
      call clock(gs,0,3)
c  Open Runcbis.cnl top control file
      open(1,file='Runcbis.cnl',err=2,iostat=io2,status='old')
      ip=0
      kt=0
      ia=0
      dowhile(kt<=7.and.ia==0)
       kt=kt+1
c  Clear buf
       do n=1,80
        buf(n:n)=' '
       enddo
c  Read line
       read(1,'(a80)',end=3,err=4,iostat=io3) buf
       ia=index(buf,'#')
       ig=index(buf,'*')
       if(ig==0) then
        write(6,*)'Place comment marker * after each line in Runcbis.cnl
     *l'
        write(9,*)'Place comment marker * after each line in Runcbis.cnl
     *'
        write(6,'(a70)')buf
        write(9,'(a70)')buf
        stop
       endif
       ig=ig-1
       name(ip+1:ip+ig)=buf(:ig)
       ip=ip+ig
      enddo
c  Name & path of the control file
      delim(1:2)='= '
      par='rol file'
      cnl=charead(par,name,delim,ia,ier,3)
      if(ier>0) then
       write(3,*)'Error in reading the cnl file name in Runcbis.cnl'
       write(6,*)'Error in reading the cnl file name in Runcbis.cnl'
       write(9,*)'Error in reading the cnl file name in Runcbis.cnl'
       stop
      endif
      inquire(file=cnl,exist=exs)
      if(.not.exs) then
       write(6,'(/)')
       write(6,*)'cnl file does not exist ',cnl
       write(9,'(/)')
       write(9,*)'cnl file does not exist ',cnl
       stop
      endif
c  Open fid.cnl control file
      open(11,file=cnl,err=5,iostat=io5,status='old')
      par='ata files'
      nfs=numread(par,name,delim,ia,ier,3)
      if(ier>0) then
       write(3,*)'Error in reading the no. of data files in Runcbis.cnl'
       write(6,*)'Error in reading the no. of data files in Runcbis.cnl'
       write(9,*)'Error in reading the no. of data files in Runcbis.cnl'
       stop
      endif
c-------
      call top(nfs,mt,ndy,ny,3)
c+++++++++++++++++++++++++++++++++++++++++++
c  Time and date of run & run time
      call head(mt,ndy,ny,gs,rt,3)
c+++++++++++++++++++++++++++++++++++++++++++
      call flush(3)
      stop
  1   write(6,*)'**** Error on opening Cbis.out ',io1
      stop
  2   write(6,*)'**** Error on opening Runcbis.cnl ',io2
      stop
  3   write(6,*)'End of file on runcbis.cnl reads'
      write(9,*)'End of file on runcbis.cnl reads'
  4   write(6,*)'Place the symbol # after the cnl file read ',io3
      write(9,*)'Place the symbol # after the cnl file read ',io3
      stop
  5   write(6,*)'**** Error on opening fid.cnl ',io5
      write(6,*)cnl
      write(9,*)cnl
      stop
c-------
      end
c
c===========================================
c
      subroutine top(nfs,mt,ndy,ny,io)
c********************************************************
c  Open data file(s) and read data
c========================================================
      allocatable::a,a2,br,bi,c1,c2,f,g,ik,sx,sy,ts,v,w,xr,xi,xy
      real a(:,:),a2(:,:),br(:),bi(:),sx(:),sy(:),w(:),xr(:),xi(:),xy(:)
      real*8 v(:),rlb
      complex*16 f(:),g(:)
      integer c1(:),c2(:),ik(:),div,q,nfs
      character*700 name
      character*80 buf,id
      character*50 charead,fmt,fmt2,fname,idx,idy
      character*25 ts(:)
      character*24 dt
      character*20 par
      character*10 su,su2
      character*6 divide
      character*2 delim
      character mk
      logical exs
      intent(in) nfs,mt,ndy,ny,io
c********************************************
      delim(1:2)='= '
      ip=0
      do k=1,17
c  Clear buf
       do n=1,80
        buf(n:n)=' '
       enddo
c  Read line
       read(11,'(a80)',end=1,err=2,iostat=io1) buf
       ia=index(buf,'#')
       if(ia>0) then
        exit
       endif
       ig=index(buf,'*')
       if(ig==0) then
        write(io,*)'Place comment marker * after parameters on fid.cnl'
        write(io,'(a70)')buf
        write(6,*)'Place comment marker * after parameters on fid.cnl'
        write(6,'(a70)')buf
        write(9,*)'Place comment marker * after parameters on fid.cnl'
        write(9,'(a70)')buf
        stop
       endif
       ig=ig-1
       name(ip+1:ip+ig)=buf(:ig)
       ip=ip+ig
      enddo
c-------
      if(nfs==1) then
       call fopen(name,buf,nfs,ia,ib,io)
      endif
      if(nfs==2) then
       call fopen(name,buf,1,ia,ib,io)
       call fopen(name,buf,2,ia,ib,io)
      endif
c-------
      fname=buf(ib:ia)
      nfn=ia-ib+2
      fname(nfn-1:nfn+11)='-cbis.out'
c  Open cbis output list file
      open(7,file=fname,err=3,iostat=io3)
      call fdate(dt)
      write(7,'('' Version '',i2,''-'',i2,''-'',i4,7x,''Run '',a24)')mt
     *,ndy,ny,dt
      fname(nfn-1:nfn+12)='.graph'
c  Open cbis probability graph file
      open(4,file=fname,err=4,iostat=io4)
      call fdate(dt)
      write(4,'('' Version '',i2,''-'',i2,''-'',i4,7x,''Run '',a24)')mt
     *,ndy,ny,dt
      fname(nfn-1:nfn+11)='.gain'
c  Open cbis probability graph file
      open(8,file=fname,err=5,iostat=io5)
      call fdate(dt)
      write(8,'('' Version '',i2,''-'',i2,''-'',i4,7x,''Run '',a24)')mt
     *,ndy,ny,dt
c-------
c  x column name
      par='ame of x'
      idx=charead(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(io,'(4x,''ERROR reading x column name '',a50)')idx
       write(6,'(4x,''ERROR reading x column name '',a50)')idx
       write(9,'(4x,''ERROR reading x column name '',a50)')idx
       stop
      endif
c  y column name
      par='ame of y'
      idy=charead(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(io,'(4x,''ERROR reading y column name '',a50)')idy
       write(6,'(4x,''ERROR reading y column name '',a50)')idy
       write(9,'(4x,''ERROR reading y column name '',a50)')idy
       stop
      endif
c  x column number
      kcox=0
      par='ata x'
      kcox=numread(par,name,delim,ia,ier,io)
      if(ier>0.or.kcox<=0) then
       write(io,'(7x,''ERROR reading x column number '',i3)')kcox
       write(6,'(7x,''ERROR reading x column number '',i3)')kcox
       write(9,'(7x,''ERROR reading x column number '',i3)')kcox
       stop
      endif
c-------
c  y column number
      kcoy=0
      par='ata y'
      kcoy=numread(par,name,delim,ia,ier,io)
      if(ier>0.or.kcoy<=0) then
       write(io,'(7x,''ERROR reading y column number '',i3)')kcoy
       write(6,'(7x,''ERROR reading y column number '',i3)')kcoy
       write(9,'(7x,''ERROR reading y column number '',i3)')kcoy
       stop
      endif
c  Resolution bandwidth
      par='idth'
      rb=rnumread(par,name,delim,ia,ier,io)
      if(ier==4) then
       write(io,*)'Place a period in the resolution bandwidth value'
       ia=index(name,'idth')
       write(io,*) name(ia+5:ia+15)
       stop
      endif
      if(rb<=0.) then
       write(io,*)'Resolution bandwidth < 0'
       write(6,*)'Resolution bandwidth < 0'
       stop
      endif
c  Sampling unit
      par='unit'
      buf=charead(par,name,delim,ia,ier,3)
      if(ier==1) then
       write(3,'(/7x,''Sampling unit in fid.cnl file is incorrect'')')
       write(6,'(/7x,''Sampling unit in fid.cnl file is incorrect'')')
       write(9,'(/7x,''Sampling unit in fid.cnl file is incorrect'')')
       stop
      endif
      su=buf(:10)
      ia=index(buf,'kHz')
      ib=index(buf,'khz')
      ic=index(buf,'Hz')
      ij=index(buf,'hz')
      if(ia>0.or.ib>0) then
c  kHz
       iunit=1
      endif
      if((ic>0.or.ij>0).and.(ia==0.and.ib==0)) then
c  Hz
       iunit=2
      endif
      if(ia==0.and.ib==0.and.ic==0.and.ij==0) then
c  Unit set in cnl file
       iunit=0
      endif
c-------
      if(iunit>0) then
c  Lower frequency
       par='wer freq'
       fl=rnumread(par,name,delim,ia,ier,io)
       if(ier==1) then
        write(io,'(/'' Type lower frequency rather than longest period f
     *or the passband in %.cnl'')')
        write(6,'(/'' Type lower frequency rather than longest period fo
     *r the passband in %.cnl'')')
        stop
       endif
c  Upper frequency
       par='per freq'
       fu=rnumread(par,name,delim,ia,ier,io)
       if(ier==1) then
        write(io,'(/'' Type upper frequency rather than shortest period
     *for the passband in %.cnl'')')
        write(6,'(/'' Type upper frequency rather than shortest period f
     *or the passband in %.cnl'')')
        stop
       endif
        else
c  Longest period
         par='gest period'
         fl=rnumread(par,name,delim,ia,ier,io)
         if(ier==1) then
          write(io,'(/'' Type longest period rather than lower frequency
     * for the passband in %.cnl'')')
          write(6,'(/'' Type longest period rather than lower frequency
     *for the passband in %.cnl'')')
          stop
         endif
c  Shortest period
         par='test period'
         fu=rnumread(par,name,delim,ia,ier,io)
         if(ier==1) then
          write(io,'(/'' Type shortest period rather than upper frequenc
     *y for the passband in %.cnl'')')
          write(6,'(/'' Type shortest period rather than upper frequency
     * for the passband in %.cnl'')')
          stop
         endif
      endif
c-------
c  % taper
      par='taper'
      pt=rnumread(par,name,delim,ia,ier,io)
c-------
      if(iunit==0) then
c  Trap fl=0 for low band period
       if(fl<=0.0) then
        fl=1.e6
       endif
c  Invert period to frequency
       fl=1.0/fl
       if(fu>0.0) then
        fu=1.0/fu
         else
          fu=0.5
       endif
      endif
c-------
      if(fu>0.0.and. fl>fu) then
       write(io,*)'Upper frequency of band ',fu,' < lower frequency ',fl
       if(iunit==1) then
        write(io,*)'Check to see the (fl,fu) units are kHz'
         elseif(iunit==2) then
          write(io,*)'Check to see the (fl,fu) units are Hz'
           else
            write(io,*)'Check to see the (fl,fu) units are ',su
       endif
       stop
      endif
c-------
c  Trap taper %
      if(pt<0.or.pt>25.) then
       write(io,*)'Taper % not in range (0,25) ',pt
       stop
      endif
c  File read start
      par='start'
      is=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(io,*) 'Error ',ier,' in read of start datum'
       write(io,*) par
       stop
      endif
c  File read end
      par='end'
      ie=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(io,*)'Error in read of end datum'
       write(io,*)par
       write(6,*)'Error in read of end datum'
       write(6,*)par
       stop
      endif
      if(is<1) then
       write(io,*)'IS in CNL file is < 1 ',is
       write(6,*)'IS in CNL file is < 1 ',is
       write(9,*)'IS in CNL file is < 1 ',is
       stop
        elseif(ie/=0.and.ie<is) then
         write(io,*)'IE < IS in CNL file ',is,ie
         write(6,*)'IE < IS in CNL file ',is,ie
         write(9,*)'IE < IS in CNL file ',is,ie
          stop
      endif
c-------
      div=0
c  Normalization=divide/no
      par='zation'
      divide=charead(par,name,delim,ia,ier,3)
      idv=index(divide,'ivide')
      if(idv>0) then
       div=1
      endif
c-------
c  Read data files header
      if(nfs==1) then
       call readhead(id,nfl,nc,sr,fmt,iunit,nfs,io)
       write(3,'(1x,a79)')id
      endif
      if(nfs==2) then
       call readhead(id,nfl,nc,sr,fmt,iunit,1,io)
       write(3,'(1x,a79)')id
       call readhead(id,nfl2,nc2,sr2,fmt2,iunit,2,io)
       write(3,'(1x,a79)')id
      endif
      write(io,'(2x,''X Column Number = '',i3,2x,''Y Column Number = '',
     *i3,2x,''No. of Columns = '',i3/)')kcox,kcoy,nc
c  Trap parameters in data files
       if(nfl2/=nfl) then
        write(6,'(/2x,''File size in file 1 /= file size in file 2 '',i6
     *,1x,i6)')nfl,nfl2
        write(9,'(/2x,''File size in file 1 /= file size in file 2 '',i6
     *,1x,i6)')nfl,nfl2
        stop
       endif
       if(nc2/=nc) then
        write(6,'(/2x,''No. of columns in file 1 =/ no. in file 2 '',i3,
     *1x,i3)')nc,nc2
        write(9,'(/2x,''No. of columns in file 1 =/ no. in file 2 '',i3,
     *1x,i3)')nc,nc2
        stop
       endif
       if(sr2/=sr) then
        write(6,'(2x,''Sampling interval in file 1 /= interval in file 2
     *'')')
        write(6,'(2x,f10,4x,2x,f104)')sr,sr2
        write(9,'(2x,''Sampling interval in file 1 /= interval in file 2
     *'')')
        write(9,'(2x,f10,4x,2x,f104)')sr,sr2
        stop
       endif
       if(fmt2/=fmt) then
        write(6,'(/2x,''Format in file 1 /= format in file 2'')')
        write(6,'(1x,a50,1x,a50)')fmt,fmt2
        write(9,'(/2x,''Format in file 1 /= format in file 2'')')
        write(9,'(1x,a50,1x,a50)')fmt,fmt2
        stop
       endif
c-------
c  Parse format line in buf to see if there is a 'a' for time marks
      mk='a'
      ia=index(fmt,'a')
      if(ia>0) then
       mk='m'
      endif
      nmk=1
      if(mk=='m') then
       nmk=nfl
      endif
c-------
c  Trap top band
      if(fu>sr/2.0) then
       write(io,*)'Upper frequency = ',fu,' > Nyquist freq. ',sr/2.
       write(6,*)'Upper frequency = ',fu,' > Nyquist freq. ',sr/2.
       write(9,*)'Upper frequency = ',fu,' > Nyquist freq. ',sr/2.
       if(iunit==1) then
        write(io,*)'Check to see the (fl,fu) units are kHz'
        write(6,*)'Check to see the (fl,fu) units are kHz'
        write(9,*)'Check to see the (fl,fu) units are kHz'
         elseif(iunit==2) then
          write(io,*)'Check to see the (fl,fu) units are Hz'
          write(6,*)'Check to see the (fl,fu) units are Hz'
          write(9,*)'Check to see the (fl,fu) units are Hz'
           else
            write(io,*)'Check to see the (fl,fu) units are ',su
            write(6,*)'Check to see the (fl,fu) units are ',su
            write(9,*)'Check to see the (fl,fu) units are ',su
       endif
       write(io,*)' Sampling unit ',su,' Sampling rate/interval ',sr
       write(io,*)'iunit = ',iunit
       write(6,*)' Sampling unit ',su,' Sampling rate/interval ',sr
       write(6,*)'iunit = ',iunit
       write(9,*)' Sampling unit ',su,' Sampling rate/interval ',sr
       write(9,*)'iunit = ',iunit
       stop
      endif
c  Convert to normalized values (0,.5)
      if(iunit>0) then
       gl=fl/sr
       if(fu>0.0) then
        gu=fu/sr
         else
          gu=0.5
       endif
        else
         gl=fl*sr
         gu=fu*sr
      endif
c-------
c  Trap file length
      if(ie>nfl) then
       write(io,*) 'IE ',ie, '> N in data file',nfl
       stop
        elseif(ie==0) then
         ie=nfl
      endif
      n=ie-is+1
c  Find frame length lb
      if(iunit==1) then
c  kHz
       lb=nint(sr*1.e3/rb)
        elseif(iunit==2) then
c  Hz
         lb=nint(sr/rb)
          else
           lb=nint(rb/sr)
      endif
c-------
      if(lb<=9) then
       write(io,*) 'Frame length = ',lb,' < 10'
       write(io,*) 'Decrease resolution bandwidth = ',rb
       stop
c  Trap n<lb
        elseif(n/lb<4) then
         write(io,*) 'Frame length > n/4'
         write(io,*) 'Increase resolution bandwidth = ',rb
         stop
      endif
c-------
c  Number of values tapered from each end
      nt=nint(pt*float(lb)/200.)
c  If no. of tapered values at ends < 2 set pt=0
      if(nt<2) then
       pt=0.
      endif
      rlb=dfloat(lb)
      n2=lb/2
c  Set il - lower band integer limit
      il=nint(sngl(dble(gl)*rlb))
      il=max(il,1)
c  Set iu - upper band integer limit
      if(fu>0.) then
       iu=nint(sngl(dble(gu)*rlb))
       iu=min(iu,n2)
        else
         iu=n2
      endif
      if(il==iu) then
       write(6,*)'fu-fl is too small. Increase fu-fl'
       write(9,*)'fu-fl is too small. Increase fu-fl'
       stop
        elseif(il==iu .and. iu==n2) then
         write(6,*)'fl ',fl,'=',fu,' fu: ERROR in band settings in cnl'
         write(9,*)'fl ',fl,'=',fu,' fu: ERROR in band settings in cnl'
        stop
      endif
c-------
c  No. of frequencies used
      jf=iu-il+1
      e=alog(float(lb))/alog(float(n))
      if(e<=0.01.or.e>=.99) then
       write(io,*)'e = ',e,' is out of bound'
       write(6,*)'e = ',e,' is out of bound'
       write(9,*)'e = ',e,' is out of bound'
       stop
        elseif(e>=0.33) then
         write(io,'(/2x,''Condition for Consistency of Estimates is Viol
     *ated - e ='',f7.4)')e
         write(io,'(7x,''Resolution Bandwidth = '',f10.2)')rb
      endif
c-------
c  Get array length for cross bispec
      call csize(il,iu,nbis)
      if(nbis<3) then
       write(3,*) 'No. of cross bispec values < 3'
       stop
      endif
c-------
      nn=nbis-1
      na=4*lb+14
      q=lb-1
      n2=lb/2-1
c  Allocate space
      allocate(br(0:nn),bi(0:nn),c1(0:nn),c2(0:nn),f(0:q),g(0:q),ik(0:nn
     *),sx(0:n2),sy(0:n2),v(0:na),xr(0:nn),xi(0:nn),xy(0:q),w(0:q),stat=
     *ibad)
      if(ibad/=0) then
       write(io,*)'Unable to allocate work space for arrays in top'
       write(6,*)'Unable to allocate work space for arrays in top'
       write(9,*)'Unable to allocate work space for arrays in top'
       stop
      endif
      allocate(a(0:n-1,0:nc-1),a2(0:n-1,0:nc-1),stat=ibad)
      if(ibad/=0) then
       write(3,*)'Unable to allocate work space for data arrays'
       write(6,*)'Unable to allocate work space for data arrays'
       write(9,*)'Unable to allocate work space for data arrays'
       stop
      endif
      allocate(ts(0:nmk-1),stat=ibad)
      if(ibad/=0) then
       write(io,*)'Unable to allocate work space for ts'
       write(6,*)'Unable to allocate work space for ts'
       write(9,*)'Unable to allocate work space for ts'
       stop
      endif
c-------
      call datc(a,ts,fmt,mk,nmk,is,ie,n,nc,1,io)
      if(nfs==2) then
       call datc(a2,ts,fmt,mk,nmk,is,ie,n,nc,2,io)
      endif
c-------
c  Write header for summary stats in file 3
      write(io,'(70(''*''))')
      call wr(idx,idy,su,sr,e,n,lb,pt,is,ie,il,iu,iunit,div,3)
c  Write header for xxy cbis plot file 7
      write(7,'(70(''*''))')
      call wr(idx,idy,su,sr,e,n,lb,pt,is,ie,il,iu,iunit,div,7)
c-------
      if(nfs==1) then
       call cb(a(0,kcox-1),a(0,kcoy-1),br,bi,c1,c2,f,g,ik,sx,sy,v,w,xr,x
     *i,xy,idx,idy,n,lb,il,iu,jf,iunit,kcox,kcoy,nbis,sr,pt,e,su,rb,div,
     *io)
      endif
      if(nfs==2) then
       call cb(a(0,kcox-1),a2(0,kcoy-1),br,bi,c1,c2,f,g,ik,sx,sy,v,w,xr,
     *xi,xy,idx,idy,n,lb,il,iu,jf,iunit,kcox,kcoy,nbis,sr,pt,e,su,rb,div
     *,io)
      endif
c-------
      deallocate(a,a2,br,bi,c1,c2,f,g,ik,sx,sy,ts,v,w,xr,xi,xy)
      return
  1   write(6,*)'End of file on runcbis.cnl reads'
      write(9,*)'End of file on runcbis.cnl reads'
  2   write(6,*)'Place the symbol # after the cnl file read ',io1
      write(9,*)'Place the symbol # after the cnl file read ',io1
      stop
  3   write(6,*)'**** Error on opening ',fname,' ',io3
      write(9,*)'**** Error on opening ',fname,' ',io3
      stop
  4   write(6,*)'**** Error on opening graph file',io4
      write(9,*)'**** Error on opening graph file',io4
      stop
  5   write(6,*)'**** Error on opening gain file',io5
      write(9,*)'**** Error on opening gain file',io5
      stop
      end
c
c===========================================
c
      subroutine cb(x,y,br,bi,c1,c2,f,g,ik,sx,sy,v,w,xr,xi,xy,idx,idy,n,
     *lb,il,iu,jf,iunit,kcox,kcoy,nbis,sr,pt,e,su,rb,div,io)
c**********************************************************************
c  xxy cross bispectrum
c   x,y - data arrays, su - sampling unit, sr - sampling rate
c   il - bin no. of lower freq, iu - bin no. of upper freq of passband
c   lb - frame length, div=1 divide by spectra normalization
c   nbis - no. of cross bispec values, pt - % taper
c======================================================================
      allocatable::u
      real x(n),y(n),sx(lb/2),sy(lb/2),xr(nbis),xi(nbis),br(nbis),bi(nbi
     *s),w(lb),xy(lb),u(:,:)
      real*8 v(4*lb+15)
      complex*16 f(0:lb-1),g(0:lb-1)
      integer c1(nbis),c2(nbis),ik(nbis),div,t
      character*50 fmt,idx,idy
      character*10 su
      intent(in)idx,idy,n,lb,il,iu,jf,iunit,kcox,kcoy,nbis,sr,pt,e,su,rb
     *,div,io
c****************************
      allocate(u(0:iu-1,0:iu-1),stat=ibad)
      if(ibad/=0) then
       write(io,*)'Unable to allocate work space for u in cb'
       write(6,*)'Unable to allocate work space for u in cb'
       write(9,*)'Unable to allocate work space for u in cb'
       stop
      endif
c  Compute statistics for input data
      call stat(n,x,ax,sdx,skx,c4x,c6x,smaxx,sminx,io)
      call stat(n,y,ay,sdy,sky,c4y,c6y,smaxy,sminy,io)
      write(io,'(60(''=''))')
      write(io,'(10x,''Descriptive Statistics of X Data'')')
      call wrstat(io,ax,sdx,skx,c4x,c6x,smaxx,sminx)
      write(io,'(/60(''=''))')
      write(io,'(10x,''Descriptive Statistics of Y Data'')')
      call wrstat(io,ay,sdy,sky,c4y,c6y,smaxy,sminy)
c-------
c  Initialize v for lb
      call cffti(lb,v)
c-------
c  Compute spectra and cross bispecta loops
      i2=nbis+1
      call xxy(x,y,xr,xi,br,bi,c1,c2,w,sx,sy,xy,f,g,v,n,lb,il,iu,nbis,pt
     *,div,x1,x2,ax,ay,sdx,sdy,io)
c-------
c  Output xxy cross bispec list
      call wrtcis(c1,c2,xr,xi,br,nbis,lb,su,sr,iunit,7)
      write(io,'(70(''#'')/)')
      call outcis(c1,c2,xr,xi,br,ik,nbis,lb,su,sr,x1,x2,iunit,io)
c  Output xxy probabilities for graphing
      call wrtcisn(c1,c2,br,u,iu,nbis,rb,iunit,4)
c  Output xxy gains for graphing
      call wrtcisn(c1,c2,xr,u,iu,nbis,rb,iunit,8)
      write(io,'(70(''*'')/)')
      call wrtspxy(sx,sy,xy,lb,il,iu,su,sr,iunit,io)
c--------
      deallocate(u)
      return
      end
c
c===========================================
c
      subroutine xxy(x,y,xr,xi,br,bi,c1,c2,w,sx,sy,xy,f,g,v,n,lb,il,iu,n
     *bis,pt,div,x1,x2,ax,ay,sdx,sdy,io)
c********************************************
c  Input: il - low freq bin of passband, iu - high freq bin of passband
c   nbis - no. of cbis values, div - divide switch
c   ax, ay - means, sdx,sdy - standard deviations of x & y
c  Output: xr, xi - xxy crosbis, c1, c2 - freq values, br - prob values
c   x,y - data arrays, n - sample size, lb - frame size, pt - taper %
c   (br,bi) xxx bispectrum with symmetries
c   sx,sy - spectra of inputs, xy - cross spectrum
c============================================================
      real x(n),y(n),sx(lb/2),sy(lb/2),xr(nbis),xi(nbis),br(nbis),bi(nbi
     *s),w(lb),xy(lb)
      real*8 v(4*lb+15),dlb
      complex*16 f(0:lb-1),g(0:lb-1),cs2,cs3,cxy
      integer c1(nbis),c2(nbis),bp,div,t,tb,tp
      intent(in)x,y,v,n,lb,il,iu,nbis,div,ax,ay,sdx,sdy,pt,io
      intent(out)xr,xi,br,bi,c1,c2,sx,sy,xy,x1,x2
c********************************************
      n2=lb/2
      jf=iu-il+1
      np=n/lb
      sq2=sqrt(2.)
      bk=float(np)
c  Scale factor
      dlb=dfloat(lb)
      cs2=dcmplx(dlb)
      cs3=dcmplx(dsqrt(dlb)*dlb)
c--------
c  Initialize arrays
      do k=1,n2
       sx(k)=0.
       sy(k)=0.
       xy(k)=0.
       xy(k+n2)=0.
      enddo
      do k=1,nbis
       xr(k)=0.
       xi(k)=0.
       br(k)=0.
       bi(k)=0.
       c1(k)=0
       c2(k)=0
      enddo
c-------
c  Initialize taper
      if(pt>0.) then
       call taper(w,lb,pt,sw,nt,io)
c  If no. of tapered values at ends < 2 set pt=0
       if(nt<2) then
        sw=1.
       endif
        else
         sw=1.
      endif
c-------
c  Start loop over frames
      do tb=1,np
c  Time-1 of tb block
       tp=(tb-1)*lb
       do k=1,lb
        t=tp+k
        xt=(x(t)-ax)/sdx
        yt=(y(t)-ay)/sdy
        if(pt>0.) then
c  Taper frame kb
          xt=w(k)*xt
          yt=w(k)*yt
        endif
c  Move x,y tb frame to complex f,g
        f(k-1)=dcmplx(dble(xt))
        g(k-1)=dcmplx(dble(yt))
       enddo
c  FFT of frames
       call cfftf(lb,f,v)
       call cfftf(lb,g,v)
c  Initialize cross bispec counter
       ic=0
c  Triangle il<k1<iu, k2<k1
       do k1=il,iu
c  Periodograms for k1
        sx(k1)=sx(k1)+sngl(dreal(f(k1)/cs2*f(lb-k1)))
        sy(k1)=sy(k1)+sngl(dreal(g(k1)/cs2*g(lb-k1)))
c  Cross spectrum
        cxy=f(k1)/cs2*g(lb-k1)
        xy(k1)=xy(k1)+sngl(dreal(cxy))
        xy(n2+k1)=xy(n2+k1)+sngl(dimag(cxy))
c--------
        if(tb==np) then
c  Set indices if kb=np
         sx(k1)=sx(k1)/bk
         sy(k1)=sy(k1)/bk
         xy(k1)=xy(k1)/cmplx(bk)
         xy(n2+k1)=xy(n2+k1)/cmplx(bk)
         if(pt>0) then
          sx(k1)=sx(k1)/sw
          sy(k1)=sy(k1)/sw
          xy(k1)=xy(k1)/cmplx(sw)
          xy(n2+k1)=xy(n2+k1)/cmplx(sw)
         endif
         if(sx(k1)<1.e-6) then
          sx(k1)=1.e-6
         endif
         if(sy(k1)<1.e-6) then
          sy(k1)=1.e-6
         endif
        endif
c--------
c  Loop on k2 if k1>1
        if(k1>1) then
         do k2=il,k1-1
          ks=k1-k2
          if(ks>=il .and. ks<=iu) then
           ic=ic+1
c  xxy spectrum
           xr(ic)=xr(ic)+sngl(dreal(f(k1)/cs3*f(lb-ks)*g(lb-k2)))
           xi(ic)=xi(ic)+sngl(dimag(f(k1)/cs3*f(lb-ks)*g(lb-k2)))
c  xxx spectrum
           br(ic)=br(ic)+sngl(dreal(f(k1)/cs3*f(lb-ks)*f(lb-k2)))
           bi(ic)=bi(ic)+sngl(dimag(f(k1)/cs3*f(lb-ks)*f(lb-k2)))
           if(tb==np) then
c  Set indices
            c1(ic)=k1
            c2(ic)=k2
            xr(ic)=xr(ic)/bk
            xi(ic)=xi(ic)/bk
            br(ic)=br(ic)/bk
            bi(ic)=bi(ic)/bk
           endif
          endif
         enddo
        endif
c-------
c  Triangle il<k1<iu, k1<k2<min(2*k1,iu)
        do k2=k1+1,min(2*k1,iu)
         ks=k2-k1
         if(ks>=il .and. ks<=iu) then
          ic=ic+1
c  xxy spectrum
          xr(ic)=xr(ic)+sngl(dreal(f(k1)/cs3*f(ks)*g(lb-k2)))
          xi(ic)=xi(ic)+sngl(dimag(f(k1)/cs3*f(ks)*g(lb-k2)))
c  xxx spectrum
          br(ic)=br(ic)+sngl(dreal(f(k1)/cs3*f(ks)*f(lb-k2)))
          bi(ic)=bi(ic)+sngl(dimag(f(k1)/cs3*f(ks)*f(lb-k2)))
          if(tb==np) then
c  Set indices
           c1(ic)=k1
           c2(ic)=k2
           if(k2<2*k1) then
            xr(ic)=xr(ic)/bk
            xi(ic)=xi(ic)/bk
            br(ic)=br(ic)/bk
            bi(ic)=bi(ic)/bk
             else
              xr(ic)=xr(ic)/(sq2*bk)
              xi(ic)=xi(ic)/(sq2*bk)
              br(ic)=br(ic)/(sq2*bk)
              bi(ic)=bi(ic)/(sq2*bk)
           endif
          endif
         endif
        enddo
c  End of k2 loop
       enddo
c  End of k1 loop
      enddo
c  End of loop on frames
c--------
      call sn(xr,xi,br,bi,c1,c2,sx,sy,np,lb,nbis,x1,x2,sw,div,io)
c-------
      return
      end
c
c===========================================
c
      subroutine sn(xr,xi,br,bi,c1,c2,sx,sy,np,lb,nbis,x1,x2,sw,div,io)
c********************************************************
c  Computes periodogram & cross bispectrum after loop in xxy
c   (xr,xi) complex xxy cross bispectrum, c1 & c2 - freq values
c   (br,bi) xxx bispectrum with symmetries
c   sx, sy - spectra
c   x1 - sum of squares of real normalized cross bispectrum
c   x2 - sum of squares of imag normalized cross bispectrum
c========================================================
      real sx(lb/2),sy(lb/2),xr(nbis),xi(nbis),br(nbis),bi(nbis),cr,ci
      integer c1(nbis),c2(nbis),div
      intent(in)sx,sy,np,lb,nbis,sw,div,io
      intent(out)x1,x2
      intent(in out)xr,xi,br,bi
c********************************************
      bk2=float(np)*2.
c  Initialize count for spec=0 values
      n0=0
c  Accumulate squared real & imag xxy
      a1=0.
      a2=0.
c  Divide cross bispectrum by spectra or variances
      do k=1,nbis
       if(div==1) then
        h=sx(c1(k))*sx(abs(c2(k)-c1(k)))*sy(c2(k))
        hb=sx(c1(k))*sx(abs(c2(k)-c1(k)))*sx(c2(k))
         else
          h=sw*sw*sw
          hb=sw*sw*sw
       endif
c--------
       h=sqrt(h/bk2)
       hb=sqrt(hb/bk2)
       tt=.001/float(lb)
       if(h>tt .and. hb>tt) then
c  Denominator is not too small
        br(k)=br(k)/hb
        bi(k)=bi(k)/hb
        xr(k)=xr(k)/h
        xi(k)=xi(k)/h
c  Sum squares of real and imaginary parts of cross bispec
        ar=xr(k)*xr(k)
        ai=xi(k)*xi(k)
        a1=a1+ar
        a2=a2+ai
c  Divide by the bispectrum of x
        ab=br(k)*br(k)+bi(k)*bi(k)
        cr=(xr(k)*br(k)+xi(k)*bi(k))/ab
        ci=(xi(k)*br(k)-xr(k)*bi(k))/ab
c  Gain & phase of cross bispec
        ch2=cr*cr+ci*ci
        xr(k)=sqrt(ch2)
        xi(k)=atan2(ci,cr)*57.2958
c  Uniform variates
        br(k)=p2chi(ar+ai,0.,4)
         else
          xr(k)=0.
          xi(k)=0.
          br(k)=1.e-9
          bi(k)=1.e-9
          n0=n0+1
       endif
      enddo
c-------
c  No. of terms with non-zero spectral products
      nx=nbis-n0
      rnx=float(nx)
c  xxy crossbis stats
c  p-value of test statistic U(0,1) of normalized a1
      call cdchi(a1,rnx,x1,3)
      x1=1.-x1
c  p-value of test statistic U(0,1) of normalized a2
      call cdchi(a2,rnx,x2,3)
      x2=1.-x2
c-------
      return
      end
c
c===========================================
c
      subroutine wr(idx,idy,su,sr,e,n,lb,pt,is,ie,il,iu,iunit,div,io)
c********************************************************
c  Writes run information for main output in output files
c========================================================
      integer div
      character*50 idx,idy
      character*10 su
      intent(in) idx,idy,su,sr,e,n,lb,pt,is,ie,il,iu,iunit,div,io
c*****************************************
      jf=iu-il+1
      rl=float(lb)
      write(io,'(3x,''X Data Name'',2x,a35,1x,''Y Data Name'',2x,a35)')i
     *dx,idy
      write(io,'(7x,''First datum read'',i7,2x,''Last datum read'',i8)')
     *is,ie
c-------
      write(io,'(2x,57(''='')/)')
      if(iunit==1) then
       rb=1.e3*sr/rl
c  Lower and upper frequencies for band
       a=float(il)/rl*sr
       b=float(iu)/rl*sr
       write(io,'('' Sampling rate ='',f12.3,'' kHz.'',3x,''Frame length
     * '',i6/)')sr,lb
       write(io,'('' Resolution Bandwidth = '',f7.3,'' Hz.'',2x,''No. of
     * frames '',i6/)')rb,n/lb
       write(io,'('' Passband ('',f7.4,2x,f7.4,'' ) kHz'')')a,b
       if(il>1) then
c  Bandpass filter residuals
        write(io,'(/1x,''Data Filtered using an FFT to Remove the Low Fr
     *equency Band (0, '',f7.4,'') kHz'')')a
       endif
      endif
c-------
      if(iunit==2) then
       rb=1.e3*sr/rl
c  Lower and upper frequencies for band
       a=float(il)/rl*sr*1.e3
       b=float(iu)/rl*sr*1.e3
       write(io,'('' Sampling rate ='',f12.3,'' kHz.'',3x,''Frame length
     * '',i6/)')1.e3*sr,lb
       write(io,'('' Resolution Bandwidth = '',f7.3,'' Hz.'',2x,''No. of
     * frames '',i6/)')rb,n/lb
       write(io,'('' Passband ('',f7.4,2x,f7.4,'' ) Hz'')')a,b
      endif
c-------
      if(iunit==0) then
       rb=rl*sr
       a=rb/float(il)
       b=rb/float(iu)
       write(io,'('' Sampling interval = '',g12.6,1x,a7,3x,''Frame lengt
     *h '',i6/)')sr,su,lb
       write(io,'('' Resolution Bandwidth = '',f7.3,1x,a7,2x,''No. of fr
     *ames '',i6/)')rb,su,n/lb
       write(io,'('' Passband ('',f9.2,2x,f9.2,'' )'',1x,a7)')a,b,su
      endif
      write(io,'(/7x,''No. of frequencies in band '',i6/)')jf
c  Number of values tapered from each end
      nt=nint(pt*rl/200.)
c  If no. of tapered values at ends < 2 set nt=0
      if(nt<2) then
       nt=0
      endif
      if(pt>0.0.and.nt>0) then
       write(io,'('' % taper = '',f5.2,4x,''No. tapered at each end = ''
     *,i5/)')pt,nt
      endif
      if(nt==0) then
       write(io,'('' % Taper yields no. of tapered values < 2''/)')
      endif
      write(io,'(2x,4x,''exp = '',f5.2/)')e
c-------
      if(div==1) then
       write(io,'(7x,''Cross bispectrum is divided by the products of sp
     *ectral values''/)')
        else
         write(io,'(1x,''Cross Bispectrum is divided by the cube of the
     *estimated variance''/)')
      endif
c-------
      return
      end
c
c===========================================
c
      subroutine wrstat(i,am,sig,sk,c4,c6,smax,smin)
c****************************************************
c  Data statistics
c============================================
      write(i,'(60(''='')/)')
      write(i,'('' Mean ='',g12.3,7x,''Std Dev ='',g12.3/)') am,sig
      write(i,'('' Skew ='',g10.3,4x,''Kurtosis ='',g10.3,4x,''C(6) ='',
     *f10.3/)') sk,c4,c6
      write(i,'('' Max value ='',g12.3,7x,''Min value ='',g12.3)') smax,
     *smin
      return
      end
c
c=============================================
c
      subroutine wrtspxy(sx,sy,xy,lb,il,iu,su,sr,iunit,io)
c*******************************************
c  Output x & y spectra & the xy coherence function
c===========================================
      real sx(lb/2),sy(lb/2),xy(lb)
      character*10 su
      intent(in) sx,sy,xy,lb,il,iu,su,sr,iunit,io
c********************************************
      n2=lb/2
      if(iunit==1) then
       rb=sr/float(lb)
       write(io,'('' Frequency/kHz'',2x,''x Spectrum'',6x,''y Spectrum
     *'',3x,''Gain'',3x,''Phase'',1x,''Coherence''/)')
       do k=il,iu
        gn=(xy(k)*xy(k)+xy(n2+k)*xy(n2+k))/(sx(k)*sy(k))
        gn=sqrt(gn)
        ph=atan2(xy(n2+k),xy(k))*57.2958
        coh=abs(xy(k))/sqrt(sx(k)*sy(k))
        tk=float(k)*rb
        if(tk>1.e-2) then
         write(io,'(f12.3,2x,g14.4,2x,g14.4,2x,f5.2,2x,f6.1,2x,f5.2)')
     *tk,10.*alog10(sx(k)),10.*alog10(sy(k)),gn,ph,coh
          else
           write(io,'(g12.3,2x,g14.4,2x,g14.4,2x,f5.2,2x,f6.1,2x,f5.2)')
     *tk,10.*alog10(sx(k)),10.*alog10(sy(k)),gn,ph,coh
        endif
       enddo
      endif
c--------
      if(iunit==2) then
       rb=sr/float(lb)*1.e3
       write(io,'('' Frequency/Hz'',2x,''x Spectrum'',6x,''y Spectrum ''
     *,3x,''Gain'',3x,''Phase'',1x,''Coherence''/)')
       do k=il,iu
        gn=(xy(k)*xy(k)+xy(n2+k)*xy(n2+k))/(sx(k)*sy(k))
        gn=sqrt(gn)
        ph=atan2(xy(n2+k),xy(k))*57.2958
        coh=abs(xy(k))/sqrt(sx(k)*sy(k))
        tk=float(k)*rb
        if(tk>1.e-2) then
         write(io,'(f12.3,2x,g14.4,2x,g14.4,2x,f5.2,2x,f6.1,2x,f5.2)')
     *tk,10.*alog10(sx(k)),10.*alog10(sy(k)),gn,ph,coh
          else
           write(io,'(g12.3,2x,g14.4,2x,g14.4,2x,f5.2,2x,f6.1,2x,f5.2)')
     *tk,10.*alog10(sx(k)),10.*alog10(sy(k)),gn,ph,coh
        endif
       enddo
      endif
c--------
      if(iunit==0) then
       write(io,'(4x,''Period '',a7,1x,''x Spectrum'',6x,''y Spectrum ''
     *,3x,''Gain'',2x,''Phase'',1x,''Coherence''/)') su
       do k=il,iu
        rt=float(lb)*sr
        gn=(xy(k)*xy(k)+xy(n2+k)*xy(n2+k))/(sx(k)*sy(k))
        gn=sqrt(gn)
        ph=atan2(xy(n2+k),xy(k))*57.2958
        coh=abs(xy(k))/sqrt(sx(k)*sy(k))
        tk=rt/float(k)
        write(io,'(g14.4,3x,g14.4,2x,g14.4,1x,f5.2,1x,f6.1,2x,f5.2)') tk
     *,10.*alog10(sx(k)),10.*alog10(sy(k)),gn,ph,coh
       enddo
      endif
c-------
      return
      end
c
c===========================================
c
      subroutine outcis(c1,c2,xr,xi,br,ik,nbis,lb,su,sr,x1,x2,iunit,io)
c*******************************************
c  Output large cross bispectral values 8-20-97
c===========================================
      real xr(nbis),xi(nbis),br(nbis),fi(7),fo(7),z(7)
      integer c1(nbis),c2(nbis),ik(nbis),i1(25),i2(25),i3(25)
      character*10 su
      data fi/.1,.2,.4,.6,.8,.9,.99/
      intent(in) c1,c2,xr,xi,br,ik,nbis,lb,su,sr,x1,x2,iunit,io
c********************************************
      write(io,'('' No. of cross bispectral values = '',i6)')nbis
c  Write p-values for test stats for real and imag parts of cbis
      write(io,'(/2x,''p-value = '',f7.4,'' for test statistic that Re{x
     *xy(f1,f2)} = 0 for all f1,f2''/)')x1
      write(io,'(2x,''p-value = '',f7.4,'' for test statistic that Im{xx
     *y(f1,f2)} = 0 for all f1,f2''/)')x2
      write(io,'(2x,70(''='')/)')
      call sort(ik,br,nbis,7,fi,fo,'B',17)
c  Write fractiles
      write(io,'(12x,''Seven Fractiles of Cross Bispec Probs''/)')
      write(io,'(4x,''10%'',6x,''20%'',6x,''40%'',6x,''60%'',6x,''80%'',
     *6x,''90%'',6x,''99%''/)')
      write(io,'(7(f7.2,2x)/)')(fo(k),k=1,7)
      do k=1,7
       z(k)=(fo(k)-fi(k))/sqrt(fi(k)*(1.-fi(k))/float(nbis))
      enddo
      write(io,'(7(f7.2,2x),'' Z Values'')')(z(k),k=1,7)
      write(io,'(2x,70(''='')/)')
c-------
c  Output top cross bispectral values
      nn=nint(.02*float(nbis))
      nn=min(nn,25)
      nn=nbis-max(nn,4)+1
      nout=0
      pt=10./float(nbis)
      pt=min(pt,.01)
      do j=nbis,nn,-1
       ji=ik(j)
       pr=1.-br(ji)
       if(pr<=pt) then
        i1(nbis-j+1)=c1(ji)
        i2(nbis-j+1)=ji
        nout=nout+1
       endif
      enddo
c  Sort on 1st index if nout > 2
      if(nout>2) then
       call indexi(nout,i1,i3)
       write(io,'(7x,''Top '',f5.2,'' % Nonlinear Cross Bispectal Values
     * & Phases''/)')(float(nout)/float(nbis))*100.
       else
        return
      endif
c-------
      if(iunit==1) then
       rb=sr/float(lb)
       write(io,'(12x,'' Frequencies in kHz'')')
       write(io,'(1x,60(''='')/)')
       write(io,'(6x,''f1'',8x,''f2-f1'',7x,''f2'',10x,''Gain'',3x,''Pha
     *se'',4x,''1-Pr''/)')
       do j=1,nout
        k=i3(j)
        ji=i2(k)
        t1=float(i1(k))*rb
        t2=float(c2(ji))*rb
        fd=t2-t1
        pr=1.-br(ji)
        if(t1>1.e-2.or.t2>1.e-2) then
         write(io,'(1x,f9.3,2x,f9.3,2x,f9.3,3x,f9.2,2x,f6.1,3x,f7.5)')t1
     *,fd,t2,xr(k),xi(k),pr
          else
           write(io,'(1x,g9.3,2x,g9.3,2x,g9.3,3x,f9.2,2x,f6.1,3x,f7.5)')
     *t1,fd,t2,xr(k),xi(k),pr
        endif
       enddo
      endif
c--------
      if(iunit==2) then
       rb=sr/float(lb)*1.e3
       write(io,'(12x,'' Frequencies in Hz'')')
       write(io,'(1x,60(''='')/)')
       write(io,'(6x,''f1'',8x,''f2-f1'',7x,''f2'',10x,''Gain'',3x,''Pha
     *se'',4x,''1-Pr''/)')
       do j=1,nout
        k=i3(j)
        ji=i2(k)
        t1=float(i1(k))*rb
        t2=float(c2(ji))*rb
        fd=t2-t1
        pr=1.-br(ji)
        if(t1>1.e-2.or.t2>1.e-2) then
         write(io,'(1x,f9.3,2x,f9.3,2x,f9.3,3x,f9.2,2x,f6.1,3x,f7.5)')t1
     *,fd,t2,xr(k),xi(k),pr
          else
           write(io,'(1x,g9.3,2x,g9.3,2x,g9.3,3x,f9.2,2x,f6.1,3x,f7.5)')
     *t1,fd,t2,xr(k),xi(k),pr
        endif
       enddo
      endif
c--------
      if(iunit==0) then
       rt=float(lb)*sr
       write(io,'(17x,''Periods in '',a7)') su
       write(io,'(1x,60(''='')/)')
       write(io,'(7x,''Period 1'',8x,''Per(f2-f1)'',6x,''Period 2'',10x,
     *''Gain'',4x,''Phase'',4x,''1-Pr''/)')
       do j=1,nout
        k=i3(j)
        ji=i2(k)
        t1=rt/float(i1(k))
        t2=rt/float(c2(ji))
        fd=rt/(float(c2(k))-float(c1(k)))
        pr=1.-br(ji)
        if(t1>1.e-2.or.t2>1.e-2) then
         write(io,'(f14.3,2x,f14.3,2x,f14.3,6x,f9.2,2x,f6.1,3x,f7.5)')t1
     *,fd,t2,xr(k),xi(k),pr
          else
           write(io,'(5x,g11.3,5x,g11.3,5x,g11.3,4x,f9.2,3x,f6.1,3x,f7.5
     *)') t1,fd,t2,xr(k),xi(k),br(k)
        endif
        enddo
      endif
c-------
      return
      end
c
c===========================================
c
      subroutine wrtcis(c1,c2,xr,xi,br,nbis,lb,su,sr,iunit,io)
c*******************************************
c  Output cross bispectrum 12-2-2005
c===========================================
      real xr(nbis),xi(nbis),br(nbis)
      integer c1(nbis),c2(nbis)
      character*10 su
      intent(in) c1,c2,xr,xi,br,nbis,lb,su,sr,iunit,io
c********************************************
      if(iunit==1) then
       rb=sr/float(lb)
       write(io,'(12x,'' Frequencies in kHz'')')
       write(io,'(1x,60(''='')/)')
       write(io,'(6x,''f1'',8x,''f2-f1'',7x,''f2'',10x,''Gain'',3x,''Pha
     *se'',5x,''Pr''/)')
       do k=1,nbis
        t1=float(c1(k))*rb
        t2=float(c2(k))*rb
        fd=t2-t1
        if(t1>1.e-2.or.t2>1.e-2) then
        write(io,'(1x,f9.3,2x,f9.3,2x,f9.3,3x,f9.2,2x,f6.1,3x,f6.4)')t1,
     *fd,t2,xr(k),xi(k),br(k)
          else
           write(io,'(1x,g9.3,2x,g9.3,2x,g9.3,3x,f9.2,2x,f6.1,3x,f6.4)')
     *t1,fd,t2,xr(k),xi(k),br(k)
        endif
       enddo
      endif
c--------
      if(iunit==2) then
       rb=sr/float(lb)*1.e3
       write(io,'(12x,'' Frequencies in Hz'')')
       write(io,'(1x,60(''='')/)')
       write(io,'(6x,''f1'',8x,''f2-f1'',7x,''f2'',10x,''Gain'',3x,''Pha
     *se'',5x,''Pr''/)')
       do k=1,nbis
        t1=float(c1(k))*rb
        t2=float(c2(k))*rb
        fd=t2-t1
        if(t1>1.e-2.or.t2>1.e-2) then
        write(io,'(1x,f9.3,2x,f9.3,2x,f9.3,3x,f9.2,2x,f6.1,3x,f6.4)')t1,
     *fd,t2,xr(k),xi(k),br(k)
          else
           write(io,'(1x,g9.3,2x,g9.3,2x,g9.3,3x,f9.2,2x,f6.1,3x,f6.4)')
     *t1,fd,t2,xr(k),xi(k),br(k)
        endif
       enddo
      endif
c--------
      if(iunit==0) then
       rt=float(lb)*sr
       write(io,'(17x,''Periods in '',a5)') su
       write(io,'(1x,60(''='')/)')
       write(io,'(7x,''Period 1'',8x,''Per(f2-f1)'',6x,''Period 2'',10x,
     *''Gain'',4x,''Phase'',5x,''Pr''/)')
       do k=1,nbis
        t1=rt/float(c1(k))
        t2=rt/float(c2(k))
        fd=rt/(float(c2(k))-float(c1(k)))
        if(t1>1.e-2.or.t2>1.e-2) then
         write(io,'(f14.3,2x,f14.3,2x,f14.3,6x,f9.2,2x,f6.1,3x,f6.4)')
     *t1,fd,t2,xr(k),xi(k),br(k)
          else
           write(io,'(5x,g11.3,5x,g11.3,5x,g11.3,4x,f9.2,3x,f6.1,3x,f6.4
     *)')t1,fd,t2,xr(k),xi(k),br(k)
        endif
       enddo
      endif
c-------
      return
      end
c
c===========================================
c
      subroutine wrtcisn(b1,b2,br,u,iu,nbis,rb,iunit,io)
c*******************************************
c  Output cross bispectrum for plotting 12-3-2005
c===========================================
      real br(nbis),u(iu,iu)
      integer b1(nbis),b2(nbis)
      intent(in) b1,b2,br,iu,nbis,rb,iunit,io
c********************************************
      do k1=1,iu
       do k2=1,iu
        u(k1,k2)=0.
       enddo
      enddo
      do k=1,nbis
       u(b1(k),b2(k))=br(k)
      enddo
      do k2=1,iu
       if(iunit>0) then
        write(io,'(:,10000(f9.4))') (u(k1,k2),k1=1,iu)
         else
          write(io,'(:,10000(f14.3))') (u(k1,k2),k1=1,iu)
       endif
      enddo
      if(iunit==1) then
       write(io,'(:,10000(f9.3))') (float(k)*rb*1.e-3,k=1,iu)
        elseif(iunit==2) then
         write(io,'(:,10000(f9.3))') (float(k)*rb,k=1,iu)
          else
           write(io,'(:,10000(f14.3))') (rb/float(k),k=1,iu)
      endif
c-------
      return
      end
c
c===========================================
c
      subroutine datc(x,ts,fmt,mk,nmk,is,ie,n,nc,ii,io)
c********************************************************
c  Reads n data values from nc columns
c   x - input data matrix, ts - time marks if indicated by 'm'
c========================================================
      real x(n,nc)
      character*50 fmt
      character*25 ts(nmk)
      character mk
      intent(in) fmt,mk,nmk,is,ie,n,nc,ii,io
      intent(out) x,ts
c********************************************
      nt=0
      if(mk=='m') then
       ia=index(fmt,'(')
       ib=index(fmt,'a')
       if(ib==ia+1) then
c  Time mark column on left
        nt=1
         elseif(ib>ia+2) then
c  Time mark column on right
          nt=2
           elseif(nt==0) then
            write(io,'(7x,''/Data & time format is incorrect''/)')
            write(io,'(a50)')fmt
            write(6,'(7x,''/Data & time format is incorrect''/)')
            write(6,'(a50)')fmt
            write(109,'(7x,''/Data & time format is incorrect''/)')
            write(109,'(a50)')fmt
            stop
       endif
      endif
c-------
      rewind(ii)
c  Skip header and parameter records
      read(ii,*)
      read(ii,*)
c  Read all the data
      if(mk=='m') then
       if(is==1) then
        do i=1,ie
         if(nt==1) then
          read(ii,fmt,end=1,err=2,iostat=ko)ts(i),(x(i,k),k=1,nc)
         endif
         if(nt==2) then
          read(ii,fmt,end=1,err=2,iostat=ko)(x(i,k),k=1,nc),ts(i)
         endif
        enddo
       endif
       if(is>1) then
c  Read is:ie
        do i=1,is-1
         read(ii,fmt,end=1,err=2,iostat=ko)
        enddo
        do i=is,ie
         if(nt==1) then
          read(ii,fmt,end=1,err=2,iostat=ko)ts(i-is+1),(x(i-is+1,k),k=1,
     *nc)
         endif
         if(nt==2) then
          read(ii,fmt,end=1,err=2,iostat=ko)(x(i-is+1,k),k=1,nc),ts(i-is
     *+1)
         endif
        enddo
       endif
      endif
      if(mk/='m') then
       if(is==1) then
        do i=1,n
         read(ii,fmt,end=1,err=2,iostat=ko)(x(i,k),k=1,nc)
        enddo
       endif
       if(is>1) then
        do i=1,is-1
         read(ii,fmt,end=1,err=2,iostat=ko)
        enddo
        do i=is,ie
         read(ii,fmt,end=1,err=2,iostat=ko)(x(i-is+1,k),k=1,nc)
        enddo
       endif
      endif
      return
c--------
  1   write(io,*)'**** eof in read from data file '
      write(io,*)'i = ',i
      return
      write(6,*)'**** eof in read from data file '
      write(6,*)'i = ',i
      write(109,*)'**** eof in read from data file '
      write(109,*)'i = ',i
c  Error in data read
  2   write(io,*)'** Error in read of data file ',ko
      write(io,*)'Check format',fmt
      write(6,*)'** Error in read of data file ',ko
      write(6,*)'Check format',fmt
      write(109,*)'** Error in read of data file ',ko
      write(109,*)'Check format',fmt
      stop
      end
c
c===========================================
c
      subroutine fopen(name,buf,nfs,ia,ib,io)
c*******************************************
c  Opens data file, nfs - no. of data files (1 or 2)
c========================================================
      integer nfs,ia,ib,io
      character*700 name
      character*80 buf
      character*50 charead
      character*20 par
      character*2 delim
      logical exs
      intent(in) name,nfs
      intent(out) buf,ia,ib
c********************************************
c  Data file name(s) & path
      if(nfs==1) then
       par='ile name1'
      endif
      if(nfs==2) then
       par='ile name2'
      endif
      delim(1:2)='= '
      buf=charead(par,name,delim,ia,ier,3)
      if(ier>0) then
       write(io,*)'Error in reading the data file name from %.cnl'
       write(io,*)'Error number for charead = ',ier
       write(6,*)'Error in reading the data file name from %.cnl'
       write(6,*)'Error number for charead = ',ier
       write(9,*)'Error in reading the data file name from %.cnl'
       write(9,*)'Error number for charead = ',ier
       stop
      endif
c-------
      close(nfs)
c  Inquire if file exists
      inquire(file=buf,exist=exs)
      if(.not.exs) then
       write(3,*)'Data file ',buf,' does not exist'
       write(6,*)'Data file ',buf,' does not exist'
       write(9,*)'Data file ',buf,' does not exist'
       stop
      endif
c  Cbis file
      ia=index(buf,'.')
      if(ia==0) then
       write(6,*)'Data file name does not have an extension'
       write(6,*)buf
       write(9,*)'Data file name does not have an extension'
       write(9,*)buf
       stop
      endif
c-------
c  Parse for up to 5 path switches
      ib=index(buf,'\')
      if(ib>0) then
       ic=index(buf(ib:),'\')
       ib=ib+ic
       if(ic>0) then
        ic=index(buf(ib:),'\')
        ib=ib+ic
        if(ic>0) then
         ic=index(buf(ib:),'\')
         ib=ib+ic
        endif
        if(ic>0) then
         ic=index(buf(ib:),'\')
         ib=ib+ic
        endif
       endif
        else
         ib=1
      endif
c-------
c  Open data file
      open(nfs,file=buf,err=1,iostat=io1,status='old')
      rewind(nfs)
      return
c-------
  1   write(6,*)'**** Error on opening data file ',io1
      write(6,*)buf
      write(9,*)'**** Error on opening data file ',io1
      write(9,*)buf
      stop
      end
c
c===========================================
c
      subroutine csize(fl,fu,nbis)
c*******************************************
c  Computes cross bispectrum loops and gets array size
c   nbis - actual size of xr & xi arrays
c========================================================
      integer fl,fu
      intent(in) fl,fu
      intent(out) nbis
c********************************************
      nbis=0
c  Triangle fl<k2<fu, 2*k2<k1<k2
      do k2=fl,fu
       k1=k2/2
       if(2*k1==k2.and.k1>=fl) then
c  Cbis for k2=2*k1
        nbis=nbis+1
       endif
c-------
c  Rest of first triangle
       i1=k1+1
       do k1=i1,k2-1
        ks=k2-k1
        if(ks>=fl.and.ks<=fu) then
         nbis=nbis+1
        endif
       enddo
c-------
c  Cbis for triangle 2*fl<k2<fu, k2<k1<fu
       do k1=k2+1,fu
        ks=k1-k2
        if(ks>=fl.and.ks<=fu) then
         nbis=nbis+1
        endif
       enddo
      enddo
c  End of k2 loop
c-------
      return
      end
c
c===================================================
c
      subroutine indexi(n,irin,indx)
c********************************************
c  Sorts real array arrin in ascending order using an indexed Heapsort
c  Input: irin  Version: 4-11-95
c============================================
      integer irin(n),indx(n),q
c********************************************
      if(n==1) return
      do j=1,n
       indx(j)=j
      enddo
      l=n/2+1
      ir=n
  2   continue
        if(l>1)then
          l=l-1
          indxt=indx(l)
          q=irin(indxt)
        else
          indxt=indx(ir)
          q=irin(indxt)
          indx(ir)=indx(1)
          ir=ir-1
          if(ir==1)then
            indx(1)=indxt
            return
          endif
        endif
        i=l
        j=l+l
  3     if(j<=ir)then
          if(j<ir)then
           if(irin(indx(j)).lt.irin(indx(j+1)))j=j+1
          endif
          if(q.lt.irin(indx(j)))then
            indx(i)=indx(j)
            i=j
            j=j+j
          else
            j=ir+1
          endif
        goto 3
        endif
        indx(i)=indxt
      goto 2
c-------
      end