From 06d761181a6d2b7f2fd922c55b158fd9b97651cd Mon Sep 17 00:00:00 2001
From: Roberto Susino <roberto.susino@inaf.it>
Date: Tue, 18 Jan 2022 14:50:55 +0100
Subject: [PATCH] Add error calculation in flat-field corr. procedure

---
 metis_flat_field.pro | 17 ++++++++++++++---
 1 file changed, 14 insertions(+), 3 deletions(-)

diff --git a/metis_flat_field.pro b/metis_flat_field.pro
index e3a039c..fe47071 100644
--- a/metis_flat_field.pro
+++ b/metis_flat_field.pro
@@ -1,4 +1,4 @@
-function metis_flat_field, data, header, cal_pack, history = history
+function metis_flat_field, data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history
 
 	journal, 'Flat-field correction:'
 
@@ -6,11 +6,13 @@ function metis_flat_field, data, header, cal_pack, history = history
 		ff_info = cal_pack.vl_channel.flat_field
 
 		ff_image = readfits(cal_pack.path + ff_info.file_name, /silent)
+		ff_error = readfits(cal_pack.path + ff_info.file_name, exten_no = 1, /silent)
 		ff_file = ff_info.file_name
 
 		first_row = header.ref_rows eq 'Included' ? 2 : 0
 		detector_size = cal_pack.vl_channel.detector_size.value
 		ff_image = ff_image[*, first_row : first_row + detector_size - 1]
+		ff_error = ff_error[*, first_row : first_row + detector_size - 1]
 
 		journal, '  VL ref. rows = ' + header.ref_rows.tolower()
 	endif
@@ -18,19 +20,28 @@ function metis_flat_field, data, header, cal_pack, history = history
 	if header.filter.contains('UV', /fold) then begin
 		ff_info = cal_pack.uv_channel.flat_field
 
-		nbin = sqrt(header.nbin)
-		i = where(ff_info.nbin eq nbin)
+		i = where(ff_info.nbin eq sqrt(header.nbin))
 		ff_image = readfits(cal_pack.path + ff_info[i].file_name, /silent)
+		
+		; WARN - flat field for the uv channel must still be provided
+		
+		ff_error = ff_image * 0.
 		ff_file = ff_info[i].file_name
 	endif
 
+	; WARN - error for binning must be checked
+	
+	ff_nobin = ff_image
 	ff_image = rebin(ff_image, header.naxis1, header.naxis2)
+	ff_error = ff_image * sqrt(rebin((ff_error/ff_nobin)^2, header.naxis1, header.naxis2)/header.nbin)
 	
 	mask = where(ff_image le 0.)
 	ff_image[mask] = 1.
 	data = data/ff_image
 	data[mask] = 0.
 
+	if isa(error) then error += (ff_error/ff_image)^2
+
 	if ~ isa(history) then history = !null
 	history = [history, 'Flat-field correction: ', '  ' + ff_file]
 
-- 
GitLab