Jeff's functional description is commented with notes of how to do it in practice (memory use, data structure details, etc.) and where to get the necessary information from. Comments begin with a ---> ---> Variable names are denoted by single quotes. ---> This will be a IRAF task similar to stis$x1d. It will be written in C and use a CL wrapper. Task name: sscatter ? ---> Reference file names necessary for the algorithm will be input via task parameters exclusively, since no provisions exist in current STIS headers to hold then. ---> As it stands, the algorithm cannot be used with multiple-IMSET files. The reason is that it uses calstis6 iteratively, and calstis6, once initiated, always works on an entire file, processing all IMSETS at once. So it's no use to wrap the algorithm in an external, IMSET-scanning loop. One possible solution is to write the C code such that only single IMSET files can be processed, and eventually provide the sscatter Cl wrapper with the ability to unpack IMSETS and table extensions in parallel and call the C code in a loop sequence. *** Prepare to enter main loop *** Run standard calstis extrinf = 'yes' Howk & Sembach algorithm disabled ---> This can be performed by a call to CalStis6 with appropriate parameters. Output goes to a temporary file. What to do with the lengthy stdout output ? Cannot be turned off. Read SCI, ERR, and DQ from _raw image ---> A simple GetSingleGroup call into 'input' SingleGroup variable. ---> Are ERR and DQ used at all ? Read WAVELENGTH, GROSS, NET, SPORDER, and EXTRLOCY from _x1d file ---> should we read DQ as well ? Would it be useful for the algorithm ? ---> Data goes into 'extracted', an array of data structures of type: typedef struct { short sporder; short npts; double scale; double *wave; double *gross; double *net; short *dq; double *extrlocy; } RowContents; Each array element stores data from a table row. To implement this step, we need a function that can *read* a _x1d table. We only have now functions that *write* _x1d tables. Calculate NM+1 background positions [NM = number of echelle orders] [IDT uses variable name POS for our EXTRLOCY] Midpoints between EXTRLOCY order positions Linear extrapolation of last 2 points in EXTRLOCY to get last background position on each end ---> Store in 2-D double array. If hires pixel mode, simulate some lores pixel quantities: Copy SCI, ERR, and DQ to SCI_HI, ERR_HI, and DQ_HI Bin 2x2 groups of pixels in SCI, ERR, and DQ Bin adjacent pixel pairs in WAVELENGTH, GROSS, NET, and EXTRLOCY ---> AllocSingleGroup with appropriate size in 'input_image' SingleGroup'. alloc array of RowContents data structures in 'extracted_data'. alloc 1-D arrays in 'extracted_data'. bin image from 'input' into 'input_image'. bin arrays from 'extracted' into 'extracted_data'. If not hires, just copy pointers so 'input_image' and 'extracted_data' become the working data areas. This will mix up DQ bits. Define internal algorithm parameters based on OPT_ELEM from header Note that SCALE may be a function of SPORDER Niter = 3 OPT_ELEM == E140M: SCALE = (14.5 - 0.275*SPORDER + 0.001435*SPORDER*SPORDER) Enforce minimum value of 2.5 in SCALE from formula above INC_SCALE = 1.4 YYAMP = 1.0 OPT_ELEM == E140H: SCALE = 1.5 INC_SCALE = 1.3 YYAMP = 1.25 OPT_ELEM == E230M: SCALE = 1.0 INC_SCALE = 1.0 YYAMP = 0.0 OPT_ELEM == E230H: SCALE = 1.5 INC_SCALE = 1.0 YYAMP = 0.0 [YYAMP above is called YYOFF by IDT, who redefine YYOFF to be a [vector that is the linear interpolation through points: [ (0,0), (250,YYAMP), (500,0), (750,-YYAMP), (1000,0) [ sampled onto x values from 0 to 999 in steps of 1 [In main loop, IDT code computes ECHELLE_SCAT_OFFSET from YYOFF, but [ it should be more efficient for us to compute ECHELLE_SCAT_OFFSET as [ needed from YYAMP.] ---> Get OPT_ELEM from primary header. Store 'inc_scale' and 'yyamp' as double scalars, and 'scale' in 'extracted_data' structure. Remove heliocentric velocity correction from WAVELENGTH, if applied ---> This is a tricky one. If heliocentric correction took place during the _x1d extraction, the only places where info was recorded are HISTORY keywords in the _x1d file. It's no use to check for the HELCORR='perform' switch in the file, because that can be overriden with x1d.helcorr="omit" by the user. This forces us to violate the FITS standard and look for values in an HISTORY keyword (I've seen this before...). get heliocentric correction from _x1d file header; correct wavelength arrays in 'extracted_data' data structure. Lookup true Xsize and Ysize of aperture in arcseconds. Example: for APERTURE == '0.1X0.2', Xsize = 0.096 and Ysize = 0.196 ---> Use modified version of GetApDes6 from calstis6 to store Xsize and Ysize in double scalar variables 'ap_xsize' and 'ap_ysize'. Build scattering kernels for given OPT_ELEM and CENWAVE FFT1, FFT2, FFT3, FFTO1, FFTO2, and FFTO3 are 2048x2048 complex arrays that must be constructed from detector halo and telescope PSF profiles in reference files (halo_*.fits and stis*c00.fits). The arrays are used once per iteration (currently 3) of the algorithm. Thus, we could trade the 200 Mb of (virtual) memory required to preload these arrays with CPU cycles to construct the arrays 3 times instead of once. IDT Variables returned by this processing step: MSCAT INT = Array[140] (described below) NELEM INT = Array[140] (described below) SCAT FLOAT = Array[947, 140] (described below) ESCAT_PSF FLOAT = Array[121] MRIP INT = Array[53] (described below) WRIP DOUBLE = Array[500, 53] (described below) RIP DOUBLE = Array[500, 53] (described below) NWAVE INT = 3 WAVE1 FLOAT = 1150.00 WAVE2 FLOAT = 1200.00 WAVE3 FLOAT = 1275.00 FFT1 COMPLEX = Array[2048, 2048] FFT2 COMPLEX = Array[2048, 2048] FFT3 COMPLEX = Array[2048, 2048] FFTO1 COMPLEX = Array[2048, 2048] FFTO2 COMPLEX = Array[2048, 2048] FFTO3 COMPLEX = Array[2048, 2048] Some of these variables still need to be described in detail, along with the reference files from which they are constructed. ---> echelle_scat_read in IDL code. ---> halo_*.fits: single 1024 sq image, empty header. For E140 there are three possible wavelengths, for E230 just one. Names in HAL1FILE, HAL2FILE, HAL3FILE keywords in primary header. Repeat for E230 ? This file is the first parameter in make_fft. ---> stis*c00.fit: single 196 sq image, empty header. OTA PSF ? ---> Scattering function in SCATFILE keyword. Find integer row and fractional order where scattering kernels are centered. Read reference wavelengths for up to 3 kernels from reference file ---> Which reference file ? In stis*c00.fit files, wavelength appears to be coded in a COMMENT header keyword, so we have to violate the FITS standard to access it. Better to add this info to input image header (KERNWAV1, KERNWAV2, KERNWAV3) Use linear interpolation of A2CENTER vs. WAVELENGTH(512,*) to find row numbers (LINEPOS1, LINEPOS2, LINEPOS3) at reference wavelengths. ---> This is performed in the x1d table (EXTRLOCY as a function of WAVELNGTH) Use linear interpolation of SPORDER vs. A2CENTER to find fractional order numbers (MPSFPOS1, MPSFPOS2, MPSFPOS3) at reference wavelengths. ---> This step will create scalar variables 'lineposX' and 'mpsfposX'. Calculate fractional order number (MLINE) for each row in image. ---> Alloc double array 'mline'. ---> This is the same computation as in previous step, except that it scans the entire set of image lines. Results stored in 'mline'. Extend observed WAVELENGTH to NSBIG pixels creating WBIG) NSBIG 3000 implies 988 pixels of padding on each end to handle scattering onto the image from spectrum off the edge of the detector. ---> Where does NSBIG comes from ? [This model ignores physical baffling in instrument?] Dimensions of extended wavelength array (WBIG) are (NSBIG,NM). I1 = (NSBIG-1024)/2 = 988 is the first pixel in WBIG with real data I2 = I1 + 1023 = 2011 is the last pixel in WBIG with real data Loop through SPORDER Assume extrapolated wavelengths have uniform dispersion (DISP) equal to the mean dispersion across the current order (I). Insert WAVELENGTH into WBIG(I1,I) through WBIG(I2,I). Insert extrapolation into ends of WBIG End loop ---> Alloc 2-D 'wbig' double array. ---> Compute 'i1' and 'i2'. ---> Loop through SPORDER (from 'extracted_data') ---> Compute mean dispersion from wavelength array. ---> Fill center of 'wbig[*][sporder]' row with data from 'extracted_data' ---> Fill edge regions with extrapolated wavelength values. ---> End loop. Read relevant echelle ripple data from ripple_* reference files. Each of the 4 gratings has its own reference file Files are FITS binary tables with columns CENWAVE, M, WAVELENGTH, and BLAZE Only load into memory reference file rows for which: Reference file CENWAVE equals observed CENWAVE Reference file M is in observed SPORDER list MRIP is list of orders satisfying two conditions immediately above RIP contains 500 ripple points for each order in MRIP (from BLAZE column) WRIP contains wavelengths for every point in RIP (from WAVELENGTH column) Note that many orders in SPORDER will not have a ripple function Extend ripple functions to 3000 pixels, padding each end with 988 pixels. Dimensions of extended ripple array (BLAZE) are (3000,NM). Loop through orders MRIP(GOOD) is closest match to current order number. Extract ripple wavelengths (WW) from WRIP for closest match. Multiply elements in WW by { MRIP(GOOD) / current order number M }. Use linear interpolation of RIP for closest match vs. WW to determine BLAZE at each wavelength in extended wavelength scale WBIG. For 988 exprapolated points in BLAZE on both ends of each order, linearly extrapolate the two end values in WRIP and RIP. [Check for negative values of BLAZE; set to zero? Not in IDT package] Endloop ---> Alloc 'blaze' 2-D array. ---> Pretty much described above. Extracting closest match can be done by a modified InterpTrace6 routine from calstis6. Compute median rows per Angstrom (DYDW) in cross-dispersion direction. Disp=(A2CENTER(1)-A2CENTER(0))/(WAVELENGTH(512,1)-WAVELENGTH(512,0)), etc. DWDY is median of NM-1 values in Disp Create model spectrum (FMERGE). Loop through orders Multiply observed NET by SCALE (from reference file above) to make NETM Divide NETM by pixels 988 through 2011 in BLAZE orders into one long spectrum, making WMERGE and FMERGE Enforce minimum value of -20 counts/pixel in FMERGE ---> idt code uses -20 counts/pixel/second instead [Minimum of -20 counts/pixel seems awfully low. Use -1 perhaps?] End loop -- Separate routine to merge echelle orders. May exist already? Define as separate function, since called in two places. Ignore first 4 and last 3 pixels of each order, which may be bad Do not combine fluxes where wavelength from adjacent orders overlap Switch to next order at wavelength midway between last wavelength of current order and first wavelength of next order Read echelle scattering profiles from ech_scat_* reference files. Each of the 4 gratings has its own reference file Files are FITS binary tables with columns M, NELEM, and SCAT Only retain rows for echelle orders contained in observation NELEM contains number of valid pixels in SCAT for each order SCAT contains echelle grating scatter profile for each order, plus any zero padding required to fill out the rectangular data area in file. NELEM should be odd, so center of kernel is well defined Only the first NELEM points in SCAT are valid for each order All orders should have echelle scattering profiles [We will only read orders in SPORDER, so don't need separate MSCAT list.] *1* Begin Main Loop *1* Allocate arrays for scattering model IM_MOD1 is echelle scattering component 1 (wings, real*4, 1024x1024) IM_MOD2 is echelle scattering component 2 (core, real*4, 1024x1024) O_MOD is model of image without scattering (real*4, 1024x1024) These array are initialized to zero at beginning of each iteration Create model spectrum (FBIG) on extended wavelength scale (WBIG). [IDT uses FLUX for both normal and extended wavelength scale. We will use] [ FLUX for observed flux on the observed wavelength scale and define FBIG] [ to be the flux extrapolated to the extended wavelength scale WBIG.] Loop through orders Linearly interpolate FMERGE vs. WMERGE onto WBIG for each order. In first order, where WBIG < WMERGE(0), set FBIG to FLUX(0) In last order, where WBIG > WMERGE(last), set FBIG to FLUX(last) ---> Multiply each order by blaze. End loop Build image containing scattered light prediction *2* Begin Loop Through Orders *2* ..SCALE_LSF is echelle scattering profile (SCAT) for current order ---> SCALE_LSF = scf.scfunc[].values Pixel offset (LSF_OFF) for profile is NELEM/2 (truncate remainder). Recall NELEM should be odd, so SCALE_LSF should have 2*LSF_OFF+1 points. ECHELLE_SCAT_OFFSET is the the linear interpolation through points: (0,0), (0.25*NELEM,YYAMP), (0.5*NELEM,0), (0.75*NELEM,-YYAMP), (NELEM,0) sampled onto x values from 0 to NELEM-1 in steps of 1 [IDT code uses frebin(YYOFF,NELEM) to compress/expand NELEM=1000 case] [It is simpler to compute ECHELLE_SCAT_OFFSET directly from definition] ..Define [two, not three] echelle scattering functions SCALE_LSF1 = SCALE_LSF with central 11 pixels (LSF_OFF-5 to LSF_OFF+5) set equal to SCALE_LSF(LSF_OFF-5), i.e. clip peak. SCALE_LSF2 = SCALE_LSF - SCALE_LSF1, i.e. wings set to zero [SCALE_OLSF = SCALE_LSF2, identical, so no need for both] ..Extract data for current order [we may just define pointers] Recall that I1 and I2 define the set of pixels in WBIG and associated arrays that corresponding to to actual pixels in observed image. E = BLAZE array elements for current order ---> renamed eblaze to ease finding with editor. EONIMAGE = E(I1:I2) is part of blaze function actually on image FORDER = FBIG array elements for current order [IDT variable FONIMAGE is defined, but not used, so ignore] YORDER_ONIMAGE = EXTRLOCY array elements for current order WORDER_ONIMAGE = WAVE array elements for current order WORDER = WBIG(I1:I2) is entire extended wavelengths for current order [The arrays in this paragraph facilitate vector operations in IDL and and perhaps add clarity, but the arrays are never modified. Thus, we can use pointers into the original arrays instead. I will use the names in the description below to facilitate comparison with the IDL code.] ..Indexes (ISTART, ISTOP) for extended arrays (e.g. WBIG) where scattering by LSF can affect actual image region (I1 to I2). Recall LSF_OFF is half the scattering kernel width ISTART = I1 - LSF_OFF is lowest pixel (e.g., in WBIG) to consider ISTOP = I2 + LSF_OFF is highest pixel (e.g., in WBIG) to consider *3* Begin Loop Through Pixels ISTART to ISTOP in Extended Arrays *3* ....Loop info: These are all the pixels that could scatter light into image IPIX is the index of the current pixel [IDT uses variable "I"] ....Define index ranges used for addressing arrays while processing pixel IMAGE1POS and IMAGE2POS bound pixels in observed image that can receive flux from current pixel by LSF scattering ILSF1 and ILSF2 bound points in *LSF* arrays to use in calculation ILSF1 = 0 ILSF2 = NELEM - 1, recall indexes begin at zero IMAGE1POS = IPIX - LSF_OFF - I1 If (IMAGE1POS<0) { ILSF1=ILSF1-IMAGE1POS; IMAGE1POS=0; } IMAGE2POS = IPIX + LSF_OFF - I1 If (IMAGE2POS>1023) { ILSF2=1023+ILSF2-IMAGE2POS; IMAGE1POS=1023; } ^ V ---> If (IMAGE2POS>1023) { ILSF2=1023+ILSF2-IMAGE2POS; IMAGE2POS=1023; } Note: (IMAGE2POS-IMAGE1POS) equals (ILSF2-ILSF1) Examples: IPIX IMAGE1POS IMAGE2POS ILSF1 ILSF2 ===================================================================== I1-LSF_OFF 0 0 2*LSF_OFF =NELEM-1 NELEM-1 I1-LSF_OFF+1 0 1 2*LSF_OFF-1=NELEM-2 NELEM-1 I1-LSF_OFF+2 0 2 2*LSF_OFF-2=NELEM-3 NELEM-1 ------------ ---------------- ----------- ------------------- ------- I1-1 0 LSF_OFF-1 LSF_OFF-1 NELEM-1 I1 0 LSF_OFF LSF_OFF NELEM-1 I1+1 0 LSF_OFF+1 LSF_OFF+1 NELEM-1 ------------ ---------------- ----------- ------------------- ------- I1+LSF_OFF-1 0 2*LSF_OFF-1 1 NELEM-1 I1+LSF_OFF 0 2*LSF_OFF 0 NELEM-1 I1+LSF_OFF+1 1 2*LSF_OFF+1 0 NELEM-1 ===================================================================== I2-LSF_OFF-1 1023-2*LSF_OFF-1 1022 0 NELEM-1 I2-LSF_OFF 1023-2*LSF_OFF 1023 0 2*LSF_OFF =NELEM-1 I2-LSF_OFF+1 1023-2*LSF_OFF+1 1023 0 2*LSF_OFF-1=NELEM-2 ------------ ---------------- ----------- ------- ------------------- I2-1=I1+1022 1023-LSF_OFF-1 1023 0 LSF_OFF+1 I2 =I1+1023 1023-LSF_OFF 1023 0 LSF_OFF I2+1=I1+1024 1023-LSF_OFF+1 1023 0 LSF_OFF-1 ------------ ---------------- ----------- ------- ------------------- I2+LSF_OFF-2 1021 1023 0 2 I2+LSF_OFF-1 1022 1023 0 1 I2+LSF_OFF 1023 1023 0 0 ===================================================================== If IMAGE2POS < IMAGE1POS, skip rest of loop for this pixel. [This should never be true, so trap as an error if debugging] *4* [Begin loop through pixels to receive scattered light?] *4* ....[..]Loop info: IDT code avoids explictly looping by using IDL array operations ....[..]Translate FBIG in IPIX to flux elsewhere in echelle blaze function F is array defined for pixels "j" in range IMAGE1POS to IMAGE2POS F(j) = FORDER(IPIX) * EONIMAGE(j)/E(IPIX) Right term in product replaces blaze at IPIX with blaze at "j" [Conceptually clearer to multiply *LSF* by EONIMAGE(j)/E(IPIX) below?] Assumes *LSF* scattering happens before or at grating [Is this true?] ....[..]Determine where to place scattered light in the model images. XX is the list of integer columns that will recieve scattered light [IDT defines XX in advance to use vector operations in IDL. We can [ simply replace references to XX(j) with IMAGE1POS+j] Y is the list of fractional rows that follow scattered light vs. XX Scattering assumed to follow lines of constant wavelength on image, whereas positions in YORDER_ONIMAGE reflect changes in wavelength. ECHELLE_SCAT_OFFSET is presumably an empirical fudge to match data. dw=WORDER_ONIMAGE(IMAGE1POS+j)-WORDER(IPIX) Y(j)=YORDER_ONIMAGE(IMAGE1POS+j)-DYDW*dw+ECHELLE_SCAT_OFFSET(ILSF1+j) The term DYDW*dw is the change in row number required to get from wavelength WORDER in column IMAGE1POS+j to wavelength WORDER(IPIX). ....[..]Calculate fraction of scattered light for row above and below Y Y1 is the integer part of Y FR2 is the fractional part of Y [Y2 is just Y1+1, which we can calculate as needed.] [FR1 is just 1-FR2, which we can calculate as needed.] ....[..]Add scattered light contribution into model images Add into each pixel along Y vs. XX the observed flux in pixel IPIX, modified for blaze function changes, line spread function, and fractional splitting between two bracketting rows. Loop (j) through number of columns to receive scattered light Maximum value of j is IMAGE1POS-IMAGE2POS ixx = IMAGE1POS + j iy1 = Y1(j) iy2 = Y1(j) + 1 ilsf = ILSF1 + j fr1 = 1 - FR2(j) IM_MOD1(ixx,iy1) = IM_MOD1(ixx,iy1) + F(j)*SCALE_LSF1(ilsf)*fr1 IM_MOD1(ixx,iy2) = IM_MOD1(ixx,iy1) + F(j)*SCALE_LSF1(ilsf)*FR2(j) IM_MOD2(ixx,iy1) = IM_MOD2(ixx,iy1) + F(j)*SCALE_LSF2(ilsf)*fr1 IM_MOD2(ixx,iy2) = IM_MOD2(ixx,iy1) + F(j)*SCALE_LSF2(ilsf)*FR2(j) If final iteration of main loop, then build object image O_MOD(ixx,iy1) = O_MOD(ixx,iy1) + F(j)*SCALE_LSF2(ilsf)*fr1 O_MOD(ixx,iy2) = O_MOD(ixx,iy2) + F(j)*SCALE_LSF2(ilsf)*FR2(j) [Recall that SCALE_LSF2 equals the IDT SCALE_OLSF] Endif End loop through columns to receive scattered light This loop will probably be expanded to include previous 4 paragraphs *4* [End loop through pixels to receive scattered light?] *4* *3* End loop through pixels in extended arrays *3* *2* End loop through orders *2* Contents of model images at this point: IM_MOD2 is current guess for what image would look like with no scattered light and infinite resolution. Orders angle across image at positions in EXTRLOCY. All flux in order is contained in two pixel wide swath. Spectral lines are visible in image. Cores of strongly saturated line cores are negative in at least the first iteration. IM_MOD1 is current guess for location of scattered light. Maximum light is in several pixel wide swath *midway* between orders [why?], but there is a pervasive component that sits under orders as well. Strong absorption lines appear as diagonal dark swaths the points where such features would appear in adjacent orders. Convolve scattering component (IM_MOD1) with cross-disperser profile ESCAT_PSF is the cross-disperser scattering profile ESCAT_PSF is read from reference file before entering main loop Loop through 1024 columns Extract current column from IM_MOD1 Pad both ends with 50 extra pixels of zeros. This is an approximation, since there will actually be orders off the the edge, but a fair bit of code would be required to extrapolate the visible orders off the top and bottom of the image. 50 pixels is not always enough to mathematically separate top and bottom of image during convolution, since ESCAT_PSF can have more than 101 pixels. [Boost amount of padding to the number of elements in ESCAT_PSF divided by two (with integer truncation)?] Convolve padded column with ESCAT_PSF Replace original column with corresponding part of convolved column End loop After column by column convolution, contrast is reduced between pervasive component and peaks midway between orders, but peaks are still present. Sum two model image components to construct composite model image (IM_MOD) Loop through all pixels in model image IM_MOD = IM_MOD1 + IM_MOD2 End loop Note that convolution with cross-disperser scattering function can only be done once the scattered light model image (IM_MOD1) has been fully constructed. One could bypass the need for multiple 1024x1024 images at this stage by constructing IM_MOD1 and convolving with cross-disperser scattering function. Then make a second pass through orders, pixels in extended arrays, and pixels to receive scattered light, adding in the IM_MOD2 component. The savings in memory probably isn't worth the repeated code, especially since array is reused below. Convolve model image with first telescope PSF and detector halo profile Allocate 2048x2048 array (BIG) and set contents to zero Insert IM_MOD into low quadrant of BIG, i.e BIG(0:1023,0:1023) [Check that this placement is consistent with non-IDL FFT algorithm] Calculate discrete Fourier transform of BIG Multiply element by element with Fourier transform of PSF and Halo (FFT1) FFT1 is computed from reference files before entering main loop Compute real part of inverse transform of product, overwriting BIG Extract 1024x1024 image from low quadrant of BIG Multiply result by 2048*2048 [depending on normalization of FFT routine] Save result in IM_MOD1, reusing array from above with different content Convolve model image with second telescope PSF and detector halo profile Multiple kernels are used to handle wavelength variations in PSF Skip this step if STIS mode described by a single PSF and halo profile Same procedure as first component, except: BIG need not be reallocated, but set contents to zero Convolve with FFT2, rather than FFT1 Store result in IM_MOD2, instead of IM_MOD1 Convolve model image with third telescope PSF and detector halo profile Multiple kernels are used to handle wavelength variations in PSF Skip this step if STIS mode described by only two PSF and halo profiles Same procedure as first component, except: BIG need not be reallocated, but set contents to zero Convolve with FFT3, rather than FFT1 Store result in IM_MOD3, instead of IM_MOD1 If final iteration, convolve object model (O_MOD) with PSF and halo Same procedure as first, second and third components above, except: Replace FFT1,FFT2,FFT3 with FFTO1,FFTO2,FFTO3 Replace IM_MOD1,IM_MOD2,IM_MOD3 with O_MOD1,O_MOD2,O_MOD3 Coadd convolved model images with row weights based on order number. This scheme is designed to handle variations in PSF shape with row Up to 3 model images are constructed using constant PSF kernels Each kernel applies to a reference wavelength (WAVE1,WAVE2,WAVE3) Recall from above main loop: LINEPOS1,LINEPOS2,LINEPOS3 are rows corresponding to WAVE1,WAVE2,WAVE3 MPSFPOS1,MPSFPOS2,MPSFPOS3 are orders coresponding to WAVE1,WAVE2,WAVE3 MLINE is fractional order associated with each image row If NWAVE is 1 (constant PSF) then Copy all of IM_MOD1 into IM_MOD Else Only copy rows 0 through LINEPOS1-1 of IM_MOD1 into IM_MOD [IDT algorithm copies entire image in both cases, but unnecessary] Endif If NWAVE is at least 2 then Fill rows LINEPOS1 through LINEPOS2-1 with FRAC1*IM_MOD1 + FRAC2*IM_MOD2 FRAC1 is unity at LINEPOS1 and zero at LINEPOS2 FRAC2 is zero at LINEPOS1 and unity at LINEPOS2 The ramp is not quite linear with row number, however. Instead, weights are linear with order number in MLINE Thus, FRAC2(row) = (MLINE(row) - MPSFPOS1) / (MPSFPOS2 - MPSFPOS1) FRAC1(row) = 1 - FRAC2(row) If NWAVE is 2, fill rows LINEPOS2 through 1023 of IM_MOD with IM_MOD2 Endif If NWAVE is at least 3 then Fill rows LINEPOS2 through LINEPOS3-1 with FRAC2*IM_MOD2 + FRAC3*IM_MOD3 FRAC2 is unity at LINEPOS2 and zero at LINEPOS3 FRAC3 is zero at LINEPOS2 and unity at LINEPOS3 The ramp is not quite linear with row number, however. Instead, weights are linear with order number in MLINE Thus, FRAC3(row) = (MLINE(row) - MPSFPOS2) / (MPSFPOS3 - MPSFPOS2) FRAC2(row) = 1 - FRAC3(row) If NWAVE is 3, fill rows LINEPOS3 through 1023 of IM_MOD with IM_MOD3 [IDT code reuses the name FRAC1 and FRAC2, but clearer to change name] Endif If final iteration, coadd object models (O_MOD1,O_MOD2,O_MOD3) with weights Same as procedure for model images above, except: Replace IM_MOD with O_MOD Replace IM_MOD1,IM_MOD2,IM_MOD3 with O_MOD1,O_MOD2,O_MOD3 Add ghosts due to reflections from window above NUV MAMA Define transformation (rotation and expansion) coefficients for ghost image 2x2 arrays for E230M Mode: KX = / -33.9909 -0.000844246 \ = / KX(0,0) KX(1,0) \ \ 0.998378 3.51736e-06 / \ KX(0,1) KX(1,1) / KY = / 5.46450 1.00869 \ \ 9.48275e-05 -1.10486e-06 / 2x2 arrays for E230H Mode: KX = / -16.4834 -0.000280014 \ \ 0.998752 1.51807e-06 / KY = / 5.59971 1.00791 \ \ 1.58386e-05 -7.14570e-07 / Build transformed image using nearest neighbor polynomial warping of image Call new image GHOST [no name used in IDT code] Calculate which pixel (ii,jj) in IM_MOD maps to pixel (i,j) in GHOST x = KX(0,0) + j*KX(1,0) + i*KX(0,1) + j*i*KX(1,1) y = KY(0,0) + j*KY(1,0) + i*KY(0,1) + j*i*KY(1,1) Use nearest neighbor approximation: ii = round(x), jj = round(y) Use closest edge pixel for points that map off edge of IM_MOD GHOST(i,j) = IM_MOD(ii,jj) Scale ghost to 0.3% of original image intensity (multiply by 0.003) Convolve with 3x3 pixel smoothing kernel: 1 / 0.2 0.5 0.2 \ Normalized KERNEL = --- * | 0.5 1.0 0.5 | 3.8 \ 0.2 0.5 0.2 / Convolve GHOST with KERNEL Add GHOST to IM_MOD, pixel by pixel General notes about transformation scheme: KX is bilinear mapping of GHOST (i,j) indexes back to IM_MOD ii index KY is bilinear mapping of GHOST (i,j) indexes back to IM_MOD jj index Coefficients assume that all indexes begin with zero! KX = / 0 0 \ \ 1 0 / is an identity mapping for columns KY = / 0 1 \ \ 0 0 / is an identity mapping for rows Extract new FLUX spectrum from IM_MOD using standard calstis Howk & Sembach algorithm disabled Produce MOD_NET, which will be compared to NET from previous extraction ---> NET = wx1d.net, MOD_NET = wx1d2.net Apply correction to extracted flux based on change in NET NET is the net counts from previous iteration ---> No, NET is the original input data. Flux increment (INC) is INC_SCALE * (NET - MOD_NET) Recall INC_SCALE comes from an OPT_ELEM specific reference file. Update FLUX: FLUX(j) = FBIG(I1+j) + INC(j) INC_SCALE boosts the flux correction for E140M and E140H, perhaps to accelerate convergence [while hopefully avoiding overshoot?] Update merged model spectrum (FMERGE). Loop through orders Divide FLUX by pixels 988 through 2011 in BLAZE (i.e. BLAZE_ONIMAGE) orders into one long spectrum, making WMERGE and FMERGE ---> The following step is *not* performed in the idl code ! Enforce minimum value of -20 counts/pixel in FMERGE [Minimum of -20 counts/pixel seems awfully low. Use -1 perhaps?] End loop *1* End Main Loop *1* ---> Free 2-D array 'blaze'. ---> Free 2-D array 'wbig'. ---> Free 1-D array 'mline'. Construct scattered light image SCAT = IM_MOD - O_MOD Return to hires pixel mode, if necessary Reapply heliocentric velocity correction to WAVELENGTH, if removed above ---> Why ? Subtract scattered light image from original raw image IMAGE = IMAGE - SCAT [No need for IDT variable NEW_IMAGE, since IMAGE array can be used] Extract final spectrum from IMAGE using standard calstis Howk & Sembach algorithm disabled Write results in standard calstis format. ---> FreeSingleGroup 'input_image'. free 1-D arrays in 'extracted data' free 'extracted_data' data structure. _________________________________________________________________________ About 15 min on a Ultra 60 with 380 Mb memory (aura). *Way* more if virtual memory is necessary. Peak data area uses about 200 Mb. FIX THE FILE NAME (.FITS) PROBLEM !! Psets: Since x1d has now a large number of parameters, I reorganized its user interface. Parameters specific to IDT algorithm: Keywords: The IDT algorithm depends on a number of reference files whose paths/names must be supplied via primary header keywords in the input image. There is also a calibration switch that enables the pipeline version of the code to select the IDT algorithm. This switch is ignored when the task is executed via the 'x1d' script. An example taken from my test data set (o4qx04010_flt.fits) is as follows: IDTCORR = 'perform ' SCATTAB = 'reffiles/ech_scat_e140h.fits' SPSFFILE= 'reffiles/ech_xdisp_e140h.fits' XDSPFILE= 'reffiles/xdisp_e140h.fits' RIPPLTAB= 'reffiles/ripple_e140h.fits' KERNWAV1= 1150 KERNWAV2= 1200 KERNWAV3= 1275 HAL1FILE= 'reffiles/halo_e140h_1150.fits' HAL2FILE= 'reffiles/halo_e140h_1200.fits' HAL3FILE= 'reffiles/halo_e140h_1275.fits' PSF1FILE= 'reffiles/stis1150c00.fit' PSF2FILE= 'reffiles/stis1200c00.fit' PSF3FILE= 'reffiles/stis1300c00.fit' Of course, I made up the keyword names; feel free to suggest other names if these are not OK. Temporary files: Because the standard calstis6 code was designed around disk files, the IDT algorithm is constrained to use disk files, both images and tables, to pass information into and out of calstis6. These files are deleted at the end of a successful completion, but in case of premature abort they are left in place (in the current working directory). Their names begin with the prefix "TEMP" and they must be manually removed whenever the task is interrupted. Verbosity: The 1-D extraction step is performed four times. This step either prints a very summarized report with just the spectral order numbers, or a full length report (the original _trl report from calstis6), depending on the 'verbose' task parameter setting. Once selected, the format is the same in all four iterations, but we can easily change this beahvior. For example, we can have the summarized report be printed at the first three iterations and the long one only at the end. What should we do ? What to do when called from the pipeline ?