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