pro srsp_mf_plot_,filename,hrec,findex ; 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 date = (srsp_tl_decode(header)).date ; Read the observing sequence record (comes directly after ; the header) and return the OBSEQ structure. obseq = srsp_get_obseq(getdata(rec+1,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 = hrec+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 = lonarr(nmax) ; Long integer array to hold the data atndat = rawdat 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 ; use this to determine the factor to multiply the data by. atnbits = ishft(thesedat,-12) atnfac = (10^(atnbits/2.0)) ; Limit the attenuation correction to a factor of 100 (20 DB), ; setting the attenuation factor to unity for larger values. bad = where(atnbits ge 5,nbad) if (nbad gt 0) then atnfac(bad) = 1. ; Correct the data and place the result in RAWDAT rawdat[k:k+495] = long((thesedat and 'FFF'x)*atnfac) atndat[k:k+495] = long(atnbits) ; 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. ncycle = n_elements(rawdat)/480L ; Reform the single long RAWDAT array into a three-dimensional IMG array with ; the first dimension being the sample type with the following defn ; 0 = RCP, ND-OFF ; 1 = LCP, ND-OFF ; 2 = RCP, ND-ON ; 3 = LCP, ND-ON ; the second dimension is the frequency (given in *obseq.pf), and the ; last dimension being the times (given in the array TIME). imgR = reform(rawdat[0:ncycle*480L-1],240,2,ncycle) rawdat = 0 ; free memory of rawdat imgL = imgR[indgen(120)*2+1,*,*] ; Extract all of the LCP data imgR = imgR[indgen(120)*2,*,*] ; Extract all of the RCP data atnimg = reform(atndat[0:ncycle*480L-1],240,2,ncycle) if (!d.name eq 'WIN') then begin ; Set up the screen for display device,decomp=0 window,/free,xsiz=894,ysiz=400 wid = !d.window erase,255 window,/free,xsiz=894,ysiz=130,/pixmap wpix = !d.window erase,255 ymin=140 endif else begin ; Set up the printer for display device,scale=2.35 ymin=200 endelse if (!d.name eq 'WIN') then wset,wid ; Show the RCP data as an image, with logarithmic scaling tvscl,congrid(rotate(alog10((reform(imgR[*,0,*])+50-reform(imgR[*,0,3000<(ncycle-1)])#replicate(1,ncycle))>10),-4),800,120),70,ymin+130 tvscl,congrid(rotate(alog10((reform(imgL[*,0,*])+50-reform(imgL[*,0,3000<(ncycle-1)])#replicate(1,ncycle))>10),-4),800,120),70,ymin ; plots,/dev,[0,69],(130+ymin+findex)*[1,1],color=0 ; plots,/dev,[871,894],(130+ymin+findex)*[1,1],color=0 ; plots,/dev,[0,69],(ymin+findex)*[1,1],color=0 ; plots,/dev,[871,894],(ymin+findex)*[1,1],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 while (curgood ge -1) do begin cursor,x,y,0,/dev 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 findex = y-270 curgood = 1 end y gt 0 and y le 20 and x gt 0 and x le 40: curgood = -99 else: curgood = -1 endcase if (curgood ge 0) then begin iplt = lindgen(n_elements(imgR[0,0,*]))*120L + findex p1 = reform(imgR[findex,0,*]-imgR[findex,0,500<(ncycle-1)]) p2 = reform(imgL[findex,0,*]-imgL[findex,0,500<(ncycle-1)]) maxy = max(p1)>max(p2) miny = min(p1)max(imgL[findex,0,*]) utplot,time[iplt]/1000.,imgR[findex,0,*],date,color=0,xsty=1,yran=[0,maxy],ysty=1 outplot,time[iplt]/1000.,imgL[findex,0,*],date,color=250 xyouts,/norm,0.9,0.15,string((*obseq.pf)[findex],format='(I5)')+' GHz',color=0,align=1,charsize=1.0 endelse if (!d.name eq 'WIN') then begin wdelete,wpix wdelete,wid endif !p = psav loadct,0 ptr_free,obseq.pf,obseq.ppac,obseq.pndon,obseq.pcalon return end