;+ ; NAME: ; adderr ; ; PURPOSE: ; Add or subtract two arrays and propagate their statistical errors ; ; CATEGORY: ; Statistics, error propagation ; ; CALLING SEQUENCE: ; sum = adderr(array1,errarray1,array2,errarray2,[ERR=sumerr,[SUBTRACT=subtract]]) ; ; INPUTS: ; array1 -- Array to add to or subtract from ; errarray1 -- Array containing the statistical errors of array1 ; array2 -- Array to add to or subtract from array1 ; errarray1 -- Array containing the statistical errors of array1 ; ; KEYWORD PARAMETERS: ; subtract -- Set this keyword to perform subtraction instead of addition ; ; OUTPUTS: ; The function returns the sum (or difference, if the subtract keyword is set) ; of the two input arrays ; ; OPTIONAL OUTPUTS: ; err -- Set this keyword to a variable that will contain the statistical ; error associated with the summed (differenced) array ; ; COMMON BLOCKS: ; None ; ; SIDE EFFECTS: ; None ; ; RESTRICTIONS: ; None ; ; MODIFICATION HISTORY: ; 2007-02-19 AJS Added documentation ;- FUNCTION adderr, array1, errarray1, array2, errarray2, ERR=sumerr, SUBTRACT=subtract IF (NOT KEYWORD_SET(subtract)) THEN sum = array1 + array2 $ ELSE sum = array1 - array2 sumerr = SQRT(errarray1^2 + errarray2^2) RETURN, sum END ; Simple program to propagate errors through a division process function diverr,u,erru,v,errv,err=errx,mult=mult if not keyword_set(mult) then x=u/v else x=u*v u1=u v1=v uzero=where(u eq 0,uzerocnt) vzero=where(v eq 0,vzerocnt) if uzerocnt gt 0 then u1[uzero]=!values.f_nan if vzerocnt gt 0 then v1[vzero]=!values.f_nan if total(finite(x)) lt n_elements(x) then errx=sqrt(x^2)*sqrt(erru^2/u1^2+$ errv^2/v1^2) else errx=abs(x)*sqrt(erru^2/u1^2+errv^2/v1^2) return,x end ;+ ; NAME: ; lamp_mike_deadtime ; ; PURPOSE: ; Apply the deadtime correction to LAMP data ; ; CATEGORY: ; Data processing ; ; CALLING SEQUENCE: ; deadtime_correction = lamp_mike_deadtime(cntrate) ; ; INPUTS: ; cntrate -- the measured countrate (Hz) of the data ; ; OUTPUTS: ; Returns the factor by which to multiply the data to yield the ; count rate that would have been detected given infinitely fast ; electronics. C_out = C_in * EXP(d*C_in) ; ; COMMON BLOCKS: ; None. ; ; SIDE EFFECTS: ; None. ; ; NOTES: ; The deadtime correction time constant is 18 microseconds. At ; input rates below 50 kHz, the detector electronics is non-paralyzable ; (fixed deadtime per event that is not re-triggerable). To calculate ; the detector output count rate, the following formula is used: ; ; output_rate = (input_rate)/(1 + input_rate * tau) ; ; with an input_rate = 10,000 c/s and tau = 18 microseconds, the ; calculated output_rate is 8474.6 c/s. This formual can be inverted ; to calculate the input rate knowing tau and the output rate. ; ; MODIFICATION HISTORY: ; V1.0 2008-11-15 D.E. Kaufmann (adapted from Andrew Steffl's routine ; for P-Alice) ;- FUNCTION lamp_mike_deadtime, cntrate tau = !LAMP_DEADTIME ; defined in lamp_mike_sysv.pro, thus this procedure ; must be run prior in order to define the constant deadtime_correction = 1. / (1. - tau * cntrate) RETURN, deadtime_correction END ;+ ; NAME: ; lamp_mike_event_locate ; ; PURPOSE: ; Locate a pixellist event on the LAMP detector ; ; CATEGORY: ; Data processing ; ; CALLING SEQUENCE: ; location = lamp_mike_event_locate(x, y) ; ; INPUTS: ; x -- the detector x location of the event ; y -- the detector y location of the event ; ; OUTPUTS: ; Returns the nonnegative constant defined in lamp_mike_sysv.pro ; for the region of the detector in which the event was detected. ; If outside all defined regions, returns -1. ; ; COMMON BLOCKS: ; None. ; ; SIDE EFFECTS: ; None. ; ; NOTES: ; Current constants are: ; FOV = 0 ; DARK1 = 1 ; DARK2 = 2 ; DARK3 = 3 ; DARK4 = 4 ; STIM1 = 5 ; STIM2 = 6 ; STRAY1 = 7 ; STRAY2 = 8 ; ; MODIFICATION HISTORY: ; V1.0 15 Jan 2009 D.E. Kaufmann ; V1.1 01 Mar 2010 D.E. Kaufmann - added STRAY1 and STRAY2 regions ; ;- FUNCTION lamp_mike_event_locate, x, y location = -1 IF ((x GE !LAMP_FOV_X_MIN) AND (x LE !LAMP_FOV_X_MAX) AND (y GE !LAMP_FOV_Y_MIN) AND (y LE !LAMP_FOV_Y_MAX)) THEN location = !LAMP_FOV_ID $ ELSE IF ((x GE !LAMP_DARK1_X_MIN) AND (x LE !LAMP_DARK1_X_MAX) AND (y GE !LAMP_DARK1_Y_MIN) AND (y LE !LAMP_DARK1_Y_MAX)) THEN location = !LAMP_DARK1_ID $ ELSE IF ((x GE !LAMP_DARK2_X_MIN) AND (x LE !LAMP_DARK2_X_MAX) AND (y GE !LAMP_DARK2_Y_MIN) AND (y LE !LAMP_DARK2_Y_MAX)) THEN location = !LAMP_DARK2_ID $ ELSE IF ((x GE !LAMP_DARK3_X_MIN) AND (x LE !LAMP_DARK3_X_MAX) AND (y GE !LAMP_DARK3_Y_MIN) AND (y LE !LAMP_DARK3_Y_MAX)) THEN location = !LAMP_DARK3_ID $ ELSE IF ((x GE !LAMP_DARK4_X_MIN) AND (x LE !LAMP_DARK4_X_MAX) AND (y GE !LAMP_DARK4_Y_MIN) AND (y LE !LAMP_DARK4_Y_MAX)) THEN location = !LAMP_DARK4_ID $ ELSE IF ((x GE !LAMP_STIM1_X_MIN) AND (x LE !LAMP_STIM1_X_MAX) AND (y EQ !LAMP_STIM1_Y)) THEN location = !LAMP_STIM1_ID $ ELSE IF ((x GE !LAMP_STIM2_X_MIN) AND (x LE !LAMP_STIM2_X_MAX) AND (y EQ !LAMP_STIM2_Y)) THEN location = !LAMP_STIM2_ID $ ELSE IF ((x GE !LAMP_STRAY1_X_MIN) AND (x LE !LAMP_STRAY1_X_MAX) AND (y GE !LAMP_STRAY1_Y_MIN) AND (y LE !LAMP_STRAY1_Y_MAX)) THEN location = !LAMP_STRAY1_ID $ ELSE IF ((x GE !LAMP_STRAY2_X_MIN) AND (x LE !LAMP_STRAY2_X_MAX) AND (y GE !LAMP_STRAY2_Y_MIN) AND (y LE !LAMP_STRAY2_Y_MAX)) THEN location = !LAMP_STRAY2_ID RETURN, location END ;+ ; NAME: ; lamp_mike_log ; ; PURPOSE: ; Write status message to log, and to screen if desired ; ; CATEGORY: ; Status logging, I/O ; ; CALLING SEQUENCE: ; lamp_mike_log, message, [ERRNUM=errnum, INITFLAG=initflag, ; STARTFLAG=startflag, TIMEFLAG=timeflag, ; VERBOSELEVEL=verboselevel] ; ; OPTIONAL INPUTS: ; None ; ; KEYWORD PARAMETERS: ; ERRNUM = value: 0 or not specified = normal msg, 1 = warning, ; 2 = fatal error. ; /INITFLAG: Flag to (re)initialize message log ; /STARTFLAG: Flag to place initial start message in log ; /TIMEFLAG: Flag to insert the current time in log message ; VERBOSELEVEL = value: Sets the level of verbosity at which ; this message gets printed to the screen. ; ; OUTPUTS: ; None ; ; COMMON BLOCKS: ; mike_data ; ; SIDE EFFECTS: ; None ; ; RESTRICTIONS: ; None ; ; PROCEDURES CALLED: ; None ; ; EXAMPLE: ; lamp_mike_log, message, /STARTFLAG, /INITFLAG, /TIMEFLAG, VERBOSELEVEL = 1 ; ; MODIFICATION HISTORY: ; V1.0 18 March 2009 DEK initial implementation ; ;- PRO lamp_mike_log, message, $ ERRNUM = errnum, $ INITFLAG = initflag, $ STARTFLAG = startflag, $ TIMEFLAG = timeflag, $ VERBOSELEVEL = verboselevel COMMON mike_data ; Return immediately if logging disabled IF (mikeflags.log EQ 0) THEN RETURN ; Initialize if requested IF (KEYWORD_SET(initflag)) THEN BEGIN logmessages = '' mikeflags.error = 0 ENDIF ; If verboselevel not set, assume message UNimportant and set vlevel to 3 IF (KEYWORD_SET(verboselevel)) THEN vlevel = verboselevel ELSE vlevel = 3 ; Initialize log entry and adjust vlevel according to errnum IF (N_ELEMENTS(errnum) EQ 0) THEN errnum = 0 CASE (errnum) OF 0 : log_entry = '' 1 : log_entry = '# WARNING: ' 2 : log_entry = '## FATAL ERROR: ' ENDCASE IF (errnum NE 0) THEN mikeflags.error = errnum > mikeflags.error IF (errnum EQ 1) THEN vlevel = vlevel < 2 IF (errnum EQ 2) THEN vlevel = vlevel < 1 ; Process startflag, timeflag, and message IF (KEYWORD_SET(startflag)) THEN BEGIN HELP, CALLS = calls log_entry = log_entry + 'Starting: ' + (STRSPLIT(calls[1], /EXTRACT))[0] ENDIF IF (KEYWORD_SET(message)) THEN log_entry = log_entry + message IF (KEYWORD_SET(timeflag)) THEN $ log_entry[timeflag-1] = log_entry[timeflag-1] + ' {' + SYSTIME() + '}' ; Add entry to log logmessages = [logmessages, log_entry] ; Print entry to screen if vlevel is low enough IF (mikeflags.verbose GE vlevel) THEN PRINT, log_entry, FORMAT = '(A)' END ;+ ; NAME: ; lamp_mike ; ; PURPOSE: ; Process CODMAC Level 2 (Engineering Level) LAMP data ; into CODMAC Level 3 (Science Level) LAMP data ; ; CATEGORY: ; Data processing ; ; CALLING SEQUENCE: ; lamp_mike, INFILES=infiles, OUTFILES=outfiles, STATFILES=statfiles, ; CALIBRATION_DIR=calibration_dir, SPICE_DIR=spice_dir, ; TEMP_DIR=temp_dir, AEFFFLAG=aeffflag, AEFFFILE=aefffile, ; BADFLAG=badflag, BADFILE=badfile, DARKFLAG=darkflag, ; DARKFILE=darkfile, DEADFLAG=deadflag, FLATFLAG=flatflag, ; FLATFILE=flatfile, RECTFLAG=rectflag, METAKERNEL=metakernel, ; LOGFILE=logfile, VERBOSE=verbose ; ; OPTIONAL INPUTS: ; None. ; ; KEYWORD PARAMETERS: ; INFILES = value: String array of input CODMAC Level 2 file names to ; process. If not supplied, procedure will attempt to ; get this input from environment variable MIKE_IN_FILE. ; If not supplied in either way, procedure errors out. ; OUTFILES = value: String array of output CODMAC Level 3 file names to ; write. If not supplied, procedure will attempt to get ; this input from environment variable MIKE_OUT_FILE. ; Number of outfiles must match number of infiles. ; STATFILES = value: String array of output processing status file names ; to write. If not supplied, procedure will attempt to ; get this input from environment variable MIKE_OUT_STATUS. ; If specified, number of output status files must match ; number of infiles. ; CALIBRATION_DIR = value: String specifying location of root directory ; for calibration data. If not specified, procedure will ; attempt to use default value. ; SPICE_DIR = value: String specifying location of root directory for SPICE ; kernel data. If not specified, procedure will attempt ; to use default value. ; TEMP_DIR = value: String specifying directory location to use for ; generating temporary files, if needed. If not specified, ; procecure will attempt to use default value. ; /AEFFFLAG: Flag to apply effective area calibration (default = 1). ; AEFFFILE = value: String specifying effective area calibration file. ; This must be supplied if AEFFFLAG is set. ; /BADFLAG: Flag to apply bad pixel filtering (default = 0). [NOT IMPLEMENTED] ; BADFILE = value: String specifying bad pixel mask. This must be supplied ; if BADFLAG is set. [NOT IMPLEMENTED] ; /DARKFLAG: Flag to apply dark count correction to data taken in histogram ; mode (default = 1). For data taken in pixel list mode, ; no dark count correction is applied at this level. ; Instead, the dark count level is characterized for later ; correction at the stage of map construction. ; DARKFILE = value: String specifying dark image data file to use for this ; correction. This must be supplied if DARKFLAG is set. ; /DEADFLAG: Flag to apply deadtime correction (default = 1). ; /FLATFLAG: Flag to apply flatfield correction (default = 1). ; FLATFILE = value: String specifying flatfield image data file to use for ; this correction. This must be specified if FLATFLAG is set. ; /RECTFLAG: Flag to apply geometrical distortion correction (default = 1). ; [NOT IMPLEMENTED] ; METAKERNEL = value: String specifying SPICE metakernel to use for performing ; geospatial location of the LAMP pixel list data. This ; must be specified, else the procedure exits with an error. ; LOGFILE = value: String specifying the name of a file to which processing ; log messages are written. ; VERBOSE = value: Integer (0-3) specifying verbose level in messaging. The ; higher the value, the more verbose the messaging. ; ; OUTPUTS: ; None. ; ; COMMON BLOCKS: ; mike_data: Contains processing flags and log message array ; ; SIDE EFFECTS: ; None. ; ; RESTRICTIONS: ; None. ; ; PROCEDURES CALLED: ; ; INTERNAL ; ======== ; lamp_mike_sysv ; lamp_mike_log ; poisson_errorbar ; lamp_mike_wavecal ; lamp_mike_event_locate ; lamp_mike_deadtime ; adderr ; ; IDLASTRO ; ======== ; rdfits_struct ; readfits ; readcol ; mrdfits ; fxpar ; sxaddpar ; mkhdr ; mwrfits ; fxbhmake ; fxbaddcol ; fxbcreate ; fxbwritm ; fxbfinish ; ; ICY (SPICE) ; =========== ; cspice_furnsh ; cspice_ktotal ; cspice_kdata ; cspice_pxform ; cspice_getfov ; cspice_mxv ; cspice_gdpool ; cspice_scs2e ; cspice_subpt ; cspice_srfxpt ; cspice_reclat ; cspice_illum ; cspice_recrad ; cspice_spkssb ; cspice_spkezp ; cspice_subpnt ; cspice_vsep ; cspice_subslr ; ; EXAMPLE: ; lamp_mike, /VERBOSE, LOGFILE = 'log', DARKFILE = 'pa_dark_001.fit', ; FLATFILE = 'pa_flat_000.fit', AEFFFILE = 'lamp_aeff_000.tab', ; METAKERNEL = 'meta_kernellro.tm' ; ; MODIFICATION HISTORY: ; 17 Mar 2009, Initial implementation. David E. Kaufmann ; 02 Mar 2010, Added code for stray factor. David E. Kaufmann ; 27 Apr 2010, Updated extension 6 contents from "Calibrated Calculated ; Count Rate" to "Reduced Count Rate" based on Kurt's email ; of 08 March 2010. David E. Kaufmann ; 01 Feb 2011, Improved performance by adding "BUFFERSIZE=1048576L" to ; fxbwritm calls. ; ;- PRO lamp_mike, $ INFILES = infiles, $ OUTFILES = outfiles, $ STATFILES = statfiles, $ CALIBRATION_DIR = calibration_dir, $ SPICE_DIR = spice_dir, $ TEMP_DIR = temp_dir, $ AEFFFLAG = aeffflag, $ AEFFFILE = aefffile, $ BADFLAG = badflag, $ BADFILE = badfile, $ DARKFLAG = darkflag, $ DARKFILE = darkfile, $ DEADFLAG = deadflag, $ FLATFLAG = flatflag, $ FLATFILE = flatfile, $ RECTFLAG = rectflag, $ METAKERNEL = metakernel, $ LOGFILE = logfile, $ VERBOSE = verbose COMMON mike_data, mikeflags, logmessages ; Initialize Mike system variables (program constants) lamp_mike_sysv ; Initialize error reporting start_time = SYSTIME(1); seconds since January 1, 1970 error_status = 0 IF (N_ELEMENTS(mikeflags) EQ 0) THEN $ mikeflags = {version: !MIKE_SL2VER + '[' + !MIKE_SL2DATE + ']', $ mode: '', $ verbose: KEYWORD_SET(verbose) ? verbose : 0, $ log: 1, $ error: 0, $ starttime: start_time, $ statflag: 0} ; Initialize message logging lamp_mike_log, message, /STARTFLAG, /INITFLAG, /TIMEFLAG, VERBOSELEVEL = 1 ; Set default processing flags (1 = apply, 0 = do NOT apply) IF (N_ELEMENTS(aeffflag) EQ 0) THEN aeffflag = 1 IF (N_ELEMENTS(badflag) EQ 0) THEN badflag = 0 IF (N_ELEMENTS(darkflag) EQ 0) THEN darkflag = 1 IF (N_ELEMENTS(deadflag) EQ 0) THEN deadflag = 1 IF (N_ELEMENTS(flatflag) EQ 0) THEN flatflag = 1 IF (N_ELEMENTS(rectflag) EQ 0) THEN rectflag = 1 ;; [[DEK REMINDER: Add code below to implement the geometric distortion correction ;; controlled by rectflag whenever we define the needed algorithm.]] ; Validate input (command line arguments take precedence over environment variables) ; Infile(s) IF (N_ELEMENTS(infiles) EQ 0) THEN BEGIN env_infile = GETENV('MIKE_IN_FILE') IF (KEYWORD_SET(env_infile)) THEN result = EXECUTE('infiles = ' + env_infile) ENDIF nfiles = N_ELEMENTS(infiles) IF (nfiles EQ 0) THEN BEGIN message = 'No input file(s) specified.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF ; Outfile(s) IF (N_ELEMENTS(outfiles) EQ 0) THEN BEGIN env_outfile = GETENV('MIKE_OUT_FILE') IF (KEYWORD_SET(env_outfile)) THEN result = EXECUTE('outfiles = ' + env_outfile) ENDIF ntmp = N_ELEMENTS(outfiles) IF (ntmp NE nfiles) THEN BEGIN message = 'Number of output files does not match number of input files.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF ; Statfile(s) IF (N_ELEMENTS(statfiles) EQ 0) THEN BEGIN env_statfile = GETENV('MIKE_OUT_STATUS') IF (KEYWORD_SET(env_statfile)) THEN result = EXECUTE('statfiles = ' + env_statfile) ENDIF ntmp = N_ELEMENTS(statfiles) IF (ntmp GT 0) THEN BEGIN mikeflags.statflag = 1 IF (ntmp NE nfiles) THEN BEGIN message = 'Number of output status files is nonzero and does not match number of input files.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF ENDIF ; Calibration directory IF (N_ELEMENTS(calibration_dir) EQ 0) THEN BEGIN env_calibration_dir = GETENV('MIKE_CALIBRATION_DIR') IF (KEYWORD_SET(env_calibration_dir)) THEN calibration_dir = env_calibration_dir ENDIF IF (N_ELEMENTS(calibration_dir) EQ 0) THEN calibration_dir = !MIKE_DEFAULT_CALIBRATION_DIR ; SPICE directory IF (N_ELEMENTS(spice_dir) EQ 0) THEN BEGIN env_spice_dir = GETENV('MIKE_SPICE_DIR') IF (KEYWORD_SET(env_spice_dir)) THEN spice_dir = env_spice_dir ENDIF IF (N_ELEMENTS(spice_dir) EQ 0) THEN spice_dir = !MIKE_DEFAULT_SPICE_DIR ; Temporary file directory IF (N_ELEMENTS(temp_dir) EQ 0) THEN BEGIN env_temp_dir = GETENV('MIKE_TEMP_DIR') IF (KEYWORD_SET(env_temp_dir)) THEN temp_dir = env_temp_dir ENDIF IF (N_ELEMENTS(temp_dir) EQ 0) THEN temp_dir = !MIKE_DEFAULT_TEMP_DIR ; Dark correction file IF (KEYWORD_SET(darkflag)) THEN BEGIN IF (KEYWORD_SET(darkfile)) THEN BEGIN darkfile = calibration_dir + PATH_SEP() + 'dark' + PATH_SEP() + darkfile rdfits_struct, darkfile, darkstruct, /SILENT darkimage = darkstruct.im0 darkerror = darkstruct.im1 ENDIF ELSE BEGIN message = 'Dark correction file not specified.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDELSE ENDIF ; Flatfield correction file ; flatAO : HVPS = A , aperture door = Open ; flatBO : HVPS = B , aperture door = Open ; flatABO: HVPS = both, aperture door = Open ; flatAC : HVPS = A , aperture door = Closed ; flatBC : HVPS = B , aperture door = Closed ; flatABC: HVPS = both, aperture door = Closed IF (KEYWORD_SET(flatflag)) THEN BEGIN IF (KEYWORD_SET(flatfile)) THEN BEGIN flatfile = calibration_dir + PATH_SEP() + 'flat' + PATH_SEP() + flatfile flatAO = mrdfits(flatfile, 1, hdr1, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading ' + flatfile + ', status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF flatAC = mrdfits(flatfile, 2, hdr2, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading ' + flatfile + ', status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF flatBO = mrdfits(flatfile, 3, hdr3, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading ' + flatfile + ', status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF flatBC = mrdfits(flatfile, 4, hdr4, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading ' + flatfile + ', status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF flatABO = mrdfits(flatfile, 5, hdr5, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading ' + flatfile + ', status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF flatABC = mrdfits(flatfile, 6, hdr6, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading ' + flatfile + ', status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF ENDIF ELSE BEGIN message = 'Flatfield correction file not specified.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDELSE ENDIF ; Effective area calibration file ; aeffAO : HVPS = A , aperture door = Open ; aeffBO : HVPS = B , aperture door = Open ; aeffABO: HVPS = both, aperture door = Open ; aeffAC : HVPS = A , aperture door = Closed ; aeffBC : HVPS = B , aperture door = Closed ; aeffABC: HVPS = both, aperture door = Closed IF (KEYWORD_SET(aeffflag)) THEN BEGIN IF (KEYWORD_SET(aefffile)) THEN BEGIN aefffile = calibration_dir + PATH_SEP() + 'aeff' + PATH_SEP() + aefffile tab1 = mrdfits(aefffile, 1, hdr1, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading ' + aefffile + ', status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF aeffwave = tab1.WAVELENGTH aeffAO = tab1.EFFECTIVE_AREA_A aeffBO = tab1.EFFECTIVE_AREA_B aeffABO = tab1.EFFECTIVE_AREA_BOTH apratio = fxpar(hdr1, 'APRATIO'); airglow-to-pinhole aperture ratio aeffAC = aeffAO / apratio aeffBC = aeffBO / apratio aeffABC = aeffABO / apratio ENDIF ELSE BEGIN message = 'Effective area calibration file not specified.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDELSE ENDIF ;; [[DEK REMINDER: Be prepared to add code later to handle possible temperature ;; dependence of these calibration steps, once the nature of the ;; dependence is better understood and characterized.]] ;; [[DEK REMINDER: Dark signals outside FOV have been showing extra signal dependent ;; on the source in the FOV for R- and P-Alices. This is difficult ;; to decipher, but most likely some combination of instrumental ;; scattering, electronic noise, etc. Make sure we get a robust set ;; of scattered light measurements during commissioning and be ;; prepared to add code to correct for these effects later on. ; SPICE metakernel file IF (KEYWORD_SET(metakernel)) THEN BEGIN PUSHD, spice_dir cspice_furnsh, metakernel POPD cspice_ktotal, 'ALL', nkernels kernel_list = STRARR(nkernels) FOR k = 0, nkernels - 1 DO BEGIN cspice_kdata, k, 'ALL', file, type, source, handle, found IF (found) THEN kernel_list[k] = FILE_BASENAME(file) ENDFOR ENDIF ELSE BEGIN message = 'SPICE metakernel not specified.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDELSE ; Initialize instrument pointing vectors ; get transformation matrix from LAMP frame to SC_BUS frame cspice_pxform, !SPICE_LAMP_FRAME, !SPICE_SC_BUS_FRAME, 0.0D0, rmatrix; transform assumed time-independent, so set input et = 0 ; get boresight vector cspice_getfov, !SPICE_LRO_LAMP_ID, 66, shape, frame, bsight, bounds ; transform boresight vector to SC_BUS frame cspice_mxv, rmatrix, bsight, bsight ; get pixel center vectors cspice_gdpool, 'INS-85300_PIXEL_CENTERS', 0, 96, values, found IF (found) THEN BEGIN nvalues = N_ELEMENTS(values) pixel_centers = FLTARR(3, nvalues/3, /NOZERO) FOR k = 0, nvalues - 1 DO pixel_centers[k] = values[k] ENDIF ELSE BEGIN message = 'In cspice_gdpool(), unable to find pixel centers.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDELSE ; transform pixel center vectors to SC_BUS frame pixel_centers = TRANSPOSE(MATRIX_MULTIPLY(pixel_centers, rmatrix, /ATRANSPOSE)) ; get pixel corner vectors cspice_gdpool, 'INS-85300_FOV_BOUNDARY_CORNERS', 0, 198, values, found IF (found) THEN BEGIN nvalues = N_ELEMENTS(values) pixel_corners = FLTARR(3, nvalues/3, /NOZERO) FOR k = 0, nvalues - 1 DO pixel_corners[k] = values[k] ENDIF ELSE BEGIN message = 'In cspice_gdpool(), unable to find pixel corners.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDELSE ; transform pixel corner vectors to SC_BUS frame pixel_corners = TRANSPOSE(MATRIX_MULTIPLY(pixel_corners, rmatrix, /ATRANSPOSE)) ; Get vector of Poisson errors poisson_errors = DBLARR(2501, /NOZERO) FOR i = 0, 2500 DO poisson_errors[i] = poisson_errorbar(i, /DOUBLE) poisson_errors = FLOAT(poisson_errors) ; Get LRO orbit number and SCUT time vectors from ASCII database file readcol, !MIKE_LRO_ORBIT_DATABASE_FILE, lro_orbit_number, lro_orbit_scut, FORMAT = 'I,D', /SILENT n_lro_orbits = N_ELEMENTS(lro_orbit_number) ; Determine the total number of detector elements in the 4 regions where ; dark signal is being estimated n_dark_elements = (!LAMP_DARK1_X_MAX - !LAMP_DARK1_X_MIN + 1)*(!LAMP_DARK1_Y_MAX - !LAMP_DARK1_Y_MIN + 1) + $ (!LAMP_DARK2_X_MAX - !LAMP_DARK2_X_MIN + 1)*(!LAMP_DARK2_Y_MAX - !LAMP_DARK2_Y_MIN + 1) + $ (!LAMP_DARK3_X_MAX - !LAMP_DARK3_X_MIN + 1)*(!LAMP_DARK3_Y_MAX - !LAMP_DARK3_Y_MIN + 1) + $ (!LAMP_DARK4_X_MAX - !LAMP_DARK4_X_MIN + 1)*(!LAMP_DARK4_Y_MAX - !LAMP_DARK4_Y_MIN + 1) ; Determine the total number of detector elements in the 2 regions where ; the stray light signal is being estimated n_stray_elements = (!LAMP_STRAY1_X_MAX - !LAMP_STRAY1_X_MIN + 1)*(!LAMP_STRAY1_Y_MAX - !LAMP_STRAY1_Y_MIN + 1) + $ (!LAMP_STRAY2_X_MAX - !LAMP_STRAY2_X_MIN + 1)*(!LAMP_STRAY2_Y_MAX - !LAMP_STRAY2_Y_MIN + 1) message = STRING(nfiles) + ' file(s) to process.' lamp_mike_log, message, VERBOSELEVEL = 1 ; Loop over input files and process FOR ifile = 0, nfiles - 1 DO BEGIN istat_first = N_ELEMENTS(logmessages) ; index into logmessages of first entry for this file ; Verify file exists and report processing start found = FILE_TEST(infiles[ifile]) IF (NOT found) THEN BEGIN message = 'File ' + infiles[ifile] + ' not found.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF filenum = STRING(ifile + 1, STRTRIM(nfiles, 2), FORMAT = '(I5,"/",A)') IF (KEYWORD_SET(verbose)) THEN $ PRINT, '% ' + !MIKE_SL2PROC + ': Now processing file ' + infiles[ifile] + ' ' + filenum IF (STRPOS(STRUPCASE(infiles[ifile]), 'SCI') GE 0) THEN BEGIN message = 'File has already been processed!' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Make sure file is valid FITS buffer = BYTARR(10, /NOZERO) OPENR, lun, infiles[ifile], /GET_LUN READU, lun, buffer FREE_LUN, lun IF (STRING(buffer) NE 'SIMPLE = ') THEN BEGIN message = 'File ' + infiles[ifile] + ' is not a valid FITS file.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Read in data ; Door open histogram, primary extension, image im0 = mrdfits(infiles[ifile], 0, hdr0, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading extension 0 with mrdfits, status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Door closed histogram, extension 1, image im1 = mrdfits(infiles[ifile], 1, hdr1, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading extension 1 with mrdfits, status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Acquisition list, extension 2, ASCII table tab2 = mrdfits(infiles[ifile], 2, hdr2, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading extension 2 with mrdfits, status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Raw science frame data, extension 3, binary table tab3 = mrdfits(infiles[ifile], 3, hdr3, /UNSIGNED, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading extension 3 with mrdfits, status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Calculated countrate, extension 4, binary table tab4 = mrdfits(infiles[ifile], 4, hdr4, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading extension 4 with mrdfits, status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; High-rate LTS data, extension 5, binary table tab5 = mrdfits(infiles[ifile], 5, hdr5, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading extension 5 with mrdfits, status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Housekeeping table, extension 6, binary table tab6 = mrdfits(infiles[ifile], 6, hdr6, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading extension 6 with mrdfits, status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Check file HK status good_housekeeping = 1 nhkpackets = fxpar(hdr0, 'PACKETS') hkgaps = fxpar(hdr0, 'HKGAPS') IF (nhkpackets EQ 0) THEN BEGIN good_housekeeping = 0 message = 'No housekeeping dataset detected! Processing of current file aborted.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ELSE IF (hkgaps EQ 1) THEN BEGIN message = 'Partial housekeeping dataset detected! Attempting normal processing.' lamp_mike_log, message, ERRNUM = 1, VERBOSELEVEL = verbose ENDIF ; Create wavelength calibration image for output FITS using average detector housing temperature; ; STIM positions, if available, will be used instead of this for calibration of pixellist data avgdetetemp = 0.0D0 IF (good_housekeeping) THEN avgdetetemp = fxpar(hdr0, 'AVRDETHT') ; waveimage = lamp_mike_wavecal(DETETEMP = avgdetetemp) waveimage = lamp_mike_wavecal() ; DEK, use "fiducial" wavelength solution here ;; [[DEK REMINDER: Note that the above image is in some sense "time-averaged" and this could be ;; misleading to end users. We need to be very clear about this when delivering ;; the RDR product files and instructing on their use.]] ; SCUT = hack_clock * hack_time (wrap-extended) + hack_offset hack_clock = fxpar(hdr0, 'HACKCLCK') hack_offset = fxpar(hdr0, 'HACKOFFS') ; Generate some acquisition statistics nframes = fxpar(hdr3, 'NFRAMES') IF (nframes EQ 0) THEN BEGIN message = 'No science dataset detected! Processing of current file aborted.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF npixellists = 0L nhistograms = 0L nfovevents = 0L ; photon events in detector FOV ncrentries = 0L ; number of countrate table entries hvps_state_sc = -1 ; HVPS state for all pixellist frames in this file; treat as error if it changes first_hack = LONARR(nframes) - 1L skip_frame = BYTARR(nframes) ncrentries_per_frame = LONARR(nframes) FOR iframe = 0L, nframes - 1 DO BEGIN ; iframe indexes into the raw science frame data table (extension 3) ; and the acquisition list table (extension 2) frame_data = (tab3.FRAME_DATA)[*, iframe] frame_number = ISHFT(ISHFT(frame_data[0], 4), -4) ; make sure HK data for frame is present in Acquisition List IF (((tab2.APD_STATE)[iframe] EQ 0) && ((tab2.LTS_STATE)[iframe] EQ 0) && $ ((tab2.STIM_STATE)[iframe] EQ 0) && ((tab2.HVPS_STATE)[iframe] EQ 0)) THEN BEGIN message = 'No housekeeping data for frame ' + STRING(frame_number) + ', skipping...' lamp_mike_log, message, ERRNUM = 1, VERBOSELEVEL = verbose skip_frame[iframe] = 1 CONTINUE ENDIF ; make sure this is a data frame IF ((tab2.FR_PATTERN)[iframe] NE !LAMP_DATA_PATTERN) THEN BEGIN message = 'Frame ' + STRING(frame_number) + ' is not a data frame, skipping...' lamp_mike_log, message, ERRNUM = 1, VERBOSELEVEL = verbose skip_frame[iframe] = 1 CONTINUE ENDIF ; make sure STIMs are on IF ((tab2.STIM_STATE)[iframe] NE !LAMP_STIM_ON) THEN BEGIN message = 'STIM pixels not ON in frame ' + STRING(frame_number) + ', skipping...' lamp_mike_log, message, ERRNUM = 1, VERBOSELEVEL = verbose skip_frame[iframe] = 1 CONTINUE ENDIF ; make sure selected HVPS state is ok if aeff calibration is requested IF (KEYWORD_SET(aeffflag)) THEN BEGIN hvps_state_frame = (tab2.HVPS_STATE)[iframe] IF ((hvps_state_frame NE !LAMP_HVPS_A_ON) && (hvps_state_frame NE !LAMP_HVPS_B_ON) && $ (hvps_state_frame NE !LAMP_HVPS_AB_ON)) THEN BEGIN message = 'Bad HVPS state = ' + STRING(hvps_state_frame) + ' in frame ' + STRING(frame_number) + ', skipping...' lamp_mike_log, message, ERRNUM = 1, VERBOSELEVEL = verbose skip_frame[iframe] = 1 CONTINUE ENDIF IF ((tab2.ACQ_MODE)[iframe] EQ !LAMP_PIXELLIST_MODE) THEN BEGIN IF (hvps_state_sc LT 0) THEN hvps_state_sc = hvps_state_frame $ ELSE IF (hvps_state_frame NE hvps_state_sc) THEN BEGIN message = 'HVPS state not constant for all pixellist frames in this file, aborting...' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose GOTO, END_OF_LOOP ENDIF ENDIF ENDIF IF ((tab2.ACQ_MODE)[iframe] EQ !LAMP_PIXELLIST_MODE) THEN BEGIN ; pixellist mode npixellists += 1 nhacks = 0L ; total number of time hacks in this frame FOR iword = 1L, 32767L DO BEGIN IF (frame_data[iword] NE 0) THEN BEGIN ; frame data word is not zero fill IF (ISHFT(frame_data[iword], -15) EQ 0) THEN BEGIN ; photon event xtmp = ISHFT(ISHFT(frame_data[iword], 6), -6) ytmp = ISHFT(ISHFT(frame_data[iword], 1), -11) location = lamp_mike_event_locate(xtmp, ytmp) IF (location EQ !LAMP_FOV_ID) THEN nfovevents += 1L ENDIF ELSE BEGIN ; time hack IF (first_hack[iframe] LT 0L) THEN first_hack[iframe] = ISHFT(ISHFT(frame_data[iword], 1), -1) nhacks += 1 ENDELSE ENDIF ENDFOR IF ((tab2.HACK_RATE)[iframe] GT 0) THEN BEGIN nhacks_per_cr_interval = MAX([2^(6 - (tab2.HACK_RATE)[iframe]), 1]) ; interval is >= 128 ms ncrentries_per_frame[iframe] = (nhacks - 1) / nhacks_per_cr_interval ncrentries += ncrentries_per_frame[iframe] ENDIF estimated_first_hack = ROUND(((tab2.START_TIME)[iframe] - hack_offset)/hack_clock) IF (first_hack[iframe] GE 0L) THEN BEGIN nrollovers = ROUND((estimated_first_hack - first_hack[iframe])/32768.0D0) first_hack[iframe] += nrollovers*32768L ENDIF ELSE $ ; only when no hacks were detected in frame, thus should never occur first_hack[iframe] = estimated_first_hack ENDIF ELSE $ ; histogram mode nhistograms += 1L ENDFOR IF ((npixellists GT 0L) AND (nfovevents EQ 0L)) THEN BEGIN message = 'File contains pixellist mode data but no events within LAMP FOV' lamp_mike_log, message, ERRNUM = 1, VERBOSELEVEL = verbose ENDIF ; Allocate arrays for fast pixellist table IF (nfovevents GT 0L) THEN BEGIN hack_time = LONARR(nfovevents, /NOZERO) ; wrap-extended values detector_x = INTARR(nfovevents, /NOZERO) detector_y = INTARR(nfovevents, /NOZERO) spatial_row = INTARR(nfovevents, /NOZERO) wavelength = FLTARR(nfovevents, /NOZERO) scut_time = DBLARR(nfovevents, /NOZERO) weighted_count = FLTARR(nfovevents, /NOZERO) weighted_variance = FLTARR(nfovevents, /NOZERO) longitude = FLTARR(nfovevents, /NOZERO) latitude = FLTARR(nfovevents, /NOZERO) longitude_ll = FLTARR(nfovevents, /NOZERO) latitude_ll = FLTARR(nfovevents, /NOZERO) longitude_ul = FLTARR(nfovevents, /NOZERO) latitude_ul = FLTARR(nfovevents, /NOZERO) longitude_ur = FLTARR(nfovevents, /NOZERO) latitude_ur = FLTARR(nfovevents, /NOZERO) longitude_lr = FLTARR(nfovevents, /NOZERO) latitude_lr = FLTARR(nfovevents, /NOZERO) altitude = FLTARR(nfovevents, /NOZERO) incidence_angle = FLTARR(nfovevents, /NOZERO) emission_angle = FLTARR(nfovevents, /NOZERO) deadtime_factor = FLTARR(nfovevents, /NOZERO) pix_data_qual = BYTARR(nfovevents, /NOZERO) apdoor_state = BYTARR(nfovevents, /NOZERO) ; auxiliary array hvps_state = BYTARR(nfovevents, /NOZERO) ; auxiliary array cr_index = LONARR(nfovevents, /NOZERO) ; for cross referencing with the countrate arrays ENDIF ELSE BEGIN ; allocate dummy arrays for FITS file hack_time = LONARR(1, /NOZERO) ; wrap-extended values detector_x = INTARR(1, /NOZERO) detector_y = INTARR(1, /NOZERO) spatial_row = INTARR(1, /NOZERO) wavelength = FLTARR(1, /NOZERO) scut_time = DBLARR(1, /NOZERO) weighted_count = FLTARR(1, /NOZERO) weighted_variance = FLTARR(1, /NOZERO) longitude = FLTARR(1, /NOZERO) latitude = FLTARR(1, /NOZERO) longitude_ll = FLTARR(1, /NOZERO) latitude_ll = FLTARR(1, /NOZERO) longitude_ul = FLTARR(1, /NOZERO) latitude_ul = FLTARR(1, /NOZERO) longitude_ur = FLTARR(1, /NOZERO) latitude_ur = FLTARR(1, /NOZERO) longitude_lr = FLTARR(1, /NOZERO) latitude_lr = FLTARR(1, /NOZERO) altitude = FLTARR(1, /NOZERO) incidence_angle = FLTARR(1, /NOZERO) emission_angle = FLTARR(1, /NOZERO) deadtime_factor = FLTARR(1, /NOZERO) pix_data_qual = BYTARR(1, /NOZERO) ENDELSE IF (ncrentries GT 0L) THEN BEGIN cr_hack_time = LONARR(ncrentries, /NOZERO) cr_scut_time = DBLARR(ncrentries, /NOZERO) ; cr_countrate = FLTARR(ncrentries) cr_countrate = DBLARR(ncrentries) cr_hack_incr = INTARR(ncrentries, /NOZERO) ENDIF ; Allocate data cube for calibrated histogram and error image pairs IF (nhistograms GT 0L) THEN histogram_cube = FLTARR(1024, 32, 2, nhistograms) ; Allocate arrays for slow pixellist table (1 minute intervals) ; scut_scifile_start = fxpar(hdr0, 'SCSTRSEC') scut_scifile_start = MIN([fxpar(hdr0, 'SCSTRSEC'), (tab2.START_TIME)[0]]) ; scut_scifile_end = fxpar(hdr0, 'SCENDSEC') scut_scifile_end = MAX([fxpar(hdr0, 'SCENDSEC'), (tab2.STOP_TIME)[nframes - 1]]) nslow = FLOOR((scut_scifile_end - scut_scifile_start)/60.0D0) + 1 IF (nslow GT 0L) THEN BEGIN scut_time_slow = 60.0D0*(DINDGEN(nslow) + 1.0D0) + scut_scifile_start orbit_number = INTARR(nslow, /NOZERO) solar_zenith_angle = FLTARR(nslow, /NOZERO) sc_shadow_flag = BYTARR(nslow) subsolar_latitude = FLTARR(nslow, /NOZERO) subsolar_longitude = FLTARR(nslow, /NOZERO) subsc_latitude = FLTARR(nslow, /NOZERO) subsc_longitude = FLTARR(nslow, /NOZERO) sc_zenith_angle = FLTARR(nslow, /NOZERO) earth_phase_angle = FLTARR(nslow, /NOZERO) phase_angle = FLTARR(nslow, /NOZERO) azimuth_angle = FLTARR(nslow, /NOZERO) ra = FLTARR(nslow, /NOZERO) dec = FLTARR(nslow, /NOZERO) ecliptic_x = FLTARR(nslow, /NOZERO) ecliptic_y = FLTARR(nslow, /NOZERO) ecliptic_z = FLTARR(nslow, /NOZERO) stim_drift1 = FLTARR(nslow) stim_drift2 = FLTARR(nslow) other_inst_interf = BYTARR(nslow) global_data_qual = INTARR(nslow) heuristic_quality = INTARR(nslow) dark_factor = FLTARR(nslow) stray_factor = FLTARR(nslow) off_nadir_pointing = BYTARR(nslow) ; Auxiliary arrays stim1_data = FLTARR(!LAMP_STIM1_X_MAX - !LAMP_STIM1_X_MIN + 1, nslow) stim2_data = FLTARR(!LAMP_STIM2_X_MAX - !LAMP_STIM2_X_MIN + 1, nslow) slow_table_waveimg = FLTARR(1024, 32, nslow, /NOZERO) IF (KEYWORD_SET(aeffflag)) THEN slow_table_aeffimg = FLTARR(1024, 32, nslow, 3, 2, /NOZERO) stim1x = FINDGEN(!LAMP_STIM1_X_MAX - !LAMP_STIM1_X_MIN + 1) + !LAMP_STIM1_X_MIN stim2x = FINDGEN(!LAMP_STIM2_X_MAX - !LAMP_STIM2_X_MIN + 1) + !LAMP_STIM2_X_MIN ENDIF ELSE BEGIN ; allocate dummy arrays for FITS file scut_time_slow = DINDGEN(1, /NOZERO) orbit_number = INTARR(1, /NOZERO) solar_zenith_angle = FLTARR(1, /NOZERO) sc_shadow_flag = BYTARR(1, /NOZERO) subsolar_latitude = FLTARR(1, /NOZERO) subsolar_longitude = FLTARR(1, /NOZERO) subsc_latitude = FLTARR(1, /NOZERO) subsc_longitude = FLTARR(1, /NOZERO) sc_zenith_angle = FLTARR(1, /NOZERO) earth_phase_angle = FLTARR(1, /NOZERO) phase_angle = FLTARR(1, /NOZERO) azimuth_angle = FLTARR(1, /NOZERO) ra = FLTARR(1, /NOZERO) dec = FLTARR(1, /NOZERO) ecliptic_x = FLTARR(1, /NOZERO) ecliptic_y = FLTARR(1, /NOZERO) ecliptic_z = FLTARR(1, /NOZERO) stim_drift1 = FLTARR(1, /NOZERO) stim_drift2 = FLTARR(1, /NOZERO) other_inst_interf = BYTARR(1, /NOZERO) global_data_qual = INTARR(1, /NOZERO) heuristic_quality = INTARR(1, /NOZERO) dark_factor = FLTARR(1, /NOZERO) stray_factor = FLTARR(1, /NOZERO) off_nadir_pointing = BYTARR(1, /NOZERO) ENDELSE ; Load basic frame data into arrays j = -1L ; FOV photon event index k = -1L ; histogram index m = -1L ; countrate table entry index FOR iframe = 0L, nframes - 1 DO BEGIN IF (skip_frame[iframe]) THEN CONTINUE ; current frame was skipped above frame_data = (tab3.FRAME_DATA)[*, iframe] IF ((tab2.ACQ_MODE)[iframe] EQ !LAMP_PIXELLIST_MODE) THEN BEGIN ; pixellist mode hack_rate = (tab2.HACK_RATE)[iframe] hack_incr = 2^(hack_rate - 1) current_hack = first_hack[iframe] - hack_incr ; back off one hack increment for initial photon events scut = hack_clock * (current_hack + 0.5D0*hack_incr) + hack_offset ; offset by 1/2 hack_incr for midpoint slow_idx = FLOOR((scut - scut_scifile_start)/60.0D0) IF (slow_idx LT 0L) THEN slow_idx = 0L IF (slow_idx GT nslow - 1) THEN slow_idx = nslow - 1 totcounts_this_hack = 0L fovevents_this_hack = 0L drkcounts_this_hack = 0L straycnts_this_hack = 0L IF ((hack_rate GT 0) AND (ncrentries_per_frame[iframe] GT 0)) THEN BEGIN ncrentries_this_frame = 1L nhacks_per_cr_interval = MAX([2^(6 - hack_rate), 1]) m += 1 cr_hack_time[m] = first_hack[iframe] + nhacks_per_cr_interval * hack_incr / 2 cr_scut_time[m] = hack_clock * cr_hack_time[m] + hack_offset cr_hack_incr[m] = nhacks_per_cr_interval * hack_incr cr_hack_counter = -1 ENDIF FOR iword = 1L, 32767L DO BEGIN IF (frame_data[iword] NE 0) THEN BEGIN ; frame data word is not zero fill IF (ISHFT(frame_data[iword], -15) EQ 0) THEN BEGIN ; photon event totcounts_this_hack += 1 xtmp = ISHFT(ISHFT(frame_data[iword], 6), -6) ytmp = ISHFT(ISHFT(frame_data[iword], 1), -11) location = lamp_mike_event_locate(xtmp, ytmp) SWITCH (location) OF !LAMP_FOV_ID : BEGIN fovevents_this_hack += 1 j += 1 hack_time[j] = current_hack detector_x[j] = xtmp detector_y[j] = ytmp scut_time[j] = scut weighted_count[j] = 1.0 weighted_variance[j] = 0.0 deadtime_factor[j] = 1.0 pix_data_qual[j] = 0 hvps_state[j] = (tab2.HVPS_STATE)[iframe] apdoor_state[j] = (tab2.APD_STATE)[iframe] cr_index[j] = m BREAK END !LAMP_STIM1_ID : BEGIN stim1_data[xtmp - !LAMP_STIM1_X_MIN, slow_idx] += 1.0 BREAK END !LAMP_STIM2_ID : BEGIN stim2_data[xtmp - !LAMP_STIM2_X_MIN, slow_idx] += 1.0 BREAK END !LAMP_DARK1_ID : ; falls through !LAMP_DARK2_ID : ; falls through !LAMP_DARK3_ID : ; falls through !LAMP_DARK4_ID : BEGIN drkcounts_this_hack += 1 BREAK END !LAMP_STRAY1_ID : ; falls through !LAMP_STRAY2_ID : BEGIN straycnts_this_hack += 1 BREAK END ENDSWITCH ENDIF ELSE BEGIN ; time hack IF (KEYWORD_SET(deadflag)) THEN BEGIN count_rate = totcounts_this_hack / (hack_clock * hack_incr) ; Hz deadtime_correction = lamp_mike_deadtime(count_rate) FOR i = j - fovevents_this_hack + 1, j DO BEGIN weighted_count[i] *= deadtime_correction deadtime_factor[i] = deadtime_correction ENDFOR dark_factor[slow_idx] += deadtime_correction * drkcounts_this_hack stray_factor[slow_idx] += deadtime_correction * straycnts_this_hack ENDIF ELSE BEGIN dark_factor[slow_idx] += FLOAT(drkcounts_this_hack) stray_factor[slow_idx] += FLOAT(straycnts_this_hack) ENDELSE current_hack += hack_incr expected_rolled_hack = current_hack MOD 32768L IF (expected_rolled_hack LT 0L) THEN expected_rolled_hack += 32768L actual_rolled_hack = ISHFT(ISHFT(frame_data[iword], 1), -1) IF (actual_rolled_hack NE expected_rolled_hack) THEN BEGIN message = ' unexpected hack value (' + STRING(actual_rolled_hack) + ') ' + $ 'in word ' + STRING(iword) + ' of frame ' + STRING(iframe) + '; ' + $ 'expected value (' + STRING(expected_rolled_hack) + ')' lamp_mike_log, message, ERRNUM = 1, VERBOSELEVEL = verbose ENDIF scut = hack_clock * (current_hack + 0.5D0*hack_incr) + hack_offset ; offset by 1/2 hack_incr for midpoint slow_idx = FLOOR((scut - scut_scifile_start)/60.0D0) IF (slow_idx LT 0L) THEN slow_idx = 0L IF (slow_idx GT nslow - 1) THEN slow_idx = nslow - 1 totcounts_this_hack = 0L fovevents_this_hack = 0L drkcounts_this_hack = 0L straycnts_this_hack = 0L IF ((hack_rate GT 0) AND (ncrentries_per_frame[iframe] GT 0)) THEN BEGIN cr_hack_counter += 1 IF ((cr_hack_counter EQ nhacks_per_cr_interval) AND $ (ncrentries_this_frame LT ncrentries_per_frame[iframe])) THEN BEGIN ncrentries_this_frame += 1L m += 1 cr_hack_time[m] = current_hack + nhacks_per_cr_interval * hack_incr / 2 cr_scut_time[m] = hack_clock * cr_hack_time[m] + hack_offset cr_hack_incr[m] = nhacks_per_cr_interval * hack_incr cr_hack_counter = 0 ENDIF ENDIF ENDELSE ENDIF ENDFOR ENDIF ELSE BEGIN ; histogram mode - perform radiometric calibrations and store in histogram data cube image = FLOAT(REFORM(frame_data, 1024, 32)) exptime = (tab2.STOP_TIME)[iframe] - (tab2.START_TIME)[iframe] IF (good_housekeeping) THEN $ indices = WHERE((tab6.SCUT_TIME GE (tab2.START_TIME)[iframe]) AND $ (tab6.SCUT_TIME LE (tab2.STOP_TIME)[iframe])) ; create the wavelength image for this histogram waveimg = waveimage ; use default if STIM and HK data not available IF ((tab2.STIM_STATE)[iframe] EQ !LAMP_STIM_ON) THEN BEGIN ; stim1y = image[!LAMP_STIM1_X_MIN : !LAMP_STIM1_X_MAX, !LAMP_STIM1_Y] stim1y = image[stim1x, !LAMP_STIM1_Y] ; start of Andrew's code for stim1 stim1err = SQRT(stim1y) stim1tot = TOTAL(stim1y) stim1toterr = SQRT(stim1tot) stim1_index1 = stim1y * stim1x stim1err_index1 = stim1err * stim1x stim1indextot = totalerr(stim1_index1, stim1err_index1, err = stim1indextoterr) stimpos1 = diverr(stim1indextot, stim1indextoterr, stim1tot, stim1toterr, err = stimpos1err_i) stimpos1err = stimpos1err_i ; end of Andrew's code stim1pos = (MATRIX_MULTIPLY(stim1x, stim1y, /ATRANSPOSE) / TOTAL(stim1y))[0] ; my value ; stim2y = image[!LAMP_STIM2_X_MIN : !LAMP_STIM2_X_MAX, !LAMP_STIM2_Y] stim2y = image[stim2x, !LAMP_STIM2_Y] ; start of Andrew's code for stim2 stim2err = SQRT(stim2y) stim2tot = TOTAL(stim2y) stim2toterr = SQRT(stim2tot) stim2_index2 = stim2y * stim2x stim2err_index2 = stim2err * stim2x stim2indextot = totalerr(stim2_index2, stim2err_index2, err = stim2indextoterr) stimpos2 = diverr(stim2indextot, stim2indextoterr, stim2tot, stim2toterr, err = stimpos2err_i) stimpos2err = stimpos2err_i ; end of Andrew's code stim2pos = (MATRIX_MULTIPLY(stim2x, stim2y, /ATRANSPOSE) / TOTAL(stim2y))[0] ; my value ; stimpos = [stim1pos, stim2pos] stimpos = [stimpos1, stimpos2] ; Andrew's values stimposerr = [stimpos1err, stimpos2err] ; Andrew's values ; waveimg = lamp_mike_wavecal(STIMPOS = stimpos) ; DEK - do not recompute the wave image to go along with the ; histogram; rather, resample the histogram image to match the ; fiducial wavelength image x_corrected = lamp_mike_temp_correction(stimpos, /REVERSE) yindex = FINDGEN(32) imgtmp = INTERPOLATE(image, x_corrected, yindex, CUBIC = -0.5, /GRID) > 0.0 input_total = TOTAL(image[!LAMP_FOV_X_MIN - 50:!LAMP_FOV_X_MAX + 50, *]) output_total = TOTAL(imgtmp[!LAMP_FOV_X_MIN - 50:!LAMP_FOV_X_MAX + 50, *]) IF (input_total GT 0.0 AND output_total GT 0.0) THEN imgtmp *= input_total / output_total image = imgtmp ENDIF ELSE IF (good_housekeeping AND (N_ELEMENTS(indices) GT 0)) THEN BEGIN detetemp = (tab6.DET_HOUS_TMP)[indices] avgdetetemp = AVG(detetemp) ; waveimg = lamp_mike_wavecal(DETETEMP = avgdetetemp) ENDIF ; zero out PHD and frame header word image[!LAMP_PHD_X_MIN : !LAMP_PHD_X_MAX, !LAMP_PHD_Y_MIN : !LAMP_PHD_Y_MAX] = 0.0 ; PHD image[0, 0] = 0.0 ; frame header word ; get countrate count_rate = TOTAL(image) / exptime ; zero out STIM pixels image[!LAMP_STIM1_X_MIN : !LAMP_STIM1_X_MAX, !LAMP_STIM1_Y] = 0.0 image[!LAMP_STIM2_X_MIN : !LAMP_STIM2_X_MAX, !LAMP_STIM2_Y] = 0.0 ; set saturated values to NaN saturated = WHERE(image EQ 65535, nsaturated) IF (nsaturated GT 0) THEN image[saturated] = !VALUES.F_NAN ; create error image; use gaussian approximation where image > 2500 counts, else poissonian error_image = FLTARR(SIZE(image, /DIMENSIONS), /NOZERO) over2500 = WHERE(image GT 2500.0, nover2500, COMPLEMENT = under2500, NCOMPLEMENT = nunder2500) IF (nover2500 GT 0) THEN error_image[over2500] = SQRT(image[over2500]) IF (nunder2500 GT 0) THEN error_image[under2500] = poisson_errors[image[under2500]] ; divide by exposure time image /= exptime error_image /= exptime ; dead time correction IF (KEYWORD_SET(deadflag)) THEN BEGIN deadtime_correction = lamp_mike_deadtime(count_rate) image *= deadtime_correction error_image *= deadtime_correction message = ' dead time correction = YES' lamp_mike_log, message ENDIF ; dark subtraction IF (KEYWORD_SET(darkflag)) THEN BEGIN image = adderr(image, error_image, darkimage, darkerror, ERR = error_image, /SUBTRACT) message = ' dark subtraction = YES' lamp_mike_log, message ENDIF ; flatfield correction IF (KEYWORD_SET(flatflag)) THEN BEGIN CASE ((tab2.HVPS_STATE)[iframe]) OF !LAMP_HVPS_A_ON : IF ((tab2.APD_STATE)[iframe] EQ !LAMP_APDOOR_OPEN) THEN $ flatimage = flatAO $ ELSE $ flatimage = flatAC !LAMP_HVPS_B_ON : IF ((tab2.APD_STATE)[iframe] EQ !LAMP_APDOOR_OPEN) THEN $ flatimage = flatBO $ ELSE $ flatimage = flatBC !LAMP_HVPS_AB_ON : IF ((tab2.APD_STATE)[iframe] EQ !LAMP_APDOOR_OPEN) THEN $ flatimage = flatABO $ ELSE $ flatimage = flatABC ENDCASE image /= flatimage error_image /= flatimage message = ' flatfield correction = YES' lamp_mike_log, message ENDIF ; effective area calibration IF (KEYWORD_SET(aeffflag)) THEN BEGIN CASE ((tab2.HVPS_STATE)[iframe]) OF !LAMP_HVPS_A_ON : IF ((tab2.APD_STATE)[iframe] EQ !LAMP_APDOOR_OPEN) THEN $ aeffimage = INTERPOL(aeffAO, aeffwave, waveimg) $ ELSE $ aeffimage = INTERPOL(aeffAC, aeffwave, waveimg) !LAMP_HVPS_B_ON : IF ((tab2.APD_STATE)[iframe] EQ !LAMP_APDOOR_OPEN) THEN $ aeffimage = INTERPOL(aeffBO, aeffwave, waveimg) $ ELSE $ aeffimage = INTERPOL(aeffBC, aeffwave, waveimg) !LAMP_HVPS_AB_ON : IF ((tab2.APD_STATE)[iframe] EQ !LAMP_APDOOR_OPEN) THEN $ aeffimage = INTERPOL(aeffABO, aeffwave, waveimg) $ ELSE $ aeffimage = INTERPOL(aeffABC, aeffwave, waveimg) ENDCASE image /= aeffimage error_image /= aeffimage message = ' effective area calibration = YES' lamp_mike_log, message ENDIF ; zero entries outside LAMP spectral range and set error = 1 out_of_range = WHERE((waveimg LT !LAMP_WAVELENGTH_LO) OR $ (waveimg GT !LAMP_WAVELENGTH_HI), count) IF (count GT 0) THEN BEGIN image[out_of_range] = 0.0 error_image[out_of_range] = 1.0 ENDIF ; store results in histogram cube k += 1 histogram_cube[*, *, 0, k] = image histogram_cube[*, *, 1, k] = error_image ENDELSE ENDFOR ; Divide dark_factor array by (60.0 * n_dark_elements) to get average counts/sec/element ; in our dark signal measuring regions dark_factor /= (60.0 * n_dark_elements) ; Divide stray_factor array by (60.0 * n_stray_elements) to get average counts/sec/element ; in our stray signal measuring regions stray_factor /= (60.0 * n_stray_elements) ; Compute wavelength and effective area calibration images for pixellist entries, ; but at the slow table cadence FOR i = 0L, nslow - 1L DO BEGIN waveimg = waveimage ; use default if STIM and HK data not available IF ((MAX(stim1_data[*, i]) GT 0.0) AND (MAX(stim2_data[*, i]) GT 0.0)) THEN BEGIN ; start of Andrew's code for stim1 stim1err = SQRT(stim1_data[*, i]) stim1tot = TOTAL(stim1_data[*, i]) stim1toterr = SQRT(stim1tot) stim1_index1 = stim1_data[*, i] * stim1x stim1err_index1 = stim1err * stim1x stim1indextot = totalerr(stim1_index1, stim1err_index1, err = stim1indextoterr) stimpos1 = diverr(stim1indextot, stim1indextoterr, stim1tot, stim1toterr, err = stimpos1err_i) stimpos1err = stimpos1err_i ; end of Andrew's code stim1pos = (MATRIX_MULTIPLY(stim1x, stim1_data[*, i], /ATRANSPOSE) / TOTAL(stim1_data[*, i]))[0] ; my value ; start of Andrew's code for stim2 stim2err = SQRT(stim2_data[*, i]) stim2tot = TOTAL(stim2_data[*, i]) stim2toterr = SQRT(stim2tot) stim2_index2 = stim2_data[*, i] * stim2x stim2err_index2 = stim2err * stim2x stim2indextot = totalerr(stim2_index2, stim2err_index2, err = stim2indextoterr) stimpos2 = diverr(stim2indextot, stim2indextoterr, stim2tot, stim2toterr, err = stimpos2err_i) stimpos2err = stimpos2err_i ; end of Andrew's code stim2pos = (MATRIX_MULTIPLY(stim2x, stim2_data[*, i], /ATRANSPOSE) / TOTAL(stim2_data[*, i]))[0] ; my value ; stim_drift1[i] = stim1pos ; stim_drift2[i] = stim2pos stim_drift1[i] = stimpos1 ; Andrew's value stim_drift2[i] = stimpos2 ; Andrew's value stimpos = [stim_drift1[i], stim_drift2[i]] stimposerr = [stimpos1err, stimpos2err] ; Andrew's values waveimg = lamp_mike_wavecal(STIMPOS = stimpos) ENDIF ELSE IF (good_housekeeping) THEN BEGIN result = MIN(ABS(tab6.SCUT_TIME - scut_time_slow[i]), idx) IF (result LE 1.0D0) THEN BEGIN ; HK entry within a second of desired time detetemp = (tab6.DET_HOUS_TMP)[idx] waveimg = lamp_mike_wavecal(DETETEMP = detetemp) ENDIF ENDIF slow_table_waveimg[*, *, i] = waveimg IF (KEYWORD_SET(aeffflag)) THEN BEGIN slow_table_aeffimg[*, *, i, 0, 0] = INTERPOL(aeffAO, aeffwave, waveimg) slow_table_aeffimg[*, *, i, 1, 0] = INTERPOL(aeffBO, aeffwave, waveimg) slow_table_aeffimg[*, *, i, 2, 0] = INTERPOL(aeffABO, aeffwave, waveimg) slow_table_aeffimg[*, *, i, 0, 1] = INTERPOL(aeffAC, aeffwave, waveimg) slow_table_aeffimg[*, *, i, 1, 1] = INTERPOL(aeffBC, aeffwave, waveimg) slow_table_aeffimg[*, *, i, 2, 1] = INTERPOL(aeffABC, aeffwave, waveimg) ENDIF ENDFOR ; Scan pixellist array and perform remaining radiometric calibrations FOR j = 0L, nfovevents - 1L DO BEGIN slow_idx = FLOOR((scut_time[j] - scut_scifile_start)/60.0D0) IF (slow_idx LT 0L) THEN slow_idx = 0L IF (slow_idx GT nslow - 1) THEN slow_idx = nslow - 1 ; assign spatial_row value, performing geometric distortion correction if requested spatial_row[j] = detector_y[j] ;; IF (KEYWORD_SET(rectflag)) THEN BEGIN ;; ENDIF ; wavelength wavelength[j] = slow_table_waveimg[detector_x[j], detector_y[j], slow_idx] ; flatfield correction IF (KEYWORD_SET(flatflag)) THEN BEGIN CASE (hvps_state[j]) OF !LAMP_HVPS_A_ON : IF (apdoor_state[j] EQ !LAMP_APDOOR_OPEN) THEN $ weighted_count[j] /= flatAO[detector_x[j], detector_y[j]] $ ELSE $ weighted_count[j] /= flatAC[detector_x[j], detector_y[j]] !LAMP_HVPS_B_ON : IF (apdoor_state[j] EQ !LAMP_APDOOR_OPEN) THEN $ weighted_count[j] /= flatBO[detector_x[j], detector_y[j]] $ ELSE $ weighted_count[j] /= flatBC[detector_x[j], detector_y[j]] !LAMP_HVPS_AB_ON : IF (apdoor_state[j] EQ !LAMP_APDOOR_OPEN) THEN $ weighted_count[j] /= flatABO[detector_x[j], detector_y[j]] $ ELSE $ weighted_count[j] /= flatABC[detector_x[j], detector_y[j]] ENDCASE ENDIF ; add to reduced countrate totals IF ((detector_x[j] GE !LAMP_FOV_X_MIN_CR) AND (detector_x[j] LE !LAMP_FOV_X_MAX_CR) AND (ncrentries GT 0L)) THEN $ cr_countrate[cr_index[j]] += weighted_count[j] ; effective area calibration IF (KEYWORD_SET(aeffflag)) THEN BEGIN CASE (hvps_state[j]) OF !LAMP_HVPS_A_ON : IF (apdoor_state[j] EQ !LAMP_APDOOR_OPEN) THEN $ weighted_count[j] /= slow_table_aeffimg[detector_x[j], detector_y[j], slow_idx, 0, 0] $ ELSE $ weighted_count[j] /= slow_table_aeffimg[detector_x[j], detector_y[j], slow_idx, 0, 1] !LAMP_HVPS_B_ON : IF (apdoor_state[j] EQ !LAMP_APDOOR_OPEN) THEN $ weighted_count[j] /= slow_table_aeffimg[detector_x[j], detector_y[j], slow_idx, 1, 0] $ ELSE $ weighted_count[j] /= slow_table_aeffimg[detector_x[j], detector_y[j], slow_idx, 1, 1] !LAMP_HVPS_AB_ON : IF (apdoor_state[j] EQ !LAMP_APDOOR_OPEN) THEN $ weighted_count[j] /= slow_table_aeffimg[detector_x[j], detector_y[j], slow_idx, 2, 0] $ ELSE $ weighted_count[j] /= slow_table_aeffimg[detector_x[j], detector_y[j], slow_idx, 2, 1] ENDCASE ENDIF ENDFOR ; Compute reduced count rates IF (ncrentries GT 0L) THEN BEGIN cr_dark_factor = INTERPOL(dark_factor, scut_time_slow, cr_scut_time) n_elements_fov_cr = (!LAMP_FOV_X_MAX_CR - !LAMP_FOV_X_MIN_CR + 1) * (!LAMP_FOV_Y_MAX - !LAMP_FOV_Y_MIN + 1) cr_dark_factor *= n_elements_fov_cr cr_countrate /= (hack_clock * cr_hack_incr) cr_countrate -= cr_dark_factor ; FOR m = 0L, ncrentries - 1L DO BEGIN ; cr_countrate[m] /= (hack_clock * cr_hack_incr[m]) ; cr_countrate[m] -= cr_dark_factor[m] ; ENDFOR ENDIF ; Scan pixellist array and perform SPICE-based geometric calibrations (vectorized method) IF (nfovevents GT 0L) THEN BEGIN FOR y = !LAMP_FOV_Y_MIN, !LAMP_FOV_Y_MAX DO BEGIN ; Get array indices for this value of spatial row index1 = WHERE(spatial_row EQ y, count1) IF (count1 NE 0) THEN BEGIN ; Get ET from SCUT scut_seconds = FLOOR(scut_time[index1]) scut_fractions = ROUND(65536.0D0 * (scut_time[index1] - scut_seconds)) cspice_scs2e, !SPICE_LRO_ID, STRCOMPRESS(STRING(scut_seconds) + '.' + STRING(scut_fractions), /REMOVE_ALL), $ et ; Get altitude cspice_subpt, 'NEARPOINT', !SPICE_TARGET, et, !SPICE_ABCORR, !SPICE_OBSRVR, $ spoint, alt altitude[index1] = alt ; Allocate arrays for longitude and latitude values and illumination angles dlon = DBLARR(count1, /NOZERO) dlat = DBLARR(count1, /NOZERO) dinc = DBLARR(count1, /NOZERO) demi = DBLARR(count1, /NOZERO) ; Get latitude and longitude of pixel center, plus incidence and emission angles dvec = pixel_centers[*, y] ;; DEK: SPICE version N0063 dies on the call below with a system memory error; had to revert to N0062 to get it to work cspice_srfxpt, !SPICE_METHOD[0], !SPICE_TARGET, et, !SPICE_ABCORR, !SPICE_OBSRVR, !SPICE_SC_BUS_FRAME, dvec, $ spoint, dist, trgepc, obspos, found index2 = WHERE(found, count2, COMPLEMENT = index3, NCOMPLEMENT = count3) IF (count2 NE 0) THEN BEGIN cspice_reclat, spoint[*, index2], radius, lon, lat dlon[index2] = lon * !SPICE_DPR index4 = WHERE(dlon LT 0.0D0, count4) IF (count4 NE 0) THEN dlon[index4] += 360.0D0 dlat[index2] = lat * !SPICE_DPR cspice_illum, !SPICE_TARGET, et[index2], !SPICE_ABCORR, !SPICE_OBSRVR, spoint[*, index2], $ phase, solar, emissn dinc[index2] = solar * !SPICE_DPR demi[index2] = emissn * !SPICE_DPR ENDIF IF (count3 NE 0) THEN BEGIN dlon[index3] = !SPICE_BADANGLE dlat[index3] = !SPICE_BADANGLE dinc[index3] = !SPICE_BADANGLE demi[index3] = !SPICE_BADANGLE ENDIF longitude[index1] = dlon latitude[index1] = dlat incidence_angle[index1] = dinc emission_angle[index1] = demi ; Get latitude and longitude for lower left pixel corner dvec = pixel_corners[*, 64 - y] cspice_srfxpt, !SPICE_METHOD[0], !SPICE_TARGET, et, !SPICE_ABCORR, !SPICE_OBSRVR, !SPICE_SC_BUS_FRAME, dvec, $ spoint, dist, trgepc, obspos, found index2 = WHERE(found, count2, COMPLEMENT = index3, NCOMPLEMENT = count3) IF (count2 NE 0) THEN BEGIN cspice_reclat, spoint[*, index2], radius, lon, lat dlon[index2] = lon * !SPICE_DPR index4 = WHERE(dlon LT 0.0D0, count4) IF (count4 NE 0) THEN dlon[index4] += 360.0D0 dlat[index2] = lat * !SPICE_DPR ENDIF IF (count3 NE 0) THEN BEGIN dlon[index3] = !SPICE_BADANGLE dlat[index3] = !SPICE_BADANGLE ENDIF longitude_ll[index1] = dlon latitude_ll[index1] = dlat ; Get latitude and longitude for upper left pixel corner dvec = pixel_corners[*, 65 - y] cspice_srfxpt, !SPICE_METHOD[0], !SPICE_TARGET, et, !SPICE_ABCORR, !SPICE_OBSRVR, !SPICE_SC_BUS_FRAME, dvec, $ spoint, dist, trgepc, obspos, found index2 = WHERE(found, count2, COMPLEMENT = index3, NCOMPLEMENT = count3) IF (count2 NE 0) THEN BEGIN cspice_reclat, spoint[*, index2], radius, lon, lat dlon[index2] = lon * !SPICE_DPR index4 = WHERE(dlon LT 0.0D0, count4) IF (count4 NE 0) THEN dlon[index4] += 360.0D0 dlat[index2] = lat * !SPICE_DPR ENDIF IF (count3 NE 0) THEN BEGIN dlon[index3] = !SPICE_BADANGLE dlat[index3] = !SPICE_BADANGLE ENDIF longitude_ul[index1] = dlon latitude_ul[index1] = dlat ; Get latitude and longitude for upper right pixel corner dvec = pixel_corners[*, y] cspice_srfxpt, !SPICE_METHOD[0], !SPICE_TARGET, et, !SPICE_ABCORR, !SPICE_OBSRVR, !SPICE_SC_BUS_FRAME, dvec, $ spoint, dist, trgepc, obspos, found index2 = WHERE(found, count2, COMPLEMENT = index3, NCOMPLEMENT = count3) IF (count2 NE 0) THEN BEGIN cspice_reclat, spoint[*, index2], radius, lon, lat dlon[index2] = lon * !SPICE_DPR index4 = WHERE(dlon LT 0.0D0, count4) IF (count4 NE 0) THEN dlon[index4] += 360.0D0 dlat[index2] = lat * !SPICE_DPR ENDIF IF (count3 NE 0) THEN BEGIN dlon[index3] = !SPICE_BADANGLE dlat[index3] = !SPICE_BADANGLE ENDIF longitude_ur[index1] = dlon latitude_ur[index1] = dlat ; Get latitude and longitude for lower right pixel corner dvec = pixel_corners[*, y + 1] cspice_srfxpt, !SPICE_METHOD[0], !SPICE_TARGET, et, !SPICE_ABCORR, !SPICE_OBSRVR, !SPICE_SC_BUS_FRAME, dvec, $ spoint, dist, trgepc, obspos, found index2 = WHERE(found, count2, COMPLEMENT = index3, NCOMPLEMENT = count3) IF (count2 NE 0) THEN BEGIN cspice_reclat, spoint[*, index2], radius, lon, lat dlon[index2] = lon * !SPICE_DPR index4 = WHERE(dlon LT 0.0D0, count4) IF (count4 NE 0) THEN dlon[index4] += 360.0D0 dlat[index2] = lat * !SPICE_DPR ENDIF IF (count3 NE 0) THEN BEGIN dlon[index3] = !SPICE_BADANGLE dlat[index3] = !SPICE_BADANGLE ENDIF longitude_lr[index1] = dlon latitude_lr[index1] = dlat ENDIF ; count1 NE 0 ENDFOR ; y = !LAMP_FOV_Y_MIN, !LAMP_FOV_Y_MAX ENDIF ; (nfovevents GT 0L) ; Compute and store remaining slow pixellist table entries FOR i = 0L, nslow - 1L DO BEGIN ; Get ET from SCUT scut_seconds = FLOOR(scut_time_slow[i]) scut_fractions = ROUND(65536.0D0 * (scut_time_slow[i] - scut_seconds)) cspice_scs2e, !SPICE_LRO_ID, STRCOMPRESS(STRING(scut_seconds, '.', scut_fractions), /REMOVE_ALL), et ; Get lunar surface crossing point of the LAMP boresight cspice_srfxpt, !SPICE_METHOD[0], !SPICE_TARGET, et, !SPICE_ABCORR, !SPICE_OBSRVR, !SPICE_SC_BUS_FRAME, bsight, $ spoint, dist, trgepc, obspos, found ; From this compute solar and S/C zenith angles and phase IF (found) THEN BEGIN cspice_illum, !SPICE_TARGET, et, !SPICE_ABCORR, !SPICE_OBSRVR, spoint, $ phase, solar, emissn azimuth_angle[i] = ACOS((COS(phase) - COS(solar)*COS(emissn))/(SIN(solar)*SIN(emissn))) * !SPICE_DPR ; azimuth angle computed via spherical trig identity, range = [0, PI] solar_zenith_angle[i] = solar * !SPICE_DPR sc_zenith_angle[i] = emissn * !SPICE_DPR phase_angle[i] = phase * !SPICE_DPR ENDIF ELSE BEGIN ; insert bad values to flag failure azimuth_angle[i] = !SPICE_BADANGLE solar_zenith_angle[i] = !SPICE_BADANGLE sc_zenith_angle[i] = !SPICE_BADANGLE phase_angle[i] = !SPICE_BADANGLE ENDELSE ; Get equatorial inertial frame coordinates of LAMP boresight and transform to RA and DEC cspice_pxform, !SPICE_SC_BUS_FRAME, !SPICE_EQUATORIAL_FRAME, et, rotate cspice_mxv, rotate, bsight, bsight_inertial cspice_recrad, bsight_inertial, radius, ra_tmp, dec_tmp ra[i] = ra_tmp * !SPICE_DPR dec[i] = dec_tmp * !SPICE_DPR ; Get the ecliptic coordinates of LRO wrt the solar system barycenter cspice_spkssb, !SPICE_LRO_ID, et, !SPICE_ECLIPTIC_FRAME, starg ; starg also returns velocity components ecliptic_x[i] = starg[0] ecliptic_y[i] = starg[1] ecliptic_z[i] = starg[2] ; Get the direction in the equatorial inertial frame from LRO to the Sun cspice_spkezp, !SPICE_SUN_ID, et, !SPICE_EQUATORIAL_FRAME, !SPICE_ABCORR, !SPICE_LRO_ID, $ ptarg, ltime cspice_srfxpt, !SPICE_METHOD[0], !SPICE_TARGET, et, !SPICE_ABCORR, !SPICE_OBSRVR, !SPICE_EQUATORIAL_FRAME, ptarg, $ spoint, dist, trgepc, obspos, found IF (found) THEN sc_shadow_flag[i] = 1 ; vector from LRO to Sun intersects the Moon, thus LRO is in shadow ; NOTE: this algorithm does NOT account for a finite-size Sun, so there is some inaccuracy near sunrise/sunset as seen from LRO ; Get sub-spacecraft latitude and longitude cspice_subpnt, !SPICE_METHOD[1], !SPICE_TARGET, et, !SPICE_MOON_FRAME, !SPICE_ABCORR, !SPICE_OBSRVR, $ spoint, trgepc, srfvec cspice_reclat, spoint, radius, lon, lat lon = lon * !SPICE_DPR lat = lat * !SPICE_DPR IF (lon LT 0.0D0) THEN lon = lon + 360.0D0 subsc_longitude[i] = lon subsc_latitude[i] = lat ; Get the boresight vector in the Moon frame and compute the angle between boresight and nadir cspice_pxform, !SPICE_SC_BUS_FRAME, !SPICE_MOON_FRAME, et, rotate cspice_mxv, rotate, bsight, bsight_moonfixed separation = cspice_vsep(srfvec, bsight_moonfixed) * !SPICE_DPR IF (separation GT !LAMP_NADIR_PTG_THRESH) THEN off_nadir_pointing[i] = 1 ; Get sub-solar latitude and longitude cspice_subslr, !SPICE_METHOD[1], !SPICE_TARGET, et, !SPICE_MOON_FRAME, !SPICE_ABCORR, !SPICE_OBSRVR, $ spoint, trgepc, srfvec cspice_reclat, spoint, radius, lon, lat lon = lon * !SPICE_DPR lat = lat * !SPICE_DPR IF (lon LT 0.0D0) THEN lon = lon + 360.0D0 subsolar_longitude[i] = lon subsolar_latitude[i] = lat ; Get Earth phase angle cspice_subpnt, !SPICE_METHOD[1], STRING(!SPICE_EARTH_ID), et, !SPICE_EARTH_FRAME, !SPICE_ABCORR, !SPICE_OBSRVR, $ spoint, trgepc, srfvec cspice_illum, STRING(!SPICE_EARTH_ID), et, !SPICE_ABCORR, !SPICE_OBSRVR, spoint, $ phase, solar, emissn earth_phase_angle[i] = phase * !SPICE_DPR ; Get LRO orbit number orbit_number[i] = -1 IF ((scut_time_slow[i] GE lro_orbit_scut[0]) && (scut_time_slow[i] LE lro_orbit_scut[n_lro_orbits-1])) THEN BEGIN IF (scut_time_slow[i] EQ lro_orbit_scut[n_lro_orbits-1]) THEN orbit_number[i] = lro_orbit_number[n_lro_orbits-1] ELSE BEGIN jhi = 1L WHILE (scut_time_slow[i] GE lro_orbit_scut[jhi]) DO jhi = jhi + 1L jlo = jhi - 1L IF ((scut_time_slow[i] LT lro_orbit_scut[jhi]) && (scut_time_slow[i] GE lro_orbit_scut[jlo]) && $ (lro_orbit_number[jhi] EQ (lro_orbit_number[jlo] + 1))) THEN orbit_number[i] = lro_orbit_number[jlo] ENDELSE ENDIF ENDFOR ; i = 0L, nslow - 1L ; Prepare data and headers for writing to RDR FITS ; Calibrated door open histogram, primary extension, image im0_rdr = FLTARR(1024, 32) im1_rdr = im0_rdr ; generate image for door closed histogram at same time IF (nfovevents GT 0L) THEN BEGIN ; use pixellist data preferentially if available FOR j = 0L, nfovevents - 1L DO BEGIN IF (apdoor_state[j] EQ !LAMP_APDOOR_OPEN) THEN $ im0_rdr[detector_x[j], spatial_row[j]] += weighted_count[j] $ ELSE IF (apdoor_state[j] EQ !LAMP_APDOOR_CLOSED) THEN $ im1_rdr[detector_x[j], spatial_row[j]] += weighted_count[j] ENDFOR ENDIF ELSE IF (nhistograms GT 0L) THEN BEGIN ; use histograms next if available im0_rdr = histogram_cube[*, *, 0, 0] IF (nhistograms GT 1L) THEN im1_hdr = histogram_cube[*, *, 0, 1] ENDIF ELSE BEGIN ; reuse EDR images as a last resort when no other data are present im0_rdr = im0 im1_rdr = im1 ENDELSE ; reuse EDR header with the following modifications sxaddpar, hdr0, 'BITPIX' , -32 sxaddpar, hdr0, 'COMMENT' , 'LAMP Level 2 (CODMAC Level 3) processing information ------------', AFTER = 'DOHVPS' sxaddpar, hdr0, 'SL2PROC' , !MIKE_SL2PROC, 'SOC pipeline level 2 software' sxaddpar, hdr0, 'SL2VER' , !MIKE_SL2VER, 'SOC pipeline level 2 software version' sxaddpar, hdr0, 'SL2DATE' , !MIKE_SL2DATE, 'SOC pipeline level 2 software date' sxaddpar, hdr0, 'FILE_IN' , FILE_BASENAME(infiles[ifile]), 'Input file for processing' sxaddpar, hdr0, 'FILE_OUT', FILE_BASENAME(outfiles[ifile]), 'Output file after processing' sxaddpar, hdr0, 'DIR_CAL' , calibration_dir, 'Directory for calibration data' sxaddpar, hdr0, 'BADFILE' , KEYWORD_SET(badfile) ? FILE_BASENAME(badfile) : '', 'File containing bad pixel mask' sxaddpar, hdr0, 'BADFLAG' , KEYWORD_SET(badflag) ? 'T' : 'F', 'Bad pixel mask applied to data' sxaddpar, hdr0, 'BADVALUE', 'IEEE NaN', 'Bad pixel value' sxaddpar, hdr0, 'DEADFLAG', KEYWORD_SET(deadflag) ? 'T' : 'F', 'Dead-time correction applied' sxaddpar, hdr0, 'DARKFILE', KEYWORD_SET(darkfile) ? FILE_BASENAME(darkfile) : '', 'Detector dark count rate file' sxaddpar, hdr0, 'DARKFLAG', KEYWORD_SET(darkflag) ? 'T' : 'F', 'Dark count correction applied' sxaddpar, hdr0, 'FLATFILE', KEYWORD_SET(flatfile) ? FILE_BASENAME(flatfile) : '', 'Detector flatfield file' sxaddpar, hdr0, 'FLATFLAG', KEYWORD_SET(flatflag) ? 'T' : 'F', 'Flatfield correction applied' sxaddpar, hdr0, 'WAVEPRO' , 'lamp_mike_wavecal.pro', 'IDL wavelength calibration procedure' sxaddpar, hdr0, 'WAVEFLAG', 'T', 'Wavelength calibration performed' sxaddpar, hdr0, 'AEFFFILE', KEYWORD_SET(aefffile) ? FILE_BASENAME(aefffile) : '', 'Detector effective area file' sxaddpar, hdr0, 'AEFFFLAG', KEYWORD_SET(aeffflag) ? 'T' : 'F', 'Effective area correction applied' sxaddpar, hdr0, 'LOG_FILE', KEYWORD_SET(logfile) ? FILE_BASENAME(logfile) : '', 'LAMP-MIKE processing logfile' sxaddpar, hdr0, 'MIKE_ERR', error_status, 'Exit status for LAMP-MIKE' sxaddpar, hdr0, 'COMMENT' , 'SPICE files used in geometry calculations -----------------------', AFTER = 'MIKE_ERR' sxaddpar, hdr0, 'NUMBSPCK', nkernels, 'Number of loaded SPICE kernels' FOR k = 0, nkernels - 1 DO BEGIN strkp1 = STRING(k + 1, FORMAT = '(I04)') sxaddpar, hdr0, 'SPCK' + strkp1, kernel_list[k], 'SPICE kernel ' + strkp1 ENDFOR ; Calibrated door closed histogram, extension 1, image ; image generated above with primary extension image ; reuse EDR header with the following modifications sxaddpar, hdr1, 'BITPIX', -32 sxaddpar, hdr1, 'EXTNAME', 'Calibrated Spectral Image Door Closed' ; Acquisition list, extension 2, ASCII table ; no modifications of EDR required, straight copy ; Calibrated pixel list mode data (fast), extension 3, binary table ; no modifications of EDR, this is a new RDR extension ; just write out new bintable below ; Calibrated pixel list mode data (slow), extension 4, binary table ; no modifications of EDR, this is a new RDR extension ; just write out new bintable below ; Calibrated histogram mode data, extension 5, image (stacked image cube) mkhdr, hdr5_rdr, '', /IMAGE sxaddpar, hdr5_rdr, 'XTENSION', 'IMAGE', ' Extension 5: Calibrated Histogram Mode Data' sxaddpar, hdr5_rdr, 'BITPIX', -32 sxaddpar, hdr5_rdr, 'BZERO', 0 sxaddpar, hdr5_rdr, 'BSCALE', 1 sxaddpar, hdr5_rdr, 'EXTNAME', 'Calibrated Histogram Mode Data' sxaddpar, hdr5_rdr, 'EXTVER', 1, ' Extension version number' ; Calibrated calculated countrate, extension 6, binary table ; IF (ncrentries GT 0L) THEN BEGIN ; hdr6_rdr = hdr4 ; sxaddpar, hdr6_rdr, 'EXTNAME', 'Calibrated Calculated Countrate' ; sxaddpar, hdr6_rdr, 'NAXIS2', ncrentries ; ENDIF ; High-rate LTS data, extension 7, binary table ; no modifications of EDR required, straight copy ; Housekeeping table, extension 8, binary table ; no modifications of EDR required, straight copy ; Wavelength lookup image, extension 9, image mkhdr, hdr9_rdr, waveimage, /IMAGE sxaddpar, hdr9_rdr, 'XTENSION', 'IMAGE', ' Extension 9: Wavelength Lookup Image' sxaddpar, hdr9_rdr, 'EXTNAME' , 'Wavelength Lookup Image' sxaddpar, hdr9_rdr, 'EXTVER' , 1, ' Extension version number' ; Write the RDR FITS mwrfits, im0_rdr, outfiles[ifile], hdr0, /CREATE ; calibrated door open histogram, primary extension, image mwrfits, im1_rdr, outfiles[ifile], hdr1 ; calibrated door closed histogram, extension 1, image mwrfits, tab2 , outfiles[ifile], hdr2, /ASCII ; acquisition list, extension 2, ASCII table ; calibrated pixel list mode data (fast), extension 3, binary table fxbhmake, hdr3_rdr, nfovevents, 'Calibrated Pixel List Mode Data', EXTVER = 1 fxbaddcol, 1 , hdr3_rdr, hack_time[0] , 'HACK_TIME' fxbaddcol, 2 , hdr3_rdr, detector_x[0] , 'DETECTOR_X' fxbaddcol, 3 , hdr3_rdr, detector_y[0] , 'DETECTOR_Y' fxbaddcol, 4 , hdr3_rdr, spatial_row[0] , 'SPATIAL_ROW' fxbaddcol, 5 , hdr3_rdr, wavelength[0] , 'WAVELENGTH' fxbaddcol, 6 , hdr3_rdr, scut_time[0] , 'SCUT_TIME' fxbaddcol, 7 , hdr3_rdr, weighted_count[0] , 'WEIGHTED_COUNT' fxbaddcol, 8 , hdr3_rdr, weighted_variance[0], 'WEIGHTED_VARIANCE' fxbaddcol, 9 , hdr3_rdr, latitude[0] , 'LATITUDE' fxbaddcol, 10, hdr3_rdr, longitude[0] , 'LONGITUDE' fxbaddcol, 11, hdr3_rdr, latitude_ll[0] , 'LATITUDE_LL' fxbaddcol, 12, hdr3_rdr, longitude_ll[0] , 'LONGITUDE_LL' fxbaddcol, 13, hdr3_rdr, latitude_ul[0] , 'LATITUDE_UL' fxbaddcol, 14, hdr3_rdr, longitude_ul[0] , 'LONGITUDE_UL' fxbaddcol, 15, hdr3_rdr, latitude_ur[0] , 'LATITUDE_UR' fxbaddcol, 16, hdr3_rdr, longitude_ur[0] , 'LONGITUDE_UR' fxbaddcol, 17, hdr3_rdr, latitude_lr[0] , 'LATITUDE_LR' fxbaddcol, 18, hdr3_rdr, longitude_lr[0] , 'LONGITUDE_LR' fxbaddcol, 19, hdr3_rdr, altitude[0] , 'ALTITUDE' fxbaddcol, 20, hdr3_rdr, incidence_angle[0] , 'INCIDENCE_ANGLE' fxbaddcol, 21, hdr3_rdr, emission_angle[0] , 'EMISSION_ANGLE' fxbaddcol, 22, hdr3_rdr, deadtime_factor[0] , 'DEADTIME_FACTOR' fxbaddcol, 23, hdr3_rdr, pix_data_qual[0] , 'PIX_DATA_QUAL' fxbcreate, lun, outfiles[ifile], hdr3_rdr IF (nfovevents GT 0L) THEN BEGIN column_names = ['HACK_TIME','DETECTOR_X','DETECTOR_Y','SPATIAL_ROW','WAVELENGTH','SCUT_TIME','WEIGHTED_COUNT', $ 'WEIGHTED_VARIANCE','LATITUDE','LONGITUDE','LATITUDE_LL','LONGITUDE_LL','LATITUDE_UL','LONGITUDE_UL',$ 'LATITUDE_UR','LONGITUDE_UR','LATITUDE_LR','LONGITUDE_LR','ALTITUDE','INCIDENCE_ANGLE', $ 'EMISSION_ANGLE','DEADTIME_FACTOR','PIX_DATA_QUAL'] fxbwritm, lun, column_names, hack_time, detector_x, detector_y, spatial_row, wavelength, scut_time, weighted_count, $ weighted_variance, latitude, longitude, latitude_ll, longitude_ll, latitude_ul, $ longitude_ul, latitude_ur, longitude_ur, latitude_lr, longitude_lr, altitude, $ incidence_angle, emission_angle, deadtime_factor, pix_data_qual, BUFFERSIZE=1048576L ENDIF fxbfinish, lun ; calibrated pixel list mode data (slow), extension 4, binary table fxbhmake, hdr4_rdr, nslow, 'Ancillary Data', EXTVER = 1 fxbaddcol, 1 , hdr4_rdr, scut_time_slow[0] , 'SCUT_TIME_SLOW' fxbaddcol, 2 , hdr4_rdr, orbit_number[0] , 'ORBIT_NUMBER' fxbaddcol, 3 , hdr4_rdr, solar_zenith_angle[0], 'SOLAR_ZENITH_ANGLE' fxbaddcol, 4 , hdr4_rdr, sc_shadow_flag[0] , 'SC_SHADOW_FLAG' fxbaddcol, 5 , hdr4_rdr, subsolar_latitude[0] , 'SUBSOLAR_LATITUDE' fxbaddcol, 6 , hdr4_rdr, subsolar_longitude[0], 'SUBSOLAR_LONGITUDE' fxbaddcol, 7 , hdr4_rdr, subsc_latitude[0] , 'SUBSC_LATITUDE' fxbaddcol, 8 , hdr4_rdr, subsc_longitude[0] , 'SUBSC_LONGITUDE' fxbaddcol, 9 , hdr4_rdr, sc_zenith_angle[0] , 'SC_ZENITH_ANGLE' fxbaddcol, 10, hdr4_rdr, earth_phase_angle[0] , 'EARTH_PHASE_ANGLE' fxbaddcol, 11, hdr4_rdr, phase_angle[0] , 'PHASE_ANGLE' fxbaddcol, 12, hdr4_rdr, azimuth_angle[0] , 'AZIMUTH_ANGLE' fxbaddcol, 13, hdr4_rdr, ra[0] , 'RA' fxbaddcol, 14, hdr4_rdr, dec[0] , 'DEC' fxbaddcol, 15, hdr4_rdr, ecliptic_x[0] , 'ECLIPTIC_X' fxbaddcol, 16, hdr4_rdr, ecliptic_y[0] , 'ECLIPTIC_Y' fxbaddcol, 17, hdr4_rdr, ecliptic_z[0] , 'ECLIPTIC_Z' fxbaddcol, 18, hdr4_rdr, stim_drift1[0] , 'STIM_DRIFT1' fxbaddcol, 19, hdr4_rdr, stim_drift2[0] , 'STIM_DRIFT2' fxbaddcol, 20, hdr4_rdr, other_inst_interf[0] , 'OTHER_INST_INTERF' fxbaddcol, 21, hdr4_rdr, global_data_qual[0] , 'GLOBAL_DATA_QUAL' fxbaddcol, 22, hdr4_rdr, heuristic_quality[0] , 'HEURISTIC_QUALITY' fxbaddcol, 23, hdr4_rdr, dark_factor[0] , 'DARK_FACTOR' fxbaddcol, 24, hdr4_rdr, stray_factor[0] , 'STRAY_FACTOR' fxbaddcol, 25, hdr4_rdr, off_nadir_pointing[0], 'OFF_NADIR_POINTING' fxbcreate, lun, outfiles[ifile], hdr4_rdr IF (nslow GT 0L) THEN BEGIN column_names = ['SCUT_TIME_SLOW','ORBIT_NUMBER','SOLAR_ZENITH_ANGLE','SC_SHADOW_FLAG','SUBSOLAR_LATITUDE', $ 'SUBSOLAR_LONGITUDE','SUBSC_LATITUDE','SUBSC_LONGITUDE','SC_ZENITH_ANGLE','EARTH_PHASE_ANGLE', $ 'PHASE_ANGLE','AZIMUTH_ANGLE','RA','DEC','ECLIPTIC_X','ECLIPTIC_Y','ECLIPTIC_Z','STIM_DRIFT1', $ 'STIM_DRIFT2','OTHER_INST_INTERF','GLOBAL_DATA_QUAL','HEURISTIC_QUALITY','DARK_FACTOR', $ 'STRAY_FACTOR','OFF_NADIR_POINTING'] fxbwritm, lun, column_names, scut_time_slow, orbit_number, solar_zenith_angle, sc_shadow_flag, $ subsolar_latitude, subsolar_longitude, subsc_latitude, subsc_longitude, $ sc_zenith_angle, earth_phase_angle, phase_angle, azimuth_angle, ra, dec, $ ecliptic_x, ecliptic_y, ecliptic_z, stim_drift1, stim_drift2, other_inst_interf, $ global_data_qual, heuristic_quality, dark_factor, stray_factor, off_nadir_pointing, BUFFERSIZE=1048576L ENDIF fxbfinish, lun mwrfits, histogram_cube, outfiles[ifile], hdr5_rdr ; calibrated histogram mode data, extension 5, image cube ; calibrated calculated countrate, extension 6, binary table IF (ncrentries GT 0L) THEN BEGIN fxbhmake, hdr6_rdr, ncrentries, 'Reduced Count Rate', EXTVER = 1 fxbaddcol, 1, hdr6_rdr, cr_hack_time[0] , 'HACK_TIME' fxbaddcol, 2, hdr6_rdr, cr_scut_time[0] , 'SCUT_TIME' fxbaddcol, 3, hdr6_rdr, LONG(cr_countrate[0]), 'COUNT_RATE' fxbcreate, lun, outfiles[ifile], hdr6_rdr column_names = ['HACK_TIME','SCUT_TIME','COUNT_RATE'] fxbwritm, lun, column_names, cr_hack_time, cr_scut_time, LONG(cr_countrate), BUFFERSIZE=1048576L fxbfinish, lun ENDIF ELSE $ ; if we have no calibrated counts, just write out old EDR extension mwrfits, tab4, outfiles[ifile], hdr4 ; calculated countrate, extension 6, binary table mwrfits, tab5, outfiles[ifile], hdr5 ; high-rate LTS data, extension 7, binary table mwrfits, tab6, outfiles[ifile], hdr6 ; housekeeping table, extension 8, binary table mwrfits, waveimage, outfiles[ifile], hdr9_rdr ; wavelength lookup image, extension 9, image END_OF_LOOP: ; Check for non-zero error status before ending loop. ; If so, then an error has occurred and the calibrated file could not be created. IF (error_status NE 0) THEN BEGIN HELP, CALLS = calls lamp_mike_log, '%% Fatal error at ' + calls, ERRNUM = 2 IF (STRLEN(!ERROR_STATE.MSG) GT 0) THEN $ lamp_mike_log, !ERROR_STATE.MSG, ERRNUM = 2 lamp_mike_log, 'File ' + infiles[ifile] + ' not processed!', ERRNUM = 2 ENDIF ; Write log messages for this file to output status file, if requested IF (KEYWORD_SET(mikeflags.statflag)) THEN BEGIN istat_last = N_ELEMENTS(logmessages) - 1 OPENW, lun, statfiles[ifile], /GET_LUN IF (istat_last GE istat_first) THEN PRINTF, lun, logmessages[istat_first : istat_last] FREE_LUN, lun ENDIF ENDFOR ; loop over input files ; Program clean-up END_OF_PROGRAM: ; Unload spice kernels cspice_unload, metakernel ; Compute wall clock time elapsed during processing stop_time = SYSTIME(1); seconds since January 1, 1970 elapsed_time = stop_time - start_time; seconds elapsed_time = STRTRIM(STRING(elapsed_time/60.0, FORMAT = '(F20.2)'), 2); minutes message = 'Total processing time = ' + elapsed_time + ' minutes.' lamp_mike_log, message, VERBOSELEVEL = 1 ; Write internal log to file IF (KEYWORD_SET(logfile)) THEN BEGIN OPENW, lun, logfile, /GET_LUN PRINTF, lun, logmessages FREE_LUN, lun ENDIF END ;+ ; NAME: ; lamp_mike_sysv ; ; PURPOSE: ; Define global system variables used by lamp_mike ; ; CATEGORY: ; Data definition ; ; CALLING SEQUENCE: ; lamp_mike_sysv ; ; OPTIONAL INPUTS: ; None ; ; KEYWORD PARAMETERS: ; None ; ; OUTPUTS: ; None ; ; COMMON BLOCKS: ; None ; ; SIDE EFFECTS: ; None ; ; RESTRICTIONS: ; None ; ; PROCEDURES CALLED: ; cspice_dpr() ; ; EXAMPLE: ; lamp_mike_sysv ; ; MODIFICATION HISTORY: ; V0.10 19 November 2008 DEK - initial implementation ; V0.11 28 September 2009 DEK - changed default directories and files to point to "nominal" phase values instead of "commissioning" ; - changed nadir-pointing threshold to 6.0 degrees from 3.0 degrees ; V0.12 02 October 2009 DEK - changed LAMP deadtime from 18 usec to 17.3 usec based on commissioning data ; V0.13 10 November 2009 DEK - updated version to 0.13 for new wavelength and flatfield calibrations ; V0.14 16 February 2010 DEK - updated version to 0.14 for new (v0.0.2) LAMP instrument kernel ; V0.15 02 March 2010 DEK - updated definitions to include Kurt's "B" dark regions and also his stray factor regions ; V0.16 27 April 2010 DEK - updated extension 6 contents from "Calibrated Calculated Count Rate" to "Reduced Count Rate" ; based on Kurt's email of 08 March 2010. ; ;- PRO lamp_mike_sysv ; Make all Mike system variables read-only by default read_only = 1 ; Code identification DEFSYSV, '!MIKE_SL2PROC', 'LAMP-MIKE' , read_only DEFSYSV, '!MIKE_SL2VER' , '0.16' , read_only DEFSYSV, '!MIKE_SL2DATE', '2010-04-27', read_only ; Default directories and files DEFSYSV, '!MIKE_DEFAULT_CALIBRATION_DIR', '/home/soc/pipeline/nominal/data/calibration' , read_only DEFSYSV, '!MIKE_DEFAULT_TEMP_DIR' , '/var/tmp' , read_only DEFSYSV, '!MIKE_DEFAULT_SPICE_DIR' , '/home/soc/pipeline/nominal/data/spice/lamp/mk' , read_only DEFSYSV, '!MIKE_LRO_ORBIT_DATABASE_FILE', '/home/soc/pipeline/nominal/data/support/LRO_orbit_database.txt', read_only ; Instrument quantities DEFSYSV, '!LAMP_PIXELLIST_MODE', 0 , read_only DEFSYSV, '!LAMP_HISTOGRAM_MODE', 1 , read_only DEFSYSV, '!LAMP_DATA_PATTERN' , 0 , read_only DEFSYSV, '!LAMP_HVPS_A_ON' , 2 , read_only DEFSYSV, '!LAMP_HVPS_B_ON' , 3 , read_only DEFSYSV, '!LAMP_HVPS_AB_ON' , 4 , read_only DEFSYSV, '!LAMP_FOV_X_MIN' , 125 , read_only DEFSYSV, '!LAMP_FOV_X_MAX' , 900 , read_only DEFSYSV, '!LAMP_FOV_X_MIN_CR' , 140 , read_only DEFSYSV, '!LAMP_FOV_X_MAX_CR' , 875 , read_only DEFSYSV, '!LAMP_FOV_Y_MIN' , 4 , read_only DEFSYSV, '!LAMP_FOV_Y_MAX' , 24 , read_only DEFSYSV, '!LAMP_FOV_ID' , 0 , read_only DEFSYSV, '!LAMP_DARK1_X_MIN' , 150 , read_only DEFSYSV, '!LAMP_DARK1_X_MAX' , 450 , read_only DEFSYSV, '!LAMP_DARK1_Y_MIN' , 25 , read_only DEFSYSV, '!LAMP_DARK1_Y_MAX' , 30 , read_only DEFSYSV, '!LAMP_DARK1_ID' , 1 , read_only DEFSYSV, '!LAMP_DARK2_X_MIN' , 550 , read_only DEFSYSV, '!LAMP_DARK2_X_MAX' , 850 , read_only DEFSYSV, '!LAMP_DARK2_Y_MIN' , 25 , read_only DEFSYSV, '!LAMP_DARK2_Y_MAX' , 30 , read_only DEFSYSV, '!LAMP_DARK2_ID' , 2 , read_only DEFSYSV, '!LAMP_DARK3_X_MIN' , 150 , read_only DEFSYSV, '!LAMP_DARK3_X_MAX' , 450 , read_only DEFSYSV, '!LAMP_DARK3_Y_MIN' , 1 , read_only DEFSYSV, '!LAMP_DARK3_Y_MAX' , 3 , read_only DEFSYSV, '!LAMP_DARK3_ID' , 3 , read_only DEFSYSV, '!LAMP_DARK4_X_MIN' , 550 , read_only DEFSYSV, '!LAMP_DARK4_X_MAX' , 850 , read_only DEFSYSV, '!LAMP_DARK4_Y_MIN' , 1 , read_only DEFSYSV, '!LAMP_DARK4_Y_MAX' , 3 , read_only DEFSYSV, '!LAMP_DARK4_ID' , 4 , read_only DEFSYSV, '!LAMP_STIM1_X_MIN' , 0 , read_only DEFSYSV, '!LAMP_STIM1_X_MAX' , 31 , read_only DEFSYSV, '!LAMP_STIM1_Y' , 21 , read_only DEFSYSV, '!LAMP_STIM1_ID' , 5 , read_only DEFSYSV, '!LAMP_STIM2_X_MIN' , 992 , read_only DEFSYSV, '!LAMP_STIM2_X_MAX' , 1023 , read_only DEFSYSV, '!LAMP_STIM2_Y' , 9 , read_only DEFSYSV, '!LAMP_STIM2_ID' , 6 , read_only DEFSYSV, '!LAMP_STRAY1_X_MIN' , 450 , read_only DEFSYSV, '!LAMP_STRAY1_X_MAX' , 550 , read_only DEFSYSV, '!LAMP_STRAY1_Y_MIN' , 3 , read_only DEFSYSV, '!LAMP_STRAY1_Y_MAX' , 3 , read_only DEFSYSV, '!LAMP_STRAY1_ID' , 7 , read_only DEFSYSV, '!LAMP_STRAY2_X_MIN' , 450 , read_only DEFSYSV, '!LAMP_STRAY2_X_MAX' , 550 , read_only DEFSYSV, '!LAMP_STRAY2_Y_MIN' , 25 , read_only DEFSYSV, '!LAMP_STRAY2_Y_MAX' , 25 , read_only DEFSYSV, '!LAMP_STRAY2_ID' , 8 , read_only DEFSYSV, '!LAMP_STIM_ON' , 2 , read_only DEFSYSV, '!LAMP_WAVELENGTH_LO' , 575.0 , read_only ; Angstroms DEFSYSV, '!LAMP_WAVELENGTH_HI' , 1965.0 , read_only ; Angstroms DEFSYSV, '!LAMP_APDOOR_OPEN' , 3 , read_only DEFSYSV, '!LAMP_APDOOR_CLOSED' , 2 , read_only DEFSYSV, '!LAMP_FRAME_DTMAX' , 1.0 , read_only ; seconds DEFSYSV, '!LAMP_DEADTIME' , 1.73D-5 , read_only ; seconds DEFSYSV, '!LAMP_PHD_X_MIN' , 0 , read_only DEFSYSV, '!LAMP_PHD_X_MAX' , 31 , read_only DEFSYSV, '!LAMP_PHD_Y_MIN' , 1 , read_only DEFSYSV, '!LAMP_PHD_Y_MAX' , 2 , read_only ; SPICE quantities DEFSYSV, '!SPICE_LRO_ID' , -85 , read_only DEFSYSV, '!SPICE_LRO_SC_BUS_ID' , -85000 , read_only DEFSYSV, '!SPICE_LRO_LAMP_ID' , -85300 , read_only DEFSYSV, '!SPICE_SUN_ID' , 10 , read_only DEFSYSV, '!SPICE_MOON_ID' , 301 , read_only DEFSYSV, '!SPICE_EARTH_ID' , 399 , read_only DEFSYSV, '!SPICE_MOON_FRAME' , 'MOON_ME' , read_only DEFSYSV, '!SPICE_EARTH_FRAME' , 'IAU_EARTH' , read_only DEFSYSV, '!SPICE_METHOD' , ['ELLIPSOID', 'NEARPOINT:ELLIPSOID'] , read_only DEFSYSV, '!SPICE_OBSRVR' , STRING(!SPICE_LRO_ID) , read_only DEFSYSV, '!SPICE_TARGET' , STRING(!SPICE_MOON_ID) , read_only DEFSYSV, '!SPICE_ABCORR' , 'LT+S' , read_only DEFSYSV, '!SPICE_LAMP_FRAME' , 'LRO_LAMP' , read_only DEFSYSV, '!SPICE_SC_BUS_FRAME' , 'LRO_SC_BUS' , read_only DEFSYSV, '!SPICE_ECLIPTIC_FRAME' , 'ECLIPJ2000' , read_only DEFSYSV, '!SPICE_EQUATORIAL_FRAME', 'J2000' , read_only DEFSYSV, '!SPICE_BADANGLE' , -999.0D0 , read_only DEFSYSV, '!SPICE_DPR' , cspice_dpr() , read_only ; Miscellaneous DEFSYSV, '!LAMP_NADIR_PTG_THRESH', 6.0D0, read_only ; degrees END ;+ ; NAME: ; lamp_mike_sysv_v3 ; ; PURPOSE: ; Define global system variables used by lamp_mike ; ; CATEGORY: ; Data definition ; ; CALLING SEQUENCE: ; lamp_mike_sysv_v3 ; ; OPTIONAL INPUTS: ; None ; ; KEYWORD PARAMETERS: ; None ; ; OUTPUTS: ; None ; ; COMMON BLOCKS: ; None ; ; SIDE EFFECTS: ; None ; ; RESTRICTIONS: ; None ; ; PROCEDURES CALLED: ; cspice_dpr() ; ; EXAMPLE: ; lamp_mike_sysv_v3 ; ; MODIFICATION HISTORY: ; V0.10 19 November 2008 DEK - initial implementation ; V0.11 28 September 2009 DEK - changed default directories and files to point to "nominal" phase values instead of "commissioning" ; - changed nadir-pointing threshold to 6.0 degrees from 3.0 degrees ; V0.12 02 October 2009 DEK - changed LAMP deadtime from 18 usec to 17.3 usec based on commissioning data ; V0.13 10 November 2009 DEK - updated version to 0.13 for new wavelength and flatfield calibrations ; V0.14 16 February 2010 DEK - updated version to 0.14 for new (v0.0.2) LAMP instrument kernel ; V0.15 02 March 2010 DEK - updated definitions to include Kurt's "B" dark regions and also his stray factor regions ; V0.16 27 April 2010 DEK - updated extension 6 contents from "Calibrated Calculated Count Rate" to "Reduced Count Rate" ; based on Kurt's email of 08 March 2010. ; V0.17 24 January 2010 DEK - implemented use of existing RDRs for geometrical information, plus implemented Randy's ; Lya gain correction ; ;- PRO lamp_mike_sysv_v3 ; Make all Mike system variables read-only by default read_only = 1 ; Code identification DEFSYSV, '!MIKE_SL2PROC', 'LAMP-MIKE' , read_only DEFSYSV, '!MIKE_SL2VER' , '0.17' , read_only DEFSYSV, '!MIKE_SL2DATE', '2011-01-24', read_only ; Default directories and files DEFSYSV, '!MIKE_DEFAULT_CALIBRATION_DIR', '/home/soc/pipeline/nominal/data/calibration' , read_only DEFSYSV, '!MIKE_DEFAULT_TEMP_DIR' , '/var/tmp' , read_only DEFSYSV, '!MIKE_DEFAULT_SPICE_DIR' , '/home/soc/pipeline/nominal/data/spice/lamp/mk' , read_only DEFSYSV, '!MIKE_LRO_ORBIT_DATABASE_FILE', '/home/soc/pipeline/nominal/data/support/LRO_orbit_database.txt', read_only ; Instrument quantities DEFSYSV, '!LAMP_PIXELLIST_MODE', 0 , read_only DEFSYSV, '!LAMP_HISTOGRAM_MODE', 1 , read_only DEFSYSV, '!LAMP_DATA_PATTERN' , 0 , read_only DEFSYSV, '!LAMP_HVPS_A_ON' , 2 , read_only DEFSYSV, '!LAMP_HVPS_B_ON' , 3 , read_only DEFSYSV, '!LAMP_HVPS_AB_ON' , 4 , read_only DEFSYSV, '!LAMP_FOV_X_MIN' , 125 , read_only DEFSYSV, '!LAMP_FOV_X_MAX' , 900 , read_only DEFSYSV, '!LAMP_FOV_X_MIN_CR' , 140 , read_only DEFSYSV, '!LAMP_FOV_X_MAX_CR' , 875 , read_only DEFSYSV, '!LAMP_FOV_Y_MIN' , 4 , read_only DEFSYSV, '!LAMP_FOV_Y_MAX' , 24 , read_only DEFSYSV, '!LAMP_FOV_ID' , 0 , read_only DEFSYSV, '!LAMP_DARK1_X_MIN' , 150 , read_only DEFSYSV, '!LAMP_DARK1_X_MAX' , 450 , read_only DEFSYSV, '!LAMP_DARK1_Y_MIN' , 25 , read_only DEFSYSV, '!LAMP_DARK1_Y_MAX' , 30 , read_only DEFSYSV, '!LAMP_DARK1_ID' , 1 , read_only DEFSYSV, '!LAMP_DARK2_X_MIN' , 550 , read_only DEFSYSV, '!LAMP_DARK2_X_MAX' , 850 , read_only DEFSYSV, '!LAMP_DARK2_Y_MIN' , 25 , read_only DEFSYSV, '!LAMP_DARK2_Y_MAX' , 30 , read_only DEFSYSV, '!LAMP_DARK2_ID' , 2 , read_only DEFSYSV, '!LAMP_DARK3_X_MIN' , 150 , read_only DEFSYSV, '!LAMP_DARK3_X_MAX' , 450 , read_only DEFSYSV, '!LAMP_DARK3_Y_MIN' , 1 , read_only DEFSYSV, '!LAMP_DARK3_Y_MAX' , 3 , read_only DEFSYSV, '!LAMP_DARK3_ID' , 3 , read_only DEFSYSV, '!LAMP_DARK4_X_MIN' , 550 , read_only DEFSYSV, '!LAMP_DARK4_X_MAX' , 850 , read_only DEFSYSV, '!LAMP_DARK4_Y_MIN' , 1 , read_only DEFSYSV, '!LAMP_DARK4_Y_MAX' , 3 , read_only DEFSYSV, '!LAMP_DARK4_ID' , 4 , read_only DEFSYSV, '!LAMP_STIM1_X_MIN' , 0 , read_only DEFSYSV, '!LAMP_STIM1_X_MAX' , 31 , read_only DEFSYSV, '!LAMP_STIM1_Y' , 21 , read_only DEFSYSV, '!LAMP_STIM1_ID' , 5 , read_only DEFSYSV, '!LAMP_STIM2_X_MIN' , 992 , read_only DEFSYSV, '!LAMP_STIM2_X_MAX' , 1023 , read_only DEFSYSV, '!LAMP_STIM2_Y' , 9 , read_only DEFSYSV, '!LAMP_STIM2_ID' , 6 , read_only DEFSYSV, '!LAMP_STRAY1_X_MIN' , 450 , read_only DEFSYSV, '!LAMP_STRAY1_X_MAX' , 550 , read_only DEFSYSV, '!LAMP_STRAY1_Y_MIN' , 3 , read_only DEFSYSV, '!LAMP_STRAY1_Y_MAX' , 3 , read_only DEFSYSV, '!LAMP_STRAY1_ID' , 7 , read_only DEFSYSV, '!LAMP_STRAY2_X_MIN' , 450 , read_only DEFSYSV, '!LAMP_STRAY2_X_MAX' , 550 , read_only DEFSYSV, '!LAMP_STRAY2_Y_MIN' , 25 , read_only DEFSYSV, '!LAMP_STRAY2_Y_MAX' , 25 , read_only DEFSYSV, '!LAMP_STRAY2_ID' , 8 , read_only DEFSYSV, '!LAMP_STIM_ON' , 2 , read_only DEFSYSV, '!LAMP_WAVELENGTH_LO' , 575.0 , read_only ; Angstroms DEFSYSV, '!LAMP_WAVELENGTH_HI' , 1965.0 , read_only ; Angstroms DEFSYSV, '!LAMP_APDOOR_OPEN' , 3 , read_only DEFSYSV, '!LAMP_APDOOR_CLOSED' , 2 , read_only DEFSYSV, '!LAMP_FRAME_DTMAX' , 1.0 , read_only ; seconds DEFSYSV, '!LAMP_DEADTIME' , 1.73D-5 , read_only ; seconds DEFSYSV, '!LAMP_PHD_X_MIN' , 0 , read_only DEFSYSV, '!LAMP_PHD_X_MAX' , 31 , read_only DEFSYSV, '!LAMP_PHD_Y_MIN' , 1 , read_only DEFSYSV, '!LAMP_PHD_Y_MAX' , 2 , read_only DEFSYSV, '!LAMP_LYA_X_MIN' , 470 , read_only DEFSYSV, '!LAMP_LYA_X_MAX' , 515 , read_only ; SPICE quantities DEFSYSV, '!SPICE_LRO_ID' , -85 , read_only DEFSYSV, '!SPICE_LRO_SC_BUS_ID' , -85000 , read_only DEFSYSV, '!SPICE_LRO_LAMP_ID' , -85300 , read_only DEFSYSV, '!SPICE_SUN_ID' , 10 , read_only DEFSYSV, '!SPICE_MOON_ID' , 301 , read_only DEFSYSV, '!SPICE_EARTH_ID' , 399 , read_only DEFSYSV, '!SPICE_MOON_FRAME' , 'MOON_ME' , read_only DEFSYSV, '!SPICE_EARTH_FRAME' , 'IAU_EARTH' , read_only DEFSYSV, '!SPICE_METHOD' , ['ELLIPSOID', 'NEARPOINT:ELLIPSOID'] , read_only DEFSYSV, '!SPICE_OBSRVR' , STRING(!SPICE_LRO_ID) , read_only DEFSYSV, '!SPICE_TARGET' , STRING(!SPICE_MOON_ID) , read_only DEFSYSV, '!SPICE_ABCORR' , 'LT+S' , read_only DEFSYSV, '!SPICE_LAMP_FRAME' , 'LRO_LAMP' , read_only DEFSYSV, '!SPICE_SC_BUS_FRAME' , 'LRO_SC_BUS' , read_only DEFSYSV, '!SPICE_ECLIPTIC_FRAME' , 'ECLIPJ2000' , read_only DEFSYSV, '!SPICE_EQUATORIAL_FRAME', 'J2000' , read_only DEFSYSV, '!SPICE_BADANGLE' , -999.0D0 , read_only DEFSYSV, '!SPICE_DPR' , cspice_dpr() , read_only ; Miscellaneous DEFSYSV, '!LAMP_NADIR_PTG_THRESH', 6.0D0, read_only ; degrees END ;+ ; NAME: ; lamp_mike_temp_correction ; ; PURPOSE: ; Return a 1024 length vector of wavelength bin numbers corrected ; for stim position or, if the REVERSE keyword is set, a 1024 ; length vector of adjusted wavelength bin numbers such that they ; would correct back to integral bin numbers [0..1023] based on the ; input stim positions if the forward correction were applied ; ; CATEGORY: ; Data processing ; ; CALLING SEQUENCE: ; x_corrected = lamp_mike_temp_correction(stimpos, [REVERSE=reverse]) ; ; INPUTS: ; stimpos = [stim1, stim2], where ; stim1 -- position (centroid) of the left stim pulse ; stim2 -- position (centroid) of the right stim pulse ; reverse = keyword that controls forward or reverse correction ; ; OUTPUTS: ; Returns a 1024 length vector of wavelength bin numbers corrected ; for stim position or, if the REVERSE keyword is set, a 1024 ; length vector of adjusted wavelength bin numbers such that they ; would correct back to integral bin numbers [0..1023] based on ; the input stim positions if the forward correction were applied; ; to obtain the actual wavelengths associated with these bin ; numbers, for the forward correction, one needs to call: ; ; lambda = INTERPOL(lambda, x, x_corrected) ; ; to obtain an adjusted image registered to the fiducial ; wavelength scale, one needs to do the following after getting ; the x_corrected array with the /REVERSE keyword specified: ; ; correction = x_corrected # (FLTARR(32) + 1.0) ; yindex = (FLTARR(1024) + 1.0) # FINDGEN(32) ; image = INTERPOLATE(image, correction, yindex, CUBIC = ; -0.5) > 0.0 ; ; or, equivalently, ; ; yindex = FINDGEN(32) ; image = INTERPOLATE(image, x_corrected, yindex, CUBIC = ; -0.5, /GRID) > 0.0 ; ; COMMON BLOCKS: ; None. ; ; SIDE EFFECTS: ; None. ; ; NOTES: ; Algorithm taken from Wilkenson et al. 2001. Code is a direct ; copy of Andrew Steffl's code for P-Alice. Waiting on a ; reply from D. Slater in order to update for LAMP. ; ; MODIFICATION HISTORY: ; V1.0 18 Nov 2008 D.E. Kaufmann ; V2.0 06 Nov 2009 D.E. Kaufmann, updated with Andrew's fiducial ; values, which are taken from the first minute's worth of data ; from the file LAMP_ENG_0273802038_01.fit taken on 2009-248 while ; observing the lunar dayside through the pinhole aperture. ;- FUNCTION lamp_mike_temp_correction, stimpos, REVERSE=reverse S_x1 = 10.102996 S_x2 = 1008.1874 Sp_x1 = stimpos[0] Sp_x2 = stimpos[1] IF (KEYWORD_SET(reverse)) THEN BEGIN RETURN, Sp_x1 + (FINDGEN(1024) - S_x1) * (Sp_x2 - Sp_x1) / (S_x2 - S_x1) ENDIF ELSE BEGIN RETURN, S_x1 + (FINDGEN(1024) - Sp_x1) * (S_x2 - S_x1) / (Sp_x2 - Sp_x1) ENDELSE END ;+ ; NAME: ; lamp_mike_v3 ; ; PURPOSE: ; Process CODMAC Level 2 (Engineering Level) LAMP data ; into CODMAC Level 3 (Science Level) LAMP data ; ; CATEGORY: ; Data processing ; ; CALLING SEQUENCE: ; lamp_mike, INFILES=infiles, GEOMFILES=geomfiles, OUTFILES=outfiles, ; STATFILES=statfiles, CALIBRATION_DIR=calibration_dir, ; SPICE_DIR=spice_dir, TEMP_DIR=temp_dir, AEFFFLAG=aeffflag, ; AEFFFILE=aefffile, BADFLAG=badflag, BADFILE=badfile, ; DARKFLAG=darkflag, DARKFILE=darkfile, DEADFLAG=deadflag, ; FLATFLAG=flatflag, FLATFILE=flatfile, PHDFLAG=phdflag, ; PHDFILE=phdfile, RECTFLAG=rectflag, METAKERNEL=metakernel, ; LOGFILE=logfile, VERBOSE=verbose ; ; OPTIONAL INPUTS: ; None. ; ; KEYWORD PARAMETERS: ; INFILES = value: String array of input CODMAC Level 2 file names to ; process. If not supplied, procedure will attempt to ; get this input from environment variable MIKE_IN_FILE. ; If not supplied in either way, procedure errors out. ; GEOMFILES = value: String array of input CODMAC Level 3 file names to ; use as the respective sources of all geometrical data ; rather than using the SPICE toolkit. This greatly ; speeds up reprocessing when only radiometric updates ; are needed to existing RDR files. If this argument is ; not present nor defined by means of the environment ; variable MIKE_GEOM_FILE, SPICE will be used by default. ; OUTFILES = value: String array of output CODMAC Level 3 file names to ; write. If not supplied, procedure will attempt to get ; this input from environment variable MIKE_OUT_FILE. ; Number of outfiles must match number of infiles. ; STATFILES = value: String array of output processing status file names ; to write. If not supplied, procedure will attempt to ; get this input from environment variable MIKE_OUT_STATUS. ; If specified, number of output status files must match ; number of infiles. ; CALIBRATION_DIR = value: String specifying location of root directory ; for calibration data. If not specified, procedure will ; attempt to use default value. ; SPICE_DIR = value: String specifying location of root directory for SPICE ; kernel data. If not specified, procedure will attempt ; to use default value. ; TEMP_DIR = value: String specifying directory location to use for ; generating temporary files, if needed. If not specified, ; procecure will attempt to use default value. ; /AEFFFLAG: Flag to apply effective area calibration (default = 1). ; AEFFFILE = value: String specifying effective area calibration file. ; This must be supplied if AEFFFLAG is set. ; /BADFLAG: Flag to apply bad pixel filtering (default = 0). [NOT IMPLEMENTED] ; BADFILE = value: String specifying bad pixel mask. This must be supplied ; if BADFLAG is set. [NOT IMPLEMENTED] ; /DARKFLAG: Flag to apply dark count correction to data taken in histogram ; mode (default = 1). For data taken in pixel list mode, ; no dark count correction is applied at this level. ; Instead, the dark count level is characterized for later ; correction at the stage of map construction. ; DARKFILE = value: String specifying dark image data file to use for this ; correction. This must be supplied if DARKFLAG is set. ; /DEADFLAG: Flag to apply deadtime correction (default = 1). ; /FLATFLAG: Flag to apply flatfield correction (default = 1). ; FLATFILE = value: String specifying flatfield image data file to use for ; this correction. This must be specified if FLATFLAG is set. ; /PHDFLAG: Flag to apply PHD (Lya gain) correction (default = 1). ; PHDFILE = value: String specifying file containing PHD grand sum data in ; columns (SCUT time [sec], Lya accumulated counts, cumulative ; exposure time [sec], detector housing temperature [deg C]) ; to use for this correction. This must be supplied if ; PHDFLAG is set. ; /RECTFLAG: Flag to apply geometrical distortion correction (default = 1). ; [NOT IMPLEMENTED] ; METAKERNEL = value: String specifying SPICE metakernel to use for performing ; geospatial location of the LAMP pixel list data. This ; must be specified if GEOMFILES is not used, else the procedure ; exits with an error. ; LOGFILE = value: String specifying the name of a file to which processing ; log messages are written. ; VERBOSE = value: Integer (0-3) specifying verbose level in messaging. The ; higher the value, the more verbose the messaging. ; ; OUTPUTS: ; None. ; ; COMMON BLOCKS: ; mike_data: Contains processing flags and log message array ; ; SIDE EFFECTS: ; None. ; ; RESTRICTIONS: ; None. ; ; PROCEDURES CALLED: ; ; INTERNAL ; ======== ; lamp_mike_sysv_v3 ; lamp_mike_log ; poisson_errorbar ; lamp_mike_wavecal ; lamp_mike_event_locate ; lamp_mike_deadtime ; adderr ; ; IDLASTRO ; ======== ; rdfits_struct ; readfits ; readcol ; mrdfits ; fxpar ; sxaddpar ; mkhdr ; mwrfits ; fxbhmake ; fxbaddcol ; fxbcreate ; fxbwritm ; fxbfinish ; ; ICY (SPICE) ; =========== ; cspice_furnsh ; cspice_ktotal ; cspice_kdata ; cspice_pxform ; cspice_getfov ; cspice_mxv ; cspice_gdpool ; cspice_scs2e ; cspice_subpt ; cspice_srfxpt ; cspice_reclat ; cspice_illum ; cspice_recrad ; cspice_spkssb ; cspice_spkezp ; cspice_subpnt ; cspice_vsep ; cspice_subslr ; ; EXAMPLE: ; lamp_mike_v3, /VERBOSE, LOGFILE = 'log', DARKFILE = 'darkLAMP_001.fit', ; FLATFILE = 'flatLAMP_001.fit', PHDFILE = 'lyafluence.dat', ; AEFFFILE = 'aeffLAMP_001.fit', METAKERNEL = 'meta_kernellro.tm' ; ; MODIFICATION HISTORY: ; 17 March 2009, Initial implementation. David E. Kaufmann ; 02 March 2010, Added code for stray factor. David E. Kaufmann ; 27 April 2010, Updated extension 6 contents from "Calibrated Calculated ; Count Rate" to "Reduced Count Rate" based on Kurt's email ; of 08 March 2010. David E. Kaufmann ; 15 December 2010, Updated to use geomfiles instead of SPICE, if present; ; Implemented Lya correction factor for photons with ; detector column (X) values between 470 and 515, inclusive, ; or, equivalently, with nominal wavelengths between 118.8 ; and 127.0 nm. David E. Kaufmann ; ;- PRO lamp_mike_v3, $ INFILES = infiles, $ GEOMFILES = geomfiles, $ OUTFILES = outfiles, $ STATFILES = statfiles, $ CALIBRATION_DIR = calibration_dir, $ SPICE_DIR = spice_dir, $ TEMP_DIR = temp_dir, $ AEFFFLAG = aeffflag, $ AEFFFILE = aefffile, $ BADFLAG = badflag, $ BADFILE = badfile, $ DARKFLAG = darkflag, $ DARKFILE = darkfile, $ DEADFLAG = deadflag, $ FLATFLAG = flatflag, $ FLATFILE = flatfile, $ PHDFLAG = phdflag, $ PHDFILE = phdfile, $ RECTFLAG = rectflag, $ METAKERNEL = metakernel, $ LOGFILE = logfile, $ VERBOSE = verbose COMMON mike_data, mikeflags, logmessages ; Initialize Mike system variables (program constants) lamp_mike_sysv_v3 ; Initialize error reporting start_time = SYSTIME(1); seconds since January 1, 1970 error_status = 0 IF (N_ELEMENTS(mikeflags) EQ 0) THEN $ mikeflags = {version: !MIKE_SL2VER + '[' + !MIKE_SL2DATE + ']', $ mode: '', $ verbose: KEYWORD_SET(verbose) ? verbose : 0, $ log: 1, $ error: 0, $ starttime: start_time, $ statflag: 0} ; Initialize message logging lamp_mike_log, message, /STARTFLAG, /INITFLAG, /TIMEFLAG, VERBOSELEVEL = 1 ; Set default processing flags (1 = apply, 0 = do NOT apply) IF (N_ELEMENTS(aeffflag) EQ 0) THEN aeffflag = 1 IF (N_ELEMENTS(badflag) EQ 0) THEN badflag = 0 IF (N_ELEMENTS(darkflag) EQ 0) THEN darkflag = 1 IF (N_ELEMENTS(deadflag) EQ 0) THEN deadflag = 1 IF (N_ELEMENTS(flatflag) EQ 0) THEN flatflag = 1 IF (N_ELEMENTS(phdflag) EQ 0) THEN phdflag = 1 IF (N_ELEMENTS(rectflag) EQ 0) THEN rectflag = 1 ;; [[DEK REMINDER: Add code below to implement the geometric distortion correction ;; controlled by rectflag whenever we define the needed algorithm.]] ; Validate input (command line arguments take precedence over environment variables) ; Infile(s) IF (N_ELEMENTS(infiles) EQ 0) THEN BEGIN env_infile = GETENV('MIKE_IN_FILE') IF (KEYWORD_SET(env_infile)) THEN result = EXECUTE('infiles = ' + env_infile) ENDIF nfiles = N_ELEMENTS(infiles) IF (nfiles EQ 0) THEN BEGIN message = 'No input file(s) specified.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF ; Geomfile(s) use_spice = 0 IF (N_ELEMENTS(geomfiles) EQ 0) THEN BEGIN env_geomfile = GETENV('MIKE_GEOM_FILE') IF (KEYWORD_SET(env_geomfile)) THEN result = EXECUTE('geomfiles = ' + env_geomfile) ENDIF ntmp = N_ELEMENTS(geomfiles) IF (ntmp EQ 0) THEN BEGIN use_spice = 1 ENDIF ELSE IF (ntmp NE nfiles) THEN BEGIN message = 'Number of geometry files does not match number of input files.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF ; Outfile(s) IF (N_ELEMENTS(outfiles) EQ 0) THEN BEGIN env_outfile = GETENV('MIKE_OUT_FILE') IF (KEYWORD_SET(env_outfile)) THEN result = EXECUTE('outfiles = ' + env_outfile) ENDIF ntmp = N_ELEMENTS(outfiles) IF (ntmp NE nfiles) THEN BEGIN message = 'Number of output files does not match number of input files.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF ; Statfile(s) IF (N_ELEMENTS(statfiles) EQ 0) THEN BEGIN env_statfile = GETENV('MIKE_OUT_STATUS') IF (KEYWORD_SET(env_statfile)) THEN result = EXECUTE('statfiles = ' + env_statfile) ENDIF ntmp = N_ELEMENTS(statfiles) IF (ntmp GT 0) THEN BEGIN mikeflags.statflag = 1 IF (ntmp NE nfiles) THEN BEGIN message = 'Number of output status files is nonzero and does not match number of input files.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF ENDIF ; Calibration directory IF (N_ELEMENTS(calibration_dir) EQ 0) THEN BEGIN env_calibration_dir = GETENV('MIKE_CALIBRATION_DIR') IF (KEYWORD_SET(env_calibration_dir)) THEN calibration_dir = env_calibration_dir ENDIF IF (N_ELEMENTS(calibration_dir) EQ 0) THEN calibration_dir = !MIKE_DEFAULT_CALIBRATION_DIR ; SPICE directory IF (use_spice) THEN BEGIN IF (N_ELEMENTS(spice_dir) EQ 0) THEN BEGIN env_spice_dir = GETENV('MIKE_SPICE_DIR') IF (KEYWORD_SET(env_spice_dir)) THEN spice_dir = env_spice_dir ENDIF IF (N_ELEMENTS(spice_dir) EQ 0) THEN spice_dir = !MIKE_DEFAULT_SPICE_DIR ENDIF ; Temporary file directory IF (N_ELEMENTS(temp_dir) EQ 0) THEN BEGIN env_temp_dir = GETENV('MIKE_TEMP_DIR') IF (KEYWORD_SET(env_temp_dir)) THEN temp_dir = env_temp_dir ENDIF IF (N_ELEMENTS(temp_dir) EQ 0) THEN temp_dir = !MIKE_DEFAULT_TEMP_DIR ; Dark correction file IF (KEYWORD_SET(darkflag)) THEN BEGIN IF (KEYWORD_SET(darkfile)) THEN BEGIN darkfile = calibration_dir + PATH_SEP() + 'dark' + PATH_SEP() + darkfile rdfits_struct, darkfile, darkstruct, /SILENT darkimage = darkstruct.im0 darkerror = darkstruct.im1 ENDIF ELSE BEGIN message = 'Dark correction file not specified.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDELSE ENDIF ; Flatfield correction file ; flatAO : HVPS = A , aperture door = Open ; flatBO : HVPS = B , aperture door = Open ; flatABO: HVPS = both, aperture door = Open ; flatAC : HVPS = A , aperture door = Closed ; flatBC : HVPS = B , aperture door = Closed ; flatABC: HVPS = both, aperture door = Closed IF (KEYWORD_SET(flatflag)) THEN BEGIN IF (KEYWORD_SET(flatfile)) THEN BEGIN flatfile = calibration_dir + PATH_SEP() + 'flat' + PATH_SEP() + flatfile flatAO = mrdfits(flatfile, 1, hdr1, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading ' + flatfile + ', status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF flatAC = mrdfits(flatfile, 2, hdr2, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading ' + flatfile + ', status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF flatBO = mrdfits(flatfile, 3, hdr3, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading ' + flatfile + ', status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF flatBC = mrdfits(flatfile, 4, hdr4, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading ' + flatfile + ', status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF flatABO = mrdfits(flatfile, 5, hdr5, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading ' + flatfile + ', status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF flatABC = mrdfits(flatfile, 6, hdr6, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading ' + flatfile + ', status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF ENDIF ELSE BEGIN message = 'Flatfield correction file not specified.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDELSE ENDIF ; Effective area calibration file ; aeffAO : HVPS = A , aperture door = Open ; aeffBO : HVPS = B , aperture door = Open ; aeffABO: HVPS = both, aperture door = Open ; aeffAC : HVPS = A , aperture door = Closed ; aeffBC : HVPS = B , aperture door = Closed ; aeffABC: HVPS = both, aperture door = Closed IF (KEYWORD_SET(aeffflag)) THEN BEGIN IF (KEYWORD_SET(aefffile)) THEN BEGIN aefffile = calibration_dir + PATH_SEP() + 'aeff' + PATH_SEP() + aefffile tab1 = mrdfits(aefffile, 1, hdr1, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading ' + aefffile + ', status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF aeffwave = tab1.WAVELENGTH aeffAO = tab1.EFFECTIVE_AREA_A aeffBO = tab1.EFFECTIVE_AREA_B aeffABO = tab1.EFFECTIVE_AREA_BOTH apratio = fxpar(hdr1, 'APRATIO'); airglow-to-pinhole aperture ratio aeffAC = aeffAO / apratio aeffBC = aeffBO / apratio aeffABC = aeffABO / apratio ENDIF ELSE BEGIN message = 'Effective area calibration file not specified.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDELSE ENDIF ;; [[DEK REMINDER: Be prepared to add code later to handle possible temperature ;; dependence of these calibration steps, once the nature of the ;; dependence is better understood and characterized.]] ;; [[DEK REMINDER: Dark signals outside FOV have been showing extra signal dependent ;; on the source in the FOV for R- and P-Alices. This is difficult ;; to decipher, but most likely some combination of instrumental ;; scattering, electronic noise, etc. Make sure we get a robust set ;; of scattered light measurements during commissioning and be ;; prepared to add code to correct for these effects later on. ; PHD file (for Lya gain correction) IF (KEYWORD_SET(phdflag)) THEN BEGIN READCOL, phdfile, phd_scut_time, phd_ac_lya, FORMAT = 'D,D', SKIPLINE = 1 IF (!ERROR_STATE.CODE LT 0) THEN BEGIN message = 'Error reading columns from PHD correction file.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDIF ENDIF ; SPICE metakernel file IF (use_spice) THEN BEGIN IF (KEYWORD_SET(metakernel)) THEN BEGIN PUSHD, spice_dir cspice_furnsh, metakernel POPD cspice_ktotal, 'ALL', nkernels kernel_list = STRARR(nkernels) FOR k = 0, nkernels - 1 DO BEGIN cspice_kdata, k, 'ALL', file, type, source, handle, found IF (found) THEN kernel_list[k] = FILE_BASENAME(file) ENDFOR ENDIF ELSE BEGIN message = 'SPICE metakernel not specified.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDELSE ENDIF ; Initialize instrument pointing vectors IF (use_spice) THEN BEGIN ; get transformation matrix from LAMP frame to SC_BUS frame cspice_pxform, !SPICE_LAMP_FRAME, !SPICE_SC_BUS_FRAME, 0.0D0, rmatrix; transform assumed time-independent, so set input et = 0 ; get boresight vector cspice_getfov, !SPICE_LRO_LAMP_ID, 66, shape, frame, bsight, bounds ; transform boresight vector to SC_BUS frame cspice_mxv, rmatrix, bsight, bsight ; get pixel center vectors cspice_gdpool, 'INS-85300_PIXEL_CENTERS', 0, 96, values, found IF (found) THEN BEGIN nvalues = N_ELEMENTS(values) pixel_centers = FLTARR(3, nvalues/3, /NOZERO) FOR k = 0, nvalues - 1 DO pixel_centers[k] = values[k] ENDIF ELSE BEGIN message = 'In cspice_gdpool(), unable to find pixel centers.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDELSE ; transform pixel center vectors to SC_BUS frame pixel_centers = TRANSPOSE(MATRIX_MULTIPLY(pixel_centers, rmatrix, /ATRANSPOSE)) ; get pixel corner vectors cspice_gdpool, 'INS-85300_FOV_BOUNDARY_CORNERS', 0, 198, values, found IF (found) THEN BEGIN nvalues = N_ELEMENTS(values) pixel_corners = FLTARR(3, nvalues/3, /NOZERO) FOR k = 0, nvalues - 1 DO pixel_corners[k] = values[k] ENDIF ELSE BEGIN message = 'In cspice_gdpool(), unable to find pixel corners.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_PROGRAM ENDELSE ; transform pixel corner vectors to SC_BUS frame pixel_corners = TRANSPOSE(MATRIX_MULTIPLY(pixel_corners, rmatrix, /ATRANSPOSE)) ENDIF ; use_spice ; Get vector of Poisson errors poisson_errors = DBLARR(2501, /NOZERO) FOR i = 0, 2500 DO poisson_errors[i] = poisson_errorbar(i, /DOUBLE) poisson_errors = FLOAT(poisson_errors) ; Get LRO orbit number and SCUT time vectors from ASCII database file readcol, !MIKE_LRO_ORBIT_DATABASE_FILE, lro_orbit_number, lro_orbit_scut, FORMAT = 'I,D', /SILENT n_lro_orbits = N_ELEMENTS(lro_orbit_number) ; Determine the total number of detector elements in the 4 regions where ; dark signal is being estimated n_dark_elements = (!LAMP_DARK1_X_MAX - !LAMP_DARK1_X_MIN + 1)*(!LAMP_DARK1_Y_MAX - !LAMP_DARK1_Y_MIN + 1) + $ (!LAMP_DARK2_X_MAX - !LAMP_DARK2_X_MIN + 1)*(!LAMP_DARK2_Y_MAX - !LAMP_DARK2_Y_MIN + 1) + $ (!LAMP_DARK3_X_MAX - !LAMP_DARK3_X_MIN + 1)*(!LAMP_DARK3_Y_MAX - !LAMP_DARK3_Y_MIN + 1) + $ (!LAMP_DARK4_X_MAX - !LAMP_DARK4_X_MIN + 1)*(!LAMP_DARK4_Y_MAX - !LAMP_DARK4_Y_MIN + 1) ; Determine the total number of detector elements in the 2 regions where ; the stray light signal is being estimated n_stray_elements = (!LAMP_STRAY1_X_MAX - !LAMP_STRAY1_X_MIN + 1)*(!LAMP_STRAY1_Y_MAX - !LAMP_STRAY1_Y_MIN + 1) + $ (!LAMP_STRAY2_X_MAX - !LAMP_STRAY2_X_MIN + 1)*(!LAMP_STRAY2_Y_MAX - !LAMP_STRAY2_Y_MIN + 1) message = STRING(nfiles) + ' file(s) to process.' lamp_mike_log, message, VERBOSELEVEL = 1 ; Loop over input files and process FOR ifile = 0, nfiles - 1 DO BEGIN istat_first = N_ELEMENTS(logmessages) ; index into logmessages of first entry for this file ; Verify file exists and report processing start found = FILE_TEST(infiles[ifile]) IF (NOT found) THEN BEGIN message = 'File ' + infiles[ifile] + ' not found.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF filenum = STRING(ifile + 1, STRTRIM(nfiles, 2), FORMAT = '(I5,"/",A)') IF (KEYWORD_SET(verbose)) THEN $ PRINT, '% ' + !MIKE_SL2PROC + ': Now processing file ' + infiles[ifile] + ' ' + filenum IF (STRPOS(STRUPCASE(infiles[ifile]), 'SCI') GE 0) THEN BEGIN message = 'File has already been processed!' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Make sure file is valid FITS buffer = BYTARR(10, /NOZERO) OPENR, lun, infiles[ifile], /GET_LUN READU, lun, buffer FREE_LUN, lun IF (STRING(buffer) NE 'SIMPLE = ') THEN BEGIN message = 'File ' + infiles[ifile] + ' is not a valid FITS file.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Check geomfile status if not using SPICE IF (NOT use_spice) THEN BEGIN ; Verify geomfile exists found = FILE_TEST(geomfiles[ifile]) IF (NOT found) THEN BEGIN message = 'File ' + geomfiles[ifile] + ' not found.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Make sure geomfile is valid FITS buffer = BYTARR(10, /NOZERO) OPENR, lun, geomfiles[ifile], /GET_LUN READU, lun, buffer FREE_LUN, lun IF (STRING(buffer) NE 'SIMPLE = ') THEN BEGIN message = 'File ' + geomfiles[ifile] + ' is not a valid FITS file.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Read primary geomfile header geom_hdr0 = headfits(geomfiles[ifile]) ENDIF ; NOT use_spice ; Read in data ; Door open histogram, primary extension, image im0 = mrdfits(infiles[ifile], 0, hdr0, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading extension 0 with mrdfits, status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Door closed histogram, extension 1, image im1 = mrdfits(infiles[ifile], 1, hdr1, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading extension 1 with mrdfits, status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Acquisition list, extension 2, ASCII table tab2 = mrdfits(infiles[ifile], 2, hdr2, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading extension 2 with mrdfits, status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Raw science frame data, extension 3, binary table tab3 = mrdfits(infiles[ifile], 3, hdr3, /UNSIGNED, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading extension 3 with mrdfits, status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Calculated countrate, extension 4, binary table tab4 = mrdfits(infiles[ifile], 4, hdr4, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading extension 4 with mrdfits, status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; High-rate LTS data, extension 5, binary table tab5 = mrdfits(infiles[ifile], 5, hdr5, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading extension 5 with mrdfits, status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Housekeeping table, extension 6, binary table tab6 = mrdfits(infiles[ifile], 6, hdr6, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading extension 6 with mrdfits, status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Check file HK status good_housekeeping = 1 nhkpackets = fxpar(hdr0, 'PACKETS') hkgaps = fxpar(hdr0, 'HKGAPS') IF (nhkpackets EQ 0) THEN BEGIN good_housekeeping = 0 message = 'No housekeeping dataset detected! Processing of current file aborted.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ELSE IF (hkgaps EQ 1) THEN BEGIN message = 'Partial housekeeping dataset detected! Attempting normal processing.' lamp_mike_log, message, ERRNUM = 1, VERBOSELEVEL = verbose ENDIF ; Create wavelength calibration image for output FITS using average detector housing temperature; ; STIM positions, if available, will be used instead of this for calibration of pixellist data avgdetetemp = 0.0D0 IF (good_housekeeping) THEN avgdetetemp = fxpar(hdr0, 'AVRDETHT') ; waveimage = lamp_mike_wavecal(DETETEMP = avgdetetemp) waveimage = lamp_mike_wavecal() ; DEK, use "fiducial" wavelength solution here ;; [[DEK REMINDER: Note that the above image is in some sense "time-averaged" and this could be ;; misleading to end users. We need to be very clear about this when delivering ;; the RDR product files and instructing on their use.]] ; SCUT = hack_clock * hack_time (wrap-extended) + hack_offset hack_clock = fxpar(hdr0, 'HACKCLCK') hack_offset = fxpar(hdr0, 'HACKOFFS') ; Generate some acquisition statistics nframes = fxpar(hdr3, 'NFRAMES') IF (nframes EQ 0) THEN BEGIN message = 'No science dataset detected! Processing of current file aborted.' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF npixellists = 0L nhistograms = 0L nfovevents = 0L ; photon events in detector FOV ncrentries = 0L ; number of countrate table entries hvps_state_sc = -1 ; HVPS state for all pixellist frames in this file; treat as error if it changes first_hack = LONARR(nframes) - 1L skip_frame = BYTARR(nframes) ncrentries_per_frame = LONARR(nframes) FOR iframe = 0L, nframes - 1 DO BEGIN ; iframe indexes into the raw science frame data table (extension 3) ; and the acquisition list table (extension 2) frame_data = (tab3.FRAME_DATA)[*, iframe] frame_number = ISHFT(ISHFT(frame_data[0], 4), -4) ; make sure HK data for frame is present in Acquisition List IF (((tab2.APD_STATE)[iframe] EQ 0) && ((tab2.LTS_STATE)[iframe] EQ 0) && $ ((tab2.STIM_STATE)[iframe] EQ 0) && ((tab2.HVPS_STATE)[iframe] EQ 0)) THEN BEGIN message = 'No housekeeping data for frame ' + STRING(frame_number) + ', skipping...' lamp_mike_log, message, ERRNUM = 1, VERBOSELEVEL = verbose skip_frame[iframe] = 1 CONTINUE ENDIF ; make sure this is a data frame IF ((tab2.FR_PATTERN)[iframe] NE !LAMP_DATA_PATTERN) THEN BEGIN message = 'Frame ' + STRING(frame_number) + ' is not a data frame, skipping...' lamp_mike_log, message, ERRNUM = 1, VERBOSELEVEL = verbose skip_frame[iframe] = 1 CONTINUE ENDIF ; make sure STIMs are on IF ((tab2.STIM_STATE)[iframe] NE !LAMP_STIM_ON) THEN BEGIN message = 'STIM pixels not ON in frame ' + STRING(frame_number) + ', skipping...' lamp_mike_log, message, ERRNUM = 1, VERBOSELEVEL = verbose skip_frame[iframe] = 1 CONTINUE ENDIF ; make sure selected HVPS state is ok if aeff calibration is requested IF (KEYWORD_SET(aeffflag)) THEN BEGIN hvps_state_frame = (tab2.HVPS_STATE)[iframe] IF ((hvps_state_frame NE !LAMP_HVPS_A_ON) && (hvps_state_frame NE !LAMP_HVPS_B_ON) && $ (hvps_state_frame NE !LAMP_HVPS_AB_ON)) THEN BEGIN message = 'Bad HVPS state = ' + STRING(hvps_state_frame) + ' in frame ' + STRING(frame_number) + ', skipping...' lamp_mike_log, message, ERRNUM = 1, VERBOSELEVEL = verbose skip_frame[iframe] = 1 CONTINUE ENDIF IF ((tab2.ACQ_MODE)[iframe] EQ !LAMP_PIXELLIST_MODE) THEN BEGIN IF (hvps_state_sc LT 0) THEN hvps_state_sc = hvps_state_frame $ ELSE IF (hvps_state_frame NE hvps_state_sc) THEN BEGIN message = 'HVPS state not constant for all pixellist frames in this file, aborting...' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose GOTO, END_OF_LOOP ENDIF ENDIF ENDIF IF ((tab2.ACQ_MODE)[iframe] EQ !LAMP_PIXELLIST_MODE) THEN BEGIN ; pixellist mode npixellists += 1 nhacks = 0L ; total number of time hacks in this frame FOR iword = 1L, 32767L DO BEGIN IF (frame_data[iword] NE 0) THEN BEGIN ; frame data word is not zero fill IF (ISHFT(frame_data[iword], -15) EQ 0) THEN BEGIN ; photon event xtmp = ISHFT(ISHFT(frame_data[iword], 6), -6) ytmp = ISHFT(ISHFT(frame_data[iword], 1), -11) location = lamp_mike_event_locate(xtmp, ytmp) IF (location EQ !LAMP_FOV_ID) THEN nfovevents += 1L ENDIF ELSE BEGIN ; time hack IF (first_hack[iframe] LT 0L) THEN first_hack[iframe] = ISHFT(ISHFT(frame_data[iword], 1), -1) nhacks += 1 ENDELSE ENDIF ENDFOR IF ((tab2.HACK_RATE)[iframe] GT 0) THEN BEGIN nhacks_per_cr_interval = MAX([2^(6 - (tab2.HACK_RATE)[iframe]), 1]) ; interval is >= 128 ms ncrentries_per_frame[iframe] = (nhacks - 1) / nhacks_per_cr_interval ncrentries += ncrentries_per_frame[iframe] ENDIF estimated_first_hack = ROUND(((tab2.START_TIME)[iframe] - hack_offset)/hack_clock) IF (first_hack[iframe] GE 0L) THEN BEGIN nrollovers = ROUND((estimated_first_hack - first_hack[iframe])/32768.0D0) first_hack[iframe] += nrollovers*32768L ENDIF ELSE $ ; only when no hacks were detected in frame, thus should never occur first_hack[iframe] = estimated_first_hack ENDIF ELSE $ ; histogram mode nhistograms += 1L ENDFOR IF ((npixellists GT 0L) AND (nfovevents EQ 0L)) THEN BEGIN message = 'File contains pixellist mode data but no events within LAMP FOV' lamp_mike_log, message, ERRNUM = 1, VERBOSELEVEL = verbose ENDIF ; Allocate arrays for fast pixellist table IF (nfovevents GT 0L) THEN BEGIN hack_time = LONARR(nfovevents, /NOZERO) ; wrap-extended values detector_x = INTARR(nfovevents, /NOZERO) detector_y = INTARR(nfovevents, /NOZERO) spatial_row = INTARR(nfovevents, /NOZERO) wavelength = FLTARR(nfovevents, /NOZERO) scut_time = DBLARR(nfovevents, /NOZERO) weighted_count = FLTARR(nfovevents, /NOZERO) weighted_variance = FLTARR(nfovevents, /NOZERO) longitude = FLTARR(nfovevents, /NOZERO) latitude = FLTARR(nfovevents, /NOZERO) longitude_ll = FLTARR(nfovevents, /NOZERO) latitude_ll = FLTARR(nfovevents, /NOZERO) longitude_ul = FLTARR(nfovevents, /NOZERO) latitude_ul = FLTARR(nfovevents, /NOZERO) longitude_ur = FLTARR(nfovevents, /NOZERO) latitude_ur = FLTARR(nfovevents, /NOZERO) longitude_lr = FLTARR(nfovevents, /NOZERO) latitude_lr = FLTARR(nfovevents, /NOZERO) altitude = FLTARR(nfovevents, /NOZERO) incidence_angle = FLTARR(nfovevents, /NOZERO) emission_angle = FLTARR(nfovevents, /NOZERO) deadtime_factor = FLTARR(nfovevents, /NOZERO) pix_data_qual = BYTARR(nfovevents, /NOZERO) apdoor_state = BYTARR(nfovevents, /NOZERO) ; auxiliary array hvps_state = BYTARR(nfovevents, /NOZERO) ; auxiliary array cr_index = LONARR(nfovevents, /NOZERO) ; for cross referencing with the countrate arrays ENDIF ELSE BEGIN ; allocate dummy arrays for FITS file hack_time = LONARR(1, /NOZERO) ; wrap-extended values detector_x = INTARR(1, /NOZERO) detector_y = INTARR(1, /NOZERO) spatial_row = INTARR(1, /NOZERO) wavelength = FLTARR(1, /NOZERO) scut_time = DBLARR(1, /NOZERO) weighted_count = FLTARR(1, /NOZERO) weighted_variance = FLTARR(1, /NOZERO) longitude = FLTARR(1, /NOZERO) latitude = FLTARR(1, /NOZERO) longitude_ll = FLTARR(1, /NOZERO) latitude_ll = FLTARR(1, /NOZERO) longitude_ul = FLTARR(1, /NOZERO) latitude_ul = FLTARR(1, /NOZERO) longitude_ur = FLTARR(1, /NOZERO) latitude_ur = FLTARR(1, /NOZERO) longitude_lr = FLTARR(1, /NOZERO) latitude_lr = FLTARR(1, /NOZERO) altitude = FLTARR(1, /NOZERO) incidence_angle = FLTARR(1, /NOZERO) emission_angle = FLTARR(1, /NOZERO) deadtime_factor = FLTARR(1, /NOZERO) pix_data_qual = BYTARR(1, /NOZERO) ENDELSE IF (ncrentries GT 0L) THEN BEGIN cr_hack_time = LONARR(ncrentries, /NOZERO) cr_scut_time = DBLARR(ncrentries, /NOZERO) ; cr_countrate = FLTARR(ncrentries) cr_countrate = DBLARR(ncrentries) cr_hack_incr = INTARR(ncrentries, /NOZERO) ENDIF ; Allocate data cube for calibrated histogram and error image pairs IF (nhistograms GT 0L) THEN histogram_cube = FLTARR(1024, 32, 2, nhistograms) ; Allocate arrays for slow pixellist table (1 minute intervals) ; scut_scifile_start = fxpar(hdr0, 'SCSTRSEC') scut_scifile_start = MIN([fxpar(hdr0, 'SCSTRSEC'), (tab2.START_TIME)[0]]) ; scut_scifile_end = fxpar(hdr0, 'SCENDSEC') scut_scifile_end = MAX([fxpar(hdr0, 'SCENDSEC'), (tab2.STOP_TIME)[nframes - 1]]) nslow = FLOOR((scut_scifile_end - scut_scifile_start)/60.0D0) + 1 IF (nslow GT 0L) THEN BEGIN scut_time_slow = 60.0D0*(DINDGEN(nslow) + 1.0D0) + scut_scifile_start orbit_number = INTARR(nslow, /NOZERO) solar_zenith_angle = FLTARR(nslow, /NOZERO) sc_shadow_flag = BYTARR(nslow) subsolar_latitude = FLTARR(nslow, /NOZERO) subsolar_longitude = FLTARR(nslow, /NOZERO) subsc_latitude = FLTARR(nslow, /NOZERO) subsc_longitude = FLTARR(nslow, /NOZERO) sc_zenith_angle = FLTARR(nslow, /NOZERO) earth_phase_angle = FLTARR(nslow, /NOZERO) phase_angle = FLTARR(nslow, /NOZERO) azimuth_angle = FLTARR(nslow, /NOZERO) ra = FLTARR(nslow, /NOZERO) dec = FLTARR(nslow, /NOZERO) ecliptic_x = FLTARR(nslow, /NOZERO) ecliptic_y = FLTARR(nslow, /NOZERO) ecliptic_z = FLTARR(nslow, /NOZERO) stim_drift1 = FLTARR(nslow) stim_drift2 = FLTARR(nslow) other_inst_interf = BYTARR(nslow) global_data_qual = INTARR(nslow) heuristic_quality = INTARR(nslow) dark_factor = FLTARR(nslow) stray_factor = FLTARR(nslow) off_nadir_pointing = BYTARR(nslow) ; Auxiliary arrays stim1_data = FLTARR(!LAMP_STIM1_X_MAX - !LAMP_STIM1_X_MIN + 1, nslow) stim2_data = FLTARR(!LAMP_STIM2_X_MAX - !LAMP_STIM2_X_MIN + 1, nslow) slow_table_waveimg = FLTARR(1024, 32, nslow, /NOZERO) IF (KEYWORD_SET(aeffflag)) THEN slow_table_aeffimg = FLTARR(1024, 32, nslow, 3, 2, /NOZERO) stim1x = FINDGEN(!LAMP_STIM1_X_MAX - !LAMP_STIM1_X_MIN + 1) + !LAMP_STIM1_X_MIN stim2x = FINDGEN(!LAMP_STIM2_X_MAX - !LAMP_STIM2_X_MIN + 1) + !LAMP_STIM2_X_MIN ENDIF ELSE BEGIN ; allocate dummy arrays for FITS file scut_time_slow = DINDGEN(1, /NOZERO) orbit_number = INTARR(1, /NOZERO) solar_zenith_angle = FLTARR(1, /NOZERO) sc_shadow_flag = BYTARR(1, /NOZERO) subsolar_latitude = FLTARR(1, /NOZERO) subsolar_longitude = FLTARR(1, /NOZERO) subsc_latitude = FLTARR(1, /NOZERO) subsc_longitude = FLTARR(1, /NOZERO) sc_zenith_angle = FLTARR(1, /NOZERO) earth_phase_angle = FLTARR(1, /NOZERO) phase_angle = FLTARR(1, /NOZERO) azimuth_angle = FLTARR(1, /NOZERO) ra = FLTARR(1, /NOZERO) dec = FLTARR(1, /NOZERO) ecliptic_x = FLTARR(1, /NOZERO) ecliptic_y = FLTARR(1, /NOZERO) ecliptic_z = FLTARR(1, /NOZERO) stim_drift1 = FLTARR(1, /NOZERO) stim_drift2 = FLTARR(1, /NOZERO) other_inst_interf = BYTARR(1, /NOZERO) global_data_qual = INTARR(1, /NOZERO) heuristic_quality = INTARR(1, /NOZERO) dark_factor = FLTARR(1, /NOZERO) stray_factor = FLTARR(1, /NOZERO) off_nadir_pointing = BYTARR(1, /NOZERO) ENDELSE ; Load basic frame data into arrays j = -1L ; FOV photon event index k = -1L ; histogram index m = -1L ; countrate table entry index FOR iframe = 0L, nframes - 1 DO BEGIN IF (skip_frame[iframe]) THEN CONTINUE ; current frame was skipped above frame_data = (tab3.FRAME_DATA)[*, iframe] IF ((tab2.ACQ_MODE)[iframe] EQ !LAMP_PIXELLIST_MODE) THEN BEGIN ; pixellist mode hack_rate = (tab2.HACK_RATE)[iframe] hack_incr = 2^(hack_rate - 1) current_hack = first_hack[iframe] - hack_incr ; back off one hack increment for initial photon events scut = hack_clock * (current_hack + 0.5D0*hack_incr) + hack_offset ; offset by 1/2 hack_incr for midpoint slow_idx = FLOOR((scut - scut_scifile_start)/60.0D0) IF (slow_idx LT 0L) THEN slow_idx = 0L IF (slow_idx GT nslow - 1) THEN slow_idx = nslow - 1 totcounts_this_hack = 0L fovevents_this_hack = 0L drkcounts_this_hack = 0L straycnts_this_hack = 0L IF ((hack_rate GT 0) AND (ncrentries_per_frame[iframe] GT 0)) THEN BEGIN ncrentries_this_frame = 1L nhacks_per_cr_interval = MAX([2^(6 - hack_rate), 1]) m += 1 cr_hack_time[m] = first_hack[iframe] + nhacks_per_cr_interval * hack_incr / 2 cr_scut_time[m] = hack_clock * cr_hack_time[m] + hack_offset cr_hack_incr[m] = nhacks_per_cr_interval * hack_incr cr_hack_counter = -1 ENDIF FOR iword = 1L, 32767L DO BEGIN IF (frame_data[iword] NE 0) THEN BEGIN ; frame data word is not zero fill IF (ISHFT(frame_data[iword], -15) EQ 0) THEN BEGIN ; photon event totcounts_this_hack += 1 xtmp = ISHFT(ISHFT(frame_data[iword], 6), -6) ytmp = ISHFT(ISHFT(frame_data[iword], 1), -11) location = lamp_mike_event_locate(xtmp, ytmp) SWITCH (location) OF !LAMP_FOV_ID : BEGIN fovevents_this_hack += 1 j += 1 hack_time[j] = current_hack detector_x[j] = xtmp detector_y[j] = ytmp scut_time[j] = scut weighted_count[j] = 1.0 weighted_variance[j] = 0.0 deadtime_factor[j] = 1.0 pix_data_qual[j] = 0 hvps_state[j] = (tab2.HVPS_STATE)[iframe] apdoor_state[j] = (tab2.APD_STATE)[iframe] cr_index[j] = m BREAK END !LAMP_STIM1_ID : BEGIN stim1_data[xtmp - !LAMP_STIM1_X_MIN, slow_idx] += 1.0 BREAK END !LAMP_STIM2_ID : BEGIN stim2_data[xtmp - !LAMP_STIM2_X_MIN, slow_idx] += 1.0 BREAK END !LAMP_DARK1_ID : ; falls through !LAMP_DARK2_ID : ; falls through !LAMP_DARK3_ID : ; falls through !LAMP_DARK4_ID : BEGIN drkcounts_this_hack += 1 BREAK END !LAMP_STRAY1_ID : ; falls through !LAMP_STRAY2_ID : BEGIN straycnts_this_hack += 1 BREAK END ENDSWITCH ENDIF ELSE BEGIN ; time hack IF (KEYWORD_SET(deadflag)) THEN BEGIN count_rate = totcounts_this_hack / (hack_clock * hack_incr) ; Hz deadtime_correction = lamp_mike_deadtime(count_rate) FOR i = j - fovevents_this_hack + 1, j DO BEGIN weighted_count[i] *= deadtime_correction deadtime_factor[i] = deadtime_correction ENDFOR dark_factor[slow_idx] += deadtime_correction * drkcounts_this_hack stray_factor[slow_idx] += deadtime_correction * straycnts_this_hack ENDIF ELSE BEGIN dark_factor[slow_idx] += FLOAT(drkcounts_this_hack) stray_factor[slow_idx] += FLOAT(straycnts_this_hack) ENDELSE current_hack += hack_incr expected_rolled_hack = current_hack MOD 32768L IF (expected_rolled_hack LT 0L) THEN expected_rolled_hack += 32768L actual_rolled_hack = ISHFT(ISHFT(frame_data[iword], 1), -1) IF (actual_rolled_hack NE expected_rolled_hack) THEN BEGIN message = ' unexpected hack value (' + STRING(actual_rolled_hack) + ') ' + $ 'in word ' + STRING(iword) + ' of frame ' + STRING(iframe) + '; ' + $ 'expected value (' + STRING(expected_rolled_hack) + ')' lamp_mike_log, message, ERRNUM = 1, VERBOSELEVEL = verbose ENDIF scut = hack_clock * (current_hack + 0.5D0*hack_incr) + hack_offset ; offset by 1/2 hack_incr for midpoint slow_idx = FLOOR((scut - scut_scifile_start)/60.0D0) IF (slow_idx LT 0L) THEN slow_idx = 0L IF (slow_idx GT nslow - 1) THEN slow_idx = nslow - 1 totcounts_this_hack = 0L fovevents_this_hack = 0L drkcounts_this_hack = 0L straycnts_this_hack = 0L IF ((hack_rate GT 0) AND (ncrentries_per_frame[iframe] GT 0)) THEN BEGIN cr_hack_counter += 1 IF ((cr_hack_counter EQ nhacks_per_cr_interval) AND $ (ncrentries_this_frame LT ncrentries_per_frame[iframe])) THEN BEGIN ncrentries_this_frame += 1L m += 1 cr_hack_time[m] = current_hack + nhacks_per_cr_interval * hack_incr / 2 cr_scut_time[m] = hack_clock * cr_hack_time[m] + hack_offset cr_hack_incr[m] = nhacks_per_cr_interval * hack_incr cr_hack_counter = 0 ENDIF ENDIF ENDELSE ENDIF ENDFOR ENDIF ELSE BEGIN ; histogram mode - perform radiometric calibrations and store in histogram data cube image = FLOAT(REFORM(frame_data, 1024, 32)) exptime = (tab2.STOP_TIME)[iframe] - (tab2.START_TIME)[iframe] IF (good_housekeeping) THEN $ indices = WHERE((tab6.SCUT_TIME GE (tab2.START_TIME)[iframe]) AND $ (tab6.SCUT_TIME LE (tab2.STOP_TIME)[iframe])) ; create the wavelength image for this histogram waveimg = waveimage ; use default if STIM and HK data not available IF ((tab2.STIM_STATE)[iframe] EQ !LAMP_STIM_ON) THEN BEGIN ; stim1y = image[!LAMP_STIM1_X_MIN : !LAMP_STIM1_X_MAX, !LAMP_STIM1_Y] stim1y = image[stim1x, !LAMP_STIM1_Y] ; start of Andrew's code for stim1 stim1err = SQRT(stim1y) stim1tot = TOTAL(stim1y) stim1toterr = SQRT(stim1tot) stim1_index1 = stim1y * stim1x stim1err_index1 = stim1err * stim1x stim1indextot = totalerr(stim1_index1, stim1err_index1, err = stim1indextoterr) stimpos1 = diverr(stim1indextot, stim1indextoterr, stim1tot, stim1toterr, err = stimpos1err_i) stimpos1err = stimpos1err_i ; end of Andrew's code stim1pos = (MATRIX_MULTIPLY(stim1x, stim1y, /ATRANSPOSE) / TOTAL(stim1y))[0] ; my value ; stim2y = image[!LAMP_STIM2_X_MIN : !LAMP_STIM2_X_MAX, !LAMP_STIM2_Y] stim2y = image[stim2x, !LAMP_STIM2_Y] ; start of Andrew's code for stim2 stim2err = SQRT(stim2y) stim2tot = TOTAL(stim2y) stim2toterr = SQRT(stim2tot) stim2_index2 = stim2y * stim2x stim2err_index2 = stim2err * stim2x stim2indextot = totalerr(stim2_index2, stim2err_index2, err = stim2indextoterr) stimpos2 = diverr(stim2indextot, stim2indextoterr, stim2tot, stim2toterr, err = stimpos2err_i) stimpos2err = stimpos2err_i ; end of Andrew's code stim2pos = (MATRIX_MULTIPLY(stim2x, stim2y, /ATRANSPOSE) / TOTAL(stim2y))[0] ; my value ; stimpos = [stim1pos, stim2pos] stimpos = [stimpos1, stimpos2] ; Andrew's values stimposerr = [stimpos1err, stimpos2err] ; Andrew's values ; waveimg = lamp_mike_wavecal(STIMPOS = stimpos) ; DEK - do not recompute the wave image to go along with the ; histogram; rather, resample the histogram image to match the ; fiducial wavelength image x_corrected = lamp_mike_temp_correction(stimpos, /REVERSE) yindex = FINDGEN(32) imgtmp = INTERPOLATE(image, x_corrected, yindex, CUBIC = -0.5, /GRID) > 0.0 input_total = TOTAL(image[!LAMP_FOV_X_MIN - 50:!LAMP_FOV_X_MAX + 50, *]) output_total = TOTAL(imgtmp[!LAMP_FOV_X_MIN - 50:!LAMP_FOV_X_MAX + 50, *]) IF (input_total GT 0.0 AND output_total GT 0.0) THEN imgtmp *= input_total / output_total image = imgtmp ENDIF ELSE IF (good_housekeeping AND (N_ELEMENTS(indices) GT 0)) THEN BEGIN detetemp = (tab6.DET_HOUS_TMP)[indices] avgdetetemp = AVG(detetemp) ; waveimg = lamp_mike_wavecal(DETETEMP = avgdetetemp) ENDIF ; zero out PHD and frame header word image[!LAMP_PHD_X_MIN : !LAMP_PHD_X_MAX, !LAMP_PHD_Y_MIN : !LAMP_PHD_Y_MAX] = 0.0 ; PHD image[0, 0] = 0.0 ; frame header word ; get countrate count_rate = TOTAL(image) / exptime ; zero out STIM pixels image[!LAMP_STIM1_X_MIN : !LAMP_STIM1_X_MAX, !LAMP_STIM1_Y] = 0.0 image[!LAMP_STIM2_X_MIN : !LAMP_STIM2_X_MAX, !LAMP_STIM2_Y] = 0.0 ; set saturated values to NaN saturated = WHERE(image EQ 65535, nsaturated) IF (nsaturated GT 0) THEN image[saturated] = !VALUES.F_NAN ; create error image; use gaussian approximation where image > 2500 counts, else poissonian error_image = FLTARR(SIZE(image, /DIMENSIONS), /NOZERO) over2500 = WHERE(image GT 2500.0, nover2500, COMPLEMENT = under2500, NCOMPLEMENT = nunder2500) IF (nover2500 GT 0) THEN error_image[over2500] = SQRT(image[over2500]) IF (nunder2500 GT 0) THEN error_image[under2500] = poisson_errors[image[under2500]] ; divide by exposure time image /= exptime error_image /= exptime ; dead time correction IF (KEYWORD_SET(deadflag)) THEN BEGIN deadtime_correction = lamp_mike_deadtime(count_rate) image *= deadtime_correction error_image *= deadtime_correction message = ' dead time correction = YES' lamp_mike_log, message ENDIF ; dark subtraction IF (KEYWORD_SET(darkflag)) THEN BEGIN image = adderr(image, error_image, darkimage, darkerror, ERR = error_image, /SUBTRACT) message = ' dark subtraction = YES' lamp_mike_log, message ENDIF ; flatfield correction IF (KEYWORD_SET(flatflag)) THEN BEGIN CASE ((tab2.HVPS_STATE)[iframe]) OF !LAMP_HVPS_A_ON : IF ((tab2.APD_STATE)[iframe] EQ !LAMP_APDOOR_OPEN) THEN $ flatimage = flatAO $ ELSE $ flatimage = flatAC !LAMP_HVPS_B_ON : IF ((tab2.APD_STATE)[iframe] EQ !LAMP_APDOOR_OPEN) THEN $ flatimage = flatBO $ ELSE $ flatimage = flatBC !LAMP_HVPS_AB_ON : IF ((tab2.APD_STATE)[iframe] EQ !LAMP_APDOOR_OPEN) THEN $ flatimage = flatABO $ ELSE $ flatimage = flatABC ENDCASE image /= flatimage error_image /= flatimage message = ' flatfield correction = YES' lamp_mike_log, message ENDIF ; effective area calibration IF (KEYWORD_SET(aeffflag)) THEN BEGIN CASE ((tab2.HVPS_STATE)[iframe]) OF !LAMP_HVPS_A_ON : IF ((tab2.APD_STATE)[iframe] EQ !LAMP_APDOOR_OPEN) THEN $ aeffimage = INTERPOL(aeffAO, aeffwave, waveimg) $ ELSE $ aeffimage = INTERPOL(aeffAC, aeffwave, waveimg) !LAMP_HVPS_B_ON : IF ((tab2.APD_STATE)[iframe] EQ !LAMP_APDOOR_OPEN) THEN $ aeffimage = INTERPOL(aeffBO, aeffwave, waveimg) $ ELSE $ aeffimage = INTERPOL(aeffBC, aeffwave, waveimg) !LAMP_HVPS_AB_ON : IF ((tab2.APD_STATE)[iframe] EQ !LAMP_APDOOR_OPEN) THEN $ aeffimage = INTERPOL(aeffABO, aeffwave, waveimg) $ ELSE $ aeffimage = INTERPOL(aeffABC, aeffwave, waveimg) ENDCASE image /= aeffimage error_image /= aeffimage message = ' effective area calibration = YES' lamp_mike_log, message ENDIF ; zero entries outside LAMP spectral range and set error = 1 out_of_range = WHERE((waveimg LT !LAMP_WAVELENGTH_LO) OR $ (waveimg GT !LAMP_WAVELENGTH_HI), count) IF (count GT 0) THEN BEGIN image[out_of_range] = 0.0 error_image[out_of_range] = 1.0 ENDIF ; store results in histogram cube k += 1 histogram_cube[*, *, 0, k] = image histogram_cube[*, *, 1, k] = error_image ENDELSE ENDFOR ; Divide dark_factor array by (60.0 * n_dark_elements) to get average counts/sec/element ; in our dark signal measuring regions dark_factor /= (60.0 * n_dark_elements) ; Divide stray_factor array by (60.0 * n_stray_elements) to get average counts/sec/element ; in our stray signal measuring regions stray_factor /= (60.0 * n_stray_elements) ; Compute wavelength and effective area calibration images for pixellist entries, ; but at the slow table cadence FOR i = 0L, nslow - 1L DO BEGIN waveimg = waveimage ; use default if STIM and HK data not available IF ((MAX(stim1_data[*, i]) GT 0.0) AND (MAX(stim2_data[*, i]) GT 0.0)) THEN BEGIN ; start of Andrew's code for stim1 stim1err = SQRT(stim1_data[*, i]) stim1tot = TOTAL(stim1_data[*, i]) stim1toterr = SQRT(stim1tot) stim1_index1 = stim1_data[*, i] * stim1x stim1err_index1 = stim1err * stim1x stim1indextot = totalerr(stim1_index1, stim1err_index1, err = stim1indextoterr) stimpos1 = diverr(stim1indextot, stim1indextoterr, stim1tot, stim1toterr, err = stimpos1err_i) stimpos1err = stimpos1err_i ; end of Andrew's code stim1pos = (MATRIX_MULTIPLY(stim1x, stim1_data[*, i], /ATRANSPOSE) / TOTAL(stim1_data[*, i]))[0] ; my value ; start of Andrew's code for stim2 stim2err = SQRT(stim2_data[*, i]) stim2tot = TOTAL(stim2_data[*, i]) stim2toterr = SQRT(stim2tot) stim2_index2 = stim2_data[*, i] * stim2x stim2err_index2 = stim2err * stim2x stim2indextot = totalerr(stim2_index2, stim2err_index2, err = stim2indextoterr) stimpos2 = diverr(stim2indextot, stim2indextoterr, stim2tot, stim2toterr, err = stimpos2err_i) stimpos2err = stimpos2err_i ; end of Andrew's code stim2pos = (MATRIX_MULTIPLY(stim2x, stim2_data[*, i], /ATRANSPOSE) / TOTAL(stim2_data[*, i]))[0] ; my value ; stim_drift1[i] = stim1pos ; stim_drift2[i] = stim2pos stim_drift1[i] = stimpos1 ; Andrew's value stim_drift2[i] = stimpos2 ; Andrew's value stimpos = [stim_drift1[i], stim_drift2[i]] stimposerr = [stimpos1err, stimpos2err] ; Andrew's values waveimg = lamp_mike_wavecal(STIMPOS = stimpos) ENDIF ELSE IF (good_housekeeping) THEN BEGIN result = MIN(ABS(tab6.SCUT_TIME - scut_time_slow[i]), idx) IF (result LE 1.0D0) THEN BEGIN ; HK entry within a second of desired time detetemp = (tab6.DET_HOUS_TMP)[idx] waveimg = lamp_mike_wavecal(DETETEMP = detetemp) ENDIF ENDIF slow_table_waveimg[*, *, i] = waveimg IF (KEYWORD_SET(aeffflag)) THEN BEGIN slow_table_aeffimg[*, *, i, 0, 0] = INTERPOL(aeffAO, aeffwave, waveimg) slow_table_aeffimg[*, *, i, 1, 0] = INTERPOL(aeffBO, aeffwave, waveimg) slow_table_aeffimg[*, *, i, 2, 0] = INTERPOL(aeffABO, aeffwave, waveimg) slow_table_aeffimg[*, *, i, 0, 1] = INTERPOL(aeffAC, aeffwave, waveimg) slow_table_aeffimg[*, *, i, 1, 1] = INTERPOL(aeffBC, aeffwave, waveimg) slow_table_aeffimg[*, *, i, 2, 1] = INTERPOL(aeffABC, aeffwave, waveimg) ENDIF ENDFOR ; Scan pixellist array and perform remaining radiometric calibrations ; PHD correction for Lyman alpha photons from Randy Gladstone's formula ac_lya = INTERPOL(phd_ac_lya, phd_scut_time, scut_time) ; accumulated Lyman alpha counts ( = fluence measure) det_hous_tmp = INTERPOL(tab6.DET_HOUS_TMP, tab6.SCUT_TIME, scut_time) ; phd_factor = -1.104 + 7.862E-9 * scut_time + 3.391E-3 * det_hous_tmp ; formula using time and temperature (second choice) phd_factor = 1.024 + 3.970E-11 * ac_lya + 3.258E-3 * det_hous_tmp ; formula using Lya fluence and temperature (preferable) phd_nonlya_indices = WHERE((detector_x LT !LAMP_LYA_X_MIN) OR (detector_x GT !LAMP_LYA_X_MAX)) phd_factor[phd_nonlya_indices] = 1.0 weighted_count *= FLOAT(phd_factor) ;; DEK DEBUG ; OPENW, luntmp, '/home/soc/tmp/phd.dat', /GET_LUN ; FOR j = 0L, nfovevents - 1L DO BEGIN ; PRINTF, luntmp, scut_time[j], phd_factor[j], detector_x[j], FORMAT = '(F20.5,F20.5,I5)' ; ENDFOR ; FREE_LUN, luntmp ;; DEK DEBUG FOR j = 0L, nfovevents - 1L DO BEGIN slow_idx = FLOOR((scut_time[j] - scut_scifile_start)/60.0D0) IF (slow_idx LT 0L) THEN slow_idx = 0L IF (slow_idx GT nslow - 1) THEN slow_idx = nslow - 1 ; assign spatial_row value, performing geometric distortion correction if requested spatial_row[j] = detector_y[j] ;; IF (KEYWORD_SET(rectflag)) THEN BEGIN ;; ENDIF ; wavelength wavelength[j] = slow_table_waveimg[detector_x[j], detector_y[j], slow_idx] ; flatfield correction IF (KEYWORD_SET(flatflag)) THEN BEGIN CASE (hvps_state[j]) OF !LAMP_HVPS_A_ON : IF (apdoor_state[j] EQ !LAMP_APDOOR_OPEN) THEN $ weighted_count[j] /= flatAO[detector_x[j], detector_y[j]] $ ELSE $ weighted_count[j] /= flatAC[detector_x[j], detector_y[j]] !LAMP_HVPS_B_ON : IF (apdoor_state[j] EQ !LAMP_APDOOR_OPEN) THEN $ weighted_count[j] /= flatBO[detector_x[j], detector_y[j]] $ ELSE $ weighted_count[j] /= flatBC[detector_x[j], detector_y[j]] !LAMP_HVPS_AB_ON : IF (apdoor_state[j] EQ !LAMP_APDOOR_OPEN) THEN $ weighted_count[j] /= flatABO[detector_x[j], detector_y[j]] $ ELSE $ weighted_count[j] /= flatABC[detector_x[j], detector_y[j]] ENDCASE ENDIF ; add to reduced countrate totals IF ((detector_x[j] GE !LAMP_FOV_X_MIN_CR) AND (detector_x[j] LE !LAMP_FOV_X_MAX_CR) AND (ncrentries GT 0L)) THEN $ cr_countrate[cr_index[j]] += weighted_count[j] ; effective area calibration IF (KEYWORD_SET(aeffflag)) THEN BEGIN CASE (hvps_state[j]) OF !LAMP_HVPS_A_ON : IF (apdoor_state[j] EQ !LAMP_APDOOR_OPEN) THEN $ weighted_count[j] /= slow_table_aeffimg[detector_x[j], detector_y[j], slow_idx, 0, 0] $ ELSE $ weighted_count[j] /= slow_table_aeffimg[detector_x[j], detector_y[j], slow_idx, 0, 1] !LAMP_HVPS_B_ON : IF (apdoor_state[j] EQ !LAMP_APDOOR_OPEN) THEN $ weighted_count[j] /= slow_table_aeffimg[detector_x[j], detector_y[j], slow_idx, 1, 0] $ ELSE $ weighted_count[j] /= slow_table_aeffimg[detector_x[j], detector_y[j], slow_idx, 1, 1] !LAMP_HVPS_AB_ON : IF (apdoor_state[j] EQ !LAMP_APDOOR_OPEN) THEN $ weighted_count[j] /= slow_table_aeffimg[detector_x[j], detector_y[j], slow_idx, 2, 0] $ ELSE $ weighted_count[j] /= slow_table_aeffimg[detector_x[j], detector_y[j], slow_idx, 2, 1] ENDCASE ENDIF ENDFOR ; Compute reduced count rates IF (ncrentries GT 0L) THEN BEGIN cr_dark_factor = INTERPOL(dark_factor, scut_time_slow, cr_scut_time) n_elements_fov_cr = (!LAMP_FOV_X_MAX_CR - !LAMP_FOV_X_MIN_CR + 1) * (!LAMP_FOV_Y_MAX - !LAMP_FOV_Y_MIN + 1) cr_dark_factor *= n_elements_fov_cr cr_countrate /= (hack_clock * cr_hack_incr) cr_countrate -= cr_dark_factor ; FOR m = 0L, ncrentries - 1L DO BEGIN ; cr_countrate[m] /= (hack_clock * cr_hack_incr[m]) ; cr_countrate[m] -= cr_dark_factor[m] ; ENDFOR ENDIF ; Fill in geometric quantities for any pixellist events IF (nfovevents GT 0L) THEN BEGIN IF (use_spice) THEN BEGIN ; Scan pixellist array and perform SPICE-based geometric calibrations (vectorized method) FOR y = !LAMP_FOV_Y_MIN, !LAMP_FOV_Y_MAX DO BEGIN ; Get array indices for this value of spatial row index1 = WHERE(spatial_row EQ y, count1) IF (count1 NE 0) THEN BEGIN ; Get ET from SCUT scut_seconds = FLOOR(scut_time[index1]) scut_fractions = ROUND(65536.0D0 * (scut_time[index1] - scut_seconds)) cspice_scs2e, !SPICE_LRO_ID, STRCOMPRESS(STRING(scut_seconds) + '.' + STRING(scut_fractions), /REMOVE_ALL), $ et ; Get altitude cspice_subpt, 'NEARPOINT', !SPICE_TARGET, et, !SPICE_ABCORR, !SPICE_OBSRVR, $ spoint, alt altitude[index1] = alt ; Allocate arrays for longitude and latitude values and illumination angles dlon = DBLARR(count1, /NOZERO) dlat = DBLARR(count1, /NOZERO) dinc = DBLARR(count1, /NOZERO) demi = DBLARR(count1, /NOZERO) ; Get latitude and longitude of pixel center, plus incidence and emission angles dvec = pixel_centers[*, y] ;; DEK: SPICE version N0063 dies on the call below with a system memory error; had to revert to N0062 to get it to work cspice_srfxpt, !SPICE_METHOD[0], !SPICE_TARGET, et, !SPICE_ABCORR, !SPICE_OBSRVR, !SPICE_SC_BUS_FRAME, dvec, $ spoint, dist, trgepc, obspos, found index2 = WHERE(found, count2, COMPLEMENT = index3, NCOMPLEMENT = count3) IF (count2 NE 0) THEN BEGIN cspice_reclat, spoint[*, index2], radius, lon, lat dlon[index2] = lon * !SPICE_DPR index4 = WHERE(dlon LT 0.0D0, count4) IF (count4 NE 0) THEN dlon[index4] += 360.0D0 dlat[index2] = lat * !SPICE_DPR cspice_illum, !SPICE_TARGET, et[index2], !SPICE_ABCORR, !SPICE_OBSRVR, spoint[*, index2], $ phase, solar, emissn dinc[index2] = solar * !SPICE_DPR demi[index2] = emissn * !SPICE_DPR ENDIF IF (count3 NE 0) THEN BEGIN dlon[index3] = !SPICE_BADANGLE dlat[index3] = !SPICE_BADANGLE dinc[index3] = !SPICE_BADANGLE demi[index3] = !SPICE_BADANGLE ENDIF longitude[index1] = dlon latitude[index1] = dlat incidence_angle[index1] = dinc emission_angle[index1] = demi ; Get latitude and longitude for lower left pixel corner dvec = pixel_corners[*, 64 - y] cspice_srfxpt, !SPICE_METHOD[0], !SPICE_TARGET, et, !SPICE_ABCORR, !SPICE_OBSRVR, !SPICE_SC_BUS_FRAME, dvec, $ spoint, dist, trgepc, obspos, found index2 = WHERE(found, count2, COMPLEMENT = index3, NCOMPLEMENT = count3) IF (count2 NE 0) THEN BEGIN cspice_reclat, spoint[*, index2], radius, lon, lat dlon[index2] = lon * !SPICE_DPR index4 = WHERE(dlon LT 0.0D0, count4) IF (count4 NE 0) THEN dlon[index4] += 360.0D0 dlat[index2] = lat * !SPICE_DPR ENDIF IF (count3 NE 0) THEN BEGIN dlon[index3] = !SPICE_BADANGLE dlat[index3] = !SPICE_BADANGLE ENDIF longitude_ll[index1] = dlon latitude_ll[index1] = dlat ; Get latitude and longitude for upper left pixel corner dvec = pixel_corners[*, 65 - y] cspice_srfxpt, !SPICE_METHOD[0], !SPICE_TARGET, et, !SPICE_ABCORR, !SPICE_OBSRVR, !SPICE_SC_BUS_FRAME, dvec, $ spoint, dist, trgepc, obspos, found index2 = WHERE(found, count2, COMPLEMENT = index3, NCOMPLEMENT = count3) IF (count2 NE 0) THEN BEGIN cspice_reclat, spoint[*, index2], radius, lon, lat dlon[index2] = lon * !SPICE_DPR index4 = WHERE(dlon LT 0.0D0, count4) IF (count4 NE 0) THEN dlon[index4] += 360.0D0 dlat[index2] = lat * !SPICE_DPR ENDIF IF (count3 NE 0) THEN BEGIN dlon[index3] = !SPICE_BADANGLE dlat[index3] = !SPICE_BADANGLE ENDIF longitude_ul[index1] = dlon latitude_ul[index1] = dlat ; Get latitude and longitude for upper right pixel corner dvec = pixel_corners[*, y] cspice_srfxpt, !SPICE_METHOD[0], !SPICE_TARGET, et, !SPICE_ABCORR, !SPICE_OBSRVR, !SPICE_SC_BUS_FRAME, dvec, $ spoint, dist, trgepc, obspos, found index2 = WHERE(found, count2, COMPLEMENT = index3, NCOMPLEMENT = count3) IF (count2 NE 0) THEN BEGIN cspice_reclat, spoint[*, index2], radius, lon, lat dlon[index2] = lon * !SPICE_DPR index4 = WHERE(dlon LT 0.0D0, count4) IF (count4 NE 0) THEN dlon[index4] += 360.0D0 dlat[index2] = lat * !SPICE_DPR ENDIF IF (count3 NE 0) THEN BEGIN dlon[index3] = !SPICE_BADANGLE dlat[index3] = !SPICE_BADANGLE ENDIF longitude_ur[index1] = dlon latitude_ur[index1] = dlat ; Get latitude and longitude for lower right pixel corner dvec = pixel_corners[*, y + 1] cspice_srfxpt, !SPICE_METHOD[0], !SPICE_TARGET, et, !SPICE_ABCORR, !SPICE_OBSRVR, !SPICE_SC_BUS_FRAME, dvec, $ spoint, dist, trgepc, obspos, found index2 = WHERE(found, count2, COMPLEMENT = index3, NCOMPLEMENT = count3) IF (count2 NE 0) THEN BEGIN cspice_reclat, spoint[*, index2], radius, lon, lat dlon[index2] = lon * !SPICE_DPR index4 = WHERE(dlon LT 0.0D0, count4) IF (count4 NE 0) THEN dlon[index4] += 360.0D0 dlat[index2] = lat * !SPICE_DPR ENDIF IF (count3 NE 0) THEN BEGIN dlon[index3] = !SPICE_BADANGLE dlat[index3] = !SPICE_BADANGLE ENDIF longitude_lr[index1] = dlon latitude_lr[index1] = dlat ENDIF ; count1 NE 0 ENDFOR ; y = !LAMP_FOV_Y_MIN, !LAMP_FOV_Y_MAX ENDIF ELSE BEGIN ; Obtain geometric quantities from RDR geomfile ; Read in calibrated pixel list mode data, extension 3, binary table geom_tab3 = mrdfits(geomfiles[ifile], 3, geom_hdr3, /UNSIGNED, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading geomfile extension 3 with mrdfits, status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Make sure the number of geomfile entries equals the number of current FOV events IF (N_ELEMENTS(geom_tab3.SCUT_TIME) NE nfovevents) THEN BEGIN message = 'Number mismatch between geomfile table 3 entries and current FOV events' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Assign the values altitude[0:nfovevents-1] = (geom_tab3.ALTITUDE)[0:nfovevents-1] incidence_angle[0:nfovevents-1] = (geom_tab3.INCIDENCE_ANGLE)[0:nfovevents-1] emission_angle[0:nfovevents-1] = (geom_tab3.EMISSION_ANGLE)[0:nfovevents-1] longitude[0:nfovevents-1] = (geom_tab3.LONGITUDE)[0:nfovevents-1] latitude[0:nfovevents-1] = (geom_tab3.LATITUDE)[0:nfovevents-1] longitude_ll[0:nfovevents-1] = (geom_tab3.LONGITUDE_LL)[0:nfovevents-1] latitude_ll[0:nfovevents-1] = (geom_tab3.LATITUDE_LL)[0:nfovevents-1] longitude_ul[0:nfovevents-1] = (geom_tab3.LONGITUDE_UL)[0:nfovevents-1] latitude_ul[0:nfovevents-1] = (geom_tab3.LATITUDE_UL)[0:nfovevents-1] longitude_ur[0:nfovevents-1] = (geom_tab3.LONGITUDE_UR)[0:nfovevents-1] latitude_ur[0:nfovevents-1] = (geom_tab3.LATITUDE_UR)[0:nfovevents-1] longitude_lr[0:nfovevents-1] = (geom_tab3.LONGITUDE_LR)[0:nfovevents-1] latitude_lr[0:nfovevents-1] = (geom_tab3.LATITUDE_LR)[0:nfovevents-1] ENDELSE ENDIF ; (nfovevents GT 0L) ; Compute and store remaining slow pixellist table entries IF (nslow GT 0L) THEN BEGIN IF (use_spice) THEN BEGIN FOR i = 0L, nslow - 1L DO BEGIN ; Get ET from SCUT scut_seconds = FLOOR(scut_time_slow[i]) scut_fractions = ROUND(65536.0D0 * (scut_time_slow[i] - scut_seconds)) cspice_scs2e, !SPICE_LRO_ID, STRCOMPRESS(STRING(scut_seconds, '.', scut_fractions), /REMOVE_ALL), et ; Get lunar surface crossing point of the LAMP boresight cspice_srfxpt, !SPICE_METHOD[0], !SPICE_TARGET, et, !SPICE_ABCORR, !SPICE_OBSRVR, !SPICE_SC_BUS_FRAME, bsight, $ spoint, dist, trgepc, obspos, found ; From this compute solar and S/C zenith angles and phase IF (found) THEN BEGIN cspice_illum, !SPICE_TARGET, et, !SPICE_ABCORR, !SPICE_OBSRVR, spoint, $ phase, solar, emissn azimuth_angle[i] = ACOS((COS(phase) - COS(solar)*COS(emissn))/(SIN(solar)*SIN(emissn))) * !SPICE_DPR ; azimuth angle computed via spherical trig identity, range = [0, PI] solar_zenith_angle[i] = solar * !SPICE_DPR sc_zenith_angle[i] = emissn * !SPICE_DPR phase_angle[i] = phase * !SPICE_DPR ENDIF ELSE BEGIN ; insert bad values to flag failure azimuth_angle[i] = !SPICE_BADANGLE solar_zenith_angle[i] = !SPICE_BADANGLE sc_zenith_angle[i] = !SPICE_BADANGLE phase_angle[i] = !SPICE_BADANGLE ENDELSE ; Get equatorial inertial frame coordinates of LAMP boresight and transform to RA and DEC cspice_pxform, !SPICE_SC_BUS_FRAME, !SPICE_EQUATORIAL_FRAME, et, rotate cspice_mxv, rotate, bsight, bsight_inertial cspice_recrad, bsight_inertial, radius, ra_tmp, dec_tmp ra[i] = ra_tmp * !SPICE_DPR dec[i] = dec_tmp * !SPICE_DPR ; Get the ecliptic coordinates of LRO wrt the solar system barycenter cspice_spkssb, !SPICE_LRO_ID, et, !SPICE_ECLIPTIC_FRAME, starg ; starg also returns velocity components ecliptic_x[i] = starg[0] ecliptic_y[i] = starg[1] ecliptic_z[i] = starg[2] ; Get the direction in the equatorial inertial frame from LRO to the Sun cspice_spkezp, !SPICE_SUN_ID, et, !SPICE_EQUATORIAL_FRAME, !SPICE_ABCORR, !SPICE_LRO_ID, $ ptarg, ltime cspice_srfxpt, !SPICE_METHOD[0], !SPICE_TARGET, et, !SPICE_ABCORR, !SPICE_OBSRVR, !SPICE_EQUATORIAL_FRAME, ptarg, $ spoint, dist, trgepc, obspos, found IF (found) THEN sc_shadow_flag[i] = 1 ; vector from LRO to Sun intersects the Moon, thus LRO is in shadow ; NOTE: this algorithm does NOT account for a finite-size Sun, so there is some inaccuracy near sunrise/sunset as seen from LRO ; Get sub-spacecraft latitude and longitude cspice_subpnt, !SPICE_METHOD[1], !SPICE_TARGET, et, !SPICE_MOON_FRAME, !SPICE_ABCORR, !SPICE_OBSRVR, $ spoint, trgepc, srfvec cspice_reclat, spoint, radius, lon, lat lon = lon * !SPICE_DPR lat = lat * !SPICE_DPR IF (lon LT 0.0D0) THEN lon = lon + 360.0D0 subsc_longitude[i] = lon subsc_latitude[i] = lat ; Get the boresight vector in the Moon frame and compute the angle between boresight and nadir cspice_pxform, !SPICE_SC_BUS_FRAME, !SPICE_MOON_FRAME, et, rotate cspice_mxv, rotate, bsight, bsight_moonfixed separation = cspice_vsep(srfvec, bsight_moonfixed) * !SPICE_DPR IF (separation GT !LAMP_NADIR_PTG_THRESH) THEN off_nadir_pointing[i] = 1 ; Get sub-solar latitude and longitude cspice_subslr, !SPICE_METHOD[1], !SPICE_TARGET, et, !SPICE_MOON_FRAME, !SPICE_ABCORR, !SPICE_OBSRVR, $ spoint, trgepc, srfvec cspice_reclat, spoint, radius, lon, lat lon = lon * !SPICE_DPR lat = lat * !SPICE_DPR IF (lon LT 0.0D0) THEN lon = lon + 360.0D0 subsolar_longitude[i] = lon subsolar_latitude[i] = lat ; Get Earth phase angle cspice_subpnt, !SPICE_METHOD[1], STRING(!SPICE_EARTH_ID), et, !SPICE_EARTH_FRAME, !SPICE_ABCORR, !SPICE_OBSRVR, $ spoint, trgepc, srfvec cspice_illum, STRING(!SPICE_EARTH_ID), et, !SPICE_ABCORR, !SPICE_OBSRVR, spoint, $ phase, solar, emissn earth_phase_angle[i] = phase * !SPICE_DPR ; Get LRO orbit number orbit_number[i] = -1 IF ((scut_time_slow[i] GE lro_orbit_scut[0]) && (scut_time_slow[i] LE lro_orbit_scut[n_lro_orbits-1])) THEN BEGIN IF (scut_time_slow[i] EQ lro_orbit_scut[n_lro_orbits-1]) THEN orbit_number[i] = lro_orbit_number[n_lro_orbits-1] ELSE BEGIN jhi = 1L WHILE (scut_time_slow[i] GE lro_orbit_scut[jhi]) DO jhi = jhi + 1L jlo = jhi - 1L IF ((scut_time_slow[i] LT lro_orbit_scut[jhi]) && (scut_time_slow[i] GE lro_orbit_scut[jlo]) && $ (lro_orbit_number[jhi] EQ (lro_orbit_number[jlo] + 1))) THEN orbit_number[i] = lro_orbit_number[jlo] ENDELSE ENDIF ENDFOR ; i = 0L, nslow - 1L ENDIF ELSE BEGIN ; Obtain geometric quantities from RDR geomfile ; Read in ancillary data, extension 4, binary table geom_tab4 = mrdfits(geomfiles[ifile], 4, geom_hdr4, /UNSIGNED, STATUS = status) IF (status LT 0) THEN BEGIN message = 'Error reading geomfile extension 4 with mrdfits, status = ' + STRING(status) lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Make sure the number of geomfile entries equals the number of expected slow table entries IF (N_ELEMENTS(geom_tab4.SCUT_TIME_SLOW) NE nslow) THEN BEGIN message = 'Number mismatch between geomfile table 4 entries and expected slow table entries' lamp_mike_log, message, ERRNUM = 2, VERBOSELEVEL = verbose error_status = 1 GOTO, END_OF_LOOP ENDIF ; Assign the values azimuth_angle[0:nslow-1] = (geom_tab4.AZIMUTH_ANGLE)[0:nslow-1] solar_zenith_angle[0:nslow-1] = (geom_tab4.SOLAR_ZENITH_ANGLE)[0:nslow-1] sc_zenith_angle[0:nslow-1] = (geom_tab4.SC_ZENITH_ANGLE)[0:nslow-1] phase_angle[0:nslow-1] = (geom_tab4.PHASE_ANGLE)[0:nslow-1] ra[0:nslow-1] = (geom_tab4.RA)[0:nslow-1] dec[0:nslow-1] = (geom_tab4.DEC)[0:nslow-1] ecliptic_x[0:nslow-1] = (geom_tab4.ECLIPTIC_X)[0:nslow-1] ecliptic_y[0:nslow-1] = (geom_tab4.ECLIPTIC_Y)[0:nslow-1] ecliptic_z[0:nslow-1] = (geom_tab4.ECLIPTIC_Z)[0:nslow-1] sc_shadow_flag[0:nslow-1] = (geom_tab4.SC_SHADOW_FLAG)[0:nslow-1] subsc_longitude[0:nslow-1] = (geom_tab4.SUBSC_LONGITUDE)[0:nslow-1] subsc_latitude[0:nslow-1] = (geom_tab4.SUBSC_LATITUDE)[0:nslow-1] off_nadir_pointing[0:nslow-1] = (geom_tab4.OFF_NADIR_POINTING)[0:nslow-1] subsolar_longitude[0:nslow-1] = (geom_tab4.SUBSOLAR_LONGITUDE)[0:nslow-1] subsolar_latitude[0:nslow-1] = (geom_tab4.SUBSOLAR_LATITUDE)[0:nslow-1] earth_phase_angle[0:nslow-1] = (geom_tab4.EARTH_PHASE_ANGLE)[0:nslow-1] orbit_number[0:nslow-1] = (geom_tab4.ORBIT_NUMBER)[0:nslow-1] ENDELSE ENDIF ; (nslow GT 0L) ; Prepare data and headers for writing to RDR FITS ; Calibrated door open histogram, primary extension, image im0_rdr = FLTARR(1024, 32) im1_rdr = im0_rdr ; generate image for door closed histogram at same time IF (nfovevents GT 0L) THEN BEGIN ; use pixellist data preferentially if available FOR j = 0L, nfovevents - 1L DO BEGIN IF (apdoor_state[j] EQ !LAMP_APDOOR_OPEN) THEN $ im0_rdr[detector_x[j], spatial_row[j]] += weighted_count[j] $ ELSE IF (apdoor_state[j] EQ !LAMP_APDOOR_CLOSED) THEN $ im1_rdr[detector_x[j], spatial_row[j]] += weighted_count[j] ENDFOR ENDIF ELSE IF (nhistograms GT 0L) THEN BEGIN ; use histograms next if available im0_rdr = histogram_cube[*, *, 0, 0] IF (nhistograms GT 1L) THEN im1_hdr = histogram_cube[*, *, 0, 1] ENDIF ELSE BEGIN ; reuse EDR images as a last resort when no other data are present im0_rdr = im0 im1_rdr = im1 ENDELSE ; reuse EDR header with the following modifications sxaddpar, hdr0, 'BITPIX' , -32 sxaddpar, hdr0, 'COMMENT' , 'LAMP Level 2 (CODMAC Level 3) processing information ------------', AFTER = 'DOHVPS' sxaddpar, hdr0, 'SL2PROC' , !MIKE_SL2PROC, 'SOC pipeline level 2 software' sxaddpar, hdr0, 'SL2VER' , !MIKE_SL2VER, 'SOC pipeline level 2 software version' sxaddpar, hdr0, 'SL2DATE' , !MIKE_SL2DATE, 'SOC pipeline level 2 software date' sxaddpar, hdr0, 'FILE_IN' , FILE_BASENAME(infiles[ifile]), 'Input file for processing' sxaddpar, hdr0, 'FILE_OUT', FILE_BASENAME(outfiles[ifile]), 'Output file after processing' sxaddpar, hdr0, 'DIR_CAL' , calibration_dir, 'Directory for calibration data' sxaddpar, hdr0, 'BADFILE' , KEYWORD_SET(badfile) ? FILE_BASENAME(badfile) : '', 'File containing bad pixel mask' sxaddpar, hdr0, 'BADFLAG' , KEYWORD_SET(badflag) ? 'T' : 'F', 'Bad pixel mask applied to data' sxaddpar, hdr0, 'BADVALUE', 'IEEE NaN', 'Bad pixel value' sxaddpar, hdr0, 'DEADFLAG', KEYWORD_SET(deadflag) ? 'T' : 'F', 'Dead-time correction applied' sxaddpar, hdr0, 'DARKFILE', KEYWORD_SET(darkfile) ? FILE_BASENAME(darkfile) : '', 'Detector dark count rate file' sxaddpar, hdr0, 'DARKFLAG', KEYWORD_SET(darkflag) ? 'T' : 'F', 'Dark count correction applied' sxaddpar, hdr0, 'FLATFILE', KEYWORD_SET(flatfile) ? FILE_BASENAME(flatfile) : '', 'Detector flatfield file' sxaddpar, hdr0, 'FLATFLAG', KEYWORD_SET(flatflag) ? 'T' : 'F', 'Flatfield correction applied' sxaddpar, hdr0, 'WAVEPRO' , 'lamp_mike_wavecal.pro', 'IDL wavelength calibration procedure' sxaddpar, hdr0, 'WAVEFLAG', 'T', 'Wavelength calibration performed' sxaddpar, hdr0, 'AEFFFILE', KEYWORD_SET(aefffile) ? FILE_BASENAME(aefffile) : '', 'Detector effective area file' sxaddpar, hdr0, 'AEFFFLAG', KEYWORD_SET(aeffflag) ? 'T' : 'F', 'Effective area correction applied' sxaddpar, hdr0, 'LOG_FILE', KEYWORD_SET(logfile) ? FILE_BASENAME(logfile) : '', 'LAMP-MIKE processing logfile' sxaddpar, hdr0, 'MIKE_ERR', error_status, 'Exit status for LAMP-MIKE' sxaddpar, hdr0, 'COMMENT' , 'SPICE files used in geometry calculations -----------------------', AFTER = 'MIKE_ERR' IF (use_spice) THEN BEGIN sxaddpar, hdr0, 'NUMBSPCK', nkernels, 'Number of loaded SPICE kernels' FOR k = 0, nkernels - 1 DO BEGIN strkp1 = STRING(k + 1, FORMAT = '(I04)') sxaddpar, hdr0, 'SPCK' + strkp1, kernel_list[k], 'SPICE kernel ' + strkp1 ENDFOR ENDIF ELSE BEGIN nkernels = sxpar(geom_hdr0, 'NUMBSPCK') sxaddpar, hdr0, 'NUMBSPCK', nkernels, 'Number of loaded SPICE kernels' FOR k = 0, nkernels - 1 DO BEGIN strkp1 = STRING(k + 1, FORMAT = '(I04)') kernel_name = sxpar(geom_hdr0, 'SPCK' + strkp1) sxaddpar, hdr0, 'SPCK' + strkp1, kernel_name, 'SPICE kernel ' + strkp1 ENDFOR ENDELSE ; Calibrated door closed histogram, extension 1, image ; image generated above with primary extension image ; reuse EDR header with the following modifications sxaddpar, hdr1, 'BITPIX', -32 sxaddpar, hdr1, 'EXTNAME', 'Calibrated Spectral Image Door Closed' ; Acquisition list, extension 2, ASCII table ; no modifications of EDR required, straight copy ; Calibrated pixel list mode data (fast), extension 3, binary table ; no modifications of EDR, this is a new RDR extension ; just write out new bintable below ; Calibrated pixel list mode data (slow), extension 4, binary table ; no modifications of EDR, this is a new RDR extension ; just write out new bintable below ; Calibrated histogram mode data, extension 5, image (stacked image cube) mkhdr, hdr5_rdr, '', /IMAGE sxaddpar, hdr5_rdr, 'XTENSION', 'IMAGE', ' Extension 5: Calibrated Histogram Mode Data' sxaddpar, hdr5_rdr, 'BITPIX', -32 sxaddpar, hdr5_rdr, 'BZERO', 0 sxaddpar, hdr5_rdr, 'BSCALE', 1 sxaddpar, hdr5_rdr, 'EXTNAME', 'Calibrated Histogram Mode Data' sxaddpar, hdr5_rdr, 'EXTVER', 1, ' Extension version number' ; Calibrated calculated countrate, extension 6, binary table ; IF (ncrentries GT 0L) THEN BEGIN ; hdr6_rdr = hdr4 ; sxaddpar, hdr6_rdr, 'EXTNAME', 'Calibrated Calculated Countrate' ; sxaddpar, hdr6_rdr, 'NAXIS2', ncrentries ; ENDIF ; High-rate LTS data, extension 7, binary table ; no modifications of EDR required, straight copy ; Housekeeping table, extension 8, binary table ; no modifications of EDR required, straight copy ; Wavelength lookup image, extension 9, image mkhdr, hdr9_rdr, waveimage, /IMAGE sxaddpar, hdr9_rdr, 'XTENSION', 'IMAGE', ' Extension 9: Wavelength Lookup Image' sxaddpar, hdr9_rdr, 'EXTNAME' , 'Wavelength Lookup Image' sxaddpar, hdr9_rdr, 'EXTVER' , 1, ' Extension version number' ; Write the RDR FITS mwrfits, im0_rdr, outfiles[ifile], hdr0, /CREATE ; calibrated door open histogram, primary extension, image mwrfits, im1_rdr, outfiles[ifile], hdr1 ; calibrated door closed histogram, extension 1, image mwrfits, tab2 , outfiles[ifile], hdr2, /ASCII ; acquisition list, extension 2, ASCII table ; calibrated pixel list mode data (fast), extension 3, binary table fxbhmake, hdr3_rdr, nfovevents, 'Calibrated Pixel List Mode Data', EXTVER = 1 fxbaddcol, 1 , hdr3_rdr, hack_time[0] , 'HACK_TIME' fxbaddcol, 2 , hdr3_rdr, detector_x[0] , 'DETECTOR_X' fxbaddcol, 3 , hdr3_rdr, detector_y[0] , 'DETECTOR_Y' fxbaddcol, 4 , hdr3_rdr, spatial_row[0] , 'SPATIAL_ROW' fxbaddcol, 5 , hdr3_rdr, wavelength[0] , 'WAVELENGTH' fxbaddcol, 6 , hdr3_rdr, scut_time[0] , 'SCUT_TIME' fxbaddcol, 7 , hdr3_rdr, weighted_count[0] , 'WEIGHTED_COUNT' fxbaddcol, 8 , hdr3_rdr, weighted_variance[0], 'WEIGHTED_VARIANCE' fxbaddcol, 9 , hdr3_rdr, latitude[0] , 'LATITUDE' fxbaddcol, 10, hdr3_rdr, longitude[0] , 'LONGITUDE' fxbaddcol, 11, hdr3_rdr, latitude_ll[0] , 'LATITUDE_LL' fxbaddcol, 12, hdr3_rdr, longitude_ll[0] , 'LONGITUDE_LL' fxbaddcol, 13, hdr3_rdr, latitude_ul[0] , 'LATITUDE_UL' fxbaddcol, 14, hdr3_rdr, longitude_ul[0] , 'LONGITUDE_UL' fxbaddcol, 15, hdr3_rdr, latitude_ur[0] , 'LATITUDE_UR' fxbaddcol, 16, hdr3_rdr, longitude_ur[0] , 'LONGITUDE_UR' fxbaddcol, 17, hdr3_rdr, latitude_lr[0] , 'LATITUDE_LR' fxbaddcol, 18, hdr3_rdr, longitude_lr[0] , 'LONGITUDE_LR' fxbaddcol, 19, hdr3_rdr, altitude[0] , 'ALTITUDE' fxbaddcol, 20, hdr3_rdr, incidence_angle[0] , 'INCIDENCE_ANGLE' fxbaddcol, 21, hdr3_rdr, emission_angle[0] , 'EMISSION_ANGLE' fxbaddcol, 22, hdr3_rdr, deadtime_factor[0] , 'DEADTIME_FACTOR' fxbaddcol, 23, hdr3_rdr, pix_data_qual[0] , 'PIX_DATA_QUAL' fxbcreate, lun, outfiles[ifile], hdr3_rdr IF (nfovevents GT 0L) THEN BEGIN column_names = ['HACK_TIME','DETECTOR_X','DETECTOR_Y','SPATIAL_ROW','WAVELENGTH','SCUT_TIME','WEIGHTED_COUNT', $ 'WEIGHTED_VARIANCE','LATITUDE','LONGITUDE','LATITUDE_LL','LONGITUDE_LL','LATITUDE_UL','LONGITUDE_UL',$ 'LATITUDE_UR','LONGITUDE_UR','LATITUDE_LR','LONGITUDE_LR','ALTITUDE','INCIDENCE_ANGLE', $ 'EMISSION_ANGLE','DEADTIME_FACTOR','PIX_DATA_QUAL'] fxbwritm, lun, column_names, hack_time, detector_x, detector_y, spatial_row, wavelength, scut_time, weighted_count, $ weighted_variance, latitude, longitude, latitude_ll, longitude_ll, latitude_ul, $ longitude_ul, latitude_ur, longitude_ur, latitude_lr, longitude_lr, altitude, $ incidence_angle, emission_angle, deadtime_factor, pix_data_qual, BUFFERSIZE=1048576L ENDIF fxbfinish, lun ; calibrated pixel list mode data (slow), extension 4, binary table fxbhmake, hdr4_rdr, nslow, 'Ancillary Data', EXTVER = 1 fxbaddcol, 1 , hdr4_rdr, scut_time_slow[0] , 'SCUT_TIME_SLOW' fxbaddcol, 2 , hdr4_rdr, orbit_number[0] , 'ORBIT_NUMBER' fxbaddcol, 3 , hdr4_rdr, solar_zenith_angle[0], 'SOLAR_ZENITH_ANGLE' fxbaddcol, 4 , hdr4_rdr, sc_shadow_flag[0] , 'SC_SHADOW_FLAG' fxbaddcol, 5 , hdr4_rdr, subsolar_latitude[0] , 'SUBSOLAR_LATITUDE' fxbaddcol, 6 , hdr4_rdr, subsolar_longitude[0], 'SUBSOLAR_LONGITUDE' fxbaddcol, 7 , hdr4_rdr, subsc_latitude[0] , 'SUBSC_LATITUDE' fxbaddcol, 8 , hdr4_rdr, subsc_longitude[0] , 'SUBSC_LONGITUDE' fxbaddcol, 9 , hdr4_rdr, sc_zenith_angle[0] , 'SC_ZENITH_ANGLE' fxbaddcol, 10, hdr4_rdr, earth_phase_angle[0] , 'EARTH_PHASE_ANGLE' fxbaddcol, 11, hdr4_rdr, phase_angle[0] , 'PHASE_ANGLE' fxbaddcol, 12, hdr4_rdr, azimuth_angle[0] , 'AZIMUTH_ANGLE' fxbaddcol, 13, hdr4_rdr, ra[0] , 'RA' fxbaddcol, 14, hdr4_rdr, dec[0] , 'DEC' fxbaddcol, 15, hdr4_rdr, ecliptic_x[0] , 'ECLIPTIC_X' fxbaddcol, 16, hdr4_rdr, ecliptic_y[0] , 'ECLIPTIC_Y' fxbaddcol, 17, hdr4_rdr, ecliptic_z[0] , 'ECLIPTIC_Z' fxbaddcol, 18, hdr4_rdr, stim_drift1[0] , 'STIM_DRIFT1' fxbaddcol, 19, hdr4_rdr, stim_drift2[0] , 'STIM_DRIFT2' fxbaddcol, 20, hdr4_rdr, other_inst_interf[0] , 'OTHER_INST_INTERF' fxbaddcol, 21, hdr4_rdr, global_data_qual[0] , 'GLOBAL_DATA_QUAL' fxbaddcol, 22, hdr4_rdr, heuristic_quality[0] , 'HEURISTIC_QUALITY' fxbaddcol, 23, hdr4_rdr, dark_factor[0] , 'DARK_FACTOR' fxbaddcol, 24, hdr4_rdr, stray_factor[0] , 'STRAY_FACTOR' fxbaddcol, 25, hdr4_rdr, off_nadir_pointing[0], 'OFF_NADIR_POINTING' fxbcreate, lun, outfiles[ifile], hdr4_rdr IF (nslow GT 0L) THEN BEGIN column_names = ['SCUT_TIME_SLOW','ORBIT_NUMBER','SOLAR_ZENITH_ANGLE','SC_SHADOW_FLAG','SUBSOLAR_LATITUDE', $ 'SUBSOLAR_LONGITUDE','SUBSC_LATITUDE','SUBSC_LONGITUDE','SC_ZENITH_ANGLE','EARTH_PHASE_ANGLE', $ 'PHASE_ANGLE','AZIMUTH_ANGLE','RA','DEC','ECLIPTIC_X','ECLIPTIC_Y','ECLIPTIC_Z','STIM_DRIFT1', $ 'STIM_DRIFT2','OTHER_INST_INTERF','GLOBAL_DATA_QUAL','HEURISTIC_QUALITY','DARK_FACTOR', $ 'STRAY_FACTOR','OFF_NADIR_POINTING'] fxbwritm, lun, column_names, scut_time_slow, orbit_number, solar_zenith_angle, sc_shadow_flag, $ subsolar_latitude, subsolar_longitude, subsc_latitude, subsc_longitude, $ sc_zenith_angle, earth_phase_angle, phase_angle, azimuth_angle, ra, dec, $ ecliptic_x, ecliptic_y, ecliptic_z, stim_drift1, stim_drift2, other_inst_interf, $ global_data_qual, heuristic_quality, dark_factor, stray_factor, off_nadir_pointing, BUFFERSIZE=1048576L ENDIF fxbfinish, lun mwrfits, histogram_cube, outfiles[ifile], hdr5_rdr ; calibrated histogram mode data, extension 5, image cube ; calibrated calculated countrate, extension 6, binary table IF (ncrentries GT 0L) THEN BEGIN fxbhmake, hdr6_rdr, ncrentries, 'Reduced Count Rate', EXTVER = 1 fxbaddcol, 1, hdr6_rdr, cr_hack_time[0] , 'HACK_TIME' fxbaddcol, 2, hdr6_rdr, cr_scut_time[0] , 'SCUT_TIME' fxbaddcol, 3, hdr6_rdr, LONG(cr_countrate[0]), 'COUNT_RATE' fxbcreate, lun, outfiles[ifile], hdr6_rdr column_names = ['HACK_TIME','SCUT_TIME','COUNT_RATE'] fxbwritm, lun, column_names, cr_hack_time, cr_scut_time, LONG(cr_countrate), BUFFERSIZE=1048576L fxbfinish, lun ENDIF ELSE $ ; if we have no calibrated counts, just write out old EDR extension mwrfits, tab4, outfiles[ifile], hdr4 ; calculated countrate, extension 6, binary table mwrfits, tab5, outfiles[ifile], hdr5 ; high-rate LTS data, extension 7, binary table mwrfits, tab6, outfiles[ifile], hdr6 ; housekeeping table, extension 8, binary table mwrfits, waveimage, outfiles[ifile], hdr9_rdr ; wavelength lookup image, extension 9, image END_OF_LOOP: ; Check for non-zero error status before ending loop. ; If so, then an error has occurred and the calibrated file could not be created. IF (error_status NE 0) THEN BEGIN HELP, CALLS = calls lamp_mike_log, '%% Fatal error at ' + calls, ERRNUM = 2 IF (STRLEN(!ERROR_STATE.MSG) GT 0) THEN $ lamp_mike_log, !ERROR_STATE.MSG, ERRNUM = 2 lamp_mike_log, 'File ' + infiles[ifile] + ' not processed!', ERRNUM = 2 ENDIF ; Write log messages for this file to output status file, if requested IF (KEYWORD_SET(mikeflags.statflag)) THEN BEGIN istat_last = N_ELEMENTS(logmessages) - 1 OPENW, lun, statfiles[ifile], /GET_LUN IF (istat_last GE istat_first) THEN PRINTF, lun, logmessages[istat_first : istat_last] FREE_LUN, lun ENDIF ENDFOR ; loop over input files ; Program clean-up END_OF_PROGRAM: ; Unload spice kernels cspice_unload, metakernel ; Compute wall clock time elapsed during processing stop_time = SYSTIME(1); seconds since January 1, 1970 elapsed_time = stop_time - start_time; seconds elapsed_time = STRTRIM(STRING(elapsed_time/60.0, FORMAT = '(F20.2)'), 2); minutes message = 'Total processing time = ' + elapsed_time + ' minutes.' lamp_mike_log, message, VERBOSELEVEL = 1 ; Write internal log to file IF (KEYWORD_SET(logfile)) THEN BEGIN OPENW, lun, logfile, /GET_LUN PRINTF, lun, logmessages FREE_LUN, lun ENDIF END ;+ ; NAME: ; lamp_mike_wavecal ; ; PURPOSE: ; Calculate the wavelength at each pixel location in a LAMP data file. ; ; CATEGORY: ; Data processing ; ; CALLING SEQUENCE: ; wavelength = lamp_mike_wavecal(PREFLIGHT=preflight,STIMPOS=stimpos,DETETEMP=detetemp) ; ; INPUTS: ; preflight -- use nominal pre-flight wavelength solution ; stimpos = [stim1, stim2], where ; stim1 -- position (centroid) of the left stim pulse ; stim2 -- position (centroid) of the right stim pulse, ; detetemp -- detector temperature ; ; OUTPUTS: ; Returns a 1024x32 element array containing the wavelength (in Angstroms) ; at each pixel. ; ; COMMON BLOCKS: ; None. ; ; SIDE EFFECTS: ; None. ; ; RESTRICTIONS: ; None. ; ; NOTES: ; - Adaptation of Andrew Steffl's palice_wavecal ; - As of v4, there is no independent solution that uses only the ; detector housing temperature; if the user specifies only the ; DETETEMP keyword (with any value), the routine will return the ; nominal pre-flight wavelength solution, just as if the user ; did not specify either keyword (STIMPOS or DETETEMP) ; ; MODIFICATION HISTORY: ; v1: 2007-03-07 AJS ; v2: 2007-09-05 AJS Switched to a grating equation-based ; non-linear solution for the wavelength. The constants c1 and ; c2 are derived from a fit to a spectrum of the Io Plasma Torus ; obtained in ali_0034698959_0x4b2_eng_1.fit ; v3: 2008-08-20 DEK Adapted for use with LAMP ; v4: 2009-11-06 DEK Updated to use Andrew's LAMP wavelength ; solution as determined by in-flight commissioning data. ;- FUNCTION lamp_mike_wavecal, PREFLIGHT=preflight, STIMPOS=stimpos, DETETEMP=detetemp IF (KEYWORD_SET(preflight)) THEN BEGIN ; Use Mike Davis' linear fit lambda (Angstroms) = 1.82798 * X (detector position) + 326.13241 slope = 1.82798 intercept = 326.13241 lambda = FINDGEN(1024) lambda = slope * lambda + intercept lambda #= (FLTARR(32) + 1.0) RETURN, lambda ENDIF ; If both the stim position and T keywords are set, use the stim ; position to determine the wavelength scale correction. IF (KEYWORD_SET(stimpos) AND KEYWORD_SET(detetemp)) THEN detetemp = 0 ; Grating equation: ; sin(alpha) - sin(beta) = lambda / d ; 2 * R * beta = c1 * x + c2 c1 = -0.043573533D0 ; mm/pixel c2 = 32.867475D0 ; mm dtor = !DPI/180.0D0 x = FINDGEN(1024) d = 6250.0 ; Angstroms (from 1600 lines/mm) alpha = 15.819 * dtor R = 75.0 ; mm (rowland circle) lambda = d * (sin(alpha) - sin((c1 * x + c2) / (2.0 * R))) ; Use of the detector housing temperature as the basis for the ; temperature correction is TBD. ;IF (KEYWORD_SET(detetemp)) THEN BEGIN ;; correction TBD ;ENDIF ; Use the observed position of the stim pulses to derive the ; temperature correction. This is the preferred way of doing things. IF (KEYWORD_SET(stimpos)) THEN BEGIN x_correction = lamp_mike_temp_correction(stimpos) lambda = INTERPOL(lambda, x, x_correction) ENDIF ; As of v4, there is no correction for the wavelength solution as a ; function of row. Since the slit is slightly bow-shaped (i.e. not ; straight) this might not be a good assumption. lambda = FLOAT(lambda) lambda #= (FLTARR(32) + 1.) RETURN, lambda END ;+ ; NAME: ; poisson_errorbar ; ; PURPOSE: ; For a given number of counts, determine the statistical uncertainty delta, ; based on Poisson statistics, such that the interval of: N - delta <-> N + delta ; has a confidence level of X percent (default is 68.3%) ; ; CATEGORY: ; statistics, error propagation ; ; CALLING SEQUENCE: ; error = poisson_errorbar(counts, [CONFIDENCE_LIMIT=confidence_limit,[DOUBLE=double]]) ; ; INPUTS: ; n_events -- number of detetcted counts ; ; KEYWORD PARAMETERS: ; confidence_limit -- the desired level of confidence for the interval ; double -- set this keyword to force double precision calculations ; ; OUTPUTS: ; returns the error bar associated with the specified level of confidence ; ; COMMON BLOCKS: ; poisson_errorbar_common -- contains variables cl, the specified confidence ; level and n, the number of counts ; ; RESTRICTIONS: ; Since the code relies on FX_ROOT to find the root of confidence limits ; function, the n_events must be a scalar ; ; MODIFICATION HISTORY: ; 2008-11-24 V1. AJS Initial version ;- FUNCTION confidence_limits_func, delta_n COMMON poisson_errorbar_common, cl, n upper_limit = n + delta_n lower_limit = n - delta_n ; These are based on Eq. 4 & 5 of Geherls 1986 ApJ 303: 336-346 prob_above = 1. - CHISQR_PDF(2.*upper_limit, 2.*n + 2.) prob_below = CHISQR_PDF(2.*lower_limit, 2.*n) RETURN, prob_above + prob_below + cl - 1. END FUNCTION poisson_errorbar, n_events, CONFIDENCE_LIMIT=confidence_limit, DOUBLE=double FORWARD_FUNCTION confidence_limits_func COMMON poisson_errorbar_common, cl, n IF (N_PARAMS() NE 1) THEN BEGIN PRINT, '% POISSON_ERRORBAR: errorbar = poisson_errorbar(n_events, CONFIDENCE_LIMIT=confidence_limit)' RETURN, -1 ENDIF IF (N_ELEMENTS(n_events) NE 1) THEN BEGIN PRINT, '% POISSON_ERRORBAR: n_events must be a scalar.' PRINT, 'Returning...' RETURN, -1 ENDIF ; Default level of confidence is one sigma (68.3%) IF (NOT KEYWORD_SET(confidence_limit)) THEN confidence_limit = GAUSS_PDF(1.) - GAUSS_PDF(-1.) ; Use double-precision if requested IF (KEYWORD_SET(double)) THEN n = DOUBLE(n_events) ELSE n = n_events cl = confidence_limit ; Use the IDL routine FX_ROOT to solve for the Poissonian error bar errorbar = FX_ROOT([(SQRT(n)-1.) > .01, SQRT(n) > .1, SQRT(n)+1.], 'confidence_limits_func') RETURN, errorbar END ; dimension is defined the same way as for the total.pro function function totalerr,array,errarray,dim,err=toterrarr,nan=nan nan = 1;why wouldn't you want to just ignore NaN's instead of returning a NaN total? errarr=errarray notfinite=where(finite(array,/nan),cnt) if cnt gt 0 then errarr[where(finite(array) eq 0)]=!values.f_nan if n_params() eq 3 then begin totarray=total(array,dim,nan=nan) finitetotal=total(finite(array),dim) notfinite=where(finitetotal eq 0,nancnt) toterrarr=sqrt(total(errarr^2,dim,nan=nan)) if nancnt gt 0 then begin totarray[notfinite]=!values.f_nan toterrarr[notfinite]=!values.f_nan endif endif else begin totarray=total(array,nan=nan) toterrarr=sqrt(total(errarr^2,nan=nan)) if total(finite(array)) eq 0 then begin totarray=!values.fnan toterrarray=!values.fnan endif endelse return,totarray end