Version 0.0

decode_obt.pro

+function decode_obt, t_coarse, t_fine, to_decimal = to_decimal, from_decimal = from_decimal
+	if keyword_set(to_decimal) then return, double(t_coarse) + double(t_fine)/65536.D0
+	if keyword_set(from_decimal) then obt = [ulong(t_coarse), uint((t_coarse - ulong(t_coarse)) * 65536)] else obt = [t_coarse, t_fine]
+	return, (obt.tostring()).join(':')
+end

fits_hdr2struct.pro

+function fits_hdr2struct, string_hdr
+
+	struct = !null
+	n = n_elements(string_hdr)
+	for i = 0, n - 1 do begin
+		tag = string_hdr[i].substring(0, 7)
+		tag = tag.trim()
+		if tag eq '' then continue
+		if tag.startswith('end', /fold) or tag.contains('continue', /fold) then continue
+		if isa(struct) then if where(tag_names(struct) eq tag) ge 0 then continue
+		struct = create_struct(struct, tag.replace('-', ' '), fxpar(string_hdr, tag))
+	endfor
+
+	return, struct
+end

get_light_time.pro

+function get_light_time, utc, body, target, radvel = radvel
+
+	; convert the requested date into ephemeris time
+
+	cspice_str2et, utc, et
+
+	; coordinates of the body wrt the target in the j2000 system
+
+	cspice_spkezr, body, et, 'J2000', 'NONE', target, state, ltime
+
+	; radial component of the velocity
+
+	radvel = 1000.d0 * total(state[0 : 2] * state[3 : 5]) / sqrt(total(state[0 : 2]^2))
+
+	return, ltime
+end

json_write.pro

+pro json_write, struct, filename
+	openw, unit, filename, /get_lun
+	printf, unit, struct, /implied_print
+	free_lun, unit
+end

load_spice_kernels.pro

+pro load_spice_kernels, path, unload = unload, kernel_list = kernel_list, kernel_version = kernel_version
+	if keyword_set(unload) then cspice_kclear else begin
+		metakernel = path + '/mk/solo_ANC_soc-flown-mk.tm'
+		openr, in, metakernel, /get_lun
+		kernel_list = list()
+		while ~ eof(in) do begin
+			line = ''
+			readf, in, line
+			if line.contains('$KERNELS', /fold) then begin
+				line = line.replace('$KERNELS', path)
+				line = line.replace("'", '')
+				kernel_list.add, line.trim()
+			endif
+			if line.contains('SKD_VERSION', /fold) then begin
+				fields = stregex(line, ".*'([a-z_0-9]*)'", /extract, /sub)
+				kernel_version = fields[1]
+			endif
+		endwhile
+		close, /all
+		foreach item, kernel_list do cspice_furnsh, item
+	endelse
+end

metis_dark_uvda.pro

+function metis_dark_uvda, data, header, cal_pack, history = history
+
+	dark = cal_pack.uv_channel.dark
+
+	dit = header.dit/1000.d0
+	ndit1 = header.ndit1
+	ndit2 = header.ndit2
+	nbin = sqrt(header.nbin)
+	tsensor = header.tsensor
+	obj_cnt = header.obj_cnt
+	obt_beg = header.obt_beg
+
+	for i = 0, n_elements(dark) - 1 do begin
+		if (dark[i].dit eq dit) and (dark[i].nbin eq nbin) and (dark[i].ndit1 eq ndit1) and (dark[i].ndit2 eq ndit2) and (abs(dark[i].tsensor - tsensor) lt 5) and (abs(obt_beg - dark[i].obt_beg) lt (3600 * 24 * 10)) then begin
+			if obj_cnt eq 1 then begin
+				dark_file = dark[i].file_name.cnt_1
+				dark_image = float(readfits(cal_pack.path + dark[i].file_name.cnt_1, /silent))
+				dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin^2
+				mask = where(data eq 0., count)
+				data = data - dark_image * ndit1 * ndit2
+				data[mask] = 0.
+
+				goto, jump
+			endif else begin
+				dark_file = dark[i].file_name.cnt_2
+				dark_image = float(readfits(cal_pack.path + dark[i].file_name.cnt_2, /silent))
+				dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin^2
+				mask = where(data eq 0., count)
+				data = data - dark_image * ndit1 * ndit2
+				data[mask] = 0.
+
+				goto, jump
+			endelse
+		endif
+	endfor
+
+	obt_available = !null
+	for i = 0, n_elements(dark) - 1 do begin
+		if (dark[i].dit eq dit) and (abs(dark[i].tsensor - tsensor) lt 5) then begin
+			obt_available = [obt_available, dark[i].obt_beg]
+		endif
+	endfor
+
+	delta_obt = min(abs(obt_beg - obt_available), j)
+	i = where(dark.obt_beg eq obt_available[j])
+
+	if obj_cnt eq 1 then begin
+		if (ndit1 eq dark[i].ndit1) and (ndit2 eq dark[i].ndit2) then begin
+			dark_image = float(readfits(cal_pack.path + dark[i].file_name.cnt_1, /silent))
+			dark_file = dark[i].file_name.cnt_1
+		endif else begin
+			q = ndit2/dark[i].ndit2
+			dark_1 = float(readfits(cal_pack.path + dark[i].file_name.cnt_1, /silent))
+			dark_2 = float(readfits(cal_pack.path + dark[i].file_name.cnt_2, /silent))
+			dark_image = (dark_1 + (q - 1) * dark_2)/q
+			dark_file = dark[i].file_name.cnt_1 + ' - ' + dark[i].file_name.cnt_2
+		endelse
+	endif else begin
+		dark_image = float(readfits(cal_pack.path + dark[i].file_name.cnt_2, /silent))
+		dark_file = dark[i].file_name.cnt_2
+	endelse
+
+	dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin^2
+
+	mask = where(data eq 0., count)
+	data = data - dark_image * ndit1 * ndit2
+	data[mask] = 0.
+
+	jump:
+
+	if ~ isa(history) then history = !null
+	history = [history, 'Dark correction: ', '  ' + dark_file]
+
+	return, data
+end

metis_dark_vlda.pro

