From cc820a6c0007d1d01470d349a67c85ca7d1163af Mon Sep 17 00:00:00 2001
From: Roberto Susino <roberto.susino@inaf.it>
Date: Wed, 10 Mar 2021 17:55:05 +0100
Subject: [PATCH] Version 2.3

---
 CHANGELOG.md      |   6 +++
 metis_l1_prep.pro | 114 ++++++++++++++++++++++++++++++----------------
 2 files changed, 80 insertions(+), 40 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6c1d64a..fccee5a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,6 +1,12 @@
 # CHANGELOG
 All changes to the Metis L1 pipeline are documented in this file.
 
+## 2.3 - 2021-03-10
+- Fix a bug that caused incorrect pixel counts in re-binned images. Pixel counts are now correctly multiplied by the square of the binning factor after re-binning.
+- Introduce the correction for the total effective exposure time: pixel counts are now multiplied by the number of frames that were averaged during the acquisition, so as not to change the real signal-to-noise statistics in each pixel, related to counts integrated over the total effective exposure time. The XPOSURE and NSUMEXP keywords of the FITS, which are affected, are properly updated using parameters DIT, NDIT, NDIT1, and NDIT2.
+- Add the HV_TEMP keyword in the FITS header of UV images.
+- Fix a bug that caused some keywords in the FITS header of light curves to be incorrectly ordered.
+
 ## 2.2.1 – 2021-02-23
 - Minor changes and optimisations
 
diff --git a/metis_l1_prep.pro b/metis_l1_prep.pro
index 76703d8..ed93f2f 100644
--- a/metis_l1_prep.pro
+++ b/metis_l1_prep.pro
@@ -35,7 +35,7 @@ pro metis_l1_prep
 
 	journal, 'File ' + input.file_name
 
