Skip to content
Snippets Groups Projects
Commit 84140d2e authored by Roberto Susino's avatar Roberto Susino
Browse files

Code optimizations

parent 5f6290e3
Branches
Tags
No related merge requests found
...@@ -7,7 +7,7 @@ function fits_hdr2struct, string_hdr ...@@ -7,7 +7,7 @@ function fits_hdr2struct, string_hdr
tag = tag.trim() tag = tag.trim()
if tag eq '' then continue if tag eq '' then continue
if tag.startswith('end', /fold) or tag.contains('continue', /fold) 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)) struct = create_struct(struct, tag.replace('-', ' '), fxpar(string_hdr, tag))
endfor endfor
......
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 ; convert the requested date into ephemeris time
cspice_str2et, utc, et 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 ; 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 end
...@@ -126,6 +126,7 @@ pro metis_l2_prep_uv ...@@ -126,6 +126,7 @@ pro metis_l2_prep_uv
fxaddpar, primary_header, 'BUNIT', bunit fxaddpar, primary_header, 'BUNIT', bunit
fxaddpar, primary_header, 'DATAMIN', min(data, /nan) fxaddpar, primary_header, 'DATAMIN', min(data, /nan)
fxaddpar, primary_header, 'DATAMAX', max(data, /nan) fxaddpar, primary_header, 'DATAMAX', max(data, /nan)
fxaddpar, primary_header, 'WAVEBAND', cal_pack.uv_channel.name
fxaddpar, primary_header, 'XPOSURE', header.xposure fxaddpar, primary_header, 'XPOSURE', header.xposure
fxaddpar, primary_header, 'NSUMEXP', header.nsumexp fxaddpar, primary_header, 'NSUMEXP', header.nsumexp
...@@ -136,10 +137,10 @@ pro metis_l2_prep_uv ...@@ -136,10 +137,10 @@ pro metis_l2_prep_uv
; append solar ephemeris keywords ; append solar ephemeris keywords
ephemerides = solo_get_ephemerides(header, constants = cal_pack.constants) ephemeris = solo_get_ephemeris(header, cal_pack)
foreach element, ephemerides do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE' foreach element, ephemeris do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
history = [history, 'Solar ephemerides and WCS:', ' SKD version = ' + kernel_version] history = [history, 'Update WCS and solar ephemeris:', ' SKD version = ' + kernel_version]
; delete useless keywords ; delete useless keywords
......
...@@ -120,6 +120,7 @@ pro metis_l2_prep_vl_generic ...@@ -120,6 +120,7 @@ pro metis_l2_prep_vl_generic
fxaddpar, primary_header, 'BUNIT', bunit fxaddpar, primary_header, 'BUNIT', bunit
fxaddpar, primary_header, 'DATAMIN', min(data, /nan) fxaddpar, primary_header, 'DATAMIN', min(data, /nan)
fxaddpar, primary_header, 'DATAMAX', max(data, /nan) fxaddpar, primary_header, 'DATAMAX', max(data, /nan)
fxaddpar, primary_header, 'WAVEBAND', cal_pack.vl_channel.name
fxaddpar, primary_header, 'XPOSURE', header.xposure fxaddpar, primary_header, 'XPOSURE', header.xposure
fxaddpar, primary_header, 'NSUMEXP', header.nsumexp fxaddpar, primary_header, 'NSUMEXP', header.nsumexp
...@@ -130,10 +131,10 @@ pro metis_l2_prep_vl_generic ...@@ -130,10 +131,10 @@ pro metis_l2_prep_vl_generic
; append solar ephemeris keywords ; append solar ephemeris keywords
ephemerides = solo_get_ephemerides(header, constants = cal_pack.constants) ephemeris = solo_get_ephemeris(header, cal_pack)
foreach element, ephemerides do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE' foreach element, ephemeris do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
history = [history, 'Solar ephemerides and WCS:', ' SKD version = ' + kernel_version] history = [history, 'Update WCS and solar ephemeris:', ' SKD version = ' + kernel_version]
; delete useless keywords ; delete useless keywords
......
...@@ -224,6 +224,7 @@ pro metis_l2_prep_vl_polariz ...@@ -224,6 +224,7 @@ pro metis_l2_prep_vl_polariz
fxaddpar, primary_header, 'OBT_BEG', header.obt_beg, 'Start acquisition time in on-board time', format = 'F0.5' 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, '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, 'TELAPSE', telapse, '[s] Elapsed time between beginning and end of observation'
fxaddpar, primary_header, 'WAVEBAND', cal_pack.vl_channel.name
fxaddpar, primary_header, 'XPOSURE', header.xposure fxaddpar, primary_header, 'XPOSURE', header.xposure
fxaddpar, primary_header, 'NSUMEXP', header.nsumexp fxaddpar, primary_header, 'NSUMEXP', header.nsumexp
fxaddpar, primary_header, 'TSENSOR', header.tsensor fxaddpar, primary_header, 'TSENSOR', header.tsensor
...@@ -236,10 +237,11 @@ pro metis_l2_prep_vl_polariz ...@@ -236,10 +237,11 @@ pro metis_l2_prep_vl_polariz
; append solar ephemeris keywords ; append solar ephemeris keywords
ephemerides = solo_get_ephemerides(header, constants = cal_pack.constants) ephemeris = solo_get_ephemeris(header, cal_pack)
foreach element, ephemerides do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE' foreach element, ephemeris do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
history = ['Update WCS and solar ephemeris:', ' SKD version = ' + kernel_version]
history = ['Solar ephemerides and WCS:', ' SKD version = ' + kernel_version]
tb_history = [tb_history, history] tb_history = [tb_history, history]
pb_history = [pb_history, history] pb_history = [pb_history, history]
......
...@@ -17,7 +17,7 @@ function metis_rad_cal, data, header, cal_pack, polarimetric = polarimetric, his ...@@ -17,7 +17,7 @@ function metis_rad_cal, data, header, cal_pack, polarimetric = polarimetric, his
rad_info = cal_pack.uv_channel.radiometry[0] rad_info = cal_pack.uv_channel.radiometry[0]
units = cal_pack.uv_channel.cal_units units = cal_pack.uv_channel.cal_units
pmp_factor = 1.D0 pmp_factor = 1.D0
unit_factor = 1. unit_factor = 1.D0
endif endif
rad_factor = rad_info.rad_response.value * pmp_factor * $ rad_factor = rad_info.rad_response.value * pmp_factor * $
......
function metis_rectify, image, filter function metis_rectify, image, filter
if filter.toupper() eq 'VL' then image = rotate(image, 1) if not keyword_set(filter) then filter = 'VL'
if filter.toupper() eq 'UV' then begin if filter.contains('VL', /fold) then image = rotate(image, 1) else begin
image = reverse(image, 1) image = reverse(image, 1)
image = rotate(image, 1) image = rotate(image, 1)
endif endelse
return, image return, image
end end
function solo_get_carrot, utc function solo_get_carringrot, utc
; initial estimate of the carrington rotation number ; initial estimate of the carrington rotation number
jul_day = date_conv(utc, 'julian') 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 ; convert the requested date into ephemeris time
...@@ -12,7 +12,7 @@ function solo_get_carrot, utc ...@@ -12,7 +12,7 @@ function solo_get_carrot, utc
; get the carrington longitude of earth ; 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 ; convert rectangular coordinates into angular coordinates
......
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 ; convert the requested date into ephemeris time
...@@ -21,7 +21,7 @@ function solo_get_coords, utc, frame, obs, velocity = velocity, spherical = sphe ...@@ -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 ; 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] coord = state[0 : 2]
......
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 utc = header.date_avg
rsun = constants.rsun.value ; get useful physical constant
au = constants.au.value
rsun = cal_pack.constants.rsun.value
au = cal_pack.constants.au.value
; compute and add the wcs keyword ; compute and add the wcs keyword
...@@ -47,7 +51,7 @@ function solo_get_ephemerides, header, constants = constants ...@@ -47,7 +51,7 @@ function solo_get_ephemerides, header, constants = constants
value: solar_angles[2], $ value: solar_angles[2], $
comment: '[deg] S/C ecliptic North to solar North angle'} comment: '[deg] S/C ecliptic North to solar North angle'}
carrot = solo_get_carrot(utc) carrot = solo_get_carringrot(utc)
ephemerides.add, { $ ephemerides.add, { $
name: 'CAR_ROT', $ name: 'CAR_ROT', $
...@@ -91,7 +95,8 @@ function solo_get_ephemerides, header, constants = constants ...@@ -91,7 +95,8 @@ function solo_get_ephemerides, header, constants = constants
; coordinates of the S/C in the hee frame ; 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, { $ ephemerides.add, { $
name: 'HEEX_OBS', $ name: 'HEEX_OBS', $
...@@ -108,9 +113,9 @@ function solo_get_ephemerides, header, constants = constants ...@@ -108,9 +113,9 @@ function solo_get_ephemerides, header, constants = constants
; coordinates of the S/C in the hci frame ; 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. coord = coord * 1000.
veloc = veloc * 1000. vel = vel * 1000.
ephemerides.add, { $ ephemerides.add, { $
name: 'HCIX_OBS', $ name: 'HCIX_OBS', $
...@@ -127,15 +132,15 @@ function solo_get_ephemerides, header, constants = constants ...@@ -127,15 +132,15 @@ function solo_get_ephemerides, header, constants = constants
ephemerides.add, { $ ephemerides.add, { $
name: 'HCIX_VOB', $ name: 'HCIX_VOB', $
value: veloc[0], $ value: vel[0], $
comment: '[m/s] S/C Heliocentric Inertial X velocity'} comment: '[m/s] S/C Heliocentric Inertial X velocity'}
ephemerides.add, { $ ephemerides.add, { $
name: 'HCIY_VOB', $ name: 'HCIY_VOB', $
value: veloc[1], $ value: vel[1], $
comment: '[m/s] S/C Heliocentric Inertial Y velocity'} comment: '[m/s] S/C Heliocentric Inertial Y velocity'}
ephemerides.add, { $ ephemerides.add, { $
name: 'HCIZ_VOB', $ name: 'HCIZ_VOB', $
value: veloc[2], $ value: vel[2], $
comment: '[m/s] S/C Heliocentric Inertial Z velocity'} comment: '[m/s] S/C Heliocentric Inertial Z velocity'}
; coordinates of the S/C in the hae frame ; coordinates of the S/C in the hae frame
...@@ -195,16 +200,16 @@ function solo_get_ephemerides, header, constants = constants ...@@ -195,16 +200,16 @@ function solo_get_ephemerides, header, constants = constants
; light travel times and radial velocity of the S/C ; light travel times and radial velocity of the S/C
earth_time = get_light_time(utc, 'EARTH', 'SUN') earth_time = get_light_time(utc, 'EARTH', 'SUN')
sun_time = get_light_time(utc, 'SOLO', 'SUN', radvel = radvel) sun_time = get_light_time(utc, 'SOLO', 'SUN', rad_vel = rad_vel)
tdel = earth_time - sun_time t_del = earth_time - sun_time
ephemerides.add, { $ ephemerides.add, { $
name: 'OBS_VR', $ name: 'OBS_VR', $
value: radvel, $ value: rad_vel, $
comment: '[m/s] Radial velocity of S/C relative to Sun'} comment: '[m/s] Radial velocity of S/C relative to Sun'}
ephemerides.add, { $ ephemerides.add, { $
name: 'EAR_TDEL', $ name: 'EAR_TDEL', $
value: tdel, $ value: t_del, $
comment: '[s] Time(Sun to Earth) - Time(Sun to S/C)'} comment: '[s] Time(Sun to Earth) - Time(Sun to S/C)'}
ephemerides.add, { $ ephemerides.add, { $
name: 'SUN_TIME', $ name: 'SUN_TIME', $
...@@ -216,12 +221,12 @@ function solo_get_ephemerides, header, constants = constants ...@@ -216,12 +221,12 @@ function solo_get_ephemerides, header, constants = constants
utc = header.date_beg utc = header.date_beg
jul_utc = date_conv(utc, 'julian') jul_utc = date_conv(utc, 'julian')
date_earth = date_conv(jul_utc + tdel / 86400.d0, 'fits') date_ear = date_conv(jul_utc + t_del/86400., 'fits')
date_sun = date_conv(jul_utc - sun_time / 86400.d0, 'fits') date_sun = date_conv(jul_utc - sun_time/86400., 'fits')
ephemerides.add, { $ ephemerides.add, { $
name: 'DATE_EAR', $ name: 'DATE_EAR', $
value: date_earth, $ value: date_ear, $
comment: '[UTC] Start time of observation, corrected to Earth'} comment: '[UTC] Start time of observation, corrected to Earth'}
ephemerides.add, { $ ephemerides.add, { $
name: 'DATE_SUN', $ name: 'DATE_SUN', $
......
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 ; convert the requested date into ephemeris time
...@@ -6,7 +6,7 @@ function solo_get_pointing, utc, degrees = degrees, radians = radians, celestial ...@@ -6,7 +6,7 @@ function solo_get_pointing, utc, degrees = degrees, radians = radians, celestial
; convert the ephemeris time into s/c internal time ; 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 ; set the proper solar orbiter frame
...@@ -14,7 +14,7 @@ function solo_get_pointing, utc, degrees = degrees, radians = radians, celestial ...@@ -14,7 +14,7 @@ function solo_get_pointing, utc, degrees = degrees, radians = radians, celestial
if keyword_set(celestial) then frame = 'J2000' 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) 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 ...@@ -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 ; 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 ; 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) yaw = yaw - !dpi * signum(yaw)
pitch = !dpi * signum(pitch) - pitch pitch = !dpi * signum(pitch) - pitch
roll = roll - !dpi * signum(roll) roll = roll - !dpi * signum(roll)
...@@ -42,8 +42,8 @@ function solo_get_pointing, utc, degrees = degrees, radians = radians, celestial ...@@ -42,8 +42,8 @@ function solo_get_pointing, utc, degrees = degrees, radians = radians, celestial
; apply correct units to the pointing vector ; apply correct units to the pointing vector
if keyword_set(radians) then rad_factor = 1. else rad_factor = 180.d0 / !dpi if keyword_set(degrees) or keyword_set(arcsec) then rad_factor = 180. / !dpi else rad_factor = 1.
if keyword_set(radians) or keyword_set(degrees) then arc_factor = 1. else arc_factor = 3600. if keyword_set(arcsec) then arc_factor = 3600. else arc_factor = 1.
yaw = yaw * rad_factor * arc_factor yaw = yaw * rad_factor * arc_factor
pitch = pitch * rad_factor * arc_factor pitch = pitch * rad_factor * arc_factor
......
...@@ -7,11 +7,12 @@ function solo_get_solar_angles, utc ...@@ -7,11 +7,12 @@ function solo_get_solar_angles, utc
cspice_pxform, 'IAU_SUN', 'J2000', et, rot cspice_pxform, 'IAU_SUN', 'J2000', et, rot
sun_north = rot[2, *] 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] coord = state[0 : 2]
rad = coord/sqrt(total(coord^2, 1)) 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_north - total(rad * sun_north) * rad
sun_proj = sun_proj/sqrt(total(sun_proj^2)) sun_proj = sun_proj/sqrt(total(sun_proj^2))
...@@ -25,13 +26,14 @@ function solo_get_solar_angles, utc ...@@ -25,13 +26,14 @@ function solo_get_solar_angles, utc
x_proj = total(sun_proj * cel_proj) x_proj = total(sun_proj * cel_proj)
y_proj = total(sun_proj * vec_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 cspice_pxform, 'IAU_SUN', 'ECLIPJ2000', et, rot
sun_north = rot[2, *] 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] coord = state[0 : 2]
rad = coord/sqrt(total(coord^2, 1)) rad = coord/sqrt(total(coord^2, 1))
eclip_north = cel_north eclip_north = cel_north
...@@ -48,14 +50,14 @@ function solo_get_solar_angles, utc ...@@ -48,14 +50,14 @@ function solo_get_solar_angles, utc
x_proj = total(sun_proj * ecl_proj) x_proj = total(sun_proj * ecl_proj)
y_proj = total(sun_proj * vec_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] coord = state[0 : 2]
cspice_reclat, coord, dist, lon, lat cspice_reclat, coord, dist, lon, lat
b0 = -lat * 180.D0/!dpi b0 = -lat * 180./!dpi
return, [b0, p0, ep] return, [b0, p0, ep]
end end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment