Code optimizations

check_quality.pro

-function check_quality, data, cal_pack, filter = filter
+function check_quality, data, cal_pack, filter
 	
-	if filter.contains('VL', /fold) then $
-		channel = cal_pack.vl_channel else $
-		channel = cal_pack.uv_channel
+	if not keyword_set(filter) then filter = 'VL'
+	
+	if filter.contains('VL', /fold) then channel = cal_pack.vl_channel else channel = cal_pack.uv_channel
 
 	levels = channel.sat_level.value
 

fits_hdr2struct.pro

@@ -7,7 +7,7 @@ function fits_hdr2struct, string_hdr
 		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
+		if isa(struct) then if (tag_names(struct)).hasvalue(tag) then continue
 		struct = create_struct(struct, tag.replace('-', ' '), fxpar(string_hdr, tag))
 	endfor
 

get_light_time.pro

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

interpol_param.pro

 function interpol_param, table, par_name, date, empty_params = empty_params
-	sample = where(table.par_name eq par_name, n)
+	if not keyword_set(empty_params) then empty_params = !null
+	s = where(table.par_name eq par_name, n)
+	if n eq 0 then return, 0.0
+	par_time = table.gen_time[s]
+	par_val = table.eng_val[s]
 	par_date = dblarr(n)
-	for k = 0, n - 1 do par_date[k] = date_conv(table.gen_time[sample[k]], 'JULIAN')
-	if total(table.eng_val[sample].contains('N/A')) then begin
-		table.eng_val[sample] = table.raw_val[sample]
+	for i = 0, n - 1 do par_date[i] = date_conv(par_time[i], 'JULIAN')
+	if max(par_val.contains('N/A')) then begin
 		empty_params = [empty_params, par_name]
+		par_val = table.raw_val[s]
 	endif
-	par_val = float(table.eng_val[sample])
 	s = sort(par_date)
 	par_val = par_val[s]
 	par_date = par_date[s]
 	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 jul_date ge min(par_date) and jul_date le max(par_date) then value = interpol(float(par_val), par_date, jul_date) else value = 0.0
 	if finite(value) then return, value else return, 0.0
 end

metis_l1_prep.pro

@@ -44,6 +44,10 @@ pro metis_l1_prep
 
 	cal_pack = create_struct('path', input.cal_pack_path + path_sep(), cal_pack)
 
+	journal, 'Calibration package correctly imported:'
+	journal, '  version = ' + cal_pack.version
+	journal, '  validity range = ' + string(cal_pack.validity_range.start, cal_pack.validity_range._end, format = '(A, "-", A)')
+
 	; read the primary hdu
 
 	data = mrdfits(input.file_name, 0, primary_header, /silent)
@@ -100,7 +104,8 @@ pro metis_l1_prep
 	date_beg_string = date_beg_string.replace('-', '')
 	date_beg_string = date_beg_string.replace(':', '')
 
-	journal, 'UTC time of start acquisition = ' + date_beg
+	journal, 'Calculated UTC time of start acquisition:'
+	journal, '  time = ' + date_beg
 
 	; name of the fits file
 
@@ -130,11 +135,11 @@ pro metis_l1_prep
 	; instrument/telescope/detector keywords
 
 	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]
+	waveband = channel.name
 	
 	; campaign keywords
 
@@ -177,7 +182,7 @@ pro metis_l1_prep
 
 			; evaluate the quality matrix
 
-			quality_matrix = check_quality(data, cal_pack, filter = filter)
+			quality_matrix = check_quality(data, cal_pack, filter)
 
 			; correct values including binning factor only for vl and uv images and accumulation matrices
 
@@ -191,7 +196,7 @@ pro metis_l1_prep
 			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)
+			quality_matrix = check_quality(data, cal_pack, filter)
 
 			journal, 'Data was not binned during acquistion.'
 		endelse
@@ -222,16 +227,18 @@ pro metis_l1_prep
 
 	if datatype le 2 then begin
 		nsumexp = ndit * ndit1 * ndit2
-		xposure = dit/1000. * nsumexp
+		xposure = dit/1000.D0 * nsumexp
 		data = long(data) * nsumexp
 
 		if ~ isa(comment) then comment = !null
 		comment = [comment, 'Image values were corrected for the total exposure time.']
 
