pro srsp_mf_plot,filename,hrec,findex,doprint=doprint,auto=auto doprint = 0 ; Open file. Note that NREC is returned, which is the ; total number of records in the file. lun = srsp_openarc(filename,a,nrec) if (lun le 0) then begin return endif ; Read header record rec = hrec-1 header = getdata(rec,a) ; Determine date from the header record headtls = srsp_tl_decode(header) date = headtls.date ; Read the observing sequence record (comes directly after ; the header) and return the OBSEQ structure. rec = rec+1 obseq = srsp_get_obseq(rec,a) ; Save the plot environment so that changes to environment ; can be restored later psav = !p loadct,39 ; We can call this with the WIN device (output to screen) ; or the PRINTER device (output to printer). Make changes ; appropriate to screen if device is WIN if (!d.name ne 'WIN') then begin !p.multi=[1,1,3] !p.charsize=2 endif ; Read the first data record (IREC will be incremented below ; to point to current data record). irec = rec+1 data = uint(getdata(irec,a)) ; Read the TIME/LABEL part of data record (first 32 bytes) ; and return structure in TLS variable. tls = srsp_tl_decode(data) ; Read each record in turn, putting all of the data into ; a single, very long array called RAWDAT. ; NMAX is the total expected length of data elements. There ; are 248 pairs of data measurements (total of 496 2-byte integers), ; each pair is RCP, then LCP. This makes 992 bytes, which, when ; added to the 32 byte TIME/LABEL field accounts for the 1024 bytes ; of each data record. nmax = (nrec-hrec)*496L time = lonarr(nmax/2) ; Time array only half as long because the data are in pairs rawdat = uintarr(nmax) ; Long integer array to hold the data atndat = rawdat xtrafac = [1.,1.] istart = irec tstart = tls.msec ; Start time is time in TIME/LABEL structure for record while (tls.segmentcode eq !seg.data) do begin if (tls.msec lt tstart) then begin ; Looks like the time is out of sequence--probably a zero meaning we have ; come to the end of the data, so fudge an end-of-sequence (EOS) record. tls.segmentcode = !seg.eos endif else begin ; Set the start time for this record to the time in the TIME/LABEL structure tstart = tls.msec k = (irec-istart)*496L ; This is a pointer to the current point in RAWDAT kt = k/2 ; Only half as many times (because data are in pairs) time[kt:kt+247] = lindgen(248)*10 + tstart ; Fill in next 248 times with 10 msec ; increment to TSTART ; Read data (skipping the 32 bytes of the TIME/LABEL field) ; into a temporary array called THESEDAT thesedat = data[indgen(496)+16] ; Any data with hex *FFF are bad data. Find the bad data bad = where((thesedat and '0fff'x) eq 'fff'x,nbad) ; Set any bad data to zero if (nbad ne 0) then thesedat[bad] = 0 ; The first four bits of each data sample are a code for the ; current attenuation setting (an integer in units of 5 DB), so ; strip these off and save them. atndat[k:k+495] = uint(ishft(thesedat,-12)) ; The last 12 bits represent the data, so isolate these bits ; and put into RAWDAT array rawdat[k:k+495] = uint((thesedat and 'FFF'x)) ; Get the next record of data and repeat irec = irec + 1 data = uint(getdata(irec,a)) if (data[0] eq -1) then tls.segmentcode = !seg.eos else tls = srsp_tl_decode(data) endelse endwhile ; We have reached the EOS record (or blank or otherwise bad data) ; so close the file and free then logical unit free_lun,lun ; Eliminate the first measurement sample, which is always not real data ; and is not part of the observing sequence time = time[1:kt+247] rawdat = rawdat[2:k+495] atndat = atndat[2:k+495] ; We know that there are always 240 pairs of samples in each cycle ; through the frequencies, so determine how many cycles there are. nfp = obseq.nx ncycle = n_elements(rawdat)/(nfp*2L) ; Reform the single long RAWDAT array into a three-dimensional IMG ; array, where the first dimension is the frequency (given in ; *obseq.pf), the second dimension is the noise diode state ; 0 = ND-OFF, 1 = ND-ON, and the last dimension is the times ; (given in the array TIME). imgR = reform(rawdat[0:ncycle*nfp*2L-1],nfp,2,ncycle) rawdat = 0 ; free memory of rawdat imgL = imgR[indgen(nfp/2)*2+1,*,*] ; Extract all of the LCP data imgR = imgR[indgen(nfp/2)*2,*,*] ; Extract all of the RCP data atnimgR = reform(atndat[0:ncycle*nfp*2L-1],nfp,2,ncycle) atndat = 0 ; free memory of atndat atnimgL = atnimgR[indgen(nfp/2)*2+1,*,*] atnimgR = atnimgR[indgen(nfp/2)*2,*,*] atnimgR = shift(atnimgR,1,0,0);<4 atnimgL = shift(atnimgL,1,0,0);<4 roff = 10.0 loff = 7.0 imgR = (imgR-roff)*10^(atnimgR/2.0) imgL = (imgL-loff)*10^(atnimgL/2.0) if (!d.name eq 'WIN') then begin ; Set up the screen for display device,decomp=0 ytop=440 window,/free,xsiz=894,ysiz=ytop ; ysiz changed from 400 (cgm) wid = !d.window erase,255 window,/free,xsiz=894,ysiz=130,/pixmap wpix = !d.window erase,255 xmin=0 ymin=140 endif else begin ; Set up the printer for display ;device,scale=2.35,/portrait ; xmin=5 ; ymin=200 device,scale=3.29,/landscape xmin=-20 ymin=210 endelse if (!d.name eq 'WIN') then wset,wid ; open calibration file and read cal data openr,/get_lun,lun,'c:\working\fluxcal.fit' cal = fltarr(6,500) cal1 = fltarr(6) i = 0 while (not eof(lun)) do begin readf,lun,cal1 cal[*,i] = cal1 i = i + 1 endwhile cal = cal[*,0:i-1] free_lun,lun ; Get observing frequencies and find indexes of these frequencies ; corresponding to the cal data nf = obseq.nx/2 f = (*obseq.pf)[0:nf-1] ; Get frequency list from OBSEQ structure indx = intarr(nf) for i = 0, nf-1 do begin indx[i] = where(cal[0,*] eq f[i]) endfor ; If the calibration file does not contain a frequency in F, then ; INDX is set to -1, so find any of those and set the CAL data ; to NaN bad = where(indx eq -1,nbad) if (nbad ne 0) then cal[1:2,indx[bad]] = !values.f_nan ; IMGR has size (NF,2,NTIMES), where the middle index =0 for ND-OFF, =1 for ND-ON imgR = reform(imgR[*,0,*]) imgL = reform(imgL[*,0,*]) ; apply calibration imgR = imgR*(reform(cal[2,indx])#replicate(1,ncycle)) imgL = imgL*(reform(cal[1,indx])#replicate(1,ncycle)) ; Show the RCP data as an image, with logarithmic scaling ; take log tempR=alog10((imgR+50-imgR[*,3000<(ncycle-1)]#replicate(1,ncycle))>10) tempL=alog10((imgL+50-imgL[*,3000<(ncycle-1)]#replicate(1,ncycle))>10) ; transpose and sample tempR=congrid(rotate(tempR,-4),800,120) tempL=congrid(rotate(tempL,-4),800,120) ; set max and min, and plot using tv (for plots with same max and min for LH and RH if keyword_set(auto) then begin minv=min(tempR) < min(tempL) ; use these 2 lines for automatic scaling maxv=max(tempR) > max(tempL) endif else begin minv=1.0 ; these 2 lines set min to 10. and max to 3162. maxv=3.5 endelse rngv=maxv-minv tv,!D.TABLE_SIZE*(tempR-minv)/rngv,70+xmin,ymin+130 tv,!D.TABLE_SIZE*(tempL-minv)/rngv,70+xmin,ymin ; next 2 lines are for automatic independent scaling of LH and RH panels ; tvscl,congrid(rotate(alog10((imgR+50-imgR[*,3000<(ncycle-1)]#replicate(1,ncycle))>10),-4),800,120),70+xmin,ymin+130 ; tvscl,congrid(rotate(alog10((imgL+50-imgL[*,3000<(ncycle-1)]#replicate(1,ncycle))>10),-4),800,120),70+xmin,ymin ; this block added by cgm 5/15/2003 (misc labels) flab=intarr(4) for i=0,2 do flab(i)=min(where(f gt 1000*2^(i+1))) flab(2)=flab(2)-1 flab(3)=119 for i=0,3 do begin xyouts,/dev,25,flab(i)+ymin,string(f(flab(i)),format='(i5)'),color=0,align=0,charsize=1.0 xyouts,/dev,25,flab(i)+ymin+130,string(f(flab(i)),format='(i5)'),color=0,align=0,charsize=1.0 endfor xyouts,/dev,400,ymin+281,'SRSP DATA PLOTS, Bell Labs, Murray Hill, NJ',color=0,charsize=1.5,align=0.5 xyouts,/dev,400,ymin+261,date+' [Day '+strmid(tls.yrday,5,3)+']',color=0,charsize=1,align=0.5 xyouts,/dev,20,ymin+10,'Frequency, MHz',color=0,orient=90,charsize=1. xyouts,/dev,20,ymin+140,'Frequency, MHz',color=0,orient=90,charsize=1. xyouts,/dev,75,ymin+251,'Right hand channel',color=0,charsize=1. xyouts,/dev,75,ymin+121,'Left hand channel',color=250,charsize=1. colorbar,XRANGE=[10^minv,10^maxv],/xlog,XTIT='SFU',color=0 if (!d.name eq 'WIN') then begin wset,wid curgood = 1 polyfill,/dev,[0,0,40,40,0],[0,20,20,0,0],color=0 xyouts,/dev,20,5,'Exit',color=255,align=0.5,charthick=2 polyfill,/dev,[0,0,40,40,0],[0,20,20,0,0]+ytop-21,color=0 xyouts,/dev,20,ytop-16,'Print',color=255,align=0.5,charthick=2 mousestate = 0 xyouts,20,42,/dev,'Mouse',align=0.5,color=0 xyouts,10,30,/dev,'off',color=0 while (curgood ge -1) do begin cursor,x,y,0,/dev if (!mouse.button eq 1) then begin mousestate = 1-mousestate if (mousestate eq 0) then begin xyouts,10,30,/dev,'on' xyouts,10,30,/dev,'off',color=0 endif else begin xyouts,10,30,/dev,'off' xyouts,10,30,/dev,'on',color=0 endelse wait,0.5 endif else if (!mouse.button eq 4 and x ne -1) then begin imid = nint((x - 70)*ncycle/800.) istart = (imid-250)>0<(ncycle-1) iend = (imid+249)<(ncycle-1)>0 if (iend-istart ge 100) then begin wdisplay,headtls,f,time[istart*240L:iend*240L+239],imgr[*,istart:iend],imgl[*,istart:iend] endif else print,'SRSP_MF_PLOT: Error, plot range too small.' wait,0.5 wset,wid endif if (mousestate) then begin case 1 of y lt 259 and y ge 140: begin findex = y-140 curgood = 1 end y ge 270 and y lt 389: begin if (x gt 40) then begin findex = y-270 curgood = 1 endif end y gt 0 and y le 20 and x gt 0 and x le 40: curgood = -99 y gt (ytop-20) and y le (ytop-1) and x gt 0 and x le 40: begin doprint = 1 goto,bail end else: curgood = -1 endcase if (curgood ge 0) then begin iplt = lindgen(n_elements(imgR[0,*]))*240L + findex p1 = reform(imgR[findex,*]-imgR[findex,3000<(ncycle-1)]) p2 = reform(imgL[findex,*]-imgL[findex,3000<(ncycle-1)]) good = where(finite(p1),ngood) if (ngood gt 2) then begin maxy = max(p1)>max(p2) miny = min(p1)max(p2) miny = min(p1)