PRO get2hrs,fname,inpath,outpath
; Get a signal array (beam-formed) from a single 2-hour block of PSM data.
; Write the array in a file.
;	r0124506.27m -> r0124506.27b (beam-formed signal)
; Then carry out quadrature mixing at 260 Hz, down sample, 
; and write to a file:
;		     -> r0124506.27w
; Beam forming and downsampling done by averaging (not summing), and
; mixing done with a sine wave of amplitude 2, as of Dec. 30, 2000.
; RWB, November 3, 2001.
; Changed to cos(omega*t) and sin(omega*t), RWB, January 10, 2002.
  str1=''; buffer
  day0=double(0); starting day for file (double-precision floating)
  IF n_params() lt 3 THEN BEGIN; need to get paths and file name
    fname='r0124506.27m'; name of a PSM 2-hr raw data file
    inpath='/usr/data/mild/pioneer/2001/'; path for "m" files
    outpath='/usr/data/bland/pioneer/2001/'; path for "b" and "w" files
    ;***** Ask for file name *****
     read,'input file: '+fname+' : ',str1
     IF strlen(str1) gt 0 THEN fname=str1
    ;***** Ask for input path *****
     read,'input path: '+inpath+' : ',str1
     IF strlen(str1) gt 0 THEN inpath=str1
    ;***** Ask for output path *****
     read,'output path: '+outpath+' : ',str1
     IF strlen(str1) gt 0 THEN outpath=str1
  ENDIF
  infile=inpath+fname; full input file name
  outfileb=outpath+fname; beam-formed-signal file name
  strput,outfileb,'b',strlen(outfileb)-1
  outfilew=outfileb; Webb-filtered-signal file name
  strput,outfilew,'w',strlen(outfilew)-1
  
  print,'infile:   ',infile
  print,'outfileb: ',outfileb
  print,'outfilew: ',outfilew

  nblk=long(1801); This is the number of blocks on standard PSM files.
  lblk=long(4096); samples per block
  nsamp=nblk*lblk; number of samples per 2-hr file (=7,376,896)
  idat=intarr(nsamp)
  idat1blk=intarr(lblk)
  openr, lun1, infile, /get_lun
  FOR iblk=0,nblk-1 DO BEGIN
    IF iblk eq 0 THEN BEGIN
      blk=assoc( lun1, bytarr(32786) )
      b=blk[iblk]
      day0=dofyear(b(0:17))
      print,string(b(0:16)),' day of year = ',day0,format='(2a,f14.9)'
    ENDIF
    get1blk,lun1,iblk,idat1blk,dayofyr; Read in one block, form beam.
    idat(iblk*lblk:(iblk+1)*lblk-1)=idat1blk; Store in array.
  ENDFOR
  close,lun1
  free_lun,lun1
  openw,lunout,outfileb,/get_lun
  writeu,lunout,day0,idat
  close,lunout
  free_lun,lunout
;***** Mix and filter for Webb sources.

;   s1=2.*cos(260.*2.*!pi*(lindgen(nsamp) mod 50)/1000.); carrier, phase 1
;   s2=2.*sin(260.*2.*!pi*(lindgen(nsamp) mod 50)/1000.); carrier, phase 1

;   s1=s1*idat; Mix data with phase-1 carrier.
;   s2=s2*idat; Mix data with phase-2 carrier.
   lblk=long(128); averaging block, about 1/8 second
   nblk=nsamp/lblk; number of blocks
;***** Carry out mixing and downsampling for the phase-1 carrier.
   sig=2.*cos(260.*2.*!pi*(lindgen(nsamp) mod 50)/1000.); carrier, phase 1
   sig=sig*idat; Mix data with carrier.
   webs1=fltarr(nblk)
   FOR iblk=long(0),nblk-1 DO BEGIN
      webs1(iblk)=total(sig(iblk*lblk:(iblk+1)*lblk-1))/lblk
   ENDFOR
;***** Carry out mixing and downsampling for the phase-2 carrier.
   sig=2.*sin(260.*2.*!pi*(lindgen(nsamp) mod 50)/1000.); carrier, phase 2
   sig=sig*idat; Mix data with carrier.
   webs2=fltarr(nblk)
   FOR iblk=long(0),nblk-1 DO BEGIN
      webs2(iblk)=total(sig(iblk*lblk:(iblk+1)*lblk-1))/lblk
   ENDFOR
;***** Write out the mixed, downsampled time series.
   openw,lunout,outfilew,/get_lun
   writeu,lunout,day0,lblk,nblk,webs1,webs2
   close,lunout
   free_lun,lunout
;   save
END

PRO get1blk,lun1,iblk,vlabeam,dayofyr
; Read one block of 4096 samples; average the four channels.
   blk=assoc( lun1, bytarr(32786) )
   b=blk[iblk]
   data = fix(b, 18, 4, 4096 )
   data = transpose( temporary(data))
   byteorder, data, /sswap
   vlabeam=fix((long(data(*,0))+long(data(*,1))+long(data(*,2))+ $
      long(data(*,3)))/4)
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