PRO scn1raf,infile,hr1,hr2,iverb=iverb ; Identify RAFOS regions, find peaks, determine height and arrival time. ; RWB, January 8, 2002. ; Write out a 1-line result to file scn1raf1.tmm; RWB, January 23, 2002. 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) !p.background='ffffff'x; white !p.color=0; black IF n_params() eq 0 THEN BEGIN; if no file on command line, get one infile='/usr/data/bland/pioneer/2001/r0127403.26d'; 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 fname=strmid(infile,strlen(infile)-12,12); just the file name IF n_params() lt 3 THEN BEGIN; If no time window given, prompt for it. hr1=4.19; hr2=4.22; read,'RAFOS arrival-time window starts at '+string(hr1)+': ',str1 IF strlen(str1) gt 0 THEN hr1=float(str1) read,'RAFOS arrival-time window ends at '+string(hr2)+': ',str1 IF strlen(str1) gt 0 THEN hr2=float(str1) ENDIF IF n_elements(iverb) eq 0 THEN iverb=3; verbosity ; iverb=0: scan silently, writing to output file. ; iverb=1: scan continuously, writing lines but without graphs. ; iverb=2: scan continuously, writing lines and displaying graphs. ; iverb=3: scan and wait for prompt from console to continue. iyear = 2000 + fix(strmid(infile,strlen(infile)-11,2)) iday = fix(strmid(infile,strlen(infile)-9,3)) ; print,'iyear, iday = ',iyear, iday 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 rsig = fltarr(npt); signal array readu, lun1, rsig close,lun1 free_lun,lun1 dt = double(.001)*lblk/86400.; time per point tsig = dindgen(npt)*dt+day0; time array (sec) hrsig=(tsig mod 1.)*24. IF iverb ge 2 THEN plot,hrsig,rsig,title=fname IF iverb ge 3 THEN read,'Hit return to continue',str1 ; **** Find out if we have a usable time series. hra=hrsig(0) hrb=hrsig(npt-1) IF hra lt hr1 and hrb gt hr2 THEN BEGIN hstart=hr1 hend=hr2 ENDIF ELSE IF hra lt hr1+12. and hrb gt hr2+12. THEN BEGIN hstart=hr1+12. hend=hr2+12. ENDIF ELSE BEGIN print,'The search block is not included in this time series.' return ENDELSE ; **** Search block is included. dummy=min((hrsig-hstart)^2,ihstart) dummy=min((hrsig-hend)^2,ihend) IF iverb ge 2 THEN BEGIN ; **** Plot +/- 2 min about center of search block, light gray. hrplt1=(hstart+hend)/2.-0.0333333333; lower bouud hrplt2=hrplt1+0.06666666666 dummy=min((hrsig-hrplt1)^2,iplt1) dummy=min((hrsig-hrplt2)^2,iplt2) plot,hrsig(iplt1:iplt2),rsig(iplt1:iplt2),title=fname, $ /nodata,xstyle=1 oplot,hrsig(iplt1:iplt2),rsig(iplt1:iplt2),color='c0c0c0'x ; **** Plot +/- 80 sec about center of search block, darker gray. hrplt1=(hstart+hend)/2.-0.0222222222; lower bouud hrplt2=hrplt1+0.04444444444 dummy=min((hrsig-hrplt1)^2,iplt1) dummy=min((hrsig-hrplt2)^2,iplt2) oplot,hrsig(iplt1:iplt2),rsig(iplt1:iplt2),color='808080'x ; **** Plot search band, black. oplot,hrsig(ihstart:ihend),rsig(ihstart:ihend) IF iverb ge 3 THEN read,'Hit return to continue',str1 ENDIF ts = (tsig(ihstart:ihend) mod 1.)*24. rs = rsig(ihstart:ihend) npts=ihend-ihstart+1 ; **** Find the highest peak in the selected region. rsmax=max(rs,irsmax); find index of peak in rs array. irsigmax=irsmax+ihstart; index of peak in rsig array. ; **** Determine the average noise near this peak. This will be the average ; **** of the signal in two 80-second blocks , located to the right and left ; **** of the peak, at the edges of a +/-80 sec band around the peak which ; **** might be affected by the RAFOS signal itself. dirafos=long(80)*1000/lblk; number of points in one RAFOS signal ib11=max([0,irsigmax-2*dirafos]); start of left-hand band ib12=irsigmax-dirafos; end of left-hand band ib21=irsigmax+dirafos; start of right-hand band ib22=min([npt-1,irsigmax+2*dirafos]); end of right-hand band ; print,'ib11,ib12,ib21,ib22 = ',ib11,ib12,ib21,ib22 npband=0 IF ib12 ge ib11 THEN npband=npband+ib12-ib11+1; points from left-hand band IF ib22 ge ib21 THEN npband=npband+ib22-ib21+1; points from right-hand band sigband=0. IF ib12 ge ib11 THEN sigband=sigband+total(rsig(ib11:ib12)); LH band IF ib22 ge ib21 THEN sigband=sigband+total(rsig(ib21:ib22)); RH band rsigns = sigband/max([npband,1]); average noise rsigns = max([rsigns,.0001]); Don't let the noise vanish. SoverN = rsmax/rsigns; signal-to-noise ratio ; **** Find the half-height points and so estimate the rms, assuming a ; **** Gaussian. hval = (rsmax+rsigns)/2. ;** Walk out to the right-hand side looking for the half-height point. FOR istep=1,20 DO BEGIN IF rsig(irsigmax+istep) lt hval THEN BREAK ENDFOR ;** Interpolate. thalfr=tsig(irsigmax+istep-1)+dt*(rsig(irsigmax+istep-1)-hval)/ $ (rsig(irsigmax+istep-1)-rsig(irsigmax+istep)) ;** Walk out to the left-hand side looking for the half-height point. FOR istep=1,20 DO BEGIN IF rsig(irsigmax-istep) lt hval THEN BREAK ENDFOR ;** Interpolate. thalfl=tsig(irsigmax-istep+1)+dt*(rsig(irsigmax-istep+1)-hval)/ $ (rsig(irsigmax-istep+1)-rsig(irsigmax-istep)) thalfsec=(thalfr-thalfl)/2.*86400. tsigma = thalfsec/1.18; estimator of the Gaussian sigma. ; print,'thalfr, thalfl, thalfsec ,tsigma = ',thalfr,thalfl,thalfsec,tsigma, $ ; format = '(A,2f20.12,2f10.6)' ; **** Determine the arrival time for this signal. nout=5 nout = fix(3.*thalfsec/(.001*lblk)+0.5); number of points out i1dt=max([0,irsigmax-nout]); start of first trial band for tpeak and twidth i2dt=min([npt-1,irsigmax+nout]); end of first trial band ndt = i2dt-i1dt+1; number of points in first trial band rsigsum=total(double(rsig(i1dt:i2dt))); r sum tav=total((tsig(i1dt:i2dt)-iday)*double(rsig(i1dt:i2dt)))/rsigsum ; tsqav=total((tsig(i1dt:i2dt)-iday)^2*double(rsig(i1dt:i2dt)))/rsigsum ; trms=sqrt(tsqav-tav^2); rms width of the peak ; print,'ndt,rsigsum,tav,tsqav,trms = ',ndt,rsigsum,tav,tsqav,trms ; print,'tav, tsqav, trms = ',tav,tsqav,trms,format='(A,3f20.10)' ; twidth=trms*86400.; Convert to seconds. ipeakday=fix(tav); day of peak arrival tpeakhr=(tav mod 1.)*24.; peak arrival time in hours ipeakhr=fix(tpeakhr); hour of peak arrival tpeakmin=(tpeakhr mod 1.)*60. ipeakmin=fix(tpeakmin); minute of peak arrival tpeaksec=(tpeakmin mod 1.)*60. ; print, 'Peak arrived on day ',strtrim(string(iday),2), $ ; ' of ',strtrim(string(iyear),2), $ ; ' at ',strtrim(string(ipeakhr),2), $ ; ':',strtrim(string(ipeakmin),2), $ ; ':',strtrim(string(tpeaksec),2), ' ; HWHM = ', $ ; strtrim(string(thalfsec),2), ' sec.' scnline=string(FORMAT='(i5,i4,i3,i3,f8.4,f8.4,2f10.4,f14.8)', $ iyear,iday,ipeakhr,ipeakmin,tpeaksec,min([99,thalfsec]), $ min([9999.,rsmax]),min([9999.,rsigns]), $ tav+iday) ; print,iyear,iday,ipeakhr,ipeakmin,tpeaksec,min([99,thalfsec]), $ ; min([9999.,rsmax]),min([9999.,rsigns]), $ ; tav+iday,format='(i5,i4,i3,i3,f8.4,f8.4,2f10.4,f14.8)' ; print,'For comparison, the scanline is' print, scnline ;***** Write the scan line to a temporary communication file. openw,lunscan,'scn1raf1.tmp',/get_lun printf,lunscan,scnline close,lunscan free_lun,lunscan IF iverb ge 2 THEN BEGIN di=50; spread for blowup i1=max([0,irsmax-di]) i2=min([irsmax+di,npts-1]) tss=ts(i1:i2); blowup block (time in hr) tss=(tss-tss(0))*3600.; convert to sec, starting at zero rss=rs(i1:i2) plot,tss,rss,title=fname+' '+strtrim(string(ipeakhr),2)+':'+ $ strtrim(string(ipeakmin),2)+':'+strtrim(string(tpeaksec),2)+ $ ' HWHM: '+strtrim(string(thalfsec),2), $ xtitle='time (sec) (add '+strtrim(string(ts(i1)*3600.),2)+')', $ /nodata oplot,tss,rss,color='aaaaaa'x; plot in gray i3=max([0,irsmax-nout-i1]) i4=min([irsmax+nout-i1,i2-i1]) oplot,tss(i3:i4),rss(i3:i4),thick=2; highlight the peak (+/- 3 sigmas) IF iverb ge 3 THEN read,'Hit return to continue',str1 ENDIF END