PRO msu_xrt_prep, input1, input2, index_out, image_out, $ image_out2, $ normalize=normalize, $ float=float, $ clean_type=clean_type, $ bsub_type=bsub_type, $ dark_type=dark_type, $ nsigma=nsigma, $ nmed=nmed, $ stop_bleed=stop_bleed, $ despike_despot=despike_despot, $ only_despike=only_despike, $ only_despot=only_despot, $ sens_despike=sens_despike, $ spotthresh=spotthresh, $ miss_map=miss_map, $ n_miss_pixels=n_miss_pixels, $ grade_type=grade_type, $ grade_map=grade_map, $ n_grade_pixels=n_grade_pixels, $ coalign=coalign, $ spike_map=spike_map, $ n_spike_pixels=n_spike_pixels, $ uncert_map=uncert_map, $ nosaturation=nosaturation, $ ; MSU local darkff=darkff, $ ; MSU local run_time=run_time, $ quiet=quiet, $ verbose=verbose, $ qabort=qabort, $ qstop=qstop ; ========================================================================= ; ;+ ; PROJECT: ; Solar-B / XRT ; ; NAME: ; XRT_PREP ; ; CATEGORY: ; Data calibration ; ; PURPOSE: ; Process XRT image(s) from Level-0 (raw) to Level-1 (calibrated). ; ; ; CALLING SEQUENCE: ; XRT_PREP, input1, input2, index_out, image_out, ; [image_out2] ; [,/normalize] [,/float] ; [,clean_type=clean_type] [,bsub_type=bsub_type] ; [,dark_type=dark_type] [,nsigma=nsigma] [,nmed=nmed] ; [,/stop_bleed] [,despike_despot=despike_despot] ; [,only_despike=only_despike] ; [,/only_despot] [,sens_despike=sens_despike] ; [,spotthresh=spotthresh] ; [,miss_map=miss_map] [,n_miss_pixels=n_miss_pixels] ; [,grade_type=grade_type] [,grade_map=grade_map] ; [,n_grade_pixels=n_grade_pixels] [,coalign=coalign] ; [,spike_map=spike_map] [,n_spike_pixels=n_spike_pixels] ; [,uncert_map=uncert_map] [,run_time=run_time] ; [,/quiet] [,/verbose] [,qabort=qabort] [,/qstop] ; ; ; INPUTS: ; *** There are 2 methods for calling XRT_PREP: *** ; ; | INPUT1 - [Mandatory] (string scalar or array [Nimg]) ; Case 1.| The input XRT FITS filelist name(s). ; | INPUT2 - [Mandatory] (long array, [Nimg]) ; The dataset number(s) to extract and process. ; ; | INPUT1 - [Mandatory] (structure array, [Nimg]) ; Case 2.| The index structure for each of the input images ; | (e.g., index from output of "read_xrt.pro".) ; | INPUT2 - [Mandatory] (integer array, [Nx,Ny,Nimg]) ; The input data array (cube). ; ; OUTPUTS: ; INDEX_OUT - [Mandatory] (structure array, [Nimg]) ; The updated index structure of the input images. ; DATA_OUT - [Mandatory] (integer array, [Nx,Ny,Nimg]) ; Processed output XRT images (data cube). ; Default data type is I*2. (See /FLOAT keyword.) ; ; OPTIONAL OUTPUT: ; DATA_OUT2 - [Optional] (integer array, [Nx,Ny,Nimg]) ; Processed output XRT images (data cube) like DATA_OUT ; but *WITHOUT* cosmetic correction for contamination ; spots, dust, and dust growth. ; Default data type is I*2. (See /FLOAT keyword.) ; ; OPTIONAL INPUT KEYWORD PARAMETERS: ; /NORMALIZE - (Boolean) Set to normalize output image to ; DN per sec. (Warning: data output will *not* ; default to floating point. Use the /FLOAT keyword. ; See Note #8.) ; /FLOAT - [Recommended] (Boolean) Set if you want to ; return floating point. [Default is I*2] ; (Users should *always* use this keyword for better ; precision, unless they have memory constraints. ; See Note #8.) ; CLEAN_TYPE - [Optional] Type of Fourier cleaning to use ; (see xrt_fourier_vacuum.pro for details): ; 0 = prefilter on semi-fixed streaks, then remove ; remaining Fourier streaks, stars, & smudges ; [default] ; 1 = remove semi-fixed streaks only ; 2 = remove Fourier streaks, stars, & smudges ; 3 = don't Fourier clean! ; BSUB_TYPE - [Optional] Type of (model) background subtraction to use: ; 0 = remove Nyquist ringing and large-scale background ; "ramp" [default] ; 1 = remove Nyquist ringing only ; 2 = remove large-scale background "ramp" only ; 3 = no (model) background subraction! ; (note: dark_type=1 resets bsub_type=3) ; DARK_TYPE - [Optional] Type of dark subtraction to use: ; 0 = Use model dark with median adjusted to median of ; cleaned darks nearby in time (if possible) [default] ; 1 = Use median of cleaned darks nearby in time ; (if possible) ; 2 = Use model darks without median adjustment ; (note: dark_type=1 resets bsub_type=3) ; NSIGMA - [Optional] Number of standard deviations ; beyond which an FFT pixel is ; considered bad in "Fourier star" removal. ; (see xrt_fourier_vacuum.pro for details): ; Default is 4.5 ; (Note: much below 4 not recommended) ; NMED - [Optional] Number of standard deviations ; above smoothed central minimum background ; to begin shielding (presumed) data from correction ; (see xrt_fourier_vacuum.pro for details): ; Default is 3.5 ; (Note: outside range of 2.0 to 4.5 not recommended) ; /STOP_BLEED - [Optional] (Boolean) Set = 1 if pixels corrected for ; saturation "bloom/bleed" at the edge of saturated ; region(s) are to be retained (otherwise they are ; reverted to their input values) ; /DESPIKE_DESPOT - [Optional] (Boolean or integer) Set = N (1<=N<~3) for N ; passes of removal radiation-belt/cosmic-ray spikes (two ; passes works well even for heavy spiking). This method ; uses convolution and thresholding to remove spikes and ; may remove small real features. ALSO automatically ; makes a cosmetic correction to contamination spots ; (spline-based; see for an alternative, ; median-cap cosmetic correction.) ; Default = no despiking or despotting if the keyword is ; not set. Combines /only_despike and /only_despot ; functions; overrides /only_despike and /only_despot ; if set. See also Note #4. ; /ONLY_DESPIKE - [Optional] (Boolean or integer) Set = N (1<=N<~3) for N ; passes of removal radiation-belt/cosmic-ray spikes (two ; passes works well even for heavy spiking). This method ; uses convolution and thresholding to remove spikes and ; may remove small real features. (Two passes works ; well even for heavy spiking.) This method uses ; convolution and thresholding to remove spikes and ; may remove small real features. DOES NOT ; make a cosmetic correction to contamination spots. ; Default = no despiking if the keyword is ; not set. ; See also Note #4. ; /ONLY_DESPOT - [Optional] (Boolean) Set to automatically ; makes a cosmetic correction to contamination spots ; (spline-based; see for an alternative, ; median-cap cosmetic correction.) ; Default = no despiking or despotting if the keyword is ; not set. DOES NOT remove radiation-belt/cosmic-ray ; spikes. Default = no despotting if the keyword is ; not set. ; SENS_DESPIKE - [Optional] (Float scalar) This number controls ; the aggressiveness of the despiking routine. ; If ("pixel value" GE SENS_DESPIKE*"mean of neighbors") ; then the pixel is treated as a spike. Values in the ; range of 1.1 to 1.3 work well. Values LT 1.0 are ; rejected for the default. Default = 1.1. ; (See for more information.) ; SPOTTHRESH - [Optional] (Float scalar 0 <= spotthresh <= 1, or -1) ; If between 0 and 1, the threshold fraction of a binned ; pixel which is spotted at which ; the binned pixel is treated as spotted. [default = 0.5] ; (Special case) if = -1, force treatment of *all* spots. ; (Normally the program only corrects spots if it ; determines they are sufficiently "dark" for speed.) ; GRADE_TYPE - [Optional] integer. Type of grade_map array to output ; (see grade_map below, and ) ; 0 = [default]: bytearr containing coded grade of pixel ; indicating pixel affected by saturation, saturation ; bloom/bleed, contamination spot, dust, has a "hot" ; response, and/or shows dust growth. ; 1 = double array, same as grade_type=0 but with added ; encoding for the number of 1x1 pixels affected within ; each binned pixel (useful only for binned data). ; 2 = do not return a grade array ; COALIGN - [Optional] Type of co-alignment to apply to XRT data ; 0 = [default]: use UFSS data ; 1 = use UFSS data + cross-corelation with AIA 335 ; 2 = no adjustment ; /QUIET - (Boolean) Set for no messages or errors. Default is ; some messages and all errors. ; /VERBOSE - (Boolean) Set for all messages, errors, and ; intermediate data listings. Suppresses quiet. ; /QSTOP - (Boolean) For debugging. ; ; OPTIONAL OUTPUT KEYWORD PARAMETERS: ; MISS_MAP - (byte array, [Nx,Ny,Nimg]) This is a 2D Boolean ; map of each image: ; 0: Image pixel had data. ; 1: Image pixel was missing data and was replaced ; with the image average. ; N_MISS_PIXELS - (long array, [Nimg]) Number of missing pixels ; found in each image. ; GRADE_MAP - (byte array, or double array if grade_type=1 ; [Nx,Ny,Nimg]) This is a 2D map of each image, ; with a value given by the sum of: ; 0: Image pixel was OK. ; 1: Image pixel was saturated and was replaced ; with a constant. ; 2: Pixel was affected by saturation bloom/bleed. ; (If stop_bleed=1 it was replaced by a local median) ; 4: Pixel was affected by a contamination spot. ; 8: Pixel was affected by a dust speck. ; 16: Pixel was a hot pixel. ; 32: Pixel was a potential dust growth pixel. ; If grade_type=1, a fraction is added encoding the ; number of 1x1 spot, dust, and hot pixels per binned ; pixel. Note that fractions are relevant only for spot, ; dust and hot grades. For each of these grade types, two ; digits to the right of the decimal place are reserved ; for the number of 1x1 pixels within the binned pixel ; which are of the given type. The digits are assigned in ; the same order as listed above, eg, for a 4x4 binned ; pixel which had 1x1 pixels of which 8 were OK, 3 were ; spot, 2 were dust, and 4 were hot, the grade would be ; 2+8+16 = 26, plus 0.030204, or 26.030204 ; (This is potentially useful for thresholding pixels ; by the degree of "damage" - eg, if only one 1x1 pixel ; in an 8x8 pixel is affected by a contamination spot, ; (ie, only 1/64th spotted) the user may wish to treat it ; as "OK" for purposes of their further analysis (see ; notes below and for more ; details). ; The grade map may be decoded into separate arrays ; (one map per grade-type) using . ; N_GRADE_PIXELS - (long array, [6,Nimg]) Number of graded pixels ; of each type (saturation, bloom/bleed, spot, dust, and ; hot respectively) found in each image. ; SPIKE_MAP - (byte array, [Nx,Ny,Nimg]) This is a 2D Boolean ; map of each image: ; 0: Image pixel was OK ; 1: Image pixel was found to be a particle hit and . ; replaced by local average ; N_SPIKE_PIXELS - (long array, [Nimg]) Number of spike pixels ; found in each image. ; RUN_TIME - (float scalar) The run time in seconds for ; . ; QABORT - (Boolean) Indicates that the program exited ; gracefully without completing. (Might be useful ; for calling programs.) ; 0: Program ran to completion. ; 1: Program aborted before completion. ; UNCERT_MAP - Returns an array of the photometric uncertainties ; within the XRT. These are NON-STATISTICAL errors ; caused by JPEG compression, and uncertainties in ; vignetting, dark correction and filtering. NO ; statistical/photometric errors are included. ; See the routine xrt_cvfact.pro to add in the photon statistics ; ; ; ; EXAMPLES (see note 9): ; ; Read in and prep a subset of XRT FITS files: ; IDL> xrt_prep, infiles, dset_arr, index_out, image_out, /float ; ; Prep an index/data pair: ; IDL> xrt_prep, index, data, index_out, image_out, /float ; ; Same as above, plus normalize all data to DN/sec: ; IDL> xrt_prep, index, data, index_out, image_out, /norm, /float ; ; Same as above, except allow output as I*2 to save memory: ; (See Note #8.) ; IDL> xrt_prep, index, data, index_out, image_out, /norm ; ; Same as above, return grade map ; IDL> grade = 1 ; note - set to non-zero value! ; IDL> xrt_prep, index, data, index_out, image_out, /norm, /float, $ ; IDL> grade_map=grade ; ; Same as above, return grade map with fractional contributions ; of spot, dust, & hot pixels per binned pixel ; IDL> grade = 1 ; note - set to non-zero value! ; IDL> xrt_prep, index, data, index_out, image_out, /norm, /float, $ ; IDL> grade_type=1, grade_map=grade ; ; Same as above, use a stronger Fourier filter ; IDL> grade = 1 ; note - set to non-zero value! ; IDL> xrt_prep, index, data, index_out, image_out, /norm, /float, $ ; IDL> grade_type=1, grade_map=grade, nsigma=4.0, nmed=2.5 ; ; Same as above, use a stronger Fourier filter, and pure model dark ; subtraction (somewhat faster) ; IDL> grade = 1 ; note - set to non-zero value! ; IDL> xrt_prep, index, data, index_out, image_out, /norm, /float, $ ; grade_type=1, grade_map=grade, nsigma=4.0, nmed=2.5, dark_type=2 ; ; Return the uncertainties as well as normalize floating point data ; IDL> uncmap = 1 ; note - set to non-zero value! ; IDL> xrt_prep, index, data, index_out, image_out, /norm, /float, $ ; uncert_map=uncmap ; ; Same as above, also do cosmetic correction for spots, despike 2 times, ; and return grade map and spike map ; IDL> uncmap = 1 ; note - set to non-zero value! ; IDL> spike = 1 ; note - set to non-zero value! ; IDL> grade = 1 ; note - set to non-zero value! ; IDL> xrt_prep, index, data, index_out, image_out, /norm, /float, $ ; uncert_map=uncmap, grade_map=grade, despike_despot=2, $ ; spike_map=spike ; ; Same as above, also do cosmetic correction for spots, no despike, ; and return grade map ; IDL> uncmap = 1 ; note - set to non-zero value! ; IDL> grade = 1 ; note - set to non-zero value! ; IDL> xrt_prep, index, data, index_out, image_out, /norm, /float, $ ; uncert_map=uncmap, grade_map=grade, /only_despot ; ; Same as above, also do cosmetic correction for ALL spots, no despike, ; and return grade map ; IDL> uncmap = 1 ; note - set to non-zero value! ; IDL> grade = 1 ; note - set to non-zero value! ; IDL> xrt_prep, index, data, index_out, image_out, /norm, /float, $ ; uncert_map=uncmap, grade_map=grade, /only_despot, spotthresh=-1 ; ; ; ; ; COMMON BLOCKS: ; None. ; ; ; PROCEDURE: ; ; A processing routine for preparing the XRT image data for analysis. ; The steps performed are: ; 1. Read in raw FITS image(s) from a filelist ; or read in a datacube and structure. ; 2. Fill pixels of value = 0 (missing data) with a "missing data ; value" = -999 ; 3. Replace near-saturated pixels for values > some threshold. ; 4. Option to remove radiation-belt/cosmic-ray hits and streaks. ; 5. Calibrate for read-out signals. ; 6. Remove the CCD bias (pedestal), and dark current ; (This is done in the routine which also ; calibrates for read-out signals.) ; 7. Remove vignetting ; 8. Option to normalize each image for exposure. ; 9. Option to compute map of calibration uncertainties. ; 10. Option to cosmetically correct for contamination spots and dust. ; 11. Output the corrected image(s) in an updated structure ; and data cube. ; ; In the future this routine may also: ; 12. Option to extract a subimage from each image. ; 13. Calibrate for flat-field. ; 14. Instrument response normalization to physical units. ; 15. Output to a FITS file. (See .) ; ; ; NOTES: ; ; 1) There are three levels of verbosity. ; a) "verbose" = highest priority. All errors and messages are ; displayed. ("if q_vb") ; b) "quiet" = lower priority. No errors or messages are ; displayed. ("if q_qt") ; c) neither = lowest priority. All errors and some messages are ; displayed. ("if not q_qt") ; ; 2) The "index_out.HISTORY" is updated to reflect the various ; calibrations performed on the data. ; ; 3) (v2007.May.16) There is a subroutine called ; that searches for certain known read-out signals to erase. ; If a given feature is not present in an image, you may see the ; following error messages: ; % CURVEFIT: Failed to converge ; and/or ; % Program caused arithmetic error: Floating divide by 0 ; % Program caused arithmetic error: Floating underflow ; % Program caused arithmetic error: Floating overflow ; % Program caused arithmetic error: Floating illegal operand ; This does not pose a problem for the calibration, but may be ; disconcerting to the end-user. This error reporting will be ; improved in the future. ; ; 4) The "despike" module identifies "spikey" pixels in the image ; and applies a median filter process to flatten them, ; under the assumption that these represent high-energy ; particle hits in the detector. Therefore, this keyword ; should be used cautiously when prepping display images, ; and possibly should not be used at all for quantitative ; data analysis. See for more information. ; ; 5) The dust, dust growth, and hot pixel information output in GRADE_MAP ; is preliminary and will be updated as better data become available. ; Note that not all "hot" pixels are equally "problematic" ; and their degree of deviation appears to be temperature and exposure ; time related (higher CCD temperatures and longer exposures ; cause more "hot" pixels). The dust growth pixels are based on ; a model of the growth rate, and are also somewhat uncertain. ; Most recent data may not have growth estimates yet tabulated; ; consult the index_out.history to check if this is the case. ; For more details, see xrt_pixel_grade.pro. ; ; 6) The spot locations given here indicate the maximum extent of their ; effect (i.e., in Al-mesh or G-band); thicker filters show ; progressively weaker effects, with some spots having little or ; no effect. Since the contaminant is primarily a long wavelength ; absorber, the spots are often less prominent (fractionally) when ; illuminated by hotter (active region/flare) plasma. ; The spots are basically transparent for filters thicker than ; Be-thin. ; ; 7) The dark subtraction will be incrementally improved in ; the near future as further calibration is completed. ; ; 8) The data output from XRT_PREP.pro (i.e., "image_out" parameter) ; defaults to I*2 for historical and memory-conservation reasons. ; Users are recommended to apply the /FLOAT keyword always unless ; they are constrained by available memory. ; ; 9) For the moment, users wanting to use the "best" (default) ; dark subtraction should have the environment variable, XRT_ARCHIVE, ; properly set to find the dark files needed for the zero point ; correction of the model dark. XRT_ARCHIVE should be set to the base ; of the directory tree with the level0 hinode xrt data. Example: ; ; if the path to your level0 data is: ; ; /archive/hinode/xrt/level0 ; ; then set XRT_ARCHIVE to: ; ; /archive/hinode ; ; To avoid need for archival darks, ; xrt_prep can be run in pure "model dark" mode, e.g., ; ; IDL> xrt_prep, index, data, index_out, image_out, /norm, $ ; /float, dark_type=2 ; ; (see also examples above). ; If the user does NOT have a local XRT data archive, a new program ; can be used to set proper environment ; variables to permit copying over the network of the dark files ; needed for best dark correction, eg, ; ; IDL> xrt_search_network,/enable ; IDL> xrt_prep, index, data, index_out, image_out, /norm, /float ; ; see for more details. ; ; 10) 2048x2048 darks became rare after the high-gain antenna failure; ; 512x2048 "strip darks" are now used to approximate them for ; use in the default dark correction (median-dark adjusted model, ; dark_type=0) ; ; 11) Changed how program treats missing data, now resetting it to a ; clearly bogus large negative value (-999). ; ; 12) Code now can also provide calibration uncertainties (uncert_map), ; and a cosmetic correction for contamination spots (/despike_despot ; [default], or /only_despot). In some cases it is useful to force ; correction of all spots; to do this set spotthresh=-1. ; ; 13) Code can now improve alignment by coaligning using UFSS (Ultra Fine ; Sun Sensor) data or UFSS + cross-correlation with AIA 335. Note ; that the latter choice may work best in some situations and not ; in others, depending on the average plasma temperature the ; code is focusing on compared to AIA 335 (so use with care!). ; Default is to use UFSS only (*not* to coalign to AIA 335 as well). ; ; ; CONTACT: ; ; Comments, feedback, and bug reports regarding this routine may be ; directed to this email address: ; xrt_manager ~at~ head.cfa.harvard.edu ; ; MODIFICATION HISTORY: ; progver = 'v2006.Nov.23' ;--- (JWC) Written. (Based heavily on ; .) progver = 'v2006.Nov.30' ;--- (MW) Systemic changes. progver = 'v2007.Feb.14' ;--- (MW) Added RO correction block. Modified ; dark subtraction to use constant DC rate. ; Added "no_rectify" keyword and function. ; Removed functionality related to dark files. progver = 'v2007.Mar.06' ;--- (MW) Added functionality for "sat_map". ; Changed saturation limit to 2500. ; Reintroduced update to exptime if normalized. progver = 'v2007.Mar.30' ;--- (MW) Fixed "input2" for Input Case #1. ; Wasn't selecting correct files, or accepting ; "input2" as a Byte or Integer type array. progver = 'v2007.Apr.11' ;--- (MW) Fixed functionality for "miss_map". ; Also adjusted "dc_rate" to 0.1/59. progver = 'v2007.May.07' ;--- (MW) Fixed functionality for "n_sat_pixels" ; and "n_miss_pixels". progver = 'v2007.May.10' ;--- (MW) Added Note #3. progver = 'v2007.May.12' ;--- (MW) Removed the constant pedestal kludge. ; Now the CCD pedestal is corrected as part of ; the "ramp" model done be . ; Added RF1 tags. Switched exposure references ; from EXPTIME to E_ETIM. progver = 'v2007.May.16' ;--- (MW) Removed the rectify/no_rectify ; functionality. progver = 'v2007.May.17' ;--- (MW) Fixed /NORMALIZE bug that rounded off ; too much. progver = 'v2007.May.19' ;--- (MW) Changed the DESPIKE module to ; , which is based on the ; way that ANA despikes. Added ; SENS_DESPIKE keyword, too. progver = 'v2008.May.18' ;--- (SSaar) Improved large-scale read-out signal ; removal. Included a vignetting correction. ; Added more options/functionality at the top ; level for background and Fourier cleaning. ; Removed older, redundant dark correction. progver = 'v2008-Aug-20' ;--- (K.Reeves (SAO)) Added update to EXCCDEX keyword ; for normalized dark frames progver = 'v2009-Jul-11' ;--- (SSaar (SAO)) Dark/read-noise improved by using ; median of contemporaneous darks to adjust model ; dark zero point. Orbit and CCD temperature ; term in dark model merged. Saturation map expanded ; to include saturation bloom/bleed, contamination ; spot, dust, and hot pixel information (now ; called grade_map). More flexibility added at the ; top level (dark_type and grade_type options). progver = 'v2009-Sep-14' ;--- (SSaar (SAO)) Fixed some bugs related to dark ; subtraction in cases when there aren't sufficient ; darks nearby in time to make a good median image. progver = 'v2009-Oct-10' ;--- (SSaar (SAO)) Fixed bug related to calculation ; of n_grade_pixel. Public release of Jul-11 ; and later upgrades. progver = 'v2009-Oct-19' ;--- (SSaar (SAO)) Added some extra header information ; especially with reference to XRT_ARCHIVE ; environment variable (note #9). progver = 'v2009-Oct-22' ;--- (M.Weber (SAO)) Updated header documentation to ; better explain the usage of /FLOAT. progver = 'v2009-Nov-15' ;--- (SSaar (SAO)) Fixed bug related to operating ; with dark_type=2 progver = 'v2010-Aug-19' ;--- (SSaar (SAO)) Implemented use of ; 512x2048 "strip darks" for more recent data. ; Added note #10 and a history notation if ; strip darks are used. Fixed bug related ; to operation when the catalog is not updated. ; Reduced error messages when appropriate dark ; data files aren't found. progver = 'v2010-Aug-19' ;--- (PGrigis,SSaar (SAO)) Implemented ability to run ; without a local XRT data archive, using ; to permit copying of ; files from a remote archive. ; Modified note #9 to explain use when running ; without a local archive. progver = 'v2011-Jan-10' ;--- (SSaar (SAO)) Changed how missing data is flagged ; (now set = -999) and changed corresponding ; procedure note #2 above. Improved Fourier ; filtering when there is missing data. Modified ; note 9. Corrected sens_spike in header. progver = 'v2013-Aug-25' ;--- (SSaar (SAO)) Improved use of medianed darks. ; Added default cosmetic correction of contamination ; spots and dust, new keywords /despike_despot, ; /only_despot, /only_despike, and ; spotthresh, replacing old /despike keyword. progver = 'v2013-Aug-25' ;--- (AKobelski (MSU), SSaar (SAO)) Implemented ; ability to obtain the photometric calibration ; uncertainty via the uncert keyword. progver = 'v2013-Oct-31' ;--- (SSaar (SAO)) Added estimate for dust growth; ; modified Note 5. Added ability to output ; both unmodified and cosmetically corrected output. progver = 'v2013-Nov-08' ;--- (SSaar (SAO)) Added a fix for correcting all ; spots, enhanced procedure writeup, examples, and ; note 12 above progver = 'v2013-Nov-19' ;--- (SSaar (SAO)) Added a fix for readport = 'L' ; cases progver = 'v2014-Jan-15' ;--- (SSaar (SAO)) Added improved spike removal, and ; info above on spike_map. Added quiet to ; uncertainty calculation. progver = 'v2014-Oct-20' ;--- (SSaar (SAO)) Added coalignment ability and ; metadata index update. progver = 'v2017-Jun-27' ;--- (KY (MSU)) Replaced xrt_med_dark.pro with ; msu_xrt_med_dark.pro to get around the problem ; in finding dark data. Added "nosaturation" option ; to turn off "sat_threshold" cut-off. Added ; "darkff" option to re-use dark frame. ; ; ;- ; ========================================================================= ;=== Initial setup ==================================================== ;===== Track program runtime. start_time = systim(1) ;===== Initialize program constants. prognam = 'XRT_PREP' ; sat_threshold = 2500 ; sat_value = 2500 despike_limit = 3 miss_data_val = -999 ; value to reset pixels w/missing data, post-prep dc_rate = 0.1 / 59 ; [12-bit ADU per pixel per sec @ -65C] base_dir= getenv('SSWDB') + '/hinode/xrt/xrt_msu_coalign/' ;===== Set Booleans which control print statements. ;===== Keyword "verbose" overrules "quiet". q_vb = keyword_set(verbose) q_qt = keyword_set(quiet) and (not q_vb) ;===== Announce version of . if (not q_qt) then box_message, 'XRT_PREP: Running ' + prognam + ' ' $ + progver + '.' ;===== Set some keyword Booleans, plus miscellaneous. q_despike_despot = keyword_set(despike_despot) q_despot = keyword_set(only_despot) or q_despike_despot q_despike = keyword_set(only_despike) if q_despike then despike=only_despike else despike=0 ; set despike if q_despike_despot then despike=despike_despot ; despike_despot overrides q_despike = q_despike or q_despike_despot ; generalize q_despike q_norm = keyword_set(normalize) q_float = keyword_set(float) q_miss_map = keyword_set(miss_map) q_spike_map = keyword_set(spike_map) and q_despike q_spotthresh = keyword_set(spotthresh) q_coal = keyword_set(coalign) qstop = keyword_set(qstop) qabort = 0B allspot=0 ;;;; MSU localization (start) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; d_time_thre=3600.*24*15 d_time=0 ; Time of Dark Frame (initialize) if keyword_set(nosaturation) then begin sat_threshold = 4096 sat_value = 4096 endif else begin sat_threshold = 2500 sat_value = 2500 endelse if keyword_set(darkff) then dark_type=5 ;;;; MSU localization (end) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if n_elements(uncert_map) eq 0 then begin uncert_map=0 q_uncert=0 endif else q_uncert=1 if not keyword_set(outdir) then outdir = curdir() if n_elements(grade_map) eq 0 then grade_map=0 if keyword_set(grade_type) eq 0 then grade_type=0 ; killed blank space[SS] if grade_type ne 2 then begin grade_frac=grade_type q_grade_map=1 endif else q_grade_map=0 if keyword_set(dark_type) eq 0 then dark_type=0 if dark_type eq 1 then bsub_type=3 if q_spotthresh eq 1 then begin if spotthresh lt 0 or spotthresh gt 1 then begin if spotthresh ne -1 then begin if q_vb then print,'spot threshold out of range .. resetting to 0.5' endif else allspot = 1 spotthresh = 0.5 endif endif if q_coal eq 1 then coal = coalign else coal = 0 ;=== Check inputs ===================================================== ;===== This block checks to see whether the parameters "input1" and ;===== "input2" are recognized as Case 1 or as Case 2. If yes, set the ;===== Boolean "q_incase1" appropriately. If no, then abort. if q_vb then print, 'XRT_PREP: Checking input parameters...' if (n_params() lt 4) then begin if (not q_qt) then box_message,'XRT_PREP: Requires >= 4 parameters.' $ + ' Aborting...' qabort = 1B return ;=== If n_params() ge 4. endif else begin in1_type = datatype(input1) in2_type = datatype(input2) ;=== Is input Case 1? if ((in1_type eq 'STR') and ((in2_type eq 'BYT') or $ (in2_type eq 'INT') or (in2_type eq 'LON'))) then begin q_incase1 = 1B if q_vb then print, 'XRT_PREP: Recognize inputs as ' $ + 'filenames + indices.' ;=== Do files exist? Nexist = total(file_exist(input1[input2])) if (Nexist ne n_elements(input2)) then begin if (not q_qt) then box_message, $ ['XRT_PREP: Not all of the input files exist. Aborting...', $ ' Try checking your inputs with something like this:', $ ' IDL> print, where(file_exist(input1[input2]) eq 0)'] qabort = 1B return endif ;=== Establish useful quantities. Nimg = n_elements(input2) read_xrt, input1[input2], index0, /nodata, verbose=verbose, $ quiet=quiet, qabort=qbort, qstop=0 if qabort then begin if (not q_qt) then box_message, 'XRT_PREP: Problem reading ' $ + 'the input files. Aborting...' return endif ;=== Is input Case 2? endif else begin ;=== Ensure input2 is of a number type. num_check = total(is_number(input2)) eq n_elements(input2) ;=== Ensure input2 has 2 or 3 dimensions. dim_check = (size(input2,/n_dim) eq 2) or (size(input2,/n_dim) eq 3) if ((in1_type eq 'STC') and num_check and dim_check) then begin q_incase1 = 0B if q_vb then print, 'XRT_PREP: Recognize inputs as ' $ + 'index + data array.' ;=== Do "input1" and "input2" match on the number of images? Nimg = n_elements(input1) if (size(input2,/n_dim) eq 3) $ then Nimg2 = n_elements(input2[0,0,*]) $ else Nimg2 = 1 if (Nimg ne Nimg2) then begin if (not q_qt) then box_message, ['XRT_PREP: Input ' $ + 'parameters do not match on the number of images.', $ 'Aborting...'] qabort = 1B return endif ;=== Check whether the data array is at least 8x8 in size. Nx = n_elements(input2[*,0,0]) Ny = n_elements(input2[0,*,0]) if ((Nx lt 8) or (Ny lt 8)) then begin if (not q_qt) then box_message, 'XRT_PREP: Input data ' $ + 'images must be at least 8x8 in size. Aborting...' qabort = 1B return endif ;=== Establish useful quantities. index0 = input1 ;=== Not input Case 1 or 2. endif else begin if (not q_qt) then box_message, 'XRT_PREP: Input parameters ' $ + 'not recognized. Aborting...' qabort = 1B return endelse endelse endelse if q_vb then print, ' ...OK' ;===== Verify that index0 has necessary tags. if q_vb then print, 'XRT_PREP: Checking whether this data has ' $ + 'necessary XRT header tags...' if not required_tags(index0, 'HISTORY,NAXIS1,NAXIS2,E_ETIM,' + $ 'PLATESCL', missing=missing) then begin if (not q_qt) then box_message, ['XRT_PREP: This input header ' $ + 'does not have necessary XRT tags:', $ ' ' + missing, $ ' Aborting...'] qabort = 1B return endif if q_vb then print, ' ' $ + ' ...OK' ;===== This block checks whether the data has already been prepped. if q_vb then print, 'XRT_PREP: Checking whether this data has ' $ + 'already been prepped...' hist0 = get_history(index0, 'XRT_PREP', found=found) if found then begin if (not q_qt) then box_message, ['XRT_PREP: Some of this data ' $ + 'has already been run through', $ ' . ' $ + 'Check header tag "history". Aborting...'] prepped=1 endif else prepped=0 if q_vb then print, ' ' $ + ' ...OK' if prepped eq 1 then begin if (not q_qt) then box_message, ['XRT_PREP: ' $ + 'requires unprepped input.', $ + 'Aborting grade mapping...'] grade_type=2 ; turn off grade mapping q_grade_map=0 endif if max(anytim(index0.date_obs)) ge anytim('2007-JUL-28 00:00') and $ (not q_qt) and (not q_despot) then $ box_message, ['XRT_PREP: This data has contamination' $ + ' spots. It is strongly recommended that you ' $ + 'run with /despike_despot or ' $ + '/only_despot (for cosmetic correction) and ' $ + 'grade_map options (to get a spot map).'] ;===== This block checks the "despike" keywords. if q_despike then begin if q_vb then print, 'XRT_PREP: Checking "despike" keywords ...' ;=== If "despike" isn't within the valid range: if ((fix(despike) lt 1) or $ (fix(despike) gt despike_limit)) then begin if (not q_qt) then box_message, ['XRT_PREP: Keyword ' $ + '"despike" has an invalid value.', $ ' Setting to default value = 1'] despike = 1 ;=== Else "despike" is valid: endif else begin if q_vb then print, ' ...OK' endelse endif ;=== Calibration loop ================================================= ;=== MAIN LOOP begins here. if q_vb then print, 'XRT_PREP: Starting calibration loop...' mdark0=0. ; initializing median dark grade_map_old=0. ; initializing prev. grade map for ii = 0, (Nimg-1) do begin loop_start = systim(1) if q_vb then print, 'XRT_PREP: Loop #' + strcompress(ii,/rem) $ + ' of ' + strcompress(Nimg-1,/rem) + '.' tagval0 = '' ;===== This block establishes a single index0/image0 to be processed. ;=== Is input Case 1? if q_incase1 then begin read_xrt, input1[input2[ii]], index0, image0, verbose=verbose, $ quiet=quiet, qabort=qabort, qstop=0 if qabort then begin aa0 = strcompress(ii,/rem) if (not q_qt) then box_message, 'XRT_PREP: Problem reading ' $ + 'Image #' + aa0 + '. Aborting...' return endif Nx = gt_tagval(index0, /naxis1) Ny = gt_tagval(index0, /naxis2) ;=== Input Case 2. endif else begin index0 = input1[ii] Nx = gt_tagval(index0, /naxis1) Ny = gt_tagval(index0, /naxis2) image0 = input2[0:Nx-1,0:Ny-1,ii] endelse ;=== Do calibrations with floats. image0 = float(temporary(image0)) if q_uncert then image_init=image0 ; this is necessary for xrt_unc ; === readport = L fix if index0.readport eq 'L' then begin index0.p2col = index0.p2col - index0.naxis2 +1 index0.p1col = index0.p1col - index0.naxis2 +1 endif ;===== This block replaces missing pixels with a constant miss_data_val ;=== First check to see if image0 has already been processed. previous = get_history(index0, 'XRT_MISSING_PIXELS', found=found) if not found then begin if q_vb then print, 'XRT_PREP: Replacing missing pixels...' if q_miss_map then miss_map0 = byte(image0) * 0B ;=== Data values < 0? if (min(image0) lt 0) then begin if (not q_qt) then box_message, $ ['XRT_PREP: Image ' + strcompress(ii,/rem) + ': Data values < 0', $ ' ==> cannot determine if there are missing pixels.', $ ' Skipping this image for missing pixel replacement.'] ;=== Make note for the index0.HISTORY tag. tagval0 = [tagval0, '(XRT_MISSING_PIXELS) Negative data values ' $ + 'precluded determining presence of missing pixels.' $ + ' Skipped.'] n_miss_pixel0 = 0 ;=== Else all data ge 0, so can determine missing pixels. endif else begin ;=== Find missing pixels. ss_miss = where(image0 eq 0, n_miss_pixel0) ;=== temporarily replace missing pixels w/linear patch if not *all* missing. if ((n_miss_pixel0 gt 0) and $ (n_miss_pixel0 lt n_elements(image0))) then begin xyi0=array_indices(image0,ss_miss) ; find unique rows with 0s iuy0=uniq(xyi0(1,*)) ; (assumption: 0s come in horiz iy0=reform(xyi0(1,iuy0)) ; blocks [*,n]) - find 0 blocks iseg=[where(abs(iy0-shift(iy0,1)) gt 1,ns),n_elements(iy0)] ; avim=avg(image0,0) for j=0,ns-1 do begin ; loop thru blocks of 0s iy0a=(iy0(iseg(j))-1)>0 iy0b=(iy0(iseg(j+1)-1)+1)<(Ny-1) fa=image0(*,iy0a) fb=image0(*,iy0b) dy=iy0b-iy0a-1 + (iy0b eq (Ny-1)) + (iy0a eq 0) ; slab length in y df = (fb - fa)*(fa ne 0)*(fb ne 0) ; df, insuring fa & fb<>0 sl= df/dy ; slope @ each x f0 = fa*(fa ne 0) + fb*(fa eq 0) ; fa, insuring fa<>0 y=findgen(dy) ; construct a linear patch slab = sl#y + f0#(fltarr(dy)+1.) ; over y @ each x in block image0(0,iy0a+1-(iy0a eq 0))=slab ; insert interpolated patch endfor ; the patch improves fourier ; filter performance and is later ; replaced by miss_data_val if q_miss_map then miss_map0[ss_miss] = 1B endif ;=== Make note for the index0.HISTORY tag. aa0 = strcompress(n_miss_pixel0,/rem) aa1 = strcompress(miss_data_val,/rem) tagval0 = [tagval0, '(XRT_MISSING_PIXELS) Replaced ' + aa0 + $ ' missing pixels with ' + aa1 + '.'] endelse if q_vb then print, ' ...OK' ;==== Else there is a prior history of fixing missing pixels. endif else begin aa0 = strcompress(ii,/rem) if (not q_qt) then box_message, 'XRT_PREP: Image ' + aa0 + $ ': Already processed for XRT_MISSING_PIXELS. Skipping.' n_miss_pixel0 = 0 endelse ;===== This block replaces saturated pixels. ;=== First check to see if image0 has already been processed. previous = get_history(index0, 'XRT_SATURATED_PIXELS', found=found) if not found then begin if q_vb then print, 'XRT_PREP: Replacing saturated pixels...' grade_map0 = byte(image0) * 0B ;=== Data values > sat_threshold? ss_sat = where(image0 gt sat_threshold, n_sat_pixel0) if (n_sat_pixel0 gt 0) then begin ;=== Replace saturated pixels with large value, so that ;=== 1x1, 2x2 & 4x4 are still >4095 after dark subtraction.(?!?) image0[ss_sat] = sat_value grade_map0[ss_sat] = 1B endif ;=== Make note for the index0.HISTORY tag. aa0 = strcompress(n_sat_pixel0,/rem) aa1 = strcompress(sat_value,/rem) tagval0 = [tagval0, '(XRT_SATURATED_PIXELS) Replaced ' + aa0 + $ ' saturated pixels with value = ' + aa1 + '.'] if q_vb then print, ' ...OK' ;==== Else there is a prior history of fixing saturated pixels. endif else begin aa0 = strcompress(ii,/rem) if (not q_qt) then box_message, 'XRT_PREP: Image ' + aa0 + $ ': Already processed for XRT_SATURATED_PIXELS. Skipping.' endelse ;===== This block does despiking of radiation hits. ;=== Despike an x-ray image to remove bright bad pixels due to ;=== rad-belt/cosmic-ray hits. ;=== (Note that there is no check against the despike HISTORY.) if q_despike then begin if prepped eq 0 then begin if q_vb then print, 'XRT_PREP: Performing despiking...' image1 = image0 spike_map0 = bytarr(Nx, Ny) for jj = 1,despike do begin image1 = xrt_despike2(temporary(image1), index0, sens=sens_despike, $ hist=hist, spike_map=spike_map1, $ verbose=q_vb, quiet=q_qt, qabort=qabort) ;=== If there was a problem with on this loop: if qabort then begin aa0 = strcompress(ii,/rem) if (not q_qt) then box_message, 'XRT_PREP: Image ' + aa0 + $ ': Unable to perform despiking. Skipping.' ;=== Break out of the loop. jj = despike + 10 ;=== Else there was no problem on this loop: endif else begin if q_spike_map and spike_map1[0] ne -1 then $ spike_map0[spike_map1] = 1B tagval0 = [tagval0, hist] endelse endfor ;=== If there was no "qabort" problem: if (not qabort) then begin image0 = image1 ss = where(spike_map0 eq 1, n_spike_pixel0) ;=== Make note for the index0.HISTORY tag. aa0 = strcompress(n_spike_pixel0,/rem) aa1 = strcompress(despike,/rem) tagval0 = [tagval0, '(XRT_DESPIKING) Successfully despiked ' + $ aa0 + ' pixels in ' + aa1 + ' passes.'] if q_vb then print, ' ...OK' ;=== Else there was a "qabort" problem: endif else begin qabort = 0B ;=== Make note for the index0.HISTORY tag. tagval0 = [tagval0, '(XRT_DESPIKING) Was unable to ' $ + 'complete despiking. Skipped.'] endelse endif else begin ; data already prepped aa0 = strcompress(ii,/rem) if (not q_qt) then box_message, 'XRT_PREP: Image ' + aa0 + $ ': Already prepped. Run to despike. Skipping.' endelse ;=== Else no despiking requested. endif else begin endelse ;===== This block deals with flat-fielding. if qstop then stop ;===== This block corrects various known read-out signals. ;=== First check to see if image0 has already been processed. previous = get_history(index0, 'XRT_RO_CORRECTION', found=found) if not found then begin stripdrk=0 ; set flag for strip darks ;;;; MSU localization (start) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; if dark_type ne 2 then mdark = $ ; compute medianed contemporaneous dark ; msu_xrt_med_dark(index0, ddates=dates, darkff=darkff, stripdark=stripdrk, $ ; dfile=dfile, quiet=q_qt) else mdark=0 ; if dark_type<>2, else mdark=0 if dark_type lt 2 then begin if n_elements(index1) ne 0 then begin if abs(d_time-anytim(index0.date_obs)) le d_time_thre and $ index0.chip_sum eq index1.chip_sum and $ index0.p1col eq index1.p1col and $ index0.p2col eq index1.p2col and $ index0.p1row eq index1.p1row and $ index0.p2row eq index1.p2row then begin mdark=darkff endif endif else begin mdark = msu_xrt_med_dark(index0, ddates=dates, darkff=darkff, $ nday=2,stripdark=stripdrk,dfile=dfile, quiet=q_qt) if n_elements(dates) ne 0 then $ ; Try wider search mdark = msu_xrt_med_dark(index0, ddates=dates, darkff=darkff, $ stripdark=stripdrk,dfile=dfile, quiet=q_qt) if n_elements(dates) ne 0 then begin d_time=average(anytim(dates)) ; Update d_time index1=index0 delvarx,dates endif endelse endif else begin if dark_type eq 5 then begin i_sz=size(image0) d_sz=size(darkff) if max(abs(i_sz(1:2)-d_sz(1:2))) eq 0 then mdark=darkff $ else mdark=0 endif else mdark=0 endelse ;;;; MSU localization (end) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if q_vb then print, 'XRT_PREP: Correcting for known', $ ; contemp. dark + ' readout signals...' if total(mdark) eq 0 then begin dark_typ=2 ; if no median dark, use model if dark_type eq 1 then bsub_type=0 ; and reset bsub_type endif else begin if dark_type eq 5 then dark_typ=0 else dark_typ=dark_type endelse xrt_clean_ro, image0, index0, mdark, image1, grade_map0, model, $ hist=hist_ro, clean_type=clean_type, bsub_type=bsub_type, $ dark_type=dark_typ, nsigma=nsigma, nmed=nmed, $ stop_bleed=stop_bleed, verbose=verbose, quiet=q_qt, $ qabort=qabort ;stop ; ***HERE if grade_type ne 2 then begin ; if grade is of interest if total(mdark-mdark0) ne 0 or (ii eq 0) then begin ; and if mdark ; has changed or is new xrt_pixel_grade,index0,grade_map0,g_interp,ngrade0, $ ; redo maps grade_frac=grade_frac ; for spot,dust,hot pix endif else begin ; otherwise reuse old one del = 1d-8 ; bump to fix round-off err xrt_pixel_grade,index_old,grade_map_old,g_interp,/interpret ; get sat_old = g_interp(*,*,0) ; sat,bloom from last bloom_old = g_interp(*,*,1) ; grade_map, ngrade0([0,1])=[total(abs(grade_map0-1) le del), $ total(abs(grade_map0-2) le del)] ; modify old ngrade0 grade_map0 = grade_map0 + grade_map_old $ ; add old spot,dust,hot - sat_old - bloom_old*2 ; but remove old sat,bloom if grade_type eq 1 then grade_map0 = grade_map0 + $ ; & add fracs (grade_map_old - fix(grade_map_old + del)) else $ grade_map0 = fix(grade_map0 + del) ; or not endelse endif mdark0=mdark ; set current medianed dark="standard" grade_map_old=grade_map0 ; set old grade map index_old=index0 ; set old index ;=== If there was no "qabort" problem: if (not qabort) then begin image0 = image1 ;=== Make note for the index0.HISTORY tag. if stripdrk eq 1 then hist_ro=[hist_ro, $ ; note use of strip dark '512x2048 strip darks used -- zero-point less well determined.'] tagval0 = [tagval0, hist_ro, '(XRT_RO_CORRECTION) Successful.'] if q_vb then print, ' ' $ + ' ...OK' ;=== Else there was a "qabort" problem: endif else begin qabort = 0B model=1. ;=== Make note for the index0.HISTORY tag. tagval0 = [tagval0, '(XRT_RO_CORRECTION) Was unable to ' $ + 'complete correction of readout signals. Skipped.'] endelse ;==== Else there is a prior history of correcting the readout signals. endif else begin aa0 = strcompress(ii,/rem) model=1. q_grade_map=0 if (not q_qt) then box_message, 'XRT_PREP: Image ' + aa0 + $ ': Already processed for XRT_RO_CORRECTION. Skipping.' endelse ;===== This block does the vignetting correction. previous = get_history(index0, 'XRT_VIGNETTING', found=found) if (not found) then begin if q_vb then print, 'XRT_VIGNETTING: Removing vignetting ...' image0 = nono_vignette(image0, index0, vmodel, hist=hist) model = model/vmodel ;=== Make note for the index0.HISTORY tag. tagval0 = [tagval0, hist] if q_vb then print, ' ...OK' ;=== Else there is a prior history of vignetting correction. endif else begin aa0 = strcompress(ii,/rem) if (not q_qt) then box_message, 'XRT_PREP: Image ' + aa0 + $ ': Already processed for vignetting. Skipping.' endelse ;===== This block calculates the photometric uncertainty if q_uncert then begin if q_vb then print, 'XRT_PREP: Calculating uncertainties...' uncert0=sqrt(xrt_unc_rel(index0, temporary(image_init), image0,vmodel, $ dark_type=dark_type,four_type=clean_type,strip_dark=stripdrk,$ quiet=q_qt)*image0^2) tagval0 = [tagval0, '(XRT_UNC_REL) Calculated photometric uncertainty.'] if q_vb then print, ' ...OK' endif ;===== This block does renormalization to DN / sec. ;=== Renormalize to DN per sec if "/normalize" is set. ;=== (Note that there is no check against the normalize HISTORY.) if q_norm then begin previous = get_history(index0, 'XRT_RENORMALIZE', found=found) if (not found) then begin if q_vb then print, 'XRT_PREP: Performing normalization...' if qstop then stop img_exp = get_xrt_expdur(index0, /sec) image0 = temporary(image0) / (img_exp) if q_uncert then uncert0 = temporary(uncert0) / (img_exp) ;;=== Dark exposure times are recorded in the exccdex keyword index0.exccdex = 1000000L ; [microsec]$ index0.e_etim = 1000000L ; [microsec] ;=== Make note for the index0.HISTORY tag. aa0 = strcompress(img_exp,/rem) tagval0 = [tagval0, '(XRT_RENORMALIZE) Normalized from ' + $ aa0 + ' sec --> 1.00 sec.'] if q_vb then print, ' ...OK' endif else begin aa0 = strcompress(ii,/rem) if (not q_qt) then box_message, 'XRT_PREP: Image ' + aa0 + $ ': Already normalized. Skipping.' endelse ;=== Else keyword "/normalize" NOT set. endif else begin endelse ;=== perform spot correction if q_despot then begin if q_vb then print, 'XRT_PREP: Performing cosmetic spot correction...' bin=index0.chip_sum if q_spotthresh eq 0 then begin case 1 of (bin eq 2): spotthresh = 0.5 (bin eq 4): spotthresh = 0.5 (bin eq 8): spotthresh = 0.5 else: spotthresh = 1. endcase endif if allspot eq 0 then begin xrt_spotcor,index0,image0,ind_out,dat_out,nspots,spotmap, $ thresh=spotthresh,quiet=q_qt endif else begin xrt_spotcor,index0,image0,ind_out,dat_out,nspots,spotmap, $ thresh=spotthresh,quiet=q_qt,/all endelse ; xrt_tup_contam,index0,image0,ind_out,dat_out,thresh=spotthresh index0 = ind_out image1 = image0 ; ***NEW image0 = dat_out if q_vb then print, ' ...OK' endif ;=== Coalign case 1 of (coal eq 0): begin ind_co0=xrt_read_coaldb(index0,base_dir=base_dir,cal=caltyp) aa0 = strcompress(caltyp,/rem) tagval0 = [tagval0, '(XRT_READ_COALDB) Coaligned using UFSS data.' +$ ' Calib. type: ' + aa0] end (coal eq 1): begin ind_co0=xrt_read_coaldb(index0,base_dir=base_dir,/aia_cc,cal=caltyp) aa0 = strcompress(caltyp,/rem) tagval0 = [tagval0, '(XRT_READ_COALDB) Coaligned using UFSS data' + $ '+ cross-correlation with AIA 335 data. Calib. type: ' + aa0] end else: tagval0 = [tagval0, '(XRT_READ_COALDB) Not coaligned.'] endcase if coal ne 2 then begin index0.xcen = ind_co0.xcen index0.ycen = ind_co0.ycen index0.crota1 = ind_co0.crota1 index0.crota2 = ind_co0.crota2 endif ;=== Replace missing pixels (if any) with "bogus" value miss_data_val if n_miss_pixel0 gt 0 then image0(ss_miss) = miss_data_val ;===== Now is the time to update the index0.HISTORY tag. if q_vb then print, 'XRT_PREP: Updating HISTORY...' hist_prefix = prognam + ' ' + progver + ': (' + systim() + ') ' hist_entry = hist_prefix + tagval0[1:*] update_history, index0, hist_entry, /noroutine ;===== Match up the number of HISTORY records so that index0 ;===== can be concatenated into index_out. ;=== If first image processed: if (ii eq 0) then begin index_out = index0 ;=== Else if NOT the first image processed: endif else begin hist_out_len = n_elements(index_out[0].history) hist0_len = n_elements(index0[0].history) ;=== If index_out is bigger: if (hist_out_len gt hist0_len) then begin aa0 = strarr(hist_out_len - hist0_len) update_history, index0, aa0, /noroutine endif else begin ;=== Or if index0 is bigger: if (hist_out_len lt hist0_len) then begin aa0 = strarr(hist0_len - hist_out_len) update_history, index_out, aa0, /noroutine ;=== Else index_out must be same size as index0: endif else begin endelse endelse ;=== Concatenate index0 onto index_out. index_out = concat_struct(temporary(index_out), index0) endelse ;===== Match up the size of the data arrays so that image0 ;===== can be concatenated into image_out. Also, pay attention ;===== to the "/float" keyword. ;=== Return image_out as Integer*2 unless "/float" is set. if (not q_float) then begin image0 = fix(round(temporary(image0))) ; This is the default. if q_despot eq 1 then image1 = fix(round(temporary(image1))) endif else begin image0 = float(temporary(image0)) if q_despot eq 1 then image1 = float(temporary(image1)) endelse ;=== If first image processed: if (ii eq 0) then begin image_out = image0 if q_despot eq 1 then image_out2 = image1 if q_uncert then uncert_map = uncert0 if q_grade_map then begin grade_map = grade_map0 n_grade_pixels = ngrade0 endif if q_miss_map then begin miss_map = miss_map0 n_miss_pixels = total(miss_map0) endif if q_spike_map then begin spike_map = spike_map0 n_spike_pixels = total(spike_map0) endif ;=== Else if NOT the first image processed: endif else begin Nx_out = n_elements(image_out[*,0,0]) Ny_out = n_elements(image_out[0,*,0]) Nx0 = n_elements(image0[*,0,0]) Ny0 = n_elements(image0[0,*,0]) ;=== First do x-dimension. ;=== If image_out (x) is bigger: if (Nx_out gt Nx0) then begin aa0 = Nx_out - Nx0 image0 = [temporary(image0), bytarr(aa0,Ny0)] if q_despot then image1 = [temporary(image1), bytarr(aa0,Ny0)] if q_uncert then uncert0 = [temporary(uncert0),$ bytarr(aa0,Ny0)] if q_grade_map then grade_map0 = [temporary(grade_map0), $ bytarr(aa0,Ny0)] if q_miss_map then miss_map0 = [temporary(miss_map0), $ bytarr(aa0,Ny0)] if q_spike_map then spike_map0 = [temporary(spike_map0), $ bytarr(aa0,Ny0)] Nx0 = n_elements(image0[*,0,0]) endif else begin ;=== Else if image0 (x) is bigger: if (Nx_out lt Nx0) then begin aa0 = Nx0 - Nx_out image_out = [temporary(image_out), bytarr(aa0,Ny_out,ii)] if q_despot then $ image_out2 = [temporary(image_out2), bytarr(aa0,Ny_out,ii)] if q_uncert then $ uncert_map = [temporary(uncert_map), bytarr(aa0,Ny_out,ii)] if q_grade_map then $ grade_map = [temporary(grade_map), bytarr(aa0,Ny_out,ii)] if q_miss_map then $ miss_map = [temporary(miss_map), bytarr(aa0,Ny_out,ii)] if q_spike_map then $ spike_map = [temporary(spike_map), bytarr(aa0,Ny_out,ii)] Nx_out = n_elements(image_out[*,0,0]) ;=== Else image_out (x) and image0 (x) are same size. endif else begin endelse endelse ;=== Now do y-dimension. ;=== If image_out (y) is bigger: if (Ny_out gt Ny0) then begin aa0 = Ny_out - Ny0 image0 = [[temporary(image0)], [bytarr(Nx0,aa0)]] if q_despot then $ image1 = [[temporary(image1)], [bytarr(Nx0,aa0)]] if q_uncert then $ uncert0= [[temporary(uncert0)], [bytarr(Nx0,aa0)]] if q_grade_map then $ grade_map0 = [[temporary(grade_map0)], [bytarr(Nx0,aa0)]] if q_miss_map then $ miss_map0 = [[temporary(miss_map0)], [bytarr(Nx0,aa0)]] if q_spike_map then $ spike_map0 = [[temporary(spike_map0)], [bytarr(Nx0,aa0)]] Ny0 = n_elements(image0[0,*,0]) endif else begin ;=== Else if image0 (y) is bigger: if (Ny_out lt Ny0) then begin aa0 = Ny0 - Ny_out image_out = [[temporary(image_out)], [bytarr(Nx_out,aa0,ii)]] if q_despot then $ image_out2 = [[temporary(image_out2)], [bytarr(Nx_out,aa0,ii)]] if q_uncert then $ uncert_map=[[temporary(uncert_map)], [bytarr(Nx_out,aa0,ii)]] if q_grade_map then $ grade_map = [[temporary(grade_map)], [bytarr(Nx_out,aa0,ii)]] if q_miss_map then $ miss_map = [[temporary(miss_map)], [bytarr(Nx_out,aa0,ii)]] if q_spike_map then $ spike_map = [[temporary(spike_map)], [bytarr(Nx_out,aa0,ii)]] Ny_out = n_elements(image_out[0,*,0]) ;=== Else image_out (y) and image0 (y) are same size. endif else begin endelse endelse ;=== Concatenate image0 onto image_out. image_out = [[[temporary(image_out)]], [[image0]]] if q_despot then $ image_out2 = [[[temporary(image_out2)]], [[image1]]] if q_uncert then $ uncert_map = [[[temporary(uncert_map)]], [[uncert0]]] if q_grade_map then begin grade_map = [[[temporary(grade_map)]], [[grade_map0]]] n_grade_pixels = [[n_grade_pixels],[ngrade0]] endif if q_miss_map then begin miss_map = [[[temporary(miss_map)]], [[miss_map0]]] n_miss_pixels = [n_miss_pixels,total(miss_map0)] endif if q_spike_map then begin spike_map = [[[temporary(spike_map)]], [[spike_map0]]] n_spike_pixels = [n_spike_pixels,total(spike_map0)] endif endelse ;===== Print processing information. aa0 = strcompress(ii,/rem) aa1 = string(systim(1)-loop_start, format='(F8.3)') if q_vb then print, 'XRT_PREP: Image ' + aa0 + ' took ' + aa1 + $ ' seconds to process.' ;===== End the MAIN LOOP of calibration. endfor ;===== Add metadata index_out = xrt_index2planinfo(temporary(index_out), /quiet) ;=== Finish up ======================================================== index_out.data_lev = 1 get_utc, utc, /ccsds index_out = add_tag(temporary(index_out), utc, 'DATE_RF1', $ index='DATE_RF0', /top_level) index_out = add_tag(temporary(index_out),get_logenv('$SSW_SITE'), $ 'ORIG_RF1', index='VER_RF0', /top_level) index_out = add_tag(temporary(index_out),progver,'VER_RF1', $ index='ORIG_RF1', /top_level) run_time = systim(1) - start_time aa0 = string(run_time, format='(F8.3)') aa1 = strcompress(Nimg,/rem) if (not q_qt) then box_message, 'XRT_PREP: Took ' + aa0 + $ ' seconds to process ' + aa1 + ' images.' if (not q_qt) then box_message, 'XRT_PREP: Finished.' if qstop then stop END ;======================================================================