-		journal, 'Image values were corrected for the total exposure time.'
+		journal, 'Image values were corrected for the total exposure time:'
+		journal, '  dit = ' + string(dit, format = '(I0)')
+		journal, '  nsumexp = ' + string(nsumexp, format = '(I0)')
 	endif else begin
 		nsumexp = 1
-		xposure = dit/1000.
+		xposure = dit/1000.D0
 	endelse
 
 	; adjust the primary header (it is almost the same for all data product types)
@@ -398,10 +405,10 @@ pro metis_l1_prep
 
 	; 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'
+	ephemeris = solo_get_ephemeris(header, cal_pack)
+	foreach element, ephemeris do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
 
-	fxaddpar, primary_header, 'HISTORY', 'WCS and solar ephemerides:'
+	fxaddpar, primary_header, 'HISTORY', 'WCS and solar ephemeris:'
 	fxaddpar, primary_header, 'HISTORY', '  SKD version = ' + kernel_version
 
 	; remove redundant and misleading keywords
@@ -469,6 +476,8 @@ pro metis_l1_prep
 		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
+
+		journal, 'Quality-matrix extension correctly added.'
 	endif
 
 	; build the telemetry extension
@@ -481,7 +490,7 @@ pro metis_l1_prep
 	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'
+	journal, 'HK binary-table extension correctly added.'
 
 	; unload the spice kernels
 
@@ -500,16 +509,11 @@ pro metis_l1_prep
 
 	json_write, output, 'output/contents.json'
 
-	journal, 'Output saved.'
+	journal, 'Output saved'
 
 	;close the log
 
 	journal, 'Exiting without errors.'
 	journal
 	exit, status = 0
-
-	error_handling:
-	journal, 'Errors occurred while processing.'
-	journal
-	exit, status = 1
 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
+	if not keyword_set(filter) then filter = 'VL'
+	if filter.contains('VL', /fold) then image = rotate(image, 1) else begin
 		image = reverse(image, 1)
 		image = rotate(image, 1)
-	endif
+	endelse
 	
 	return, image
 end

metis_wcs.pro

@@ -42,7 +42,7 @@ function metis_wcs, header, cal_pack, ref_detector = ref_detector
 
 	; determine spacecract pointing information in the hpc reference frame using the spice kernels
 
-	pointing = solo_get_pointing(header.date_avg)
+	pointing = solo_get_pointing(header.date_avg, /degrees, /arcsec)
 
 	; NOTE - values are defined as follows:
 	; pointing[0] = yaw (arcsec)
@@ -187,14 +187,14 @@ function metis_wcs, header, cal_pack, ref_detector = ref_detector
 		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: '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], $

solo_get_carrot.pro → solo_get_carringrot.pro

-function solo_get_carrot, utc
+function solo_get_carringrot, utc
 
 	; initial estimate of the carrington rotation number
 
 	jul_day = date_conv(utc, 'julian')
 
-	carr = (jul_day - 2398167.D0)/27.2753D0 + 1
+	carr = (jul_day - 2398167.)/27.2753D0 + 1
 
 	; convert the requested date into ephemeris time
 
@@ -12,7 +12,7 @@ function solo_get_carrot, utc
 
 	; get the carrington longitude of earth
 
-	cspice_spkezr, 'EARTH', et, 'IAU_SUN', 'NONE', 'SUN', state, ltime
+	cspice_spkezr, 'EARTH', et, 'IAU_SUN', 'NONE', 'SUN', state, light_time
 
 	; convert rectangular coordinates into angular coordinates
 

solo_get_coords.pro

-function solo_get_coords, utc, frame, obs, velocity = velocity, spherical = spherical, degrees = degrees
+function solo_get_coords, utc, frame, observer, spherical = spherical, degrees = degrees, velocity = velocity
 
 	; convert the requested date into ephemeris time
 
@@ -21,7 +21,7 @@ function solo_get_coords, utc, frame, obs, velocity = velocity, spherical = sphe
 
 	; coordinates of the s/c in the selected reference frame
 
-	cspice_spkezr, 'SOLO', et, frame, 'NONE', obs, state, ltime
+	cspice_spkezr, 'SOLO', et, frame, 'NONE', observer, state, light_time
 
 	coord = state[0 : 2]
 

