From d54fe7c1f704a63824c5bfa0ece65245572e9b27 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 4 Mar 2015 21:21:30 -0500 Subject: Initial commit --- idl/cf_obsplot.pro | 334 ++++++++++++++++++++++++++++++++++++++++++++ idl/cf_plot_extract3.pro | 210 ++++++++++++++++++++++++++++ idl/cf_plot_rate3.pro | 227 ++++++++++++++++++++++++++++++ idl/cf_plot_rate_linear.pro | 175 +++++++++++++++++++++++ idl/fuse_extract.pro | 154 ++++++++++++++++++++ idl/fuse_sic_only.pro | 159 +++++++++++++++++++++ 6 files changed, 1259 insertions(+) create mode 100644 idl/cf_obsplot.pro create mode 100644 idl/cf_plot_extract3.pro create mode 100644 idl/cf_plot_rate3.pro create mode 100644 idl/cf_plot_rate_linear.pro create mode 100644 idl/fuse_extract.pro create mode 100644 idl/fuse_sic_only.pro (limited to 'idl') diff --git a/idl/cf_obsplot.pro b/idl/cf_obsplot.pro new file mode 100644 index 0000000..9a4c8c2 --- /dev/null +++ b/idl/cf_obsplot.pro @@ -0,0 +1,334 @@ +;+ +; cf_obsplot.pro reads the final, combined *00000all* spectral file +; for a given observation and produces quick-look plots for +; each of the four FUSE detectors. A fifth plot shows the +; best three channels that span the regions 900-1000, +; 1000-1100, and 1100-1200 A. +; +; Author: Ken Sembach +; Last Update: February 23, 2001 +; +; V. Dixon, 07/02/2001 - Modified for inclusion in pipeline. +; V. Dixon, 08/16/2001 - Read and write to current directory. +; V. Dixon, 08/31/2001 - Add a sixth plot spanning FUSE band. +; Mark region of worm in LiF 1B. +; V. Dixon, 11/05/2001 - Compute mean flux with 1-sig errors. +; V. Dixon, 11/20/2001 - Modify position of mean-flux plots. +; Add sentence describing error bars. +; V. Dixon, 01/17/2002 - Deal gracefully with missing data. +; V. Dixon, 03/12/2002 - Create .gif files if possible. +; B. Godard, 04/26/2004 - Print units and CalFUSE Ver +; Don't keep exposures with exptime=0 +; or quality array=0 +; Use quality array to compute min/max +; V. Dixon, 05/04/2004 - Scale ymax by 1.1 in all plots. +; Change output file names to match +; MAST's expectations. +; V. Dixon, 08/20/2004 - Change A to Angstroms at the bottom +; of each page. +; Exclude C III 977 when setting Y +; scale. +; V. Dixon, 04/14/2005 - Rewrite the program to use the new +; *all* files produced by CalFUSE v3.1. +; V. Dixon, 08/01/2005 - If IDL version > 6.0, generate GIF file. +; V. Dixon, 03/21/2006 - If channel has no data, set mflux = -1. +; Use "Exposure" or "Exposures" as +; appropriate. +; Consider all four channels for the +; 1000-1100 A plot. +; Pass xrange and yscale to TOP_PLOT +; and BOTTOM_PLOT. +; V. Dixon, 05/19/2006 - For BKGD targets, omit flux +; comparison when selecting which +; detector segments to use. +; Compute mean flux using same wavelength +; intervals for each band. +; V. Dixon, 05/24/2006 - For PC targets, don't include O VI in +; calculation of mean flux. +; V. Dixon, 12/12/2006 - If OBSTIME < 10 s for preferred segment +; and > 100 s for other segment, use other. +; V. Dixon, 12/19/2006 - Use the same sample regions in all plots. +; V. Dixon, 04/06/2007 - Comment out calls to set_plot,'X' and +; CLEANPLOT, as they can cause problems +; for OPUS. +; V. Dixon, 04/11/2008 - If obstime < 1000, quote in seconds. +; V. Dixon, 08/08/2008 - If airglow keyword is set, image files +; are named *00900* +; +; Calling sequence: OBSPLOT, obsname +; +; Input: obsname (e.g., P1030604) +; Output: obsname00000lif1ttagf.jpg, obsname00000lif2ttagf.jpg, +; obsname00000sic1ttagf.jpg, obsname00000sic2ttagf.jpg, +; obsname00000specttagf.jpg +; +; +;- + +pro TOP_PLOT, data, header, xrange, yscale, aper, sample, mflux, title + s=sample + !x.range=xrange + obstime = SXPAR(header,'OBSTIME') + extname = SXPAR(header,'EXTNAME') + num_exp = SXPAR(header,'NSPEC') + if (num_exp eq 1) then exposures=' Exposure, ' else exposures=' Exposures, ' + IF (obstime gt 0) THEN BEGIN + w = data.wave & f = data.flux + loc = WHERE(((w GT s[0]) AND (w LE s[1])) OR ((w GT s[2]) and (w LT s[3])) OR $ + ((w GT s[4]) and (w LT s[5]))) + ftmp = SMOOTH(f[loc],6) & nftmp = N_ELEMENTS(ftmp) + ymin = MIN(f[loc]) & ymax = MAX(ftmp) + ydiv = FIX(1-ALOG10(ymax)) + ymax = yscale * ymax *10.^ydiv + !y.range=[ymin,ymax] + mflux = mean(f[loc]) + + if (obstime lt 1000) then time_str = STRING(round(obstime),'(I3)')+' sec)' else $ + time_str = STRING(obstime/1000,'(F4.1)')+' ksec)' + + PLOT,w,f*10.^ydiv,nsum=12,title=title,ytitle='Flux (10!u-'+STRTRIM(ydiv,2)+'!n)',xtitle=' ' + XYOUTS,(!x.range[1]-!x.range[0])*0.5+!x.range[0],(!y.range[1]-!y.range[0])*0.88+!y.range[0], $ + extname+' '+aper+' ('+STRTRIM(num_exp,2)+exposures+time_str,charsize=1.3,align=0.5 + ENDIF ELSE BEGIN + !y.range=[0,1] + mflux = -1. + PLOT,[0,0],[0,0],title=title,ytitle='Flux',xtitle=' ' + XYOUTS,(!x.range[1]-!x.range[0])*0.5+!x.range[0],0.88, $ + extname+' '+aper+': No data found',charsize=1.5,align=0.5 + ENDELSE + RETURN +END + +pro BOTTOM_PLOT, data, header, xrange, yscale, aper, sample, mflux, vers + s=sample + !x.range=xrange + obstime = SXPAR(header,'OBSTIME') + extname = SXPAR(header,'EXTNAME') + num_exp = SXPAR(header,'NSPEC') + if (num_exp eq 1) then exposures=' Exposure, ' else exposures=' Exposures, ' + IF (obstime gt 0) THEN BEGIN + w = data.wave & f = data.flux + loc = WHERE(((w GT s[0]) AND (w LE s[1])) OR ((w GT s[2]) AND (w LT s[3])) OR $ + ((w GT s[4]) AND (w LT s[5]))) + ftmp = SMOOTH(f(loc),6) & nftmp = N_ELEMENTS(ftmp) + ymin = MIN(f[loc]) & ymax = MAX(ftmp) + ydiv = FIX(1-ALOG10(ymax)) + ymax = yscale * ymax *10.^ydiv + !y.range=[ymin,ymax] + mflux = mean(f[loc]) + + if (obstime lt 1000) then time_str = STRING(round(obstime),'(I3)')+' sec)' else $ + time_str = STRING(obstime/1000,'(F4.1)')+' ksec)' + + PLOT,w,f*10.^ydiv,nsum=12,title=' ',ytitle='Flux (10!u-'+STRTRIM(ydiv,2)+'!n)', $ + xtitle='Wavelength ('+string("305B)+')',/noeras + XYOUTS,(!x.range[1]-!x.range[0])*0.5+!x.range[0],(!y.range[1]-!y.range[0])*0.88+!y.range[0], $ + extname+' '+aper+' ('+STRTRIM(num_exp,2)+exposures+time_str, charsize=1.3,align=0.5 + IF (STRCMP(extname,'1BLIF',5, /FOLD_CASE) EQ 1) THEN $ + XYOUTS, 1160., (!y.range[1]-!y.range[0])*0.1+!y.range[0],'---- REGION OF WORM ----',align=0.5 + ENDIF ELSE BEGIN + !y.range=[0,1] + mflux = -1. + PLOT,[0,0],[0,0],title=' ',ytitle='Flux', xtitle='Wavelength ('+string("305B)+')',/noeras + XYOUTS,(!x.range[1]-!x.range[0])*0.5+!x.range[0],0.88, $ + extname+' '+aper+': No data found',charsize=1.5,align=0.5 + ENDELSE + XYOUTS, 0.05, 0.01, 'Flux units are erg cm!E-2!N s!E-1!N '+string("305B)+'!E-1!N.', alignment=0., /NORMAL + XYOUTS, 0.95, 0.01, 'CalFUSE v'+vers, alignment=1., /NORMAL + RETURN +END + +PRO CF_OBSPLOT, obsname, airglow=airglow + + ver = float(!version.release) + + if keyword_set(airglow) then fnumber = '00900' else fnumber = '00000' + list = FINDFILE(obsname+fnumber+'all*fcal.fit', COUNT=count) + lun = fxposit(list[0], 0, /SILENT) + foo = MRDFITS(lun,0,h,/SILENT) + + root = SXPAR(h,'PRGRM_ID')+SXPAR(h,'TARG_ID')+SXPAR(h,'SCOBS_ID') + root = STRCOMPRESS(root, /REMOVE_ALL) + targ = STRCOMPRESS(SXPAR(h,'TARGNAME'), /REMOVE_ALL) + objclass = SXPAR(h,'OBJCLASS') + sp_type = STRCOMPRESS(SXPAR(h,'SP_TYPE'), /REMOVE_ALL) + src_type = STRCOMPRESS(SXPAR(h,'SRC_TYPE'), /REMOVE_ALL) + mode = STRLOWCASE(STRCOMPRESS(SXPAR(h,'INSTMODE'), /REMOVE_ALL)) + aper = STRCOMPRESS(SXPAR(h,'APERTURE'), /REMOVE_ALL) + vers = STRCOMPRESS(SXPAR(h,'CF_VERS'), /REMOVE_ALL) + ext = '4' + IF aper EQ 'MDRS' THEN ext = '2' + IF aper EQ 'HIRS' THEN ext = '3' + if keyword_set(airglow) then title = ''+root+' [AIRGLOW]' else $ + title = ''+root+' ['+targ+']' + if (src_type eq 'PE') then yscale = 1.0 else yscale = 1.3 + if (objclass eq 7 or sp_type eq 'BKGD') then bkgd_obs = 1 else bkgd_obs = 0 + +; +; Read each channel and its header. +; + L1A = mrdfits(lun,0,HL1A) + L1B = mrdfits(lun,0,HL1B) + L2B = mrdfits(lun,0,HL2B) + L2A = mrdfits(lun,0,HL2A) + S1A = mrdfits(lun,0,HS1A) + S1B = mrdfits(lun,0,HS1B) + S2B = mrdfits(lun,0,HS2B) + S2A = mrdfits(lun,0,HS2A) + +; +; Use different sample regions for point and extended sources. +; + S1B_SAMPLE = [912, 935, 955, 970, 980, 985] + S2A_SAMPLE = [912, 935, 955, 970, 980, 985] + + if (src_type eq 'PC') then begin + L1A_SAMPLE = [0, 0, 0, 0, 1045, 1070] + S1A_SAMPLE = [0, 0, 0, 0, 1045, 1070] + L2B_SAMPLE = [0, 0, 0, 0, 1045, 1070] + S2B_SAMPLE = [0, 0, 0, 0, 1045, 1070] + endif else begin + L1A_SAMPLE = [0, 0, 1030, 1039, 1045, 1070] + S1A_SAMPLE = [0, 0, 1030, 1039, 1045, 1070] + L2B_SAMPLE = [0, 0, 1030, 1039, 1045, 1070] + S2B_SAMPLE = [0, 0, 1030, 1039, 1045, 1070] + endelse + + L1B_SAMPLE = [1095, 1130, 1140, 1165, 1170, 1190] + L2A_SAMPLE = [1095, 1130, 1140, 1165, 1170, 1190] + +; +; Plot individual channels. Skip this for airglow files. +; + ;set_plot,'X' + ;CLEANPLOT + set_plot,'Z' + device,set_resolution=[650,580] + + !x.style = 1 & !y.style = 1 + !x.ticklen=0.045 + !x.charsize=1.5 + !y.charsize=1.5 + !y.minor = 2 + !p.title = ' ' + !x.title = ' ' + !y.title = ' ' + !p.color = 0 + !p.background = 255 + +; +; LiF 1A and 1B +; + !p.position = [0.13,0.57,0.95,0.95] + TOP_PLOT, L1A, HL1A, [985,1085], yscale, aper, L1A_SAMPLE, FL1A, title + + !p.position = [0.13,0.1,0.95,0.48] + BOTTOM_PLOT, L1B, HL1B, [1092,1190], yscale, aper, L1B_SAMPLE, FL1B, vers + + if (ver ge 5.4 and ver le 6.0) then write_jpeg,root+fnumber+'lif1'+mode+'f.jpg',TVRD() $ + else write_gif,root+fnumber+'lif1'+mode+'f.gif',TVRD() + +; +; SiC 1A and 1B +; + !p.position = [0.13,0.57,0.95,0.95] + TOP_PLOT, S1A, HS1A, [1000,1094], yscale, aper, S1A_SAMPLE, FS1A, title + + !p.position = [0.13,0.1,0.95,0.48] + BOTTOM_PLOT, S1B, HS1B, [910,995], yscale, aper, S1B_SAMPLE, FS1B, vers + + if (ver ge 5.4 and ver le 6.0) then write_jpeg,root+fnumber+'sic1'+mode+'f.jpg',TVRD() $ + else write_gif,root+fnumber+'sic1'+mode+'f.gif',TVRD() + +; +; LiF 2A and 2B +; + !p.position = [0.13,0.57,0.95,0.95] + TOP_PLOT, L2A, HL2A, [1085,1185], yscale, aper, L2A_SAMPLE, FL2A, title + + !p.position = [0.13,0.1,0.95,0.48] + BOTTOM_PLOT, L2B, HL2B, [975,1078], yscale, aper, L2B_SAMPLE, FL2B, vers + + if (ver ge 5.4 and ver le 6.0) then write_jpeg,root+fnumber+'lif2'+mode+'f.jpg',TVRD() $ + else write_gif,root+fnumber+'lif2'+mode+'f.gif',TVRD() + +; +; SiC 2A and 2B +; + !p.position = [0.13,0.57,0.95,0.95] + TOP_PLOT, S2A, HS2A, [913,1009], yscale, aper, S2A_SAMPLE, FS2A, title + + !p.position = [0.13,0.1,0.95,0.48] + BOTTOM_PLOT, S2B, HS2B, [1013,1106], yscale, aper, S2B_SAMPLE, FS2B, vers + + if (ver ge 5.4 and ver le 6.0) then write_jpeg,root+fnumber+'sic2'+mode+'f.jpg',TVRD() $ + else write_gif,root+fnumber+'sic2'+mode+'f.gif',TVRD() + +; +; Now plot combined fluxes +; + ;set_plot,'X' + ;CLEANPLOT + set_plot,'Z' + device,set_resolution=[650,870] + + !x.style = 1 & !y.style = 1 + !x.ticklen=0.045 + !x.charsize=1.5 + !y.charsize=1.5 + !y.minor = 2 + !p.title = ' ' + !x.title = ' ' + !y.title = ' ' + !p.color = 0 + !p.background = 255 + +; +; 900 - 1000 A +; + !p.position = [0.13,0.70,0.95,0.95] + obs2a = SXPAR(HS2A, 'OBSTIME') + obs1b = SXPAR(HS1B, 'OBSTIME') + + if ((bkgd_obs eq 1 and FS2A gt -1) or (FS2A gt 0.9*FS1B) or (obs2a gt 100 and obs1b lt 10)) then begin + TOP_PLOT, S2A, HS2A, [913,1009], yscale, aper, S2A_SAMPLE, foo, title + endif else begin + TOP_PLOT, S1B, HS1B, [910,995], yscale, aper, S1B_SAMPLE, foo, title + endelse + +; +; 1000 - 1100 A +; + !p.noerase = 1 + !p.position = [0.13,0.40,0.95,0.65] + obl1a = SXPAR(HL1A, 'OBSTIME') + obl2b = SXPAR(HL2B, 'OBSTIME') + + if ((bkgd_obs eq 1 and FL1A gt -1) or $ + ((FL1A gt 0.9*FL2B and FL1A gt 0.7*FS1A and FL1A gt 0.7*FS2B) $ + and not (obl2b gt 100 and obl1a lt 10))) then $ + TOP_PLOT, L1A, HL1A, [985,1085], yscale, aper, L1A_SAMPLE, foo, '' else $ + if ((bkgd_obs eq 1 and FL2B gt -1) or $ + ((FL2B gt 0.7*FS1A and FL2B gt 0.7*FS2B) and not (obl2b lt 10 and obl1a gt 100))) then $ + TOP_PLOT, L2B, HL2B, [975,1078], yscale, aper, L2B_SAMPLE, foo, '' else $ + if ((bkgd_obs eq 1 and FS1A gt -1) or ((FS1A gt FS2B) and not (obl1a lt 10 and obl2b gt 100))) then $ + TOP_PLOT, S1A, HS1A, [1000,1094], yscale, aper, S1A_SAMPLE, foo, '' else $ + TOP_PLOT, S2B, HS2B, [1013,1106], yscale, aper, S2B_SAMPLE, foo, '' + +; +; 1100 - 1200 A +; + !p.position = [0.13,0.1,0.95,0.35] + obl2a = SXPAR(HL2A, 'OBSTIME') + obl1b = SXPAR(HL1B, 'OBSTIME') + + if ((bkgd_obs eq 1 and FL2A gt -1) or (FL2A gt 0.9*FL1B) or (obl2a gt 100 and obl1b lt 10)) then $ + BOTTOM_PLOT, L2A, HL2A, [1085,1182], yscale, aper, L2A_SAMPLE, foo, vers else $ + BOTTOM_PLOT, L1B, HL1B, [1092,1190], yscale, aper, L1B_SAMPLE, foo, vers + + if (ver ge 5.4 and ver le 6.0) then write_jpeg,root+fnumber+'spec'+mode+'f.jpg',TVRD() $ + else write_gif,root+fnumber+'spec'+mode+'f.gif',TVRD() + + RETURN +END + diff --git a/idl/cf_plot_extract3.pro b/idl/cf_plot_extract3.pro new file mode 100644 index 0000000..c201754 --- /dev/null +++ b/idl/cf_plot_extract3.pro @@ -0,0 +1,210 @@ +;+ +; cf_plot_extract3.pro is a procedure to plot the extraction windows over a +; greyscale image for a given IDF to a standard image file. +; +; Author: Edward M. Murphy +; Written: 1999 November 24 +; +; Calling sequence: cf_plot_extract3, rootname +; +; Inputs: rootname will have 'idf.fit' added to it. +; +; 02/12/03 jch: Adapt from cf_plot_extract.pro to work for CalFUSE 3.0. +; 03/10/03 wvd: Silence calls to mrdfits. Simplify image construction. +; 04/08/03 wvd: Shift apertures in X to match FPA positions. +; 06/02/03 wvd: Don't plot non-target apertures in HIST mode. +; 10/10/03 wvd: Exclude photons from bad times and those outside PHA limits. +; 12/09/03 wvd: Use extended apertures when appropriate. +; 02/26/04 wvd: Use strcmp to compare strings. +; 05/03/04 wvd: In HIST mode, don't plot non-target apertures. +; Read expected spectral centroid from CHID_CAL file header. +; 05/14/04 wvd: Don't use strcmp to compare strings. +; Modify boundaries between LiF and SiC channels. +; 06/14/04 wvd: Don't shift apertures in X to match FPA positions. +; 07/12/05 wvd: Clean up call to HIST_2D. +; 08/01/05 wvd: If IDL version > 6.0, generate GIF file. +; 08/30/05 wvd: If there are no good photons, plot all photons. +; 11/01/05 wvd: If DAYNIGHT = NIGHT, exclude daytime photons. +; 03/14/08 wvd: Require LOC_FLAG < 16B for good data. +;- + +pro cf_plot_extract3,rootname + + !quiet=1 + filename=rootname+'idf.fit' + a=mrdfits(filename,0,ah,/SILENT) + + ; Read input file's primary header keywords + wave_cal = strcompress(fxpar(ah,'WAVE_CAL'),/rem) + chid_cal = strcompress(fxpar(ah,'CHID_CAL'),/rem) + filedate = strcompress(fxpar(ah,'DATE'),/rem) + obsdate = strcompress(fxpar(ah,'DATEOBS'),/rem) + obstime = strcompress(fxpar(ah,'TIMEOBS'),/rem) + cfvers = strcompress(fxpar(ah,'CF_VERS'),/rem) + exptime = fxpar(ah,'EXPTIME') + instmode = strcompress(fxpar(ah,'INSTMODE'),/rem) + detector = strcompress(fxpar(ah,'DETECTOR'),/rem) + aperture = strcompress(fxpar(ah,'APERTURE'),/rem) + srctype = strmid(fxpar(ah,'SRC_TYPE'),0,1) + daynight = strcompress(fxpar(ah,'DAYNIGHT'),/rem) + + ; Interpret aperture and srctype keywords. + if (aperture EQ 'MDRS') then begin + aplif = 2 + apsic = 6 + endif else if (aperture EQ 'HIRS') then begin + aplif = 1 + apsic = 5 + endif else begin ; assume LWRS + aplif = 3 + apsic = 7 + endelse + if (srctype EQ 'P') then extend = 0 else extend = 8 + + cal_path = getenv("CF_CALDIR") + slitname=["HIRS","MDRS","LWRS","PINH","HIRS","MDRS","LWRS","PINH"] + specname=["LiF","LiF","LiF","LiF","SiC","SiC","SiC","SiC"] + + ; Set the plot resolution. + set_plot,'Z' + device,set_resolution=[1124,612] + + ; Define plot attributes + !x.style=1 + !y.style=1 + !x.tickformat='(I5)' + !p.title=filename + !x.title="X pixel" + !y.title="Y pixel" + !p.charsize=1.2 + cs=0.8 + + ; Read the X and Y columns from the first extension (photon list), + ; exclude bad events, and generate an image from the data. + ftab_ext, filename, 'X,Y,TIMEFLGS,LOC_FLGS', xarr, yarr, timeflag, loc_flag + if (daynight EQ 'NIGHT') then begin + good = where (timeflag eq 0 and loc_flag lt 16, n) + endif else begin + good = where ((timeflag and not 1B) eq 0 and loc_flag lt 16, n) + endelse + if (n gt 0) then begin + xarr = xarr[good] + yarr = yarr[good] + endif + a = HIST_2D(xarr,yarr,bin1=16,bin2=2,min1=0,max1=16383,min2=0,max2=1023) + b = HIST_EQUAL(a) + + !p.charsize=1.3 + x0=85 + y0=65 + x1=1024+x0 + y1=512+y0 + dummy_x = indgen(2)*16383 + dummy_y = indgen(2)*1023 + + ; Plot the pixel axes and the photon image. + plot,dummy_x,dummy_y,/nodata,position=[x0,y0,x1,y1],/DEVICE,$ + background=255,color=0,xticks=8 + tv,255-b,x0,y0 + plot,dummy_x,dummy_y,/nodata,/noerase,position=[x0,y0,x1,y1],/DEVICE,$ + background=255,color=0,xticks=8,xminor=4 + + miny=10000 + maxy=0.0 + + ; Plot the extraction windows. + for i=1,7 do begin + ; Skip the pinhole and any non-target apertures for HIST data. + if ((i ne 4) and ((instmode EQ 'TTAG') $ + or (i eq aplif) or (i eq apsic))) then begin + + ; Read aperture limits from calibration file. + j = 8; use extended apertures for non-target apertures + if ((i eq aplif) or (i eq apsic)) then j = extend + xw=mrdfits(cal_path+'/'+chid_cal,i+j,chidhdr,/SILENT) + yhigh = xw.yhigh + ylow = xw.ylow + default_y_centroid = fxpar(chidhdr,'CENTROID') + + ; Get the YCENT keywords and adjust the windows accordingly. + key = 'YCENT' + strcompress(string(i),/rem) + ycent = fxpar(ah, key) + shift = nint(default_y_centroid - ycent) + yhigh = yhigh - shift + ylow = ylow - shift + + ; Plot the extraction windows. + oplot,yhigh,color=0 + oplot,ylow,color=0 + miny=min([miny,ylow]) + maxy=max([maxy,yhigh]) + xyouts,15600,ylow[16350]+5,specname[i-1],$ + color=0,charsize=0.8 + xyouts,15600,ylow[16350]-20,slitname[i-1],$ + color=0,charsize=0.8 + endif + endfor + + ; Write history information. + xyouts,50,20,'Observation date='+obsdate+'T'+obstime,$ + charsize=0.90,color=0,/DEVICE + xyouts,50,5,'IDF file date= '+filedate,charsize=0.90,color=0,/DEVICE + xyouts,1124/2,5,'CALFUSE version '+cfvers,$ + charsize=0.90,color=0,/DEVICE,align=0.5 + xyouts,1124-50,20,'CHID calibration file='+chid_cal,$ + charsize=0.90,color=0,/DEVICE,align=1.0 + xyouts,1124-50,5,'Exposure time='+$ + strcompress(string(exptime,format="(F10.2)"),/rem),$ + charsize=0.90,color=0,/DEVICE,align=1.0 + + ; Read in and plot the wavelength scale + wl_lif=mrdfits(cal_path+'/'+wave_cal,4,ah,/SILENT) + wl_sic=mrdfits(cal_path+'/'+wave_cal,8,ah,/SILENT) + + wlsc_lif=[wl_lif[0].wavelength,wl_lif[16383].wavelength] + wlsc_sic=[wl_sic[0].wavelength,wl_sic[16383].wavelength] + + detector=strtrim(detector,2) + + ; Determine locations of the wavelength axes. + if (detector EQ '1A') then begin + hip=maxy+25 + midp=425 + lowp=miny-15 + endif else if (detector EQ '1B') then begin + hip=maxy+25 + midp=425 + lowp=miny-10 + endif else if (detector EQ '2A') then begin + hip=maxy+100 + midp=500 + lowp=miny-100 + endif else if (detector EQ '2B') then begin + hip=maxy+100 + midp=500 + lowp=miny-100 + endif else begin + hip=maxy+25 + midp=512 + lowp=miny-15 + endelse + + ; Draw wavelength axes. + axis,0,hip,xaxis=1,xrange=wlsc_lif,xstyle=1,charsize=cs,$ + color=0,xtitle="Wavelength" + axis,0,midp,xaxis=0,xrange=wlsc_lif,xstyle=1,charsize=0.001,$ + color=0,xtitle="Wavelength" + + axis,0,midp,xaxis=1,xrange=wlsc_sic,xs=1,charsize=0.001,$ + color=0,xtitle="" + axis,0,lowp,xaxis=0,xrange=wlsc_sic,xs=1,charsize=cs,$ + color=0,xtitle="Wavelength" + + ; Create a GIF/JPEG file: + ver = float(!version.release) + if (ver ge 5.4 and ver le 6.0) then write_jpeg,rootname+'ext.jpg',TVRD() $ + else write_gif,rootname+'ext.gif',TVRD() + +EXIT: +return +end diff --git a/idl/cf_plot_rate3.pro b/idl/cf_plot_rate3.pro new file mode 100644 index 0000000..6a2669b --- /dev/null +++ b/idl/cf_plot_rate3.pro @@ -0,0 +1,227 @@ +;+ +; cf_plot_rate3.pro is a procedure to plot the count rate and screening +; as a function of time for calfuse 3.0 data. +; +; Author: Edward M. Murphy +; Written: 2000 April 21 +; +; Calling sequence: cf_plot_rate3, rootname +; +; Inputs: rootname will have 'idf.fit' added to it. +; +; 02/20/03 JC Hsu Adapt from cf_plot_rate.pro to work on version 3.0 data. +; 06/23/04 V Dixon Mark low-voltage periods with cross-hatch. +; 03/16/05 V Dixon Ignore timeline table entries that extend beyond +; the end of the exposure. +; 08/01/05 V Dixon If IDL version > 6.0, generate GIF file. +; 11/03/05 V Dixon Shade OPUS and SAA bad times before drawing rest of plot. +;- + +pro cf_plot_rate3,rootname + + !quiet=1 + + ; Open input IDF file + filename=rootname+'idf.fit' + dummy=mrdfits(filename,0,ah,/SILENT) + + ; Read primary header keywords + pdate = fxpar(ah,'DATE') + dateobs = fxpar(ah,'DATEOBS') + expstart = fxpar(ah,'EXPSTART') + expend = fxpar(ah,'EXPEND') + detector = fxpar(ah,'DETECTOR') + + exposure = (expend - expstart) * 24. * 3600. + + ; Define plot attributes. + !p.multi=[0,1,2] + set_plot,'Z' + device,set_resolution=[1024,512] + + !x.style=1 + !y.style=1 + !x.tickformat='(I6)' + + !x.title="Time (seconds)" + !y.title="Count rate" + !p.charsize=1.2 + !p.background=255 + !p.color=0 + + ; Read columns from the timeline extension (ext 3) + ftab_ext, filename, 'TIME', time0, exten_no=3 + ftab_ext, filename, 'STATUS_FLAGS', flags, exten_no=3 + if (exposure gt 55000) then begin + good = where(time0 lt 55000, n) + exposure = time0[n-1] - time0[0] + endif + good = where(time0 lt exposure) + time = [0, time0[good], exposure] + + ; Compute the day/night screenings (bit 1) + flags1 = ishft(flags, 7) + flags1 = ishft(flags1, -7) ; only keep bit 1 + + ; Compute the limb angle screening (bit 2) + flags2 = ishft(flags, 6) ; + flags2 = ishft(flags2, -7) ; only keep bit 2 + + ; Compute the SAA screening (bit 3) + flags3 = ishft(flags, 5) ; + flags3 = ishft(flags3, -7) ; only keep bit 3 + + ; Compute the HV screening (bit 4) + flags4 = ishft(flags, 4) ; + flags4 = ishft(flags4, -7) ; only keep bit 4 + + ; Compute burst screening (bit 5) + flags5 = ishft(flags, 3) ; + flags5 = ishft(flags5, -7) ; only keep bit 5 + + ; Compute OPUS screening (bit 6) + flags6 = ishft(flags, 2) ; + flags6 = ishft(flags6, -7) ; only keep bit 6 + + ; Compute jitter screening (bit 7) + flags7 = ishft(flags, 1) ; + flags7 = ishft(flags7, -7) ; only keep bit 7 + + ftab_ext, filename, 'BKGD_CNT_RATE', bkgd, exten_no=3 + len = size(time0) + packet = 1 > len[1]/100 ; determine rebinning size + wht = where(time0 eq time0) + whbin = where(wht mod packet eq 0) + + for i = 1, 2 do begin + + if (i eq 1) then text = ' LiF' else text = ' SiC' + !p.title=filename + text + + ; Draw grid, shade bad times, then over-plot count rate arrays. + if (i eq 1) then col = 'LIF' else col = 'SIC' + ftab_ext, filename, col+'_CNT_RATE', crate, exten_no=3 + ymin = 1 + ymax = max([1000,crate * 2.]) + crate = crate > ymin + !y.range=[ymin, ymax] + plot_io, time0, crate, /nodata, xr=[0,exposure] + + ; Shade OPUS + flags16 = flags6 * (ymax-ymin) + ymin + flags16 = [0, flags16[good], 0] + polyfill, time, flags16, COL = 0.75 * !D.N_COLORS + ;oplot, time, flags16, psym=10, linestyle=2 + + ; Shade SAA violations + flags13 = flags3 * (ymax/100-ymin) + ymin + flags13 = [0., flags13[good], 0.] + polyfill, time, flags13, COL = 0.25 * !D.N_COLORS + flags13 = flags3 * (ymax-ymin) + ymin + oplot, time, flags13, psym=10, linestyle=2 + + ; Shade day/night screenings + flags11 = flags1 * (0.2*ymax) + ymax + flags11 = [ymax, flags11[good], ymax] + oplot, time, flags11, psym=10 + polyfill, time, flags11 + + ; Shade limb angle + flags12 = flags2 * (ymax-ymin) + ymin + flags12 = [0, flags12[good], 0] + oplot, time, flags12, psym=10, linestyle=2 + polyfill, time, flags12, /line_fill, orientation=0 + + ; Shade low-voltage periods. + flags14 = flags4 * (ymax-ymin) + ymin + flags14 = [0, flags14[good], 0] + oplot, time, flags14, psym=10, linestyle=2 + polyfill, time, flags14, /line_fill, orientation=45 + polyfill, time, flags14, /line_fill, orientation=135 + + ; Shade burst + flags15 = flags5 * (ymax-ymin) + ymin + flags15 = [0, flags15[good], 0] + oplot, time, flags15, psym=10, linestyle=2 + polyfill, time, flags15, /line_fill, orientation=45 + + ; Shade jitter + flags17 = flags7 * (ymax-ymin) + ymin + flags17 = [0, flags17[good], 0] + oplot, time, flags17, psym=10, linestyle=2 + polyfill, time, flags17, /line_fill, orientation=135 + + ; Overplot count-rate arrays. + oplot, time0, crate, linestyle=0, psym=10 + oplot, time0(whbin), bkgd(whbin)+0.001, linestyle=1 + endfor + + ; Draw legend + plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[345-20,5,345+21,25],$ + /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$ + yrange=[0,1],/DEVICE + polyfill,[0,1,1,0],[0,0,1,1], COL = 0.75 * !D.N_COLORS + oplot,[0,1,1,0,0],[0,0,1,1,0],psym=10 + xyouts,345-17,10,'OPUS',charsize=1,color=0,/DEVICE + + plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[400-20,5,400+21,25],$ + /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$ + yrange=[0,1],/DEVICE + polyfill,[0,1,1,0],[0,0,0.36,0.36], COL = 0.25 * !D.N_COLORS + oplot,[0,1,1,0,0],[0,0,1,1,0],psym=10 + xyouts,400-12,13,'SAA',charsize=1,color=0,/DEVICE + + plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[450-15,5,450+15,25],$ + /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$ + yrange=[0,1],/DEVICE + polyfill,[0,1,1,0],[0.8,0.8,1,1] + xyouts,450-11,10,'Day',charsize=1,color=0,/DEVICE + + plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[710-15,5,710+15,25],$ + /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$ + yrange=[0,1],/DEVICE + polyfill,[0,1,1,0],[0,0,1,1],/line_fill,orientation=0 + oplot,[0,1,1,0,0],[0,0,1,1,0],psym=10 + xyouts,710+20,10,'Limb',charsize=1,color=0,/DEVICE + + plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[790-15,5,790+15,25],$ + /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$ + yrange=[0,1],/DEVICE + polyfill,[0,1,1,0],[0,0,1,1],/line_fill,orientation=45 + oplot,[0,1,1,0,0],[0,0,1,1,0],psym=10 + xyouts,790+20,10,'Burst',charsize=1,color=0,/DEVICE + + plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[870-15,5,870+15,25],$ + /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$ + yrange=[0,1],/DEVICE + polyfill,[0,1,1,0],[0,0,1,1],/line_fill,orientation=135 + oplot,[0,1,1,0,0],[0,0,1,1,0],psym=10 + xyouts,870+20,10,'Jitter',charsize=1,color=0,/DEVICE + + plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[950-15,5,950+15,25],$ + /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$ + yrange=[0,1],/DEVICE + polyfill,[0,1,1,0],[0,0,1,1],/line_fill,orientation=45 + polyfill,[0,1,1,0],[0,0,1,1],/line_fill,orientation=135 + oplot,[0,1,1,0,0],[0,0,1,1,0],psym=10 + xyouts,950+20,10,'HV',charsize=1,color=0,/DEVICE + + xyouts,835,27,'Screenings :',charsize=.9,color=0,/DEVICE,alignment=.5 + xyouts,20,16,'RAW file written at '+pdate,charsize=0.90,color=0,$ + /DEVICE + xyouts,925,495,'DATEOBS: '+dateobs,charsize=.9,color=0,$ + /DEVICE,alignment=.5 + + ; Create a GIF/JPEG file: + ver = float(!version.release) + if (ver ge 5.4 and ver le 6.0) then begin + xyouts,20,6,'JPG file written at '+!stime,charsize=0.90,$ + color=0,/DEVICE + write_jpeg,rootname+'rat.jpg',TVRD() + endif else begin + xyouts,20,6,'GIF file written at '+!stime,charsize=0.90,$ + color=0,/DEVICE + write_gif,rootname+'rat.gif',TVRD() + endelse +return +end diff --git a/idl/cf_plot_rate_linear.pro b/idl/cf_plot_rate_linear.pro new file mode 100644 index 0000000..e2e3aa7 --- /dev/null +++ b/idl/cf_plot_rate_linear.pro @@ -0,0 +1,175 @@ +;+ +; cf_plot_rate_linear.pro is a procedure to plot the count rate and screening +; as a function of time for calfuse 3.0 data. +; THIS VERSION PRODUCES A LINEAR PLOT OF COUNT RATE +; +; Author: Edward M. Murphy +; Written: 2000 April 21 +; +; Calling sequence: cf_plot_rate3, rootname +; +; Inputs: rootname will have 'idf.fit' added to it. +; +; 02/20/03 JC Hsu Adapt from cf_plot_rate.pro to work on version 3.0 data. +; 05/17/04 V Dixon Adapt from cf_plot_rate3.pro +;- + +pro cf_plot_rate_linear, rootname + + !quiet=1 + + ; Open input IDF file + filename=rootname+'idf.fit' + dummy=mrdfits(filename,0,ah,/SILENT) + + ; Read primary header keywords + pdate = fxpar(ah,'DATE') + dateobs = fxpar(ah,'DATEOBS') + expstart = fxpar(ah,'EXPSTART') + expend = fxpar(ah,'EXPEND') + detector = fxpar(ah,'DETECTOR') + + exposure = (expend - expstart) * 24. * 3600. + + ; Define plot attributes. + !p.multi=[0,1,2] + set_plot,'Z' + device,set_resolution=[1024,512] + + !x.style=1 + !y.style=1 + !x.tickformat='(I6)' + + !x.title="Time (seconds)" + !y.title="Count rate" + !p.charsize=1.2 + !p.background=255 + !p.color=0 + + ; Read columns from the timeline extension (ext 3) + ftab_ext, filename, 'TIME', time0, exten_no=3 + ftab_ext, filename, 'STATUS_FLAGS', flags, exten_no=3 + time = [0, time0, exposure] + + ; Compute the day/night screenings (bit 1) + flags1 = ishft(flags, 7) + flags1 = ishft(flags1, -7) ; only keep bit 1 + + ; Compute the limb angle screening (bit 2) + flags2 = ishft(flags, 6) ; + flags2 = ishft(flags2, -7) ; only keep bit 2 + + ; Compute the SAA screening (bit 3) + flags3 = ishft(flags, 5) ; + flags3 = ishft(flags3, -7) ; only keep bit 3 + + ; Compute burst screening (bit 5) + flags5 = ishft(flags, 3) ; + flags5 = ishft(flags5, -7) ; only keep bit 5 + + ; Compute jitter screening (bit 7) + flags7 = ishft(flags, 1) ; + flags7 = ishft(flags7, -7) ; only keep bit 7 + + ftab_ext, filename, 'BKGD_CNT_RATE', bkgd, exten_no=3 + len = size(time0) + packet = 1 > len[1]/100 ; determine rebinning size + wht = where(time0 eq time0) + whbin = where(wht mod packet eq 0) + + for i = 1, 2 do begin + + if (i eq 1) then text = ' LiF' else text = ' SiC' + !p.title=filename + text + + ; Plot the count rate + if (i eq 1) then col = 'LIF' else col = 'SIC' + ftab_ext, filename, col+'_CNT_RATE', crate, exten_no=3 + ymax = max(crate)*1.2 + ymin = 0. + !y.range=[ymin, ymax] + plot, time0, crate, xr=[0,exposure], linestyle=0, psym=10 + oplot, time0(whbin), bkgd(whbin)*(10./i), linestyle=1 + if (i eq 1) then text = '5' else text = '10' + xyouts,exposure/50,ymax/20,'BKGD scaled by '+text,charsize=1 + + + ; Shade day/night screenings + ;;;flags11 = flags1 * (-0.2*ymin) + ymin + ;;;flags11 = [0, flags11, 0] + flags11 = flags1 * (0.025*ymax) + ymax + flags11 = [ymax, flags11, ymax] + oplot, time, flags11, psym=10 + polyfill, time, flags11 + + ; Shade limb angle + flags12 = flags2 * (ymax-ymin) + ymin + flags12 = [ymin, flags12, ymin] + oplot, time, flags12, psym=10, linestyle=2 + polyfill, time, flags12, /line_fill, orientation=0 + + ; Shade SAA + ;;;flags13 = flags3 * (ymax-ymin) + ymin + ;;;flags13 = [ymin, flags13, ymin] + ;;;oplot, time, flags13, psym=10, linestyle=2 + + ; Shade burst + flags15 = flags5 * (ymax-ymin) + ymin + flags15 = [ymin, flags15, ymin] + oplot, time, flags15, psym=10, linestyle=2 + polyfill, time, flags15, /line_fill, orientation=45 + + ; Shade jitter + flags17 = flags7 * (ymax-ymin) + ymin + flags17 = [ymin, flags17, ymin] + oplot, time, flags17, psym=10, linestyle=2 + polyfill, time, flags17, /line_fill, orientation=135 + endfor + + ; Draw legend + plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[710-15,5,710+15,25],$ + /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$ + yrange=[0,1],/DEVICE + polyfill,[0,1,1,0],[0.8,0.8,1,1] + xyouts,710+20,10,'Day',charsize=1,color=0,/DEVICE + + plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[790-15,5,790+15,25],$ + /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$ + yrange=[0,1],/DEVICE + polyfill,[0,1,1,0],[0,0,1,1],/line_fill,orientation=0 + oplot,[0,1,1,0,0],[0,0,1,1,0],psym=10 + xyouts,790+20,10,'Limb',charsize=1,color=0,/DEVICE + + plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[870-15,5,870+15,25],$ + /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$ + yrange=[0,1],/DEVICE + polyfill,[0,1,1,0],[0,0,1,1],/line_fill,orientation=45 + oplot,[0,1,1,0,0],[0,0,1,1,0],psym=10 + xyouts,870+20,10,'Burst',charsize=1,color=0,/DEVICE + + plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[950-15,5,950+15,25],$ + /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$ + yrange=[0,1],/DEVICE + polyfill,[0,1,1,0],[0,0,1,1],/line_fill,orientation=135 + oplot,[0,1,1,0,0],[0,0,1,1,0],psym=10 + xyouts,950+20,10,'Jitter',charsize=1,color=0,/DEVICE + + xyouts,835,27,'Screenings :',charsize=.9,color=0,/DEVICE,alignment=.5 + xyouts,20,16,'RAW file written at '+pdate,charsize=0.90,color=0,$ + /DEVICE + xyouts,925,495,'DATEOBS: '+dateobs,charsize=.9,color=0,$ + /DEVICE,alignment=.5 + + ; Create a GIF/JPEG file: + ver = float(!version.release) + if (ver ge 5.4) then begin + xyouts,20,6,'JPG file written at '+!stime,charsize=0.90,$ + color=0,/DEVICE + write_jpeg,rootname+'lin.jpg',TVRD() + endif else begin + xyouts,20,6,'GIF file written at '+!stime,charsize=0.90,$ + color=0,/DEVICE + write_gif,rootname+'lin.gif',TVRD() + endelse +return +end diff --git a/idl/fuse_extract.pro b/idl/fuse_extract.pro new file mode 100644 index 0000000..0b784d7 --- /dev/null +++ b/idl/fuse_extract.pro @@ -0,0 +1,154 @@ +pro fuse_extract, input +;+ +; NAME: +; FUSE_EXTRACT +; +;*PURPOSE: +; To extract full-resolution spectra from the primary LiF and SiC +; apertures of an intermediate data file. +; +;*CATEGORY: +; INPUT/OUTPUT +; +;*CALLING SEQUENCE: +; FUSE_EXTRACT, input +; +;*INPUT: +; input: intermediate data file +; +;*OUTPUT: +; LiF and SiC spectral files in v2.x format. +; +; +;*PROCEDURES USED: +; FXADDPAR, FXHMAKE, FXBADDCOL, FXBCREATE, FXBWRITE, FXBFINISH, MRDFITS, +; SXPAR +; +;*HISTORY: +; Written by: V. Dixon March 2003 +; 03/23/2003 V. Dixon Correct X values for Doppler shift. +; 03/24/2003 V. Dixon Fix bug in call to where(). +; 06/27/2003 R. Robinson Update to new version of IDF +; - added fscale keyword to mrdfits call +; - divide flux by exposure time +; 06/08/2004 V. Dixon Use EXPTIME rather than RAWTIME. +; Correct calculation of COUNTS +; 07/01/2004 V. Dixon Don't read headers of extensions. +; Employ faster scheme for calculating +; Doppler shift. +; 07/02/2004 V. Dixon Don't compute dl per pixel. It suffers +; from truncation errors. Instead, use +; value from file header. +;- +;----------------------------------------------------------------------------- +; +; print calling sequence if no parameters supplied +; + if n_params(0) lt 1 then begin + print,'CALLING SEQUENCE: fuse_extract, input' + return + end +; +; Physical constants +; + C = 2.99792458e5 ; km/s +; +; Open input file and read keywords from header. +; + dummy = mrdfits(input, 0, hdr, status=status) + rootname = SXPAR(hdr, 'ROOTNAME') + vhelio = SXPAR(hdr, 'V_HELIO') + detector = SXPAR(hdr, 'DETECTOR') + aperture = SXPAR(hdr, 'APERTURE') + wavecal = SXPAR(hdr, 'WAVE_CAL') + exptime = SXPAR(hdr, 'EXPTIME') + + ap_array = ['HIRS', 'MDRS', 'LWRS', 'PINH'] + ap_code = ['3','2','4','1'] + for k = 0, 3 do if (strcmp(aperture, ap_array[k], 1)) then ap = k + channel = ap+1 + + GET_DATE, FITS_date +; +; Read photon information in the first extension of the IDF. +; + a = mrdfits(input, 1, status=status, /fscale) +; +; Read timeline table from the third extension of the IDF. +; + b = mrdfits(input, 3, status=status, /fscale) + nseconds = n_elements(b.time) + if (nint(b.time[nseconds-1] - b.time[0]) gt nseconds) then $ + print, 'WARNING: '+input+' timeline table not continuous.' + + for k = channel, channel+4, 4 do begin + print, 'Analyzing extension ', k + ; + ; Read corresponding wavelength array from WAVE_CAL file. + ; + wcalfile=strtrim('/data1/fuse/calfuse/v3.0/calfiles/'+wavecal,2) + print, 'Reading wavelength data from ', wcalfile + w=mrdfits(wcalfile,k,w_hdr) + wave = w.wavelength + wpc = SXPAR(w_hdr, 'DISPAPIX') + ; + ; Select photons from aperture of interest. + ; + g = where (a.channel eq k and (a.timeflgs and not 1B) eq 0 and $ + (a.loc_flgs and 16B) eq 0, n) + print, 'Number of photons in spectrum = ', n + x = a.x[g] + ; + ; Correct X values for Doppler shift. + ; + print, 'Applying Doppler shift' + t = nint(a.time[g] - b.time[0]) + ovc = (b.ORBITAL_VEL + vhelio)/C + for i = 0L, n-1L do x[i] = x[i] + ovc[t[i]] * a.lambda[g[i]]/wpc + ; + ; Extract spectrum and populate output arrays. + ; + print, 'Extracting spectrum' + x = nint(x) + flux = fltarr(16384) + error = fltarr(16384) + cntserr = fltarr(16384) + for i = 0L, n-1L do flux[x[i]] = flux[x[i]] + a.ergcm2[g[i]] + flux = flux / exptime / abs(wpc) + counts = histogram(x, min=0, max=16383) + g = where (counts gt 0) + cntserr[g] = sqrt(counts[g]) + error[g] = flux[g]/cntserr[g] + quality = indgen(16384) + ; + ; Create primary HDU of output file + ; + if (k lt 5) then ch = 'lif' else ch = 'sic' + filename = strcompress(rootname+strlowcase(detector)+ch+ap_code[ap]+'ttagfcal.fit',/remove_all) + FXADDPAR, hdr, 'DATE', FITS_date + FXADDPAR, hdr, 'FILENAME', filename + FXADDPAR, hdr, 'FILETYPE', 'CALIBRATED EXTRACTED SPECTRUM' + FXADDPAR, hdr, 'APER_ACT', strupcase(strcompress(ap_array[ap]+'_'+ch,/remove_all)) + FXHMAKE, hdr, /EXTEND + FXWRITE, filename, hdr + + FXBHMAKE, HEADER, 1, 'FUSE 1D Spectrum' + FXBADDCOL, iwave, HEADER, wave, 'WAVE', TUNIT='ANGSTROMS' + FXBADDCOL, iflux, HEADER, flux, 'FLUX', TUNIT='ERG/CM2/S/A' + FXBADDCOL, ierror, HEADER, error, 'ERROR', TUNIT='ERG/CM2/S/A' + FXBADDCOL, iquality, HEADER, quality, 'QUALITY', TUNIT='UNITLESS' + FXBADDCOL, icounts, HEADER, counts, 'COUNTS', TUNIT='COUNTS' + FXBADDCOL, icntserr, HEADER, cntserr, 'CNTSERR', TUNIT='COUNTS' + + FXBCREATE, out, filename, HEADER + FXBWRITE, out, wave, iwave, 1 + FXBWRITE, out, flux, iflux, 1 + FXBWRITE, out, error, ierror, 1 + FXBWRITE, out, quality, iquality, 1 + FXBWRITE, out, counts, icounts, 1 + FXBWRITE, out, cntserr, icntserr, 1 + FXBFINISH, out + endfor + close, /all + return + end diff --git a/idl/fuse_sic_only.pro b/idl/fuse_sic_only.pro new file mode 100644 index 0000000..f0eadc0 --- /dev/null +++ b/idl/fuse_sic_only.pro @@ -0,0 +1,159 @@ +pro fuse_sic_only, input +;+ +; NAME: +; FUSE_EXTRACT +; +;*PURPOSE: +; To extract full-resolution spectra from the primary LiF and SiC +; apertures of an intermediate data file. +; +; THIS VERSION EXTRACTS ONLY A SIC SPECTRUM. +; +;*CATEGORY: +; INPUT/OUTPUT +; +;*CALLING SEQUENCE: +; FUSE_EXTRACT, input +; +;*INPUT: +; input: intermediate data file +; +;*OUTPUT: +; LiF and SiC spectral files in v2.x format. +; +; +;*PROCEDURES USED: +; FXADDPAR, FXHMAKE, FXBADDCOL, FXBCREATE, FXBWRITE, FXBFINISH, MRDFITS, +; SXPAR +; +;*HISTORY: +; Written by: V. Dixon March 2003 +; 03/23/2003 V. Dixon Correct X values for Doppler shift. +; 03/24/2003 V. Dixon Fix bug in call to where(). +; 06/27/2003 R. Robinson Update to new version of IDF +; - added fscale keyword to mrdfits call +; - divide flux by exposure time +; 06/08/2004 V. Dixon Use EXPTIME rather than RAWTIME. +; Correct calculation of COUNTS +; 07/01/2004 V. Dixon Don't read headers of extensions. +; Employ faster scheme for calculating +; Doppler shift. +; 07/02/2004 V. Dixon Don't compute dl per pixel. It suffers +; from truncation errors. Instead, use +; value from file header. +; 07/15/2004 V. Dixon Extract only SiC spectrum. +;- +;----------------------------------------------------------------------------- +; +; print calling sequence if no parameters supplied +; + if n_params(0) lt 1 then begin + print,'CALLING SEQUENCE: fuse_sic_only, input' + return + end +; +; Physical constants +; + C = 2.99792458e5 ; km/s +; +; Open input file and read keywords from header. +; + dummy = mrdfits(input, 0, hdr, status=status) + rootname = SXPAR(hdr, 'ROOTNAME') + vhelio = SXPAR(hdr, 'V_HELIO') + detector = SXPAR(hdr, 'DETECTOR') + aperture = SXPAR(hdr, 'APERTURE') + wavecal = SXPAR(hdr, 'WAVE_CAL') + exptime = SXPAR(hdr, 'EXPTIME') + + ap_array = ['HIRS', 'MDRS', 'LWRS', 'PINH'] + ap_code = ['3','2','4','1'] + for k = 0, 3 do if (strcmp(aperture, ap_array[k], 1)) then ap = k + channel = ap+1 + + GET_DATE, FITS_date +; +; Read photon information in the first extension of the IDF. +; + a = mrdfits(input, 1, status=status, /fscale) +; +; Read timeline table from the third extension of the IDF. +; + b = mrdfits(input, 3, status=status, /fscale) + nseconds = n_elements(b.time) + if (nint(b.time[nseconds-1] - b.time[0]) gt nseconds) then $ + print, 'WARNING: '+input+' timeline table not continuous.' + + ; for k = channel, channel+4, 4 do begin + k = channel + 4 + print, 'Analyzing extension ', k + ; + ; Read corresponding wavelength array from WAVE_CAL file. + ; + wcalfile=strtrim('/data1/fuse/calfuse/v3.0/calfiles/'+wavecal,2) + print, 'Reading wavelength data from ', wcalfile + w=mrdfits(wcalfile,k,w_hdr) + wave = w.wavelength + wpc = SXPAR(w_hdr, 'DISPAPIX') + ; + ; Select photons from aperture of interest. + ; + g = where (a.channel eq k and (a.timeflgs and not 1B) eq 0 and $ + (a.loc_flgs and 16B) eq 0, n) + print, 'Number of photons in spectrum = ', n + x = a.x[g] + ; + ; Correct X values for Doppler shift. + ; + print, 'Applying Doppler shift' + t = nint(a.time[g] - b.time[0]) + ovc = (b.ORBITAL_VEL + vhelio)/C + for i = 0L, n-1L do x[i] = x[i] + ovc[t[i]] * a.lambda[g[i]]/wpc + ; + ; Extract spectrum and populate output arrays. + ; + print, 'Extracting spectrum' + x = nint(x) + flux = fltarr(16384) + error = fltarr(16384) + cntserr = fltarr(16384) + for i = 0L, n-1L do flux[x[i]] = flux[x[i]] + a.ergcm2[g[i]] + flux = flux / exptime / abs(wpc) + counts = histogram(x, min=0, max=16383) + g = where (counts gt 0) + cntserr[g] = sqrt(counts[g]) + error[g] = flux[g]/cntserr[g] + quality = indgen(16384) + ; + ; Create primary HDU of output file + ; + if (k lt 5) then ch = 'lif' else ch = 'sic' + filename = strcompress(rootname+strlowcase(detector)+ch+ap_code[ap]+'ttagfcal.fit',/remove_all) + FXADDPAR, hdr, 'DATE', FITS_date + FXADDPAR, hdr, 'FILENAME', filename + FXADDPAR, hdr, 'FILETYPE', 'CALIBRATED EXTRACTED SPECTRUM' + FXADDPAR, hdr, 'APER_ACT', strupcase(strcompress(ap_array[ap]+'_'+ch,/remove_all)) + FXHMAKE, hdr, /EXTEND + FXWRITE, filename, hdr + + FXBHMAKE, HEADER, 1, 'FUSE 1D Spectrum' + FXBADDCOL, iwave, HEADER, wave, 'WAVE', TUNIT='ANGSTROMS' + FXBADDCOL, iflux, HEADER, flux, 'FLUX', TUNIT='ERG/CM2/S/A' + FXBADDCOL, ierror, HEADER, error, 'ERROR', TUNIT='ERG/CM2/S/A' + FXBADDCOL, iquality, HEADER, quality, 'QUALITY', TUNIT='UNITLESS' + FXBADDCOL, icounts, HEADER, counts, 'COUNTS', TUNIT='COUNTS' + FXBADDCOL, icntserr, HEADER, cntserr, 'CNTSERR', TUNIT='COUNTS' + + FXBCREATE, out, filename, HEADER + FXBWRITE, out, wave, iwave, 1 + FXBWRITE, out, flux, iflux, 1 + FXBWRITE, out, error, ierror, 1 + FXBWRITE, out, quality, iquality, 1 + FXBWRITE, out, counts, icounts, 1 + FXBWRITE, out, cntserr, icntserr, 1 + FXBFINISH, out + ;endfor + + close, /all + return + end -- cgit