Version 3.0

CHANGELOG.md

 # CHANGELOG
 All changes to the Metis L1 pipeline are documented in this file.
 
+## 3.0 - 2021-07-19
+- Add the WCS (World Coordinate System) keywords to the FITS header and a set of routines to perform calculation of the instrument boresight parameters and S/C coordinates.
+- Modify the pipeline to handle the new input I/F which includes the full set of Solar Orbiter SPICE kernels and the Metis Calibration Package.
+- Add the derivation of the quality matrix which is now added in an extension to the L1 FITS file.
+- Add several code optimisations.
+
 ## 2.3.1 - 2021-04-19
-- Fixed a bug that caused an incorrect FITS file format: images with unsigned or long integer data type were saved as signed integers with the wrong BSCALE and BZERO scale keywords.
+- Fix a bug that caused an incorrect FITS file format: images with unsigned or long integer data type were saved as signed integers with the wrong BSCALE and BZERO scale keywords.
 
 ## 2.3 - 2021-03-10
 - Fix a bug that caused incorrect pixel counts in re-binned images. Pixel counts are now correctly multiplied by the square of the binning factor after re-binning.
@@ -41,7 +47,7 @@ All changes to the Metis L1 pipeline are documented in this file.
 - Revise the structure of the binary table containing the house-keeping parameters to make it more user-friendly.
 
 ## 1.0 – 2020-03-06
-- Add several optimisations to the code.
+- Add several code optimisations.
 - Add the decode_obt.pro routine.
 
 ## 0.0 – 2020-01-16

check_quality.pro

+function check_quality, data, cal_pack, filter = filter
+
+	if filter.contains('VL', /fold) then $
+		channel = cal_pack.vl_channel else $
+		channel = cal_pack.uv_channel
+
+	levels = channel.sat_level.value
+
+	quality_matrix = float(data) * 0.
+
+	s1 = where(data le levels[0], count)
+	if count gt 0 then quality_matrix[s1] = !values.f_nan
+
+	s2 = where(data gt levels[0] and data le levels[1], count)
+	if count gt 0 then quality_matrix[s2] = 0
+
+	s3 = where(data gt levels[1] and data lt levels[2], count)
+	if count gt 0 then quality_matrix[s3] = 1
+
+	s4 = where(data ge levels[2] and data lt levels[3], count)
+	if count gt 0 then quality_matrix[s4] = 0
+
+	s5 = where(data ge levels[3], count)
+	if count gt 0 then quality_matrix[s5] = !values.f_nan
+
+	return, quality_matrix
+end

decode_obt.pro

-function decode_obt, t_coarse, t_fine, to_double = to_double, from_double = from_double
-	if keyword_set(to_double) then return, double(t_coarse) + double(t_fine)/65536.0D
-	if keyword_set(from_double) then obt = [ulong(t_coarse), uint((t_coarse - ulong(t_coarse)) * 65536)] else obt = [t_coarse, t_fine]
+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

interpol_param.pro

@@ -10,6 +10,7 @@ function interpol_param, table, par_name, date, empty_params = empty_params
 	s = sort(par_date)
 	par_val = par_val[s]
 	par_date = par_date[s]
-	value = interpol(par_val, par_date, date_conv(date, 'JULIAN'))
+	jul_date = date_conv(date, 'JULIAN')
+	if jul_date ge min(par_date) and jul_date le max(par_date) then value = interpol(par_val, par_date, jul_date) else value = 0.0
 	if finite(value) then return, value else return, 0.0
 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_l1_prep.pro

 pro metis_l1_prep
 
+	; keyword defining if the detector reference frame must be used for the output
+
+	ref_detector = 1
+
 	; start the log
 
 	journal, 'output/metis_l1_prep_log.txt'
@@ -19,16 +23,6 @@ pro metis_l1_prep
 		'VL light-curve' $					; 9
 	)
 
-	metis_obj_size = list( $
-		2048L, $					; 0
-		1024L, $					; 1
-		1024L, $					; 2
-		2048L, $					; 3
-		1024L, $					; 4
-		2048L, $					; 5
-		1024L $						; 6
-	)
-
 	; read the auxiliary input file
 
 	input = json_parse('input/contents.json', /toarray, /tostruct)
