PRO rafdet,infile
; Carry out cross-correlation detection for RAFOS sources.
; Read in a 2-hour file of beam-formed, mixed, down-sampled data, e.g.
;	 r0124506.27w
; Then carry out quadrature cross-correlation with a RAFOS signal.
; RWB, December 29, 2001.
  str1=''; buffer
  lblk = long(0); length of averaging block used in downsampling
  npt = long(0); number of points in data arrays
  day0=double(0); starting day for file (double-precision floating)
  IF n_params() eq 0 THEN BEGIN; if no file on command line, get one
    infile='/usr/data/bland/pioneer/2001/r0124506.27w'; name of a PSM 2-hr data
					; file (beam-formed, mixed, downsampled)
    ;***** Ask for file name *****
     read,'input file: '+strmid(infile,strlen(infile)-12,12)+' : ',str1
     IF strlen(str1) gt 0 THEN strput,infile,str1,strlen(infile)-12
  ENDIF
  print,'infile:   ',infile
  outfiled=infile; detected-signal file name
  strput,outfiled,'d',strlen(outfiled)-1
  print,'outfiled: ',outfiled

  openr, lun1, infile, /get_lun
  readu, lun1, day0; time of first point
  print,'day0 = ',day0
  readu, lun1, lblk; number of 1000-Hz points per averaged point
  print,'lblk = ',lblk
  readu, lun1, npt; number of points
  print,'npt = ',npt
  dt = double(.001)*lblk; time per point
  rsig = fltarr(npt); signal array
  tsig = dindgen(npt)*dt; time array (sec)
  readu, lun1, rsig

;***** Calculate the quadrature replicas for detection.
  f0=double(260.); reference frequency
  f1=double(259.375); starting frequency
  phi1=double(0.); phase offset
  tramp=double(80.); ramp time in seconds
  delf=1.523; frequency sweep
  rafsig,f0,f1,phi1,delf/tramp,tramp,dt,traf,q1raf,q2raf; quadrature signals

  nraf=n_elements(traf)
  ncc=npt-nraf+1; number of cross-correlation points
  cc1=fltarr(ncc)
  cc2=cc1
  ccquad=cc1
  tcc=findgen(ncc)*dt; time in seconds
  FOR icc=long(0),ncc-1 DO BEGIN
    ic1=icc
    ic2=icc+nraf-1
    cc1(icc) =total(q1raf*rsig(ic1:ic2))/nraf
    cc2(icc)=total(q2raf*rsig(ic1:ic2))/nraf
  ENDFOR
  ccquad=sqrt(cc1^2+cc2^2)
;***** Write out the result array.
  print,'outfiled = ',outfiled
  openw,lund,outfiled,/get_lun
  writeu,lund,day0,lblk,ncc,ccquad
  close,lund
  free_lun,lund

  save
  return
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