C+ C NAME: C smei_base C PURPOSE: C Determines pedestal & dark current, finds measles C CALLING SEQUENCE: program smei_base C INPUTS: C SMEI data frames read from specified source (usually the C SMEI frame data base) C C smei_base -start_time= -stop_time= C -camera=<1,2,3> [NOT YET: -mode=<0,1,2>] C -source= -destination= C C and bracket the time period to be processed. C The times are specified either as orbit numbers or as a UT time in the C format YYYY_DOY_hhmmss (i.e. the standard SMEI time format). C MUST be specified. If is not specified than a C single orbit starting at is processed. C C The camera number (1,2 or 3) MUST be specified. C C is the source directory where SMEI frames are picked up. C C is the destination directory for the updated SMEI frames C OUTPUTS: C Updated SMEI frames as Fits files in C INCLUDE: include 'openfile.h' include 'smei_frm_layout.h' include 'smei_frm_basepar.h' include 'smei_frm_hdr.h' include 'math.h' include 'dirspec.h' include 'filparts.h' C CALLS: C smei_frm_measle, Str2Str, Time2Str, Int2Str, bOpenFile C BadR4, BadI4, smei_Time2Split, iFilePath, iSetFileSpec C iGetFileSpec, ArrR4Mask, ArrR4Mean, ArrR4AddConstant, smei_frm_read C smei_frm_base, smei_frm_getlist, smei_frm_pickup, smei_frm_clean C smei_orbit_time2, smei_orbit2, smei_hdr_str, smei_frm_saturated C smei_frm_stats, smei_frm_c3fudge, smei_orb_mkname, ForeignI4Arg C ForeignR4Arg, ForeignArgSet, smei_frm_write, smei_foreign, BadR4Test C Say, iOSDeleteFile, smei_cal_init, smei_cal_get, iFreeLun, smei_frm_fts C iScratchLun, cDbl2Str, cHideLogical C PROCEDURE: C smei_base is run on every individual SMEI frame to fill in three groups C of header entries in the Fits header: C C 1) base-related entries (pedestal, dark current, pattern names, etc.) C 2) optical axis entries (RA, Dec, and PA) C These depend on the Coriolis sky-to-spacecraft quaternions and C on the fixed camera-to-spacecraft quaternions C 3) unit vectors for Sun, Moon, Venus in the UCSD camera frame, and C eclipse information. These depend on the Coriolis orbital C elements and the quaternions. C C If smei_base is first run on a .nic.gz file (which does not contain C any of the above information) all header entries are calculated C and a new .fts.gz file is created (i.e. both header and data frame C are written). This will happen in the conversion of the SMEI C data base from .nic.gz to .fts.gz files. C C As of Mar-2005 we do not produce .nic.gz files from the C real-time SMEI data, but create .fts.gz files instead with the C flag SMEI__HDR_BASE_DONE set to zero. These Fits files C do not yet have the base headers filled (group 1 above), while C the unit vectors of Sun, Moon and Venus are probably based C on orbital elements that could be weeks old, and hence may C not be very accurate. If smei_base is first run on these files C then the zero value in SMEI__HDR_BASE_DONE will trigger the C calculation of all three groups of headers. If the source is the C same as the destination, then the existing Fits file is updated C with the header information, i.e. the header is only partially C updated, and the data array remains untouched. If the C destination is different from the origin then a new Fits C file is written. In both cases the flag SMEI__HDR_BASE_DONE C is set to one. C C If smei_base is run on a Fits file which already has the C SMEI__HDR_BASE_DONE flag set, then by default no action is taken C (not even when destination and source are different) unless the C version number of smei_base has been increased. C Overide this behavior by setting the -force keyword. This will C trigger a recalculation of all three groups of header entries. C C The -savefrm option will do the base calculation and save three C fts files into the destination directory: C - *_0.fts.gz: frame with constant pedestal+dark_current subtracted C - *_1.fts.gz: frame with pedestal and pattern subtracted C - *_2.fts.gz: frame with pedestal and pattern subtracted and measles removed C The original Fits file is never modified when -savefrm is set. C MODIFICATION HISTORY: C FEB-2005, Paul Hick (UCSD/CASS) C Complete overhaul of Andys version to allow direct processing C from the SMEI data base. C JUN-2005, Paul Hick (UCSD/CASS), V1.02. C nped_min for cam 1 and 2 in mode 2: from 115 to 105 C ped_median_excess for cam 1 and 2 in mode 2: from 2 to 4 C Added bSaveFrm option C JUL-2005, Paul Hick (UCSD/CASS), V1.03 C Changed ped_right_offset_cam(3,1:2) from bad to zero C Fixed bug in cal_name determination. Left-right pedestal offset C in ped_aux(SMEI__PED_RIGHT_OFFSET) cleared for each frame instead C of only when calibration pattern changes. C JUL-2005, Paul Hick (UCSD/CASS), V1.04 C Modification in smei_frm_ped. ped_aux(SMEI__PED_RIGHT_OFFSET) C is now updated only by frames that pass the pedestal test on C ped_aux(SMEI__PED_MEAN_REF)+ped_aux(SMEI__PED_MEAN_EXCESS). C There also is a warning printed if ped_mean_diff (difference C between left and right mean pedestal) rises above C ped_aux(SMEI__PED_MEAN_EXCESS). If this happens C ped_aux(SMEI__PED_RIGHT_OFFSET) is not updated, and the pedestal C calculation might get stuck in a bad state. C AUG-2005, Paul Hick (UCSD/CASS), V1.05 C smei_frm_read contained a bug that affected the cam3/mode1 frames C only. The frame time was not picked up correctly, and as a result C the onboard factor 2 gain change in Feb 2005 was not compensated C for correctly. C AUG-2005, Paul Hick (UCSD/CASS), V1.06 C Changed ndark_min_cam(3,1) from 100 to 450. C Added dark_max_ratio to list of command line arguments. C AUG-2005, Paul Hick (UCSD/CASS), V1.07 C In smei_frm_saturated the saturation value for frames where C the onboard flatfield is not enabled (primarily camera 3 mode 1 C data) was changed from 65535 to 65533. C SEP-2005, Paul Hick (UCSD/CASS), V2.00 C Residual values for squares and center group are now a mean C (instead of a median). C NOV-2005, Paul Hick (UCSD/CASS), V2.01 C dark_scale_ratio_cam is now initialized to BadInitFlag for C camera 3. This triggers calculation of the ratio using C the function smei_frm_ratio (called in smei_frm_base). C NOV-2005, Paul Hick (UCSD/CASS), V2.02 C dark_scale_ratio_cam is set to BadInitFlag only for mode 1 C (mode 0 and 2 are set to one). C dark_scale_fudge_cam is also set to BadInitFlag for mode 1. C This triggers calculation of fudge ratio with smei_frm_c3fudge. C NOV-2005, Paul Hick (UCSD/CASS), V2.03 C The fudge scale and offset are now applied whenever a dark C current value is available (even a rejected value), and not C only to dark currents of accepted frames. Note that this affects C camera 3 only (cam 1 and 2 do not have fudge factors). C MAR-2007, Paul Hick (UCSD/CASS), V3.10 C Set ped_aux(SMEI__PED_LEFT_SIDE_ONLY) to BadR4 to compensate C for changes to smei_frm_ped. C For camera 3 only after 2005, Doy 159 00 UT, the pedestal is C now derived from the left side (underscan) pedestal columns C only. Note that this also means that the left-side offset for C the pedestal is no longer tracked. C MAR-2007, Paul Hick (UCSD/CASS), V4.00 C For camera 3 only a check is made whether the "bad pixel" C mask is in effect. If it is then the calibration pattern with C the mask taken into account is used instead of the regular C mask. C For camera 3 the header entry for the orbital (on-the-fly) C pattern always used to be set to a name based on the orbit C start time for each frame. Now this entry is cleared (if C present) if the mode is not 1, or if the mode is 1 but the C the corresponding orbital pattern file does not exist. C (href=smei_orb= is now responsible for setting this entry). C JUN-2007, Paul Hick (UCSD/CASS), V4.10 C Bug fixes to smei_frm_read introduced in V4.00. C For time periods with a "bad pixel" mask onboard C (smei_frm_c3mask_active(0,ctime,-1,cmask)=1) C two errors were made: C - For camera 3, mode 1, the gain change of a factor two was C not applied to frames with the flat_enabled flag off (zero), C i.e. the pedestal was 2x the correct value (and all frames C were rejected as bad). C - For cameras 1 and 2 the headroom 0.75 was not corrected for, C when flat_enabled=1, i.e. the pedestal and dark current C were 0.75 times the correct value. C V4.00 was run for all cameras from 2006_001 to 2007_160. C JUL-2007, Paul Hick (UCSD/CASS) C Modified use of c3mask cmd line argument (used for camera 3 C in mode 1 only). C C If -c3mask=0 is specified then the regular calibration pattern C (with no "bad-pixel" mask taken into account) is used for all C frames. C If -c3mask=1 is specified then the calibration pattern with C the "bad-pixel" mask taken into account is used. This is only C useful for periods were a "bad-pixel" mask was in effect C (only for these periods does a "bad-pixel" calibration C pattern exist). C If -c3mask=-1 then the program decides which pattern to use C on a frame-by-frame basis (see href=smei_frm_c3mask_active=) C C If no explicit pattern is specified (i.e. the pattern from the C frame headers is used then -c3mask=-1 is the default. C C If an explicit pattern file is specified on the cmd line, then C only -c3mask=0 or -c3mask=1 are valid. The default is -c3mask=0. C JUL-2007, Paul Hick (UCSD/CASS; pphick@ucsd.edu), V4.11 C For camera 2 the pedestal is now calculated only from the C underscan (left-side) columns for frames after 2006_125. C JUN-2008, Paul Hick (UCSD/CASS) C Added -test cmd line argument (no frames written) C Addde -nofudge keyword (suppress fudge factors for dark current) C JUL-2008, Paul Hick (UCSD/CASS), V4.12 C Added Fits keyword DARK_AVE for mean of dark current. C (Keyword DARK_CUR contains the fudged median) C JUL-2008, Paul Hick (UCSD/CASS), V4.13 C Fudge factors for camera 3 with "bad-pixel" mask active C are now calculated differently (single fudge factor C for each mask). C OCT-2008, Paul Hick (UCSD/CASS), V4.14 C Significant changes to href=smei_frm_c3mask_active to fix some C problems with transitions between masks. Added new mask. C Added cmd line args -ped_median_deficit and -ped_median_excess C NOV-2008, Paul Hick (UCSD/CASS), V4.15 C Added cmd line args -ped_mean_ref and -ped_mean_excess. C Especially the default ped_mean_ref is not adequate anymore for C camera 3. The default ped_mean_excess is probably still OK. C NOV-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu), V5.00 C Removed cmd line arg -ped_mean_ref again. C Whenever some estimate for the pedestal is needed this C is now obtained by a call to function smei_frm_ped_mean. C This replaces hardcoded values in ped_mean_ref_cam (copied C to ped_aux(SMEI__PED_MEAN_REF)) and ped_init_cam (copied C to base_aux(SMEI__BAS_PED_INIT+[0,1,2]). C- character cSay*4 /'base'/ double precision version /5.00d0/ !=============== ! Function types logical ForeignArgSet logical bOpenFile integer Str2Str integer Time2Str character cDbl2Str*14 character cHideLogical*(FIL__LENGTH) logical smei_frm_read logical smei_frm_base integer smei_frm_getlist integer smei_frm_pickup integer smei_frm_saturated integer smei_frm_c3mask_active integer BadI4 logical BadR4Test !=============== character cArg *(FIL__LENGTH) character cStr *(FIL__LENGTH) character cFile*(FIL__LENGTH) character cBase*(FIL__LENGTH) character cDest*(FIL__LENGTH) character cList*(FIL__LENGTH) character cCalPattern*(FIL__LENGTH) character cal_name*(FIL__LENGTH) /' '/ character orb_name*(FIL__LENGTH) /' '/ ! Only set for cam3/mode1 real ccd_frame (SMEI__FRM_NPIX*3) real cal_pattern(SMEI__FRM_NPIX) real orb_pattern(SMEI__FRM_NPIX) integer mask (SMEI__FRM_NPIX) integer indx (SMEI__FRM_NPIX) real dark_pixs (SMEI__DRK_NPIX) real sqre_pixs (SMEI__SQR_NPIX) real cntr_pixs (SMEI__PIX_NPIX) equivalence (dark_pixs, sqre_pixs, cntr_pixs) real*8 hdr(SMEI__HDR_N) !========================= ! In-/output arrays for smei_frm_ped, smei_frm_dark, smei_frm_base real ped_aux (SMEI__PED_NUM) real dark_aux(SMEI__DRK_NUM) real base_aux(SMEI__BAS_NUM) ! Control variables logical bOverwrite logical bFrameOK logical bDigBase logical bDigDest logical bNoDups logical bForce logical bSaveFrm logical bNewCal logical bTest logical bNoFudge integer TBeg(2) integer TEnd(2) integer TFrm(2) integer TOrb(2) integer hms (4) double precision dorbit integer OPEN_ parameter (OPEN_ = OPN__REOPEN+OPN__STOP) integer PRINT_ parameter (PRINT_ = 100) !============================================================ ! Constraints on pedestal and dark current (camera dependent) ! For eng mode (no binning): 222.22 adu's = 1000 electrons real BadInitFlag parameter (BadInitFlag=-9999.0) ! ========================= ! Constants determining the pattern calculation. ! ped_mean_ref_cam: Andy's pedinit ! Reference pedestal used in smei_frm_ped in combination with ! ped_mean_excess_cam to detect 'crazy frames' with unrealistic ! high pedestal. !real ped_mean_ref_cam (3,0:2) /996.0, 1013.0, 1046.0, ! = ped_init_cam(1,*) !& 996.0, 1016.0, 1046.0, !& 996.0, 1016.0, 1046.0/ ! ped_mean_excess_cam: hardcoded in Andy's programs ! Excess above naive mean pedestal used in smei_frm_ped to detect 'crazy ! frames' with unrealistic high pedestal. real ped_mean_excess_cam (3,0:2) /3*50.0, 3*50.0, 3*50.0/ ! ped_median_deficit_cam, ped_median_excess_cam ! Deficit and excess, respectively, relative to naive median, used in ! smei_frm_ped to exclude low and high outliers from contributing to the ! pedestal calculation. The deficit is not applied to camera 3 or mode 0. real ped_median_deficit_cam (3,0:2) /BadInitFlag,BadInitFlag,BadInitFlag, & -10.0, -10.0,BadInitFlag, & -10.0, -10.0,BadInitFlag/ real ped_median_excess_cam (3,0:2) /4.0, 4.0, 4.0, & 2.0, 2.0, 4.0, & 4.0, 4.0, 4.0/ !========================= ! Constants that do not affect the pattern calculation, ! and only affect the performance test. ! ped_init_cam(0,*) = Andys pedinit2 (slightly bigger than pedinit) ! ped_init_cam(1,*) = Andys pedinit ! Used to reinitialize ped_last in smei_frm_base after more than ! 12 rejected frames. ! ped_init_cam(2,*) = Andys pedinit (cam 1 and 2) or pedinit2 (cam 3) ! Used to initialize ped_last in smei_frm_base. C real ped_init_cam (0:2,3,0:2) /1000.0, 996.0, 996.0, C & 1020.0, 1013.0, 1013.0, C & 1049.0, 1046.0, 1049.0, C & 998.0, 996.0, 996.0, C & 1020.0, 1016.0, 1016.0, C & 1051.0, 1046.0, 1046.0, C & 998.0, 996.0, 996.0, C & 1020.0, 1016.0, 1016.0, C & 1051.0, 1046.0, 1046.0/ ! dark_cut_cam: Andys patcr. In Andys program for cam 3 patcr is multiplied ! by 2 when the cut is applied. This factor of two is absorbed here in the value ! for dark_cut_cam itself. ! Used in smei_frm_dark to exclude pixels from contributing to the ! calculated dark current once the pattern is known. real dark_cut_cam (3,0:2) /44.44, 44.44, 88.88, & 6.94, 6.94, 11.10, & 6.94, 6.94, 11.10/ ! dark_power_cam: hardcoded in Andys programs. ! When the pattern is subtracted the dark current ratio of frame and pattern ! is scaled with this power. Note that this power is applied inconsistently: ! it is not used when pixels are excluded from the dark current calculation ! by applying dark_cut_cam. real dark_power_cam (3,0:2) / 2*1.5,1.0, 3*1.0, 3*1.0/ ! For camera 3 only, the pattern is used in the pedestal calculation to correct ! for a left/right asymmetry (see smei_frm_ped). Only applied for camera 3 in mode 0?? real ped_right_offset_cam (3,0:2) /2*BadInitFlag,0.0, & 2*BadInitFlag,0.0, & 2*BadInitFlag,0.0/ ! nped_min_cam, ndark_min_cam: hardcoded in Andys programs ! Threshold on minimum number of pixels required for validating a pedestal ! and dark current, respectively. integer nped_min_cam (3,0:2) /1800,1800,1800, & 450, 450, 450, & 105, 105, 115/ integer ndark_min_cam (3,0:2) / 1, 1, 684, & 100, 100, 450, & 100, 100, 1/ ! dark_scale_ratio_cam and dark_max_ratio_cam are used in smei_frm_base ! to control the ratio of the naive mean dark current and the pattern dark ! current. If set to BadInitFlag the ratio is calculated by calling ! smei_frm_ratio. real dark_scale_ratio_cam (3,0:2) /1.0, 1.0, 1.0, & 1.0, 1.0, BadInitFlag, & 1.0, 1.0, 1.0/ real dark_max_ratio_cam (3,0:2) /1.0, 1.0, 2.5, & 1.0, 1.0, 2.5, & 1.0, 1.0, 2.5/ ! Cam 3: 1.0 in November ! Ad hoc adjustments to dark current returned from smei_frm_base ! for camera 3 only. ! - scale change for 2x2 binned data: 0.76 ! - shift because Mark Cooke's 2x2 binning must differ from mine... real dark_offset_fudge_cam(3,0:2) /0.0, 0.0, 0.00, & 0.0, 0.0, -0.75, & 0.0, 0.0, -0.75/ real dark_scale_fudge_cam (3,0:2) /1.0, 1.0, 1.00, & 1.0, 1.0, BadInitFlag, & 1.0, 1.0, 0.76/ integer c3mask integer c3mask_ref double precision cal_version ! Pick up command line arguments call smei_foreign(3,cCalPattern,cBase,cDest,0,TBeg,TEnd,0d0,icam,mode,cArg) if (ForeignArgSet(cArg,'help')) then write (*,'(3X,A,T47,A)') cSay,' ', & ' ',' ', & ' '//cSwitch(:iSwitch)//'start_time=<'//SMEI__UT_FORMAT//'>' ,' ', & ' '//cSwitch(:iSwitch)//'stop_time=<'//SMEI__UT_FORMAT//'>' ,' ', & ' '//cSwitch(:iSwitch)//'camera=<1,2,3>' ,' ', & ' '//cSwitch(:iSwitch)//'mode=<0,1,2>' ,'default: 2 for cam 1,2; 1 for cam 3', & ' '//cSwitch(:iSwitch)//'source=' ,'default: SMEIDB?' , & ' '//cSwitch(:iSwitch)//'digsource' ,'default: off (on for SMEIDB?)', & ' '//cSwitch(:iSwitch)//'destination=' ,'default: $TUB' , & ' '//cSwitch(:iSwitch)//'digdestination' ,'default: off (on for SMEIDB?)', & ' '//cSwitch(:iSwitch)//'force' ,'default: off', & ' '//cSwitch(:iSwitch)//'nodups' ,'default: off', & ' '//cSwitch(:iSwitch)//'quiet=<0|1>' ,'default: 0', & ' '//cSwitch(:iSwitch)//'nped_min=' ,'default: depends on camera/mode', & ' '//cSwitch(:iSwitch)//'ndark_min=' ,'default: depends on camera/mode', & ' '//cSwitch(:iSwitch)//'dark_cut=' ,'default: depends on camera/mode', & ' '//cSwitch(:iSwitch)//'dark_scale_ratio=' ,'default: depends on camera/mode', & ' '//cSwitch(:iSwitch)//'dark_max_ratio=' ,'default: depends on camera/mode', & ' '//cSwitch(:iSwitch)//'ped_mean_excess=' ,'default: depends on camera/mode', & ' '//cSwitch(:iSwitch)//'ped_median_deficit=','default: depends on camera/mode', & ' '//cSwitch(:iSwitch)//'ped_median_excess=' ,'default: depends on camera/mode', & ' '//cSwitch(:iSwitch)//'help' ,'command line syntax', & ' '//cSwitch(:iSwitch)//'savefrm' ,'save processed frames only', & ' '//cSwitch(:iSwitch)//'test' ,'test mode (no frames written)', & ' '//cSwitch(:iSwitch)//'parameters' ,'print values of all free parameters', & ' '//cSwitch(:iSwitch)//'version' ,'print version number', & ' '//cSwitch(:iSwitch)//'c3mask=<-1,0,1>' ,'"bad-pixel" mask: auto/off/on', & ' '//cSwitch(:iSwitch)//'nofudge' ,'skip fudge factors for dark current', & ' '//cSwitch(:iSwitch)//'near_pattern' ,'select nearest calibration pattern' call Say(cSay,'S','Stop','end syntax') end if call ForeignI4Arg(cArg, 'quiet', 0, iQuiet) if (ForeignArgSet(cArg,'version')) ! Print version number & call Say(cSay,'S','Stop','Version '//cDbl2Str(version,2)// & 'Paul Hick (UCSD/CASS; pphick@ucsd.edu)') if (ForeignArgSet(cArg,'parameters')) then ! Print all free parameters write (*,'(1X,A23)') 'ped_mean_excess' write (*,'(3X,3A7)' ) & 'c1', 'c2', 'c3' write (*,'(1X,A1,I1,3F7.1)') & ('m', j ,(ped_mean_excess_cam(i,j),i=1,3),j=0,2) C write (*,'(1X,A23,7X,A23)') 'ped_mean_ref','ped_mean_excess' C write (*,'(3X,3A7,9X,3A7)' ) C & 'c1', 'c2', 'c3', 'c1', 'c2', 'c3' C write (*,'(1X,A1,I1,3F7.1,7X,A1,I1,3F7.1)') C & ('m', j ,(ped_mean_ref_cam (i,j),i=1,3), C & 'm', j ,(ped_mean_excess_cam(i,j),i=1,3),j=0,2) write (*,*) ' ' write (*,'(1X,A23,7X,A23)') 'ped_median_deficit','ped_median_excess' write (*,'(3X,3A7,9X,3A7)' ) & 'c1', 'c2', 'c3', 'c1', 'c2', 'c3' write (*,'(1X,A1,I1,3F7.1,7X,A1,I1,3F7.1)') & ('m', j ,(ped_median_deficit_cam(i,j),i=1,3), & 'm', j ,(ped_median_excess_cam (i,j),i=1,3),j=0,2) C write (*,*) ' ' C write (*,'(1X,A23,7X,A23)') 'ped_init_cam' C write (*,'(3X,3A7,5X,3A7,5X,3A7)') C & 'c1', 'c2', 'c3', 'c1', 'c2', 'c3', 'c1', 'c2', 'c3' C write (*,'(1X,A1,I1,3F7.1,3X,A1,I1,3F7.1,3X,A1,I1,3F7.1)') C & ('m', j ,(ped_init_cam(0,i,j),i=1,3), C & 'm', j ,(ped_init_cam(1,i,j),i=1,3), C & 'm', j ,(ped_init_cam(2,i,j),i=1,3),j=0,2) write (*,*) ' ' write (*,'(1X,A23,7X,A23)') 'dark_cut','dark_power' write (*,'(3X,3A7,9X,3A7)' ) 'c1', 'c2', 'c3', 'c1', 'c2', 'c3' write (*,'(1X,A1,I1,3F7.1,7X,A1,I1,3F7.1)') & ('m', j ,(dark_cut_cam (i,j),i=1,3), & 'm', j ,(dark_power_cam(i,j),i=1,3),j=0,2) write (*,*) ' ' write (*,'(1X,A23)') 'ped_right_offset' write (*,'(3X,3A7)') 'c1', 'c2', 'c3' write (*,'(1X,A1,I1,3F7.1)') & ('m', j ,(ped_right_offset_cam(i,j),i=1,3),j=0,2) write (*,*) ' ' write (*,'(1X,A23,7X,A23)') 'nped_min','ndark_min' write (*,'(3X,3A7,9X,3A7)' ) 'c1', 'c2', 'c3', 'c1', 'c2', 'c3' write (*,'(1X,A1,I1,3I7,7X,A1,I1,3I7)') & ('m', j ,(nped_min_cam (i,j),i=1,3), & 'm', j ,(ndark_min_cam(i,j),i=1,3),j=0,2) write (*,*) ' ' write (*,'(1X,A23,7X,A23)') 'dark_scale_ratio','dark_max_ratio' write (*,'(3X,3A7,9X,3A7)' ) 'c1', 'c2', 'c3', 'c1', 'c2', 'c3' write (*,'(1X,A1,I1,3F7.1,7X,A1,I1,3F7.1)') & ('m', j ,(dark_scale_ratio_cam(i,j),i=1,3), & 'm', j ,(dark_max_ratio_cam (i,j),i=1,3),j=0,2) write (*,*) ' ' write (*,'(1X,A23,7X,A23)') 'dark_offset_fudge','dark_scale_fudge' write (*,'(3X,3A7,9X,3A7)' ) 'c1', 'c2', 'c3', 'c1', 'c2', 'c3' write (*,'(1X,A1,I1,3F7.1,7X,A1,I1,3F7.1)') & ('m', j ,(dark_offset_fudge_cam(i,j),i=1,3), & 'm', j ,(dark_scale_fudge_cam (i,j),i=1,3),j=0,2) call Say(cSay,'S','Stop','end parameters') end if bDigBase = ForeignArgSet(cArg, 'digsource' ) bDigDest = ForeignArgSet(cArg, 'digdestination') bNoDups = ForeignArgSet(cArg, 'nodups' ) bForce = ForeignArgSet(cArg, 'force' ) bSaveFrm = ForeignArgSet(cArg, 'savefrm' ) bTest = ForeignArgSet(cArg, 'test' ) bNoFudge = ForeignArgSet(cArg, 'nofudge' ) ipick = 0 if (ForeignArgSet(cArg, 'near_pattern')) ipick = 1 !========================== ibad = BadI4() rbad = BadR4() ! Replace BadInitFlag by rbad call ArrR4Mask(9,ped_median_deficit_cam,BadInitFlag,rbad,0.0,0.0,1.0,ped_median_deficit_cam) call ArrR4Mask(9,ped_right_offset_cam ,BadInitFlag,rbad,0.0,0.0,1.0,ped_right_offset_cam ) call ArrR4Mask(9,dark_scale_ratio_cam ,BadInitFlag,rbad,0.0,0.0,1.0,dark_scale_ratio_cam ) call ArrR4Mask(9,dark_scale_fudge_cam ,BadInitFlag,rbad,0.0,0.0,1.0,dark_scale_fudge_cam ) !========================== if (cCalPattern .ne. ' ') then ! Cal pattern specified on cmd line ! Cam 3 only: specify whether to use regular calibration pattern (0) or ! pattern with "bad-pixel" mask applied (1). if (icam .eq. 3) then call ForeignI4Arg(cArg, 'c3mask', 0, c3mask_ref) if (c3mask_ref .ne. 1) c3mask_ref = 0 end if call smei_cal_get(cCalPattern,0,icam,mode,c3mask_ref,TFrm,cal_name,cal_pattern,pattern_dark,bNewCal,cal_version) base_aux(SMEI__BAS_PATTERN_DARK) = pattern_dark i = iSetFileSpec(cal_name) i = iGetFileSpec(FIL__NAME,FIL__NAME,cal_name) else ! Cam 3 only: specify whether to use regular calibration pattern (0), ! pattern with "bad-pixel" mask applied (1), or decide automatically ! on frame-by-frame basis (-1). ! For cam 3, mode 0 data, c3mask=1 is useful only for pseudo-frames ! created by smeidb_cal (mode 0 frames rebinned to mode 1). if (icam .eq. 3) then call ForeignI4Arg(cArg, 'c3mask', -1, c3mask_ref) if (mode .eq. 1) then if (c3mask_ref .ne. 1 .and. c3mask_ref .ne. 0) c3mask_ref = -1 else if (c3mask_ref .ne. 1) c3mask_ref = 0 end if end if if (TBeg(1) .eq. iBad .or. TEnd(1) .eq. iBad) call Say(cSay,'E','times', & 'specify begin and end time'// & '#type "'//cSay//' '//cSwitch(:iSwitch)//'help" for more info') pattern_dark = rbad end if if (icam .eq. iBad) call Say(cSay,'E','camera', & 'specify camera'// & '#type "'//cSay//' '//cSwitch(:iSwitch)//'help" for more info') i = 0 i = i+Time2Str(SMEI__UT_FORMAT,TBeg ,cStr(i+1:))+1 i = i+Str2Str ('-' ,cStr(i+1:))+1 i = i+Time2Str(SMEI__UT_FORMAT,TEnd ,cStr(i+1:)) i = i+Str2Str (', camera' ,cStr(i+1:))+1 i = i+Int2Str (icam ,cStr(i+1:)) i = i+Str2Str (', mode' ,cStr(i+1:))+1 i = i+Int2Str (mode ,cStr(i+1:)) i = i+Str2Str('#read from '//cHideLogical(cBase),cStr(i+1:)) if (bDigBase) i = i+Str2Str(' (digging)' ,cStr(i+1:)) if (bTest) then i = i+Str2Str('#test mode: no frames written',cStr(i+1:)) else i = i+Str2Str('#write to '//cHideLogical(cDest),cStr(i+1:)) end if if (bDigDest) i = i+Str2Str(' (digging)' ,cStr(i+1:)) if (bSaveFrm) then i = i+Str2Str('#saving processed frames only',cStr(i+1:)) else if (bNoDups) i = i+Str2Str('#replace existing copies' ,cStr(i+1:)) if (bForce ) i = i+Str2Str('#force new base calculation',cStr(i+1:)) end if if (bNoFudge) i = i+Str2Str('#no dark current fudging',cStr(i+1:)) if (icam .eq. 3 .and. mode .eq. 1) then if (c3mask_ref .eq. -1) then i = i+Str2Str('#automatic calibration pattern selection',cStr(i+1:)) else if (c3mask_ref .eq. 0) then i = i+Str2Str('#force regular calibration pattern',cStr(i+1:)) else if (c3mask_ref .eq. 1) then i = i+Str2Str('#force "bad-pixel" calibration pattern',cStr(i+1:)) end if end if if (ipick .ne. 0) i = i+Str2Str('#always use nearest cal pattern',cStr(i+1:)) call Say(cSay,'I','do',cStr) nFrm = smei_frm_getlist(TBeg,TEnd,icam,mode,bDigBase,cBase,cList) if (nFrm .eq. 0) call Say(cSay,'E','do','no frames found in this time range') i = 0 i = i+Int2Str(nFrm ,cStr(i+1:))+1 i = i+Str2Str('frames',cStr(i+1:)) call Say(cSay,'I','do',cStr) call smei_cal_init() !ped_aux(SMEI__PED_MEAN_REF ) = ped_mean_ref_cam (icam,mode) ped_aux(SMEI__PED_MEAN_EXCESS ) = ped_mean_excess_cam (icam,mode) call ForeignR4Arg(cArg,'ped_median_deficit',-ped_median_deficit_cam(icam,mode),ped_aux(SMEI__PED_MEDIAN_DEFICIT)) ped_aux(SMEI__PED_MEDIAN_DEFICIT) = -ped_aux(SMEI__PED_MEDIAN_DEFICIT) call ForeignR4Arg(cArg,'ped_median_excess' ,ped_median_excess_cam (icam,mode),ped_aux(SMEI__PED_MEDIAN_EXCESS )) ! For camera 3 data after 2005, Doy 159 the next two parameters ! are adjusted later on. ped_aux(SMEI__PED_RIGHT_OFFSET ) = ped_right_offset_cam (icam,mode) ped_aux(SMEI__PED_LEFT_SIDE_ONLY) = rbad call ForeignR4Arg(cArg,'dark_cut',dark_cut_cam(icam,mode),dark_aux(SMEI__DRK_CUT)) !base_aux(SMEI__BAS_PED_INIT+0 ) = ped_init_cam(0,icam,mode) !base_aux(SMEI__BAS_PED_INIT+1 ) = ped_init_cam(1,icam,mode) !base_aux(SMEI__BAS_PED_INIT+2 ) = ped_init_cam(2,icam,mode) call ForeignR4Arg(cArg,'nped_min' ,float(nped_min_cam (icam,mode)), base_aux(SMEI__BAS_NPED_MIN )) call ForeignR4Arg(cArg,'ndark_min',float(ndark_min_cam(icam,mode)), base_aux(SMEI__BAS_NDARK_MIN)) base_aux(SMEI__BAS_PATTERN_DARK ) = pattern_dark ! The ratio to scale the pattern with uses a straight mean. base_aux(SMEI__BAS_DARK_GUESS ) = rbad call ForeignR4Arg(cArg,'dark_scale_ratio',dark_scale_ratio_cam(icam,mode),base_aux(SMEI__BAS_DARK_SCALE_RATIO)) call ForeignR4Arg(cArg,'dark_max_ratio',dark_max_ratio_cam(icam,mode),base_aux(SMEI__BAS_DARK_MAX_RATIO)) base_aux(SMEI__BAS_NPED_BAD) = -1 ! Intializes smei_frm_base dark_scale_fudge = dark_scale_fudge_cam (icam,mode) dark_offset_fudge = dark_offset_fudge_cam(icam,mode) nframe = 0 iFirst = 1 ! Process all frames stored in cList. iRecl = 0 if (.not. bOpenFile(OPN__TEXT+OPN__REOPEN,iU,cList,iRecl)) then i = iOSDeleteFile(cList) call Say(cSay,'E','#'//cList,'error opening list with frame names') end if i = iScratchLun(iU) ! Mark cList as scratch file read (iU,'(A)',iostat=i) cFile do while (i .eq. 0) nframe = nframe+1 bFrameOK = smei_frm_read(OPEN_+min(mod(nframe,PRINT_),1)*OPN__NOMESSAGE, & cFile, SMEI__FRM_NPIX, nx,ny,nb,ccd_frame,hdr,headroom) ! If smei_frm_getlist works as advertised, this never happens if (nint(hdr(SMEI__HDR_CAMERA)) .ne. icam) stop 'bad camera' if (nint(hdr(SMEI__HDR_MODE )) .ne. mode) stop 'bad mode' if (bFrameOK) then ! If the version number has increased, redo the base calculation if (version .gt. hdr(SMEI__HDR_VERSION)) hdr(SMEI__HDR_BASE_DONE) = 0 bFrameOK = bForce .or. bSaveFrm .or. nint(hdr(SMEI__HDR_BASE_DONE)) .eq. 0 if (.not. bFrameOK) call Say(cSay,'W','#'//cFile,'already done') end if if (bFrameOK) then call smei_Time2Split(1,cFile,TFrm) ! After 2005, Doy 159 (camera 3 only) the pedestal is ! calculated only from the underscan region (leftside of frame) ! Frames are processed chronologically, so the next if-then ! block starts to be executed when the indicated date is passed. if (icam .eq. 2) then i = Time2Str(SMEI__UT_FORMAT,TFrm,cStr) if (cStr(:i) .gt. SMEI__UT_C2_LEFT_SIDE_ONLY) then ped_aux(SMEI__PED_LEFT_SIDE_ONLY) = 1.0 end if else if (icam .eq. 3) then i = Time2Str(SMEI__UT_FORMAT,TFrm,cStr) if (cStr(:i) .gt. SMEI__UT_C3_LEFT_SIDE_ONLY) then ped_aux(SMEI__PED_RIGHT_OFFSET ) = rbad ped_aux(SMEI__PED_LEFT_SIDE_ONLY) = 1.0 end if end if if (cCalPattern .eq. ' ') then ! Read new calibration pattern if necessary if (icam .eq. 3) then if (c3mask_ref .eq. -1) then c3mask = smei_frm_c3mask_active(0,cStr(:i), & nint(hdr(SMEI__HDR_FLAT_ENABLED)),cStr,fudge) else c3mask = c3mask_ref end if else c3mask = 0 end if cStr = cal_name ! Name of current calibration pattern call smei_cal_get(' ',ipick,icam,mode,c3mask,TFrm,cal_name, & cal_pattern,pattern_dark,bNewCal,cal_version) i = iSetFileSpec(cal_name) ! New calibration pattern i = iGetFileSpec(FIL__NAME,FIL__NAME,cal_name) if (cStr .ne. cal_name .and. ped_aux(SMEI__PED_RIGHT_OFFSET) .ne. rbad) then ped_aux(SMEI__PED_RIGHT_OFFSET) = ped_right_offset_cam(icam,mode) call Say(cSay,'I','reset','left-right pedestal offset to zero') end if base_aux(SMEI__BAS_PATTERN_DARK) = pattern_dark else ! Cal pattern specified on cmd line if (icam .eq. 3) then ! Camera 3 ! Flatfield flag for frame and pattern must match bFrameOK = c3mask_ref .eq. nint(hdr(SMEI__HDR_FLAT_ENABLED)) if (.not. bFrameOK) then! Flatfield mismatch ! If this happens the ped/dark calculation below is ! skipped if (c3mask_ref .eq. 0) then call Say(cSay,'W','#'//cCalPattern,'flatfield "off"; frm "on"') else call Say(cSay,'W','#'//cCalPattern,'flatfield "on"; frm "off"') end if call Say(cSay,'W','#'//cFile,'flatfield mismatch') end if ! c3mask will be wrong if there is a flatfield mismatch, ! but then c3mask is not used, because the call to ! smei_frm_c3fudge is bypassed. c3mask = c3mask_ref else c3mask = 0 end if end if if (icam .eq. 3) then ! Check orbital pattern name ! The orbital pattern is created by smei_orb. ! smei_orb has a command line argument -update that forces ! an update of the name of the orbital in the frame headers. ! Setting the orbital pattern name here avoids having to ! use -update, resulting in a considerable speedup of smei_orb. ! The disadvantage is that the name constructed here ! MUST be consistent with the names created by smei_orb. if (mode .eq. 1) then ! Construct orb_name from start time of orbit call smei_orb_mkname(TFrm,orb_name) else call smei_hdr_str(SMEI__HDR_ORB_PATTERN,hdr,orb_name) if (itrim(orb_name) .ne. 0) ! Orbital pattern already set & call Say(cSay,'W','hdr','erase orbital pattern (not mode 1)') orb_name = ' ' end if end if call smei_frm_clean(nx,ny,ccd_frame,ccd_frame) dark_aux(SMEI__DRK_HEADROOM) = headroom if (bFrameOK) then ! Pedestal and dark current calculation bFrameOK = smei_frm_base(cFile, nx, ny, ccd_frame, cal_pattern, & nped, ped_val, ndark, dark_val, dark_pixs, ped_aux, dark_aux, base_aux) else ! This happens only for cam 3 if the pattern (and c3mask_ref) ! were specified on the cmd line, and there is a mismatch ! between the flatfield flag on frame and pattern. nped = 0 ndark = 0 ped_val = rbad dark_val = rbad end if hdr(SMEI__HDR_BASE_DONE) = 1 ! Base subtraction tried ! Apply the fudge correction to the dark current (even if the ! base calculation is rejected). if (.not. BadR4Test(dark_val)) then ! There is a dark current value fudge = dark_scale_fudge if (BadR4Test(fudge)) fudge = smei_frm_c3fudge(TFrm,c3mask) if (.not. bNoFudge) dark_val = fudge*dark_val+dark_offset_fudge end if if (bFrameOK) then dark_ratio = dark_val/pattern_dark power = dark_power_cam(icam,mode) ! Raise to power 1.5 for cam 1 and 2 if (power .ne. 1.0) dark_ratio = dark_ratio**power n_ante = 0 +nx*ny ! Position of frame before measle removal n_post = n_ante+nx*ny ! Position of frame after measle removal ! The pattern also has non-zero values in the pedestal columns. ! We only want to subtract the pattern in the dark current and uncovered pixels ! (or at least not scale with the dark current). do n=1,nx*ny if (BadR4Test(ccd_frame(n))) then ccd_frame(n_ante+n) = ccd_frame(n) else ccd_frame(n_ante+n) = ccd_frame(n)-ped_val i = mod(n-1,nx)+1 if (SMEI__DRK_COL(1) .le. i .and. i .le. SMEI__DRK_COL(SMEI__DRK_NCOL)) & ccd_frame(n_ante+n) = ccd_frame(n_ante+n)-dark_ratio*cal_pattern(n) end if ccd_frame(n_post+n) = ccd_frame(n_ante+n) end do hdr(SMEI__HDR_BASE_OK) = 1 ! Base subtraction OK call smei_frm_measle(.TRUE.,icam,nx,ny,ccd_frame(n_post+1),npos_measles,nbig_measles,mask,indx) call smei_frm_stats(nx,ny,ccd_frame(n_post+1),pixsum,ipixsum,pixdif,ipixdif) hdr(SMEI__HDR_N_POS_MEASLES) = npos_measles hdr(SMEI__HDR_N_BIG_MEASLES) = nbig_measles hdr(SMEI__HDR_PIXSUM ) = pixsum hdr(SMEI__HDR_N_PIXSUM) = ipixsum hdr(SMEI__HDR_PIXDIF ) = pixdif hdr(SMEI__HDR_N_PIXDIF) = ipixdif ! sqre_pixs and cntr_pixs are equivalenced, so don't change order of statements n = smei_frm_pickup(2,nx,ny,ccd_frame(n_ante+1),sqre_pixs) hdr(SMEI__HDR_SQUARE) = ArrR4Mean(-n,sqre_pixs,m) n = smei_frm_pickup(3,nx,ny,ccd_frame(n_ante+1),cntr_pixs) hdr(SMEI__HDR_CENTER) = ArrR4Mean(-n,cntr_pixs,m) else hdr(SMEI__HDR_BASE_OK) = 0 ! Base subtraction failed hdr(SMEI__HDR_N_POS_MEASLES) = hdr(SMEI__HDR_BAD_DATA) hdr(SMEI__HDR_N_BIG_MEASLES) = hdr(SMEI__HDR_BAD_DATA) hdr(SMEI__HDR_PIXSUM ) = hdr(SMEI__HDR_BAD_DATA) hdr(SMEI__HDR_N_PIXSUM) = hdr(SMEI__HDR_BAD_DATA) hdr(SMEI__HDR_PIXDIF ) = hdr(SMEI__HDR_BAD_DATA) hdr(SMEI__HDR_N_PIXDIF) = hdr(SMEI__HDR_BAD_DATA) hdr(SMEI__HDR_SQUARE) = hdr(SMEI__HDR_BAD_DATA) hdr(SMEI__HDR_CENTER) = hdr(SMEI__HDR_BAD_DATA) end if call smei_hdr_str(-SMEI__HDR_CAL_PATTERN, hdr, cal_name) call smei_hdr_str(-SMEI__HDR_ORB_PATTERN, hdr, orb_name) if (BadR4Test(ped_val )) then hdr(SMEI__HDR_PEDESTAL) = hdr(SMEI__HDR_BAD_DATA) else hdr(SMEI__HDR_PEDESTAL) = ped_val end if if (BadR4Test(dark_val)) then hdr(SMEI__HDR_DARK_MEDIAN) = hdr(SMEI__HDR_BAD_DATA) hdr(SMEI__HDR_DARK_MEAN ) = hdr(SMEI__HDR_BAD_DATA) else hdr(SMEI__HDR_DARK_MEDIAN) = dark_val ! Fudged median hdr(SMEI__HDR_DARK_MEAN ) = dark_aux(SMEI__DRK_MEAN)! Unfudged mean end if if (BadR4Test(ped_aux (SMEI__PED_STDEV))) then hdr(SMEI__HDR_PED_SIGMA) = hdr(SMEI__HDR_BAD_DATA) else hdr(SMEI__HDR_PED_SIGMA) = ped_aux(SMEI__PED_STDEV) end if if (BadR4Test(dark_aux(SMEI__DRK_STDEV))) then hdr(SMEI__HDR_DARK_SIGMA) = hdr(SMEI__HDR_BAD_DATA) else hdr(SMEI__HDR_DARK_SIGMA) = dark_aux(SMEI__DRK_STDEV) end if hdr(SMEI__HDR_N_PEDESTAL ) = nped hdr(SMEI__HDR_N_DARK_CURRENT) = ndark hdr(SMEI__HDR_N_SATURATED ) = smei_frm_saturated(hdr,ccd_frame) !if (smei_frm_bad_quaternion(hdr)) then ! hdr(SMEI__HDR_BAD_QUAT) = 1 !else ! hdr(SMEI__HDR_BAD_QUAT) = 0 !endif ! Write Fits file ! NOTE: hdr(SMEI__HDR_SATURATED) was set in smei_frm_read hdr(SMEI__HDR_VERSION) = version ! Set current version number if (bSaveFrm) then call ArrR4AddConstant(-nx*ny,ccd_frame,-ped_val-dark_val,ccd_frame) i = iSetFileSpec(cFile) i = iGetFileSpec(FIL__NAME,FIL__NAME,cStr) i = iFilePath(cDest,0,' ',cStr,cStr) call ArrR4Mask(3*nx*ny,ccd_frame,rbad,0.0,0.0,0.0,1.0,ccd_frame) call smei_frm_fts(cStr(:i)//'_0.fts.gz', nx, ny, -4, ccd_frame( 1), hdr) call smei_frm_fts(cStr(:i)//'_1.fts.gz', nx, ny, -4, ccd_frame(n_ante+1), hdr) call smei_frm_fts(cStr(:i)//'_2.fts.gz', nx, ny, -4, ccd_frame(n_post+1), hdr) else if (.not. bTest) then call smei_frm_write(cFile,hdr,cDest,bDigDest,bNoDups) end if end if read (iU,'(A)',iostat=i) cFile end do i = iFreeLun(-iU) call Say(cSay,'S','Stop','Done') end