@@ -37,64 +31,47 @@ pro metis_l1_prep
 
 	; load the spice kernels
 
-	kernels = [ $
-		input.tls_file_name, $
-		input.tsc_file_name $
-	]
+	load_spice_kernels, input.spice_kernels, kernel_list = kernel_list, kernel_version = kernel_version
 
-	; load the spice kernels
+	journal, 'SPICE kernel files correctly loaded:'
+	journal, '  SDK version= ' + kernel_version
 
-	cspice_furnsh, kernels
+	; import the calibration package index file
 
-	journal, 'SPICE kernel files correctly loaded:'
-	journal, '  sclk kernel = ' + input.tsc_file_name
-	journal, '  leap seconds kernel = ' + input.tls_file_name
+	cal_pack = json_parse(input.cal_pack_path + '/index.json', /toarray, /tostruct)
 
-	; read l0 fits structure
+	; include the calibration package path into the cal_pack structure
 
-	fits_info, input.file_name, n_ext = n_ext, extname = extname, /silent
+	cal_pack = create_struct('path', input.cal_pack_path + path_sep(), cal_pack)
 
 	; read the primary hdu
 
-	ext_no = 0
-	data = mrdfits(input.file_name, ext_no, primary_header, /silent)
-
-	; if the data is not an image, read the data binary-table extension
+	data = mrdfits(input.file_name, 0, primary_header, /silent)
 
-	if fxpar(primary_header, 'NAXIS') eq 0 then begin
-		ext_no += 1
-		data_bin_table = mrdfits(input.file_name, ext_no, data_extension_header, /silent)
-	endif else data_bin_table = !null
+	naxis = fxpar(primary_header, 'NAXIS')
+	naxis1 = fxpar(primary_header, 'NAXIS1', missing = 0)
+	naxis2 = fxpar(primary_header, 'NAXIS2', missing = 0)
 
 	; read the metadata extension
 
-	ext_no += 1
-	metadata_bin_table = mrdfits(input.file_name, ext_no, metadata_extension_header, /silent)
+	metadata_bin_table = mrdfits(input.file_name, 'metis metadata', metadata_extension_header, /silent)
 
 	; identify the data product and its nominal size
 
 	datatype = fxpar(metadata_extension_header, 'DATATYPE')
 
-	if datatype le 6 then begin
-		naxis = 2
-		naxis1 = metis_obj_size[datatype]
-		naxis2 = naxis1
-	endif else begin
-		naxis = 0
-	endelse
+	journal, 'L0 FITS file correctly read:'
+	journal, '  datatype = ' + string(datatype, format = '(I0)')
 
-	; if the data product is an accumulation matrix, look for the accumulation vector extension and read it if it exists
+	; if the data product is a light curve or a pcu-event list, read the data binary hk_table
 
-	if datatype eq 2 then begin
-		if fxpar(metadata_extension_header, 'M') eq 256 then begin
-			ext_no += 1
-			accumul_vector = mrdfits(input.file_name, ext_no, vector_extension_header, /silent)
-		endif else accumul_vector = !null
-	endif
+	if datatype gt 6 then data_bin_table = mrdfits(input.file_name, 1, data_extension_header, /silent) else data_bin_table = !null
 
-	journal, 'L0 FITS file correctly read:'
-	journal, '  datatype = ' + string(datatype, format = '(I0)')
-	if datatype le 6 then journal, '  size = ' + string(naxis1, format = '(I0)') + '×' + string(naxis2, format = '(I0)')
+	; if the data product is an accumulation matrix, look for the accumulation vector extension and read it if it exists
+
+	if datatype eq 2 then $
+		if fxpar(metadata_extension_header, 'M') eq 256 then $
+			accumul_vector = mrdfits(input.file_name, 'accumulation vector', vector_extension_header, /silent) else accumul_vector = !null
 
 	; NOTE - the radialization extension is ignored
 
@@ -115,9 +92,9 @@ pro metis_l1_prep
 	obt_end = fxpar(primary_header, 'OBT_END')
 	obt_avg = (obt_beg + obt_end)/2.0D
 
