Add error calculation in radiometric cal. procedure

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))

	data = data * cal_factor

		eff_nobin = efficiency

		efficiency = rebin(efficiency, header.naxis1, header.naxis2)

		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)')

	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

 No newline at end of file