-	; prepare the spice kernels
+	; load the spice kernels
 
 	kernels = [ $
 		input.tls_file_name, $
@@ -81,7 +81,6 @@ pro metis_l1_prep
 		naxis2 = naxis1
 	endif else begin
 		naxis = 0
-		bitpix = 8
 	endelse
 
 	; if the data product is an accumulation matrix, look for the accumulation vector extension and read it if it exists
@@ -133,12 +132,6 @@ pro metis_l1_prep
 	file_name = 'solo_L1_' + file_name_fields[2] + '_' + date_beg_string + '_V' + version + '.fits'
 	out_file_name = 'output/' + file_name
 
-	; exposure times
-
-	dit = fxpar(metadata_extension_header, 'DIT')
-	telapse = obt_end - obt_beg
-	xposure = dit/1000.
-
 	; instrument keywords
 
 	if datatype eq 0 or datatype eq 3 or datatype eq 5 or datatype eq 9 then begin
@@ -167,36 +160,73 @@ pro metis_l1_prep
 
 	; join the metadata extension header to the primary header removing unwanted keywords
 
-	i = where(strmid(primary_header, 0, 8) eq 'DATASUM ')
+	i = where(strmid(primary_header, 0, 8) eq 'COMP_RAT')
 	j = where(strmid(metadata_extension_header, 0, 8) eq 'EXTNAME ')
 	k = where(strmid(metadata_extension_header, 0, 8) eq 'TFORM1  ')
 
 	primary_header = [ $
-		primary_header[0 : i - 1], $
+		primary_header[0 : i], $
 		metadata_extension_header[j + 1 : k - 1], $
-		primary_header[i : *] $
+		primary_header[i + 1: *] $
 	]
 
 	; rebin the image if binning was applied during the acquisition
 
-	if fxpar(primary_header, 'COMPR') then bin_type = fxpar(primary_header, 'B0_BIN') else bin_type = 0
+	if datatype le 6 then begin
+		if fxpar(primary_header, 'COMPR') then $
+			bin_type = fxpar(primary_header, 'B0_BIN') else bin_type = 0
+
+		pol_id = fxpar(primary_header, 'POL_ID')
+		fpfilt = fxpar(primary_header, 'VLFPFILT', count = count)
+		if count gt 0 and pol_id eq 0 and fpfilt eq 0 and bin_type eq 0 then bin_type = 1
+
+		bin_fact = 2. ^ bin_type
+		if bin_fact gt 1. then begin
+			data = rebin(data, naxis1 / bin_fact, naxis2 / bin_fact)
+			data = data * bin_fact^2
+
+			if ~ isa(comment) then comment = !null
+			comment = [comment, 'Image was rebinned according to the commanded binning factor.']
+
+			journal, 'Data was binned during acquisition:'
+			journal, '  binning = ' + string(bin_fact, format = '(I1)') + '×' + string(bin_fact, format = '(I1)')
+			journal, '  data was correctly rebinned:'
+			journal, '    new size = ' + string(naxis1 / bin_fact, format = '(I0)') + '×' + string(naxis1 / bin_fact, format = '(I0)')
+		endif else begin
+			journal, 'Data was not binned during acquistion.'
+		endelse
+	endif
+
+	; correct for the effective exposure time
+
+	telapse = float(obt_end - obt_beg)
 
-	pol_id = fxpar(primary_header, 'POL_ID')
-	fpfilt = fxpar(primary_header, 'VLFPFILT', count = count)
-	if count gt 0 and pol_id eq 0 and fpfilt eq 0 and bin_type eq 0 then bin_type = 1
+	dit = fxpar(metadata_extension_header, 'DIT')
+	ndit = fxpar(metadata_extension_header, 'NDIT', missing = 1.)
+	ndit1 = fxpar(metadata_extension_header, 'NDIT1', missing = 1.)
+	ndit2 = fxpar(metadata_extension_header, 'NDIT2', missing = 1.)
+
+	if datatype le 6 then begin
+		xposure = dit/1000. * ndit * ndit1 * ndit2
+		nsumexp = fix(ndit * ndit1 * ndit2)
+		data = data * nsumexp
+
+		if max(data, /nan) gt 32768 then data = long(data) else data = fix(data)
+
+		datamin = min(data, /nan)
+		datamax = max(data, /nan)
+
+		fxaddpar, primary_header, 'DATAMIN', datamin
+		fxaddpar, primary_header, 'DATAMAX', datamax
 
-	bin_fact = 2. ^ bin_type
-	if bin_fact gt 1. then begin
-		data = rebin(data, naxis1 / bin_fact, naxis2 / bin_fact, /sample)
 		if ~ isa(comment) then comment = !null
-		comment = [comment, 'Image was rebinned according to the commanded binning factor.']
+		comment = [comment, 'Image values were corrected for the total exposure time.']
 
-		journal, 'Data was binned during acquisition:'
-		journal, '  binning = ' + string(bin_fact, format = '(I1)') + '×' + string(bin_fact, format = '(I1)')
-		journal, '  data was correctly rebinned:'
-		journal, '    new size = ' + string(naxis1 / bin_fact, format = '(I0)') + '×' + string(naxis1 / bin_fact, format = '(I0)')
+		if nsumexp gt 1 then journal, 'Image values were corrected for the total exposure time.'
 	endif else begin
-		journal, 'Data was not binned during acquistion.'
+		xposure = dit/1000.
+		nsumexp = 1
+		data = !null
 	endelse
 
 	; adjust the primary header (it is almost the same for all data product types)
@@ -226,27 +256,28 @@ pro metis_l1_prep
 	fxaddpar, primary_header, 'WAVEMAX', wavemax, '[nm] Max. wavelength where response > 0.05 of max.', before = 'DATAMIN'
 	fxaddpar, primary_header, 'WAVEBAND', waveband, 'Bandpass description', before = 'DATAMIN'
 	fxaddpar, primary_header, 'XPOSURE', xposure, '[s] Total effective exposure time', before = 'DATAMIN'
-	fxaddpar, primary_header, 'NSUMEXP', 1, 'Number of detector readouts summed together', before = 'DATAMIN'
+	fxaddpar, primary_header, 'NSUMEXP', nsumexp, 'Number of detector readouts summed together', before = 'DATAMIN'
 	fxaddpar, primary_header, 'TELAPSE', telapse, '[s] Elapsed time between beginning and end of observation', before = 'DATAMIN'
 	fxaddpar, primary_header, 'SOOPNAME', planning_data.soop_name, 'Name of the SOOP campaign that the data belong to', before = 'DATAMIN'
 	fxaddpar, primary_header, 'SOOPTYPE', soop_type, 'Campaign ID(s) that the data belong to', before = 'DATAMIN'
 	fxaddpar, primary_header, 'OBS_ID', planning_data.obs_id, 'Unique ID of the individual observation', before = 'DATAMIN'
 	fxaddpar, primary_header, 'TARGET', 'TBD', 'Type of target from planning', before = 'DATAMIN'
-	fxaddpar, primary_header, 'BSCALE', 1., 'Ratio of physical to array value at 0 offset', before = 'DATAMIN'
-	fxaddpar, primary_header, 'BZERO', 0., 'Physical value for the array value 0', before = 'DATAMIN'
+	fxaddpar, primary_header, 'BSCALE', 1, 'Ratio of physical to array value at 0 offset', before = 'DATAMIN'
+	fxaddpar, primary_header, 'BZERO', 0, 'Physical value for the array value 0', before = 'DATAMIN'
 	fxaddpar, primary_header, 'BTYPE', metis_datatype[datatype], 'Science data object type', before = 'DATAMIN'
 	fxaddpar, primary_header, 'BUNIT', 'DN', 'Units of physical value, after application of BSCALE and BZERO', before = 'DATAMIN'
+	fxaddpar, primary_header, 'IDB_VERS', input.idb_version, '', before = 'HDR_VERS'
+	fxaddpar, primary_header, 'INFO_URL', 'http://metis.oato.inaf.it', 'Link to more information on the instrument data', before = 'HISTORY'
+
 	if datatype le 6 then begin
+		fxaddpar, primary_header, 'NBIN1', bin_fact, 'Binning factor in the dimension 1', before = 'COMPRESS'
+		fxaddpar, primary_header, 'NBIN2', bin_fact, 'Binning factor in the dimension 2', before = 'COMPRESS'
+		fxaddpar, primary_header, 'NBIN', bin_fact * bin_fact, 'Product of all NBIN values above', before = 'COMPRESS'
 		fxaddpar, primary_header, 'PXBEG1', 1, 'First pixel that has been read out in dimension 1', before = 'COMPRESS'
 		fxaddpar, primary_header, 'PXBEG2', 1, 'First pixel that has been read out in dimension 2', before = 'COMPRESS'
 		fxaddpar, primary_header, 'PXEND1', naxis1, 'Last pixel that has been read out in dimension 1', before = 'COMPRESS'
 		fxaddpar, primary_header, 'PXEND2', naxis2, 'Last pixel that has been read out in dimension 2', before = 'COMPRESS'
 	endif
-	fxaddpar, primary_header, 'NBIN1', bin_fact, 'Binning factor in the dimension 1', before = 'COMPRESS'
-	fxaddpar, primary_header, 'NBIN2', bin_fact, 'Binning factor in the dimension 2', before = 'COMPRESS'
-	fxaddpar, primary_header, 'NBIN', bin_fact * bin_fact, 'Product of all NBIN values above', before = 'COMPRESS'
-	fxaddpar, primary_header, 'IDB_VERS', input.idb_version, '', before = 'HDR_VERS'
-	fxaddpar, primary_header, 'INFO_URL', 'http://metis.oato.inaf.it', 'Link to more information on the instrument data', before = 'HISTORY'
 
 	; read the house-keeping telemetry
 
@@ -269,8 +300,8 @@ pro metis_l1_prep
 		; fxaddpar, primary_header, 'DAC1POL4', interpol_param(hk_table, 'NIT0E06A', date_avg, empty_params = empty_params)
 		; fxaddpar, primary_header, 'DAC2POL4', interpol_param(hk_table, 'NIT0E06B', date_avg, empty_params = empty_params)
 
-		fxaddpar, primary_header, 'TSENSOR ', interpol_param(hk_table, 'NIT0E0E0', date_avg, empty_params = empty_params), '[degC] VLDA temperature'
-		fxaddpar, primary_header, 'PMPTEMP ', interpol_param(hk_table, 'NIT0L00D', date_avg, empty_params = empty_params), '[degC] PMP temperature'
+		fxaddpar, primary_header, 'TSENSOR', interpol_param(hk_table, 'NIT0E0E0', date_avg, empty_params = empty_params), '[degC] VLDA temperature'
+		fxaddpar, primary_header, 'PMPTEMP', interpol_param(hk_table, 'NIT0L00D', date_avg, empty_params = empty_params), '[degC] PMP temperature'
 
 		journal, 'Header keywords were calibrated using HK parameters:'
 		journal, '  TSENSOR = ' + string(fxpar(primary_header, 'TSENSOR'), format = '(F0)')
@@ -278,11 +309,13 @@ pro metis_l1_prep
 	endif
 
 	if datatype eq 1 or datatype eq 4 or datatype eq 6 then begin
-		fxaddpar, primary_header, 'HVU_SCR ', interpol_param(hk_table, 'NIT0E070', date_avg, empty_params = empty_params), '[raw] HVU Screen commanded voltage'
-		fxaddpar, primary_header, 'HVU_MCP ', interpol_param(hk_table, 'NIT0E071', date_avg, empty_params = empty_params), '[raw] HVU MCP commanded voltage'
-		fxaddpar, primary_header, 'HV_SCR_V', interpol_param(hk_table, 'NIT0E0B7', date_avg, empty_params = empty_params), '[V] HVU MCP applied voltage', after = 'HVU_MCP'
-		fxaddpar, primary_header, 'HV_MCP_V', interpol_param(hk_table, 'NIT0E0B6', date_avg, empty_params = empty_params), '[V] HVU Screen + MCP voltage (TBD)', after = 'HV_SCR_V'
+		fxaddpar, primary_header, 'HVU_SCR', interpol_param(hk_table, 'NIT0E070', date_avg, empty_params = empty_params), '[raw] HVU Screen commanded voltage'
+		; CMDHVSCR e CMDHVMCP
+		fxaddpar, primary_header, 'HVU_MCP', interpol_param(hk_table, 'NIT0E071', date_avg, empty_params = empty_params), '[raw] HVU MCP commanded voltage'
+		fxaddpar, primary_header, 'HV_SCR_V', interpol_param(hk_table, 'NIT0E0B7', date_avg, empty_params = empty_params), '[V] HVU Screen + MCP read voltage', after = 'HVU_MCP'
+		fxaddpar, primary_header, 'HV_MCP_V', interpol_param(hk_table, 'NIT0E0B6', date_avg, empty_params = empty_params), '[V] HVU MCP read voltage', after = 'HV_SCR_V'
 		fxaddpar, primary_header, 'HV_MCP_I', interpol_param(hk_table, 'NIT0E0BF', date_avg, empty_params = empty_params), '[uA] HVU MCP current', after = 'HV_MCP_V'
+		fxaddpar, primary_header, 'HV_TEMP', interpol_param(hk_table, 'NIT0E0B5', date_avg, empty_params = empty_params), '[degC] HVU temperature', after = 'HV_MCP_I'
 		fxaddpar, primary_header, 'TSENSOR ', interpol_param(hk_table, 'NIT0E050', date_avg, empty_params = empty_params), '[degC] UVDA temperature'
 
 		journal, 'Header keywords were calibrated using HK parameters:'
@@ -291,6 +324,7 @@ pro metis_l1_prep
 		journal, '  HV_SCR_V = ' + string(fxpar(primary_header, 'HV_SCR_V'), format = '(F0)')
 		journal, '  HV_MCP_V = ' + string(fxpar(primary_header, 'HV_MCP_V'), format = '(F0)')
 		journal, '  HV_MCP_I = ' + string(fxpar(primary_header, 'HV_MCP_I'), format = '(F0)')
+		journal, '  HV_TEMP = ' + string(fxpar(primary_header, 'HV_TEMP'), format = '(F0)')
 		journal, '  TSENSOR  = ' + string(fxpar(primary_header, 'TSENSOR'), format = '(F0)')
 	endif
 
@@ -315,7 +349,7 @@ pro metis_l1_prep
 	key = fxpar(primary_header, 'VLFPFILT', count = count)
 	if count gt 0 then begin
 		pol_id = fxpar(primary_header, 'POL_ID')
-		if pol_id ne 0 then fxaddpar, primary_header, 'VLFPFILT', 'None' else begin
+		if pol_id ne 0 then fxaddpar, primary_header, 'VLFPFILT', 'Not applicable' else begin
 			if key eq 0 then fxaddpar, primary_header, 'VLFPFILT', 'Binning'
 			if key eq 1 then fxaddpar, primary_header, 'VLFPFILT', 'Masking'
 			if key eq 2 then fxaddpar, primary_header, 'VLFPFILT', 'No filter'
@@ -379,7 +413,7 @@ pro metis_l1_prep
 
 	; add checksum and datasum to the fits header
 
-	fits_add_checksum, primary_header, image
+	fits_add_checksum, primary_header, data
 
 	mwrfits, data, out_file_name, primary_header, /no_comment, /create, /silent
 
-- 
GitLab