-	date_beg = solo_obt2utc(decode_obt(obt_beg, /from_double))
-	date_end = solo_obt2utc(decode_obt(obt_end, /from_double))
-	date_avg = solo_obt2utc(decode_obt(obt_avg, /from_double))
+	date_beg = solo_obt2utc(obt_beg)
+	date_end = solo_obt2utc(obt_end)
+	date_avg = solo_obt2utc(obt_avg)
 
 	date_beg_string = strmid(date_beg, 0, 19)
 	date_beg_string = date_beg_string.replace('-', '')
@@ -128,27 +105,36 @@ pro metis_l1_prep
 	; name of the fits file
 
 	file_name_fields = strsplit(fxpar(primary_header, 'FILENAME'), '_', /extract)
-
 	file_name = 'solo_L1_' + file_name_fields[2] + '_' + date_beg_string + '_V' + version + '.fits'
 	out_file_name = 'output/' + file_name
 
-	; instrument keywords
+	; filter keyword
+
+	case datatype of
+		0: filter = 'VL'
+		1: filter = 'UV'
+		2: filter = 'UV'
+		3: filter = 'VL'
+		4: filter = 'UV'
+		5: filter = 'VL'
+		6: filter = 'UV'
+		7: filter = 'UV'
+		8: filter = 'UV'
+		9: filter = 'VL'
+	endcase
+
+	; choose the correct calibration-package structure based on the image filter
+
+	if filter eq 'VL' then channel = cal_pack.vl_channel else channel = cal_pack.uv_channel
+
+	; instrument/telescope/detector keywords
 
-	if datatype eq 0 or datatype eq 3 or datatype eq 5 or datatype eq 9 then begin
-		filter = 'VL'
-		wavelnth = 610.0
-		wavemin = 580.0
-		wavemax = 640.0
-		waveband = 'Visible light 580-640 nm'
-	endif else begin
-		filter = 'UV'
-		wavelnth = 121.6
-		wavemin = 111.6
-		wavemax = 131.6
-		waveband = 'H I Lyman-alpha 121.6 nm'
-	endelse
 	detector = filter + 'D'
+	waveband = cal_pack.vl_channel.name
 	telescope = 'SOLO/Metis/' + detector
+	wavelnth = channel.bandpass.value[1]
+	wavemin = channel.bandpass.value[0]
+	wavemax = channel.bandpass.value[2]
 
 	; campaign keywords
 
@@ -158,7 +144,7 @@ pro metis_l1_prep
 
 	; build the fits file extensions
 
-	; join the metadata extension header to the primary header removing unwanted keywords
+	; join the metadata extension header to the primary header removing useless keywords
 
 	i = where(strmid(primary_header, 0, 8) eq 'COMP_RAT')
 	j = where(strmid(metadata_extension_header, 0, 8) eq 'EXTNAME ')
@@ -170,69 +156,88 @@ pro metis_l1_prep
 		primary_header[i + 1: *] $
 	]
 
-	; rebin the image if binning was applied during the acquisition
+	; rebin the image if binning was applied during the acquisition and check for the data quality
+	; NOTE - this is done only for image data products
 
 	if datatype le 6 then begin
-		if fxpar(primary_header, 'COMPR') then $
-			bin_type = fxpar(primary_header, 'B0_BIN') else bin_type = 0
+		bin_type = fxpar(primary_header, 'B0_BIN', missing = 0)
 
-		pol_id = fxpar(primary_header, 'POL_ID')
-		fpfilt = fxpar(primary_header, 'VLFPFILT', count = count)
-		if count gt 0 and pol_id eq 0 and fpfilt eq 0 and bin_type eq 0 then bin_type = 1
+		; check for wrong binning type in fixed-polarization filtered images
+
+		pol_id = fxpar(primary_header, 'POL_ID', missing = -1)
+		vlfpfilt = fxpar(primary_header, 'VLFPFILT', missing = -1)
+		if vlfpfilt eq 0 and pol_id eq 0 and bin_type eq 0 then bin_type = 1
+
+		bin_fact = 2^bin_type ; binning factor = 1, 2, or 4
+
+		if bin_fact gt 1 then begin
+			naxis1 = naxis1/bin_fact
+			naxis2 = naxis2/bin_fact
+			data = rebin(data, naxis1, naxis2)
+
+			; evaluate the quality matrix
+
+			quality_matrix = check_quality(data, cal_pack, filter = filter)
 
