metis_l1_prep.pro

pro metis_l1_prep

	; start the log 

	journal, 'output/metis_l1_prep_log.txt'

	; some definitions

	metis_datatype = list( $
		'vl-image', $				; 0
		'uv-image', $				; 1
		'pcu-accumul', $			; 2
		'vl-temp-matrix', $			; 3
		'uv-temp-matrix', $			; 4
		'vl-c-ray-mat', $			; 5
		'uv-c-ray-mat', $			; 6
		'pcu-event-list', $			; 7
		'pcu-test-event-list', $	; 8
		'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)

	journal, 'File ' + input.file_name

	; prepare the spice kernels

	kernels = [ $
		input.tls_file_name, $
		input.tsc_file_name $
	]

	; load the spice kernels

	cspice_furnsh, kernels

	journal, 'SPICE kernel files correctly loaded:'
	journal, '  sclk kernel = ' + input.tsc_file_name
	journal, '  leap seconds kernel = ' + input.tls_file_name

	; read l0 fits structure

	fits_info, input.file_name, n_ext = n_ext, extname = extname, /silent

	; 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

	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

	; read the metadata extension

	ext_no += 1
	metadata_bin_table = mrdfits(input.file_name, ext_no, 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
		bitpix = 8
	endelse

	; 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 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

	journal, 'L0 FITS file correctly read:'
	journal, '  datatype = ' + string(datatype, format = '(I0)')
	journal, '  size = ' + string(naxis1, format = '(I0)') + '×' + string(naxis2, format = '(I0)')

	; NOTE - the radialization extension is ignored

	; read the planning data

	planning_data = json_parse(input.planning_file_name, /toarray, /tostruct)

	; definitions for the primary header
	; version of the fits file

	version = string(input.l1_version + 1, format = '(I02)')

	; creation and acquisition times in utc

	date = date_conv(systime(/julian, /utc), 'FITS')

	obt_beg = fxpar(primary_header, 'OBT_BEG')
	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_string = strreplace(strreplace(strmid(date_beg, 0, 19), '-', ''), ':', '')

	journal, 'UTC time of start acquisition = ' + date_beg

	; 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

	; exposure times

	dit = float(fxpar(metadata_extension_header, 'DIT'))
	telapse = float(obt_end - obt_beg)
	xposure = dit/1000.

	; instrument keywords
	; TODO - complete with filter information

	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'
	telescope = 'SOLO/Metis/' + detector

	; campaign keywords

	obs_id_fields = strsplit(planning_data.obs_id, '_', /extract)
	soop_type = obs_id_fields[2]
	obs_type = obs_id_fields[4]

	; build the fits file extensions

	; join the metadata extension header to the primary header removing unwanted keywords

	i = where(strmid(primary_header, 0, 8) eq 'DATASUM ')
	j = where(strmid(metadata_extension_header, 0, 8) eq 'EXTNAME ')
	k = where(strmid(metadata_extension_header, 0, 8) eq 'TFORM1  ')

	primary_header = [ $
		primary_header[0 : i - 1], $
		metadata_extension_header[j + 1 : k - 1], $
		primary_header[i : *] $
	]
	
	; rebin the image if binning was applied during the acquisition

	if fxpar(primary_header, 'COMPR') then bin_type = fxpar(primary_header, 'B0_BIN') else bin_type = 0
	
	fpfilt = fxpar(primary_header, 'VLFPFILT', count = count)
	if count gt 0 and fpfilt eq 0 and bin_type eq 0 then bin_type = 1
	
	bin_fact = 2. ^ bin_type
	if bin_fact gt 1. then begin
		data = rebin(data, naxis1 / bin_fact, naxis2 / bin_fact, /sample)
		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)')
	endif else begin
		journal, 'Data was not binned during acquistion.'
	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, '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'
	fxaddpar, primary_header, 'DATE-AVG', date_avg, 'Average time of observation', before = 'OBT_BEG'
	fxaddpar, primary_header, 'DATE-END', date_end, 'End time of observation', before = 'OBT_BEG'
	fxaddpar, primary_header, 'TIMESYS', 'UTC', 'System used for time keywords', before = 'OBT_BEG'
	fxaddpar, primary_header, 'TIMRDER', 0.0, 'Estimated random error in time values', before = 'OBT_BEG'
	fxaddpar, primary_header, 'TIMSYER', 0.0, 'Estimated systematic error in time values', before = 'OBT_BEG'
	fxaddpar, primary_header, 'LEVEL', 'L1'
	fxaddpar, primary_header, 'CREATOR', 'metis_l1_prep.pro'
	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'
	fxaddpar, primary_header, 'OBJECT', 'TBD', 'The use of the keyword OBJECT is [TBD]', before = 'DATAMIN'
	fxaddpar, primary_header, 'OBS_MODE', planning_data.obs_mode, 'Observation mode or study that has been used to acquire this image', before = 'DATAMIN'
	fxaddpar, primary_header, 'OBS_TYPE', obs_type, 'Subfield of OBS_ID that only contains an encoded version of OBS_MODE', before = 'DATAMIN'
	fxaddpar, primary_header, 'FILTER', filter, 'Filter used to acquire this image', before = 'DATAMIN'
	fxaddpar, primary_header, 'WAVELNTH', wavelnth, 'Characteristic wavelength at which the observation was taken', before = 'DATAMIN'
	fxaddpar, primary_header, 'WAVEMIN', wavemin, 'The shortest wavelength at which the net (approximate) response function becomes 0.05 times the maximum response. ', before = 'DATAMIN'
	fxaddpar, primary_header, 'WAVEMAX', wavemax, 'The longest wavelength at which the net (approximate) response function becomes 0.05 times the maximum response', before = 'DATAMIN'
	fxaddpar, primary_header, 'WAVEBAND', waveband, 'Description of the wavelength band', before = 'DATAMIN'
	fxaddpar, primary_header, 'XPOSURE', xposure, 'Total effective exposure time of the observation, in seconds', before = 'DATAMIN'
	fxaddpar, primary_header, 'NSUMEXP', 1, 'Number of images summed together to form the observation', before = 'DATAMIN'
	fxaddpar, primary_header, 'TELAPSE', telapse, 'Total elapsed time between the beginning and end of the complete observation in seconds, including any dead times between exposures', before = 'DATAMIN'
	fxaddpar, primary_header, 'SOOPNAME', planning_data.soop_name, 'SOOP(s) that this observation belongs to', before = 'DATAMIN'
	fxaddpar, primary_header, 'SOOPTYPE', soop_type, before = 'DATAMIN'
	fxaddpar, primary_header, 'OBS_ID', planning_data.obs_id, 'Unique identifier for the observation that is associated with the data acquisition', before = 'DATAMIN' 
	fxaddpar, primary_header, 'TARGET', 'TBD', 'Taget as defined in the SOOP description', before = 'DATAMIN'
	fxaddpar, primary_header, 'BSCALE', 1.0, before = 'DATAMIN'
	fxaddpar, primary_header, 'BZERO', 0.0, before = 'DATAMIN'
	fxaddpar, primary_header, 'BTYPE', strupcase(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'
	if datatype le 6 then begin
		fxaddpar, primary_header, 'PXBEG1', 1, 'First pixel that has been read out in dimension 1', before = 'COMPRESS'
		fxaddpar, primary_header, 'PXBEG2', 1, 'First pixel that has been read out in dimension 2', before = 'COMPRESS'
		fxaddpar, primary_header, 'PXEND1', naxis1, 'Last pixel that has been read out in dimension 1', before = 'COMPRESS'
		fxaddpar, primary_header, 'PXEND2', naxis2, 'Last pixel that has been read out in dimension 2', before = 'COMPRESS'
	endif
	fxaddpar, primary_header, 'NBIN1', bin_fact, 'Binning factor in the dimension 1', before = 'COMPRESS'
	fxaddpar, primary_header, 'NBIN2', bin_fact, 'Binning factor in the dimension 2', before = 'COMPRESS'
	fxaddpar, primary_header, 'NBIN', bin_fact * bin_fact, 'Product of all NBIN values above', before = 'COMPRESS'
	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'

	; read the house-keeping telemetry

	hk_table = json_parse(input.hk_file_name, /toarray, /tostruct)

	; replace raw values with calibrated values in the primary header

	if datatype eq 0 or datatype eq 3 or datatype eq 5 then begin
		
		; WARN - DACPOL parameters are not calibrated since a calibration curve does not exist in the IDB. Their calibration in physical units (e.g., voltages or angles) should be done later

		; fxaddpar, primary_header, 'DAC1POL1', interpol_param(hk_table, 'NIT0E061', date_avg, empty_param = empty_param)
		; fxaddpar, primary_header, 'DAC2POL1', interpol_param(hk_table, 'NIT0E062', date_avg, empty_param = empty_param)
		; fxaddpar, primary_header, 'DAC1POL2', interpol_param(hk_table, 'NIT0E064', date_avg, empty_param = empty_param)
		; fxaddpar, primary_header, 'DAC2POL2', interpol_param(hk_table, 'NIT0E065', date_avg, empty_param = empty_param)
		; fxaddpar, primary_header, 'DAC1POL3', interpol_param(hk_table, 'NIT0E067', date_avg, empty_param = empty_param)
		; fxaddpar, primary_header, 'DAC2POL3', interpol_param(hk_table, 'NIT0E068', date_avg, empty_param = empty_param)
		; fxaddpar, primary_header, 'DAC1POL4', interpol_param(hk_table, 'NIT0E06A', date_avg, empty_param = empty_param)
		; fxaddpar, primary_header, 'DAC2POL4', interpol_param(hk_table, 'NIT0E06B', date_avg, empty_param = empty_param)
		
		fxaddpar, primary_header, 'TSENSOR ', interpol_param(hk_table, 'NIT0E0E0', date_avg, empty_params = empty_params)
		fxaddpar, primary_header, 'PMPTEMP ', interpol_param(hk_table, 'NIT0L00D', date_avg, empty_params = empty_params)

		journal, 'Header keywords were calibrated using HK parameters:'
		journal, '  TSENSOR = ' + string(fxpar(primary_header, 'TSENSOR'), format = '(F0)')
		journal, '  PMPTEMP = ' + string(fxpar(primary_header, 'PMPTEMP'), format = '(F0)')
	endif

	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)
		fxaddpar, primary_header, 'HVU_MCP ', interpol_param(hk_table, 'NIT0E071', date_avg, empty_params = empty_params)
		fxaddpar, primary_header, 'TSENSOR ', interpol_param(hk_table, 'NIT0E050', date_avg, empty_params = empty_params)

		journal, 'Header keywords were calibrated using HK parameters:'
		journal, '  HVU_SCR = ' + string(fxpar(primary_header, 'HVU_SCR'), format = '(F0)')
		journal, '  HVU_MCP = ' + string(fxpar(primary_header, 'HVU_MCP'), format = '(F0)')
		journal, '  TSENSOR = ' + string(fxpar(primary_header, 'TSENSOR'), format = '(F0)')
	endif

	if ~ isa(empty_params) then empty_params = ''

	; replace verbose keyword values
	
	key = fxpar(primary_header, 'MEASKIND', count = count)
	if count gt 0 then begin
		pol_id = fxpar(primary_header, 'POL_ID')
		if key eq 0 and pol_id eq 0 then fxaddpar, primary_header, 'MEASKIND', 'Fixed pol.'
		if key eq 0 and pol_id ne 0 then fxaddpar, primary_header, 'MEASKIND', 'pB'
		if key eq 1 then fxaddpar, primary_header, 'MEASKIND', 'tB'
	endif

	key = fxpar(primary_header, 'FRAMEMOD', count = count)
	if count gt 0 then begin
		if key eq 0 then fxaddpar, primary_header, 'FRAMEMOD', 'Single'
		if key eq 1 then fxaddpar, primary_header, 'FRAMEMOD', 'Multiple'
	endif

	key = fxpar(primary_header, 'VLFPFILT', count = count)
	if count gt 0 then begin
		pol_id = fxpar(primary_header, 'POL_ID')
		if pol_id ne 0 then fxaddpar, primary_header, 'VLFPFILT', 'None' else begin
			if key eq 0 then fxaddpar, primary_header, 'VLFPFILT', 'Binning'
			if key eq 1 then fxaddpar, primary_header, 'VLFPFILT', 'Masking'
			if key eq 2 then fxaddpar, primary_header, 'VLFPFILT', 'No filter'
		endelse
	endif

	key = fxpar(primary_header, 'PMPSTAB', count = count)
	if count gt 0 then begin
		if key eq 0 then fxaddpar, primary_header, 'PMPSTAB', 'Not stable'
		if key eq 1 then fxaddpar, primary_header, 'PMPSTAB', 'Stable'
	endif

	key = fxpar(primary_header, 'REF_ROWS', count = count)
	if count gt 0 then begin
		if key eq 0 then fxaddpar, primary_header, 'REF_ROWS', 'Not included'
		if key eq 1 then fxaddpar, primary_header, 'REF_ROWS', 'Included'
	endif

	keys = ['PRESUM', 'CME_OBS', 'SUNDISK', 'CR_SEP', 'SP_NOISE', 'COMPR', 'MASKING', 'RADIAL']
	foreach key_name, keys do begin
		key = fxpar(primary_header, key_name, count = count)
		if count gt 0 then begin
			if key eq 0 then fxaddpar, primary_header, key_name, 'Disabled'
			if key eq 1 then fxaddpar, primary_header, key_name, 'Enabled'
		endif
	endforeach

	; remove unused keywords

	sxdelpar, primary_header, 'WIDTH'
	sxdelpar, primary_header, 'HEIGHT'
	sxdelpar, primary_header, 'X_SIZE'
	sxdelpar, primary_header, 'Y_SIZE'
	sxdelpar, primary_header, 'Z_SIZE'
	sxdelpar, primary_header, 'P_BANDS'
	sxdelpar, primary_header, 'N_BANDS'
	sxdelpar, primary_header, 'ORIG_X'
	sxdelpar, primary_header, 'ORIG_Y'

	; modify keywords for file history

	date = date_conv(systime(/julian, /utc), 'FITS')
	
	old_history = fxpar(primary_header, 'HISTORY')
	sxdelpar, primary_header, 'HISTORY'
	history = ['L1 FITS file created on ' + date, old_history]
	for k = 0, n_elements(history) - 1 do $
		fxaddpar, primary_header, 'HISTORY', history[k]

	; modify keywords for comments

	old_comment = fxpar(primary_header, 'COMMENT', count = count)
	if count eq 0 then old_comment = ''
	sxdelpar, primary_header, 'COMMENT'
	if max(old_comment.startswith('Warning: correction of the OBT for pB acquisitions was applied.')) then begin
		if ~ isa(comment) then comment = !null
		comment = ['Correction of the OBT for pB acquisitions was applied at L0 level.', comment]
	endif
	if max(old_comment.startswith('Warning: the OBT of the acquisition is missing.')) then begin
		if ~ isa(comment) then comment = !null
		comment = ['Warning: the OBT of the acquisition is missing.','  The generation time of the first science packet was used', '  to compute OBT_BEG. Polarimetric keywords DAC*POL*','  may be wrong.', comment]
	endif
	if isa(comment) then begin
		for k = 0, n_elements(comment) - 1 do $
			fxaddpar, primary_header, 'COMMENT', comment[k]
	endif

	; add checksum and datasum to the fits header
	; WARN - should this be done at the end? i don't know

	fits_add_checksum, primary_header, image

	mwrfits, data, out_file_name, primary_header, /no_comment, /create, /silent

	journal, 'Fits file created:'
	journal, '  file name = ' + out_file_name


	; 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

	; 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'

	for i = 0, n_elements(hk_table.par_name) - 1 do begin
		param = create_struct( $
			'PAR_NAME', hk_table.par_name[i], $
			'PACKET', hk_table.packet[i], $
			'GEN_TIME', hk_table.gen_time[i], $
			'REC_TIME', hk_table.rec_time[i], $
			'RAW_VAL', hk_table.raw_val[i], $
			'ENG_VAL', hk_table.eng_val[i], $
			'UNIT', hk_table.unit[i], $
			'DESCR', hk_table.desc[i] $
		)
		if ~ isa(hk_bin_table) then hk_bin_table = param else hk_bin_table = [hk_bin_table, param]
	endfor
 
	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

	journal, 'SPICE kernel files unloaded.'

	; write the auxiliary information file

	output = { $
		file_name: out_file_name, $
		l0_file_name: input.file_name, $
		log_file_name: 'output/metis_l1_prep_log.txt', $
		empty_params : empty_params[uniq(empty_params)]  $
 	}
	
	json_write, output, 'output/contents.json'

	journal, 'Output saved.'

	;close the log

	journal, 'Exiting without errors.'
	journal
	exit, status = 0

	error_handling:
	journal, 'Errors occurred while processing.', /continue
	journal
	exit, status = 1
end