From 89cc7f426967ff91cb57dbce4322c2c4d46016be Mon Sep 17 00:00:00 2001
From: Roberto Susino <roberto.susino@inaf.it>
Date: Tue, 18 Jan 2022 14:51:32 +0100
Subject: [PATCH] Add error calculation in radiometric cal. procedure
---
metis_rad_cal.pro | 114 ++++++++++++++++++++++++++++++++++++----------
1 file changed, 90 insertions(+), 24 deletions(-)
diff --git a/metis_rad_cal.pro b/metis_rad_cal.pro
index 5b3b8b6..06409c0 100644
--- a/metis_rad_cal.pro
+++ b/metis_rad_cal.pro
@@ -1,47 +1,113 @@
-function metis_rad_cal, data, header, cal_pack, polarimetric = polarimetric, history = history
+function metis_rad_cal, data, header, cal_pack, polarimetric = polarimetric, norad = norad, error = error, quality_matrix = quality_matrix, history = history
if header.filter.contains('VL', /fold) then begin
+ channel = cal_pack.vl_channel
+ radiometry = channel.radiometry
- ; NOTE - this is to take into account that calibration of polarimetric images or their combination, must not include the pmp efficiency (0.5) since this is included in the demodulation formulae
+ unit_factor = radiometry.msb.value
+ unit_error = radiometry.msb.error
- if keyword_set(polarimetric) then pmp_factor = 2.D0 else pmp_factor = 1.D0
+ ; NOTE - this is to take into account that calibration of polarimetric images or their combination, must not include the pmp efficiency (0.5) since this is already included in the demodulation
+
+ if keyword_set(polarimetric) then pmp_factor = 2. else pmp_factor = 1.
- angular_pixel = cal_pack.vl_channel.angular_pixel.value * header.nbin
- rad_info = cal_pack.vl_channel.radiometry[0]
- units = cal_pack.vl_channel.cal_units
- unit_factor = rad_info.msb.value
ndit = header.ndit
end
if header.filter.contains('UV', /fold) then begin
- angular_pixel = cal_pack.uv_channel.angular_pixel.value * header.nbin
- rad_info = cal_pack.uv_channel.radiometry[0]
- units = cal_pack.uv_channel.cal_units
- unit_factor = 1.D0
- pmp_factor = 1.D0
+ channel = cal_pack.uv_channel
+ radiometry = channel.radiometry
+
+ unit_factor = 1.
+ unit_error = 0.
+
+ ; the pmp factor does not apply to the uv channel
+
+ pmp_factor = 1.
+
ndit = header.ndit1 * header.ndit2
- endif
- rad_factor = rad_info.rad_response.value * pmp_factor * $
- cal_pack.instrument.pupil_area.value * $
- angular_pixel * $
- unit_factor
+ ; read the uv/vl ratio of the door retro-illumination, to be used as correction for the optical efficiency
- xposure = header.dit/1000.D0 * ndit
- cal_factor = 1./rad_factor/xposure
+ efficiency_file = radiometry.opt_efficiency.file_name
+ efficiency = float(readfits(cal_pack.path + efficiency_file, /silent))
+
+ eff_nobin = efficiency
+ efficiency = rebin(efficiency, header.naxis1, header.naxis2)
- data = data * cal_factor
+ eff_error_file = radiometry.opt_efficiency.error
+ eff_error = float(readfits(cal_pack.path + eff_error_file, /silent))
+ eff_error = efficiency * sqrt(rebin((eff_error/eff_nobin)^2, header.naxis1, header.naxis2)/header.nbin)
+
+ s = where(efficiency le 0., count)
+ quality_matrix[s] = 0
+
+ ; WARN - temporary patch to bypass efficiency correction of the uv channel
+ ; efficiency_file = 'not applied'
+ ; efficiency = 1.
+ endif
+
+ angular_pixel = channel.angular_pixel.value * header.nbin
+ units = channel.cal_units
if ~ isa(history) then history = !null
- history = [history, 'Radiometric calibration:', ' cal. factor = ' + string(cal_factor, format = '(E8.2)') + ' ' + units + '/DN']
+ history = [history, 'Radiometric calibration:']
journal, 'Radiometric calibration:'
- journal, ' rad. response = ' + string(rad_info.rad_response.value, format = '(F0)') + ' DN/phot.'
+ journal, ' rad. response = ' + string(radiometry.rad_response.value, format = '(F0)') + ' DN/phot.'
journal, ' pupil area = ' + string(cal_pack.instrument.pupil_area.value, format = '(F0)') + ' cm²'
journal, ' pixel solid angle = ' + string(angular_pixel, format = '(E8.2)') + ' sr'
+
+ ; exposure normalisation
+
+ xposure = header.dit/1000.D0 * ndit
+ cal_factor = 1./xposure
+ cal_error = 0.
+
journal, ' exposure time = ' + string(xposure, format = '(F0)') + ' s'
journal, ' pmp corr. factor = ' + string(pmp_factor, format = '(F0)')
- journal, ' total cal. factor = ' + string(cal_factor, format = '(E8.2)') + ' ' + units + '/DN'
+
+ if keyword_set(norad) then begin
+ history = [history, ' exposure normalisation only']
+
+ journal, ' total cal. factor = exposure normalisation only'
+ endif else begin
+ ; radiometric factor
+ rad_factor = radiometry.rad_response.value * $
+ cal_pack.instrument.pupil_area.value * $
+ angular_pixel * $
+ unit_factor * $
+ pmp_factor
+
+ cal_factor = cal_factor/rad_factor
+
+ cal_error = sqrt((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)
+
+ history = [history, ' cal. factor = ' + string(cal_factor, format = '(E8.2)') + ' ' + units + '/DN']
+
+ journal, ' total cal. factor = ' + string(cal_factor, format = '(E8.2)') + ' ' + units + '/DN'
+ journal, ' cal. factor error = ' + string(cal_error, format = '(E8.2)') + ' ' + units + '/DN'
+ endelse
+
+ ; radiometric calibration
+
+ data = data * cal_factor
+
+ if isa(error) then error += cal_error^2
+
+ if header.filter.contains('UV', /fold) then begin
+
+ ; efficiency correction
+
+ data = data/abs(efficiency)
+
+ if isa(error) then error += (eff_error/efficiency)^2
+
+ history = [history, ' efficiency corr. map = ' + efficiency_file]
+
+ journal, ' UV efficiency corr. map = ' + efficiency_file
+ endif
+
return, data
-end
+end
\ No newline at end of file
--
GitLab