-		bin_fact = 2. ^ bin_type
-		if bin_fact gt 1. then begin
-			data = rebin(data, naxis1 / bin_fact, naxis2 / bin_fact)
-			data = data * bin_fact^2
+			; correct values including binning factor only for vl and uv images and accumulation matrices
+
+			if datatype le 2 then data = long(data) * bin_fact^2
 
 			if ~ isa(comment) then comment = !null
 			comment = [comment, 'Image was rebinned according to the commanded binning factor.']
 
 			journal, 'Data was binned during acquisition:'
 			journal, '  binning = ' + string(bin_fact, format = '(I1)') + '×' + string(bin_fact, format = '(I1)')
-			journal, '  data was correctly rebinned:'
-			journal, '    new size = ' + string(naxis1 / bin_fact, format = '(I0)') + '×' + string(naxis1 / bin_fact, format = '(I0)')
+			journal, '  data was rebinned:'
+			journal, '    new size = ' + string(naxis1, format = '(I0)') + '×' + string(naxis2, format = '(I0)')
 		endif else begin
+			quality_matrix = check_quality(data, cal_pack, filter = filter)
+
 			journal, 'Data was not binned during acquistion.'
 		endelse
-	endif
+
+		; read the bad-pixel map
+
+		bad_pixel_map = readfits(cal_pack.path + channel.bad_pixel_map[0].file_name, /silent)
+
+		; rebin the bad-pixel map to the size of the image to be corrected
+
+		bad_pixel_map = rebin(bad_pixel_map, naxis1, naxis2) * bin_fact^2
+
+		; correct the quality matrix according to the bad-pixel map
+
+		bad_pixels = where(bad_pixel_map gt 0, count)
+		if count gt 0 then quality_matrix[bad_pixels] = 0
+	endif else quality_matrix = !null
 
 	; correct for the effective exposure time
+	; NOTE - this is done only for vl and uv images and accumulation matrices
 
 	telapse = float(obt_end - obt_beg)
 
 	dit = fxpar(metadata_extension_header, 'DIT')
-	ndit = fxpar(metadata_extension_header, 'NDIT', missing = 1L)
-	ndit1 = fxpar(metadata_extension_header, 'NDIT1', missing = 1L)
-	ndit2 = fxpar(metadata_extension_header, 'NDIT2', missing = 1L)
+	ndit = fxpar(metadata_extension_header, 'NDIT', missing = 1)
+	ndit1 = fxpar(metadata_extension_header, 'NDIT1', missing = 1)
+	ndit2 = fxpar(metadata_extension_header, 'NDIT2', missing = 1)
 
-	if datatype le 6 then begin
+	if datatype le 2 then begin
 		nsumexp = ndit * ndit1 * ndit2
 		xposure = dit/1000. * nsumexp
-		data = data * nsumexp
-
-		if max(data, /nan) lt 32768 then data = fix(data)
-
-		datamin = min(data, /nan)
-		fxaddpar, primary_header, 'DATAMIN', datamin
-
-		datamax = max(data, /nan)
-		fxaddpar, primary_header, 'DATAMAX', datamax
+		data = long(data) * nsumexp
 
 		if ~ isa(comment) then comment = !null
 		comment = [comment, 'Image values were corrected for the total exposure time.']
 
-		if nsumexp gt 1 then journal, 'Image values were corrected for the total exposure time.'
+		journal, 'Image values were corrected for the total exposure time.'
 	endif else begin
-		xposure = dit/1000.
 		nsumexp = 1
-		data = !null
+		xposure = dit/1000.
 	endelse
 
 	; adjust the primary header (it is almost the same for all data product types)
 
 	fxaddpar, primary_header, 'FILENAME', file_name
-	fxaddpar, primary_header, 'PARENT', file_basename(input.file_name), 'Name of the parent file that got processed to the current one', before = 'APID'
+	fxaddpar, primary_header, 'PARENT', file_basename(input.file_name), 'Name of the parent file', before = 'APID'
 	fxaddpar, primary_header, 'DATE', date, 'Date and time of FITS file creation', before = 'OBT_BEG'
 	fxaddpar, primary_header, 'DATE-OBS', date_beg, 'Same as DATE-BEG', before = 'OBT_BEG'
 	fxaddpar, primary_header, 'DATE-BEG', date_beg, 'Start time of observation', before = 'OBT_BEG'
