Select Git revision
metis_l2_prep_vl_polariz.pro
metis_l2_prep_vl_polariz.pro 22.71 KiB
pro metis_l2_prep_vl_polariz
; keyword defining if the detector reference frame must be used for the output
ref_detector = 1
; start the log
journal,'output/metis_l2_prep_log.txt'
; read the auxiliary file - here we have all the inputs we need
input = json_parse('input/contents.json', /toarray, /tostruct)
; load the spice kernels
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
; read the calibration package
cal_pack = json_parse(input.cal_pack_path + '/index.json', /toarray, /tostruct)
; include the calibration package path into the cal_pack structure
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)')
n = n_elements(input.file_name)
; 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
; calibration block
data = !null
data_header = !null
data_subdark = !null
quality_matrix = 1.
for k = 0, 3 do begin
; read the input image
image = mrdfits(input.file_name[k], 0, primary_header, /silent)
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
; ====================================
; 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.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 pmp raw voltages (dacpol)
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
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)')
; pile up images and headers
data = [[[data]], [[image]]]
data_header = [data_header, header]
; read and update the quality matrix
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, 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]]]
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)
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
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
; 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]
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
; compute the tb from the dark-subtracted stokes i and apply other calibrations
journal, 'Calibrating total brightness...'
tb_image = reform(stokes_subdark[*, *, 0])
tb_image = metis_flat_field(tb_image, header, cal_pack, history = tb_history)
tb_image = metis_vignetting(tb_image, header, cal_pack, history = tb_history)
tb_image = metis_rad_cal(tb_image, header, 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(reform(stokes[*, *, 1])^2 + reform(stokes[*, *, 2])^2)
pb_image = metis_flat_field(pb_image, header, cal_pack, history = pb_history)
pb_image = metis_vignetting(pb_image, header, cal_pack, history = pb_history)
pb_image = metis_rad_cal(pb_image, header, 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]) * !radeg
journal, 'Polarization angle correctly computed.'
; definitions for the primary header
; version of the fits file
version = string(input.l2_version + 1, 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, 'ORIGIN', ''
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, after = 'VERS_SW'
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
; 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.'
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
extension_header = !null
fxaddpar, extension_header, 'PCOUNT', 0, 'Parameter count'
fxaddpar, extension_header, 'GCOUNT', 1, 'Group count'
fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'Extension name'
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[0], extension_header, /no_comment, /silent
journal, 'Quality-matrix extension correctly added.'
; add the extension with the error matrix
extension_header = !null
fxaddpar, extension_header, 'PCOUNT', 0, 'Parameter count'
fxaddpar, extension_header, 'GCOUNT', 1, 'Group count'
fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'Extension name'
error_matrix = intarr(header.naxis1, header.naxis2)
if not ref_detector then error_matrix = metis_rectify(error_matrix, 'VL')
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
extension_header = !null
fxaddpar, extension_header, 'PCOUNT', 0, 'Parameter count'
fxaddpar, extension_header, 'GCOUNT', 1, 'Group count'
fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'Extension name'
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[1], extension_header, /no_comment, /silent
journal, 'Quality-matrix extension correctly added.'
; add the extension with the error matrix
extension_header = !null
fxaddpar, extension_header, 'PCOUNT', 0, 'Parameter count'
fxaddpar, extension_header, 'GCOUNT', 1, 'Group count'
fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'Extension name'
error_matrix = intarr(header.naxis1, header.naxis2)
if not ref_detector then error_matrix = metis_rectify(error_matrix, 'VL')
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
journal, 'Polarization-angle FITS file created:'
journal, ' file name = ' + file_basename(out_file_name[2])
; add the extension with the quality matrix
extension_header = !null fxaddpar, extension_header, 'PCOUNT', 0, 'Parameter count'
fxaddpar, extension_header, 'GCOUNT', 1, 'Group count'
fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'Extension name'
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[2], extension_header, /no_comment, /silent
journal, 'Quality-matrix extension correctly added.'
; add the extension with the error matrix
extension_header = !null
fxaddpar, extension_header, 'PCOUNT', 0, 'Parameter count'
fxaddpar, extension_header, 'GCOUNT', 1, 'Group count'
fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'Extension name'
error_matrix = intarr(header.naxis1, header.naxis2)
if not ref_detector then error_matrix = metis_rectify(error_matrix, 'VL')
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 = !null
fxaddpar, extension_header, 'PCOUNT', 0, 'Parameter count'
fxaddpar, extension_header, 'GCOUNT', 1, 'Group count'
fxaddpar, extension_header, 'EXTNAME', 'Stokes Q', 'Extension name'
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)
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 = !null
fxaddpar, extension_header, 'PCOUNT', 0, 'Parameter count'
fxaddpar, extension_header, 'GCOUNT', 1, 'Group count'
fxaddpar, extension_header, 'EXTNAME', 'Stokes U', 'Extension name'
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)
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
extension_header = !null
fxaddpar, extension_header, 'PCOUNT', 0, 'Parameter count'
fxaddpar, extension_header, 'GCOUNT', 1, 'Group count'
fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'Extension name'
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
extension_header = !null
fxaddpar, extension_header, 'PCOUNT', 0, 'Parameter count'
fxaddpar, extension_header, 'GCOUNT', 1, 'Group count'
fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'Extension name'
error_matrix = intarr(header.naxis1, header.naxis2)
if not ref_detector then error_matrix = metis_rectify(error_matrix, 'VL')
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