pro mkmulti, input_file, help=help, norefsub=norefsub, numsamrd=numsamrd, numutr=numutr, frmtime=frmtime, output_file=output_file, mimic_multiaccum=mimic_multiaccum, sca_type=sca_type ;; NAME: ;; pro mkmulti.pro ;; ;; PURPOSE: ;; This program calibrates MULTIACCUM exposures. Calibration ;; includes fitting the slope and performing reference pixel correction. ;; ;; CALLING SEQUENCE: ;; mkmulti, input_file ;; ;; INPUTS: ;; input_file - A FITS data cube. The file must use IDTL header keywords. ;; numutr - Number of up the ramp sample groups. ;; frmtime - Time required to sample one frame in seconds ;; output_file - The name of the output file ;; sca_type - Set this to specify the type of ;; SCA on the command line ;; ;; ;; KEYWORD PARAMETERS: ;; help - Print a help message ;; norefsub - Do not perform reference pixel correction ;; mimic_multiaccum - Mimic MULTIACCUM readout. This updates ;; NUMSAMRD and NUMUTR to the values supplied ;; on the command line (if supplied) and ;; adjusts the exposure time. ;; ;; EXAMPLE ;; ;; mkmulti, my_multiaccum_image, my_resultant_image ;; ;; ;; REFERENCE: ;; ;; ;; MODIFICATION HISTORY: ;; Written by: ;; Modified: B.J. Rauscher, IDTL, May 6, 2002 ;; Initial release ;; ;; B.J. Rauscher, 21 January 2003 ;; --- Major revisions. Automated calibration of general images ;; ;; B.J. Rauscher, 17 February 2003 ;; --- Read large fits files by slice. This avoids the problem ;; of running out of memory on large files. ;; --- Add ability to mimic multiaccum exposures ;; ;; B.J.Rauscher, 8 April 2003 ;; -- Add on the fly reference pixel correction ;; ;; B.J. Rauscher, 15 April 2003 ;; -- Added sca_type parameter to command line ;;- ;; ;; This version of mkmulti version = '2.2' ;; Print a help message if requested if (not keyword_set(help)) then help = 0 if (help) then begin print, 'MKMULTI version '+strtrim(string(version),2) stop, 'USAGE: mkmulti, input_file [, /help, /norefsub, numsamrd=, numutr=, frmtime=, sca_type=, /mimic_multiaccum]' endif ;; Set keyword defauts if (not keyword_set(norefsub)) then norefsub = 0; Default is to perform reference pixel correction if (not keyword_set(mimic_multiaccum)) then mimic_multiaccum = 0; Default is not to mimic ; multiaccum ;; Open the FITS image and read the header. We will be reading ;; the FITS file slice-by-slice. Close it now so that the next ;; slice will work fits_open, input_file, fcb fits_read, fcb, data, header, /header_only fits_close, fcb ;; Load variables using header information print, 'Reading header information' xsize = sxpar(header, 'NAXIS1', 'Keyword NAXIS1 not found') ; x-dimension of image ysize = sxpar(header, 'NAXIS2', 'Keyword NAXIS2 not found') ; y-dimension of image zsize = sxpar(header, 'NAXIS3', 'Keyword NAXIS3 not found') ; z-dimension of image exp_time = sxpar(header,'EXP_TIME', 'Keyword EXP_TIME not found') ;; If user has not specified the sca_type, get it from the header if (not keyword_set(sca_type)) then begin sca_type = '' ; Create an empty string to hold the SCA type getscatype, header, sca_type endif ;; Try to get the following variables from the command ;; line. If not found there, try to get them from the header ;; NUMSAMRD is the number of samples per group if (not keyword_set(numsamrd)) then $ numsamrd = sxpar(header, 'NUMSAMRD', 'Keyword NUMSAMRD not found') ;; FRMTIME is the time required to read one frame ;; This is used only when mimicing MULTIACCUM if (mimic_multiaccum) then begin if (not keyword_set(frmtime)) then begin frmtime = sxpar(header, 'FRMTIME', 'Keyword FRMTIME not found') endif endif ;; NUMUTR is the number of groups going up the ramp if (not keyword_set(numutr)) then $ numutr = sxpar(header, 'NUMUTR', 'Keyword NUMUTR not found') ;; Set up to mimic MULTIACCUM readout if desired if (mimic_multiaccum) then begin ;; MULTIACCUM exposure time is the difference between the ;; first sample of the 0th group and the first sample of ;; the last group. This will in general be shorter than ;; the SUTR exposure time. exp_time = numsamrd * (numutr - 1) * frmtime ;; Update header fields sxdelpar, header, 'NUMUTR' sxdelpar, header, 'NUMSAMRD' sxdelpar, header, 'EXP_TIME' sxaddpar, header, 'NUMUTR', numutr sxaddpar, header, 'NUMSAMRD', numsamrd sxaddpar, header, 'EXP_TIME', exp_time endif ;; Calculate the time between groups of samples group_time = exp_time / (numutr - 1) ; Time required to read a group ;; Create and initialize an array to hold the time of the read ;; We only need one of these if we update it on the fly. time_of_read = dblarr(xsize, ysize) ;; To keep it simple, use the standard formula from Numerical ;; Recipes. ;; ;; slope = N * S[x_i * y_i] - S[x_i] * S[y_i] / ( N * S[x_i^2] - ;; S[x_i]^2 ) ;; First initialize all sums that we will need. To keep it ;; simple, I'll keep some redundant sums (e.g., the ;; information already is available to the program) N = 0 Sxy = dblarr(xsize, ysize) Sx = dblarr(xsize, ysize) Sy = dblarr(xsize, ysize) Sxx = dblarr(xsize, ysize) frame0 = fltarr(xsize, ysize); Need this to subtract residual bias during reference pixel correction. ;; Compute sums print, 'Fitting least squares' for i = 0, numutr - 1 do begin ;; Update the time of read for the current slice time_of_read[*,*] = i * group_time ;; Average current group of samples group_mean = dblarr(xsize, ysize); Initialize an array to ; hold running sum for j = 0, numsamrd -1 do begin ;; Print status information current_slice = i * numsamrd + j print, 'Reading in frame # '+strtrim(string(current_slice),2) ;; Get data for the current slice i1 = (xsize * ysize) * numsamrd * i + (xsize * ysize) * j i2 = i1 + xsize * ysize - 1 ;; Unfortunately, it seems like you have to open the ;; file every read. I tried leaving it open, but this ;; caused errors. fits_open, input_file, fcb fits_read, fcb, data, header, first=i1, last=i2, /data_only fits_close, fcb data = reform(data, xsize, ysize, /overwrite) ;; Do reference pixel correction if requested frame-by-frame if (not norefsub) then begin if (i eq 0) then begin ;; Update FITS header and get necessary information refsub_version = ' ' ; Create a string to hold refsub version refsub, version=refsub_version sxaddhist, 'Processed by refsub_'+refsub_version, header sxaddhist, 'Refsub using Spatial Averaging', header ;; Save the 0th frame away to allow for subtracting residual bias frame0 = data endif else begin data = data - frame0 ;Subtract residual bias refsub, data, sca=sca_type; Perform reference pixel correction data = data + frame0 ; Add residual bias back in endelse endif ;; Update group mean group_mean = group_mean + data endfor group_mean[*,*] = group_mean[*,*] / numsamrd ; Take the average ;; Compute sums for least square fit line N = N + 1 Sxy = Sxy + time_of_read[*,*] * group_mean[*,*] Sx = Sx + time_of_read[*,*] Sy = Sy + group_mean[*,*] Sxx = Sxx + time_of_read[*,*] * time_of_read[*,*] ;;;print, time_of_read[0,0], group_mean[0,0] endfor ;; Output file should be in ADU, not ADU/s signal = fltarr(xsize, ysize) signal = exp_time * (N * Sxy - Sx * Sy) / ( N * Sxx - Sx * Sx ) ;; Update FITS header for output file sxaddhist, 'Calibrated using mkmulti_'+version, header ;; Perform reference pixel correction, if requested, using ;; Spatial Averaging technique if (not norefsub) then begin refsub, signal, sca=sca_type; This over-writes signal endif ;; Create output filename if one is not supplied on command line if (not keyword_set(output_file)) then $ output_file = strmid(input_file, 0, strpos(input_file, '.fits', /REVERSE_SEARCH))+'.cal.fits' ;; Write the FITS file print, 'Writing Fits file' fits_write, output_file, signal, header end