@@ -244,6 +249,7 @@ pro metis_l1_prep
 	fxaddpar, primary_header, 'LEVEL', 'L1'
 	fxaddpar, primary_header, 'CREATOR', 'metis_l1_prep.pro'
 	fxaddpar, primary_header, 'VERS_SW', input.sw_version
+	fxaddpar, primary_header, 'VERS_CAL', cal_pack.version, 'Version of the calibration package', after = 'VERS_SW'
 	fxaddpar, primary_header, 'OBSRVTRY', 'Solar Orbiter', 'Satellite name', before = 'INSTRUME'
 	fxaddpar, primary_header, 'TELESCOP', telescope, 'Telescope that took the measurement', before = 'INSTRUME'
 	fxaddpar, primary_header, 'DETECTOR', detector, 'Subunit/sensor', before = 'DATAMIN'
@@ -266,6 +272,9 @@ pro metis_l1_prep
 	fxaddpar, primary_header, 'BZERO', 0, 'Physical value for the array value 0', before = 'DATAMIN'
 	fxaddpar, primary_header, 'BTYPE', metis_datatype[datatype], 'Science data object type', before = 'DATAMIN'
 	fxaddpar, primary_header, 'BUNIT', 'DN', 'Units of physical value, after application of BSCALE and BZERO', before = 'DATAMIN'
+	fxaddpar, primary_header, 'BLANK', 0
+	fxaddpar, primary_header, 'DATAMIN', min(data, /nan)
+	fxaddpar, primary_header, 'DATAMAX', max(data, /nan)
 	fxaddpar, primary_header, 'IDB_VERS', input.idb_version, '', before = 'HDR_VERS'
 	fxaddpar, primary_header, 'INFO_URL', 'http://metis.oato.inaf.it', 'Link to more information on the instrument data', before = 'HISTORY'
 
@@ -310,7 +319,6 @@ pro metis_l1_prep
 
 	if datatype eq 1 or datatype eq 4 or datatype eq 6 then begin
 		fxaddpar, primary_header, 'HVU_SCR', interpol_param(hk_table, 'NIT0E070', date_avg, empty_params = empty_params), '[raw] HVU Screen commanded voltage'
-		; CMDHVSCR e CMDHVMCP
 		fxaddpar, primary_header, 'HVU_MCP', interpol_param(hk_table, 'NIT0E071', date_avg, empty_params = empty_params), '[raw] HVU MCP commanded voltage'
 		fxaddpar, primary_header, 'HV_SCR_V', interpol_param(hk_table, 'NIT0E0B7', date_avg, empty_params = empty_params), '[V] HVU Screen + MCP read voltage', after = 'HVU_MCP'
 		fxaddpar, primary_header, 'HV_MCP_V', interpol_param(hk_table, 'NIT0E0B6', date_avg, empty_params = empty_params), '[V] HVU MCP read voltage', after = 'HV_SCR_V'
@@ -330,7 +338,7 @@ pro metis_l1_prep
 
 	if empty_params eq !null then empty_params = ''
 
-	; replace verbose keyword values
+	; replace with text keyword values
 
 	key = fxpar(primary_header, 'MEASKIND', count = count)
 	if count gt 0 then begin
@@ -377,7 +385,26 @@ pro metis_l1_prep
 		endif
 	endforeach
 
-	; remove unused keywords
+	; convert the fits header into an idl structure
+
+	header = fits_hdr2struct(primary_header)
+
+	; append wcs keywords
+
+	if datatype le 6 then begin
+		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'
+	endif
+
+	; append solar keywords
+
+	ephemerides = solo_get_ephemerides(header, constants = cal_pack.constants)
+	foreach element, ephemerides do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
+
+	fxaddpar, primary_header, 'HISTORY', 'WCS and solar ephemerides:'
+	fxaddpar, primary_header, 'HISTORY', '  SKD version = ' + kernel_version
+
+	; remove redundant and misleading keywords
 
 	sxdelpar, primary_header, 'WIDTH'
 	sxdelpar, primary_header, 'HEIGHT'
