diff --git a/check_quality.pro b/check_quality.pro
index c495f783d99fbf3b755fed937c5242ac98232002..55006761b35de077374a6ae34c69e71bd283caad 100644
--- a/check_quality.pro
+++ b/check_quality.pro
@@ -1,8 +1,8 @@
-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
+function check_quality, data, cal_pack, filter
+
+ 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
diff --git a/fits_hdr2struct.pro b/fits_hdr2struct.pro
index db742027a484474c4add1bfdc8dfc2213d0c695f..e55516d39d52b5a183bc2ea1148d98061f21b414 100644
--- a/fits_hdr2struct.pro
+++ b/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
diff --git a/get_light_time.pro b/get_light_time.pro
index 4146dcf52a94fe5b330e6fd9163294606a116449..3e3c794b6089b484bac40fe4f64ed72866ac0765 100644
--- a/get_light_time.pro
+++ b/get_light_time.pro
@@ -1,16 +1,16 @@
-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
diff --git a/interpol_param.pro b/interpol_param.pro
index 839cf0c8d8eda3b124a7e8e9c02bc027b0321288..350aaf96cd1ddf9838309e39591f726902be6ae7 100644
--- a/interpol_param.pro
+++ b/interpol_param.pro
@@ -1,16 +1,19 @@
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
diff --git a/metis_l1_prep.pro b/metis_l1_prep.pro
old mode 100644
new mode 100755
index f967ae0f3386bde087482ffbf6838f59d0450c60..6e1492977bccb24d5cf0131e8ea41083e1cfb5ed
--- a/metis_l1_prep.pro
+++ b/metis_l1_prep.pro
@@ -34,7 +34,7 @@ pro metis_l1_prep
load_spice_kernels, input.spice_kernels, kernel_list = kernel_list, kernel_version = kernel_version
journal, 'SPICE kernel files correctly loaded:'
- journal, ' SDK version= ' + kernel_version
+ journal, ' SDK version = ' + kernel_version
; import the calibration package index file
@@ -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,12 +135,12 @@ 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
obs_id_fields = strsplit(planning_data.obs_id, '_', /extract)
@@ -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
diff --git a/metis_rectify.pro b/metis_rectify.pro
index 327b06025ea3a622041f6e38d825cf8de88798c6..23db709fbcdc16a692a548ad62018a33f188009f 100644
--- a/metis_rectify.pro
+++ b/metis_rectify.pro
@@ -1,9 +1,9 @@
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
diff --git a/metis_wcs.pro b/metis_wcs.pro
index a7e1edf796b1d73cab9f9bd68455c2ac03b7dc96..74b4d427e25f808e3f8d08f090de4a7793e482e4 100644
--- a/metis_wcs.pro
+++ b/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)
@@ -58,7 +58,7 @@ function metis_wcs, header, cal_pack, ref_detector = ref_detector
; 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)]]
@@ -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], $
diff --git a/solo_get_carrot.pro b/solo_get_carringrot.pro
similarity index 84%
rename from solo_get_carrot.pro
rename to solo_get_carringrot.pro
index 86127297f2094f7b17f632e04fbc8b3a597c064d..7ce924886fefff87dcd56f229b9f4ab01114d6f7 100644
--- a/solo_get_carrot.pro
+++ b/solo_get_carringrot.pro
@@ -1,10 +1,10 @@
-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
diff --git a/solo_get_coords.pro b/solo_get_coords.pro
index b8d2ec462928acb678666cad075d9d58691e8301..dab61d4c22fa0d8c8df553daadad9b3051c0b7fd 100644
--- a/solo_get_coords.pro
+++ b/solo_get_coords.pro
@@ -1,4 +1,4 @@
-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]
diff --git a/solo_get_ephemerides.pro b/solo_get_ephemeris.pro
similarity index 88%
rename from solo_get_ephemerides.pro
rename to solo_get_ephemeris.pro
index ca3d91ce842f18a38002f01711d01552af326e8f..58a18f6e598ca65430e3e77ab1d6e2f17ab1a86b 100644
--- a/solo_get_ephemerides.pro
+++ b/solo_get_ephemeris.pro
@@ -1,9 +1,13 @@
-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', $
diff --git a/solo_get_pointing.pro b/solo_get_pointing.pro
index e18ef4601650d755aa60662746556dbc37ea0da0..b3449268b2aeeab6d74713f3a6064b3818d9ab41 100644
--- a/solo_get_pointing.pro
+++ b/solo_get_pointing.pro
@@ -1,4 +1,4 @@
-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
diff --git a/solo_get_solar_angles.pro b/solo_get_solar_angles.pro
index 2538d26c4069688ecebffd0bf74fd1c99a1b76ff..217e744cfdf1e176d2b2ccd5c36e7121ef2cf438 100644
--- a/solo_get_solar_angles.pro
+++ b/solo_get_solar_angles.pro
@@ -7,55 +7,57 @@ 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]
+ rad = coord/sqrt(total(coord^2, 1))
+
+ cel_north = [0., 0., 1.]
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))
cel_proj = cel_north - total(rad * cel_north) * rad
- cel_proj = cel_proj / sqrt(total(cel_proj^2))
+ cel_proj = cel_proj/sqrt(total(cel_proj^2))
vec_proj = crossp(cel_north, rad)
- vec_proj = vec_proj / sqrt(total(vec_proj^2))
+ 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
+ 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))
+
+ 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))
+ 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))
+ ecl_proj = ecl_proj/sqrt(total(ecl_proj^2))
vec_proj = crossp(eclip_north, rad)
- vec_proj = vec_proj / sqrt(total(vec_proj^2))
+ 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
+ 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