+function metis_dark_vlda, data, header, cal_pack, history = history
+
+	bias = cal_pack.vl_channel.bias
+	dark = cal_pack.vl_channel.dark
+	bias_dark = cal_pack.vl_channel.bias_dark
+
+	dit = header.dit
+	ndit = header.ndit
+	nbin = sqrt(header.nbin)
+	tsensor = header.tsensor
+	obt_beg = header.obt_beg/1000d0
+
+	for i = 0, n_elements(bias_dark) - 1 do begin
+		if bias_dark[i].dit eq dit/1000. and bias_dark[i].nbin eq nbin and abs(bias_dark[i].tsensor - tsensor) lt 5 then begin
+			dark_image = float(readfits(cal_pack.path + bias_dark[i].file_name, /silent))
+			dark_file = bias_dark[i].file_name
+
+			mask = where(data eq 0.)
+			dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin^2 * ndit
+			data = data - dark_image
+			data[mask] = 0.
+
+			goto, jump
+		endif
+	endfor
+
+	bias_obt_available = !null
+	for i = 0, n_elements(bias) - 1 do begin
+		if abs(bias[i].tsensor - tsensor) lt 5 then $
+		bias_obt_available = [bias_obt_available, bias[i].obt_beg]
+	endfor
+
+	delta_obt = min(abs(obt_beg - bias_obt_available), j)
+	i = where(bias.obt_beg eq bias_obt_available[j])
+	bias_image = float(readfits(cal_pack.path + bias[i].file_name, /silent))
+	dark_file = bias[i].file_name
+
+	dark_obt_available = !null
+	for i = 0, n_elements(dark) - 1 do begin
+		if abs(dark[i].tsensor - tsensor) lt 5 then $
+		dark_obt_available = [dark_obt_available, dark[i].obt_beg]
+	endfor
+
+	delta_obt = min(abs(obt_beg - dark_obt_available), j)
+	i = where(dark.obt_beg eq dark_obt_available[j])
+	dark_image = float(readfits(cal_pack.path + dark[i].file_name, /silent))
+	dark_file = [dark_file, dark[i].file_name]
+
+	bias_image = rebin(bias_image, header.naxis1, header.naxis2) * nbin^2
+	dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin^2
+
+	mask = where(data eq 0.)
+	data = data - (bias_image * ndit + header.xposure * dark_image)
+	data[mask] = 0.
+
+	jump:
+
+	if ~ isa(history) then history = !null
+	history = [history, 'Bias and dark-current corrections: ', '  ' + dark_file]
+
+	return, data
+end

metis_flat_field.pro

+function metis_flat_field, data, header, cal_pack, history = history
+
+	if header.filter.contains('VL', /fold) then begin
+		ff_info = cal_pack.vl_channel.flat_field
+
+		ff_image = readfits(cal_pack.path + ff_info.file_name, /silent)
+		ff_file = ff_info.file_name
+
+		first_row = header.ref_rows eq 'Included' ? 2 : 0
+		detector_size = cal_pack.vl_channel.detector_size.value
+		ff_image = ff_image[*, first_row : first_row + detector_size - 1]
+	endif
+
+	if header.filter.contains('UV', /fold) then begin
+		ff_info = cal_pack.uv_channel.flat_field
+
+		nbin = sqrt(header.nbin)
+		i = where(ff_info.nbin eq nbin)
+		ff_image = readfits(cal_pack.path + ff_info[i].file_name, /silent)
+		ff_file = ff_info[i].file_name
+	endif
+
+	ff_image = rebin(ff_image, header.naxis1, header.naxis2)
+	
+	mask = where(ff_image le 0.)
+	ff_image[mask] = 1.
+	data = data / ff_image
+	data[mask] = 0.
+
+	if ~ isa(history) then history = !null
+	history = [history, 'Flat-field correction: ', '  ' + ff_file]
+
+	return, data
+end

metis_l2_prep_uv.pro

