Loading calculate_error.pro 0 → 100644 +46 −0 Original line number Diff line number Diff line function calculate_error, cal_pack, err_info, detector if detector.toupper() eq 'UV' then channel = cal_pack.uv_channel if detector.toupper() eq 'VL' then channel = cal_pack.vl_channel radiometry = channel.radiometry if detector.toupper() eq 'VL' then begin unit_factor = radiometry.msb.value unit_error = radiometry.msb.error endif else begin unit_factor = 1. unit_error = 0. endelse data = err_info.data.image s = where(data eq 0., count) ; calculate the initial error from dark/bias subtraction (relative, squared) if detector.toupper() eq 'UV' then $ error = (data * radiometry.iaps_gain.value * radiometry.f_noise.value ^ 2 + $ 2. * err_info.dark.error ^ 2) / data ^ 2 if detector.toupper() eq 'VL' then $ error = (data * radiometry.gain.value + 2. * err_info.bias.error ^ 2 + $ 2. * err_info.dark.error ^ 2) / data ^ 2 ; add error for flat-field correction error += (err_info.flat.error / err_info.flat.image) ^ 2 ; add error for vignetting correction error += (err_info.vign.error / err_info.vign.image) ^ 2 ; add error for radiometric calibration error += (radiometry.rad_response.error / radiometry.rad_response.value) ^ 2 $ + (cal_pack.instrument.pupil_area.error / cal_pack.instrument.pupil_area.value) ^ 2 $ + (channel.angular_pixel.error / channel.angular_pixel.value) ^ 2 $ + (unit_error / unit_factor) ^ 2 if detector.toupper() eq 'UV' then $ error += (err_info.eff.error / err_info.eff.image) ^ 2 if count gt 0 then error[s] = 0. return, error end No newline at end of file calculate_error_polariz.pro 0 → 100644 +81 −0 Original line number Diff line number Diff line function calculate_error_polariz, cal_pack, err_info channel = cal_pack.vl_channel radiometry = channel.radiometry data = err_info.data.image error = data * 0. i_error = 0. q_error = 0. u_error = 0. ; error over stokes i for k = 0, 3 do $ ; initial error for dark/bias subtraction (absolute, squared) error[*, *, k] = data[*, *, k] * radiometry.gain.value $ + 2. * err_info.bias.error^2 $ + 2. * err_info.dark.error^2 for k = 0, 3 do $ ; error for demodulation (relative, squared) i_error += (data[*, *, k]^2 * err_info.demod.error[*, *, k, 0]^2 $ + err_info.demod.image[*, *, k, 0]^2 * error[*, *, k])/err_info.stokes.image[*, *, 0]^2 ; add error for flat-field correction i_error += (err_info.flat.error/err_info.flat.image)^2 ; add error for vignetting correction i_error += (err_info.vign.error/err_info.vign.image)^2 ; add error for radiometric calibration i_error += (radiometry.rad_response.error/radiometry.rad_response.value)^2 $ + (cal_pack.instrument.pupil_area.error/cal_pack.instrument.pupil_area.value)^2 $ + (channel.angular_pixel.error/channel.angular_pixel.value)^2 $ + (radiometry.msb.error/radiometry.msb.value)^2 ; error over stokes q for k = 0, 3 do $ ; initial error for dark/bias subtraction (absolute, squared) error[*, *, k] = data[*, *, k] * radiometry.gain.value $ + err_info.bias.error^2 $ + err_info.dark.error^2 for k = 0, 3 do $ ; error for demodulation (relative, squared) q_error += (data[*, *, k]^2 * err_info.demod.error[*, *, k, 1]^2 $ + err_info.demod.image[*, *, k, 1]^2 * error[*, *, k])/err_info.stokes.image[*, *, 1]^2 ; add error for flat-field correction q_error += (err_info.flat.error/err_info.flat.image)^2 ; add error for vignetting correction q_error += (err_info.vign.error/err_info.vign.image)^2 ; add error for radiometric calibration q_error += (radiometry.rad_response.error/radiometry.rad_response.value)^2 $ + (cal_pack.instrument.pupil_area.error/cal_pack.instrument.pupil_area.value)^2 $ + (channel.angular_pixel.error/channel.angular_pixel.value)^2 $ + (radiometry.msb.error/radiometry.msb.value)^2 ; error over stokes u for k = 0, 3 do $ ; initial error for dark/bias subtraction (absolute, squared) error[*, *, k] = data[*, *, k] * radiometry.gain.value $ + err_info.bias.error^2 $ + err_info.dark.error^2 for k = 0, 3 do $ ; error for demodulation (relative, squared) u_error += (data[*, *, k]^2 * err_info.demod.error[*, *, k, 2]^2 $ + err_info.demod.image[*, *, k, 2]^2 * error[*, *, k])/err_info.stokes.image[*, *, 2]^2 ; add error for flat-field correction u_error += (err_info.flat.error/err_info.flat.image)^2 ; add error for vignetting correction u_error += (err_info.vign.error/err_info.vign.image)^2 ; add error for radiometric calibration u_error += (radiometry.rad_response.error/radiometry.rad_response.value)^2 $ + (cal_pack.instrument.pupil_area.error/cal_pack.instrument.pupil_area.value)^2 $ + (channel.angular_pixel.error/channel.angular_pixel.value)^2 $ + (radiometry.msb.error/radiometry.msb.value)^2 return, [[[i_error]], [[q_error]], [[u_error]]] end No newline at end of file check_compression.pro 0 → 100644 +41 −0 Original line number Diff line number Diff line pro check_compression, image, header bin_type = fxpar(header, 'BIN_TYPE', missing = -1) if bin_type lt 0 then return nbin = fxpar(header, 'NBIN') bin_fact = 2^bin_type if nbin ne bin_fact^2 then begin naxis1 = fxpar(header, 'NAXIS1') naxis2 = fxpar(header, 'NAXIS2') naxis1 = naxis1/bin_fact naxis2 = naxis2/bin_fact image = rebin(image, naxis1, naxis2) image = long(image) * bin_fact^2 fxaddpar, header, 'NAXIS1', naxis1 fxaddpar, header, 'NAXIS2', naxis2 fxaddpar, header, 'NBIN1', bin_fact fxaddpar, header, 'NBIN2', bin_fact fxaddpar, header, 'NBIN', bin_fact^2 fxaddpar, header, 'PXEND1', naxis1 fxaddpar, header, 'PXEND2', naxis2 fxaddpar, header, 'FIRSTROW', 0, before = 'CHECKSUM' fxaddpar, header, 'B0_BIN', bin_type, before = 'CHECKSUM' fxaddpar, header, 'B0_DQ', 0, before = 'CHECKSUM' fxaddpar, header, 'B0_STOP', 2048, before = 'CHECKSUM' fxaddpar, header, 'B1_BIN', 0, before = 'CHECKSUM' fxaddpar, header, 'B1_DQ', 0, before = 'CHECKSUM' fxaddpar, header, 'B1_STOP', 0, before = 'CHECKSUM' fxaddpar, header, 'B2_BIN', 0, before = 'CHECKSUM' fxaddpar, header, 'B2_DQ', 0, before = 'CHECKSUM' fxaddpar, header, 'B2_STOP', 0, before = 'CHECKSUM' comment = fxpar(header, 'COMMENT') comment = ['Image was rebinned according to the commanded binning factor.', comment] sxdelpar, header, 'COMMENT' foreach line, comment do fxaddpar, header, 'COMMENT', line endif end No newline at end of file metis_dark_uvda.pro +155 −139 Original line number Diff line number Diff line function m_angle_uv, matrix a = mean(matrix[10:45, 10:45]) b = mean(matrix[980:1014, 10:45]) c = mean(matrix[10:45, 980:1014]) d = mean(matrix[980:1014, 980:1014]) a = mean(matrix[8 : 43, 8 : 43]) b = mean(matrix[980 : 1015, 8 : 43]) c = mean(matrix[8 : 43, 980 : 1015]) d = mean(matrix[980 : 1015, 980 : 1015]) return, [a, b, c, d] end function m_center_uv, matrix a = mean(matrix[454:574, 460:582]) a = mean(matrix[452 : 571, 460 : 579]) return, a end function normalize_dark, data, dark dark_bin = rebin(dark, 1024, 1024, /sample) data_bin = rebin(data, 1024, 1024, /sample) Loading @@ -30,8 +27,7 @@ function normalize_dark, data, dark return, dark_norm end function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history function metis_dark_uvda, data, header, cal_pack, err_info, quality_matrix = quality_matrix, history = history dark = cal_pack.uv_channel.dark dit = header.dit Loading @@ -43,9 +39,6 @@ function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix masking = header.masking tsensor = header.tsensor ; WARN - temporary patch to handle local l1 fits files if tsensor eq 0. then tsensor = -25. journal, 'Dark-current correction:' journal, ' dit = ' + string(dit, format = '(I0)') + ' ms' journal, ' ndit1 = ' + string(ndit1, format = '(I0)') Loading @@ -53,21 +46,21 @@ function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix journal, ' obj_cnt = ' + string(obj_cnt, format = '(I0)') journal, ' tsensor = ' + string(tsensor, format = '(F0)') + ' degC' flag_use_extrapolated = boolean(1) flag_use_extrapolated = boolean(0) flag_use_notfullset = boolean(1) flag_normalize_dark = boolean(0) flag_normalize_dark = boolean(1) obt_available = !null ; check if one or more dark matrices with same dit, ndit1, ndit2, and tsensor exist ; check if one or more dark matrices with same parameters exist for i = 0, n_elements(dark) - 1 do begin if $ (dark[i].extrapol eq 'False') and $ (dark[i].dit eq dit) and $ (dark[i].ndit1 eq ndit1) and $ (dark[i].ndit2 eq ndit2) and $ (abs(dark[i].tsensor - tsensor) lt 5) and $ (flag_use_notfullset or dark[i].fullset) and $ (dark[i].extrapol eq 'False') then $ (flag_use_notfullset or dark[i].fullset) then $ obt_available = [obt_available, dark[i].obt_beg] endfor Loading @@ -75,12 +68,25 @@ function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix if ~isa(obt_available) and flag_use_extrapolated then begin for i = 0, n_elements(dark) - 1 do begin if $ (dark[i].extrapol eq 'True') and $ (dark[i].dit eq dit) and $ (dark[i].ndit1 eq ndit1) and $ (dark[i].ndit2 eq ndit2) and $ (abs(dark[i].tsensor - tsensor) lt 5) and $ (flag_use_notfullset or dark[i].fullset) and $ (dark[i].extrapol eq 'True') then $ (flag_use_notfullset or dark[i].fullset) then $ obt_available = [obt_available, dark[i].obt_beg] endfor endif ; if not, search for dark matrices with no constraint on tsensor if ~isa(obt_available) then begin for i = 0, n_elements(dark) - 1 do begin if $ (dark[i].extrapol eq 'False') and $ (dark[i].dit eq dit) and $ (dark[i].ndit1 eq ndit1) and $ (dark[i].ndit2 eq ndit2) and $ (flag_use_notfullset or dark[i].fullset) then $ obt_available = [obt_available, dark[i].obt_beg] endfor endif Loading @@ -93,7 +99,7 @@ function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix endif ; if yes, chose the temporally closest delta_obt = min(abs(obt_beg - obt_available), j) delta = min(abs(obt_beg - obt_available), j) i = where(dark.obt_beg eq obt_available[j]) transient = where(dark[i].filename.cnt_1 ne '', n) Loading @@ -119,24 +125,34 @@ function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix ; set values in masked areas equal to zero if count gt 0 then data[mask] = 0. ; save data and dark image in the error info structure err_info.data.image = data err_info.dark.image = dark_image ; read the dark error matrix dark_error_file = dark[i].error if file_test(dark_error_file) then $ dark_error = float(readfits(cal_pack.path + dark_error_file, /silent)) else $ if file_test(cal_pack.path + dark_error_file) then begin dark_error = float(readfits(cal_pack.path + dark_error_file, hdr, /silent)) dark_syst_err = fxpar(hdr, 'SYSTEMAT', missing = 0.) endif else begin dark_error = fltarr(1024, 1024) dark_syst_err = 0. endelse ; rebin it and take into account the exposure time correction dark_error = sqrt(rebin(dark_error^2, header.naxis1, header.naxis2) * nbin)/nbin ; NOTE - some systematic errors reported in the header of the dark fits files are wrong ; calculate the total uncertainty if not isa(error) then error = 0. error += ((data > 0)/cal_pack.uv_channel.radiometry.iaps_gain.value * cal_pack.uv_channel.radiometry.f_noise.value + 2. * dark_error^2)/data^2 if isa(dark_syst_err, /string) then dark_syst_err = 0. ; rebin the error matrix and take into account the exposure time correction of l1 images dark_error = sqrt(rebin(dark_error ^ 2, header.naxis1, header.naxis2) * nbin) / nbin dark_error = dark_error * ndit1 * ndit2 ; set error values in masked areas equal to zero if count gt 0 then error[mask] = 0. ; save errors in the error info structure err_info.dark.error = dark_error err_info.dark.syst_err = dark_syst_err jump: s = where(data le 0., count) s = where(data le 0, count) if isa(quality_matrix) and count gt 0 then quality_matrix[s] = 0 if ~isa(history) then history = !null Loading metis_dark_uvda_new.pro +154 −150 Original line number Diff line number Diff line function m_angle_uv, matrix a = mean(matrix[8 : 43, 8 : 43]) b = mean(matrix[980 : 1015, 8 : 43]) c = mean(matrix[8 : 43, 980 : 1015]) Loading @@ -9,14 +8,12 @@ function m_angle_uv, matrix end function m_center_uv, matrix a = mean(matrix[452 : 571, 460 : 579]) return, a end function correct_dark, data, dark, offset dark_bin = rebin(dark, 1024, 1024, /sample) data_bin = rebin(data, 1024, 1024, /sample) Loading @@ -30,8 +27,7 @@ function correct_dark, data, dark, offset return, dark_corr end function metis_dark_uvda_new, data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history function metis_dark_uvda_new, data, header, cal_pack, err_info, quality_matrix = quality_matrix, history = history dark = cal_pack.uv_channel.dark dit = header.dit Loading Loading @@ -126,26 +122,34 @@ function metis_dark_uvda_new, data, header, cal_pack, error = error, quality_mat ; set values in masked areas equal to zero if count gt 0 then data[mask] = 0. ; save data and dark image in the error info structure err_info.data.image = data err_info.dark.image = dark_image ; read the dark error matrix dark_error_file = dark[i].error if file_test(cal_pack.path + dark_error_file) then $ dark_error = float(readfits(cal_pack.path + dark_error_file, hdr, /silent)) else $ if file_test(cal_pack.path + dark_error_file) then begin dark_error = float(readfits(cal_pack.path + dark_error_file, hdr, /silent)) dark_syst_err = fxpar(hdr, 'SYSTEMAT', missing = 0.) endif else begin dark_error = fltarr(1024, 1024) dark_syst_err = 0. endelse ; NOTE - some systematic errors reported in the header of the dark fits files are wrong if isa(dark_syst_err, /string) then dark_syst_err = 0. ; rebin the error matrix and take into account the exposure time correction of l1 images dark_error = sqrt(rebin(dark_error ^ 2, header.naxis1, header.naxis2) * nbin) / nbin dark_error = dark_error * ndit1 * ndit2 if not isa(error) then error = 0. error += ((data > 0)/cal_pack.uv_channel.radiometry.iaps_gain.value * cal_pack.uv_channel.radiometry.f_noise.value + 2. * dark_error^2)/data^2 ; set error values in masked areas equal to zero if count gt 0 then error[mask] = 0. ; save errors in the error info structure err_info.dark.error = dark_error err_info.dark.syst_err = dark_syst_err jump: s = where(data le 0., count) s = where(data le 0, count) if isa(quality_matrix) and count gt 0 then quality_matrix[s] = 0 if ~isa(history) then history = !null Loading Loading
calculate_error.pro 0 → 100644 +46 −0 Original line number Diff line number Diff line function calculate_error, cal_pack, err_info, detector if detector.toupper() eq 'UV' then channel = cal_pack.uv_channel if detector.toupper() eq 'VL' then channel = cal_pack.vl_channel radiometry = channel.radiometry if detector.toupper() eq 'VL' then begin unit_factor = radiometry.msb.value unit_error = radiometry.msb.error endif else begin unit_factor = 1. unit_error = 0. endelse data = err_info.data.image s = where(data eq 0., count) ; calculate the initial error from dark/bias subtraction (relative, squared) if detector.toupper() eq 'UV' then $ error = (data * radiometry.iaps_gain.value * radiometry.f_noise.value ^ 2 + $ 2. * err_info.dark.error ^ 2) / data ^ 2 if detector.toupper() eq 'VL' then $ error = (data * radiometry.gain.value + 2. * err_info.bias.error ^ 2 + $ 2. * err_info.dark.error ^ 2) / data ^ 2 ; add error for flat-field correction error += (err_info.flat.error / err_info.flat.image) ^ 2 ; add error for vignetting correction error += (err_info.vign.error / err_info.vign.image) ^ 2 ; add error for radiometric calibration error += (radiometry.rad_response.error / radiometry.rad_response.value) ^ 2 $ + (cal_pack.instrument.pupil_area.error / cal_pack.instrument.pupil_area.value) ^ 2 $ + (channel.angular_pixel.error / channel.angular_pixel.value) ^ 2 $ + (unit_error / unit_factor) ^ 2 if detector.toupper() eq 'UV' then $ error += (err_info.eff.error / err_info.eff.image) ^ 2 if count gt 0 then error[s] = 0. return, error end No newline at end of file
calculate_error_polariz.pro 0 → 100644 +81 −0 Original line number Diff line number Diff line function calculate_error_polariz, cal_pack, err_info channel = cal_pack.vl_channel radiometry = channel.radiometry data = err_info.data.image error = data * 0. i_error = 0. q_error = 0. u_error = 0. ; error over stokes i for k = 0, 3 do $ ; initial error for dark/bias subtraction (absolute, squared) error[*, *, k] = data[*, *, k] * radiometry.gain.value $ + 2. * err_info.bias.error^2 $ + 2. * err_info.dark.error^2 for k = 0, 3 do $ ; error for demodulation (relative, squared) i_error += (data[*, *, k]^2 * err_info.demod.error[*, *, k, 0]^2 $ + err_info.demod.image[*, *, k, 0]^2 * error[*, *, k])/err_info.stokes.image[*, *, 0]^2 ; add error for flat-field correction i_error += (err_info.flat.error/err_info.flat.image)^2 ; add error for vignetting correction i_error += (err_info.vign.error/err_info.vign.image)^2 ; add error for radiometric calibration i_error += (radiometry.rad_response.error/radiometry.rad_response.value)^2 $ + (cal_pack.instrument.pupil_area.error/cal_pack.instrument.pupil_area.value)^2 $ + (channel.angular_pixel.error/channel.angular_pixel.value)^2 $ + (radiometry.msb.error/radiometry.msb.value)^2 ; error over stokes q for k = 0, 3 do $ ; initial error for dark/bias subtraction (absolute, squared) error[*, *, k] = data[*, *, k] * radiometry.gain.value $ + err_info.bias.error^2 $ + err_info.dark.error^2 for k = 0, 3 do $ ; error for demodulation (relative, squared) q_error += (data[*, *, k]^2 * err_info.demod.error[*, *, k, 1]^2 $ + err_info.demod.image[*, *, k, 1]^2 * error[*, *, k])/err_info.stokes.image[*, *, 1]^2 ; add error for flat-field correction q_error += (err_info.flat.error/err_info.flat.image)^2 ; add error for vignetting correction q_error += (err_info.vign.error/err_info.vign.image)^2 ; add error for radiometric calibration q_error += (radiometry.rad_response.error/radiometry.rad_response.value)^2 $ + (cal_pack.instrument.pupil_area.error/cal_pack.instrument.pupil_area.value)^2 $ + (channel.angular_pixel.error/channel.angular_pixel.value)^2 $ + (radiometry.msb.error/radiometry.msb.value)^2 ; error over stokes u for k = 0, 3 do $ ; initial error for dark/bias subtraction (absolute, squared) error[*, *, k] = data[*, *, k] * radiometry.gain.value $ + err_info.bias.error^2 $ + err_info.dark.error^2 for k = 0, 3 do $ ; error for demodulation (relative, squared) u_error += (data[*, *, k]^2 * err_info.demod.error[*, *, k, 2]^2 $ + err_info.demod.image[*, *, k, 2]^2 * error[*, *, k])/err_info.stokes.image[*, *, 2]^2 ; add error for flat-field correction u_error += (err_info.flat.error/err_info.flat.image)^2 ; add error for vignetting correction u_error += (err_info.vign.error/err_info.vign.image)^2 ; add error for radiometric calibration u_error += (radiometry.rad_response.error/radiometry.rad_response.value)^2 $ + (cal_pack.instrument.pupil_area.error/cal_pack.instrument.pupil_area.value)^2 $ + (channel.angular_pixel.error/channel.angular_pixel.value)^2 $ + (radiometry.msb.error/radiometry.msb.value)^2 return, [[[i_error]], [[q_error]], [[u_error]]] end No newline at end of file
check_compression.pro 0 → 100644 +41 −0 Original line number Diff line number Diff line pro check_compression, image, header bin_type = fxpar(header, 'BIN_TYPE', missing = -1) if bin_type lt 0 then return nbin = fxpar(header, 'NBIN') bin_fact = 2^bin_type if nbin ne bin_fact^2 then begin naxis1 = fxpar(header, 'NAXIS1') naxis2 = fxpar(header, 'NAXIS2') naxis1 = naxis1/bin_fact naxis2 = naxis2/bin_fact image = rebin(image, naxis1, naxis2) image = long(image) * bin_fact^2 fxaddpar, header, 'NAXIS1', naxis1 fxaddpar, header, 'NAXIS2', naxis2 fxaddpar, header, 'NBIN1', bin_fact fxaddpar, header, 'NBIN2', bin_fact fxaddpar, header, 'NBIN', bin_fact^2 fxaddpar, header, 'PXEND1', naxis1 fxaddpar, header, 'PXEND2', naxis2 fxaddpar, header, 'FIRSTROW', 0, before = 'CHECKSUM' fxaddpar, header, 'B0_BIN', bin_type, before = 'CHECKSUM' fxaddpar, header, 'B0_DQ', 0, before = 'CHECKSUM' fxaddpar, header, 'B0_STOP', 2048, before = 'CHECKSUM' fxaddpar, header, 'B1_BIN', 0, before = 'CHECKSUM' fxaddpar, header, 'B1_DQ', 0, before = 'CHECKSUM' fxaddpar, header, 'B1_STOP', 0, before = 'CHECKSUM' fxaddpar, header, 'B2_BIN', 0, before = 'CHECKSUM' fxaddpar, header, 'B2_DQ', 0, before = 'CHECKSUM' fxaddpar, header, 'B2_STOP', 0, before = 'CHECKSUM' comment = fxpar(header, 'COMMENT') comment = ['Image was rebinned according to the commanded binning factor.', comment] sxdelpar, header, 'COMMENT' foreach line, comment do fxaddpar, header, 'COMMENT', line endif end No newline at end of file
metis_dark_uvda.pro +155 −139 Original line number Diff line number Diff line function m_angle_uv, matrix a = mean(matrix[10:45, 10:45]) b = mean(matrix[980:1014, 10:45]) c = mean(matrix[10:45, 980:1014]) d = mean(matrix[980:1014, 980:1014]) a = mean(matrix[8 : 43, 8 : 43]) b = mean(matrix[980 : 1015, 8 : 43]) c = mean(matrix[8 : 43, 980 : 1015]) d = mean(matrix[980 : 1015, 980 : 1015]) return, [a, b, c, d] end function m_center_uv, matrix a = mean(matrix[454:574, 460:582]) a = mean(matrix[452 : 571, 460 : 579]) return, a end function normalize_dark, data, dark dark_bin = rebin(dark, 1024, 1024, /sample) data_bin = rebin(data, 1024, 1024, /sample) Loading @@ -30,8 +27,7 @@ function normalize_dark, data, dark return, dark_norm end function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history function metis_dark_uvda, data, header, cal_pack, err_info, quality_matrix = quality_matrix, history = history dark = cal_pack.uv_channel.dark dit = header.dit Loading @@ -43,9 +39,6 @@ function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix masking = header.masking tsensor = header.tsensor ; WARN - temporary patch to handle local l1 fits files if tsensor eq 0. then tsensor = -25. journal, 'Dark-current correction:' journal, ' dit = ' + string(dit, format = '(I0)') + ' ms' journal, ' ndit1 = ' + string(ndit1, format = '(I0)') Loading @@ -53,21 +46,21 @@ function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix journal, ' obj_cnt = ' + string(obj_cnt, format = '(I0)') journal, ' tsensor = ' + string(tsensor, format = '(F0)') + ' degC' flag_use_extrapolated = boolean(1) flag_use_extrapolated = boolean(0) flag_use_notfullset = boolean(1) flag_normalize_dark = boolean(0) flag_normalize_dark = boolean(1) obt_available = !null ; check if one or more dark matrices with same dit, ndit1, ndit2, and tsensor exist ; check if one or more dark matrices with same parameters exist for i = 0, n_elements(dark) - 1 do begin if $ (dark[i].extrapol eq 'False') and $ (dark[i].dit eq dit) and $ (dark[i].ndit1 eq ndit1) and $ (dark[i].ndit2 eq ndit2) and $ (abs(dark[i].tsensor - tsensor) lt 5) and $ (flag_use_notfullset or dark[i].fullset) and $ (dark[i].extrapol eq 'False') then $ (flag_use_notfullset or dark[i].fullset) then $ obt_available = [obt_available, dark[i].obt_beg] endfor Loading @@ -75,12 +68,25 @@ function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix if ~isa(obt_available) and flag_use_extrapolated then begin for i = 0, n_elements(dark) - 1 do begin if $ (dark[i].extrapol eq 'True') and $ (dark[i].dit eq dit) and $ (dark[i].ndit1 eq ndit1) and $ (dark[i].ndit2 eq ndit2) and $ (abs(dark[i].tsensor - tsensor) lt 5) and $ (flag_use_notfullset or dark[i].fullset) and $ (dark[i].extrapol eq 'True') then $ (flag_use_notfullset or dark[i].fullset) then $ obt_available = [obt_available, dark[i].obt_beg] endfor endif ; if not, search for dark matrices with no constraint on tsensor if ~isa(obt_available) then begin for i = 0, n_elements(dark) - 1 do begin if $ (dark[i].extrapol eq 'False') and $ (dark[i].dit eq dit) and $ (dark[i].ndit1 eq ndit1) and $ (dark[i].ndit2 eq ndit2) and $ (flag_use_notfullset or dark[i].fullset) then $ obt_available = [obt_available, dark[i].obt_beg] endfor endif Loading @@ -93,7 +99,7 @@ function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix endif ; if yes, chose the temporally closest delta_obt = min(abs(obt_beg - obt_available), j) delta = min(abs(obt_beg - obt_available), j) i = where(dark.obt_beg eq obt_available[j]) transient = where(dark[i].filename.cnt_1 ne '', n) Loading @@ -119,24 +125,34 @@ function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix ; set values in masked areas equal to zero if count gt 0 then data[mask] = 0. ; save data and dark image in the error info structure err_info.data.image = data err_info.dark.image = dark_image ; read the dark error matrix dark_error_file = dark[i].error if file_test(dark_error_file) then $ dark_error = float(readfits(cal_pack.path + dark_error_file, /silent)) else $ if file_test(cal_pack.path + dark_error_file) then begin dark_error = float(readfits(cal_pack.path + dark_error_file, hdr, /silent)) dark_syst_err = fxpar(hdr, 'SYSTEMAT', missing = 0.) endif else begin dark_error = fltarr(1024, 1024) dark_syst_err = 0. endelse ; rebin it and take into account the exposure time correction dark_error = sqrt(rebin(dark_error^2, header.naxis1, header.naxis2) * nbin)/nbin ; NOTE - some systematic errors reported in the header of the dark fits files are wrong ; calculate the total uncertainty if not isa(error) then error = 0. error += ((data > 0)/cal_pack.uv_channel.radiometry.iaps_gain.value * cal_pack.uv_channel.radiometry.f_noise.value + 2. * dark_error^2)/data^2 if isa(dark_syst_err, /string) then dark_syst_err = 0. ; rebin the error matrix and take into account the exposure time correction of l1 images dark_error = sqrt(rebin(dark_error ^ 2, header.naxis1, header.naxis2) * nbin) / nbin dark_error = dark_error * ndit1 * ndit2 ; set error values in masked areas equal to zero if count gt 0 then error[mask] = 0. ; save errors in the error info structure err_info.dark.error = dark_error err_info.dark.syst_err = dark_syst_err jump: s = where(data le 0., count) s = where(data le 0, count) if isa(quality_matrix) and count gt 0 then quality_matrix[s] = 0 if ~isa(history) then history = !null Loading
metis_dark_uvda_new.pro +154 −150 Original line number Diff line number Diff line function m_angle_uv, matrix a = mean(matrix[8 : 43, 8 : 43]) b = mean(matrix[980 : 1015, 8 : 43]) c = mean(matrix[8 : 43, 980 : 1015]) Loading @@ -9,14 +8,12 @@ function m_angle_uv, matrix end function m_center_uv, matrix a = mean(matrix[452 : 571, 460 : 579]) return, a end function correct_dark, data, dark, offset dark_bin = rebin(dark, 1024, 1024, /sample) data_bin = rebin(data, 1024, 1024, /sample) Loading @@ -30,8 +27,7 @@ function correct_dark, data, dark, offset return, dark_corr end function metis_dark_uvda_new, data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history function metis_dark_uvda_new, data, header, cal_pack, err_info, quality_matrix = quality_matrix, history = history dark = cal_pack.uv_channel.dark dit = header.dit Loading Loading @@ -126,26 +122,34 @@ function metis_dark_uvda_new, data, header, cal_pack, error = error, quality_mat ; set values in masked areas equal to zero if count gt 0 then data[mask] = 0. ; save data and dark image in the error info structure err_info.data.image = data err_info.dark.image = dark_image ; read the dark error matrix dark_error_file = dark[i].error if file_test(cal_pack.path + dark_error_file) then $ dark_error = float(readfits(cal_pack.path + dark_error_file, hdr, /silent)) else $ if file_test(cal_pack.path + dark_error_file) then begin dark_error = float(readfits(cal_pack.path + dark_error_file, hdr, /silent)) dark_syst_err = fxpar(hdr, 'SYSTEMAT', missing = 0.) endif else begin dark_error = fltarr(1024, 1024) dark_syst_err = 0. endelse ; NOTE - some systematic errors reported in the header of the dark fits files are wrong if isa(dark_syst_err, /string) then dark_syst_err = 0. ; rebin the error matrix and take into account the exposure time correction of l1 images dark_error = sqrt(rebin(dark_error ^ 2, header.naxis1, header.naxis2) * nbin) / nbin dark_error = dark_error * ndit1 * ndit2 if not isa(error) then error = 0. error += ((data > 0)/cal_pack.uv_channel.radiometry.iaps_gain.value * cal_pack.uv_channel.radiometry.f_noise.value + 2. * dark_error^2)/data^2 ; set error values in masked areas equal to zero if count gt 0 then error[mask] = 0. ; save errors in the error info structure err_info.dark.error = dark_error err_info.dark.syst_err = dark_syst_err jump: s = where(data le 0., count) s = where(data le 0, count) if isa(quality_matrix) and count gt 0 then quality_matrix[s] = 0 if ~isa(history) then history = !null Loading
