PRO mkhok,infile,outfile,idstart,ihstart,imnstart,delay
; Carry out cross-correlation detection for a new experimental swept-
; frequency RAFOS-type source on Hoke Seamount, using data from a test
; carried out on November 9, 2001.  The signals are about 134 seconds
; long, sweeping from 201 Hz to 301 Hz (estimated).
; Read in a 2-hour raw-data file; select a 268-sec section and separate
; the channels; frequency-shift, down sample and cross-correlate each channel
; separately.  Write out the results.
;	r0131319.46m --> r0131319.54H
; The format of the output file is:
;	dstart	(double)	start time for the first
;	nptd	(long)		number of detection points
; cross-correlation point (not necessarily the first point in the block
; of raw data selected
;	dc(nptd,4)	(float)  detector output, cosine phase
;	ds(nptd,4)	(float)  detector output, sine phase
;	dmag(nptd,4)	(float)  quadrature sum of dc and ds
; RWB, May 2, 2002.
  npar=n_params()
  init,npar,infile,outfile,idstart,ihstart,imnstart,delay; Startup tasks.
  getfile; Read in data, select the segment to operate on.
  calc1; Carry out the calculations.
  outfile; Write the output file.
  endall; Clean up for exit.
END


PRO init,npar,inf,outf,ids,ihs,imns,del
  common params,infile,outfile,idstart,ihstart,imnstart,delay,dstart
  common rafpar,f0,f1,phi1,dfdt,tsig,dtside
  common datin,npt,idat
  common datout,nblk,nraf,nptd,scdat,ssdat,sdscdat,sdssdat,q1raf, $
    q2raf,dc,ds,dmag
  str1=''; buffer
  IF npar eq 6 THEN BEGIN
    infile=inf; input file name
    outfile=outf; output file name
    idstart=ids; start of signal, day
    ihstart=ihs; start of signal, hour
    imnstart=imns; start of signal, minute
    delay=double(del); start of signal, seconds
  ENDIF ELSE BEGIN
    infile='/usr/data/mild/pioneer/2001/r0131319.46m'; PSM 2-hr raw-data file
    ;***** Ask for file name *****
    read,'infile:  '+strmid(infile,strlen(infile)-12,12)+' : ',str1
    IF strlen(str1) gt 0 THEN strput,infile,str1,strlen(infile)-12

    outfile='/usr/data/bland/pioneer/2001/r0131319.54H'; Hoke output file
    read,'outfile: '+strmid(outfile,strlen(outfile)-12,12)+' : ',str1
    IF strlen(str1) gt 0 THEN strput,outfile,str1,strlen(outfile)-12

; **** For the Hoke test, the day is always 313.
    idstart=313
    read,'idstart: '+string(idstart)+' : ',str1
    IF strlen(str1) gt 0 THEN idstart=fix(str1)

; **** Input the hour and minute and second of the time of arrival
; **** for the particular signal desired.

    ihstart=19
    read,'ihstart: '+string(ihstart)+' : ',str1
    IF strlen(str1) gt 0 THEN ihstart=fix(str1)

    imnstart=54
    read,'imnstart: '+string(imnstart)+' : ',str1
    IF strlen(str1) gt 0 THEN imnstart=fix(str1)

    delay=double(47.); signal arrival time, seconds
    read,'delay: '+string(delay)+' : ',str1
    IF strlen(str1) gt 0 THEN delay=double(str1)

;    tsig=double(134.); duration of signal, in seconds
;    read,'tsig (sec): '+string(tsig)+' : ',str1
;    IF strlen(str1) gt 0 THEN tsig=double(str1)
  ENDELSE
  f0=double(260.); reference frequency
  f1=double(201.2); starting frequency
  phi1=double(0.); phase offset
  dfdt=double(0.74072)
  tsig=double(131.07); length of swept-frequency signal (sec)
  dtside=double(10.); time on either side of expected arrival time
END


PRO getfile
  common params,infile,outfile,idstart,ihstart,imnstart,delay,dstart
  common rafpar,f0,f1,phi1,dfdt,tsig,dtside
  common datin,npt,idat
  nblk=long(1801); This is the number of blocks on standard PSM files.
  lblk=long(4096); number of  samples in a block
  nsamp=nblk*lblk; total number of samples in file
  idatfile=intarr(nsamp,4); arrays for the four hydrophones
  print,'infile: ',infile
  openr, lun1, infile, /get_lun
  blk=assoc( lun1, bytarr(32786) ); blocks with 18-byte header
  FOR iblk=0,nblk-1 DO BEGIN
    b=blk(iblk); get one block
    IF iblk eq 0 THEN BEGIN; Get starting time for file from header.
      day0=dofyear(b(0:17))
      print,'Header ',string(b(0:16)),' --> day of year = ',day0, $
	format='(3a,f14.9)'
    ENDIF
    vladata = fix(b, 18, 4, 4096 )
    FOR ichan=0,3 DO BEGIN
      idatfile(iblk*lblk:(iblk+1)*lblk-1,ichan)=vladata(ichan,*); Store in array.
    ENDFOR
  ENDFOR
  close,lun1
  free_lun,lun1
; **** Select the part corresponding to the RAFOS signal.
  dstsig=double(idstart)+double(ihstart)/24.+double(imnstart)/1440. + $
    double(delay)/86400.
  dstart=dstsig-dtside/double(86400.); start of section to cross correlate
  print,'day0, dstsig, tsig, dtside, dstart = ',day0,dstsig,tsig,dtside,dstart
  IF day0 gt dstart THEN STOP; Error in selecting the file.
  npt=long((tsig+2.*dtside)*1000.); number of points to take
  i1=round((dstart-day0)*86400000.); first sample to take
  i2=i1+npt-1; last sample to take
  print,'npt, i1, i2 = ',npt,i1,i2
  idat=idatfile(i1:i2,*)
  idatfile=0; convert a long array to a single number, to save space
END

PRO calc1
  common params,infile,outfile,idstart,ihstart,imnstart,delay,dstart
  common rafpar,f0,f1,phi1,dfdt,tsig,dtside
  common datin,npt,idat
  common datout,nblk,nraf,nptd,scdat,ssdat,sdscdat,sdssdat,q1raf, $
    q2raf,dc,ds,dmag
; **** Carryout the frequency shifting and downsampling.  The frequency
; shifting and each subsequent operation are done twice, for the
; cosine phase and sine phase of the detection.

  lblk=long(4); averaging block length, corresponding to 4 ms
  nblk=npt/lblk; number of blocks
  sdscdat=fltarr(nblk,4); Frequency-shifted, down-sampled waveforms (cos).
  sdssdat=fltarr(nblk,4); Frequency-shifted, down-sampled waveforms (sin).
  scdat=fltarr(npt,4); Frequency-shifted waveforms, phase 1 (cos).
  ssdat=fltarr(npt,4); Frequency-shifted waveforms, phase 1 (cos).
; **** Calculate the carriers for frequency shifting.
  csig=2.*cos(f0*2.*!dpi*double(lindgen(npt))/1000.); carrier, phase 1 (cos)
  ssig=2.*sin(f0*2.*!dpi*double(lindgen(npt))/1000.); carrier, phase 2 (sin)

  FOR ichan=0,3 DO BEGIN
    scdat(*,ichan)=csig*idat(*,ichan); Mix data with carrier, cos phase.
    ssdat(*,ichan)=ssig*idat(*,ichan); Mix data with carrier, sin phase.
    FOR iblk=long(0),nblk-1 DO BEGIN; Carry out downsampling.
      sdscdat(iblk,ichan)=total(scdat(iblk*lblk:(iblk+1)*lblk-1,ichan))/lblk
      sdssdat(iblk,ichan)=total(ssdat(iblk*lblk:(iblk+1)*lblk-1,ichan))/lblk
    ENDFOR
  ENDFOR

;***** Calculate the quadrature replicas for detection.
;***** With quadrature-detected input signals, we only need q1raf.
  dt = double(.001)*lblk; time per point
  rafsig,f0,f1,phi1,dfdt,tsig,dt,traf,q1raf,q2raf; quadrature signals

  nraf=n_elements(traf); length of the RAFOS signals.
  nptd=nblk-nraf+1
  dc=fltarr(nptd,4); cross-correlated waveforms, cosine phase
  ds=fltarr(nptd,4); cross-correlated waveforms, sine phase
  dmag=fltarr(nptd,4); cross-correlation amplitude
  td=dindgen(nptd)*dt; time array for cross-correlated waveforms, in seconds
  print,'About to start the cross-correlation.  nraf, nptd = ',nraf,nptd
  FOR ichan=0,3 DO BEGIN; Loop over hydrophones.
    FOR icc=long(0),nptd-1 DO BEGIN
      ic1=icc
      ic2=icc+nraf-1
      dc(icc,ichan) = total(q1raf*sdssdat(ic1:ic2,ichan))/nraf
      ds(icc,ichan) = total(q1raf*sdscdat(ic1:ic2,ichan))/nraf
    ENDFOR
  ENDFOR
  dmag=sqrt(dc^2+ds^2); Calculate the amplitudes.
  save
END

PRO outfile
;***** Write out the result array.
  common params,infile,outfile,idstart,ihstart,imnstart,delay,dstart
  common datin,npt,idat
  common datout,nblk,nraf,nptd,scdat,ssdat,sdscdat,sdssdat,q1raf, $
    q2raf,dc,ds,dmag
  print,'outfile = ',outfile
  openw,lund,outfile,/get_lun
  writeu,lund,dstart
  writeu,lund,nptd
  writeu,lund,dc,ds,dmag
  close,lund
  free_lun,lund
  scratchf='hoke.dat'
  openw,lunout,scratchf,/get_lun
  writeu,lunout,npt,nblk,nraf,nptd,idat,sdscdat,sdssdat,q1raf,q2raf,dc,ds,dmag
  close,lunout
  free_lun,lunout
END

PRO endall
END

PRO rafsig,f0,f1,phi1,dfdt,ton,delt,tsig,q1sig,q2sig
; Calculate a RAFOS signal.  A direct and a quadrature signal (phase
; retarded by pi/2) are calculated.
;	f0:  reference frequency
;	f1:  Sweep starts at this frequency.
;	phi1: phase offset of RAFOS signal rel. to reference
;	dfdt:  sweep rate
;	ton:  time signal is on
;	delt:  time interval between samples
;	tsig:  time array for output
;	rsig;  output signal array
; The frequency as a function of time from turnon is given by
;	f(t) = f1 + dfdt*t - f0
; and the phase as a function of time from turnon is given by
;	phi = 2.*!pi*(f1-f0)*t + !pi*dfdt*t^2
  npt=long(ton/delt); number of points
  t=dindgen(npt)*delt; times for calculated signal
  phi = phi1 + 2.*!pi*(f1-f0)*t + !pi*dfdt*t^2
  q1sig=sin(phi)
  q2sig=sin(phi-!dpi/2.)
  tsig=t
END


PRO get1blk,lun1,iblk,vladata
; Read one block of 4096 samples; average the four channels.
   blk=assoc( lun1, bytarr(32786) )
   b=blk[iblk]
   vladata = fix(b, 18, 4, 4096 )
   data = transpose( temporary(vladata))
   byteorder, vladata, /sswap
END

FUNCTION dofyear,b
; Calculate a double-precision day of the year, from a byte array b
; containing a string of the form
;    ddd:hh:mm:ss.sssS[LF]
; This is the time string from Pioneer Seamount data files.
; RWB, September 29, 2001
   return, double(string(b(0:2)))+ $
	 1./24.*(double(string(b(4:5))) + $
	 1./60.*(double(string(b(7:8))) + $
	 1./60.* double(string(b(10:15))) ) )
END