+pro metis_l2_prep_uv
+
+	; 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
+
+	; read the calibration package index
+
+	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)
+
+	; read the primary hdu
+
+	data = mrdfits(input.file_name, 0, primary_header, /silent)
+
+	; convert the string header into an idl structure
+
+	header = fits_hdr2struct(primary_header)
+
+	; read the house-keeping extension
+
+	hk_table = mrdfits(input.file_name, 'house-keeping', hk_extension_header, /silent)
+
+	; calibration block
+
+	if header.datatype eq 1 then 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
+		; ====================================
+
+		data = metis_dark_uvda(data, header, cal_pack, history = history)
+
+		; ====================================
+		; 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, history = history)
+		data = metis_vignetting(data, header, cal_pack, history = history)
+		data = metis_rad_cal(data, header, cal_pack, history = history)
+
+		; ====================================
+		; for simple radiometric calibration
+		; data = data/header.xposure
+		; cal_pack.uv_channel.cal_units = 'DN/s'
+		; history = [history, 'Radiometric calibration: ', '  exposure normalisation only']
+		; ====================================
+
+		; ====================================
+		; temporary correction based on alfa-leo data
+		; data = data * 3.4
+		; history = [history, '  additional cal. factor from stars = 3.4']
+		; ====================================
+
+		if not ref_detector then data = metis_rectify(data, header.filter)
+	endif
+
+	if header.datatype eq 4 then begin
+		; calibration of temporal noise images
+	endif
+
+	if header.datatype eq 6 then begin
+		; calibration of cr/sep log matrices
+	endif
+
+	; definitions for the primary header
+	; version of the fits file
+
+	version = string(input.l2_version + 1, 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-draft'
+	fxaddpar, primary_header, 'ORIGIN', ''
+	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, after = 'VERS_SW'
+	fxaddpar, primary_header, 'BTYPE', 'UV Lyman-alpha intensity'
+	fxaddpar, primary_header, 'BUNIT', cal_pack.uv_channel.cal_units
+	fxaddpar, primary_header, 'DATAMIN', min(data, /nan)
+	fxaddpar, primary_header, 'DATAMAX', max(data, /nan)
+	fxaddpar, primary_header, 'XPOSURE', header.xposure
+	fxaddpar, primary_header, 'NSUMEXP', header.nsumexp
+
+	; 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
+
+	ephemerides = solo_get_ephemerides(header.date_avg, constants = cal_pack.constants)
+
+	foreach element, ephemerides do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
+
+	history = [history, 'Solar ephemerides and WCS:', '  SKD version = ' + kernel_version]
+
+	; delete useless keywords
+
+	sxdelpar, primary_header, 'BLANK'
+
+	; add the history keyword
+
+	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
+
+	fits_add_checksum, primary_header, data
+
+	; add keywords for file history
+
+	data = float(data)
+	fits_add_checksum, primary_header, data
+	mwrfits, data, out_file_name, primary_header, /no_comment, /create, /silent
+
+	extension_header = !null
+	fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'Extension name'
+	fits_add_checksum, extension_header, intarr(header.naxis1, header.naxis2)
+	mwrfits, intarr(header.naxis1, header.naxis2), out_file_name, extension_header, /no_comment, /silent
+
+	extension_header = !null
+	fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'Extension name'
+	fits_add_checksum, extension_header, intarr(header.naxis1, header.naxis2)
+	mwrfits, intarr(header.naxis1, header.naxis2), out_file_name, extension_header, /no_comment, /silent
+
+	; write the auxiliary information file
+
+	output = { $
+		file_name: out_file_name, $
+		l2_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
+
+	; close the log
+
+	journal
+	;exit, status = 0
+	return
+
+	error_handling:
+	journal, 'Errors occurred while processing.'
+	journal
+	exit, status = 1
+end

metis_l2_prep_vl_generic.pro

+pro metis_l2_prep_vl_generic
+
+	; 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
+
+	; 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)
+
+	; read the primary hdu
+
+	data = mrdfits(input.file_name, 0, primary_header, /silent)
+
+	; convert the string header into an idl structure
+
+	header = fits_hdr2struct(primary_header)
+
+	; read the house-keeping extension
+
+	hk_table = mrdfits(input.file_name, 'house-keeping', hk_extension_header, /silent)
+
+	; calibration block
+
+	if header.datatype eq 0 then 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, 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]
+		; ====================================
+
+		data = metis_flat_field(data, header, cal_pack, history = history)
+		data = metis_vignetting(data, header, cal_pack, history = history)
+
+		if header.pol_id ge 1 and header.pol_id le 4 then polarimetric = 1 else polarimetric = 0
+		data = metis_rad_cal(data, header, cal_pack, polarimetric = polarimetric, history = history)
+
+		if not ref_detector then data = rotate(data, 1)
+	endif
+
+	if header.datatype eq 3 then begin
+		; calibration of temporal noise images
+	endif
+
+	if header.datatype eq 5 then begin
+		; calibration of cr/sep log matrices
+	endif
+
+	; definitions for the primary header
+	; version of the fits file
+
+	version = string(input.l2_version + 1, 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-draft'
+	fxaddpar, primary_header, 'ORIGIN', ''
+	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, after = 'VERS_SW'
+	case header.pol_id of
+		0: fxaddpar, primary_header, 'BTYPE', 'VL fixed-polarization intensity'
+		1: fxaddpar, primary_header, 'BTYPE', 'VL fixed-polarization intensity'
+		2: fxaddpar, primary_header, 'BTYPE', 'VL fixed-polarization intensity'
+		3: fxaddpar, primary_header, 'BTYPE', 'VL fixed-polarization intensity'
+		4: fxaddpar, primary_header, 'BTYPE', 'VL fixed-polarization intensity'
+		5: fxaddpar, primary_header, 'BTYPE', 'VL total brightness'
+		6: fxaddpar, primary_header, 'BTYPE', 'VL total brightness'
+	endcase
+	fxaddpar, primary_header, 'BUNIT', cal_pack.vl_channel.cal_units
+	fxaddpar, primary_header, 'DATAMIN', min(data, /nan)
+	fxaddpar, primary_header, 'DATAMAX', max(data, /nan)
+	fxaddpar, primary_header, 'XPOSURE', header.xposure
+	fxaddpar, primary_header, 'NSUMEXP', header.nsumexp
+
+	; 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
+
+	ephemerides = solo_get_ephemerides(header.date_avg, constants = cal_pack.constants)
+
+	foreach element, ephemerides do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
+
+	history = [history, 'Solar ephemerides and WCS:', '  SKD version = ' + kernel_version]
+
+	; delete useless keywords
+
+	sxdelpar, primary_header, 'BLANK'
+
+	; add the history keyword
+
+	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
+
+	data = float(data)
+	fits_add_checksum, primary_header, data
+	mwrfits, data, out_file_name, primary_header, /no_comment, /create, /silent
+
+	extension_header = !null
+	fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'Extension name'
+	fits_add_checksum, extension_header, intarr(header.naxis1, header.naxis2)
+	mwrfits, intarr(header.naxis1, header.naxis2), out_file_name, extension_header, /no_comment, /silent
+
+	extension_header = !null
+	fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'Extension name'
+	fits_add_checksum, extension_header, intarr(header.naxis1, header.naxis2)
+	mwrfits, intarr(header.naxis1, header.naxis2), out_file_name, extension_header, /no_comment, /silent
+
+	; write the auxiliary information file
+
+	output = { $
+		file_name: out_file_name, $
+		l2_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
+
+	; close the log
+
+	journal
+	;exit, status = 0
+	return
+
+	error_handling:
+	journal, 'Errors occurred while processing.', /continue
+	journal
+	; exit, status = 1
+end

metis_l2_prep_vl_polariz.pro

+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
+
+	; 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)
+
+	n = n_elements(input.file_name)
+
+	; check on the number of input images (= npol = 4)
+
+	if n ne 4 then begin
+		error_message = 'number of input images is not equal to 4'
+		goto, error_handling
+	endif
+
+	; calibration block
+
+	data = !null
+	data_header = !null
+	data_subdark = !null
+	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 consistency of polarization state
+
+		if header.pol_id lt 1 or header.pol_id gt 4 then begin
+			error_message = 'image no. ' + string(k, format = '(I0)') + $
+			' has inconsistent polarization state'
+			goto, error_handling
+		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
+			error_message = 'image no. ' + string(k, format = '(I0)') + ' has inconsistent PMP voltages'
+			goto, error_handling
+		endif
+
+		; pile up images and headers
+
+		data = [[[data]], [[image]]]
+		data_header = [data_header, header]
+
+		; 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('VL subdark/' + file, 0, /silent)
+		; image_subdark = rebin(image_subdark, header.naxis1, header.naxis2)
+		; image_subdark = image_subdark * header.nbin * 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)
+
+	; 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
+
+	demod_tensor = fltarr(header.naxis1, header.naxis2, 4)
+	stokes = dblarr(header.naxis1, header.naxis2, 3)
+	stokes_subdark = dblarr(header.naxis1, header.naxis2, 3)
+	stokes_name = ['I', 'Q', 'U']
+	for i = 0, 2 do begin
+		for j = 0, 3 do begin
+
+			; 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]
+
+			n = where(demod_info.angle eq angle[0] and demod_info.stokes eq stokes_name[i])
+			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
+		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
+
+	; compute the tb from the dark-subtracted stokes i and apply other calibrations
+
+	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)
+
+	if not ref_detector then tb_image = rotate(tb_image, 1)
+
+	; compute the pb from the stokes q and u and apply other calibrations
+
+	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)
+
+	if not ref_detector then pb_image = rotate(pb_image, 1)
+
+	; compute the polarization angle from the stokes q and u
+
+	pol_angle = 0.5d0 * atan(stokes[*, *, 2], stokes[*, *, 1]) * !radeg
+
+	if not ref_detector then pol_angle = rotate(pol_angle, 1)
+
+	; 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-draft'
+	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, '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
+
+	ephemerides = solo_get_ephemerides(header.date_avg, constants = cal_pack.constants)
+
+	foreach element, ephemerides do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
+
+	history = ['Solar ephemerides and WCS:', '  SKD version = ' + kernel_version]
+	tb_history = [tb_history, history]
+	pb_history = [pb_history, history]
+
+	; delete useless keywords
+
+	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(':', '')
+
+	; 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 = '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)
+	sxdelpar, primary_pb_header, 'BLANK'
+
+	; 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
+
+	fits_add_checksum, primary_pb_header, pb_image
+
+	mwrfits, pb_image, out_file_name, primary_pb_header, /no_comment, /create, /silent
+
+	extension_header = !null
+	fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'Extension name'
+	fits_add_checksum, extension_header, intarr(header.naxis1, header.naxis2)
+	mwrfits, intarr(header.naxis1, header.naxis2), out_file_name, extension_header, /no_comment, /silent
+
+	extension_header = !null
+	fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'Extension name'
+	fits_add_checksum, extension_header, intarr(header.naxis1, header.naxis2)
+	mwrfits, intarr(header.naxis1, header.naxis2), out_file_name, extension_header, /no_comment, /silent
+
+	; 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 = [out_file_name, '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)
+	sxdelpar, primary_tb_header, 'BLANK'
+
+	; 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
+
+	fits_add_checksum, primary_tb_header, tb_image
+
+	mwrfits, tb_image, out_file_name[1], primary_tb_header, /no_comment, /create, /silent
+
+	extension_header = !null
+	fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'Extension name'
+	fits_add_checksum, extension_header, intarr(header.naxis1, header.naxis2)
+	mwrfits, intarr(header.naxis1, header.naxis2), out_file_name[1], extension_header, /no_comment, /silent
+
+	extension_header = !null
+	fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'Extension name'
+	fits_add_checksum, extension_header, intarr(header.naxis1, header.naxis2)
+	mwrfits, intarr(header.naxis1, header.naxis2), out_file_name[1], extension_header, /no_comment, /silent
+
+	; 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 = '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)
+	sxdelpar, primary_polangle_header, 'BLANK'
+
+	; 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
+
+	fits_add_checksum, primary_polangle_header, pol_angle
+
+	mwrfits, pol_angle, out_file_name, primary_polangle_header, /no_comment, /create, /silent
+
+	extension_header = !null
+	fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'Extension name'
+	fits_add_checksum, extension_header, intarr(header.naxis1, header.naxis2)
+	mwrfits, intarr(header.naxis1, header.naxis2), out_file_name, extension_header, /no_comment, /silent
+
+	extension_header = !null
+	fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'Extension name'
+	fits_add_checksum, extension_header, intarr(header.naxis1, header.naxis2)
+	mwrfits, intarr(header.naxis1, header.naxis2), out_file_name, extension_header, /no_comment, /silent
+
+	; =================== STOKES PARAMETERS ======================
+
+	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)
+	if not ref_detector then i = rotate(i, 1)
+
+	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)
+	if not ref_detector then q = rotate(q, 1)
+
+	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)
+	if not ref_detector then u = rotate(u, 1)
+
+	; 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 = '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)
+	sxdelpar, primary_stokes_header, 'BLANK'
+
+	; 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
+
+	fits_add_checksum, primary_stokes_header, i
+	mwrfits, i, out_file_name, primary_stokes_header, /no_comment, /create, /silent
+
+	extension_header = !null
+	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)
+	fits_add_checksum, extension_header, q
+	mwrfits, q, out_file_name, extension_header, /no_comment, /silent
+
+	extension_header = !null
+	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)
+	fits_add_checksum, extension_header, u
+	mwrfits, u, out_file_name, extension_header, /no_comment, /silent
+
+	extension_header = !null
+	fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'Extension name'
+	mwrfits, intarr(header.naxis1, header.naxis2), out_file_name, extension_header, /no_comment, /silent
+
+	extension_header = !null
+	fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'Extension name'
+	mwrfits, intarr(header.naxis1, header.naxis2), out_file_name, extension_header, /no_comment, /silent
+
+	; write the auxiliary information file
+
+	output = { $
+		file_name: out_file_name, $
+		l2_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
+
+	; close the log
+
+	journal
+	;exit, status = 0
+	return
+
+	error_handling:
+	journal, 'Errors occurred while processing:'
+	journal, '  ' + error_message
+	journal
+	exit, status = 1
+end

metis_rad_cal.pro

+function metis_rad_cal, data, header, cal_pack, polarimetric = polarimetric, history = history
+
+	if header.filter.contains('VL', /fold) then begin
+
+		; NOTE - this is to take into account that calibration of polarimetric images or their combination, must not include the pmp efficiency (0.5) since this is included in the demodulation formulae
+
+		if keyword_set(polarimetric) then pmp_factor = 2.d0 else pmp_factor = 1.d0
+
+		angular_pixel = cal_pack.vl_channel.angular_pixel.value * header.nbin
+		rad_info = cal_pack.vl_channel.radiometry[0]
+		units = cal_pack.vl_channel.cal_units
+		unit_factor = rad_info.msb.value
+	end
+
+	if header.filter.contains('UV', /fold) then begin
+		angular_pixel = cal_pack.uv_channel.angular_pixel.value * header.nbin
+		rad_info = cal_pack.uv_channel.radiometry[0]
+		units = cal_pack.uv_channel.cal_units
+		pmp_factor = 1.d0
+		unit_factor = 1.
+	endif
+
+	rad_factor = rad_info.rad_response.value * pmp_factor * $
+		cal_pack.instrument.pupil_area.value * $
+		angular_pixel * $
+		unit_factor
+
+	cal_factor = 1. / rad_factor / header.xposure
+
+	data = data * cal_factor
+
+	if ~ isa(history) then history = !null
+	history = [history, 'Radiometric calibration:', '  cal. factor = ' + string(cal_factor, format = '(E8.2)') + ' ' + units + '/DN']
+
+	return, data
+end

metis_rectify.pro

+function metis_rectify, image, filter
+	if filter.toupper() eq 'VL' then image = rotate(image, 1)
+	if filter.toupper() eq 'UV' then begin
+		image = reverse(image, 1)
+		image = rotate(image, 1)
+	endif
+	return, image
+end

metis_vignetting.pro

+function metis_vignetting, data, header, cal_pack, history = history
+
+	if header.filter.contains('VL', /fold) then vig_info = cal_pack.vl_channel.vignetting
+	if header.filter.contains('UV', /fold) then vig_info = cal_pack.uv_channel.vignetting
+
+	vig_image = readfits(cal_pack.path + vig_info.file_name, /silent)
+	vig_image = median(vig_image, 5)
+	vig_file = vig_info.file_name
+
+	dx = 0.
+	dy = 0.
+
+	if header.filter.contains('UV', /fold) then begin
+		vl_boresight = cal_pack.vl_channel.boresight[0]
+		uv_boresight = cal_pack.uv_channel.boresight[0]
+
+		x1 = uv_boresight.borpix1.value
+		x2 = vl_boresight.borpix1.value
+		y1 = uv_boresight.borpix2.value
+		y2 = vl_boresight.borpix2.value
+
+		plate_scale = vl_boresight.plate_scale.value
+		dx = (x1 - x2)/plate_scale
+		dy = (y1 - y2)/plate_scale
+
+		vig_image = rotate(vig_image, 1)
+		vig_image = shift(vig_image, dx, dy)
+		vig_image[2047 - abs(dx) : 2047, *] = 0.
+		vig_image[*, 0 : abs(dy) - 1] = 0.
+		vig_image = rotate(vig_image, 3)
+		vig_image = reverse(vig_image, 1)
+	endif
+
+	vig_image = rebin(vig_image, header.naxis1, header.naxis2)
+	vig_image = (vig_image > 0.) < 1.
+
+	mask = where(vig_image eq 0.)
+	vig_image[mask] = 1.
+
+	data = data / vig_image
+	data[mask] = 0.
+
+	if ~ isa(history) then history = !null
+	history = [history, 'Vignetting correction: ', '  ' + vig_file + ' shifted by [' + string(dx, format = '(f0.1)') + ', ' + string(dy, format = '(f0.1)') + '] pixel']
+
+	return, data
+end

metis_wcs.pro

+function metis_wcs, header, cal_pack, ref_detector = ref_detector
+
+	if header.filter.contains('VL', /fold) then $
+		channel = cal_pack.vl_channel else $
+		channel = cal_pack.uv_channel
+
+	boresight = channel.boresight[0]
+	detector_size = channel.detector_size.value
+
+	; compute correct boresight parameters to account for image orientation, binning, and fits convention
+
+	; plate scale
+
+	cdelt1 = boresight.plate_scale.value * header.nbin1
+	cdelt2 = boresight.plate_scale.value * header.nbin2
+	cdelt = [cdelt1, cdelt2]
+
+	; old definitions if boresight parameters in pixel units are used
+
+	; pntpix1 = boresight.pntpix1.value + (nominal_size - 1)/2.
+	; pntpix1 = (pntpix1 + 0.5)/header.nbin1 + 0.5
+	; pntpix2 = boresight.pntpix2.value + (nominal_size - 1)/2.
+	; pntpix2 = (pntpix2 + 0.5)/header.nbin2 + 0.5
+
+	; new definitions if boresight parameters in arcsec units are used
+
+	detector_size = detector_size/sqrt(header.nbin)
+	xcen = (detector_size + 1)/2.
+	ycen = (detector_size + 1)/2.
+
+	; boresight, i.e., io center
+
+	borpix1 = boresight.borpix1.value/cdelt1 + xcen
+	borpix2 = boresight.borpix2.value/cdelt2 + ycen
+	borpix = [borpix1, borpix2]
+
+	; s/c pointing
+
+	pntpix1 = boresight.pntpix1.value/cdelt1 + xcen
+	pntpix2 = boresight.pntpix2.value/cdelt2 + ycen
+	pntpix = [pntpix1, pntpix2]
+
+	; determine spacecract pointing information in the hpc reference frame using the spice kernels
+
+	pointing = solo_get_pointing(header.date_avg)
+
+	; NOTE - values are defined as follows:
+	; pointing[0] = yaw (arcsec)
+	; pointing[1] = pitch (arcsec)
+	; pointing[2] = roll (deg)
+
+	; physical coordinates of s/c pointing (i.e., yaw and pitch) in the hpc reference frame
+
+	pntval1 = pointing[0]
+	pntval2 = pointing[1]
+	pntval = [pntval1, pntval2]
+
+	; correct the roll angle value for metis misalignment
+
+	roll = (pointing[2] + boresight.delta_roll.value) * !dtor
+
+	; wcs rotation matrix in the hpc reference frame
+
+	pc = [[cos(roll), -sin(roll)], [sin(roll), cos(roll)]]
+
+	; if requested, transform the wcs matrix to the detector reference frame and adjust the boresight and spacecraft pointing parameters
+
+	if keyword_set(ref_detector) then begin
+		ctype1 = 'HPLT-TAN'
+		ctype2 = 'HPLN-TAN'
+		if header.filter.contains('UV', /fold) then begin
+			borpix_prime = borpix
+			borpix[0] = detector_size - (borpix_prime[1] - 1.)
+			borpix[1] = detector_size - (borpix_prime[0] - 1.)
+			pntpix_prime = pntpix
+			pntpix[0] = detector_size - (pntpix_prime[1] - 1.)
+			pntpix[1] = detector_size - (pntpix_prime[0] - 1.)
+			pc = -reverse(pc, 1)
+			pc = transpose(pc)
+		endif
+		if header.filter.contains('VL', /fold) then begin
+			borpix_prime = borpix
+			borpix[0] = borpix_prime[1]
+			borpix[1] = detector_size - (borpix_prime[0] - 1.)
+			pntpix_prime = pntpix
+			pntpix[0] = pntpix_prime[1]
+			pntpix[1] = detector_size - (pntpix_prime[0] - 1.)
+			roll = roll + !dpi/2.
+			pc = [[cos(roll), -sin(roll)], [sin(roll), cos(roll)]]
+		endif
+	endif else begin
+		ctype1 = 'HPLN-TAN'
+		ctype2 = 'HPLT-TAN'
+	endelse
+
+	; get sun's center pixel
+
+	sunval = [0., 0.]
+	sunpix = (invert(pc) ## (sunval - pntval))/cdelt + pntpix
+
+	; get coordinates of the image center pixel
+
+	crpix = [xcen, ycen]
+	crval = (pc ## (crpix - pntpix)) * cdelt + pntval
+
+	; create the wcs list
+
+	wcs = list()
+
+	wcs.add, { $
+		name: 'WCSNAME', $
+		value: 'Helioprojective-Cartesian', $
+		comment: 'Name of coordinate system'}
+	wcs.add, { $
+		name: 'CTYPE1', $
+		value: ctype1, $
+		comment: ctype1 eq 'HPLT-TAN' ? 'Helioprojective latitude (Solar Y)' : 'Helioprojective longitude (Solar X)'}
+	wcs.add, { $
+		name: 'CTYPE2', $
+		value: ctype2, $
+		comment: ctype2 eq 'HPLT-TAN' ? 'Helioprojective latitude (Solar Y)' : 'Helioprojective longitude (Solar X)'}
+	wcs.add, { $
+		name: 'CUNIT1', $
+		value: 'arcsec', $
+		comment: 'Units along axis 1'}
+	wcs.add, { $
+		name: 'CUNIT2', $
+		value: 'arcsec', $
+		comment: 'Units along axis 2'}
+	wcs.add, { $
+		name: 'PC1_1', $
+		value: pc[0, 0], $
+		comment: 'WCS coordinate transformation matrix'}
+	wcs.add, { $
+		name: 'PC1_2', $
+		value: pc[1, 0], $
+		comment: 'WCS coordinate transformation matrix'}
+	wcs.add, { $
+		name: 'PC2_1', $
+		value: pc[0, 1], $
+		comment: 'WCS coordinate transformation matrix'}
+	wcs.add, { $
+		name: 'PC2_2', $
+		value: pc[1, 1], $
+		comment: 'WCS coordinate transformation matrix'}
+	wcs.add, { $
+		name: 'CDELT1', $
+		value: cdelt[0], $
+		comment: '[arcsec] Pixel scale along axis 1'}
+	wcs.add, { $
+		name: 'CDELT2', $
+		value: cdelt[1], $
+		comment: '[arcsec] Pixel scale along axis 2'}
+	wcs.add, { $
+		name: 'CROTA', $
+		value: atan(pc[0, 1], pc[0, 0]) * !radeg, $
+		comment: '[deg] Rotation angle'}
+	wcs.add, { $
+		name: 'CRVAL1', $
+		value: crval[0], $
+		comment: '[arcsec] Value of reference pixel along axis 1'}
+	wcs.add, { $
+		name: 'CRVAL2', $
+		value: crval[1], $
+		comment: '[arcsec] Value of reference pixel along axis 2'}
+	wcs.add, { $
+		name: 'CRPIX1', $
+		value: crpix[0], $
+		comment: '[pixel] Reference pixel location along axis 1'}
+	wcs.add, { $
+		name: 'CRPIX2', $
+		value: crpix[1], $
+		comment: '[pixel] Reference pixel location along axis 2'}
+	wcs.add, { $
+		name: 'SUN_XCEN', $
+		value: sunpix[0], $
+		comment: '[pixel] Sun center location along axis 1'}
+	wcs.add, { $
+		name: 'SUN_YCEN', $
+		value: sunpix[1], $
+		comment: '[pixel] Sun center location along axis 2'}
+	wcs.add, { $
+		name: 'IO_XCEN', $
+		value: borpix[0], $
+		comment: '[pixel] Metis IO center location along axis 1'}
+	wcs.add, { $
+		name: 'IO_YCEN', $
+		value: borpix[1], $
+		comment: '[pixel] Metis IO center location along axis 2'}
+	wcs.add, { $
+		name: 'FS_XCEN', $
+		value: crpix[0], $
+		comment: '[pixel] Metis field-stop center location along axis 1'}
+	wcs.add, { $
+		name: 'FS_YCEN', $
+		value: crpix[1], $
+		comment: '[pixel] Metis field-stop center location along axis 2'}
+	wcs.add, { $
+		name: 'SC_XCEN', $
+		value: pntpix[0], $
+		comment: '[pixel] S/C pointing location along axis 1'}
+	wcs.add, { $
+		name: 'SC_YCEN', $
+		value: pntpix[1], $
+		comment: '[pixel] S/C pointing location along axis 2'}
+	wcs.add, { $
+		name: 'SC_YAW', $
+		value: pointing[0], $
+		comment: '[arcsec] S/C HPC yaw'}
+	wcs.add, { $
+		name: 'SC_PITCH', $
+		value: pointing[1], $
+		comment: '[arcsec] S/C HPC pitch'}
+	wcs.add, { $
+		name: 'SC_ROLL', $
+		value: pointing[2], $
+		comment: '[deg] S/C HPC roll angle'}
+	wcs.add, { $
+		name: 'INN_FOV', $
+		value: 1.6, $
+		comment: '[deg] Inner Metis FOV'}
+	wcs.add, { $
+		name: 'OUT_FOV', $
+		value: 3.4, $
+		comment: '[deg] Outer Metis FOV'}
+
+	return, wcs
+end

solo_get_carrot.pro

+function solo_get_carrot, utc
+
+	; initial estimate of the carrington rotation number
+
+	jul_day = date_conv(utc, 'julian')
+
+	carr = (jul_day - 2398167.d0) / 27.2753d0 + 1
+
+	; convert the requested date into ephemeris time
+
+	cspice_str2et, utc, et
+
+	; get the carrington longitude of earth
+
+	cspice_spkezr, 'EARTH', et, 'IAU_SUN', 'NONE', 'SUN', state, ltime
+
+	; convert rectangular coordinates into angular coordinates
+
+	cspice_reclat, state[0 : 2], sun_dist, he_lon, he_lat
+
+	; calculate the fractional part of the decimal rotation number
+
+	if he_lon lt 0. then he_lon = he_lon + 2. * !dpi
+	frac = 1.d0 - he_lon / (2. * !dpi)
+	n_carr = round(carr - frac)
+	carr = n_carr + frac
+
+	; correction for s/c position
+
+	cspice_spkezr, 'SOLO', et, 'SUN_EARTH_CEQU', 'NONE', 'SUN', state, ltime
+
+	; convert rectangular coordinates into angular coordinates
+
+	cspice_reclat, state[0 : 2], so_dist, so_lon, so_lat
+
+	carr = carr - so_lon / (2. * !dpi)
+
+	return, carr
+end

solo_get_coords.pro

+function solo_get_coords, utc, frame, obs, velocity = velocity, spherical = spherical, degrees = degrees
+
+	; convert the requested date into ephemeris time
+
+	cspice_str2et, utc, et
+
+	; switch to the proper solar orbiter frame
+
+	if ('carrington').contains(frame, /fold) then frame = 'IAU_SUN'
+
+	case frame.toupper() of
+		'RTN' : frame = 'SOLO_SUN_RTN'
+		'HEE' : frame = 'SUN_EARTH_ECL'
+		'HCI' : frame = 'SUN_INERTIAL'
+		'HAE' : frame = 'SUN_ARIES_ECL'
+		'HEQ' : frame = 'SUN_EARTH_CEQU'
+		'GSE' : frame = 'EARTH_SUN_ECL'
+		'GEI' : frame = 'J2000'
+		else : frame = frame
+	endcase
+
+	; coordinates of the s/c in the selected reference frame
+
+	cspice_spkezr, 'SOLO', et, frame, 'NONE', obs, state, ltime
+
+	coord = state[0 : 2]
+
+	; convert rectangular coordinates into angular coordinates
+
+	if keyword_set(spherical) then begin
+		cspice_reclat, coord, solo_dist, solo_lon, solo_lat
+		if frame eq 'IAU_SUN' then solo_lon = (solo_lon + 2. * !dpi) mod (2. * !dpi)
+		if keyword_set(degrees) then begin
+			solo_lat = solo_lat * 180.d0 / !dpi
+			solo_lon = solo_lon * 180.d0 / !dpi
+		endif
+		coord = [solo_dist, solo_lon, solo_lat]
+	endif
+
+	velocity = state[3 : 5]
+
+	return, coord
+end

solo_get_ephemerides.pro

+function solo_get_ephemerides, utc, constants = constants
+
+    if ~ keyword_set(constants) then begin
+        rsun = 695508000.0 ; metres
+        au = 149597870691.0 ; metres
+    endif else begin
+        rsun = constants.rsun.value
+        au = constants.au.value
+    endelse
+
+    ; compute and add the wcs keyword
+
+    ephemerides = list()
+
+    ephemerides.add, { $
+        name: 'LONPOLE', $
+        value: 180., $
+        comment: '[deg] Native longitude of the celestial pole'}
+
+    ; spherical coordinates of the S/C in the heeq (i.e., stonyhurst) frame
+
+    coord = solo_get_coords(utc, 'HEQ', 'SUN', /spherical, /degrees)
+
+    ; solar photospheric radius as seen from the S/C
+
+    solo_dist = coord[0] * 1000.
+    rsun_arc = atan(rsun, solo_dist) * !radeg * 3600.d0
+
+    ephemerides.add, { $
+        name: 'RSUN_ARC', $
+        value: rsun_arc, $
+        comment: '[arcsec] Apparent photospheric solar radius'}
+    ephemerides.add, { $
+        name: 'RSUN_REF', $
+        value: rsun, $
+        comment: '[m] Assumed physical solar radius'}
+
+    solar_angles = solo_get_solar_angles(utc)
+
+    ephemerides.add, { $
+        name: 'SOLAR_B0', $
+        value: solar_angles[0], $
+        comment: '[deg] S/C tilt of solar North pole'}
+    ephemerides.add, { $
+        name: 'SOLAR_P0 ', $
+        value: solar_angles[1], $
+        comment: '[deg] S/C celestial North to solar North angle'}
+    ephemerides.add, { $
+        name: 'SOLAR_EP', $
+        value: solar_angles[2], $
+        comment: '[deg] S/C ecliptic North to solar North angle'}
+
+    carrot = solo_get_carrot(utc)
+
+    ephemerides.add, { $
+        name: 'CAR_ROT', $
+        value: carrot, $
+        comment: 'Carrington rotation number'}
+
+    ephemerides.add, { $
+        name: 'HGLT_OBS', $
+        value: coord[2], $
+        comment: '[deg] S/C heliographic latitude (B0 angle)'}
+    ephemerides.add, { $
+        name: 'HGLN_OBS', $
+        value: coord[1], $
+        comment: '[deg] S/C heliographic longitude'}
+
+    ; coordinates of the S/C in the carrington frame and carrington rotation number
+
+    coord = solo_get_coords(utc, 'CARRINGTON', 'SUN', /spherical, /degrees)
+
+    ephemerides.add, { $
+        name: 'CRLT_OBS', $
+        value: coord[2], $
+        comment: '[deg] S/C Carrington latitude (B0 angle)'}
+    ephemerides.add, { $
+        name: 'CRLN_OBS', $
+        value: coord[1], $
+        comment: '[deg] S/C Carrington longitude (L0 angle)'}
+
+    ephemerides.add, { $
+        name: 'DSUN_OBS', $
+        value: solo_dist, $
+        comment: '[m] S/C distance from Sun'}
+    ephemerides.add, { $
+        name: 'DSUN_AU', $
+        value: solo_dist/au, $
+        comment: '[AU] S/C distance from Sun'}
+    ephemerides.add, { $
+        name: 'AU_REF', $
+        value: au, $
+        comment: '[m] Assumed physical Astronomical Unit'}
+
+    ; coordinates of the S/C in the hee frame
+
+    coord = solo_get_coords(utc, 'HEE', 'SUN') * 1000.
+
+    ephemerides.add, { $
+        name: 'HEEX_OBS', $
+        value: coord[0], $
+        comment: '[m] S/C Heliocentric Earth Ecliptic X'}
+    ephemerides.add, { $
+        name: 'HEEY_OBS', $
+        value: coord[1], $
+        comment: '[m] S/C Heliocentric Earth Ecliptic Y'}
+    ephemerides.add, { $
+        name: 'HEEZ_OBS', $
+        value: coord[2], $
+        comment: '[m] S/C Heliocentric Earth Ecliptic Z'}
+
+    ; coordinates of the S/C in the hci frame
+
+    coord = solo_get_coords(utc, 'HCI', 'SUN', velocity = veloc)
+    coord = coord * 1000.
+    veloc = veloc * 1000.
+
+    ephemerides.add, { $
+        name: 'HCIX_OBS', $
+        value: coord[0], $
+        comment: '[m] S/C Heliocentric Inertial X'}
+    ephemerides.add, { $
+        name: 'HCIY_OBS', $
+        value: coord[1], $
+        comment: '[m] S/C Heliocentric Inertial Y'}
+    ephemerides.add, { $
+        name: 'HCIZ_OBS', $
+        value: coord[2], $
+        comment: '[m] S/C Heliocentric Inertial Z'}
+
+    ephemerides.add, { $
+        name: 'HCIX_VOB', $
+        value: veloc[0], $
+        comment: '[m/s] S/C Heliocentric Inertial X velocity'}
+    ephemerides.add, { $
+        name: 'HCIY_VOB', $
+        value: veloc[1], $
+        comment: '[m/s] S/C Heliocentric Inertial Y velocity'}
+    ephemerides.add, { $
+        name: 'HCIZ_VOB', $
+        value: veloc[2], $
+        comment: '[m/s] S/C Heliocentric Inertial Z velocity'}
+
+    ; coordinates of the S/C in the hae frame
+
+    coord = solo_get_coords(utc, 'HAE', 'SUN')
+    coord = coord * 1000.
+
+    ephemerides.add, { $
+        name: 'HAEX_OBS', $
+        value: coord[0], $
+        comment: '[m] S/C Heliocentric Aries Ecliptic X'}
+    ephemerides.add, { $
+        name: 'HAEY_OBS', $
+        value: coord[1], $
+        comment: '[m] S/C Heliocentric Aries Ecliptic Y'}
+    ephemerides.add, { $
+        name: 'HAEZ_OBS', $
+        value: coord[2], $
+        comment: '[m] S/C Heliocentric Aries Ecliptic Z'}
+
+    ; coordinates of the S/C in the heeq frame
+
+    coord = solo_get_coords(utc, 'HEQ', 'SUN')
+    coord = coord * 1000.
+
+    ephemerides.add, { $
+        name: 'HEQX_OBS', $
+        value: coord[0], $
+        comment: '[m] S/C Heliocentric Earth Equatorial X'}
+    ephemerides.add, { $
+        name: 'HEQY_OBS', $
+        value: coord[1], $
+        comment: '[m] S/C Heliocentric Earth Equatorial Y'}
+    ephemerides.add, { $
+        name: 'HEQZ_OBS', $
+        value: coord[2], $
+        comment: '[m] S/C Heliocentric Earth Equatorial Z'}
+
+    ; coordinates of the S/C in the gse frame
+
+    coord = solo_get_coords(utc, 'GSE', 'EARTH')
+    coord = coord * 1000.
+
+    ephemerides.add, { $
+        name: 'GSEX_OBS', $
+        value: coord[0], $
+        comment: '[m] S/C Geocentric Solar Ecliptic X'}
+    ephemerides.add, { $
+        name: 'GSEY_OBS', $
+        value: coord[1], $
+        comment: '[m] S/C Geocentric Solar Ecliptic Y'}
+    ephemerides.add, { $
+        name: 'GSEZ_OBS', $
+        value: coord[2], $
+        comment: '[m] S/C Geocentric Solar Ecliptic Z'}
+
+    ; light travel times and radial velocity of the S/C
+
+    earth_time = get_light_time(utc, 'EARTH', 'SUN')
+    sun_time = get_light_time(utc, 'SOLO', 'SUN', radvel = radvel)
+    tdel = earth_time - sun_time
+
+    ephemerides.add, { $
+        name: 'OBS_VR', $
+        value: radvel, $
+        comment: '[m/s] Radial velocity of S/C relative to Sun'}
+    ephemerides.add, { $
+        name: 'EAR_TDEL', $
+        value: tdel, $
+        comment: '[s] Time(Sun to Earth) - Time(Sun to S/C)'}
+    ephemerides.add, { $
+        name: 'SUN_TIME', $
+        value: sun_time, $
+        comment: '[s] Time(Sun to S/C)'}
+
+    ; corrections of the acquisition date
+
+    jul_utc = date_conv(utc, 'julian')
+    date_earth = date_conv(jul_utc + tdel / 86400.d0, 'fits')
+    date_sun = date_conv(jul_utc - sun_time / 86400.d0, 'fits')
+
+    ephemerides.add, { $
+        name: 'DATE_EAR', $
+        value: date_earth, $
+        comment: '[UTC] Start time of observation, corrected to Earth'}
+    ephemerides.add, { $
+        name: 'DATE_SUN', $
+        value: date_sun, $
+        comment: '[UTC] Start time of observation, corrected to Sun'}
+
+    return, ephemerides
+end

solo_get_pointing.pro

+function solo_get_pointing, utc, degrees = degrees, radians = radians, celestial = celestial
+
+	; convert the requested date into ephemeris time
+
+	cspice_str2et, utc, et
+
+	; convert the ephemeris time into s/c internal time
+
+	cspice_sce2c, -144L, et, clk_in
+
+	; set the proper solar orbiter frame
+
+	frame = 'SOLO_SUN_RTN'
+
+	if keyword_set(celestial) then frame = 'J2000'
+
+	cspice_ckgp, -144000L, clk_in, 1, frame, cmat, clk_out, found
+
+	if found then cspice_m2eul, cmat, 1, 2, 3, roll, pitch, yaw else return, replicate(0., 3)
+
+	; make the conversion form rtn to hpc, if necessary
+
+	if frame eq 'SOLO_SUN_RTN' then begin
+		yaw = !dpi * signum(yaw) - yaw
+		pitch = -pitch
+		roll = -roll
+	endif else begin
+		pitch = -pitch
+	endelse
+
+	; correct any cases where the roll is greater than +/- 180 degrees
+
+	if abs(roll) gt !dpi then roll = roll - 2.d0 * !dpi * signum(roll)
+
+	; correct any cases where the pitch is greater than +/- 90 degrees
+
+	if abs(pitch) gt !dpi / 2.d0 then begin
+		yaw = yaw - !dpi * signum(yaw)
+		pitch = !dpi * signum(pitch) - pitch
+		roll = roll - !dpi * signum(roll)
+	endif
+
+	; apply correct units to the pointing vector
+
+	if keyword_set(radians) then rad_factor = 1. else rad_factor = 180.d0 / !dpi
+	if keyword_set(radians) or keyword_set(degrees) then arc_factor = 1. else arc_factor = 3600.
+
+	yaw = yaw * rad_factor * arc_factor
+	pitch = pitch * rad_factor * arc_factor
+	roll = roll * rad_factor
+
+	return, [yaw, pitch, roll]
+end

solo_get_solar_angles.pro

+function solo_get_solar_angles, utc
+
+    ; convert the requested date into ephemeris time
+
+    cspice_utc2et, utc, et
+
+    cspice_pxform, 'IAU_SUN', 'J2000', et, rot
+    sun_north = rot[2, *]
+
+    cspice_spkezr, 'SUN', et, 'J2000', 'NONE', 'SOLO', state, ltime
+    coord = state[0 : 2]
+    rad = coord / sqrt(total(coord^2, 1))
+
+    cel_north = [0., 0., 1.d0]
+
+    sun_proj = sun_north - total(rad * sun_north) * rad
+    sun_proj = sun_proj / sqrt(total(sun_proj^2))
+
+    cel_proj = cel_north - total(rad * cel_north) * rad
+    cel_proj = cel_proj / sqrt(total(cel_proj^2))
+
+    vec_proj = crossp(cel_north, rad)
+    vec_proj = vec_proj / sqrt(total(vec_proj^2))
+
+    x_proj = total(sun_proj * cel_proj)
+    y_proj = total(sun_proj * vec_proj)
+
+    p0 = atan(y_proj, x_proj) * 180.d0 / !dpi
+
+    cspice_pxform, 'IAU_SUN', 'ECLIPJ2000', et, rot
+    sun_north = rot[2, *]
+
+    cspice_spkezr, 'SUN', et, 'ECLIPJ2000', 'NONE', 'SOLO', state, ltime
+    coord = state[0 : 2]
+    rad = coord / sqrt(total(coord^2, 1))
+
+    eclip_north = cel_north
+
+    sun_proj = sun_north - total(rad * sun_north) * rad
+    sun_proj = sun_proj / sqrt(total(sun_proj^2))
+	
+    ecl_proj = eclip_north - total(rad * eclip_north) * rad
+    ecl_proj = ecl_proj / sqrt(total(ecl_proj^2))
+
+    vec_proj = crossp(eclip_north, rad)
+    vec_proj = vec_proj / sqrt(total(vec_proj^2))
+
+    x_proj = total(sun_proj * ecl_proj)
+    y_proj = total(sun_proj * vec_proj)
+
+    ep = atan(y_proj, x_proj) * 180.d0 / !dpi
+
+    cspice_spkezr, 'SUN', et, 'IAU_SUN', 'NONE', 'SOLO', state, ltime
+    coord = state[0 : 2]
+
+    cspice_reclat, coord, dist, lon, lat
+
+    b0 = -lat * 180.d0 / !dpi
+
+    return, [b0, p0, ep]
+end