@@ -411,10 +438,16 @@ pro metis_l1_prep
 			fxaddpar, primary_header, 'COMMENT', comment[k]
 	endif
 
-	; add checksum and datasum to the fits header
+	; rectify the image if requested
 
-	fits_add_checksum, primary_header, data
+	if not ref_detector then begin
+		data = metis_rectify(data, filter)
+		quality_matrix = metis_rectify(quality_matrix, filter)
+	endif
+
+	; add checksum and datasum to the fits header and save
 
+	fits_add_checksum, primary_header, data
 	mwrfits, data, out_file_name, primary_header, /no_comment, /create, /silent
 
 	journal, 'Fits file created:'
@@ -422,29 +455,37 @@ pro metis_l1_prep
 
 	; if applicable, save the data binary-table extension as it is
 
-	if isa(data_bin_table) then mwrfits, data_bin_table, 'output/' + file_name, data_extension_header, /no_comment, /silent
+	if isa(data_bin_table) then begin
+		fits_add_checksum, data_extension_header, data_bin_table
+		mwrfits, data_bin_table, out_file_name, data_extension_header, /no_comment, /silent
+	endif
+
+	; save the quality matrix
+
+	if datatype le 6 and isa(quality_matrix) then begin
+		quality_matrix_header = !null
+		fxaddpar, quality_matrix_header, 'PCOUNT', 0, 'Parameter count'
+		fxaddpar, quality_matrix_header, 'GCOUNT', 1, 'Group count'
+		fxaddpar, quality_matrix_header, 'EXTNAME', 'Quality matrix', 'Extension name'
+		fits_add_checksum, quality_matrix_header, quality_matrix
+		mwrfits, quality_matrix, out_file_name, quality_matrix_header, /no_comment, /silent
+	endif
 
 	; build the telemetry extension
 
 	hk_extension_header = !null
-	fxaddpar, hk_extension_header, 'XTENSION', 'BINTABLE', 'Extension type'
-	fxaddpar, hk_extension_header, 'BITPIX', 8, 'Number of bits per data pixel'
-	fxaddpar, hk_extension_header, 'NAXIS', 2, 'Number of data axes'
-	fxaddpar, hk_extension_header, 'NAXIS1', 0, 'Number of pixels along axis 1'
-	fxaddpar, hk_extension_header, 'NAXIS2', 0, 'Number of pixels along axis 2'
 	fxaddpar, hk_extension_header, 'PCOUNT', 0, 'Parameter count'
 	fxaddpar, hk_extension_header, 'GCOUNT', 1, 'Group count'
 	fxaddpar, hk_extension_header, 'EXTNAME', 'House-keeping', 'Extension name'
 
 	hk_bin_table = make_bin_table(hk_table)
-
 	mwrfits, hk_bin_table, out_file_name, hk_extension_header, /no_comment, /silent
 
 	journal, '  HK binary-table extension correctly added'
 
 	; unload the spice kernels
 
-	cspice_unload, kernels
+	load_spice_kernels, kernel_list = kernel_list, /unload
 
 	journal, 'SPICE kernel files unloaded.'
 

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_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, header, constants = constants
+
+	utc = header.date_avg
+
+	rsun = constants.rsun.value
+	au = constants.au.value
+
+	; 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
+
+	utc = header.date_beg
+
+	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

solo_obt2utc.pro

 function solo_obt2utc, obt
-	solo_id = -144
-	cspice_scs2e, solo_id, obt, et
+
+	; check if obt is in the correct form (string, with coarse and fine times)
+
+	if isa(obt, /number) then obt_decoded = decode_obt(obt, /from_decimal) else obt_decoded = obt
+
+	; convert the requested obt into ephemeris time
+
+	cspice_scs2e, -144L, obt_decoded, et
+
+	; convert the ephemeris time into utc
+
 	cspice_et2utc, et, 'ISOC', 3, utc
+
 	return, utc
 end

solo_utc2obt.pro

+function solo_utc2obt, utc
+
+	; convert the requested utc into ephemeris time
+
+	cspice_str2et, utc, et
+
+	; convert the ephemeris time into spacecraft obt
+
+	cspice_sce2s, -144L, et, clk_str
+	obt = clk_str.substring(2)
+	return, obt
+end