Add error calculation in VLDA dark corr. procedure

metis_dark_vlda.pro

-function metis_dark_vlda, data, header, cal_pack, history = history
+function metis_dark_vlda, data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history
 
 	bias = cal_pack.vl_channel.bias
 	dark = cal_pack.vl_channel.dark
@@ -6,24 +6,36 @@ function metis_dark_vlda, data, header, cal_pack, history = history
 
 	dit = header.dit/1000.D0
 	ndit = header.ndit
+	nbin = header.nbin
 	xposure = dit * ndit
-	nbin = sqrt(header.nbin)
 	tsensor = header.tsensor
 	obt_beg = header.obt_beg/1000.D0 ; check why obt_beg must be divided by 1000.
 
+	gain = cal_pack.vl_channel.radiometry.gain.value
+
 	; WARN - temporary patch to handle local l1 fits files
-	; if tsensor eq 0. then tsensor = -30.
+	if tsensor eq 0. then tsensor = -30.
 
+	; check if one or more bias+dark matrices with same dit, binning, and tsensor exist
 	for i = 0, n_elements(bias_dark) - 1 do begin
-		if bias_dark[i].dit eq dit and bias_dark[i].nbin eq nbin and abs(bias_dark[i].tsensor - tsensor) lt 5. then begin
+		if $
+			bias_dark[i].dit eq dit and $
+			bias_dark[i].nbin^2 eq nbin and $
+			abs(bias_dark[i].tsensor - tsensor) lt 5. then begin
+			
 			dark_file = bias_dark[i].file_name
 			dark_image = float(readfits(cal_pack.path + dark_file, /silent))
-			dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin^2 * ndit
+			dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin
+			dark_image = dark_image * ndit
 			
 			mask = where(data eq 0.)
 			data = data - dark_image
 			data[mask] = 0.
 
+			; WARN - for this case the error matrix is not provided so the error is not computed
+
+			if isa(error) then error += 0.
+
 			goto, jump
 		endif
 	endfor
@@ -37,6 +49,7 @@ function metis_dark_vlda, data, header, cal_pack, history = history
 		delta_obt = min(abs(obt_beg - bias_obt_available), j)
 		i = where(bias.obt_beg eq bias_obt_available[j])
 		bias_image = float(readfits(cal_pack.path + bias[i].file_name, /silent))
+		bias_error = float(readfits(cal_pack.path + bias[i].error, /silent))
 		dark_file = bias[i].file_name
 	endif else begin
 		journal, 'Error 02: VL applicable bias file not found.'
@@ -53,6 +66,7 @@ function metis_dark_vlda, data, header, cal_pack, history = history
 		delta_obt = min(abs(obt_beg - dark_obt_available), j)
 		i = where(dark.obt_beg eq dark_obt_available[j])
 		dark_image = float(readfits(cal_pack.path + dark[i].file_name, /silent))
+		dark_error = float(readfits(cal_pack.path + dark[i].error, /silent))
 		dark_file = [dark_file, dark[i].file_name]
 	endif else begin
 		journal, 'Error 03: VL applicable dark file not found.'
@@ -60,14 +74,25 @@ function metis_dark_vlda, data, header, cal_pack, history = history
 		exit, status = 3
 	endelse
 
-	bias_image = rebin(bias_image, header.naxis1, header.naxis2) * nbin^2
-	dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin^2
+	; WARN - error for binning must be checked
+
+	bias_nobin = bias_image
+	bias_image = rebin(bias_image, header.naxis1, header.naxis2) * nbin
+	bias_error = bias_image * sqrt(rebin((bias_error/bias_nobin)^2, header.naxis1, header.naxis2)/nbin)
+	
+	dark_nobin = dark_image
+	dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin
+	dark_error = dark_image * sqrt(rebin((dark_error/dark_nobin)^2, header.naxis1, header.naxis2)/nbin)
 
 	mask = where(data eq 0.)
 	data = data - (bias_image * ndit + dark_image * xposure)
 	data[mask] = 0.
 
+	if isa(error) then error += (data * gain + 2. * (bias_error * ndit)^2 + 2. * (dark_error * xposure)^2)/data^2
+
 	jump:
+	s = where(data le 0., count)
+	if isa(quality_matrix) then quality_matrix[s] = 0
 
 	if ~ isa(history) then history = !null
 	history = [history, 'Bias and dark-current corrections: ', '  ' + dark_file]