PRO PLOTBPC,BPC_FILE,CHAN,OPTION,START ;+ ; NAME: ; PLOTBPC ; PURPOSE: ; Read the BPCxxxxxx.xxxx files written by BSDCAL and the BSD data file ; specified within it and plot the spectral fits as given by parameters ; in the BPC file over the raw data ; ; CALLING SEQUENCE: ; PLOTBPC[,bpc_file,chan,option,start] ; ; ; INPUTS: ; BPC_FILE= Equal the BPC file name or I.D. ; If BPC_FILE is undefined or the first character is ; blank, PLOTBPC will prompt user for name. ; ; CHAN = required channel; ; ; OPTION = the fit option number required to display. If supplied ; on the command line then only one option is allowed ; ; START = the time of first spectrum required ( default = all) ; ; COMMON BLOCKS: ; None. ; ; SIDE EFFECTS: ; None. ; ; MODIFICATION HISTORY: ; 0) Written May 1991 CDP ; 1) Revised to handle start time of spectra better Sep 91 CDP ; ; 2) Undo previous "bug-fix" and fix real one where rdmap variable ; was incorrectly incremented between spectra Nov 91 CDP ; ; 3) Cut out bin-width scaling (not needed once BSDCAL sorted) Nov 91 CDP ; ; 4) Generalize and UNIXify handling of file names. Mar 92 CDP ; 5) Changed BGRND call to routine POLY. June 92, RDB ; 6) Use rd_bsd_head to access BSD file header July 92 CDP ; ;-------------------------------------------------------------------------- ; set up arrays required ;----------------------- dummy = lonarr(5) rdummy = fltarr(4) time = intarr(6) ; spectrum start time (BPC file) stime = lonarr(7) ; data start time (BSD file) dw_etc = fltarr(6) ; dw, origwo, wo, wooff, sdgi, count rate num = intarr(2) ; no. options, no. back cs = intarr(2) ; channel & spectrum bin = intarr(2) ; start & end bins lparm = fltarr(6) ; line parameters backparm = fltarr(6) ; background parameters optname = ' ' ; character option name data = fltarr(26,200,10) ; array for all information ; dimensions are (parameter,spectrum,option) sarray = intarr(200) ; array for spectrum numbers loaded x = fltarr(500) ; array to hold bin values y = fltarr(500) ; array to hold data values temp = fltarr(16) ; fudge to get over read limitations bcof = fltarr(3) ; background coefficients vcof = fltarr(9) ; parameters for Voigt profile calculation fwhml = fltarr(4) ; Lorentz width for each channel wave = fltarr(4) ; wavelength of resonance line ewid = fltarr(4) ; electronic resolution of each channel tmp_val = fltarr(6) ; yet another dummy xann = fltarr(2) ; x coordinates of annotation yann = fltarr(2) ; y coordinates of annotation maxopt = intarr(4) ; max number options in each channel num_opt = intarr(10) ; array for input of options wanted tchar = ' ' ; character variable for input time charchan = ' ' ; string for channel number input ans = ' ' ; one letter answer from terminal ; initialise integers used in FOR loops (version 2 much fussier than v 1!) ; ------------------------------------------------------------------------ nreq_opt = 0 lines = 0 kspec = 0 nrecs = 0 kopt = 0 nbinsp = 0 ; load known values ;------------------ fwhml = [ 2.08e-4, 2.06e-4, 3.89e-4, 13.71e-4 ] ; Lorentzian width ewid = [ 0.0, 0.0, 0.0, 0.0 ] ; electronic width wave = [ 1.778, 1.8509, 3.1781, 5.0385 ] ; wavelength resonance maxopt = [5, 9, 3, 7] ; max number options ; in each channel clear_screen = '\033[;H\033[2J' ; clear page alpha = '\x18' ; reset to alpha screen bell = '\007' ; warning true = 1 ; psuedo logical false = 0 ; ditto nparm = n_params() ; set variable so can ; be changed to loop ; ; compile to get message away from plot ; dumx = calfun() ;;dumx = bgrnd() if nparm lt 3 then option = 0 ; if not on line, initialise if nparm eq 0 then bpc_file = ' ' ; if not even the file ; name given ;print, clear_screen ; clear screen if (strmid(bpc_file,0,1) eq ' ') or $ ; if not valid BPC file (strlen(bpc_file) eq 0) then begin ; name then manual input print,'$(a)', ' Enter the BPC file name (omit BPC prefix): ' bpc_file = ' ' & read,bpc_file endif if strlen(bpc_file) eq 0 then return ; trap exit if strpos(strupcase(bpc_file),'bpc') lt 0 then $ ; if BPC prefix not bpc_file = 'bpc'+bpc_file ; there, affix it ;------------------------------------------------------------------------ ; procedure loops back to here when request satisfied or user Quit - ; requested during plotting - ;------------------------------------------------------------------------ continue: data(*,*,*) = 0 ; for safety get_lun, lun ; get free unit number ; to attach file openr, lun, bpc_file ; open for readonly bsd_name = ' ' ; get the bsd file name from readf,lun,bsd_name ; the top line of the BPC bsd_name = strmid(bsd_name,$ ; file strpos(bsd_name,'BSD file')+8,$ strlen(bsd_name)) break_file,bsd_name,disk_log,dir_log,fil,ext,fver bsd_name = strtrim(fil + ext,2) buf = ' ' ; get the BPC file format readf,lun,buf ; code mode = fix(strtrim(strmid(buf,$ strpos(buf,'mode:')+5,$ strlen(buf)),2)) if mode eq 0 then begin ; if not 1, then error tex = ' Expected BPC format mode = 1.'$ + ' Program aborted.' print, bell print,'$(/a)',tex return endif ;------------------------------------------------------------------- ; read header of BSD file to get start time ;------------------------------------------------------------------- rd_bsd_header,bsd_name,bsd_h ; ; calculate day of the mission ; dom,bsd_h.tdf(4),bsd_h.tdf(5),bsd_h.tdf(6),day_of_mission ; time of first spectrum in mission seconds bsd_stime = long(day_of_mission)*86400 + $ long(bsd_h.tdf(0))*3600 + bsd_h.tdf(1)*60 + bsd_h.tdf(2) ;------------------------------------------------------------------- ; get channel wanted ;------------------------------------------------------------------- err_flag = 1 ; force at least one go while err_flag eq 1 do begin if (nparm lt 2) then begin ; channel not supplied print,'$(a)', ' Enter the channel number: (Q to quit)' read,charchan ; character input if strupcase(charchan) eq 'Q' then begin ; null response quits close, lun return endif else begin chan = fix(charchan) ; otherwise translate endelse endif err_flag = 0 ; reset flag if (chan lt 1 or chan gt 4) then begin ; set sensible limits err_flag = 1 endif endwhile ;----------------------------------------------------------------------- ; get options list wanted ;----------------------------------------------------------------------- err_flag = 1 ; force at least one go while err_flag eq 1 do begin if (nparm lt 3) or (option eq 0) then begin ; option not on call line nreq_opt = fix(0) free_in, chan, num_opt, nreq_opt ; free format input endif else begin nreq_opt = 1 num_opt(0) = option endelse err_flag = 0 ; reset error flag for i = 1,nreq_opt do begin ; check each is OK if (num_opt(i-1) lt 1) or $ (num_opt(i-1) gt maxopt(chan-1)) then begin ; limit errors err_flag = 1 option = 0 ; error, set flag endif endfor endwhile ;----------------------------------------------------------------------- ; get time of first spectrum wanted for plotting ;----------------------------------------------------------------------- colpos1 = -1 colpos2 = -1 if (nparm lt 4) then begin ; start time not on call line while colpos1 lt 0 do begin ; so get it here print,'$(a,i2,a,i2,a,i2,a,i2,a,i2,a,i4)',$ ' First spectrum has start time (hh:mm:ss) of ',bsd_h.tdf(0),':',$ bsd_h.tdf(1),':',bsd_h.tdf(2),' on date ',bsd_h.tdf(4),'/',$ bsd_h.tdf(5),'/',bsd_h.tdf(6) print,'$(a)',$ ' Enter time of the first spectrum required (hh:mm[:ss])',$ ' - default gives all -' read,tchar if strlen(tchar) eq 0 then begin tchar = string(bsd_h.tdf(0),':',bsd_h.tdf(1),':',bsd_h.tdf(2),$ "(i2,a,i2,a,i2)") endif colpos1 = strpos(tchar,':') endwhile colpos2 = strpos(tchar,':',colpos1+1) ; search for second colon for ; seconds field endif thr = fix(strmid(tchar,0,colpos1)) ; extract hours tmin = fix(strmid(tchar,colpos1+1,2)) ; extract minutes if colpos2 gt -1 then begin tsec = fix(strmid(tchar,colpos2+1,2)) ; extract seconds if there endif else begin tsec = 0 endelse start = long(thr)*3600 + tmin*60 + tsec ; start time in seconds start = start + long(day_of_mission)*86400 ; add start day if start lt bsd_stime then begin ; if requested start time start = start + 86400 ; is less than first spectrum endif ; assume next day is wanted print,'$(///2a)',' Spectral data from file: ',bsd_name print,'$(2a)',' Reading results file ',strupcase(bpc_file) print,'$(a20,i2,a10,10i2)',' Looking for channel ',chan,$ ' option #s',num_opt(0:nreq_opt-1) erase ; clear graphics screen print, alpha ; return to alpha readf,lun,buf ; read 1 more header record of ----'s kspec = -1 ; initialize spectrum total counter sp_skip = 0 ; number of spectra time-skipped inform_opt = 0 ; has missing options been notified? while not eof(lun) do begin ; loop while not at end of file for i=0,1 do begin readf,lun,buf ; skip 2 lines below the -----'s endfor newdata = true readf,lun,cs,time,dw_etc,num ; read header line dom,time(3),time(4),time(5),day_of_mission ; get day of mission of spectrum tsec = long(time(0))*3600 + time(1)*60 + time(2)$ + long(day_of_mission)*86400 ; time in seconds if cs(0) eq chan and (tsec ge start) then begin ; is it the desired ; channel and is it ; after specified ; start time? for i=0,2 do begin readf,lun,buf ; skip 3 lines of headers etc endfor kopt = -1 ; running total of options ; found for iopt = 1,num(0) do begin ; for all the options ; in this channel extra = false ; not yet read extra lines bin = [0,0] lparm = [0.0,0.0,0.0,0.0,0.0,0.0] ; zero out arrays backparm = [0.0,0.0,0.0,0.0,0.0,0.0] tempback = fltarr(num(1)*2) ; create temp array for backg readf,lun,opt,okflag,lines,bin,rchi,$ lparm,tempback,optname ; read data line backparm(0:num(1)*2-1) = tempback ; store backgnd for ireq=0,nreq_opt-1 do begin ; loop over requested options if opt eq num_opt(ireq) then begin ; is this a correct option? mode = 1 ; set normal mode if newdata then begin kspec = kspec+1 ; increment spectrum total newdata = false ; start of a new parameter endif kopt = kopt + 1 ; options stored data(0,kspec,kopt) = cs(1) ; save the spectrum number data(1,kspec,kopt) = bin(0) ; save the start bin number data(2,kspec,kopt) = bin(1) ; save the end bin number data(3:8,kspec,kopt) = time(0:5) ; save the time data(9, kspec,kopt) = backparm(0) ; save background constant data(10,kspec,kopt) = backparm(2) ; save background slope data(11,kspec,kopt) = backparm(4) ; save background curvature data(12,kspec,kopt) = lparm(0) ; save first line's height data(13,kspec,kopt) = lparm(2) ; save first line's width data(14,kspec,kopt) = lparm(4) ; save first line's position if lines gt 1 then begin ; if there are more lines ; in this fit extra = true ; flag: have read them for i=0,lines-2 do begin if opt eq 2 then begin ; is this opt 2 (vel)? mode = 2 readf,lun,tmp_val ; read additional lines data(15+i*3,kspec,kopt) = tmp_val(0) ; save next line's intensity data(16+i*3,kspec,kopt) = tmp_val(4) ; save next line's vel data(17+i*3,kspec,kopt) = 0.0 ; spare endif else begin readf,lun,lparm ; read additional lines data(15+i*3,kspec,kopt) = lparm(0) ; save next line's height data(16+i*3,kspec,kopt) = lparm(2) ; save next line's width data(17+i*3,kspec,kopt) = lparm(4) ; save next line's posn endelse endfor endif data(21,kspec,kopt) = mode ; save the mode for ; this option data(22,kspec,kopt) = lines ; save # of lines ireq = fix(nreq_opt - 1) ; skip the rest data(23,kspec,0) = kopt ; save # options in spectrum data(24,kspec,kopt) = okflag ; was fit OK? data(25,kspec,kopt) = dw_etc(0) ; save dw for this one endif endfor if not extra then begin for i=0,lines-2 do begin readf, lun,buf ; if necessary, skip other endfor ; lines endif endfor if inform_opt eq 0 then begin if (kopt+1) lt nreq_opt then begin print, bell, ' NOTE: not all the options requested were found.' inform_opt = 1 endif endif readf,lun,buf ; end of spectrum options so ; skip the ---'s line endif else begin ; not channel wanted or ; spectrum too early if cs(0) eq chan and (tsec lt start) then begin ; count skipped spectra sp_skip = sp_skip + 1 endif buf = ' ' minibuf = ' ' ; so skip till ------'s line while minibuf ne '----------' and not eof(lun) do begin readf,lun,buf ; search for next line of minibuf = strmid(buf,2,10) ; --------'s endwhile endelse endwhile if sp_skip gt 0 then begin ; inform of spectra skipped if sp_skip eq 1 then begin print, ' One spectrum skipped prior to start time.', bell endif else begin print, sp_skip,' spectra skipped prior to start time.', bell endelse endif if kspec eq -1 then begin ; inform if no spectra found print, bell, ' ** ERROR ** no spectra for that' print, ' channel and options-list located.', bell return endif else begin if kspec eq 0 then begin print, ' Only one spectrum loaded.' ; tidy in case only one endif else begin print, kspec+1,' spectra loaded.' ; inform of success rate endelse endelse wait,3.0 ; wait for user digestion close, lun ; finished with BPC file free_lun, lun ; release unit sarray = data(0,*,0) ; extract spectrum #s ssarray = intarr(kspec+1) ; sorted spectrum numbers ssarray = sort(sarray(0:kspec)) ; do the sort ;----------------------------------------------------------------------- ; now access the BSD data file ;----------------------------------------------------------------------- get_lun, lun ; get a unit for the BSD data file openr, lun, bsd_name ; unformatted, readonly numchan = lonarr(5) ; total spectra and number per channel ptrdmp = lonarr(4) ; roadmap pointers for each channel readu, lun,numchan,ptrdmp ; read first header record in file numspec = numchan(chan) ; total number of spectra in this channel rdmap = ptrdmp(chan-1) ; the record where roadmap for this channel is rdmap = (rdmap-1) * 64 ; byte location of start of roadmap if kspec+1 gt numspec then begin ; check number of spectra OK print, bell, ' ** ERRROR ** number of spectra referenced in BPC file' print, ' is greater than those present in BSD file', bell return endif for is=0,kspec do begin ; loop over spectra to be plotted sis = ssarray(is) ; order as sorted above numspec = fix( data(0,sis,0) ) - 1 ; spectrum to plot (zero start) nrecs = fix(numspec/16) ; there are 16 spectral pointers per record point_lun, lun,rdmap+nrecs*64 ; go to where required spectrum pointer ; is located, ready for reading pointers = lonarr(16) ; temp storage for pointers readu, lun,pointers ; read them spcrec = pointers(numspec-nrecs*16) ; finally the record where ; required spectrum is located spcbyt = (spcrec-1)*64 ; in bytes to beginning of record wanted point_lun, lun,spcbyt ; go there ready to read readu, lun,pointers ; read first housekeeping record nbinsp = pointers(14) ; to get the number bins per spectrum readu, lun,pointers ; skip second housekeeping record nrecs = nbinsp/16 ; no. of complete records used by bins etc nleft = nbinsp - nrecs*16 ; any leftovers from complete records for i=1,nrecs do begin ; read the bin values readu, lun,temp ; one record at a time start = (i-1)*16 ; where to start loading temp into x x(start:start+15) = temp ; load x endfor if nleft gt 0 then begin ; have leftovers to deal with readu, lun,temp ; read 'em start = nrecs*16 ; where to load into x x(start:start+nleft-1) = temp(0:nleft-1) ; load 'em endif for i=1,nrecs do begin ; read the spectral data values readu, lun,temp ; one record at a time start = (i-1)*16 ; where to start loading temp into y y(start:start+15) = temp ; load y endfor if nleft gt 0 then begin ; have leftovers to deal with readu, lun,temp ; read 'em start = nrecs*16 ; where to load into y y(start:start+nleft-1) = temp(0:nleft-1) ; load 'em endif ;----------------------------------------------------------------- ; attempt to plot data and overplot the fit ;----------------------------------------------------------------- kopt = fix(data(23,sis,0)) ; max # of options ; in this spectrum minbin = 255 ; inverse limits to maxbin = 0 ; find real limits for i=0,kopt do begin if data(1,sis,i) lt minbin then begin ; find min minbin = data(1,sis,i) endif if data(2,sis,i) gt maxbin then begin ; and max maxbin = data(2,sis,i) ; bins to be used endif endfor xs = minbin-1 xe = maxbin-1 ; set data limits ; -1 to account for ; zero start arrays ; ; not my idea! repeat this bit if want hardcopy after seeing terminal ; version ; hardcopy: ; set up plot title and information lines ;----------------------------------------- if nreq_opt gt 1 then begin ; more than 1 option? opt_str = string(' Options: ',num_opt(0),"(a,i1)") ; start string off for i=2,nreq_opt do begin temp_str = string(', ',num_opt(i-1),"(a,i1)") ; next option opt_str = opt_str + temp_str ; add it endfor endif else begin opt_str = string(' Option: ',num_opt(0),"(a,i1)") ; only one option endelse plt_title = string('BSDCAL results for Channel: ',chan,' Spectrum: ',$ numspec+1,opt_str,"(a,i1,a,i3,a)") time_tit = string('Time: ',data(3,sis,0),'h ',data(4,sis,0),'m ',$ data(5,sis,0),'s on ',data(6,sis,0),'/',data(7,sis,0),$ '/',data(8,sis,0),"(a,i2,a,i2,a,i2,a,i2,a,i2,a,i4)") ; plot raw data ;-------------- !p.linestyle = 0 ; solid line plot,x(xs:xe),y(xs:xe),title=plt_title,$ xtitle='Bin value',ytitle='Counts/sec/unit bin position' ; plot data xwid = (!x.crange(1) - !x.crange(0))*0.05 ; plot width xyouts, !x.crange(0)+xwid,0.95*!y.crange(1),bsd_name ; bsd file xyouts, !x.crange(0)+xwid,0.9*!y.crange(1),time_tit ; time info xann(0) = !x.crange(0)+xwid ; legend xann(1) = !x.crange(0)+xwid*2 ; x posn yann(0) = 0.8*!y.crange(1) ; and yann(1) = 0.8*!y.crange(1) ; y posn oplot, xann,yann ; output line xyouts, xann(1)+0.1*xwid,0.99*yann(1),'raw data' ; annotate ; calculate the fit for each option used ;--------------------------------------- for k=0,kopt do begin lines = fix(data(22,sis,k)) ; lines in this option bcof = data(9:11,sis,k) ; load background coefficients okflag = data(24,sis,k) ; recall fit flag this_dw = data(25,sis,k) ; recall dw for this fit for i=1,lines do begin ; the coefficient vector vcof((i-1)*3) = data(12+(i-1)*3,sis,k) ; height vcof((i-1)*3+1) = data(13+(i-1)*3,sis,k) ; width vcof((i-1)*3+2) = data(14+(i-1)*3,sis,k) ; position endfor pxmin = data(1,sis,k) - 1 ; xmin for this option pxmax = data(2,sis,k) - 1 ; xmax for this option if okflag gt 0 then begin !p.linestyle = 2 ; dashed line for OK fit endif else begin !p.linestyle = 1 ; dotted for bad fit endelse if vcof(0) gt 0.0 then begin ; check sensible mode = data(21,sis,k) ; set mode for this option yfit = calfun(x,vcof,lines,fwhml(chan-1),$ this_dw,mode,wave(chan-1),ewid(chan-1))$ + poly(x,bcof) ;;bgrnd(x,bcof) oplot, x(pxmin:pxmax),yfit(pxmin:pxmax) ; overlay calculated fit endif if(k eq 0) then begin yann(0) = 0.75*!y.crange(1) ; legend yann(1) = 0.75*!y.crange(1) ; y posn oplot, xann,yann ; output line !p.linestyle = 0 ; solid line xyouts, xann(1)+0.1*xwid,0.99*yann(1),'calculated fit' ; annotate endif if okflag eq 0 then begin yann(0) = 0.65*!y.crange(1) ; legend yann(1) = 0.65*!y.crange(1) ; y posn oplot, xann,yann ; output line !p.linestyle = 0 ; solid line xyouts, xann(1)+0.1*xwid,0.99*yann(1),'fit not converged' endif ; ; plot the background used ;------------------------- yfit = poly(x,bcof) ;;bgrnd(x,bcof) !p.linestyle = 1 oplot, x(pxmin:pxmax),yfit(pxmin:pxmax) if k eq 0 then begin yann(0) = 0.7*!y.crange(1) ; legend yann(1) = 0.7*!y.crange(1) ; y posn oplot, xann,yann ; output line !p.linestyle = 0 ; solid line xyouts, xann(1)+0.1*xwid,0.99*yann(1),'background used' ; annotate endif endfor ; end of options in this spectrum !p.linestyle = 0 ; solid line type if ans eq 'H' then begin ; back to terminal set_plot,save_dev ans = ' ' endif else begin xyouts, 0.0,0.01,$ ' Hit any key to proceed: (Q to quit: H for hard copy)',/normal ans = get_kbrd(1) ; wait for reply endelse ans = strupcase(ans) if ans eq 'Q' then begin is = kspec ; kill loop by endif ; resetting counter if ans eq 'H' then begin ; like it! want xyouts, 0.15,0.5 ,$ ; hard copy ' HARD COPY being prepared in PostScript file (IDL.PS)',/normal save_dev = !d.name set_plot,'PS' goto, hardcopy endif endfor close, lun ; close BSD file free_lun,lun ; free unit nparm = 0 ; kid it that it needs channel number etc. print, alpha ; switch back to alpha screen print, ' ' goto, continue ; naughty but nice! end