diff --git a/fix_fits_header.pro b/fix_fits_header.pro
new file mode 100644
index 0000000000000000000000000000000000000000..60503b4788d5eba08cfe40f716f822c48c3642a2
--- /dev/null
+++ b/fix_fits_header.pro
@@ -0,0 +1,12 @@
+pro fix_fits_header, filename
+ fits_info, filename, n_ext = n_ext, /silent
+ for ext = 0, n_ext do begin
+ data = readfits(filename, header, exten_no = ext, /silent)
+ for k = 0, n_elements(header) - 1 do begin
+ while header[k].matches('\/ +') do begin
+ header[k] = header[k].replace('/ ', '/ ')
+ endwhile
+ endfor
+ modfits, filename, 0, header, exten_no = ext
+ endfor
+end
\ No newline at end of file
diff --git a/metis_l2_prep_uv.pro b/metis_l2_prep_uv.pro
index 387bb4de994bc9ce61bfc35402fb789812ddc69f..40b1a9b66fc0391fb3ef7f6d9ede11184cf6560c 100644
--- a/metis_l2_prep_uv.pro
+++ b/metis_l2_prep_uv.pro
@@ -1,275 +1,280 @@
pro metis_l2_prep_uv
+ ; keyword defining if the detector reference frame must be used for the output
- ; keyword defining if the detector reference frame must be used for the output
+ ref_detector = 0
- ref_detector = 0
+ ; start the log
- ; start the log
+ journal, 'output/metis_l2_prep_log.txt'
- journal,'output/metis_l2_prep_log.txt'
+ ; read the auxiliary file - here we have all the inputs we need
- ; read the auxiliary file - here we have all the inputs we need
+ input = json_parse('input/contents.json', /toarray, /tostruct)
- input = json_parse('input/contents.json', /toarray, /tostruct)
+ journal, 'File ' + file_basename(input.file_name)
- journal, 'File ' + file_basename(input.file_name)
+ ; load the spice kernels
- ; load the spice kernels
+ load_spice_kernels, input.spice_kernels, kernel_list = kernel_list, kernel_version = kernel_version
- load_spice_kernels, input.spice_kernels, kernel_list = kernel_list, kernel_version = kernel_version
+ journal, 'SPICE kernel files correctly loaded:'
+ journal, ' SDK version= ' + kernel_version
- journal, 'SPICE kernel files correctly loaded:'
- journal, ' SDK version= ' + kernel_version
+ ; read the calibration package index
- ; read the calibration package index
+ cal_pack = json_parse(input.cal_pack_path + '/index.json', /toarray, /tostruct)
- cal_pack = json_parse(input.cal_pack_path + '/index.json', /toarray, /tostruct)
+ ; include the calibration package path into the cal_pack structure
- ; include the calibration package path into the cal_pack structure
+ cal_pack = create_struct('path', input.cal_pack_path + path_sep(), cal_pack)
- cal_pack = create_struct('path', input.cal_pack_path + path_sep(), cal_pack)
+ journal, 'Calibration package correctly imported:'
+ journal, ' version = ' + cal_pack.version
+ journal, ' validity range = ' + string(cal_pack.validity_range.start, cal_pack.validity_range._end, format = '(A, "-", A)')
- journal, 'Calibration package correctly imported:'
- journal, ' version = ' + cal_pack.version
- journal, ' validity range = ' + string(cal_pack.validity_range.start, cal_pack.validity_range._end, format = '(A, "-", A)')
+ ; read the primary hdu
- ; read the primary hdu
+ data = mrdfits(input.file_name, 0, primary_header, /silent)
- data = mrdfits(input.file_name, 0, primary_header, /silent)
+ ; read the quality matrix
- ; read the quality matrix
+ quality_matrix = mrdfits(input.file_name, 'quality matrix', /silent)
- quality_matrix = mrdfits(input.file_name, 'quality matrix', /silent)
+ ; convert the string header into an idl structure
- ; convert the string header into an idl structure
+ header = fits_hdr2struct(primary_header)
- header = fits_hdr2struct(primary_header)
+ journal, 'L1 FITS file correctly read:'
+ journal, ' datatype = ' + string(header.datatype, format = '(I0)')
+ journal, ' sess_num = ' + header.sess_num
+ journal, ' nbin = ' + string(sqrt(header.nbin), format = '(I0)')
- journal, 'L1 FITS file correctly read:'
- journal, ' datatype = ' + string(header.datatype, format = '(I0)')
- journal, ' sess_num = ' + header.sess_num
- journal, ' nbin = ' + string(sqrt(header.nbin), format = '(I0)')
+ ; error is the quadratic relative error
- ; error is the quadratic relative error
-
- error = 0.
+ error = 0.
- ; calibration block
+ ; calibration block
- case header.datatype of
- 1 : begin
- data = double(data)
+ case header.datatype of
+ 1: begin
+ data = double(data)
- ; ====================================
- ; for old l1 data
- ; data = data * header.nbin * header.ndit1 * header.ndit2
- ; header.xposure = header.dit/1000. * header.ndit1 * header.ndit2
- ; header.nsumexp = header.ndit1 * header.ndit2
- ; ====================================
+ ; ====================================
+ ; for old l1 data
+ ; data = data * header.nbin * header.ndit1 * header.ndit2
+ ; header.xposure = header.dit/1000. * header.ndit1 * header.ndit2
+ ; header.nsumexp = header.ndit1 * header.ndit2
+ ; ====================================
- if fix(header.hdr_vers) le 4 then begin
- data = metis_dark_uvda(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
- endif else begin
- data = metis_dark_uvda_new(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
- endelse
+ if fix(header.hdr_vers) le 4 then begin
+ data = metis_dark_uvda(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
+ endif else begin
+ data = metis_dark_uvda_new(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
+ endelse
- ; ====================================
- ; for data already subtracted of dark
- ; file = file_basename(header.parent, '.fits') + '_subdark.fits'
- ; data = mrdfits('UV subdark/' + file, 0, /silent)
- ; data = rebin(data, header.naxis1, header.naxis2)
- ; data = data * header.nbin * header.ndit1 * header.ndit2
- ; history = ['Dark correction: ', ' ' + file]
- ; ====================================
+ ; ====================================
+ ; for data already subtracted of dark
+ ; file = file_basename(header.parent, '.fits') + '_subdark.fits'
+ ; data = mrdfits('UV subdark/' + file, 0, /silent)
+ ; data = rebin(data, header.naxis1, header.naxis2)
+ ; data = data * header.nbin * header.ndit1 * header.ndit2
+ ; history = ['Dark correction: ', ' ' + file]
+ ; ====================================
- data = metis_flat_field(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
- data = metis_vignetting(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
- data = metis_rad_cal(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
+ data = metis_flat_field(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
+ data = metis_vignetting(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
+ data = metis_rad_cal(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
- ; ====================================
- ; for simple radiometric calibration
- ; cal_pack.uv_channel.cal_units = 'DN/s'
- ; ====================================
+ ; ====================================
+ ; for simple radiometric calibration
+ ; cal_pack.uv_channel.cal_units = 'DN/s'
+ ; ====================================
- ; ====================================
- ; temporary correction based on alfa-leo data
- ; data = data * 3.4
- ; history = [history, ' additional cal. factor from stars = 3.4']
- ; ====================================
-
- btype = 'UV Lyman-alpha intensity'
- bunit = cal_pack.uv_channel.cal_units
- end
+ ; ====================================
+ ; temporary correction based on alfa-leo data
+ ; data = data * 3.4
+ ; history = [history, ' additional cal. factor from stars = 3.4']
+ ; ====================================
- 4 : begin
- ; calibration of temporal noise images
- btype = 'UV temporal standard deviation'
- bunit = 'DN'
- history = !null
- end
+ btype = 'UV Lyman-alpha intensity'
+ bunit = cal_pack.uv_channel.cal_units
+ end
- 6 : begin
- ; calibration of cr/sep log matrices
- btype = 'UV cosmic-ray matrix'
- bunit = 'DN'
- history = !null
- end
-
- else : begin
- journal, 'Error 01: wrong input data product (expected data types 1, 4, or 6).'
- journal
- exit, status = 1
- end
- endcase
-
- ; definitions for the primary header
- ; version of the fits file
-
- version = string(input.l2_version, format = '(I02)')
-
- ; fits creation date
-
- date = date_conv(systime(/julian, /utc), 'FITS')
-
- ; name of the fits file
-
- file_name_fields = strsplit(header.filename, '_', /extract)
- file_name = 'solo_L2_' + file_name_fields[2] + '_' + file_name_fields[3] + '_V' + version + '.fits'
- out_file_name = 'output/' + file_name
-
- ; adjust the primary header
-
- fxaddpar, primary_header, 'FILENAME', file_name
- fxaddpar, primary_header, 'PARENT', file_basename(input.file_name)
- fxaddpar, primary_header, 'LEVEL', 'L2'
- fxaddpar, primary_header, 'CREATOR', 'metis_l2_prep_uv.pro'
- fxaddpar, primary_header, 'VERS_SW', input.sw_version
- fxaddpar, primary_header, 'VERS_CAL', cal_pack.version
- fxaddpar, primary_header, 'VERSION', version
- fxaddpar, primary_header, 'BTYPE', btype
- fxaddpar, primary_header, 'BUNIT', bunit
- fxaddpar, primary_header, 'DATAMIN', min(data, /nan)
- fxaddpar, primary_header, 'DATAMAX', max(data, /nan)
- fxaddpar, primary_header, 'WAVEBAND', cal_pack.uv_channel.name
- fxaddpar, primary_header, 'XPOSURE', header.xposure
- fxaddpar, primary_header, 'NSUMEXP', header.nsumexp
-
- sxdelpar, primary_header, 'BLANK'
-
- ; fix planning info keywords
-
- if header.soopname.startswith('unknown') then soopname = 'none' else soopname = header.soopname
- if header.obs_mode.startswith('unknown') then obs_mode = 'none' else obs_mode = header.obs_mode
- if soopname eq 'none' then sooptype = 'none' else sooptype = header.sooptype
- if obs_mode eq 'none' then obs_type = 'none' else obs_type = header.obs_type
- if soopname eq 'none' and obs_mode eq 'none' then obs_id = 'none' else obs_id = header.obs_id
- if obs_mode eq 'none' then obs_mode = 'METIS_GENERIC_OBS'
-
- fxaddpar, primary_header, 'SOOPNAME', soopname
- fxaddpar, primary_header, 'SOOPTYPE', sooptype
- fxaddpar, primary_header, 'OBS_MODE', obs_mode
- fxaddpar, primary_header, 'OBS_TYPE', obs_type
- fxaddpar, primary_header, 'OBS_ID', obs_id
-
- ; append wcs keywords
-
- wcs = metis_wcs(header, cal_pack, ref_detector = ref_detector)
- foreach element, wcs do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
-
- ; append solar ephemeris keywords
-
- ephemeris = solo_get_ephemeris(header, cal_pack)
- foreach element, ephemeris do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
-
- history = [history, 'Update WCS and solar ephemeris:', ' SKD version = ' + kernel_version]
-
- ; update the comment and history keywords
-
- fxaddpar, primary_header, 'COMMENT', 'WARNING: UV radiometric calibration is still preliminary.'
- fxaddpar, primary_header, 'COMMENT', 'Uncertainty matrix in the FITS extension is preliminary.'
-
- if ref_detector then begin
- fxaddpar, primary_header, 'COMMENT', 'Flip vertically and rotate CROTA degrees counter-clockwise'
- fxaddpar, primary_header, 'COMMENT', ' to have Solar North up.'
- endif else $
- fxaddpar, primary_header, 'COMMENT', 'Rotate CROTA degrees counter-clockwise to have Solar North up.'
-
- for k = 0, n_elements(history) - 1 do $
- fxaddpar, primary_header, 'HISTORY', history[k]
- fxaddpar, primary_header, 'HISTORY', 'L2 FITS file created on ' + date
-
- if not ref_detector then data = metis_rectify(data, 'UV')
- fits_add_checksum, primary_header, data
- mwrfits, float(data), out_file_name, primary_header, /no_comment, /create, /silent
-
- journal, 'Fits file created:'
- journal, ' file name = ' + file_basename(out_file_name)
-
- ; add the extension with the quality matrix
-
- base_header = primary_header
- sxdelpar, base_header, 'SIMPLE'
- sxdelpar, base_header, 'EXTEND'
- sxdelpar, base_header, 'DATASUM'
- sxdelpar, base_header, 'CHECKSUM'
- sxdelpar, base_header, 'COMMENT'
- sxdelpar, base_header, 'HISTORY'
-
- extension_header = base_header
- fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'extension name', before = 'LONGSTRN'
- fxaddpar, extension_header, 'BTYPE', 'Pixel quality'
- fxaddpar, extension_header, 'BUNIT', 'None'
- fxaddpar, extension_header, 'DATAMIN', min(quality_matrix, /nan)
- fxaddpar, extension_header, 'DATAMAX', max(quality_matrix, /nan)
- fxaddpar, extension_header, 'COMMENT', 'Quality matrix values:'
- fxaddpar, extension_header, 'COMMENT', ' NaN = saturated or null L0 pixel counts'
- fxaddpar, extension_header, 'COMMENT', ' 0 = unreliable pixel value'
- fxaddpar, extension_header, 'COMMENT', ' 1 = good pixel'
- if not ref_detector then quality_matrix = metis_rectify(quality_matrix, 'UV')
- fits_add_checksum, extension_header, quality_matrix
- mwrfits, float(quality_matrix), out_file_name, extension_header, /no_comment, /silent
-
- journal, 'Quality-matrix extension correctly added.'
-
- ; add the extension with the error matrix
-
- if not ref_detector then error = metis_rectify(error, 'UV')
- error_matrix = data * sqrt(error)
- extension_header = base_header
- fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN'
- fxaddpar, extension_header, 'BTYPE', 'Absolute error'
- fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan)
- fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan)
-
- fits_add_checksum, extension_header, float(error_matrix)
- mwrfits, float(error_matrix), out_file_name, extension_header, /no_comment, /silent
-
- journal, 'Error-matrix extension correctly added.'
-
- ; write the auxiliary information file
-
- output = { $
- file_name: out_file_name, $
- l1_file_name: input.file_name, $
- log_file_name: 'output/metis_l2_prep_log.txt' $
- }
-
- json_write, output, 'output/contents.json'
-
- ; unload the spice kernels
-
- load_spice_kernels, kernel_list = kernel_list, /unload
-
- journal, 'SPICE kernel files unloaded.'
-
- ; close the log
-
- journal, 'Exiting without errors.'
- journal
-
- exit, status = 0
-end
+ 4: begin
+ ; calibration of temporal noise images
+ btype = 'UV temporal standard deviation'
+ bunit = 'DN'
+ history = !null
+ end
+
+ 6: begin
+ ; calibration of cr/sep log matrices
+ btype = 'UV cosmic-ray matrix'
+ bunit = 'DN'
+ history = !null
+ end
+
+ else: begin
+ journal, 'Error 01: wrong input data product (expected data types 1, 4, or 6).'
+ journal
+ exit, status = 1
+ end
+ endcase
+
+ ; definitions for the primary header
+ ; version of the fits file
+
+ version = string(input.l2_version, format = '(I02)')
+
+ ; fits creation date
+
+ date = date_conv(systime(/julian, /utc), 'FITS')
+
+ ; name of the fits file
+
+ file_name_fields = strsplit(header.filename, '_', /extract)
+ file_name = 'solo_L2_' + file_name_fields[2] + '_' + file_name_fields[3] + '_V' + version + '.fits'
+ out_file_name = 'output/' + file_name
+
+ ; adjust the primary header
+
+ fxaddpar, primary_header, 'FILENAME', file_name
+ fxaddpar, primary_header, 'PARENT', file_basename(input.file_name)
+ fxaddpar, primary_header, 'LEVEL', 'L2'
+ fxaddpar, primary_header, 'CREATOR', 'metis_l2_prep_uv.pro'
+ fxaddpar, primary_header, 'VERS_SW', input.sw_version
+ fxaddpar, primary_header, 'VERS_CAL', cal_pack.version
+ fxaddpar, primary_header, 'VERSION', version
+ fxaddpar, primary_header, 'BTYPE', btype
+ fxaddpar, primary_header, 'BUNIT', bunit
+ fxaddpar, primary_header, 'DATAMIN', min(data, /nan)
+ fxaddpar, primary_header, 'DATAMAX', max(data, /nan)
+ fxaddpar, primary_header, 'WAVEBAND', cal_pack.uv_channel.name
+ fxaddpar, primary_header, 'XPOSURE', header.xposure
+ fxaddpar, primary_header, 'NSUMEXP', header.nsumexp
+
+ sxdelpar, primary_header, 'BLANK'
+
+ ; fix planning info keywords
+
+ if header.soopname.startswith('unknown') then soopname = 'none' else soopname = header.soopname
+ if header.obs_mode.startswith('unknown') then obs_mode = 'none' else obs_mode = header.obs_mode
+ if soopname eq 'none' then sooptype = 'none' else sooptype = header.sooptype
+ if obs_mode eq 'none' then obs_type = 'none' else obs_type = header.obs_type
+ if soopname eq 'none' and obs_mode eq 'none' then obs_id = 'none' else obs_id = header.obs_id
+ if obs_mode eq 'none' then obs_mode = 'METIS_GENERIC_OBS'
+
+ fxaddpar, primary_header, 'SOOPNAME', soopname
+ fxaddpar, primary_header, 'SOOPTYPE', sooptype
+ fxaddpar, primary_header, 'OBS_MODE', obs_mode
+ fxaddpar, primary_header, 'OBS_TYPE', obs_type
+ fxaddpar, primary_header, 'OBS_ID', obs_id
+
+ ; append wcs keywords
+
+ wcs = metis_wcs(header, cal_pack, ref_detector = ref_detector)
+ foreach element, wcs do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
+
+ ; append solar ephemeris keywords
+
+ ephemeris = solo_get_ephemeris(header, cal_pack)
+ foreach element, ephemeris do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
+
+ history = [history, 'Update WCS and solar ephemeris:', ' SKD version = ' + kernel_version]
+
+ ; update the comment and history keywords
+
+ fxaddpar, primary_header, 'COMMENT', 'WARNING: UV radiometric calibration is still preliminary.'
+ fxaddpar, primary_header, 'COMMENT', 'Uncertainty matrix in the FITS extension is preliminary.'
+
+ if ref_detector then begin
+ fxaddpar, primary_header, 'COMMENT', 'Flip vertically and rotate CROTA degrees counter-clockwise'
+ fxaddpar, primary_header, 'COMMENT', ' to have Solar North up.'
+ endif else $
+ fxaddpar, primary_header, 'COMMENT', 'Rotate CROTA degrees counter-clockwise to have Solar North up.'
+
+ for k = 0, n_elements(history) - 1 do $
+ fxaddpar, primary_header, 'HISTORY', history[k]
+ fxaddpar, primary_header, 'HISTORY', 'L2 FITS file created on ' + date
+
+ if not ref_detector then data = metis_rectify(data, 'UV')
+ fits_add_checksum, primary_header, float(data)
+ mwrfits, float(data), out_file_name, primary_header, /no_comment, /create, /silent
+
+ journal, 'Fits file created:'
+ journal, ' file name = ' + file_basename(out_file_name)
+
+ ; add the extension with the quality matrix
+
+ base_header = primary_header
+ sxdelpar, base_header, 'SIMPLE'
+ sxdelpar, base_header, 'EXTEND'
+ sxdelpar, base_header, 'DATASUM'
+ sxdelpar, base_header, 'CHECKSUM'
+ sxdelpar, base_header, 'COMMENT'
+ sxdelpar, base_header, 'HISTORY'
+
+ extension_header = base_header
+ fxaddpar, extension_header, 'XTENSION', 'IMAGE', 'image extension', before = 'BITPIX'
+ fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'extension name', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'BTYPE', 'Pixel quality'
+ fxaddpar, extension_header, 'BUNIT', 'None'
+ fxaddpar, extension_header, 'DATAMIN', min(quality_matrix, /nan)
+ fxaddpar, extension_header, 'DATAMAX', max(quality_matrix, /nan)
+ fxaddpar, extension_header, 'COMMENT', 'Quality matrix values:'
+ fxaddpar, extension_header, 'COMMENT', ' NaN = saturated or null L0 pixel counts'
+ fxaddpar, extension_header, 'COMMENT', ' 0 = unreliable pixel value'
+ fxaddpar, extension_header, 'COMMENT', ' 1 = good pixel'
+ if not ref_detector then quality_matrix = metis_rectify(quality_matrix, 'UV')
+ fits_add_checksum, extension_header, float(quality_matrix)
+ mwrfits, float(quality_matrix), out_file_name, extension_header, /no_comment, /silent
+
+ journal, 'Quality-matrix extension correctly added.'
+
+ ; add the extension with the error matrix
+
+ if not ref_detector then error = metis_rectify(error, 'UV')
+ error_matrix = data * sqrt(error)
+ extension_header = base_header
+ fxaddpar, extension_header, 'XTENSION', 'IMAGE', 'image extension', before = 'BITPIX'
+ fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'BTYPE', 'Absolute error'
+ fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan)
+ fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan)
+
+ fits_add_checksum, extension_header, float(error_matrix)
+ mwrfits, float(error_matrix), out_file_name, extension_header, /no_comment, /silent
+
+ journal, 'Error-matrix extension correctly added.'
+
+ ; fix some header keywords
+
+ fix_fits_header, out_file_name
+
+ ; write the auxiliary information file
+
+ output = { $
+ file_name: out_file_name, $
+ l1_file_name: input.file_name, $
+ log_file_name: 'output/metis_l2_prep_log.txt' $
+ }
+
+ json_write, output, 'output/contents.json'
+
+ ; unload the spice kernels
+
+ load_spice_kernels, kernel_list = kernel_list, /unload
+
+ journal, 'SPICE kernel files unloaded.'
+
+ ; close the log
+
+ journal, 'Exiting without errors.'
+ journal
+
+ exit, status = 0
+end
\ No newline at end of file
diff --git a/metis_l2_prep_vl_generic.pro b/metis_l2_prep_vl_generic.pro
index 6e983f2cfb745b3fff9d115fa0fc4ccdf08bf80b..70c9f5a8107180e46fa64d83b5f4030f8e833b68 100644
--- a/metis_l2_prep_vl_generic.pro
+++ b/metis_l2_prep_vl_generic.pro
@@ -1,311 +1,316 @@
pro metis_l2_prep_vl_generic
+ ; keyword defining if the detector reference frame must be used for the output
- ; keyword defining if the detector reference frame must be used for the output
+ ref_detector = 0
- ref_detector = 0
+ ; start the log
- ; start the log
+ journal, 'output/metis_l2_prep_log.txt'
- journal,'output/metis_l2_prep_log.txt'
+ ; read the auxiliary file - here we have all the inputs we need
- ; read the auxiliary file - here we have all the inputs we need
+ input = json_parse('input/contents.json', /toarray, /tostruct)
- input = json_parse('input/contents.json', /toarray, /tostruct)
+ ; load the spice kernels
- ; load the spice kernels
+ load_spice_kernels, input.spice_kernels, kernel_list = kernel_list, kernel_version = kernel_version
- load_spice_kernels, input.spice_kernels, kernel_list = kernel_list, kernel_version = kernel_version
+ journal, 'SPICE kernel files correctly loaded:'
+ journal, ' SDK version= ' + kernel_version
- journal, 'SPICE kernel files correctly loaded:'
- journal, ' SDK version= ' + kernel_version
+ ; read the calibration package
- ; read the calibration package
+ cal_pack = json_parse(input.cal_pack_path + '/index.json', /toarray, /tostruct)
- cal_pack = json_parse(input.cal_pack_path + '/index.json', /toarray, /tostruct)
+ ; include the calibration package path into the cal_pack structure
- ; include the calibration package path into the cal_pack structure
+ cal_pack = create_struct('path', input.cal_pack_path + path_sep(), cal_pack)
- cal_pack = create_struct('path', input.cal_pack_path + path_sep(), cal_pack)
-
- journal, 'Calibration package correctly imported:'
- journal, ' version = ' + cal_pack.version
- journal, ' validity range = ' + string(cal_pack.validity_range.start, cal_pack.validity_range._end, format = '(A, "-", A)')
-
- ; read the primary hdu
-
- data = mrdfits(input.file_name, 0, primary_header, /silent)
-
- ; patch to fix the lack of the keyword SEQ_NUM
-
- pol_id = fxpar(primary_header, 'POL_ID', missing = 0)
- if pol_id ge 1 and pol_id le 4 then begin
- seq_num = fxpar(primary_header, 'SEQ_NUM', missing = 0)
- if seq_num eq 0 then begin
- obj_cnt = fxpar(primary_header, 'OBJ_CNT')
- n_pol = fxpar(primary_header, 'N_POL', missing = 4)
- seq_num = ((obj_cnt - 1)/n_pol + 1)
- fxaddpar, primary_header, 'SEQ_NUM', seq_num, before = 'POL_ID'
- endif
- endif
-
- ; read the quality matrix
+ journal, 'Calibration package correctly imported:'
+ journal, ' version = ' + cal_pack.version
+ journal, ' validity range = ' + string(cal_pack.validity_range.start, cal_pack.validity_range._end, format = '(A, "-", A)')
- quality_matrix = mrdfits(input.file_name, 'quality matrix', /silent)
+ ; read the primary hdu
- ; convert the string header into an idl structure
+ data = mrdfits(input.file_name, 0, primary_header, /silent)
- header = fits_hdr2struct(primary_header)
+ ; patch to fix the lack of the keyword SEQ_NUM
- journal, 'L1 FITS file correctly read:'
- journal, ' filename = ' + file_basename(input.file_name)
- journal, ' datatype = ' + string(header.datatype, format = '(I0)')
- journal, ' sess_num = ' + header.sess_num
- journal, ' obj_cnt = ' + string(header.obj_cnt, format = '(I0)')
- journal, ' pol_id = ' + string(header.pol_id, format = '(I0)')
- journal, ' nbin = ' + string(sqrt(header.nbin), format = '(I0)')
-
- ; calibration block
-
- file_name_fields = strsplit(header.filename, '_', /extract)
-
- ; error is the quadratic relative error
-
- error = 0.
-
- case header.datatype of
- 0 : begin
- data = double(data)
-
- ; ====================================
- ; for old l1 data
- ; data = data * header.nbin * header.ndit
- ; header.xposure = header.xposure * header.ndit
- ; header.nsumexp = header.ndit
- ; ====================================
-
- data = metis_dark_vlda(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
-
- ; ====================================
- ; for data already subtracted of dark
- ; file = file_basename(header.parent, '.fits') + '_subdark.fits'
- ; data = mrdfits('VL subdark/' + file, 0, /silent)
- ; data = rebin(data, header.naxis1, header.naxis2)
- ; data = data * header.nbin * header.ndit
- ; history = ['Dark correction: ', ' ' + file]
- ; ====================================
-
- if header.pol_id eq 0 then begin
- btype = 'VL fixed-polarization intensity'
- descriptor = 'metis-vl-image'
- polarimetric = 1
- endif else if header.pol_id ge 1 and header.pol_id le 4 then begin
- btype = 'VL fixed-polarization intensity'
- descriptor = 'metis-vl-image'
- polarimetric = 1
- endif else begin
- btype = 'VL total brightness'
- descriptor = 'metis-vl-tb'
- polarimetric = 0
- endelse
-
- data = metis_flat_field(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
- data = metis_vignetting(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
- data = metis_rad_cal(data, header, cal_pack, polarimetric = polarimetric, error = error, quality_matrix = quality_matrix, history = history)
-
- bunit = cal_pack.vl_channel.cal_units
- end
-
- 3 : begin
- ; calibration of temporal noise images
- btype = 'VL temporal standard deviation'
- descriptor = file_name_fields[2]
- bunit = 'DN'
- history = !null
- end
-
- 5 : begin
- ; calibration of cr/sep log matrices
- btype = 'VL cosmic-ray matrix'
- descriptor = file_name_fields[2]
- bunit = 'DN'
- history = !null
- end
-
- else : begin
- journal, 'Error 01: wrong input data product (expected data types 0, 3, or 5).'
- journal
- exit, status = 1
- end
- endcase
-
- ; definitions for the primary header
- ; version of the fits file
-
- version = string(input.l2_version, format = '(I02)')
-
- ; fits creation date
-
- date = date_conv(systime(/julian, /utc), 'FITS')
-
- ; name of the fits file
-
- file_name = 'solo_L2_' + descriptor + '_' + file_name_fields[3] + '_V' + version + '.fits'
- out_file_name = 'output/' + file_name
-
- ; adjust the primary header
-
- fxaddpar, primary_header, 'FILENAME', file_name
- fxaddpar, primary_header, 'PARENT', file_basename(input.file_name)
- fxaddpar, primary_header, 'LEVEL', 'L2'
- fxaddpar, primary_header, 'CREATOR', 'metis_l2_prep_vl_generic.pro'
- fxaddpar, primary_header, 'VERS_SW', input.sw_version
- fxaddpar, primary_header, 'VERS_CAL', cal_pack.version
- fxaddpar, primary_header, 'VERSION', version
- fxaddpar, primary_header, 'BTYPE', btype
- fxaddpar, primary_header, 'BUNIT', bunit
- fxaddpar, primary_header, 'DATAMIN', min(data, /nan)
- fxaddpar, primary_header, 'DATAMAX', max(data, /nan)
- fxaddpar, primary_header, 'WAVEBAND', cal_pack.vl_channel.name
- fxaddpar, primary_header, 'XPOSURE', header.xposure
- fxaddpar, primary_header, 'NSUMEXP', header.nsumexp
-
- sxdelpar, primary_header, 'BLANK'
-
- ; fix planning info keywords
-
- if header.soopname.startswith('unknown') then soopname = 'none' else soopname = header.soopname
- if header.obs_mode.startswith('unknown') then obs_mode = 'none' else obs_mode = header.obs_mode
- if soopname eq 'none' then sooptype = 'none' else sooptype = header.sooptype
- if obs_mode eq 'none' then obs_type = 'none' else obs_type = header.obs_type
- if soopname eq 'none' and obs_mode eq 'none' then obs_id = 'none' else obs_id = header.obs_id
- if obs_mode eq 'none' then obs_mode = 'METIS_GENERIC_OBS'
-
- fxaddpar, primary_header, 'SOOPNAME', soopname
- fxaddpar, primary_header, 'SOOPTYPE', sooptype
- fxaddpar, primary_header, 'OBS_MODE', obs_mode
- fxaddpar, primary_header, 'OBS_TYPE', obs_type
- fxaddpar, primary_header, 'OBS_ID', obs_id
-
- ; read the calibration curve to convert pmp raw voltages (dacpol) into effective polarization angles
-
- dacpol_cal = cal_pack.vl_channel.dacpol_cal
-
- if header.pol_id ne 5 then begin
- if fix(header.hdr_vers) le 4 then begin
- case header.pol_id of
- 0: dacpol = header.dac1pol1
- 1: dacpol = header.dac1pol1
- 2: dacpol = header.dac1pol2
- 3: dacpol = header.dac1pol3
- 4: dacpol = header.dac1pol4
- endcase
- endif
-
- if fix(header.hdr_vers) ge 5 then begin
- dacpol = header.dac1pol1
- endif
-
- k = where(dacpol_cal.dacpol eq dacpol)
- angle = dacpol_cal.angle[k]
- angle = float(angle[0])
-
- fxaddpar, primary_header, 'POLANGLE', angle, '[deg] polarization angle', after = 'POL_ID'
- endif
-
- ; append wcs keywords
-
- wcs = metis_wcs(header, cal_pack, ref_detector = ref_detector)
- foreach element, wcs do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
-
- ; append solar ephemeris keywords
-
- ephemeris = solo_get_ephemeris(header, cal_pack)
- foreach element, ephemeris do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
-
- history = [history, 'Update WCS and solar ephemeris:', ' SKD version = ' + kernel_version]
-
- ; update the comment and history keywords
-
- fxaddpar, primary_header, 'COMMENT', 'Uncertainty matrix in the FITS extension is preliminary.'
- fxaddpar, primary_header, 'COMMENT', 'Rotate CROTA degrees counter-clockwise to have Solar North up.'
-
- for k = 0, n_elements(history) - 1 do $
- fxaddpar, primary_header, 'HISTORY', history[k]
- fxaddpar, primary_header, 'HISTORY', 'L2 FITS file created on ' + date
-
- ; add checksum and datasum to the fits header
-
- if not ref_detector then data = metis_rectify(data, 'VL')
- fits_add_checksum, primary_header, data
- mwrfits, float(data), out_file_name, primary_header, /no_comment, /create, /silent
-
- journal, 'Fits file created:'
- journal, ' file name = ' + file_basename(out_file_name)
-
- ; add the extension with the quality matrix
-
- base_header = primary_header
- sxdelpar, base_header, 'SIMPLE'
- sxdelpar, base_header, 'EXTEND'
- sxdelpar, base_header, 'DATASUM'
- sxdelpar, base_header, 'CHECKSUM'
- sxdelpar, base_header, 'COMMENT'
- sxdelpar, base_header, 'HISTORY'
-
- extension_header = base_header
- fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'extension name', before = 'LONGSTRN'
- fxaddpar, extension_header, 'BTYPE', 'Pixel quality'
- fxaddpar, extension_header, 'BUNIT', 'None'
- fxaddpar, extension_header, 'DATAMIN', min(quality_matrix, /nan)
- fxaddpar, extension_header, 'DATAMAX', max(quality_matrix, /nan)
- fxaddpar, extension_header, 'COMMENT', 'Quality matrix values:'
- fxaddpar, extension_header, 'COMMENT', ' NaN = saturated or null L0 pixel counts'
- fxaddpar, extension_header, 'COMMENT', ' 0 = unreliable pixel value'
- fxaddpar, extension_header, 'COMMENT', ' 1 = good pixel'
- if not ref_detector then quality_matrix = metis_rectify(quality_matrix, 'VL')
- fits_add_checksum, extension_header, quality_matrix
- mwrfits, float(quality_matrix), out_file_name, extension_header, /no_comment, /silent
-
- journal, 'Quality-matrix extension correctly added.'
-
- ; add the extension with the error matrix
-
- if not ref_detector then error = metis_rectify(error, 'VL')
- error_matrix = data * sqrt(error)
- extension_header = base_header
- fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN'
- fxaddpar, extension_header, 'BTYPE', 'Absolute error'
- fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan)
- fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan)
- fits_add_checksum, extension_header, float(error_matrix)
- mwrfits, float(error_matrix), out_file_name, extension_header, /no_comment, /silent
-
- journal, 'Error-matrix extension correctly added.'
+ pol_id = fxpar(primary_header, 'POL_ID', missing = 0)
+ if pol_id ge 1 and pol_id le 4 then begin
+ seq_num = fxpar(primary_header, 'SEQ_NUM', missing = 0)
+ if seq_num eq 0 then begin
+ obj_cnt = fxpar(primary_header, 'OBJ_CNT')
+ n_pol = fxpar(primary_header, 'N_POL', missing = 4)
+ seq_num = ((obj_cnt - 1) / n_pol + 1)
+ fxaddpar, primary_header, 'SEQ_NUM', seq_num, before = 'POL_ID'
+ endif
+ endif
- ; write the auxiliary information file
-
- output = { $
- file_name: out_file_name, $
- l1_file_name: input.file_name, $
- log_file_name: 'output/metis_l2_prep_log.txt' $
- }
-
- json_write, output, 'output/contents.json'
-
- ; unload the spice kernels
+ ; read the quality matrix
- load_spice_kernels, kernel_list = kernel_list, /unload
+ quality_matrix = mrdfits(input.file_name, 'quality matrix', /silent)
- journal, 'SPICE kernel files unloaded.'
+ ; convert the string header into an idl structure
- ; close the log
+ header = fits_hdr2struct(primary_header)
- journal, 'Exiting without errors.'
- journal
-
- exit, status = 0
+ journal, 'L1 FITS file correctly read:'
+ journal, ' filename = ' + file_basename(input.file_name)
+ journal, ' datatype = ' + string(header.datatype, format = '(I0)')
+ journal, ' sess_num = ' + header.sess_num
+ journal, ' obj_cnt = ' + string(header.obj_cnt, format = '(I0)')
+ journal, ' pol_id = ' + string(header.pol_id, format = '(I0)')
+ journal, ' nbin = ' + string(sqrt(header.nbin), format = '(I0)')
+
+ ; calibration block
+
+ file_name_fields = strsplit(header.filename, '_', /extract)
+
+ ; error is the quadratic relative error
+
+ error = 0.
+
+ case header.datatype of
+ 0: begin
+ data = double(data)
+
+ ; ====================================
+ ; for old l1 data
+ ; data = data * header.nbin * header.ndit
+ ; header.xposure = header.xposure * header.ndit
+ ; header.nsumexp = header.ndit
+ ; ====================================
+
+ data = metis_dark_vlda(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
+
+ ; ====================================
+ ; for data already subtracted of dark
+ ; file = file_basename(header.parent, '.fits') + '_subdark.fits'
+ ; data = mrdfits('VL subdark/' + file, 0, /silent)
+ ; data = rebin(data, header.naxis1, header.naxis2)
+ ; data = data * header.nbin * header.ndit
+ ; history = ['Dark correction: ', ' ' + file]
+ ; ====================================
+
+ if header.pol_id eq 0 then begin
+ btype = 'VL fixed-polarization intensity'
+ descriptor = 'metis-vl-image'
+ polarimetric = 1
+ endif else if header.pol_id ge 1 and header.pol_id le 4 then begin
+ btype = 'VL fixed-polarization intensity'
+ descriptor = 'metis-vl-image'
+ polarimetric = 1
+ endif else begin
+ btype = 'VL total brightness'
+ descriptor = 'metis-vl-tb'
+ polarimetric = 0
+ endelse
+
+ data = metis_flat_field(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
+ data = metis_vignetting(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
+ data = metis_rad_cal(data, header, cal_pack, polarimetric = polarimetric, error = error, quality_matrix = quality_matrix, history = history)
+
+ bunit = cal_pack.vl_channel.cal_units
+ end
+
+ 3: begin
+ ; calibration of temporal noise images
+ btype = 'VL temporal standard deviation'
+ descriptor = file_name_fields[2]
+ bunit = 'DN'
+ history = !null
+ end
+
+ 5: begin
+ ; calibration of cr/sep log matrices
+ btype = 'VL cosmic-ray matrix'
+ descriptor = file_name_fields[2]
+ bunit = 'DN'
+ history = !null
+ end
+
+ else: begin
+ journal, 'Error 01: wrong input data product (expected data types 0, 3, or 5).'
+ journal
+ exit, status = 1
+ end
+ endcase
+
+ ; definitions for the primary header
+ ; version of the fits file
+
+ version = string(input.l2_version, format = '(I02)')
+
+ ; fits creation date
+
+ date = date_conv(systime(/julian, /utc), 'FITS')
+
+ ; name of the fits file
+
+ file_name = 'solo_L2_' + descriptor + '_' + file_name_fields[3] + '_V' + version + '.fits'
+ out_file_name = 'output/' + file_name
+
+ ; adjust the primary header
+
+ fxaddpar, primary_header, 'FILENAME', file_name
+ fxaddpar, primary_header, 'PARENT', file_basename(input.file_name)
+ fxaddpar, primary_header, 'LEVEL', 'L2'
+ fxaddpar, primary_header, 'CREATOR', 'metis_l2_prep_vl_generic.pro'
+ fxaddpar, primary_header, 'VERS_SW', input.sw_version
+ fxaddpar, primary_header, 'VERS_CAL', cal_pack.version
+ fxaddpar, primary_header, 'VERSION', version
+ fxaddpar, primary_header, 'BTYPE', btype
+ fxaddpar, primary_header, 'BUNIT', bunit
+ fxaddpar, primary_header, 'DATAMIN', min(data, /nan)
+ fxaddpar, primary_header, 'DATAMAX', max(data, /nan)
+ fxaddpar, primary_header, 'WAVEBAND', cal_pack.vl_channel.name
+ fxaddpar, primary_header, 'XPOSURE', header.xposure
+ fxaddpar, primary_header, 'NSUMEXP', header.nsumexp
+
+ sxdelpar, primary_header, 'BLANK'
+
+ ; fix planning info keywords
+
+ if header.soopname.startswith('unknown') then soopname = 'none' else soopname = header.soopname
+ if header.obs_mode.startswith('unknown') then obs_mode = 'none' else obs_mode = header.obs_mode
+ if soopname eq 'none' then sooptype = 'none' else sooptype = header.sooptype
+ if obs_mode eq 'none' then obs_type = 'none' else obs_type = header.obs_type
+ if soopname eq 'none' and obs_mode eq 'none' then obs_id = 'none' else obs_id = header.obs_id
+ if obs_mode eq 'none' then obs_mode = 'METIS_GENERIC_OBS'
+
+ fxaddpar, primary_header, 'SOOPNAME', soopname
+ fxaddpar, primary_header, 'SOOPTYPE', sooptype
+ fxaddpar, primary_header, 'OBS_MODE', obs_mode
+ fxaddpar, primary_header, 'OBS_TYPE', obs_type
+ fxaddpar, primary_header, 'OBS_ID', obs_id
+
+ ; read the calibration curve to convert pmp raw voltages (dacpol) into effective polarization angles
+
+ dacpol_cal = cal_pack.vl_channel.dacpol_cal
+
+ if header.pol_id ne 5 then begin
+ if fix(header.hdr_vers) le 4 then begin
+ case header.pol_id of
+ 0: dacpol = header.dac1Pol1
+ 1: dacpol = header.dac1Pol1
+ 2: dacpol = header.dac1Pol2
+ 3: dacpol = header.dac1Pol3
+ 4: dacpol = header.dac1Pol4
+ endcase
+ endif
+
+ if fix(header.hdr_vers) ge 5 then begin
+ dacpol = header.dac1Pol1
+ endif
+
+ k = where(dacpol_cal.dacpol eq dacpol)
+ angle = dacpol_cal.angle[k]
+ angle = float(angle[0])
+
+ fxaddpar, primary_header, 'POLANGLE', angle, '[deg] polarization angle', after = 'POL_ID'
+ endif
+
+ ; append wcs keywords
+
+ wcs = metis_wcs(header, cal_pack, ref_detector = ref_detector)
+ foreach element, wcs do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
+
+ ; append solar ephemeris keywords
+
+ ephemeris = solo_get_ephemeris(header, cal_pack)
+ foreach element, ephemeris do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
+
+ history = [history, 'Update WCS and solar ephemeris:', ' SKD version = ' + kernel_version]
+
+ ; update the comment and history keywords
+
+ fxaddpar, primary_header, 'COMMENT', 'Uncertainty matrix in the FITS extension is preliminary.'
+ fxaddpar, primary_header, 'COMMENT', 'Rotate CROTA degrees counter-clockwise to have Solar North up.'
+
+ for k = 0, n_elements(history) - 1 do $
+ fxaddpar, primary_header, 'HISTORY', history[k]
+ fxaddpar, primary_header, 'HISTORY', 'L2 FITS file created on ' + date
+
+ ; add checksum and datasum to the fits header
+
+ if not ref_detector then data = metis_rectify(data, 'VL')
+ fits_add_checksum, primary_header, float(data)
+ mwrfits, float(data), out_file_name, primary_header, /no_comment, /create, /silent
+
+ journal, 'Fits file created:'
+ journal, ' file name = ' + file_basename(out_file_name)
+
+ ; add the extension with the quality matrix
+
+ base_header = primary_header
+ sxdelpar, base_header, 'SIMPLE'
+ sxdelpar, base_header, 'EXTEND'
+ sxdelpar, base_header, 'DATASUM'
+ sxdelpar, base_header, 'CHECKSUM'
+ sxdelpar, base_header, 'COMMENT'
+ sxdelpar, base_header, 'HISTORY'
+
+ extension_header = base_header
+ fxaddpar, extension_header, 'XTENSION', 'IMAGE', 'image extension', before = 'BITPIX'
+ fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'extension name', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'BTYPE', 'Pixel quality'
+ fxaddpar, extension_header, 'BUNIT', 'None'
+ fxaddpar, extension_header, 'DATAMIN', min(quality_matrix, /nan)
+ fxaddpar, extension_header, 'DATAMAX', max(quality_matrix, /nan)
+ fxaddpar, extension_header, 'COMMENT', 'Quality matrix values:'
+ fxaddpar, extension_header, 'COMMENT', ' NaN = saturated or null L0 pixel counts'
+ fxaddpar, extension_header, 'COMMENT', ' 0 = unreliable pixel value'
+ fxaddpar, extension_header, 'COMMENT', ' 1 = good pixel'
+ if not ref_detector then quality_matrix = metis_rectify(quality_matrix, 'VL')
+ fits_add_checksum, extension_header, float(quality_matrix)
+ mwrfits, float(quality_matrix), out_file_name, extension_header, /no_comment, /silent
+
+ journal, 'Quality-matrix extension correctly added.'
+
+ ; add the extension with the error matrix
+
+ if not ref_detector then error = metis_rectify(error, 'VL')
+ error_matrix = data * sqrt(error)
+ extension_header = base_header
+ fxaddpar, extension_header, 'XTENSION', 'IMAGE', 'image extension', before = 'BITPIX'
+ fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'BTYPE', 'Absolute error'
+ fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan)
+ fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan)
+ fits_add_checksum, extension_header, float(error_matrix)
+ mwrfits, float(error_matrix), out_file_name, extension_header, /no_comment, /silent
+
+ journal, 'Error-matrix extension correctly added.'
+
+ ; fix some header keywords
+
+ fix_fits_header, out_file_name
+
+ ; write the auxiliary information file
+
+ output = { $
+ file_name: out_file_name, $
+ l1_file_name: input.file_name, $
+ log_file_name: 'output/metis_l2_prep_log.txt' $
+ }
+
+ json_write, output, 'output/contents.json'
+
+ ; unload the spice kernels
+
+ load_spice_kernels, kernel_list = kernel_list, /unload
+
+ journal, 'SPICE kernel files unloaded.'
+
+ ; close the log
+
+ journal, 'Exiting without errors.'
+ journal
+
+ exit, status = 0
end
\ No newline at end of file
diff --git a/metis_l2_prep_vl_polariz.pro b/metis_l2_prep_vl_polariz.pro
index c41b19f0f4742496ab5bad5fecc816e472459b62..2cad5819ca6b75f9fc7552384a188d1dfadd5883 100755
--- a/metis_l2_prep_vl_polariz.pro
+++ b/metis_l2_prep_vl_polariz.pro
@@ -1,799 +1,809 @@
pro metis_l2_prep_vl_polariz
+ ; keyword defining if the detector reference frame must be used for the output
- ; keyword defining if the detector reference frame must be used for the output
+ ref_detector = 0
- ref_detector = 0
+ ; start the log
- ; start the log
+ journal, 'output/metis_l2_prep_log.txt'
- journal,'output/metis_l2_prep_log.txt'
+ ; read the auxiliary file - here we have all the inputs we need
- ; read the auxiliary file - here we have all the inputs we need
+ input = json_parse('input/contents.json', /toarray, /tostruct)
- input = json_parse('input/contents.json', /toarray, /tostruct)
+ ; load the spice kernels
- ; load the spice kernels
+ load_spice_kernels, input.spice_kernels, kernel_list = kernel_list, kernel_version = kernel_version
- load_spice_kernels, input.spice_kernels, kernel_list = kernel_list, kernel_version = kernel_version
+ journal, 'SPICE kernel files correctly loaded:'
+ journal, ' SDK version = ' + kernel_version
- journal, 'SPICE kernel files correctly loaded:'
- journal, ' SDK version = ' + kernel_version
+ ; read the calibration package
- ; read the calibration package
+ cal_pack = json_parse(input.cal_pack_path + '/index.json', /toarray, /tostruct)
- cal_pack = json_parse(input.cal_pack_path + '/index.json', /toarray, /tostruct)
+ ; include the calibration package path into the cal_pack structure
- ; include the calibration package path into the cal_pack structure
+ cal_pack = create_struct('path', input.cal_pack_path + path_sep(), cal_pack)
- cal_pack = create_struct('path', input.cal_pack_path + path_sep(), cal_pack)
+ journal, 'Calibration package correctly imported:'
+ journal, ' version = ' + cal_pack.version
+ journal, ' validity range = ' + string(cal_pack.validity_range.start, cal_pack.validity_range._end, format = '(A, "-", A)')
- journal, 'Calibration package correctly imported:'
- journal, ' version = ' + cal_pack.version
- journal, ' validity range = ' + string(cal_pack.validity_range.start, cal_pack.validity_range._end, format = '(A, "-", A)')
-
- n = n_elements(input.file_name)
+ n = n_elements(input.file_name)
- ; check on the number of input images (= npol = 4)
+ ; check on the number of input images (= npol = 4)
- if n ne 4 then begin
- journal, 'Error 07: wrong number of input images (expected 4).'
- journal
- exit, status = 7
- endif
+ if n ne 4 then begin
+ journal, 'Error 07: wrong number of input images (expected 4).'
+ journal
+ exit, status = 7
+ endif
- ; calibration block
+ ; calibration block
- data = !null
- abs_error = !null
- data_header = !null
- data_subdark = !null
+ data = !null
+ abs_error = !null
+ data_header = !null
+ data_subdark = !null
- quality_matrix = 1.
+ quality_matrix = 1.
- for k = 0, 3 do begin
+ for k = 0, 3 do begin
+ ; read the input image
- ; read the input image
+ image = mrdfits(input.file_name[k], 0, primary_header, /silent)
- image = mrdfits(input.file_name[k], 0, primary_header, /silent)
+ ; patch to fix the lack of the keyword SEQ_NUM
- ; patch to fix the lack of the keyword SEQ_NUM
-
- seq_num = fxpar(primary_header, 'SEQ_NUM', missing = 0)
- if seq_num eq 0 then begin
- obj_cnt = fxpar(primary_header, 'OBJ_CNT')
- n_pol = fxpar(primary_header, 'N_POL', missing = 4)
- seq_num = ((obj_cnt - 1)/n_pol + 1)
- fxaddpar, primary_header, 'SEQ_NUM', seq_num, before = 'POL_ID'
- endif
+ seq_num = fxpar(primary_header, 'SEQ_NUM', missing = 0)
+ if seq_num eq 0 then begin
+ obj_cnt = fxpar(primary_header, 'OBJ_CNT')
+ n_pol = fxpar(primary_header, 'N_POL', missing = 4)
+ seq_num = ((obj_cnt - 1) / n_pol + 1)
+ fxaddpar, primary_header, 'SEQ_NUM', seq_num, before = 'POL_ID'
+ endif
- if k eq 0 then comment = fxpar(primary_header, 'COMMENT') else begin
- sxdelpar, primary_header, 'COMMENT'
- foreach line, comment do $
- fxaddpar, primary_header, 'COMMENT', line
- endelse
+ if k eq 0 then comment = fxpar(primary_header, 'COMMENT') else begin
+ sxdelpar, primary_header, 'COMMENT'
+ foreach line, comment do $
+ fxaddpar, primary_header, 'COMMENT', line
+ endelse
- header = fits_hdr2struct(primary_header)
+ header = fits_hdr2struct(primary_header)
- ; ====================================
- ; for old l1 data
- ; image = image * header.nbin * header.ndit
- ; header.xposure = header.dit/1000. * header.ndit
- ; header.nsumexp = header.ndit
- ; ====================================
+ ; ====================================
+ ; for old l1 data
+ ; image = image * header.nbin * header.ndit
+ ; header.xposure = header.dit/1000. * header.ndit
+ ; header.nsumexp = header.ndit
+ ; ====================================
- ; check data type
+ ; check data type
- if header.datatype ne 0 then begin
- journal, 'Error 01: wrong input data product (expected data type 0).'
- journal
- exit, status = 1
- endif
-
- ; check consistency of polarization state
+ if header.datatype ne 0 then begin
+ journal, 'Error 01: wrong input data product (expected data type 0).'
+ journal
+ exit, status = 1
+ endif
- if header.pol_id lt 1 or header.pol_id gt 4 then begin
- journal, 'Error 04: image has inconsistent polarization state.'
- journal
- exit, status = 4
- endif
+ ; check consistency of polarization state
- ; check consistency of pmp raw voltages (dacpol)
+ if header.pol_id lt 1 or header.pol_id gt 4 then begin
+ journal, 'Error 04: image has inconsistent polarization state.'
+ journal
+ exit, status = 4
+ endif
- if $
- header.dac1pol1 ne header.dac2pol1 or $
- header.dac1pol2 ne header.dac2pol2 or $
- header.dac1pol3 ne header.dac2pol3 or $
- header.dac1pol4 ne header.dac2pol4 then begin
- journal, 'Error 05: image has inconsistent PMP voltages.'
- journal
- exit, status = 5
- endif
+ ; check consistency of pmp raw voltages (dacpol)
- journal, 'Reading L1 FITS file: ' + file_basename(input.file_name[k]) + ' ...'
- journal, ' datatype = ' + string(header.datatype, format = '(I0)')
- journal, ' sess_num = ' + header.sess_num
- journal, ' seq_num = ' + string(header.seq_num, format = '(I0)')
- journal, ' obj_cnt = ' + string(header.obj_cnt, format = '(I0)')
- journal, ' pol_id = ' + string(header.pol_id, format = '(I0)')
- journal, ' nbin = ' + string(sqrt(header.nbin), format = '(I0)')
+ if $
+ header.dac1Pol1 ne header.dac2Pol1 or $
+ header.dac1Pol2 ne header.dac2Pol2 or $
+ header.dac1Pol3 ne header.dac2Pol3 or $
+ header.dac1Pol4 ne header.dac2Pol4 then begin
+ journal, 'Error 05: image has inconsistent PMP voltages.'
+ journal
+ exit, status = 5
+ endif
- ; pile up images and headers
+ journal, 'Reading L1 FITS file: ' + file_basename(input.file_name[k]) + ' ...'
+ journal, ' datatype = ' + string(header.datatype, format = '(I0)')
+ journal, ' sess_num = ' + header.sess_num
+ journal, ' seq_num = ' + string(header.seq_num, format = '(I0)')
+ journal, ' obj_cnt = ' + string(header.obj_cnt, format = '(I0)')
+ journal, ' pol_id = ' + string(header.pol_id, format = '(I0)')
+ journal, ' nbin = ' + string(sqrt(header.nbin), format = '(I0)')
- data = [[[data]], [[image]]]
- data_header = [data_header, header]
-
- ; read and update the quality matrix
-
- image_quality_matrix = mrdfits(input.file_name[k], 'quality matrix', /silent)
-
- ; apply dark correction to compute stokes i and total brightness
-
- tb_history = !null
-
- image_subdark = metis_dark_vlda(image, header, error = image_error, quality_matrix = image_quality_matrix, cal_pack, history = tb_history)
-
- ; ====================================
- ; for data already subtracted of dark
- ; file = file_basename(header.parent, '.fits') + '_subdark.fits'
- ; image_subdark = mrdfits('input/subdark/' + file, 0, /silent)
- ; image_subdark = rebin(image_subdark, header.naxis1, header.naxis2) * header.nbin
- ; image_subdark = image_subdark * header.ndit
- ; tb_history = ['Dark correction: ', ' ' + file]
- ; ====================================
-
- data_subdark = [[[data_subdark]], [[image_subdark]]]
- abs_error = [[[abs_error]], [[image_subdark * sqrt(image_error)]]]
-
- quality_matrix *= image_quality_matrix
- endfor
-
- ; adjust header keywords
-
- obt_beg = min(data_header.obt_beg, i)
- obt_end = max(data_header.obt_end, j)
-
- header = data_header[0]
- header.date_beg = data_header[i].date_beg
- header.date_obs = data_header[i].date_beg
- header.date_end = data_header[j].date_end
-
- telapse = obt_end - obt_beg
- obt_avg = (obt_beg + obt_end)/2.
- date_avg = solo_obt2utc(decode_obt(obt_avg, /from_decimal))
-
- header.obt_beg = obt_beg
- header.obt_end = obt_end
- header.date_avg = date_avg
- header.tsensor = mean(data_header.tsensor)
- header.pmptemp = mean(data_header.pmptemp)
-
- journal, 'Average values computed:'
- journal, ' obt_beg = ' + string(header.obt_beg, format = '(F0)')
- journal, ' date_beg = ' + header.date_avg
- journal, ' tsensor = ' + string(header.tsensor, format = '(F0)')
- journal, ' pmptemp = ' + string(header.pmptemp, format = '(F0)')
-
- ; read the demodulation object of the calibration package
-
- demod_info = cal_pack.vl_channel.demodulation
-
- ; read the calibration curve to convert pmp raw voltages (dacpol) into effective polarization angles
-
- dacpol_cal = cal_pack.vl_channel.dacpol_cal
-
- ; apply the demodulation
-
- journal, 'Demodulation:'
-
- demod_tensor = fltarr(header.naxis1, header.naxis2, 4)
-
- angles = !null
- stokes_name = ['I', 'Q', 'U']
- stokes = dblarr(header.naxis1, header.naxis2, 3)
- stokes_subdark = dblarr(header.naxis1, header.naxis2, 3)
- stokes_abs_error = dblarr(header.naxis1, header.naxis2, 3)
- for i = 0, 2 do begin
-
- journal, ' stokes = ' + stokes_name[i]
-
- for j = 0, 3 do begin
-
- journal, ' pol_id = ' + string(data_header[j].pol_id, format = '(I0)')
-
- ; check the polarization state of the image and select the corresponding dacpol value
-
- if fix(data_header[j].hdr_vers) le 4 then begin
- case data_header[j].pol_id of
- 1: dacpol = data_header[j].dac1pol1
- 2: dacpol = data_header[j].dac1pol2
- 3: dacpol = data_header[j].dac1pol3
- 4: dacpol = data_header[j].dac1pol4
- endcase
- endif
-
- if fix(data_header[j].hdr_vers) ge 5 then begin
- dacpol = data_header[j].dac1pol1
- endif
-
- ; select the correct demodulation tensor element based on effective angle and stokes paramater
-
- k = where(dacpol_cal.dacpol eq dacpol)
- angle = dacpol_cal.angle[k]
- angle = angle[0]
-
- if isa(angles) then begin
- if n_elements(angles) lt 4 then angles = [angles, angle]
- endif else angles = angle
-
- journal, ' pol. angle = ' + string(angle, format = '(F0.1)')
-
- n = where(demod_info.angle eq angle and demod_info.stokes eq stokes_name[i], count)
-
- if count ne 1 then begin
- journal, 'Error 06: applicable demodulation-tensor element file not found.'
- journal
- exit, status = 6
- endif
-
- demod_file = demod_info[n].file_name
- demod_image = float(readfits(cal_pack.path + demod_file, /silent))
-
- ; rebin the demodulation tensor image to match image size
-
- demod_image = rebin(demod_image, header.naxis1, header.naxis2)
- demod_tensor[*, *, j] = demod_image
-
- journal, ' demod. tensor file = ' + demod_file
- endfor
-
- ; compute the stokes parameters
-
- stokes[*, *, i] = $
- demod_tensor[*, *, 0] * data[*, *, 0] + $
- demod_tensor[*, *, 1] * data[*, *, 1] + $
- demod_tensor[*, *, 2] * data[*, *, 2] + $
- demod_tensor[*, *, 3] * data[*, *, 3]
-
- ; compute the dark-subtracted stokes parameters
-
- stokes_subdark[*, *, i] = $
- demod_tensor[*, *, 0] * data_subdark[*, *, 0] + $
- demod_tensor[*, *, 1] * data_subdark[*, *, 1] + $
- demod_tensor[*, *, 2] * data_subdark[*, *, 2] + $
- demod_tensor[*, *, 3] * data_subdark[*, *, 3]
-
- ; compute the stokes errors
-
- stokes_abs_error[*, *, i] = sqrt( $
- demod_tensor[*, *, 0]^2 * abs_error[*, *, 0]^2 + $
- demod_tensor[*, *, 1]^2 * abs_error[*, *, 1]^2 + $
- demod_tensor[*, *, 2]^2 * abs_error[*, *, 2]^2 + $
- demod_tensor[*, *, 3]^2 * abs_error[*, *, 3]^2)
- endfor
-
- demod_history = 'Demodulation performed for angles ' + string(angles, format = '(3(f0.1, ", "), f0.1)') + ' deg'
-
- tb_history = [tb_history, demod_history]
- pb_history = demod_history
-
- i = reform(stokes_subdark[*, *, 0])
- i_abs_error = reform(stokes_abs_error[*, *, 0])
- q = reform(stokes[*, *, 1])
- q_abs_error = reform(stokes_abs_error[*, *, 1])
- u = reform(stokes[*, *, 2])
- u_abs_error = reform(stokes_abs_error[*, *, 2])
-
- ; compute the tb from the dark-subtracted stokes i and apply other calibrations
-
- journal, 'Calibrating total brightness...'
-
- tb_image = i
- tb_error = (i_abs_error/i)^2
- tb_quality_matrix = quality_matrix
- tb_image = metis_flat_field(tb_image, header, error = tb_error, quality_matrix = tb_quality_matrix, cal_pack, history = tb_history)
- tb_image = metis_vignetting(tb_image, header, error = tb_error, quality_matrix = tb_quality_matrix, cal_pack, history = tb_history)
- tb_image = metis_rad_cal(tb_image, header, error = tb_error, quality_matrix = tb_quality_matrix, cal_pack, /polarimetric, history = tb_history)
-
- ; compute the pb from the stokes q and u and apply other calibrations
-
- journal, 'Calibrating polarized brightness...'
-
- pb_image = sqrt(q^2 + u^2)
- pb_error = (q^2 * q_abs_error^2 + u^2 * u_abs_error^2)/pb_image^4
- pb_quality_matrix = quality_matrix
- pb_image = metis_flat_field(pb_image, header, error = pb_error, quality_matrix = pb_quality_matrix, cal_pack, history = pb_history)
- pb_image = metis_vignetting(pb_image, header, error = pb_error, quality_matrix = pb_quality_matrix, cal_pack, history = pb_history)
- pb_image = metis_rad_cal(pb_image, header, error = pb_error, quality_matrix = pb_quality_matrix, cal_pack, /polarimetric, history = pb_history)
-
- ; ====================================
- ; for simple radiometric calibration
- ; cal_pack.vl_channel.cal_units = 'DN/s'
- ; ====================================
-
- ; compute the polarization angle from the stokes q and u
-
- pol_angle = 0.5D0 * atan(stokes[*, *, 2], stokes[*, *, 1])
- pol_angle_error = (q^2 * u_abs_error^2 + u^2 * q_abs_error^2)/(2D0 * pb_image^4)/pol_angle^2
-
- pol_angle = pol_angle * !radeg
-
- journal, 'Polarization angle correctly computed.'
-
- ; definitions for the primary header
- ; version of the fits file
-
- version = string(input.l2_version, format = '(I02)')
-
- ; creation and acquisition times in utc
-
- date = date_conv(systime(/julian, /utc), 'FITS')
-
- ; adjust the primary header
-
- fxaddpar, primary_header, 'PARENT', strjoin(file_basename(input.file_name), ', ')
- fxaddpar, primary_header, 'LEVEL', 'L2'
- fxaddpar, primary_header, 'CREATOR', 'metis_l2_prep_vl_polariz.pro'
- fxaddpar, primary_header, 'VERS_SW', input.sw_version
- fxaddpar, primary_header, 'VERS_CAL', cal_pack.version
- fxaddpar, primary_header, 'VERSION', version
- fxaddpar, primary_header, 'DATE', date, 'date and time of FITS file creation'
- fxaddpar, primary_header, 'DATE-BEG', header.date_beg, 'start time of observation'
- fxaddpar, primary_header, 'DATE-OBS', header.date_obs, 'same as DATE-BEG'
- fxaddpar, primary_header, 'DATE-AVG', header.date_avg, 'average time of observation'
- fxaddpar, primary_header, 'DATE-END', header.date_end, 'end time of observation'
- fxaddpar, primary_header, 'OBT_BEG', header.obt_beg, 'start acquisition time in on-board time', format = 'F0.5'
- fxaddpar, primary_header, 'OBT_END', header.obt_end, 'end acquisition time in on-board time', format = 'F0.5'
- fxaddpar, primary_header, 'TELAPSE', telapse, '[s] elapsed time between beginning and end of observation'
- fxaddpar, primary_header, 'WAVEBAND', cal_pack.vl_channel.name
- fxaddpar, primary_header, 'XPOSURE', header.xposure
- fxaddpar, primary_header, 'NSUMEXP', header.nsumexp
- fxaddpar, primary_header, 'TSENSOR', header.tsensor
- fxaddpar, primary_header, 'PMPTEMP', header.pmptemp
-
- ; fix planning info keywords
-
- if header.soopname.startswith('unknown') then soopname = 'none' else soopname = header.soopname
- if header.obs_mode.startswith('unknown') then obs_mode = 'none' else obs_mode = header.obs_mode
- if soopname eq 'none' then sooptype = 'none' else sooptype = header.sooptype
- if obs_mode eq 'none' then obs_type = 'none' else obs_type = header.obs_type
- if soopname eq 'none' and obs_mode eq 'none' then obs_id = 'none' else obs_id = header.obs_id
- if obs_mode eq 'none' then obs_mode = 'METIS_GENERIC_OBS'
-
- fxaddpar, primary_header, 'SOOPNAME', soopname
- fxaddpar, primary_header, 'SOOPTYPE', sooptype
- fxaddpar, primary_header, 'OBS_MODE', obs_mode
- fxaddpar, primary_header, 'OBS_TYPE', obs_type
- fxaddpar, primary_header, 'OBS_ID', obs_id
-
- ; append wcs keywords
-
- wcs = metis_wcs(header, cal_pack, ref_detector = ref_detector)
- foreach element, wcs do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
-
- ; append solar ephemeris keywords
-
- ephemeris = solo_get_ephemeris(header, cal_pack)
- foreach element, ephemeris do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
-
- ; update the comment and history keywords
-
- fxaddpar, primary_header, 'COMMENT', 'Uncertainty matrix in the FITS extension is preliminary.'
- fxaddpar, primary_header, 'COMMENT', 'Rotate CROTA degrees counter-clockwise to have Solar North up.'
-
- history = ['Update WCS and solar ephemeris:', ' SKD version = ' + kernel_version]
-
- tb_history = [tb_history, history]
- pb_history = [pb_history, history]
-
- ; delete useless keywords
-
- sxdelpar, primary_header, 'BLANK'
- sxdelpar, primary_header, 'OBJ_CNT'
- sxdelpar, primary_header, 'POL_ID'
- sxdelpar, primary_header, 'SNRMIN'
- sxdelpar, primary_header, 'SNRMAX'
- sxdelpar, primary_header, 'NB_IMG'
- sxdelpar, primary_header, 'SN_MEAN1'
- sxdelpar, primary_header, 'SN_VAR1'
- sxdelpar, primary_header, 'SN_MEAN2'
- sxdelpar, primary_header, 'SN_VAR2'
- sxdelpar, primary_header, 'SN_MEAN3'
- sxdelpar, primary_header, 'SN_VAR3'
- sxdelpar, primary_header, 'SN_MEAN4'
- sxdelpar, primary_header, 'SN_VAR4'
- sxdelpar, primary_header, 'SN_MEAN5'
- sxdelpar, primary_header, 'SN_VAR5'
- sxdelpar, primary_header, 'N'
-
- date_beg_string = strmid(header.date_beg, 0, 19)
- date_beg_string = date_beg_string.replace('-', '')
- date_beg_string = date_beg_string.replace(':', '')
-
- ; array of output file names
-
- out_file_name = strarr(4)
-
- ; keywords specific for polarized brightness images
-
- primary_pb_header = primary_header
-
- ; name of the fits file
-
- file_name = 'solo_L2_metis-vl-pb_' + date_beg_string + '_V' + version + '.fits'
- out_file_name[0] = 'output/' + file_name
-
- fxaddpar, primary_pb_header, 'FILENAME', file_name
- fxaddpar, primary_pb_header, 'BTYPE', 'VL polarized brightness'
- fxaddpar, primary_pb_header, 'BUNIT', cal_pack.vl_channel.cal_units
- fxaddpar, primary_pb_header, 'DATAMIN', min(pb_image, /nan)
- fxaddpar, primary_pb_header, 'DATAMAX', max(pb_image, /nan)
-
- ; add the history keyword
-
- for k = 0, n_elements(pb_history) - 1 do fxaddpar, primary_pb_header, 'HISTORY', pb_history[k]
- fxaddpar, primary_pb_header, 'HISTORY', 'L2 FITS file created on ' + date
-
- ; add checksum and datasum to the fits header
-
- if not ref_detector then pb_image = metis_rectify(pb_image, 'VL')
- fits_add_checksum, primary_pb_header, pb_image
- mwrfits, float(pb_image), out_file_name[0], primary_pb_header, /no_comment, /create, /silent
-
- journal, 'Polarized-brightness FITS file created:'
- journal, ' file name = ' + file_basename(out_file_name[0])
-
- ; add the extension with the quality matrix
-
- base_header = primary_pb_header
- sxdelpar, base_header, 'SIMPLE'
- sxdelpar, base_header, 'EXTEND'
- sxdelpar, base_header, 'DATASUM'
- sxdelpar, base_header, 'CHECKSUM'
- sxdelpar, base_header, 'COMMENT'
- sxdelpar, base_header, 'HISTORY'
-
- extension_header = base_header
- fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'extension name', before = 'LONGSTRN'
- fxaddpar, extension_header, 'BTYPE', 'Pixel quality'
- fxaddpar, extension_header, 'BUNIT', 'None'
- fxaddpar, extension_header, 'DATAMIN', min(quality_matrix, /nan)
- fxaddpar, extension_header, 'DATAMAX', max(quality_matrix, /nan)
- fxaddpar, extension_header, 'COMMENT', 'Quality matrix values:'
- fxaddpar, extension_header, 'COMMENT', ' NaN = saturated or null L0 pixel counts'
- fxaddpar, extension_header, 'COMMENT', ' 0 = unreliable pixel value'
- fxaddpar, extension_header, 'COMMENT', ' 1 = good pixel'
- if not ref_detector then pb_quality_matrix = metis_rectify(pb_quality_matrix, 'VL')
- fits_add_checksum, extension_header, pb_quality_matrix
- mwrfits, float(pb_quality_matrix), out_file_name[0], extension_header, /no_comment, /silent
-
- journal, 'Quality-matrix extension correctly added.'
-
- ; add the extension with the error matrix
-
- if not ref_detector then pb_error = metis_rectify(pb_error, 'VL')
- error_matrix = pb_image * 0. ;pb_image * sqrt(pb_error)
- extension_header = base_header
- fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN'
- fxaddpar, extension_header, 'BTYPE', 'Absolute error'
- fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan)
- fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan)
- fits_add_checksum, extension_header, error_matrix
- mwrfits, float(error_matrix), out_file_name[0], extension_header, /no_comment, /silent
-
- journal, 'Error-matrix extension correctly added.'
-
- ; keywords specific for total brightness images
-
- primary_tb_header = primary_header
-
- ; name of the fits file
-
- file_name = 'solo_L2_metis-vl-tb_' + date_beg_string + '_V' + version + '.fits'
- out_file_name[1] = 'output/' + file_name
-
- fxaddpar, primary_tb_header, 'FILENAME', file_name
- fxaddpar, primary_tb_header, 'BTYPE', 'VL total brightness'
- fxaddpar, primary_tb_header, 'BUNIT', cal_pack.vl_channel.cal_units
- fxaddpar, primary_tb_header, 'DATAMIN', min(tb_image, /nan)
- fxaddpar, primary_tb_header, 'DATAMAX', max(tb_image, /nan)
-
- ; add the history keyword
-
- for k = 0, n_elements(tb_history) - 1 do fxaddpar, primary_tb_header, 'HISTORY', tb_history[k]
- fxaddpar, primary_tb_header, 'HISTORY', 'L2 FITS file created on ' + date
-
- ; add checksum and datasum to the fits header
-
- if not ref_detector then tb_image = metis_rectify(tb_image, 'VL')
- fits_add_checksum, primary_tb_header, tb_image
- mwrfits, float(tb_image), out_file_name[1], primary_tb_header, /no_comment, /create, /silent
-
- journal, 'Total-brightness FITS file created:'
- journal, ' file name = ' + file_basename(out_file_name[1])
-
- ; add the extension with the quality matrix
-
- base_header = primary_tb_header
- sxdelpar, base_header, 'SIMPLE'
- sxdelpar, base_header, 'EXTEND'
- sxdelpar, base_header, 'DATASUM'
- sxdelpar, base_header, 'CHECKSUM'
- sxdelpar, base_header, 'COMMENT'
- sxdelpar, base_header, 'HISTORY'
-
- extension_header = base_header
- fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'extension name', before = 'LONGSTRN'
- fxaddpar, extension_header, 'BTYPE', 'Pixel quality'
- fxaddpar, extension_header, 'BUNIT', 'None'
- fxaddpar, extension_header, 'DATAMIN', min(tb_quality_matrix, /nan)
- fxaddpar, extension_header, 'DATAMAX', max(tb_quality_matrix, /nan)
- fxaddpar, extension_header, 'COMMENT', 'Quality matrix values:'
- fxaddpar, extension_header, 'COMMENT', ' NaN = saturated or null L0 pixel counts'
- fxaddpar, extension_header, 'COMMENT', ' 0 = unreliable pixel value'
- fxaddpar, extension_header, 'COMMENT', ' 1 = good pixel'
- if not ref_detector then tb_quality_matrix = metis_rectify(tb_quality_matrix, 'VL')
- fits_add_checksum, extension_header, tb_quality_matrix
- mwrfits, float(tb_quality_matrix), out_file_name[1], extension_header, /no_comment, /silent
-
- journal, 'Quality-matrix extension correctly added.'
-
- ; add the extension with the error matrix
-
- if not ref_detector then tb_error = metis_rectify(tb_error, 'VL')
- error_matrix = tb_image * 0. ;tb_image * sqrt(tb_error)
- extension_header = base_header
- fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN'
- fxaddpar, extension_header, 'BTYPE', 'Absolute error'
- fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan)
- fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan)
- fits_add_checksum, extension_header, error_matrix
- mwrfits, float(error_matrix), out_file_name[1], extension_header, /no_comment, /silent
-
- journal, 'Error-matrix extension correctly added.'
-
- ; keywords specific for polarization-angle images
-
- primary_polangle_header = primary_header
-
- ; name of the fits file
-
- file_name = 'solo_L2_metis-vl-pol-angle_' + date_beg_string + '_V' + version + '.fits'
- out_file_name[2] = 'output/' + file_name
-
- fxaddpar, primary_polangle_header, 'FILENAME', file_name
- fxaddpar, primary_polangle_header, 'BTYPE', 'VL polarization angle'
- fxaddpar, primary_polangle_header, 'BUNIT', 'deg'
- fxaddpar, primary_polangle_header, 'DATAMIN', min(pol_angle, /nan)
- fxaddpar, primary_polangle_header, 'DATAMAX', max(pol_angle, /nan)
-
- ; add the history keyword
-
- for k = 0, n_elements(pb_history) - 1 do fxaddpar, primary_polangle_header, 'HISTORY', pb_history[k]
- fxaddpar, primary_polangle_header, 'HISTORY', 'L2 FITS file created on ' + date
-
- ; add checksum and datasum to the fits header
-
- if not ref_detector then pol_angle = metis_rectify(pol_angle, 'VL')
- fits_add_checksum, primary_polangle_header, pol_angle
- mwrfits, float(pol_angle), out_file_name[2], primary_polangle_header, /no_comment, /create, /silent
+ ; pile up images and headers
- journal, 'Polarization-angle FITS file created:'
- journal, ' file name = ' + file_basename(out_file_name[2])
-
- ; add the extension with the quality matrix
-
- base_header = primary_polangle_header
- sxdelpar, base_header, 'SIMPLE'
- sxdelpar, base_header, 'EXTEND'
- sxdelpar, base_header, 'DATASUM'
- sxdelpar, base_header, 'CHECKSUM'
- sxdelpar, base_header, 'COMMENT'
- sxdelpar, base_header, 'HISTORY'
-
- extension_header = base_header
- fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'extension name', before = 'LONGSTRN'
- fxaddpar, extension_header, 'BTYPE', 'Pixel quality'
- fxaddpar, extension_header, 'BUNIT', 'None'
- fxaddpar, extension_header, 'DATAMIN', min(quality_matrix, /nan)
- fxaddpar, extension_header, 'DATAMAX', max(quality_matrix, /nan)
- fxaddpar, extension_header, 'COMMENT', 'Quality matrix values:'
- fxaddpar, extension_header, 'COMMENT', ' NaN = saturated or null L0 pixel counts'
- fxaddpar, extension_header, 'COMMENT', ' 0 = unreliable pixel value'
- fxaddpar, extension_header, 'COMMENT', ' 1 = good pixel'
- if not ref_detector then pol_angle_quality_matrix = metis_rectify(quality_matrix, 'VL')
- fits_add_checksum, extension_header, pol_angle_quality_matrix
- mwrfits, float(pol_angle_quality_matrix), out_file_name[2], extension_header, /no_comment, /silent
-
- journal, 'Quality-matrix extension correctly added.'
-
- ; add the extension with the error matrix
-
- if not ref_detector then pol_angle_error = metis_rectify(pol_angle_error, 'VL')
- error_matrix = pol_angle * 0. ;pol_angle * sqrt(pol_angle_error)
- extension_header = base_header
- fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN'
- fxaddpar, extension_header, 'BTYPE', 'Absolute error'
- fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan)
- fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan)
- fits_add_checksum, extension_header, error_matrix
- mwrfits, float(error_matrix), out_file_name[2], extension_header, /no_comment, /silent
-
- journal, 'Error-matrix extension correctly added.'
-
- ; management of stokes images
-
- journal, 'Calibrating Stokes parameter I...'
-
- i = reform(stokes_subdark[*, *, 0])
- i = metis_flat_field(i, header, cal_pack)
- i = metis_vignetting(i, header, cal_pack)
- i = metis_rad_cal(i, header, cal_pack, /polarimetric)
-
- journal, 'Calibrating Stokes parameter Q...'
-
- q = reform(stokes[*, *, 1])
- q = metis_flat_field(q, header, cal_pack)
- q = metis_vignetting(q, header, cal_pack)
- q = metis_rad_cal(q, header, cal_pack, /polarimetric)
-
- journal, 'Calibrating Stokes parameter U...'
-
- u = reform(stokes[*, *, 2])
- u = metis_flat_field(u, header, cal_pack)
- u = metis_vignetting(u, header, cal_pack)
- u = metis_rad_cal(u, header, cal_pack, /polarimetric)
-
- ; keywords specific for stokes images
-
- primary_stokes_header = primary_header
-
- ; name of the fits file
-
- file_name = 'solo_L2_metis-vl-stokes_' + date_beg_string + '_V' + version + '.fits'
- out_file_name[3] = 'output/' + file_name
-
- fxaddpar, primary_stokes_header, 'FILENAME', file_name
- fxaddpar, primary_stokes_header, 'BTYPE', 'Stokes I'
- fxaddpar, primary_stokes_header, 'BUNIT', cal_pack.vl_channel.cal_units
- fxaddpar, primary_stokes_header, 'DATAMIN', min(i, /nan)
- fxaddpar, primary_stokes_header, 'DATAMAX', max(i, /nan)
-
- ; add the history keyword
-
- for k = 0, n_elements(tb_history) - 1 do fxaddpar, primary_stokes_header, 'HISTORY', tb_history[k]
- fxaddpar, primary_stokes_header, 'HISTORY', 'L2 FITS file created on ' + date
-
- ; add checksum and datasum to the fits header
- if not ref_detector then i = metis_rectify(i, 'VL')
- fits_add_checksum, primary_stokes_header, i
- mwrfits, float(i), out_file_name[3], primary_stokes_header, /no_comment, /create, /silent
-
- journal, 'Stokes parameters FITS file created:'
- journal, ' file name = ' + file_basename(out_file_name[3])
-
- ; add the extension with the stokes q image
-
- extension_header = primary_stokes_header
- sxdelpar, extension_header, 'SIMPLE'
- sxdelpar, extension_header, 'EXTEND'
- sxdelpar, extension_header, 'DATASUM'
- sxdelpar, extension_header, 'CHECKSUM'
- sxdelpar, extension_header, 'COMMENT'
- sxdelpar, extension_header, 'HISTORY'
-
- fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'EXTNAME', 'Stokes Q', 'extension name', before = 'LONGSTRN'
- fxaddpar, extension_header, 'BTYPE', 'Stokes Q'
- fxaddpar, extension_header, 'BUNIT', cal_pack.vl_channel.cal_units
- fxaddpar, extension_header, 'DATAMIN', min(q, /nan)
- fxaddpar, extension_header, 'DATAMAX', max(q, /nan)
- sxdelpar, extension_header, 'EXTEND'
- if not ref_detector then q = metis_rectify(q, 'VL')
- fits_add_checksum, extension_header, q
- mwrfits, float(q), out_file_name[3], extension_header, /no_comment, /silent
-
- journal, 'Q parameter extension correctly added.'
-
- ; add the extension with the stokes u image
-
- extension_header = primary_stokes_header
- sxdelpar, extension_header, 'SIMPLE'
- sxdelpar, extension_header, 'EXTEND'
- sxdelpar, extension_header, 'DATASUM'
- sxdelpar, extension_header, 'CHECKSUM'
- sxdelpar, extension_header, 'COMMENT'
- sxdelpar, extension_header, 'HISTORY'
-
- fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'EXTNAME', 'Stokes U', 'extension name', before = 'LONGSTRN'
- fxaddpar, extension_header, 'BTYPE', 'Stokes U'
- fxaddpar, extension_header, 'BUNIT', cal_pack.vl_channel.cal_units
- fxaddpar, extension_header, 'DATAMIN', min(u, /nan)
- fxaddpar, extension_header, 'DATAMAX', max(u, /nan)
- sxdelpar, extension_header, 'EXTEND'
- if not ref_detector then u = metis_rectify(u, 'VL')
- fits_add_checksum, extension_header, u
- mwrfits, float(u), out_file_name[3], extension_header, /no_comment, /silent
-
- journal, 'U parameter extension correctly added.'
-
- ; add the extension with the quality matrix
-
- base_header = primary_stokes_header
- sxdelpar, base_header, 'EXTEND'
- sxdelpar, base_header, 'DATASUM'
- sxdelpar, base_header, 'CHECKSUM'
- sxdelpar, base_header, 'COMMENT'
- sxdelpar, base_header, 'HISTORY'
-
- extension_header = base_header
- fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'extension name', before = 'LONGSTRN'
- fxaddpar, extension_header, 'BTYPE', 'Pixel quality'
- fxaddpar, extension_header, 'BUNIT', 'None'
- fxaddpar, extension_header, 'DATAMIN', min(quality_matrix, /nan)
- fxaddpar, extension_header, 'DATAMAX', max(quality_matrix, /nan)
- fxaddpar, extension_header, 'COMMENT', 'Quality matrix values:'
- fxaddpar, extension_header, 'COMMENT', ' NaN = saturated or null L0 pixel counts'
- fxaddpar, extension_header, 'COMMENT', ' 0 = unreliable pixel value'
- fxaddpar, extension_header, 'COMMENT', ' 1 = good pixel'
- if not ref_detector then quality_matrix = metis_rectify(quality_matrix, 'VL')
- fits_add_checksum, extension_header, quality_matrix
- mwrfits, float(quality_matrix), out_file_name[3], extension_header, /no_comment, /silent
-
- journal, 'Quality-matrix extension correctly added.'
-
- ; add the extension with the error matrix
-
- if not ref_detector then error_matrix = metis_rectify(error_matrix, 'VL')
- error_matrix = tb_image * 0. ;tb_image * sqrt(tb_error)
- extension_header = base_header
- fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
- fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN'
- fxaddpar, extension_header, 'BTYPE', 'Absolute error'
- fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan)
- fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan)
- fits_add_checksum, extension_header, error_matrix
- mwrfits, float(error_matrix), out_file_name[3], extension_header, /no_comment, /silent
-
- journal, 'Error-matrix extension correctly added.'
-
- ; write the auxiliary information file
-
- output = { $
- file_name: out_file_name, $
- l1_file_name: input.file_name, $
- log_file_name: 'output/metis_l2_prep_log.txt' $
- }
-
- json_write, output, 'output/contents.json'
-
- ; unload the spice kernels
-
- load_spice_kernels, kernel_list = kernel_list, /unload
-
- journal, 'SPICE kernel files unloaded.'
-
- ; close the log
-
- journal, 'Exiting without errors.'
- journal
-
- exit, status = 0
-end
+ data = [[[data]], [[image]]]
+ data_header = [data_header, header]
+
+ ; read and update the quality matrix
+
+ image_quality_matrix = mrdfits(input.file_name[k], 'quality matrix', /silent)
+
+ ; apply dark correction to compute stokes i and total brightness
+
+ tb_history = !null
+
+ image_subdark = metis_dark_vlda(image, header, error = image_error, quality_matrix = image_quality_matrix, cal_pack, history = tb_history)
+
+ ; ====================================
+ ; for data already subtracted of dark
+ ; file = file_basename(header.parent, '.fits') + '_subdark.fits'
+ ; image_subdark = mrdfits('input/subdark/' + file, 0, /silent)
+ ; image_subdark = rebin(image_subdark, header.naxis1, header.naxis2) * header.nbin
+ ; image_subdark = image_subdark * header.ndit
+ ; tb_history = ['Dark correction: ', ' ' + file]
+ ; ====================================
+
+ data_subdark = [[[data_subdark]], [[image_subdark]]]
+ abs_error = [[[abs_error]], [[image_subdark * sqrt(image_error)]]]
+
+ quality_matrix *= image_quality_matrix
+ endfor
+
+ ; adjust header keywords
+
+ obt_beg = min(data_header.obt_beg, i)
+ obt_end = max(data_header.obt_end, j)
+
+ header = data_header[0]
+ header.date_beg = data_header[i].date_beg
+ header.date_obs = data_header[i].date_beg
+ header.date_end = data_header[j].date_end
+
+ telapse = obt_end - obt_beg
+ obt_avg = (obt_beg + obt_end) / 2.
+ date_avg = solo_obt2utc(decode_obt(obt_avg, /from_decimal))
+
+ header.obt_beg = obt_beg
+ header.obt_end = obt_end
+ header.date_avg = date_avg
+ header.tsensor = mean(data_header.tsensor)
+ header.pmptemp = mean(data_header.pmptemp)
+
+ journal, 'Average values computed:'
+ journal, ' obt_beg = ' + string(header.obt_beg, format = '(F0)')
+ journal, ' date_beg = ' + header.date_avg
+ journal, ' tsensor = ' + string(header.tsensor, format = '(F0)')
+ journal, ' pmptemp = ' + string(header.pmptemp, format = '(F0)')
+
+ ; read the demodulation object of the calibration package
+
+ demod_info = cal_pack.vl_channel.demodulation
+
+ ; read the calibration curve to convert pmp raw voltages (dacpol) into effective polarization angles
+
+ dacpol_cal = cal_pack.vl_channel.dacpol_cal
+
+ ; apply the demodulation
+
+ journal, 'Demodulation:'
+
+ demod_tensor = fltarr(header.naxis1, header.naxis2, 4)
+
+ angles = !null
+ stokes_name = ['I', 'Q', 'U']
+ stokes = dblarr(header.naxis1, header.naxis2, 3)
+ stokes_subdark = dblarr(header.naxis1, header.naxis2, 3)
+ stokes_abs_error = dblarr(header.naxis1, header.naxis2, 3)
+ for i = 0, 2 do begin
+ journal, ' stokes = ' + stokes_name[i]
+
+ for j = 0, 3 do begin
+ journal, ' pol_id = ' + string(data_header[j].pol_id, format = '(I0)')
+
+ ; check the polarization state of the image and select the corresponding dacpol value
+
+ if fix(data_header[j].hdr_vers) le 4 then begin
+ case data_header[j].pol_id of
+ 1: dacpol = data_header[j].dac1Pol1
+ 2: dacpol = data_header[j].dac1Pol2
+ 3: dacpol = data_header[j].dac1Pol3
+ 4: dacpol = data_header[j].dac1Pol4
+ endcase
+ endif
+
+ if fix(data_header[j].hdr_vers) ge 5 then begin
+ dacpol = data_header[j].dac1Pol1
+ endif
+
+ ; select the correct demodulation tensor element based on effective angle and stokes paramater
+
+ k = where(dacpol_cal.dacpol eq dacpol)
+ angle = dacpol_cal.angle[k]
+ angle = angle[0]
+
+ if isa(angles) then begin
+ if n_elements(angles) lt 4 then angles = [angles, angle]
+ endif else angles = angle
+
+ journal, ' pol. angle = ' + string(angle, format = '(F0.1)')
+
+ n = where(demod_info.angle eq angle and demod_info.stokes eq stokes_name[i], count)
+
+ if count ne 1 then begin
+ journal, 'Error 06: applicable demodulation-tensor element file not found.'
+ journal
+ exit, status = 6
+ endif
+
+ demod_file = demod_info[n].file_name
+ demod_image = float(readfits(cal_pack.path + demod_file, /silent))
+
+ ; rebin the demodulation tensor image to match image size
+
+ demod_image = rebin(demod_image, header.naxis1, header.naxis2)
+ demod_tensor[*, *, j] = demod_image
+
+ journal, ' demod. tensor file = ' + demod_file
+ endfor
+
+ ; compute the stokes parameters
+
+ stokes[*, *, i] = $
+ demod_tensor[*, *, 0] * data[*, *, 0] + $
+ demod_tensor[*, *, 1] * data[*, *, 1] + $
+ demod_tensor[*, *, 2] * data[*, *, 2] + $
+ demod_tensor[*, *, 3] * data[*, *, 3]
+
+ ; compute the dark-subtracted stokes parameters
+
+ stokes_subdark[*, *, i] = $
+ demod_tensor[*, *, 0] * data_subdark[*, *, 0] + $
+ demod_tensor[*, *, 1] * data_subdark[*, *, 1] + $
+ demod_tensor[*, *, 2] * data_subdark[*, *, 2] + $
+ demod_tensor[*, *, 3] * data_subdark[*, *, 3]
+
+ ; compute the stokes errors
+
+ stokes_abs_error[*, *, i] = sqrt( $
+ demod_tensor[*, *, 0] ^ 2 * abs_error[*, *, 0] ^ 2 + $
+ demod_tensor[*, *, 1] ^ 2 * abs_error[*, *, 1] ^ 2 + $
+ demod_tensor[*, *, 2] ^ 2 * abs_error[*, *, 2] ^ 2 + $
+ demod_tensor[*, *, 3] ^ 2 * abs_error[*, *, 3] ^ 2)
+ endfor
+
+ demod_history = 'Demodulation performed for angles ' + string(angles, format = '(3(f0.1, ", "), f0.1)') + ' deg'
+
+ tb_history = [tb_history, demod_history]
+ pb_history = demod_history
+
+ i = reform(stokes_subdark[*, *, 0])
+ i_abs_error = reform(stokes_abs_error[*, *, 0])
+ q = reform(stokes[*, *, 1])
+ q_abs_error = reform(stokes_abs_error[*, *, 1])
+ u = reform(stokes[*, *, 2])
+ u_abs_error = reform(stokes_abs_error[*, *, 2])
+
+ ; compute the tb from the dark-subtracted stokes i and apply other calibrations
+
+ journal, 'Calibrating total brightness...'
+
+ tb_image = i
+ tb_error = (i_abs_error / i) ^ 2
+ tb_quality_matrix = quality_matrix
+ tb_image = metis_flat_field(tb_image, header, error = tb_error, quality_matrix = tb_quality_matrix, cal_pack, history = tb_history)
+ tb_image = metis_vignetting(tb_image, header, error = tb_error, quality_matrix = tb_quality_matrix, cal_pack, history = tb_history)
+ tb_image = metis_rad_cal(tb_image, header, error = tb_error, quality_matrix = tb_quality_matrix, cal_pack, /polarimetric, history = tb_history)
+
+ ; compute the pb from the stokes q and u and apply other calibrations
+
+ journal, 'Calibrating polarized brightness...'
+
+ pb_image = sqrt(q ^ 2 + u ^ 2)
+ pb_error = (q ^ 2 * q_abs_error ^ 2 + u ^ 2 * u_abs_error ^ 2) / pb_image ^ 4
+ pb_quality_matrix = quality_matrix
+ pb_image = metis_flat_field(pb_image, header, error = pb_error, quality_matrix = pb_quality_matrix, cal_pack, history = pb_history)
+ pb_image = metis_vignetting(pb_image, header, error = pb_error, quality_matrix = pb_quality_matrix, cal_pack, history = pb_history)
+ pb_image = metis_rad_cal(pb_image, header, error = pb_error, quality_matrix = pb_quality_matrix, cal_pack, /polarimetric, history = pb_history)
+
+ ; ====================================
+ ; for simple radiometric calibration
+ ; cal_pack.vl_channel.cal_units = 'DN/s'
+ ; ====================================
+
+ ; compute the polarization angle from the stokes q and u
+
+ pol_angle = 0.5d0 * atan(stokes[*, *, 2], stokes[*, *, 1])
+ pol_angle_error = (q ^ 2 * u_abs_error ^ 2 + u ^ 2 * q_abs_error ^ 2) / (2d0 * pb_image ^ 4) / pol_angle ^ 2
+
+ pol_angle = pol_angle * !radeg
+
+ journal, 'Polarization angle correctly computed.'
+
+ ; definitions for the primary header
+ ; version of the fits file
+
+ version = string(input.l2_version, format = '(I02)')
+
+ ; creation and acquisition times in utc
+
+ date = date_conv(systime(/julian, /utc), 'FITS')
+
+ ; adjust the primary header
+
+ fxaddpar, primary_header, 'PARENT', strjoin(file_basename(input.file_name), ', ')
+ fxaddpar, primary_header, 'LEVEL', 'L2'
+ fxaddpar, primary_header, 'CREATOR', 'metis_l2_prep_vl_polariz.pro'
+ fxaddpar, primary_header, 'VERS_SW', input.sw_version
+ fxaddpar, primary_header, 'VERS_CAL', cal_pack.version
+ fxaddpar, primary_header, 'VERSION', version
+ fxaddpar, primary_header, 'DATE', date, 'date and time of FITS file creation'
+ fxaddpar, primary_header, 'DATE-BEG', header.date_beg, 'start time of observation'
+ fxaddpar, primary_header, 'DATE-OBS', header.date_obs, 'same as DATE-BEG'
+ fxaddpar, primary_header, 'DATE-AVG', header.date_avg, 'average time of observation'
+ fxaddpar, primary_header, 'DATE-END', header.date_end, 'end time of observation'
+ fxaddpar, primary_header, 'OBT_BEG', header.obt_beg, 'start acquisition time in on-board time', format = 'F0.5'
+ fxaddpar, primary_header, 'OBT_END', header.obt_end, 'end acquisition time in on-board time', format = 'F0.5'
+ fxaddpar, primary_header, 'TELAPSE', telapse, '[s] elapsed time between beginning and end of observation'
+ fxaddpar, primary_header, 'WAVEBAND', cal_pack.vl_channel.name
+ fxaddpar, primary_header, 'XPOSURE', header.xposure
+ fxaddpar, primary_header, 'NSUMEXP', header.nsumexp
+ fxaddpar, primary_header, 'TSENSOR', header.tsensor
+ fxaddpar, primary_header, 'PMPTEMP', header.pmptemp
+
+ ; fix planning info keywords
+
+ if header.soopname.startswith('unknown') then soopname = 'none' else soopname = header.soopname
+ if header.obs_mode.startswith('unknown') then obs_mode = 'none' else obs_mode = header.obs_mode
+ if soopname eq 'none' then sooptype = 'none' else sooptype = header.sooptype
+ if obs_mode eq 'none' then obs_type = 'none' else obs_type = header.obs_type
+ if soopname eq 'none' and obs_mode eq 'none' then obs_id = 'none' else obs_id = header.obs_id
+ if obs_mode eq 'none' then obs_mode = 'METIS_GENERIC_OBS'
+
+ fxaddpar, primary_header, 'SOOPNAME', soopname
+ fxaddpar, primary_header, 'SOOPTYPE', sooptype
+ fxaddpar, primary_header, 'OBS_MODE', obs_mode
+ fxaddpar, primary_header, 'OBS_TYPE', obs_type
+ fxaddpar, primary_header, 'OBS_ID', obs_id
+
+ ; append wcs keywords
+
+ wcs = metis_wcs(header, cal_pack, ref_detector = ref_detector)
+ foreach element, wcs do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
+
+ ; append solar ephemeris keywords
+
+ ephemeris = solo_get_ephemeris(header, cal_pack)
+ foreach element, ephemeris do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
+
+ ; update the comment and history keywords
+
+ fxaddpar, primary_header, 'COMMENT', 'Uncertainty matrix in the FITS extension is preliminary.'
+ fxaddpar, primary_header, 'COMMENT', 'Rotate CROTA degrees counter-clockwise to have Solar North up.'
+
+ history = ['Update WCS and solar ephemeris:', ' SKD version = ' + kernel_version]
+
+ tb_history = [tb_history, history]
+ pb_history = [pb_history, history]
+
+ ; delete useless keywords
+
+ sxdelpar, primary_header, 'BLANK'
+ sxdelpar, primary_header, 'OBJ_CNT'
+ sxdelpar, primary_header, 'POL_ID'
+ sxdelpar, primary_header, 'SNRMIN'
+ sxdelpar, primary_header, 'SNRMAX'
+ sxdelpar, primary_header, 'NB_IMG'
+ sxdelpar, primary_header, 'SN_MEAN1'
+ sxdelpar, primary_header, 'SN_VAR1'
+ sxdelpar, primary_header, 'SN_MEAN2'
+ sxdelpar, primary_header, 'SN_VAR2'
+ sxdelpar, primary_header, 'SN_MEAN3'
+ sxdelpar, primary_header, 'SN_VAR3'
+ sxdelpar, primary_header, 'SN_MEAN4'
+ sxdelpar, primary_header, 'SN_VAR4'
+ sxdelpar, primary_header, 'SN_MEAN5'
+ sxdelpar, primary_header, 'SN_VAR5'
+ sxdelpar, primary_header, 'N'
+
+ date_beg_string = strmid(header.date_beg, 0, 19)
+ date_beg_string = date_beg_string.replace('-', '')
+ date_beg_string = date_beg_string.replace(':', '')
+
+ ; array of output file names
+
+ out_file_name = strarr(4)
+
+ ; keywords specific for polarized brightness images
+
+ primary_pb_header = primary_header
+
+ ; name of the fits file
+
+ file_name = 'solo_L2_metis-vl-pb_' + date_beg_string + '_V' + version + '.fits'
+ out_file_name[0] = 'output/' + file_name
+
+ fxaddpar, primary_pb_header, 'FILENAME', file_name
+ fxaddpar, primary_pb_header, 'BTYPE', 'VL polarized brightness'
+ fxaddpar, primary_pb_header, 'BUNIT', cal_pack.vl_channel.cal_units
+ fxaddpar, primary_pb_header, 'DATAMIN', min(pb_image, /nan)
+ fxaddpar, primary_pb_header, 'DATAMAX', max(pb_image, /nan)
+
+ ; add the history keyword
+
+ for k = 0, n_elements(pb_history) - 1 do fxaddpar, primary_pb_header, 'HISTORY', pb_history[k]
+ fxaddpar, primary_pb_header, 'HISTORY', 'L2 FITS file created on ' + date
+
+ ; add checksum and datasum to the fits header
+
+ if not ref_detector then pb_image = metis_rectify(pb_image, 'VL')
+ fits_add_checksum, primary_pb_header, float(pb_image)
+ mwrfits, float(pb_image), out_file_name[0], primary_pb_header, /no_comment, /create, /silent
+
+ journal, 'Polarized-brightness FITS file created:'
+ journal, ' file name = ' + file_basename(out_file_name[0])
+
+ ; add the extension with the quality matrix
+
+ base_header = primary_pb_header
+ sxdelpar, base_header, 'SIMPLE'
+ sxdelpar, base_header, 'EXTEND'
+ sxdelpar, base_header, 'DATASUM'
+ sxdelpar, base_header, 'CHECKSUM'
+ sxdelpar, base_header, 'COMMENT'
+ sxdelpar, base_header, 'HISTORY'
+
+ extension_header = base_header
+ fxaddpar, extension_header, 'XTENSION', 'IMAGE', 'image extension', before = 'BITPIX'
+ fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'extension name', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'BTYPE', 'Pixel quality'
+ fxaddpar, extension_header, 'BUNIT', 'None'
+ fxaddpar, extension_header, 'DATAMIN', min(quality_matrix, /nan)
+ fxaddpar, extension_header, 'DATAMAX', max(quality_matrix, /nan)
+ fxaddpar, extension_header, 'COMMENT', 'Quality matrix values:'
+ fxaddpar, extension_header, 'COMMENT', ' NaN = saturated or null L0 pixel counts'
+ fxaddpar, extension_header, 'COMMENT', ' 0 = unreliable pixel value'
+ fxaddpar, extension_header, 'COMMENT', ' 1 = good pixel'
+ if not ref_detector then pb_quality_matrix = metis_rectify(pb_quality_matrix, 'VL')
+ fits_add_checksum, extension_header, float(pb_quality_matrix)
+ mwrfits, float(pb_quality_matrix), out_file_name[0], extension_header, /no_comment, /silent
+
+ journal, 'Quality-matrix extension correctly added.'
+
+ ; add the extension with the error matrix
+
+ if not ref_detector then pb_error = metis_rectify(pb_error, 'VL')
+ error_matrix = pb_image * 0. ; pb_image * sqrt(pb_error)
+ extension_header = base_header
+ fxaddpar, extension_header, 'XTENSION', 'IMAGE', 'image extension', before = 'BITPIX'
+ fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'BTYPE', 'Absolute error'
+ fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan)
+ fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan)
+ fits_add_checksum, extension_header, float(error_matrix)
+ mwrfits, float(error_matrix), out_file_name[0], extension_header, /no_comment, /silent
+
+ journal, 'Error-matrix extension correctly added.'
+
+ ; fix some header keywords
+
+ fix_fits_header, out_file_name[0]
+
+ ; keywords specific for total brightness images
+
+ primary_tb_header = primary_header
+
+ ; name of the fits file
+
+ file_name = 'solo_L2_metis-vl-tb_' + date_beg_string + '_V' + version + '.fits'
+ out_file_name[1] = 'output/' + file_name
+
+ fxaddpar, primary_tb_header, 'FILENAME', file_name
+ fxaddpar, primary_tb_header, 'BTYPE', 'VL total brightness'
+ fxaddpar, primary_tb_header, 'BUNIT', cal_pack.vl_channel.cal_units
+ fxaddpar, primary_tb_header, 'DATAMIN', min(tb_image, /nan)
+ fxaddpar, primary_tb_header, 'DATAMAX', max(tb_image, /nan)
+
+ ; add the history keyword
+
+ for k = 0, n_elements(tb_history) - 1 do fxaddpar, primary_tb_header, 'HISTORY', tb_history[k]
+ fxaddpar, primary_tb_header, 'HISTORY', 'L2 FITS file created on ' + date
+
+ ; add checksum and datasum to the fits header
+
+ if not ref_detector then tb_image = metis_rectify(tb_image, 'VL')
+ fits_add_checksum, primary_tb_header, float(tb_image)
+ mwrfits, float(tb_image), out_file_name[1], primary_tb_header, /no_comment, /create, /silent
+
+ journal, 'Total-brightness FITS file created:'
+ journal, ' file name = ' + file_basename(out_file_name[1])
+
+ ; add the extension with the quality matrix
+
+ base_header = primary_tb_header
+ sxdelpar, base_header, 'SIMPLE'
+ sxdelpar, base_header, 'EXTEND'
+ sxdelpar, base_header, 'DATASUM'
+ sxdelpar, base_header, 'CHECKSUM'
+ sxdelpar, base_header, 'COMMENT'
+ sxdelpar, base_header, 'HISTORY'
+
+ extension_header = base_header
+ fxaddpar, extension_header, 'XTENSION', 'IMAGE', 'image extension', before = 'BITPIX'
+ fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'extension name', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'BTYPE', 'Pixel quality'
+ fxaddpar, extension_header, 'BUNIT', 'None'
+ fxaddpar, extension_header, 'DATAMIN', min(tb_quality_matrix, /nan)
+ fxaddpar, extension_header, 'DATAMAX', max(tb_quality_matrix, /nan)
+ fxaddpar, extension_header, 'COMMENT', 'Quality matrix values:'
+ fxaddpar, extension_header, 'COMMENT', ' NaN = saturated or null L0 pixel counts'
+ fxaddpar, extension_header, 'COMMENT', ' 0 = unreliable pixel value'
+ fxaddpar, extension_header, 'COMMENT', ' 1 = good pixel'
+ if not ref_detector then tb_quality_matrix = metis_rectify(tb_quality_matrix, 'VL')
+ fits_add_checksum, extension_header, float(tb_quality_matrix)
+ mwrfits, float(tb_quality_matrix), out_file_name[1], extension_header, /no_comment, /silent
+
+ journal, 'Quality-matrix extension correctly added.'
+
+ ; add the extension with the error matrix
+
+ if not ref_detector then tb_error = metis_rectify(tb_error, 'VL')
+ error_matrix = tb_image * 0. ; tb_image * sqrt(tb_error)
+ extension_header = base_header
+ fxaddpar, extension_header, 'XTENSION', 'IMAGE', 'image extension', before = 'BITPIX'
+ fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'BTYPE', 'Absolute error'
+ fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan)
+ fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan)
+ fits_add_checksum, extension_header, float(error_matrix)
+ mwrfits, float(error_matrix), out_file_name[1], extension_header, /no_comment, /silent
+
+ journal, 'Error-matrix extension correctly added.'
+
+ ; fix some header keywords
+
+ fix_fits_header, out_file_name[1]
+
+ ; keywords specific for polarization-angle images
+
+ primary_polangle_header = primary_header
+
+ ; name of the fits file
+
+ file_name = 'solo_L2_metis-vl-pol-angle_' + date_beg_string + '_V' + version + '.fits'
+ out_file_name[2] = 'output/' + file_name
+
+ fxaddpar, primary_polangle_header, 'FILENAME', file_name
+ fxaddpar, primary_polangle_header, 'BTYPE', 'VL polarization angle'
+ fxaddpar, primary_polangle_header, 'BUNIT', 'deg'
+ fxaddpar, primary_polangle_header, 'DATAMIN', min(pol_angle, /nan)
+ fxaddpar, primary_polangle_header, 'DATAMAX', max(pol_angle, /nan)
+
+ ; add the history keyword
+
+ for k = 0, n_elements(pb_history) - 1 do fxaddpar, primary_polangle_header, 'HISTORY', pb_history[k]
+ fxaddpar, primary_polangle_header, 'HISTORY', 'L2 FITS file created on ' + date
+
+ ; add checksum and datasum to the fits header
+
+ if not ref_detector then pol_angle = metis_rectify(pol_angle, 'VL')
+ fits_add_checksum, primary_polangle_header, float(pol_angle)
+ mwrfits, float(pol_angle), out_file_name[2], primary_polangle_header, /no_comment, /create, /silent
+
+ journal, 'Polarization-angle FITS file created:'
+ journal, ' file name = ' + file_basename(out_file_name[2])
+
+ ; add the extension with the quality matrix
+
+ base_header = primary_polangle_header
+ sxdelpar, base_header, 'SIMPLE'
+ sxdelpar, base_header, 'EXTEND'
+ sxdelpar, base_header, 'DATASUM'
+ sxdelpar, base_header, 'CHECKSUM'
+ sxdelpar, base_header, 'COMMENT'
+ sxdelpar, base_header, 'HISTORY'
+
+ extension_header = base_header
+ fxaddpar, extension_header, 'XTENSION', 'IMAGE', 'image extension', before = 'BITPIX'
+ fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'extension name', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'BTYPE', 'Pixel quality'
+ fxaddpar, extension_header, 'BUNIT', 'None'
+ fxaddpar, extension_header, 'DATAMIN', min(quality_matrix, /nan)
+ fxaddpar, extension_header, 'DATAMAX', max(quality_matrix, /nan)
+ fxaddpar, extension_header, 'COMMENT', 'Quality matrix values:'
+ fxaddpar, extension_header, 'COMMENT', ' NaN = saturated or null L0 pixel counts'
+ fxaddpar, extension_header, 'COMMENT', ' 0 = unreliable pixel value'
+ fxaddpar, extension_header, 'COMMENT', ' 1 = good pixel'
+ if not ref_detector then pol_angle_quality_matrix = metis_rectify(quality_matrix, 'VL')
+ fits_add_checksum, extension_header, float(pol_angle_quality_matrix)
+ mwrfits, float(pol_angle_quality_matrix), out_file_name[2], extension_header, /no_comment, /silent
+
+ journal, 'Quality-matrix extension correctly added.'
+
+ ; add the extension with the error matrix
+
+ if not ref_detector then pol_angle_error = metis_rectify(pol_angle_error, 'VL')
+ error_matrix = pol_angle * 0. ; pol_angle * sqrt(pol_angle_error)
+ extension_header = base_header
+ fxaddpar, extension_header, 'XTENSION', 'IMAGE', 'image extension', before = 'BITPIX'
+ fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'BTYPE', 'Absolute error'
+ fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan)
+ fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan)
+ fits_add_checksum, extension_header, float(error_matrix)
+ mwrfits, float(error_matrix), out_file_name[2], extension_header, /no_comment, /silent
+
+ journal, 'Error-matrix extension correctly added.'
+
+ ; fix some header keywords
+
+ fix_fits_header, out_file_name[2]
+
+ ; management of stokes images
+
+ journal, 'Calibrating Stokes parameter I...'
+
+ i = reform(stokes_subdark[*, *, 0])
+ i = metis_flat_field(i, header, cal_pack)
+ i = metis_vignetting(i, header, cal_pack)
+ i = metis_rad_cal(i, header, cal_pack, /polarimetric)
+
+ journal, 'Calibrating Stokes parameter Q...'
+
+ q = reform(stokes[*, *, 1])
+ q = metis_flat_field(q, header, cal_pack)
+ q = metis_vignetting(q, header, cal_pack)
+ q = metis_rad_cal(q, header, cal_pack, /polarimetric)
+
+ journal, 'Calibrating Stokes parameter U...'
+
+ u = reform(stokes[*, *, 2])
+ u = metis_flat_field(u, header, cal_pack)
+ u = metis_vignetting(u, header, cal_pack)
+ u = metis_rad_cal(u, header, cal_pack, /polarimetric)
+
+ ; keywords specific for stokes images
+
+ primary_stokes_header = primary_header
+
+ ; name of the fits file
+
+ file_name = 'solo_L2_metis-vl-stokes_' + date_beg_string + '_V' + version + '.fits'
+ out_file_name[3] = 'output/' + file_name
+
+ fxaddpar, primary_stokes_header, 'FILENAME', file_name
+ fxaddpar, primary_stokes_header, 'BTYPE', 'Stokes I'
+ fxaddpar, primary_stokes_header, 'BUNIT', cal_pack.vl_channel.cal_units
+ fxaddpar, primary_stokes_header, 'DATAMIN', min(i, /nan)
+ fxaddpar, primary_stokes_header, 'DATAMAX', max(i, /nan)
+
+ ; add the history keyword
+
+ for k = 0, n_elements(tb_history) - 1 do fxaddpar, primary_stokes_header, 'HISTORY', tb_history[k]
+ fxaddpar, primary_stokes_header, 'HISTORY', 'L2 FITS file created on ' + date
+
+ ; add checksum and datasum to the fits header
+ if not ref_detector then i = metis_rectify(i, 'VL')
+ fits_add_checksum, primary_stokes_header, float(i)
+ mwrfits, float(i), out_file_name[3], primary_stokes_header, /no_comment, /create, /silent
+
+ journal, 'Stokes parameters FITS file created:'
+ journal, ' file name = ' + file_basename(out_file_name[3])
+
+ ; add the extension with the stokes q image
+
+ base_header = primary_stokes_header
+ sxdelpar, base_header, 'SIMPLE'
+ sxdelpar, base_header, 'EXTEND'
+ sxdelpar, base_header, 'DATASUM'
+ sxdelpar, base_header, 'CHECKSUM'
+
+ extension_header = base_header
+ fxaddpar, extension_header, 'XTENSION', 'IMAGE', 'image extension', before = 'BITPIX'
+ fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'EXTNAME', 'Stokes Q', 'extension name', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'BTYPE', 'Stokes Q'
+ fxaddpar, extension_header, 'BUNIT', cal_pack.vl_channel.cal_units
+ fxaddpar, extension_header, 'DATAMIN', min(q, /nan)
+ fxaddpar, extension_header, 'DATAMAX', max(q, /nan)
+ sxdelpar, extension_header, 'EXTEND'
+ if not ref_detector then q = metis_rectify(q, 'VL')
+ fits_add_checksum, extension_header, float(q)
+ mwrfits, float(q), out_file_name[3], extension_header, /no_comment, /silent
+
+ journal, 'Q parameter extension correctly added.'
+
+ ; add the extension with the stokes u image
+
+ extension_header = base_header
+ fxaddpar, extension_header, 'XTENSION', 'IMAGE', 'image extension', before = 'BITPIX'
+ fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'EXTNAME', 'Stokes U', 'extension name', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'BTYPE', 'Stokes U'
+ fxaddpar, extension_header, 'BUNIT', cal_pack.vl_channel.cal_units
+ fxaddpar, extension_header, 'DATAMIN', min(u, /nan)
+ fxaddpar, extension_header, 'DATAMAX', max(u, /nan)
+ sxdelpar, extension_header, 'EXTEND'
+ if not ref_detector then u = metis_rectify(u, 'VL')
+ fits_add_checksum, extension_header, float(u)
+ mwrfits, float(u), out_file_name[3], extension_header, /no_comment, /silent
+
+ journal, 'U parameter extension correctly added.'
+
+ sxdelpar, base_header, 'COMMENT'
+ sxdelpar, base_header, 'HISTORY'
+
+ ; add the extension with the quality matrix
+
+ extension_header = base_header
+ fxaddpar, extension_header, 'XTENSION', 'IMAGE', 'image extension', before = 'BITPIX'
+ fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'extension name', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'BTYPE', 'Pixel quality'
+ fxaddpar, extension_header, 'BUNIT', 'None'
+ fxaddpar, extension_header, 'DATAMIN', min(quality_matrix, /nan)
+ fxaddpar, extension_header, 'DATAMAX', max(quality_matrix, /nan)
+ fxaddpar, extension_header, 'COMMENT', 'Quality matrix values:'
+ fxaddpar, extension_header, 'COMMENT', ' NaN = saturated or null L0 pixel counts'
+ fxaddpar, extension_header, 'COMMENT', ' 0 = unreliable pixel value'
+ fxaddpar, extension_header, 'COMMENT', ' 1 = good pixel'
+ if not ref_detector then quality_matrix = metis_rectify(quality_matrix, 'VL')
+ fits_add_checksum, extension_header, float(quality_matrix)
+ mwrfits, float(quality_matrix), out_file_name[3], extension_header, /no_comment, /silent
+
+ journal, 'Quality-matrix extension correctly added.'
+
+ ; add the extension with the error matrix
+
+ if not ref_detector then error_matrix = metis_rectify(error_matrix, 'VL')
+ error_matrix = tb_image * 0. ; tb_image * sqrt(tb_error)
+ extension_header = base_header
+ fxaddpar, extension_header, 'XTENSION', 'IMAGE', 'image extension', before = 'BITPIX'
+ fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN'
+ fxaddpar, extension_header, 'BTYPE', 'Absolute error'
+ fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan)
+ fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan)
+ fits_add_checksum, extension_header, float(error_matrix)
+ mwrfits, float(error_matrix), out_file_name[3], extension_header, /no_comment, /silent
+
+ journal, 'Error-matrix extension correctly added.'
+
+ ; fix some header keywords
+
+ fix_fits_header, out_file_name[3]
+
+ ; write the auxiliary information file
+
+ output = { $
+ file_name: out_file_name, $
+ l1_file_name: input.file_name, $
+ log_file_name: 'output/metis_l2_prep_log.txt' $
+ }
+
+ json_write, output, 'output/contents.json'
+
+ ; unload the spice kernels
+
+ load_spice_kernels, kernel_list = kernel_list, /unload
+
+ journal, 'SPICE kernel files unloaded.'
+
+ ; close the log
+
+ journal, 'Exiting without errors.'
+ journal
+
+ exit, status = 0
+end
\ No newline at end of file