solo_get_ephemerides.pro → solo_get_ephemeris.pro

-function solo_get_ephemerides, header, constants = constants
+function solo_get_ephemeris, header, cal_pack
+
+	; ephemerides are to be calculated for the average date of observation
 
 	utc = header.date_avg
 
-	rsun = constants.rsun.value
-	au = constants.au.value
+	; get useful physical constant
+
+	rsun = cal_pack.constants.rsun.value
+	au = cal_pack.constants.au.value
 
 	; compute and add the wcs keyword
 
@@ -47,7 +51,7 @@ function solo_get_ephemerides, header, constants = constants
 		value: solar_angles[2], $
 		comment: '[deg] S/C ecliptic North to solar North angle'}
 
-	carrot = solo_get_carrot(utc)
+	carrot = solo_get_carringrot(utc)
 
 	ephemerides.add, { $
 		name: 'CAR_ROT', $
@@ -91,7 +95,8 @@ function solo_get_ephemerides, header, constants = constants
 
 	; coordinates of the S/C in the hee frame
 
-	coord = solo_get_coords(utc, 'HEE', 'SUN') * 1000.
+	coord = solo_get_coords(utc, 'HEE', 'SUN')
+	coord = coord * 1000.
 
 	ephemerides.add, { $
 		name: 'HEEX_OBS', $
@@ -108,9 +113,9 @@ function solo_get_ephemerides, header, constants = constants
 
 	; coordinates of the S/C in the hci frame
 
-	coord = solo_get_coords(utc, 'HCI', 'SUN', velocity = veloc)
+	coord = solo_get_coords(utc, 'HCI', 'SUN', vel = vel)
 	coord = coord * 1000.
-	veloc = veloc * 1000.
+	vel = vel * 1000.
 
 	ephemerides.add, { $
 		name: 'HCIX_OBS', $
@@ -127,15 +132,15 @@ function solo_get_ephemerides, header, constants = constants
 
 	ephemerides.add, { $
 		name: 'HCIX_VOB', $
-		value: veloc[0], $
+		value: vel[0], $
 		comment: '[m/s] S/C Heliocentric Inertial X velocity'}
 	ephemerides.add, { $
 		name: 'HCIY_VOB', $
-		value: veloc[1], $
+		value: vel[1], $
 		comment: '[m/s] S/C Heliocentric Inertial Y velocity'}
 	ephemerides.add, { $
 		name: 'HCIZ_VOB', $
-		value: veloc[2], $
+		value: vel[2], $
 		comment: '[m/s] S/C Heliocentric Inertial Z velocity'}
 
 	; coordinates of the S/C in the hae frame
@@ -195,16 +200,16 @@ function solo_get_ephemerides, header, constants = constants
 	; 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
+	sun_time = get_light_time(utc, 'SOLO', 'SUN', rad_vel = rad_vel)
+	t_del = earth_time - sun_time
 
 	ephemerides.add, { $
 		name: 'OBS_VR', $
-		value: radvel, $
+		value: rad_vel, $
 		comment: '[m/s] Radial velocity of S/C relative to Sun'}
 	ephemerides.add, { $
 		name: 'EAR_TDEL', $
-		value: tdel, $
+		value: t_del, $
 		comment: '[s] Time(Sun to Earth) - Time(Sun to S/C)'}
 	ephemerides.add, { $
 		name: 'SUN_TIME', $
@@ -216,12 +221,12 @@ function solo_get_ephemerides, header, constants = constants
 	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')
+	date_ear = date_conv(jul_utc + t_del/86400., 'fits')
+	date_sun = date_conv(jul_utc - sun_time/86400., 'fits')
 
 	ephemerides.add, { $
 		name: 'DATE_EAR', $
-		value: date_earth, $
+		value: date_ear, $
 		comment: '[UTC] Start time of observation, corrected to Earth'}
 	ephemerides.add, { $
 		name: 'DATE_SUN', $

solo_get_pointing.pro

-function solo_get_pointing, utc, degrees = degrees, radians = radians, celestial = celestial
+function solo_get_pointing, utc, degrees = degrees, arcsec = arcsec
 
 	; convert the requested date into ephemeris time
 
@@ -6,7 +6,7 @@ function solo_get_pointing, utc, degrees = degrees, radians = radians, celestial
 
 	; convert the ephemeris time into s/c internal time
 
-	cspice_sce2c, -144L, et, clk_in
+	cspice_sce2c, -144, et, clk_in
 
 	; set the proper solar orbiter frame
 
@@ -14,7 +14,7 @@ function solo_get_pointing, utc, degrees = degrees, radians = radians, celestial
 
 	if keyword_set(celestial) then frame = 'J2000'
 
-	cspice_ckgp, -144000L, clk_in, 1, frame, cmat, clk_out, found
+	cspice_ckgp, -144000, 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)
 
@@ -30,11 +30,11 @@ function solo_get_pointing, utc, degrees = degrees, radians = radians, celestial
 
 	; correct any cases where the roll is greater than +/- 180 degrees
 
-	if abs(roll) gt !dpi then roll = roll - 2.d0 * !dpi * signum(roll)
+	if abs(roll) gt !dpi then roll = roll - 2. * !dpi * signum(roll)
 
 	; correct any cases where the pitch is greater than +/- 90 degrees
 
-	if abs(pitch) gt !dpi / 2.d0 then begin
+	if abs(pitch) gt !dpi/2. then begin
 		yaw = yaw - !dpi * signum(yaw)
 		pitch = !dpi * signum(pitch) - pitch
 		roll = roll - !dpi * signum(roll)
@@ -42,8 +42,8 @@ function solo_get_pointing, utc, degrees = degrees, radians = radians, celestial
 
 	; 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.
+	if keyword_set(degrees) or keyword_set(arcsec) then rad_factor = 180./!dpi else rad_factor = 1.
+	if keyword_set(arcsec) then arc_factor = 3600. else arc_factor = 1.
 
 	yaw = yaw * rad_factor * arc_factor
 	pitch = pitch * rad_factor * arc_factor

solo_get_solar_angles.pro

@@ -7,11 +7,12 @@ function solo_get_solar_angles, utc
 	cspice_pxform, 'IAU_SUN', 'J2000', et, rot
 	sun_north = rot[2, *]
 
-	cspice_spkezr, 'SUN', et, 'J2000', 'NONE', 'SOLO', state, ltime
+	cspice_spkezr, 'SUN', et, 'J2000', 'NONE', 'SOLO', state, light_time
 	coord = state[0 : 2]
+
 	rad = coord/sqrt(total(coord^2, 1))
 
-	cel_north = [0., 0., 1.d0]
+	cel_north = [0., 0., 1.]
 
 	sun_proj = sun_north - total(rad * sun_north) * rad
 	sun_proj = sun_proj/sqrt(total(sun_proj^2))
@@ -25,13 +26,14 @@ function solo_get_solar_angles, utc
 	x_proj = total(sun_proj * cel_proj)
 	y_proj = total(sun_proj * vec_proj)
 
-	p0 = atan(y_proj, x_proj) * 180.d0 / !dpi
+	p0 = atan(y_proj, x_proj) * 180./!dpi
 
 	cspice_pxform, 'IAU_SUN', 'ECLIPJ2000', et, rot
 	sun_north = rot[2, *]
 
-	cspice_spkezr, 'SUN', et, 'ECLIPJ2000', 'NONE', 'SOLO', state, ltime
+	cspice_spkezr, 'SUN', et, 'ECLIPJ2000', 'NONE', 'SOLO', state, light_time
 	coord = state[0 : 2]
+
 	rad = coord/sqrt(total(coord^2, 1))
 
 	eclip_north = cel_north
@@ -48,14 +50,14 @@ function solo_get_solar_angles, utc
 	x_proj = total(sun_proj * ecl_proj)
 	y_proj = total(sun_proj * vec_proj)
 
-	ep = atan(y_proj, x_proj) * 180.d0 / !dpi
+	ep = atan(y_proj, x_proj) * 180./!dpi
 
-	cspice_spkezr, 'SUN', et, 'IAU_SUN', 'NONE', 'SOLO', state, ltime
+	cspice_spkezr, 'SUN', et, 'IAU_SUN', 'NONE', 'SOLO', state, light_time
 	coord = state[0 : 2]
 
 	cspice_reclat, coord, dist, lon, lat
 
-	b0 = -lat * 180.d0 / !dpi
+	b0 = -lat * 180./!dpi
 
 	return, [b0, p0, ep]
 end