"""Geometry module."""
from gfinder.target import Target
from gfinder.detector import Detector
import gfinder.simulator as simulator
import numpy as np
import spiceypy as spice
import spiceypy.utils.support_types as stypes
import geojson
from geojson import Polygon, Point, Feature, MultiPolygon
import inspect
import copy
import re
# Constants
IFOV = 0.00015 # radians
MAJIS_SCAN_ANGLE_MAX = 4.0 # degrees
ABCORR = 'LT+S'
OBSRVR = 'JUICE'
JUPITER_RADIUS = 71492.0 # Jupiter equatorial radius
RINGS_RADIUS_MIN = JUPITER_RADIUS # 1.8 * JUPITER_RADIUS # main ring radius
RINGS_RADIUS_MAX = 100*JUPITER_RADIUS # 3.5 * JUPITER_RADIUS # Gossamer ring radius
# Turning SpiceyPy found check off
spice.found_check_off()
[docs]class Geometry:
"""Abstract class that represents any types of geometrical data associated to an instant acquisition (measurement), an observation, or a sequence.
Args:
name (str): name of the Geometry class
geoevt (GeometryEvent): associated GeometryEvent object (passed as optional argument to the constructor)
searchable (boolean): whether or not this class of Geometry object is searchable
required_parameters_sets (list): list of required parameters
parameters (dict): parameters required to compute geometry data
dimension (int): data array dimension (default is 0, scalar)
units (float): data units
prefix (float): prefix to be used for naming each data array element. For example, SubSC_Point, which is 2-element vector data, prefixes are defined to be ['lon', 'lat'] (see :class:`SubSC_Point`).
data (NoneType): data value, which can be a scalar, vector, or multidimentional array.
valid (boolean): whether or not Geometry object is valid
"""
def __init__(self, parameters, geoevt=None, silent=None):
"""Constructor class.
Args:
parameters (dict): parameters to be used for computation
geoevt(GeometryEvent): associated GeometryEvent object (None by default)
Returns:
Geometry: Geometry object
"""
self.name = self.__class__.__name__
self.geoevt = geoevt
self.searchable = False # TODO: check if it is used
self.required_parameters_sets = [[]]
self.parameters_set_index = -1
self.parameters = {}
self.condition = None
self.dimension = 0
self.units = ''
self.data = None
self.prefix = ''
self.valid = False
self.silent = silent
def __repr__(self):
return '<%s %r>' % (self.__class__.__name__, self.getDict())
[docs] def get_parent_class(self):
return self.__class__.__bases__[0].__name__
[docs] def check_validity(self):
"""Check overall validity of Geometry object.
Checked when not abstract objects are created.
"""
if self.valid:
if self.dimension == 0: # both prefix and units must be of type str.
if not isinstance(self.prefix, str) or not isinstance(self.units, str):
print('WARNING: `prefix` and `units` must be of type str.')
self.valid = False
elif self.dimension == 1:
if self.data is None:
print('WARNING: Default `data` not set.')
self.valid = False
if isinstance(self.prefix, str):
# set default prefixes
self.prefix = []
for i in range(len(self.data)):
self.prefix.append(str(i))
if isinstance(self.units, list):
if len(self.prefix) != len(self.units):
print(f'WARNING: Discrepancy between the number of prefixes and units, respectively: {len(self.prefix)} and {len(self.units)}.')
self.valid = False
elif isinstance(self.units, str):
# set default prefix for all list element
units = self.units
self.units = []
for i in range(len(self.prefix)):
self.units.append(units)
# else:
# print(self.name, self.dimension)
[docs] def check_parameters_set(self, parameters, required_parameters_set):
"""Check input parameters against a set of required parameters.
Args:
parameters (dict): input parameters to be checked.
required_parameters_set (list):
Returns:
valid (bool), missing_parameters (list), extra_parameters (list)
"""
# Init variables
missing_parameters = []
extra_parameters = []
valid = True
# check that each required parameter is available in input parameters
for required_parameter_name in required_parameters_set:
if required_parameter_name not in parameters.keys():
missing_parameters.append(required_parameter_name)
valid = False
# check that there are no extra parameters
for parameter in parameters:
if parameter not in required_parameters_set:
extra_parameters.append(parameter)
valid = False
return valid, missing_parameters, extra_parameters
[docs] def set_parameters(self, parameters):
"""Check that required input parameters are provided, and set input parameters.
Args:
parameters (dict): input parameters to be used for computation.
"""
valid = False
self.parameters_set_index = -1
missing_parameters_sets = []
extra_parameters_sets = []
# check against each required parameters set
for i, required_parameters_set in enumerate(self.required_parameters_sets):
valid, missing_parameters, extra_parameters = self.check_parameters_set(parameters, required_parameters_set)
missing_parameters_sets.append(missing_parameters)
extra_parameters_sets.append(extra_parameters)
if valid:
self.parameters_set_index = i
valid = True
break
# Set parameters if valid
if valid:
self.parameters = parameters
self.valid = True
else:
if not self.silent:
print(f'WARNING: Invalid input {self.name} parameters: {parameters}')
print()
for i in range(len(self.required_parameters_sets)):
print(f'Required parameters set #{i}: {self.required_parameters_sets[i]}')
if missing_parameters_sets[i]:
print(f' - missing parameters: {missing_parameters_sets[i]}')
if extra_parameters_sets[i]:
print(f' - extra parameters: {extra_parameters_sets[i]}')
print()
self.valid = False
# check overall validity of Geoemetry object
self.check_validity()
[docs] def getDict(self):
"""Returns Geometry object.
Returns:
dict: Geometry object time_inputs_dict
"""
geometry_dict = {
'name': self.name,
'parameters': self.parameters,
'value': self.data
}
return geometry_dict
[docs] def get_csv_line(self):
line = ''
if self.dimension == 0:
if self.data is None:
line = ''
else:
line = str(self.data)
elif self.dimension == 1:
values = []
data = self.data
if isinstance(self.data[0], list): # workaround when forcing dimension=1 for single-vector Director and Geolocation children classes
data = self.data[0]
for value in data:
if value is None:
values.append('')
else:
values.append(str(value))
line = ','.join(values)
elif self.dimension >= 2:
line = None
return line
[docs] def get_time_step(self):
time_step = 1.0 # default if no associated GeometryEvent object
if self.geoevt:
time_step = self.geoevt.parent_event.time_step
return time_step
[docs] def compute(self):
"""Compute this Geometry object data at ephemeris time `et`.
"""
pass
[docs] def getData(self):
return self.data
[docs] def getGeoJSON(self):
return None
[docs] def isSearchable(self):
return self.searchable
[docs] def getName(self):
return self.name
[docs] def getCondition(self):
return self.condition
[docs] def is_valid(self):
return self.valid
[docs]class Condition:
def __init__(self, condition):
self.relate = condition['relationalCondition']
self.refval = condition['referenceValue']
self.adjust = condition['adjustmentValue']
def __repr__(self):
return '<%s %r>' % (self.__class__.__name__, self.__dict__)
[docs]class GeometryFactory:
"""Class responsible for creating Geometry objects from input geometry definition.
Args:
classes_index (dict): index of all Geometry classes available.
excluded_class_names (List(str)):
lowlevel_class_names (List(str)):
highlevel_class_names (List(str)):
"""
def __init__(self):
"""Constructor method.
"""
# Make a copy of the current global symbol table (data structure maintained by a compiler which contains all
# necessary information about the program.
globals_copy = copy.copy(globals())
# Set excluded, low-level and high-level Geometry class names
self.excluded_class_names = [
'Geometry',
'Condition',
'GeometryFactory',
'Scalar',
'Vector',
'Multidimentional',
'SimulatedSCFrameError'
]
self.lowlevel_class_names = ['Scalar', 'Vector', 'Multidimentional']
self.highlevel_class_names = ['Direction', 'Geolocation']
self.classes_index = {}
for this_class_name in globals_copy.keys():
if this_class_name not in self.excluded_class_names:
this_class = globals_copy[this_class_name]
if re.search("^(<class 'gfinder.geometry.*'>)", str(this_class)):
self.classes_index.update({this_class_name: this_class})
def __repr__(self):
return '<%s %r>' % (self.__class__.__name__, self.__dict__)
[docs] def get_geometries(self, name_filter='', class_filter='', searchable=False):
"""Returns Geometry class names matching input query.
Args:
searchable:
name_filter:
class_filter:
Returns:
"""
geometry_class_names = []
for geometry_class_name in self.classes_index:
if searchable:
# print(f'is {geometry_class_name} searchable?')
print(self.get_geometry_info(geometry_class_name)['is_searchable'])
if self.get_geometry_info(geometry_class_name)['is_searchable']:
geometry_class_names.append(geometry_class_name)
elif class_filter:
if class_filter == self.get_geometry_info(geometry_class_name)['parent_class']:
geometry_class_names.append(geometry_class_name)
else:
geometry_class_names.append(geometry_class_name)
if name_filter:
geometry_class_names_copy = copy.copy(geometry_class_names)
for geometry_class_name in geometry_class_names_copy:
if name_filter.lower() not in geometry_class_name.lower():
geometry_class_names.remove(geometry_class_name)
return geometry_class_names
[docs] def get_geometry_info(self, class_name):
"""Return information about a Geometry class.
Args:
class_name:
Returns:
geometry_class_dict
"""
geometry_class_dict = {}
if class_name in self.classes_index.keys():
# Create Geometry object corresponding to input class name
geometry = self.classes_index[class_name]({}, silent=True) # could use the create() method, but geometry_def is required.
short_description = ''
if geometry.__doc__:
short_description = geometry.__doc__.split('\n')[0]
geometry_class_dict = dict(
parent_class=geometry.__class__.__bases__[0].__name__, # 'Scalar', 'Vector' or 'Multidimensional'
is_searchable=geometry.isSearchable(),
required_parameters_sets=geometry.required_parameters_sets,
prefixes=geometry.prefix,
units=geometry.units,
short_description=short_description,
source=inspect.getsource(self.classes_index[class_name].compute)
)
return geometry_class_dict
[docs] def get_geometry_required_parameters(self, class_name, set_index=0):
"""Return a Geometry class required parameters set dict."""
if class_name in self.classes_index.keys():
# Create Geometry object corresponding to input class name
geometry = self.classes_index[class_name]({}, silent=True) # could use the create() method, but geometry_def is required.
if set_index < len(geometry.required_parameters_sets):
required_parameters_set = geometry.required_parameters_sets[set_index]
return dict.fromkeys(required_parameters_set, None)
# TODO: Below is an example of what would be a much better definition for required parameters sets.
# self.required_parameters_sets = [
# {'start_angle': 0.0, 'start_scan_rate': 0.0, 'stop_scan_rate': 0.0, 'method': ''},
# {'start_angle': 0.0, 'stop_angle': 0.0, 'start_scan_rate': 0.0, 'stop_scan_rate': 0.0, 'method':''},
# {'start_angle': 0.0, 'scan_rate': 0.0}
# ]
else:
raise IndexError(f'Input set index exceeds number of required parameters sets: {len(geometry.required_parameters_sets)}')
else:
raise AttributeError(f'Input `{class_name}` Geometry class name does not exist.')
[docs] def create(self, geometry_def, geoevt=None):
name = geometry_def['name']
parameters = geometry_def['parameters']
condition = None
if 'condition' in geometry_def.keys():
condition = Condition(geometry_def['condition'])
if name in self.classes_index.keys():
geometry_class = self.classes_index[name]
geometry = geometry_class(parameters, geoevt=geoevt)
# Workaround to assign condition to Geometry object.
# TODO: Geometry classes must be re-designed/implemented to take into account for condition.
# (this should be done in concert with opportunity./sequence.search() methods)
geometry.condition = condition
return geometry
else:
geometry = Geometry({})
geometry.valid = False
print(f'[WARNING] {name} Geometry class name does not exist. Add to geometry.py modules or check for valid input name.')
return geometry
[docs]class Scalar(Geometry):
"""Abstract Geometry class that holds scalar data.
Args:
type (str): scalar type: ``numerical`` (default) or ``boolean``
"""
def __init__(self, parameters, geoevt=None, silent=None):
"""Scalar constructor class.
Returns:
Scalar: Scalar object
"""
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.type = 'numerical' # numerical or boolean
self.dimension = 0
self.valid = True
def __format__(self, geometry_fmt):
if self.data is not None:
str = '{:f}'.format(self.data)
else:
str = '{:f}'.format(0)
return str
[docs] def getData(self):
return self.data
[docs] def search(self, condition, cnfine, step): # cnfine is search_window
"""Search and return times when input geometric condition is met.
Args:
condition: input condition that relates to the scalar quantity.
cnfine (SpiceCell): SPICE window to which the search is restricted.
step (float): Step size used for locating extrema and roots.
Returns:
SpiceCell: SPICE window containing result
"""
relate = condition.relate
refval = condition.refval
adjust = condition.adjust
nintvls = 2000 # 2*n + ( m / step )
@spice.utils.callbacks.SpiceUDFUNS
def gfq(et):
value = self.compute(et)
return value
@spice.utils.callbacks.SpiceUDFUNB
def gfdecrx(udfuns, et):
return spice.uddc(udfuns, et, 1.0)
@spice.utils.callbacks.SpiceUDFUNS
def dummy_udf(et):
return et
@spice.utils.callbacks.SpiceUDFUNB
def gfq_boolean(udfuns, et):
bool_value = self.compute(et)
if refval: # true
return bool_value
else:
return not (bool_value)
result = stypes.SPICEDOUBLE_CELL(100000)
if self.type == 'boolean':
spice.gfudb(dummy_udf, gfq_boolean, step, cnfine, result)
else:
spice.gfuds(gfq, gfdecrx, relate, refval, adjust, step, nintvls, cnfine, result)
# spice.gfuds(udfuns, udqdec, relate, refval, adjust, step, nintvls, cnfine, result)
return result
[docs]class Vector(Geometry):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.dimension = 1
self.valid = True
[docs] def getData(self, prefix=None):
"""Returns Geometry data as a list or list element (given its prefix)
Args:
prefix: Name of the prefix to get data element from.
"""
prefixes = self.prefix
if prefix in prefixes:
data_element = self.data[prefixes.index(prefix)]
return data_element
else:
return self.data
[docs]class Multidimentional(Geometry):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.dimension = 2
self.valid = True
################################################################################
################################################################################
##
## GEOMETRY CLASSES
##
################################################################################
################################################################################
[docs]class Direction(Multidimentional):
"""Direction abstract class
- Any Direction children classes/objects can be used as input vector for the simulated SC attitude (when not pointing-dependent !!)
- It can be used as input to the Directions_Separation_Angle class.
- It can be represented in the "pointing"/"NOA" coordinates space (PointingViewer)
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.data = [None, None, None]
self.prefix = ['dx', 'dy', 'dz']
self.units = ''
self.n_vectors = 0
self.reference_frame = 'J2000' # default
[docs] def setDirection(self, frame, vectors): # 3-element list (vector), or list of 3-element list (vector)
# set reference frame
self.reference_frame = frame
# check inputs vectors
n_valid_vectors = 0
n_total_vectors = 1
if isinstance(vectors, np.ndarray):
vectors = vectors.tolist()
for vec in vectors:
if isinstance(vec, list):
n_total_vectors = len(vectors)
if len(vec) == 3:
n_valid_vectors += 1
else:
if len(vectors) == 3:
n_valid_vectors = 1
vectors = [vectors]
break
if n_valid_vectors == n_total_vectors:
self.n_vectors = n_valid_vectors
self.data = vectors
else:
self.n_vectors == 0
print('WARNING: Unable to set direction vectors.')
# print(vectors)
[docs] def get_pointing_angles(self, et, observer, target, frame='NOA'):
"""Computes the directions pointing angles in a two-dimensional coordinate system
based on NOA (default) or NPO SPICE dynamic frames.
"""
# Set direction vectors to be transformed to NOA angles
dir_vectors = self.data
pointing_frame = f'{observer}_{target}_{frame}'
rotmat = spice.pxform(self.reference_frame, pointing_frame, et)
dir_noa_angles = []
for dir_vec in dir_vectors:
# transform direction vector into NOA frame
dir_noa_vec = spice.mxv(rotmat, dir_vec)
# derive "pointing" angles from NOA direction vector
yz_plane = spice.nvp2pl([1.0, 0.0, 0.0], [0.0, 0.0, 0.0])
xz_plane = spice.nvp2pl([0.0, 1.0, 0.0], [0.0, 0.0, 0.0])
p_yz = spice.vprjp(dir_noa_vec, yz_plane) # projection of dir_noa_vec onto YZ_NOA plane
p_xz = spice.vprjp(dir_noa_vec, xz_plane) # projection of dir_noa_vec onto YZ_NOA plane
theta_x = np.sign(p_xz[0]) * spice.vsep([0.0, 0.0, 1.0], p_xz) * spice.dpr() # rotation angle between Z_NOA and p_xz, in degrees
theta_y = np.sign(p_yz[1]) * spice.vsep([0.0, 0.0, 1.0], p_yz) * spice.dpr() # rotation angle between Z_NOA and p_xz, in degrees
dir_noa_angles.append([theta_x, -theta_y]) # TODO: "minused" theta_y to fix issue, valid ?
return dir_noa_angles
[docs] def toJ2000(self, et, reference_frame='', vectors=[]):
# Transform a vector defined in given a reference frame into J2000 frame, taking into account for SC and MAJIS
# scanner pointing simulator.
#
# - Vectors defined in SC or SC-fixed frames are transformed using SC pointing simulator (when enabled/on):
#
# SC frame = 'JUICE_SPACECRAFT' -> Simulated_SC_Frame
# SC-fixed frames = [ 'JUICE_MAJIS_BASE', 'JUICE_MAJIS' ]
#
# - Vectors defined in Scan (SC moving structure) or Scan-fixed frames are transformed using Scan pointing
# simulator (when enabled/on):
#
# Scan frame = 'JUICE_MAJIS_SCAN' -> Simulated_Scan_Frame
# Scan-fixed frames = [ 'JUICE_MAJIS_VISNIR', 'JUICE_MAJIS_IR']
#
# - Vectors defined in any other frames are transformed using SPICE Frames sub-system pxform() function.
#
# Set input vectors and reference_frame
external_call = False
if reference_frame and vectors:
external_call = True
frame = reference_frame
vectors = vectors
else:
frame = self.reference_frame
vectors = self.data
# Init output J2000 vectors list
j2000_vectors = []
# Transform input vector in J2000 using SPICE Frames sub-system pxform() function, when:
# - input frame can not be simulated (eg: IAU_CALLISTO), or
# - simulator is disabled.
if frame not in simulator.frames() or simulator.is_off():
frame_j2000_rotmat = spice.pxform(frame, 'J2000', et)
for vec in vectors:
j2000_vectors.append(spice.mxv(frame_j2000_rotmat, vec).tolist())
if not external_call:
self.reference_frame = 'J2000'
self.data = j2000_vectors
return j2000_vectors
# Check that required simulated frame Geometry objects are accessible
if simulator.is_on():
if self.geoevt is None:
print(
'WARNING: Unable to simulate pointing because Direction object is not associated to any GeometryEvent object.')
return None
else:
# (the following could be easily scripted: for frame in simulator.simulated_frames ...)
if simulator.is_on(frame='JUICE_SPACECRAFT'):
if not self.geoevt.has_geometry('Simulated_SC_Frame'):
print(
'WARNING: Unable to simulate SC pointing because Simulated_SC_Frame not found in associated GeometryEvent.')
return None
if simulator.is_on(frame='JUICE_MAJIS_SCAN'):
if not self.geoevt.has_geometry('Simulated_Scan_Frame'):
print(
'WARNING: Unable to simulate Scan pointing because Simulated_Scan_Frame not found in associated GeometryEvent.')
return None
# (the following would be way more difficult to script)
if frame in simulator.frames(frame='JUICE_MAJIS_SCAN'): # input vector defined in a Scan/Scan-fixed frame
if simulator.is_on(frame='JUICE_MAJIS_SCAN'):
# compute transformation matrix from input frame to SC frame, using Scan frame simulator
frame_scan_rotmat = spice.pxform(frame, 'JUICE_MAJIS_SCAN', et) # fixed
#simscan_base_rotmat = self.geoevt.get_geometry_data('Simulated_Scan_Frame')
if et == self.geoevt.reference_time:
# retrieve pre-computed simulated frame at epoch et
simscan_base_rotmat = self.geoevt.get_geometry_data('Simulated_Scan_Frame')
else:
# compute simulated frame at epoch et. eg: Detector_Boresight_Ground_Drift
parameters = self.geoevt.get_geometry_parameters('Simulated_Scan_Frame')
simscan_base_rotmat = Simulated_Scan_Frame(parameters, geoevt=self.geoevt).compute(et)
frame_base_rotmat = spice.mxm(frame_scan_rotmat, simscan_base_rotmat) # hyp: scan = simscan
base_sc_rotmat = spice.pxform('JUICE_MAJIS_BASE', 'JUICE_SPACECRAFT', et) # fixed
frame_sc_rotmat = spice.mxm(frame_base_rotmat, base_sc_rotmat)
else:
# using SPICE pxform() function
frame_sc_rotmat = spice.pxform(frame, 'JUICE_SPACECRAFT', et)
# compute input vectors coordinates in SC frame.
sc_vectors = []
for vec in vectors:
sc_vectors.append(spice.mxv(frame_sc_rotmat, vec).tolist())
else: # input vector defined in a SC/SC-fixed frame
# compute input vector coordinates in SC frame, using SPICE pxform() function.
frame_sc_rotmat = spice.pxform(frame, 'JUICE_SPACECRAFT', et) # fixed
sc_vectors = []
for vec in vectors:
sc_vectors.append(spice.mxv(frame_sc_rotmat, vec).tolist())
# transform SC vector (computed in above step) into J2000
if simulator.is_on(frame='JUICE_SPACECRAFT'):
# using SC frame simulator
if et == self.geoevt.reference_time:
# retrieve pre-computed simulated frame at epoch et
simsc_j2000_rotmat = self.geoevt.get_geometry_data('Simulated_SC_Frame')
else:
# compute simulated frame at epoch et. eg: Detector_Boresight_Ground_Drift
parameters = self.geoevt.get_geometry_parameters('Simulated_SC_Frame')
simsc_j2000_rotmat = Simulated_SC_Frame(parameters, geoevt=self.geoevt).compute(et)
for sc_vec in sc_vectors:
j2000_vectors.append(spice.mxv(simsc_j2000_rotmat, sc_vec).tolist())
else:
# using SPICE pxform() function
sc_j2000_rotmat = spice.pxform('JUICE_SPACECRAFT', 'J2000', et)
for sc_vec in sc_vectors:
j2000_vectors.append(spice.mxv(sc_j2000_rotmat, sc_vec).tolist())
if not external_call:
self.reference_frame = 'J2000'
self.data = j2000_vectors
return j2000_vectors
# TODO: Properly implement a Ring_Location class and subsequent design.
# class Ring_Location(Vector):
# def __init__(self, parameters, geoevt=None, silent=None):
# super().__init__(parameters, geoevt=geoevt, silent=silent)
#
# self.n_points = 0
# self.data = []
# self.target = None
from gfinder.utils import split_polygon
[docs]class Geolocation(Multidimentional):
"""Geolocation abstract class
- MapViewer
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.data = [None, None, None]
self.prefix = ['lon', 'lat', 'alt']
self.units = 'deg'
self.n_points = 0
self.target = None
[docs] def setGeolocation(self, points): # 3-element list (point), or list of 3-element list (points)
# check input points
n_valid_points = 0
n_total_points = 1
if isinstance(points, np.ndarray):
points = points.tolist()
for point in points:
if isinstance(point, list):
n_total_points = len(points)
if len(point) == 3:
n_valid_points += 1
else:
if len(points) == 3:
n_valid_points = 1
points = [points]
break
if n_valid_points == n_total_points:
self.n_points = n_valid_points
else:
self.n_points = 0
print('WARNING: Unable to set geolocation point coordinates.')
# print(points)
points = []
#print('>> setGeolocation: n_points', self.n_points)
# print(self.n_points)
# print(type(points))
# print(points)
# convert from rectangular to planetocentric coordinates
geopoints = []
for point in points:
# compute nearest point on the surface to input observered point, which can be above (alt > 0), on (alt = 0),
# or below the surface (alt < 0)
srfpt, alt = spice.nearpt(point, self.target.radii[0], self.target.radii[1], self.target.radii[2])
if np.isclose(alt, 0.0):
alt = 0.0
# derive planetocentric coordinates.
r, lon, lat = spice.reclat(srfpt)
# add to the list of points
geopoints.append([lon * spice.dpr(), lat * spice.dpr(), alt])
self.data = geopoints
[docs] def getGeoJSON(self, surface_only=False, split=False):
self.n_points = len(self.data) # needed for when loaded opportunity data (only self.data state is preserved)
if self.n_points < 1:
feature = Feature()
elif self.n_points == 1:
point = self.data[0]
feature = Feature(geometry=Point((point[0], point[1], point[2])), properties={})
else:
coordTuples = []
for point in self.data:
if surface_only: # only add points that lie onto the surface
if point[2] <= 0.0:
# coordTuples.append((point[0], point[1], point[2]))
coordTuples.append((point[0], point[1]))
else:
coordTuples.append((point[0], point[1])) #, point[2])) # lon, lat, alt
# TODO: unpatch -> add point[2] above
# ~/workspace/git-repos/gfinder-python-cli/gfinder/utils.py in split_polygon(geojson, output_format, validate)
# 229 crossings = 0
# 230
# --> 231 for coord_index, (lon, _) in enumerate(ring[1:], start=1):
# 232 lon_prev = ring[coord_index - 1][0] # [0] corresponds to longitude coordinate
# 233 if check_crossing(lon, lon_prev, validate=validate):
#
# ValueError: too many values to unpack (expected 2)
if len(coordTuples) > 2:
coordTuples = [coordTuples + [coordTuples[0]]] # complete coordinates tuples to close Polygon
feature = Feature(geometry=Polygon(coordTuples), properties={})
if split:
if len(feature.geometry.coordinates) > 0:
try:
split_polygons = split_polygon(feature.geometry)
except:
print(f'[WARNING] Could not split polygon: {feature.geometry}')
split_polygons = []
if len(split_polygons) > 1: # assuming polygon crosses antimeridian
multi_polygon = []
for polygon in split_polygons:
multi_polygon.append(geojson.loads(polygon).coordinates)
feature = Feature(geometry=MultiPolygon(multi_polygon), properties={})
#feature = tamn.getSplitFeature(feature)
return feature
[docs] def to_rectangular(self):
"""
Returns the rectangular coordinates of the Geolocation points.
"""
recpoints = []
for geopoint in self.data:
recpoints.append(spice.georec(geopoint[0]*spice.rpd(), geopoint[1]*spice.rpd(), geopoint[2], self.target.re, self.target.f))
return recpoints
[docs] def angular_sep(self, geolocation):
"""
Returns ground angular separation between two geolocations (surface points).
Args:
geolocation:
Returns:
"""
# print('> self.n_points = {}, geolocation.n_points = {}'.format(self.n_points, geolocation.n_points))
# print('> len(self.data) = {}, len(geolocation.data) = {}'.format(len(self.data), len(geolocation.data)))
spoint1 = self.to_rectangular()[0] if self.n_points > 0 else [0.0, 0.0, 0.0]
spoint2 = geolocation.to_rectangular()[0] if geolocation.n_points > 0 else [0.0, 0.0, 0.0]
# print('> angular_sep: spoint1={}, spoint2={}'.format(spoint1, spoint2))
angular_sep = spice.vsep(spoint1, spoint2)
return angular_sep
[docs]class Surface_Point_Direction(Direction): # Direction vector from SC to a given target surface point in J2000 frame
"""Direction vector from SC to a given target surface point in J2000 frame.
Args:
required_parameters_sets (list): ['target', 'longitude', 'latitude']
data (list):
Example of ODF definition::
{
"name": "Surface_Point_Direction",
"parameters": {
"target": "JUPITER",
"longitude": 0.0,
"latitude": 0.0
}
}
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.dimension = 1 # force dimension to one; returned data is a list.
self.required_parameters_sets = [['target', 'longitude', 'latitude']]
self.set_parameters(parameters)
#self.reference_frame = Target(self.parameters['target']).frame
[docs] def compute(self, et):
# set default parameters
target = Target(self.parameters['target'])
lons = self.parameters['longitude']
lats = self.parameters['latitude']
if not isinstance(lons, list):
lons = [ lons ]
if not isinstance(lats, list):
lats = [ lats ]
alt = 0.0
spoint_dir_vectors = []
target_pos, ltime = spice.spkpos(target.name, et, target.frame, 'NONE', OBSRVR)
for lon, lat in zip(lons, lats):
lon = lon * spice.rpd()
lat = lat * spice.rpd()
# compute surface point direction from spacecraft in target's body-fixed frame
spoint = spice.georec(lon, lat, alt, target.re, target.f)
spoint_dir_vectors.append((spoint+target_pos).tolist())
self.setDirection(target.frame, spoint_dir_vectors)
self.toJ2000(et)
return self.data
[docs]class North_Pole_Direction(Direction): # SC to north-pole direction in J2000 frame
"""North_Pole_Direction.
Args:
required_parameters_sets (list): ['target']
data (list):
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.dimension = 1 # force dimension to one; returned data is a list.
self.required_parameters_sets = [['target']]
self.set_parameters(parameters)
[docs] def compute(self, et):
target = Target(self.parameters['target'])
longitude = 0.0
latitude = 90.0
self.data = Surface_Point_Direction(
{'target': target.name, 'longitude': longitude, 'latitude': latitude}).compute(et)
return self.data
[docs]class SC_Slew_Angles(Vector):
"""SC_Slew_Angles.
Attributes:
required_parameters_sets (list): [['x_start','x_stop', 'y_start', 'y_stop'], ['x_start', 'x_rate', 'y_start', 'y_rate']]
data (list):
"""
# TODO: Rename to SC_Offset_Angles and return PTR-convention X-offset and Y-offset angles.
def __init__(self, parameters, geoevt=None, silent=None, et_start=None, et_stop=None):
"""Constructs SC_Slew_Angles object.
Args:
parameters:
geoevt:
silent:
et_start (float, optional): start ephemeris time.
et_stop:
"""
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.prefix = ['x_rot', 'y_rot']
self.data = [None, None]
self.units = 'rad'
self.required_parameters_sets = [
['x_start', 'x_stop', 'y_start', 'y_stop'], # x_rate/y_rate relaxed
['x_start', 'x_rate', 'y_start', 'y_rate'], # x_stop/y_stop relaxed
['x_start', 'x_stop', 'x_rate'],
['x_start', 'y_start'], # constant x-y offset
['y_start', 'x_start', 'x_stop', 'x_rate'], # to be able to map with PTR parameters
# ['x_start', 'y_start', 'scan_delta', 'scan_speed'], # PTR approach and convention ...later.
]
self.et_start = None
self.et_stop = None
if et_start is not None:
self.et_start = et_start
if et_stop is not None:
self.et_stop = et_stop
self.set_parameters(parameters)
[docs] def compute(self, et):
# Retrieve start and stop time of the observation associated to this measurement's Geometry object
if self.geoevt:
et_start = self.geoevt.parent_event.get_start_time()
et_stop = self.geoevt.parent_event.get_stop_time()
elif (self.et_start is not None) and (self.et_stop is not None):
et_start = self.et_start
et_stop = self.et_stop
else:
print('Missing associated GeometryEvent object or `et_start` and `et_stop` optional argument.')
return [None, None]
if self.parameters_set_index == 0: # ['x_start', 'x_stop', 'y_start', 'y_stop']
# x_rot
a = (self.parameters['x_stop'] - self.parameters['x_start']) / (et_stop-et_start)
b = self.parameters['x_start'] - a * et_start
x_rot = (a*et + b) * spice.rpd()
# y_rot
a = (self.parameters['y_stop'] - self.parameters['y_start']) / (et_stop-et_start)
b = self.parameters['y_start'] - a * et_start
y_rot = (a*et + b) * spice.rpd()
elif self.parameters_set_index == 1: # ['x_start', 'x_rate', 'y_start', 'y_rate']
# x_rot angle at time et
b = self.parameters['x_start'] - self.parameters['x_rate'] * et_start
x_rot = (self.parameters['x_rate']*et + b) * spice.rpd()
# y_rot angle at time et
b = self.parameters['y_start'] - self.parameters['y_rate'] * et_start
y_rot = (self.parameters['y_rate']*et + b) * spice.rpd()
elif self.parameters_set_index == 2: # ['x_start', 'x_stop', 'x_rate']
# x_rot angle at time et
b = self.parameters['x_start'] - self.parameters['x_rate'] * et_start
x_rot = (self.parameters['x_rate'] * et + b)
if self.parameters['x_stop'] >= 0: # x_rot can't be higher than x_rot(et_stop)
x_rot = min(x_rot, self.parameters['x_stop']) * spice.rpd()
else: # x_rot can't be lower than x_rot(et_stop)
x_rot = max(x_rot, self.parameters['x_stop']) * spice.rpd()
# set fixed no rotation around Y
y_rot = 0.0
elif self.parameters_set_index == 3: # ['x_start', 'y_start']
x_rot = self.parameters['x_start'] * spice.rpd()
y_rot = self.parameters['y_start'] * spice.rpd()
elif self.parameters_set_index == 4: # ['y_start', 'x_start', 'x_stop', 'x_rate']
# x_rot angle at time et
b = self.parameters['x_start'] - self.parameters['x_rate'] * et_start
x_rot = (self.parameters['x_rate'] * et + b)
if self.parameters['x_stop'] >= 0: # x_rot can't be higher than x_rot(et_stop)
x_rot = min(x_rot, self.parameters['x_stop']) * spice.rpd()
else: # x_rot can't be lower than x_rot(et_stop)
x_rot = max(x_rot, self.parameters['x_stop']) * spice.rpd()
# set constant rotation around Y
y_rot = self.parameters['y_start'] * spice.rpd()
else:
return [None, None]
self.data = [x_rot, y_rot]
return self.data
# class MAJIS_Scan_Rate_Angle(Scalar):
# """MAJIS_Scan_Rate_Angle.
#
# Args:
# required_parameters_sets (list): ['x_start','x_rate']
# data (list):
# """
#
# def __init__(self, parameters, geoevt=None, silent=None):
# super().__init__(parameters, geoevt=geoevt, silent=silent)
#
# self.prefix = ['x_rot', 'y_rot']
# self.data = 0.0
# self.units = 'radians'
# self.required_parameters_sets = [['x_start','x_rate']] # rates in iFOV (0.00015 radians) !
#
# self.set_parameters(parameters)
#
# def compute(self, et):
# # set default parameters
# x_start = self.parameters['x_start']
# x_rate = self.parameters['x_rate']
#
# # Retrieve start and stop time of the observation associated to this measurement's Geometry object
# et_start = self.geoevt.parent_event.get_start_time()
# et_stop = self.geoevt.parent_event.get_stop_time()
#
# # x_rot angle at time et
# a = x_rate * IFOV # (x_stop - x_start) / (et_stop - et_start)
# b = x_start*spice.rpd() - a * et_start
# x_rot = (a * et + b) # * spice.rpd()
#
# self.data = x_rot
#
# return self.data
[docs]class MAJIS_Scan_Angle(Scalar):
"""MAJIS_Scan_Angle.
Args:
required_parameters_sets (list): ['start_angle','start_scan_rate', 'stop_scan_rate', 'method']
data (list):
"""
def __init__(self, parameters, geoevt=None, silent=None, scan_duration=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
# self.data = 0.0
self.units = 'rad'
self.required_parameters_sets = [
['start_angle', 'start_scan_rate', 'stop_scan_rate', 'method'],
['start_angle', 'stop_angle', 'start_scan_rate', 'stop_scan_rate', 'method'],
['start_angle', 'scan_rate']
]
# TODO: Below is an example of what would be a much better definition for required parameters sets.
# self.required_parameters_sets = [
# {'start_angle': 0.0, 'start_scan_rate': 0.0, 'stop_scan_rate': 0.0, 'method': ''},
# {'start_angle': 0.0, 'stop_angle': 0.0, 'start_scan_rate': 0.0, 'stop_scan_rate': 0.0, 'method':''},
# {'start_angle': 0.0, 'scan_rate': 0.0}
# ]
self.scan_duration = None
if scan_duration:
self.scan_duration = scan_duration
# Set MAJIS scanner parameters
self.scanner_angular_step = 360 / 8192 / 28
scanner_start_position = 2821 # "min"
scanner_stop_position = 5369 # "max"
n_scanner_total_steps = scanner_stop_position-scanner_start_position # 2548
self.scanner_angular_range = self.scanner_angular_step * n_scanner_total_steps
self.valid_scanner_angles = np.linspace(-self.scanner_angular_range/2, self.scanner_angular_range/2, n_scanner_total_steps+1)
self.valid_scanner_positions = np.arange(scanner_start_position, scanner_stop_position+1)
self.set_parameters(parameters)
if self.valid:
if self.parameters_set_index == 2:
self.method = "discrete" # default
else:
self.set_scan_parameters()
self.method = self.parameters['method'] # 'MC', for testing
[docs] def compute(self, et):
# temporary code for testing purpose only
if self.method == 'MC':
mc_angle = self.geoevt.parent_event.get_sub_events_geometry_data('Motion_Compensation_Angle')[self.geoevt.index]
majis_scan_angle_obj = MAJIS_Scan_Angle({}, silent=True)
valid_scanner_angle, scanner_position = majis_scan_angle_obj.get_valid_scanner_angle(mc_angle*spice.dpr()/2)
self.data = valid_scanner_angle*spice.rpd()*2
return self.data
et_start = 0.0
if self.geoevt:
et_start = self.geoevt.parent_event.get_start_time()
time = et-et_start
if self.parameters_set_index == 2: # ['start_angle', 'scan_rate'] (deg, deg/s)
start_angle = self.parameters['start_angle']
scan_rate = self.parameters['scan_rate']
scan_angle = (scan_rate*time + start_angle) * 2.0 * spice.rpd()
else:
scan_angle = self.get_scanner_angle(time) * 2.0 * spice.rpd()
#print('>>', scan_angle)
self.data = scan_angle # optical, in radians
return self.data
[docs] def set_scan_parameters(self):
""" Set scan parameters from input commanding parameters.
start_angle:
stop_angle:
n_steps:
positions:
times:
t_scan:
t_scan_yl:
"""
start_angle = self.parameters['start_angle']
start_scan_rate = self.parameters['start_scan_rate']
stop_scan_rate = self.parameters['stop_scan_rate']
self.start_angle = start_angle
if 'stop_angle' in self.parameters.keys():
self.stop_angle = self.parameters['stop_angle'] # asymmetrical scan
else:
self.stop_angle = -start_angle # symmetrical scan
# print('>>>', self.stop_angle)
valid_start_angle, self.start_position = self.get_valid_scanner_angle(self.start_angle)
valid_stop_angle, self.stop_position = self.get_valid_scanner_angle(self.stop_angle)
self.n_steps = np.abs(self.stop_position-self.start_position) # always 2548 for a full-range scan
self.positions = np.linspace(self.start_position, self.stop_position, self.n_steps+1, dtype=int)
self.angles = np.linspace(self.start_angle, self.stop_angle, self.n_steps+1)
self.delta_scan_rate = 0.0
if self.n_steps != 0:
self.delta_scan_rate = (stop_scan_rate - start_scan_rate) / self.n_steps
delta_times = []
for n, position in enumerate(self.positions):
scan_rate = start_scan_rate + self.delta_scan_rate * n
delta_time = 0.0
if scan_rate != 0.0:
delta_time = abs(self.scanner_angular_step / scan_rate)
delta_times.append(delta_time)
delta_times = np.array(delta_times)
self.times = np.zeros(self.n_steps+1)
for n, delta_time in enumerate(delta_times):
if n < self.n_steps:
self.times[n+1] = self.times[n] + delta_time
self.t_scan = np.sum(delta_times[:-1])
self.t_scan_yl = 0.0
if stop_scan_rate != start_scan_rate:
if stop_scan_rate != 0.0:
self.t_scan_yl = np.abs(self.scanner_angular_range / (stop_scan_rate-start_scan_rate)*np.log(stop_scan_rate/start_scan_rate))
[docs] def get_scanner_angle(self, time):
""" Returns scanner angle corresponding to a given time since scan start.
Args:
time:
"""
if self.method == 'discrete':
position = round(np.interp(time, self.times, self.positions))
scanner_angle = self.position_to_angle(position)
elif self.method == 'continuous':
scanner_angle = np.interp(time, self.times, self.angles)
else:
print('[WARNING] Unknown MAJIS_Scan_Angle method.')
scanner_angle = 0.0
return scanner_angle
[docs] def position_to_angle(self, position):
""" Returns scanner angle corresponding to a given scanner position.
Args:
position:
"""
scanner_angle = None
indexes = np.where(self.valid_scanner_positions == np.round(position))
if indexes[0].size != 0:
scanner_angle = self.valid_scanner_angles[indexes][0]
return scanner_angle
[docs] def get_valid_scanner_angle(self, angle):
""" Returns valid scanner angle corresponding to an input angle.
Args:
angle: input (mechanical) scanner angle, in degrees.
"""
index = np.abs(self.valid_scanner_angles - angle).argmin()
valid_scanner_angle = self.valid_scanner_angles[index]
scanner_position = self.valid_scanner_positions[index]
return valid_scanner_angle, scanner_position
[docs] def is_valid_scanner_angle(self, angle):
""" Returns whether or not an input angle is valid.
Args:
angle: input (mechanical) scanner angle, in degrees.
Returns:
"""
scanner_position = None
index = np.abs(self.valid_scanner_angles - angle).argmin()
is_valid_scanner_angle = abs(self.valid_scanner_angles[index]-angle) < 1e-12
if is_valid_scanner_angle:
scanner_position = self.valid_scanner_positions[index]
return is_valid_scanner_angle, scanner_position
class SimulatedSCFrameError(Exception):
"""Error raised when computing Simulated_SC_Frame data (rotation matrix) is not possible due to no primary or
secondary vector solution."""
pass
[docs]class Simulated_SC_Frame(Multidimentional):
"""Simulated_SC_Frame.
Args:
required_parameters_sets (list): ['target']
data (list): [[],[],[]]
"""
# Example of output data:
#
# [[-8.78687303e-10 -3.10418898e-05 -3.85203489e-05]
# [-5.93810778e-10 3.62125190e-05 -2.75750771e-05]
# [-3.17808890e-10 1.35058517e-05 -1.45997027e-05]]
#
# Examples of input parameters
#
# {
# "name": "Simulated_SC_Frame",
# "parameters": {
# "primary_axis": "+Z",
# "primary_vector": {
# "name": "Target_Position",
# "parameters": {
# "target": "JUPITER"
# }
# },
# "secondary_axis": "+Y",
# "secondary_vector": {
# "name": "SC_Velocity",
# "parameters": {
# "target": "JUPITER"
# }
# },
# "offset_rotations": "SC_Slew_Angles"
# }
# }
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
# self.prefix = ['x', 'y', 'z']
self.data = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
# self.units = 'km'
self.required_parameters_sets = [
['primary_axis', 'primary_vector', 'secondary_axis', 'secondary_vector', 'offset_rotations'],
['primary_axis', 'primary_vector', 'secondary_axis', 'secondary_vector', 'power_optimised', 'offset_rotations'],
['naif_sc_frame', 'offset_rotations']
]
self.set_parameters(parameters)
[docs] def compute(self, et):
power_optimised = None # default: no power-optimisation
if self.parameters_set_index == 1: # 'power_optimised' is provided
power_optimised = self.parameters['power_optimised']
elif self.parameters_set_index == 2: # 'naif_sc_frame' is provided
# Get input NAIF SC reference frame to J200 transformation matrix
naif_sc_frame = self.parameters['naif_sc_frame']
simsc_j2000_rotmat = spice.pxform(naif_sc_frame, 'J2000', et)
# Apply offset rotation angles
offset_rot_def = self.parameters['offset_rotations']
simsc_j2000_rotmat = self.apply_offset_rotations(simsc_j2000_rotmat, offset_rot_def)
self.data = simsc_j2000_rotmat.tolist()
return self.data
# set default parameters, from parameters sets 0 and 1.
pri_vector_def = self.parameters['primary_vector']
sec_vector_def = self.parameters['secondary_vector']
pri_axis_def = self.parameters['primary_axis']
sec_axis_def = self.parameters['secondary_axis']
offset_rot_def = self.parameters['offset_rotations']
# compute primary vector
pri_vector_data = GeometryFactory().create(pri_vector_def, geoevt=self.geoevt).compute(et)
if pri_vector_data:
pri_vector = pri_vector_data[0]
else:
if simulator.is_on(frame='JUICE_SPACECRAFT'):
raise SimulatedSCFrameError(f'Could not compute {pri_vector_def["name"]} as Simulated_SC_Frame primary vector at ephemeris time {et}. '
f'Try again with a different set of inputs parameters: {pri_vector_def["parameters"]}.')
else:
print(f'[WARNING] Could not compute {pri_vector_def["name"]} as Simulated_SC_Frame primary vector at ephemeris time {et}.')
return np.identity(3).tolist()
# compute secondary vector
sec_vector_data = GeometryFactory().create(sec_vector_def, geoevt=self.geoevt).compute(et)
if sec_vector_data:
sec_vector = sec_vector_data[0]
else:
if simulator.is_on(frame='JUICE_SPACECRAFT'):
raise SimulatedSCFrameError(
f'Could not compute {sec_vector_def["name"]} as Simulated_SC_Frame primary vector at ephemeris time {et}. '
f'Try again with a different set of inputs parameters: {sec_vector_def["parameters"]}.')
else:
print(f'[WARNING] Could not compute {sec_vector_def["name"]} as Simulated_SC_Frame secondary vector at ephemeris time {et}.')
return np.identity(3).tolist()
indexa, signa, framea = self.axis_str(pri_axis_def)
indexp, signp, framep = self.axis_str(sec_axis_def)
axdef = signa * np.array(pri_vector) # vector defining the primary axis
plndef = signp * np.array(sec_vector) # vector defining the secondary axis
# print('>>', axdef, indexa, signa, framea)
# print('>>', plndef, indexp, signp, framep)
# print()
# j2000_simsc_rotmat = spice.twovec(axdef, indexa, plndef, indexp)
# compute J2000 to simulated SC basic pointing frame transformation matrix based on input axis frame definition,
# either 'JUICE_SPACECRAFT' or 'JUICE_MAJIS_BASE'.
if framea == 'JUICE_SPACECRAFT':
j2000_simsc_rotmat = spice.twovec(axdef, indexa, plndef, indexp) # J2000 to simulated SC frame
if framep != framea:
raise NotImplementedError(f'Secondary axis must be defined wrt primary axis frame: {framea}.')
elif framea == 'JUICE_MAJIS_BASE':
j2000_simbase_rotmat = spice.twovec(axdef, indexa, plndef, indexp)
simbase_simsc_rotmat = spice.pxform(framea, 'JUICE_SPACECRAFT', 0.0)
j2000_simsc_rotmat = spice.xpose(spice.mxm(spice.xpose(j2000_simbase_rotmat), spice.xpose(simbase_simsc_rotmat)))
# print(f'{pri_axis_def} = {simbase_simsc_rotmat[:,2]}') # MAJIS +Z axis in SC frame
if framep == 'JUICE_SPACECRAFT': # NPO pointing case
raise NotImplementedError(f'Secondary axis must be defined wrt primary axis frame: {framea}.')
else:
raise Exception(f'Invalid primary axis frame: {framea}.')
# Power-optimised pointing: changing the sign of secondary axis depending on whether or not defined solar-array
# normal vector (-Xsc) is pointing to the Sun for power-optimisation.
#
# TODO: "Power-optimisation" method must be updated to accommodate for new primary and secondary axis definition.
if power_optimised is not None:
# print('>> "power-optimising"')
solar_array_norm_axis = -1 # (-Xsc); defines the orientation of solar panels in SC ref fram; allowed values: -3, -2, -1, 1, 2, 3
solar_array_norm = np.sign(solar_array_norm_axis) * j2000_simsc_rotmat[abs(solar_array_norm_axis) - 1] # in J2000
sunpos, ltime = spice.spkpos('SUN', et, 'J2000', ABCORR, OBSRVR)
solar_array_inc_angle = spice.vsep(solar_array_norm, sunpos)
if power_optimised:
if solar_array_inc_angle > spice.halfpi():
j2000_simsc_rotmat = spice.twovec(axdef, indexa, -plndef, indexp)
else:
if solar_array_inc_angle < spice.halfpi():
j2000_simsc_rotmat = spice.twovec(axdef, indexa, -plndef, indexp)
# Transpose to get simulated SC frame to J2000
simsc_j2000_rotmat = spice.xpose(j2000_simsc_rotmat)
# Apply offset rotation angles
simsc_j2000_rotmat = self.apply_offset_rotations(simsc_j2000_rotmat, offset_rot_def)
self.data = simsc_j2000_rotmat.tolist()
return self.data
[docs] def apply_offset_rotations(self, simsc_j2000_rotmat, offset_rot_def):
"""Returns the SC basic pointing frame rotated around X and Y by angles derived from input "offset_rotation"
geometry definition (eg: 'SC_Slew_Angles', which must be pre-computed.)
"""
# TODO: Implement similar 'offsetRefAxis' (rotation axis) parameter to be aligned with PTR definition/capability.
# Currently, rotations are always applied to input rotation matrix reference frame (basic pointing SC frame).
if offset_rot_def == '':
return simsc_j2000_rotmat
else:
offset_angles = self.geoevt.get_geometry_data(offset_rot_def) # !! again, must placed before Simulated_SC_Frame def in ODF
# parameters = self.geoevt.get_geometry_parameters(offset_rot_def)
# offset_angles = GeometryFactory().create({'name': offset_rot_def, parameters': parameters}).compute(et)
# before: not working because offset rotations geometry definition was not in `measurement_geometries` list.
xrotmat = spice.rotate(offset_angles[0], 1) # rotation around axis X; simsc to x_rot_simsc
simsc_j2000_rotmat = spice.mxm(simsc_j2000_rotmat, spice.xpose(xrotmat)) # x_rot_simsc
yrotmat = spice.rotate(offset_angles[1], 2) # rotation around axis Y; x_rot_simsc to y_rot_simsc
simsc_j2000_rotmat = spice.mxm(simsc_j2000_rotmat, spice.xpose(yrotmat)) # y_rot_simsc
return simsc_j2000_rotmat
[docs] def axis_str(self, axis_def):
"""Parse input axis definition string and returns axis parameters.
Definition string schema: `[+|-][<frame_key>_][X|Y|Z]axis`, where `<frame_key>` can be 'SC' (default) or 'MAJIS'.
"""
sign = 1
if axis_def[0] == '-':
sign = -1
axis_def = axis_def[1:] # remove sign
elif axis_def[0] == '+':
axis_def = axis_def[1:] # remove sign
valid_frames = {
'SC': 'JUICE_SPACECRAFT',
'MAJIS': 'JUICE_MAJIS_BASE'
}
# find axis definition frame and name
def_splits = axis_def.split('_')
if len(def_splits) == 1: # default is SC frame
frame = valid_frames['SC']
name = axis_def
elif len(def_splits) == 2:
if def_splits[0] in valid_frames.keys():
frame = valid_frames[def_splits[0]]
name = def_splits[1]
else:
raise Exception(f'Invalid axis frame definition: {def_splits[0]}.')
else:
raise Exception(f'Invalid axis definition: {axis_def}.')
# find axis index number
axes = {'X': 1, 'Y': 2, 'Z': 3}
if name[0] in axes.keys():
axis_number = axes[name[0]]
else:
raise Exception(f'Invalid axis definition: {axis_def}.')
return axis_number, sign, frame
# class Simulated_SC_Frame(Multidimentional):
# """Simulated_SC_Frame.
#
# Args:
# required_parameters_sets (list): ['target']
# data (list): [[],[],[]]
# """
#
# # Example of output data:
# #
# # [[-8.78687303e-10 -3.10418898e-05 -3.85203489e-05]
# # [-5.93810778e-10 3.62125190e-05 -2.75750771e-05]
# # [-3.17808890e-10 1.35058517e-05 -1.45997027e-05]]
# #
# # Examples of input parameters
# #
# # {
# # "name": "Simulated_SC_Frame",
# # "parameters": {
# # "primary_axis": "+Z",
# # "primary_vector": {
# # "name": "Target_Position",
# # "parameters": {
# # "target": "JUPITER"
# # }
# # },
# # "secondary_axis": "+Y",
# # "secondary_vector": {
# # "name": "SC_Velocity",
# # "parameters": {
# # "target": "JUPITER"
# # }
# # },
# # "offset_rotations": "SC_Slew_Angles"
# # }
# # }
#
# def __init__(self, parameters, geoevt=None, silent=None):
# super().__init__(parameters, geoevt=geoevt, silent=silent)
#
# # self.prefix = ['x', 'y', 'z']
# self.data = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
# # self.units = 'km'
# self.required_parameters_sets = [
# ['primary_axis', 'primary_vector', 'secondary_axis', 'secondary_vector', 'offset_rotations'],
# ['primary_axis', 'primary_vector', 'secondary_axis', 'secondary_vector', 'power_optimised', 'offset_rotations'],
# ['naif_sc_frame', 'offset_rotations']
# ]
#
# self.set_parameters(parameters)
#
# def compute(self, et):
#
# power_optimised = None # default: no power-optimisation
# if self.parameters_set_index == 1: # 'power_optimised' is provided
# power_optimised = self.parameters['power_optimised']
# elif self.parameters_set_index == 2: # 'naif_sc_frame' is provided
# # Get input NAIF SC reference frame to J200 transformation matrix
# naif_sc_frame = self.parameters['naif_sc_frame']
# simsc_j2000_rotmat = spice.pxform(naif_sc_frame, 'J2000', et)
#
# # Apply offset rotation angles
# offset_rot_def = self.parameters['offset_rotations']
# simsc_j2000_rotmat = self.apply_offset_rotations(simsc_j2000_rotmat, offset_rot_def)
#
# self.data = simsc_j2000_rotmat.tolist()
# return self.data
#
# # set default parameters, from parameters sets 0 and 1.
# pri_vector_def = self.parameters['primary_vector']
# sec_vector_def = self.parameters['secondary_vector']
# pri_axis_def = self.parameters['primary_axis']
# sec_axis_def = self.parameters['secondary_axis']
# offset_rot_def = self.parameters['offset_rotations']
#
# # compute primary vector
# pri_vector_data = GeometryFactory().create(pri_vector_def, geoevt=self.geoevt).compute(et)
# if pri_vector_data:
# pri_vector = pri_vector_data[0]
# else:
# if simulator.is_on(frame='JUICE_SPACECRAFT'):
# raise SimulatedSCFrameError(f'Could not compute {pri_vector_def["name"]} as Simulated_SC_Frame primary vector at ephemeris time {et}. '
# f'Try again with a different set of inputs parameters: {pri_vector_def["parameters"]}.')
# else:
# print(f'[WARNING] Could not compute {pri_vector_def["name"]} as Simulated_SC_Frame primary vector at ephemeris time {et}.')
# return np.identity(3).tolist()
#
# # compute secondary vector
# sec_vector_data = GeometryFactory().create(sec_vector_def, geoevt=self.geoevt).compute(et)
# if sec_vector_data:
# sec_vector = sec_vector_data[0]
# else:
# if simulator.is_on(frame='JUICE_SPACECRAFT'):
# raise SimulatedSCFrameError(
# f'Could not compute {sec_vector_def["name"]} as Simulated_SC_Frame primary vector at ephemeris time {et}. '
# f'Try again with a different set of inputs parameters: {sec_vector_def["parameters"]}.')
# else:
# print(f'[WARNING] Could not compute {sec_vector_def["name"]} as Simulated_SC_Frame secondary vector at ephemeris time {et}.')
# return np.identity(3).tolist()
#
# indexa, signa = self.axis_str(pri_axis_def)
# indexp, signp = self.axis_str(sec_axis_def)
#
# axdef = signa * np.array(pri_vector) # vector defining the primary axis
# plndef = signp * np.array(sec_vector) # vector defining the secondary axis
#
# #print('>>', axdef, indexa, signa, plndef, indexp, signp)
# j2000_simsc_rotmat = spice.twovec(axdef, indexa, plndef, indexp) # J2000 to simulated SC frame
#
# # Power-optimised pointing: changing the sign of secondary axis depending on whether or not defined solar-array
# # normal vector (-Xsc) is pointing to the Sun for power-optimisation.
# #
# if power_optimised is not None:
# #print('>> "power-optimising"')
# solar_array_norm_axis = -1 # (-Xsc); defines the orientation of solar panels in SC ref fram; allowed values: -3, -2, -1, 1, 2, 3
# solar_array_norm = np.sign(solar_array_norm_axis) * j2000_simsc_rotmat[abs(solar_array_norm_axis) - 1] # in J2000
# sunpos, ltime = spice.spkpos('SUN', et, 'J2000', ABCORR, OBSRVR)
# solar_array_inc_angle = spice.vsep(solar_array_norm, sunpos)
# if power_optimised:
# if solar_array_inc_angle > spice.halfpi():
# j2000_simsc_rotmat = spice.twovec(axdef, indexa, -plndef, indexp)
# else:
# if solar_array_inc_angle < spice.halfpi():
# j2000_simsc_rotmat = spice.twovec(axdef, indexa, -plndef, indexp)
#
# # Transpose to get simulated SC frame to J2000
# simsc_j2000_rotmat = spice.xpose(j2000_simsc_rotmat)
#
# # Apply offset rotation angles
# simsc_j2000_rotmat = self.apply_offset_rotations(simsc_j2000_rotmat, offset_rot_def)
#
# self.data = simsc_j2000_rotmat.tolist()
# return self.data
#
# def apply_offset_rotations(self, simsc_j2000_rotmat, offset_rot_def):
#
# if offset_rot_def == '':
# return simsc_j2000_rotmat
# else:
# offset_angles = self.geoevt.get_geometry_data(offset_rot_def) # !! again, must placed before Simulated_SC_Frame def in ODF
# # parameters = self.geoevt.get_geometry_parameters(offset_rot_def)
# # offset_angles = GeometryFactory().create({'name': offset_rot_def, parameters': parameters}).compute(et)
# # before: not working because offset rotations geometry definition was not in `measurement_geometries` list.
#
# xrotmat = spice.rotate(offset_angles[0], 1) # rotation around axis X; simsc to x_rot_simsc
# simsc_j2000_rotmat = spice.mxm(simsc_j2000_rotmat, spice.xpose(xrotmat)) # x_rot_simsc
# yrotmat = spice.rotate(offset_angles[1], 2) # rotation around axis Y; x_rot_simsc to y_rot_simsc
# simsc_j2000_rotmat = spice.mxm(simsc_j2000_rotmat, spice.xpose(yrotmat)) # y_rot_simsc
#
# return simsc_j2000_rotmat
#
#
#
# def axis_str(self, axis_def):
# sign = 1
# if axis_def[0] == '-':
# sign = -1
# axes = {'X': 1, '-X': 1, 'Y': 2, '-Y': 2, 'Z': 3, '-Z': 3}
# axis_number = axes[axis_def]
#
# return axis_number, sign
[docs]class Simulated_Scan_Frame(Multidimentional):
"""Simulated_Scan_Frame.
Defined wrt to JUICE_MAJIS_BASE frame.
scan_zero (list): list of Euler rotation angles defining the orientation of the JUICE_MAJIS_SCAN wrt the
JUICE_MAJIS_BASE when the scanner is in "zero" position (scanner rotation axis).
scan_rotation_angle:
Args:
required_parameters_sets (list): ['scan_zero', 'scan_rotation_angle']
data (list): [[],[],[]]
"""
# Example of output data:
#
# [[-8.78687303e-10 -3.10418898e-05 -3.85203489e-05]
# [-5.93810778e-10 3.62125190e-05 -2.75750771e-05]
# [-3.17808890e-10 1.35058517e-05 -1.45997027e-05]]
#
# Examples of input parameters
#
# {
# "name": "Simulated_Scan_Frame",
# "parameters": {
# "scan_zero": [0.0, 0.0, 0.0], # Eg. Euler rotation angles wrt to JUICE_MAJIS_BASE
# "scan_rotation_angle": "MAJIS_Scanner_Rate_Angles" # Function returning the scanner offset rotation angle
# } # at each acquisition time.
#
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
# self.prefix = ['x', 'y', 'z']
self.data = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
# self.units = 'km'
self.required_parameters_sets = [[
'scan_zero',
'scan_rotation_angle'
]]
self.set_parameters(parameters)
[docs] def compute(self, et):
# primary vector is +X that is intended to be co-aligned with +Xmajis_base
# We use Euler angles for now to defined the scan mirror axis
# TODO: take into accounts for maximum MAJIS scan rotation angle limit, using MAJIS_SCAN_ANGLE_MAX
# set default parameters
scan_zero_angles = self.parameters['scan_zero']
scan_rotation_angle_def = self.parameters['scan_rotation_angle']
# Compute rotation matrix from "zero"-position Scan frame to JUICE_MAJIS_BASE.
# Hypothesis on order of rotation ('Z' 'X' 'Y')
# angle3, angle2, angle1, axis3, axis2, axis1
base_simscan_rotmat = spice.eul2m(scan_zero_angles[0], scan_zero_angles[1], scan_zero_angles[2], 3, 1, 2) # JUICE_MAJIS_BASE to "zero" Scan frame
simscan_base_rotmat = spice.xpose(base_simscan_rotmat) # "zero" Scan frame to JUICE_MAJIS_BASE
# print(simscan_base_rotmat)
# Apply mirror rotation angle
scan_rotation_angle_geometry_name = scan_rotation_angle_def # Geometry class name
#scan_rotation_angle = self.geoevt.get_geometry_data(scan_rotation_angle_geometry)
if et == self.geoevt.reference_time:
# retrieve pre-computed scan rotation angle at epoch et
scan_rotation_angle = self.geoevt.get_geometry_data(scan_rotation_angle_geometry_name)
#print(f'line={self.geoevt.index+1}, time={et}, {scan_rotation_angle_geometry_name}={scan_rotation_angle*spice.dpr()}')
else:
# compute scan rotation angle at epoch et. eg: Detector_Boresight_Ground_Drift
if scan_rotation_angle_def == 'Motion_Compensation_Angle':
#print('[WARNING] MAJIS scanner motion can not be simulated at requested epoch (using reference time instead).')
scan_rotation_angle = self.geoevt.get_geometry_data(scan_rotation_angle_geometry_name)
else:
parameters = self.geoevt.get_geometry_parameters(scan_rotation_angle_def)
scan_rotation_angle_geometry_def = {'name': scan_rotation_angle_geometry_name, 'parameters':parameters}
scan_rotation_angle = GeometryFactory().create(scan_rotation_angle_geometry_def, geoevt=self.geoevt).compute(et)
#print(scan_rotation_angle)
# IMPORTANT: again, Geometry definition must placed before Simulated_Scan_Frame def in ODF
xrotmat = spice.rotate(scan_rotation_angle, 1) # apply rotation around axis X; simscan to x_rot_simscan
simscan_base_rotmat = spice.mxm(simscan_base_rotmat, spice.xpose(xrotmat)) # x_rot_simscan
self.data = simscan_base_rotmat.tolist()
return self.data
[docs]class Target_Position(Direction):
"""Target center geometric position vector with respect to spacecraft, expressed in J2000.
Args:
required_parameters_sets (list): ['target']
data (list):
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.dimension = 1 # force dimension to one; returned data is a list.
self.required_parameters_sets = [['target']]
self.set_parameters(parameters)
[docs] def compute(self, et):
ref_frame = 'J2000'
target = Target(self.parameters['target'])
target_state, ltime = spice.spkezr(target.name, et, ref_frame, 'NONE', OBSRVR)
self.setDirection(ref_frame, list(target_state[0:3]))
return self.data
[docs]class SubSC_Point_Direction(Direction):
"""SubSC_Point_Direction.
Args:
required_parameters_sets (list): ['target']
data (list):
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.dimension = 1 # force dimension to one; returned data is a list.
self.required_parameters_sets = [['target']]
self.set_parameters(parameters)
[docs] def compute(self, et):
# set default parameters (could be passed as required parameters if needed)
observer = OBSRVR # JUICE
abcorr = 'NONE' # ABCORR = 'LT+S'
target = Target(self.parameters['target'])
spoint, trgepc, srfvec = spice.subpnt('NEAR POINT/ELLIPSOID', target.name, et, target.frame, abcorr, OBSRVR)
self.setDirection(target.frame, list(srfvec))
self.toJ2000(et)
return self.data
[docs]class SC_Velocity(Direction):
"""Target geometric velocity vector with respect to target, expressed in J2000 frame.
Args:
required_parameters_sets (list): ['target']
data (list): [x, y, z]
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.dimension = 1 # force dimension to one; returned data is a list.
# self.prefix = ['vx', 'vy', 'vz']
# self.data = [0.0, 0.0, 0.0]
# self.units = 'km/s'
self.required_parameters_sets = [['target']]
self.set_parameters(parameters)
[docs] def compute(self, et):
ref_frame = 'J2000'
target = Target(self.parameters['target'])
target_state, ltime = spice.spkezr(OBSRVR, et, ref_frame, 'NONE', target.name)
self.setDirection(ref_frame, list(target_state[3:]))
return self.data
[docs]class Ground_Track_Across_Direction(Direction):
"""NOT AVAILABLE (WIP).
Ground_Track_Across_Direction.
Args:
required_parameters_sets (list): ['target']
data (list):
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.dimension = 1 # force dimension to one; returned data is a list.
self.required_parameters_sets = [['target']]
self.set_parameters(parameters)
[docs] def compute(self, et):
# set default parameters (could be passed as required parameters if needed)
observer = OBSRVR # JUICE
abcorr = 'NONE' # ABCORR = 'LT+S'
target = Target(self.parameters['target'])
spoint1, trgepc1, srfvec = spice.subpnt('NEAR POINT/ELLIPSOID', target.name, et, target.frame, abcorr, OBSRVR)
spoint2, trgepc2, srfvec = spice.subpnt('NEAR POINT/ELLIPSOID', target.name, et+1.0, target.frame, abcorr, OBSRVR)
spoint1 = spice.mxv(spice.pxform(target.frame, 'J2000', et), spoint1)
spoint2 = spice.mxv(spice.pxform(target.frame, 'J2000', et), spoint2)
dir = spice.vcrss(spoint1, spoint2)
self.setDirection('J2000', list(dir))
# self.toJ2000(et)
return self.data
[docs]class Ring_Ansae(Multidimentional):
"""Multidimensional-type data container for pointing-independent Jovian ring's ansae properties.
Geometry classes that users ``Ring_Ansae``:
- :class:`gfinder.geometry.Ring_Illumination_Angle`
- :class:`gfinder.geometry.Ring_Horizontal_Angular_Size`
- :class:`gfinder.geometry.Ring_Vertical_Angular_Size`
- :class:`gfinder.geometry.Ring_West_Ansa_Phase_Angle`
- :class:`gfinder.geometry.Ring_East_Ansa_Phase_Angle`
- :class:`gfinder.geometry.Ring_West_Ansa_Best_Resolution`
- :class:`gfinder.geometry.Ring_East_Ansa_Best_Resolution`
- :class:`gfinder.geometry.Ring_Ansae_Phase_Angle_In_Range`
Args:
required_parameters_sets (list): ['radius']
data (list): [[[sun_elev, h_ang_size, v_ang_size]],[[west_phase, west_resol]],[[east_phase, east_resol]]]
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
# set required parameters sets
self.required_parameters_sets = [['radius']]
self.set_parameters(parameters)
[docs] def compute(self, et):
radius = self.parameters['radius']
target = Target('JUPITER')
# vector from Jupiter to the observer
vec_jup_sc, ltime = spice.spkpos(OBSRVR, et, 'J2000', 'NONE', target.name)
# vector from Sun to Jupiter
vec_sun_jup, ltime = spice.spkpos(target.name, et, 'J2000', 'NONE', 'SUN')
# vector from Jupiter to Sun in Jupiter Inertial Frame ("ji") at J2000
vec_jup_sun = -vec_sun_jup
rotmat = spice.pxform('J2000', 'JUICE_JUPITER_IF_J2000', et)
vec_jup_sun_ji = spice.mxv(rotmat, vec_jup_sun)
ang_jupsun_jupaxis = spice.vsep(vec_jup_sun_ji, np.array([0., 0., 1.]))*spice.dpr()
sun_elev = 90.0-ang_jupsun_jupaxis
# vector from Jupiter to the observer in Jupiter Inertial Frame
vec_jup_sc_ji = spice.mxv(rotmat, vec_jup_sc)
# vector from Sun to Jupiter in Jupiter Inertial Frame
vec_sun_jup_ji = spice.mxv(rotmat, vec_sun_jup)
# vector from jupiter center to ansa, perpendicular to jupiter-sun line
pa = spice.vcrss(vec_sun_jup_ji, np.array([0., 0., 1.]))
pa = spice.vhat(pa)*radius*target.re
vec_sun_ansa1 = vec_sun_jup_ji + pa
vec_sun_ansa2 = vec_sun_jup_ji - pa
vec_sc_ansa1 = pa-vec_jup_sc_ji
vec_sc_ansa2 = -pa-vec_jup_sc_ji
# determine which one is the west or east from the observer point of view
ansae_orient = spice.vdot(spice.vcrss(vec_sc_ansa1, vec_sc_ansa2), np.array([0., 0., 1.]))
if ansae_orient < 0.0:
vec_sc_west_ansa = vec_sc_ansa1
vec_sun_west_ansa = vec_sun_ansa1
vec_sc_east_ansa = vec_sc_ansa2
vec_sun_east_ansa = vec_sun_ansa2
else:
vec_sc_west_ansa = vec_sc_ansa2
vec_sun_west_ansa = vec_sun_ansa2
vec_sc_east_ansa = vec_sc_ansa1
vec_sun_east_ansa = vec_sun_ansa1
# phase angle at each ansa
west_phase = spice.vsep(vec_sun_west_ansa, vec_sc_west_ansa)*spice.dpr()
east_phase = spice.vsep(vec_sun_east_ansa, vec_sc_east_ansa)*spice.dpr()
# best spatial resolution at both ansae
west_resol = spice.vnorm(vec_sc_west_ansa)*IFOV
east_resol = spice.vnorm(vec_sc_east_ansa)*IFOV
# alternative computations of horizontal and vertical angular size of the ring
## find direction vectors tangent to ring ansa.
n_points = 100
azimuths = np.linspace(0.0, 2 * spice.pi(), num=n_points)
xvec = radius*target.re * np.cos(azimuths)
yvec = radius*target.re * np.sin(azimuths)
rpoint_dir_vecs = []
rpoints = []
sep_angles = []
for x, y in zip(xvec, yvec):
rpoint = np.array([x, y, 0.0])
rpoint_dir_vec = rpoint-vec_jup_sc_ji
sep_angle = spice.vsep(rpoint, rpoint_dir_vec)
rpoints.append(rpoint)
rpoint_dir_vecs.append(rpoint_dir_vec)
sep_angles.append(sep_angle*spice.dpr())
tangent_dirs = [] # ring ansa tangent points
for i in range(len(sep_angles)-1):
if sep_angles[i+1] >= 90.0 and sep_angles[i] < 90.0:
r = (90.0-sep_angles[i]) / (sep_angles[i+1] - sep_angles[i])
rpoint0 = r*rpoints[i] + (1-r)*rpoints[i+1]
# print(i, azimuths[i]*spice.dpr(), rpoint_dir_vecs[i])
# print(spice.vsep(rpoint_dir_vecs[i], (rpoint0 - vec_jup_sc_ji)) * spice.dpr())
tangent_dirs.append(rpoint0-vec_jup_sc_ji)
if sep_angles[i+1] < 90.0 and sep_angles[i] >= 90.0:
r = (90.0-sep_angles[i+1]) / (sep_angles[i] - sep_angles[i+1])
rpoint0 = (1-r)*rpoints[i] + r*rpoints[i+1]
# print(i, azimuths[i]*spice.dpr(), rpoint_dir_vecs[i])
# print(spice.vsep(rpoint_dir_vecs[i],(rpoint0-vec_jup_sc_ji))*spice.dpr())
tangent_dirs.append(rpoint0-vec_jup_sc_ji)
h_ang_size0 = None
if len(tangent_dirs) == 2:
# print(tangent_dirs[0])
# print(tangent_dirs[1])
h_ang_size0 = spice.vsep(tangent_dirs[0], tangent_dirs[1])*spice.dpr()
# print(f'Angular size = {h_ang_size} deg')
## vertical angular size
eq_plane = spice.nvc2pl([0.0, 0.0, 1.0], 0.0) # set equatorial plane
pa = spice.vprjp(vec_jup_sc_ji,eq_plane) # projection of sc position in equatorial plane
pa = spice.vhat(pa)*radius*target.re # one of the most vertical ansa
vec_sc_ansa1 = pa-vec_jup_sc_ji
vec_sc_ansa2 = -pa-vec_jup_sc_ji
v_ang_size = spice.vsep(vec_sc_ansa1, vec_sc_ansa2)*spice.dpr()
# print(f'v_vec_sc_ansa1 = {vec_sc_ansa1}')
# print(f'v_vec_sc_ansa2 = {vec_sc_ansa2}')
## horizontal angular size
if h_ang_size0:
h_ang_size = h_ang_size0
else:
pa = spice.vcrss(vec_jup_sc_ji, np.array([0., 0., 1.])) # normal vector to plane defined by sc position and north pole
pa = spice.vhat(pa)*radius*target.re # one of the most horizontal ansa
vec_sc_ansa1 = pa-vec_jup_sc_ji
vec_sc_ansa2 = -pa-vec_jup_sc_ji
h_ang_size = spice.vsep(vec_sc_ansa1, vec_sc_ansa2) * spice.dpr()
# ### radial (horizonal) angular size
# dist = spice.vnorm(vec_jup_sc_ji)
# min_radius = target.re
# max_radius = radius*target.re
# h_ang_size0 = (np.arcsin(max_radius/dist) - np.arcsin(min_radius/dist)) * spice.dpr()
#
# ### ring opening angle (vertical angular size)
# rho = np.sqrt(vec_jup_sc_ji[0]*vec_jup_sc_ji[0] + vec_jup_sc_ji[1]*vec_jup_sc_ji[1])
# v_ang_size0 = abs(np.arctan2(vec_jup_sc_ji[2], rho)) * spice.dpr()
#
# print('Vertical angular size:')
# print(f'- v_ang_size0 = {v_ang_size0:.4} deg')
# print(f'- v_ang_size = {v_ang_size:.4} deg')
# print()
# print('Horizontal angular size:')
# print(f'- h_ang_size0 = {h_ang_size0:.4} deg')
# print(f'- h_ang_size = {h_ang_size:.4} deg')
# print()
ring_values = [sun_elev, h_ang_size, v_ang_size]
west_ansa_values = [west_phase, west_resol]
east_ansa_values = [east_phase, east_resol]
self.data = [ ring_values, west_ansa_values, east_ansa_values ]
return self.data
[docs]class Ring_Illumination_Angle(Scalar):
"""Jovian ring solar illumination angle (deg).
Args:
required_parameters_sets (list): ['radius']
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['radius']]
self.set_parameters(parameters)
[docs] def compute(self, et):
ring_ansae_data = self.geoevt.get_geometry_data('Ring_Ansae')
self.data = ring_ansae_data[0][0]
return self.data
[docs]class Ring_Horizontal_Angular_Size(Scalar):
"""Jovian ring horizontal angular size (deg).
Args:
required_parameters_sets (list): ['radius']
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['radius']]
self.set_parameters(parameters)
[docs] def compute(self, et):
ring_ansae_data = self.geoevt.get_geometry_data('Ring_Ansae')
self.data = ring_ansae_data[0][1]
return self.data
[docs]class Ring_Vertical_Angular_Size(Scalar):
"""Jovian ring vertical angular size (deg).
Args:
required_parameters_sets (list): ['radius']
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['radius']]
self.set_parameters(parameters)
[docs] def compute(self, et):
ring_ansae_data = self.geoevt.get_geometry_data('Ring_Ansae')
self.data = ring_ansae_data[0][2]
return self.data
[docs]class Ring_West_Ansa_Phase_Angle(Scalar):
"""Jovian ring West ansa phase angle (deg).
Args:
required_parameters_sets (list): ['radius']
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['radius']]
self.set_parameters(parameters)
[docs] def compute(self, et):
ring_ansae_data = self.geoevt.get_geometry_data('Ring_Ansae')
self.data = ring_ansae_data[2][0]
return self.data
[docs]class Ring_East_Ansa_Phase_Angle(Scalar):
"""Jovian ring East ansa phase angle (deg).
Args:
required_parameters_sets (list): ['radius']
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['radius']]
self.set_parameters(parameters)
[docs] def compute(self, et):
ring_ansae_data = self.geoevt.get_geometry_data('Ring_Ansae')
self.data = ring_ansae_data[1][0]
return self.data
[docs]class Ring_West_Ansa_Best_Resolution(Scalar):
"""Jovian ring West ansa best resolution (km).
Args:
required_parameters_sets (list): ['radius']
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'km'
self.required_parameters_sets = [['radius']]
self.set_parameters(parameters)
[docs] def compute(self, et):
ring_ansae_data = self.geoevt.get_geometry_data('Ring_Ansae')
self.data = ring_ansae_data[2][1]
return self.data
[docs]class Ring_East_Ansa_Best_Resolution(Scalar):
"""Jovian ring East ansa best resolution (km).
Args:
required_parameters_sets (list): ['radius']
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'km'
self.required_parameters_sets = [['radius']]
self.set_parameters(parameters)
[docs] def compute(self, et):
ring_ansae_data = self.geoevt.get_geometry_data('Ring_Ansae')
self.data = ring_ansae_data[1][1]
return self.data
[docs]class Ring_Ansae_Phase_Angle_In_Range(Scalar):
"""Whether or not at least one Jovian ring ansa phase angle is within a given range.
Args:
required_parameters_sets (list): ['radius', 'min_phase_angle', 'max_phase_angle']
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.type = 'boolean'
self.required_parameters_sets = [['radius', 'min_phase_angle', 'max_phase_angle']]
self.set_parameters(parameters)
[docs] def compute(self, et):
min_phase_angle = self.parameters['min_phase_angle']
max_phase_angle = self.parameters['max_phase_angle']
# wansa_phase_angle = self.geoevt.get_geometry_data('Ring_West_Ansa_Phase_Angle')
# eansa_phase_angle = self.geoevt.get_geometry_data('Ring_East_Ansa_Phase_Angle')
# ring_ansae_data = GeometryFactory().create(pri_vector_def, geoevt=self.geoevt).compute(et)
ring_ansae_data = Ring_Ansae({'radius': self.parameters['radius']}).compute(et)
wansa_phase_angle = ring_ansae_data[1][0]
eansa_phase_angle = ring_ansae_data[2][0]
self.data = False
# method: “at least one ansa test for phase conditions”
if ((wansa_phase_angle > min_phase_angle) and (wansa_phase_angle <= max_phase_angle) or
(eansa_phase_angle > min_phase_angle) and (eansa_phase_angle <= max_phase_angle)):
self.data = True
return self.data
# Coordinates and illumination angles at intersection point between LOS and ring/equatorial plane.
#
[docs]class Observed_Ring_Point(Vector):
"""Vector-type data container for detector's boresight observed Jovian ring's point coordinates and associated
properties.
Geometry classes that users ``Observed_Ring_Point``:
- :class:`gfinder.geometry.Observed_Ring_Point_Radius`
- :class:`gfinder.geometry.Observed_Ring_Point_Azimuth`
- :class:`gfinder.geometry.Observed_Ring_Point_Spatial_Resolution`
- :class:`gfinder.geometry.Observed_Ring_Point_Phase`
- :class:`gfinder.geometry.Observed_Ring_Point_Incidence`
- :class:`gfinder.geometry.Observed_Ring_Point_Emergence`
Args:
required_parameters_sets (list): ['detector']
data (list): [radius, azimuth, slant, phase, incdnc, emissn]
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.prefix = ['radius', 'azimuth', 'resol', 'phase', 'incdnc', 'emissn']
self.data = [None, None, None, None, None, None]
self.units = ['RJ', 'deg', 'km', 'deg', 'deg', 'deg']
self.required_parameters_sets = [['detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
target = Target('JUPITER')
ref_frame = 'JUICE_JUPITER_IF_J2000'
detector = Detector(self.parameters['detector'])
scpos, ltime = spice.spkpos(OBSRVR, et, ref_frame, ABCORR, target.name)
# retrieve and transform detector boresight in JUICE_JUPITER_IF_J2000
if self.geoevt.has_geometry('Detector_Boresight'):
dvec = self.geoevt.get_geometry_data('Detector_Boresight')[0]
else:
dvec = Detector_Boresight(self.parameters).compute(et)[0]
rotmat = spice.pxform('J2000', ref_frame, et)
dvec = spice.mxv(rotmat, dvec)
ring_normal = [0.0, 0.0, 1.0]
plane = spice.nvp2pl(ring_normal, [0.0, 0.0, 0.0])
nxpts, ring_point = spice.inrypl(scpos, dvec, plane)
# https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/inrypl_c.html
if nxpts < 1:
return self.data
radius, azimuth, elevation = spice.reclat(ring_point) # elevation should be equal to zero
if (radius >= RINGS_RADIUS_MIN) and (radius <= RINGS_RADIUS_MAX):
# compute illumination angles at ring-point
sunpos, ltime = spice.spkpos('SUN', et - ltime, ref_frame, ABCORR, target.name)
sundir = sunpos - ring_point
scdir = scpos - ring_point
phase = spice.vsep(sundir, scdir)
incdnc = spice.vsep(sundir, ring_normal)
emissn = spice.vsep(scdir, ring_normal)
resol = spice.vnorm(scdir) * IFOV
self.data = [radius / JUPITER_RADIUS, azimuth * spice.dpr(), resol, phase * spice.dpr(),
incdnc * spice.dpr(), emissn * spice.dpr()]
return self.data
[docs]class Observed_Ring_Point_Radius(Scalar):
"""Detector boresight observed Jovian ring point radius (RJ).
Args:
required_parameters_sets (list): ['detector']
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'RJ'
self.required_parameters_sets = [['detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
data = self.geoevt.get_geometry_data('Observed_Ring_Point')
self.data = data[0]
return self.data
[docs]class Observed_Ring_Point_Azimuth(Scalar):
"""Detector boresight observed Jovian ring point azimuth (deg).
Args:
required_parameters_sets (list): ['detector']
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
data = self.geoevt.get_geometry_data('Observed_Ring_Point')
self.data = data[1]
return self.data
[docs]class Observed_Ring_Point_Spatial_Resolution(Scalar):
"""Detector boresight observed Jovian ring point spatial resolution (km).
Args:
required_parameters_sets (list): ['detector']
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'km'
self.required_parameters_sets = [['detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
# data = Observed_Ring_Point(self.parameters).compute(et)
data = self.geoevt.get_geometry_data('Observed_Ring_Point')
self.data = data[2]
return self.data
[docs]class Observed_Ring_Point_Phase(Scalar):
"""Detector boresight observed Jovian ring point phase angle (deg).
Args:
required_parameters_sets (list): ['detector']
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
# data = Observed_Ring_Point(self.parameters).compute(et)
data = self.geoevt.get_geometry_data('Observed_Ring_Point')
self.data = data[3]
return self.data
[docs]class Observed_Ring_Point_Incidence(Scalar):
"""Detector boresight observed Jovian ring point incidence angle (deg).
Args:
required_parameters_sets (list): ['detector']
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
# data = Observed_Ring_Point(self.parameters).compute(et)
data = self.geoevt.get_geometry_data('Observed_Ring_Point')
self.data = data[4]
return self.data
[docs]class Observed_Ring_Point_Emergence(Scalar):
"""Detector boresight observed Jovian ring point emergence angle (deg).
Args:
required_parameters_sets (list): ['detector']
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
# data = Observed_Ring_Point(self.parameters).compute(et)
data = self.geoevt.get_geometry_data('Observed_Ring_Point')
self.data = data[5]
return self.data
[docs]class Target_Point(Vector):
"""Detector boresight target point coordinates and its associated illumination and viewing angles.
TODO: Shall become deprecated and replaced by Detector_Boresight_Target_Point.
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.prefix = ['lon', 'lat', 'alt', 'slant', 'phase', 'incdnc', 'emissn']
self.data = [None, None, None, None, None, None, None]
self.units = ['deg', 'deg', 'km', 'km', 'deg', 'deg', 'deg']
self.required_parameters_sets = [['target', 'detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
target = Target(self.parameters['target'])
detector = Detector(self.parameters['detector'])
# retrieve detector boresight vector in J2000
dref = 'J2000'
if self.geoevt.has_geometry('Detector_Boresight'):
dvec = self.geoevt.get_geometry_data('Detector_Boresight')[0]
else:
dvec = Detector_Boresight(self.parameters).compute(et)[0]
rotmat = spice.pxform('J2000', target.frame, et)
dvec = spice.mxv(rotmat, dvec)
spoint, trgepc, srfvec, found = spice.sincpt('ELLIPSOID', target.name, et, target.frame, ABCORR, OBSRVR, dref,
dvec)
if found: # observed target point is on the target body surface
r, lon, lat = spice.reclat(spoint)
alt = 0.0
slant = spice.vnorm(srfvec)
else:
targpos, ltime = spice.spkpos(target.name, et, target.frame, ABCORR, OBSRVR)
scpos = -targpos
trg_pnear, dist = spice.npedln(target.radii[0], target.radii[1], target.radii[2], scpos, dvec)
r, lon, lat = spice.reclat(trg_pnear)
spoint = trg_pnear # assign as surface point for illumination angles computation
los_pnear, alt = spice.nplnpt(scpos, dvec, trg_pnear)
slant = spice.vnorm(los_pnear - scpos)
# compute illumination angles at surface point (intercept or nearest)
trgepc, srfvec, phase, incdnc, emissn = spice.ilumin('ELLIPSOID', target.name, et, target.frame, ABCORR, OBSRVR,
spoint)
self.data = [lon * spice.dpr(), lat * spice.dpr(), alt, slant, phase * spice.dpr(), incdnc * spice.dpr(),
emissn * spice.dpr()]
return self.data
[docs]class Target_Point_Longitude(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['target', 'detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
# data = Target_Point(self.parameters).compute(et)
data = self.geoevt.get_geometry_data('Target_Point')
self.data = data[0]
return self.data
[docs]class Target_Point_Latitude(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['target', 'detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
data = self.geoevt.get_geometry_data('Target_Point')
self.data = data[1]
return self.data
[docs]class Target_Point_Altitude(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'km'
self.required_parameters_sets = [['target', 'detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
data = self.geoevt.get_geometry_data('Target_Point')
self.data = data[2]
return self.data
[docs]class Target_Point_Slant(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'km'
self.required_parameters_sets = [['target', 'detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
data = self.geoevt.get_geometry_data('Target_Point')
self.data = data[3]
return self.data
[docs]class Target_Point_Phase(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['target', 'detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
data = self.geoevt.get_geometry_data('Target_Point')
self.data = data[4]
return self.data
[docs]class Target_Point_Incidence(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['target', 'detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
data = self.geoevt.get_geometry_data('Target_Point')
self.data = data[5]
return self.data
[docs]class Target_Point_Emission(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['target', 'detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
data = self.geoevt.get_geometry_data('Target_Point')
self.data = data[6]
return self.data
[docs]class SC_Elevation(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['target', 'reference_frame']]
self.set_parameters(parameters)
[docs] def compute(self, et):
target = Target(self.parameters['target'])
ref_frame = self.parameters['reference_frame']
sc_pos = SC_Position(self.parameters).compute(et)
# sc_pos, ltime = spice.spkpos('JUICE', et, ref_frame, 'NONE', target.name)
r, lon, lat = spice.reclat(sc_pos)
self.data = lat * spice.dpr()
return self.data
[docs]class SC_Position(Vector):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.prefix = ['x', 'y', 'z']
self.data = [None, None, None]
self.units = 'km'
self.required_parameters_sets = [['target', 'reference_frame']]
self.set_parameters(parameters)
[docs] def compute(self, et):
target = Target(self.parameters['target'])
ref_frame = self.parameters['reference_frame']
sc_pos, ltime = spice.spkpos('JUICE', et, ref_frame, 'NONE', target.name)
self.data = [sc_pos[0], sc_pos[1], sc_pos[2]]
return self.data
[docs]class Ring_Ansae_Directions(Direction):
"""Direction vector(s) in J2000 from the observer to a given Jovian ring point or set of points.
TODO: Rename Ring_Ansae_Directions into Ring_Directions and create a new Ring_Point_Direction class.
Args:
required_parameters_sets (list): [['radius', 'azimuth'], ['radius', 'n_points']]
data (list):
Example of ODF definition::
{
"name": "Ring_Ansae_Directions",
"parameters": {
"radius": "@ring_radius",
"n_points": 100
}
}
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.required_parameters_sets = [['radius', 'azimuth'], ['radius', 'n_points']]
self.set_parameters(parameters)
# self.reference_frame = Target(self.parameters['target']).frame
[docs] def compute(self, et):
# set default parameters
target = Target('JUPITER')
if self.parameters_set_index == 0:
# TODO: when single 'radius', 'azimuth' ring point.
lons = self.parameters['radius']
lats = self.parameters['latitude']
if not isinstance(lons, list):
lons = [lons]
if not isinstance(lats, list):
lats = [lats]
alt = 0.0
spoint_dir_vectors = []
target_pos, ltime = spice.spkpos(target.name, et, target.frame, ABCORR, OBSRVR)
for lon, lat in zip(lons, lats):
lon = lon * spice.rpd()
lat = lat * spice.rpd()
# compute surface point direction from spacecraft in target's body-fixed frame
spoint = spice.georec(lon, lat, alt, target.re, target.f)
spoint_dir_vectors.append((spoint + target_pos).tolist())
self.setDirection(target.frame, spoint_dir_vectors)
self.toJ2000(et)
elif self.parameters_set_index == 1:
# set radius, n_points
radius = self.parameters['radius']*JUPITER_RADIUS
n_points = self.parameters['n_points']
# calculate latitudinal coordinates of ring points
azimuths = np.linspace(0.0,2*spice.pi(), num=n_points)
xvec = radius*np.cos(azimuths)
yvec = radius*np.sin(azimuths)
rpoint_dir_vectors = []
target_pos, ltime = spice.spkpos(target.name, et, target.frame, ABCORR, OBSRVR)
for x, y in zip(xvec, yvec):
rpoint= np.array([x, y, 0.0])
rpoint_dir_vectors.append((rpoint + target_pos).tolist())
self.setDirection(target.frame, rpoint_dir_vectors)
self.toJ2000(et)
return self.data
# Return the terminator points, times associated with terminator points, and terminator vectors emanating from the observer.
#
[docs]class Terminator(Multidimentional): # for the sake of a better Geometry parent class for now.
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
# self.prefix = []
self.required_parameters_sets = [['target']]
self.set_parameters(parameters)
[docs] def compute(self, et):
# set required parameter
target = Target(self.parameters['target'])
# set spice.limbpt input parameters
method = 'UMBRAL/TANGENT/ELLIPSOID'
corloc = 'CENTER'
refvec = [0., 0., 1.0]
ncuts = 100
rolstp = spice.twopi() / ncuts
schstp = 0
soltol = 0
maxn = 10000
npts, points, epochs, tangts = spice.termpt(method, 'SUN', target.name, et, target.frame, ABCORR, corloc,
OBSRVR, refvec, rolstp, ncuts, schstp, soltol, maxn);
self.data = [points, epochs, tangts]
return self.data
# Return the limb points of tangency on the target of rays emanating from the observer, in target body-fixed reference frame.
#
[docs]class Limb(Multidimentional): # for the sake of a better Geometry parent class for now.
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
# self.prefix = []
self.required_parameters_sets = [['target']]
self.set_parameters(parameters)
[docs] def compute(self, et):
# set required parameter
target = Target(self.parameters['target'])
# set spice.limbpt input parameters
method = 'TANGENT/ELLIPSOID'
corloc = 'CENTER'
refvec = [0., 0., 1.0]
ncuts = 100
rolstp = spice.twopi() / ncuts
schstp = 0
soltol = 0
maxn = 10000
npts, points, epochs, tangts = spice.limbpt(method, target.name, et, target.frame, ABCORR, corloc, OBSRVR,
refvec, rolstp, ncuts, schstp, soltol, maxn)
self.data = [points, epochs, tangts]
return self.data
# TODO: A generic Geometry.get[Sibling]GeometryData(name) should be implemented to easily retrieve
# "sibling" Geometry objects, introducing the concept of GeometryEvent "buffer" Geometry objects
# for pre-computed data. Eg. for Limb: self.get[Sibling]GeometryData('Limb') would return sibling 'Limb' geometry
# if computed (in ODF, placed before 'Limb_Disk' definition; or in buffer), or compute it (once) and
# store it in the buffer geometries list.
[docs]class Limb_Directions(Direction):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.required_parameters_sets = [['target']]
self.set_parameters(parameters)
[docs] def compute(self, et):
# set required parameter
target = Target(self.parameters['target'])
# compute limb tangent points in body-fixed frame (directions from observer)
limb_data = Limb(self.parameters).compute(et)
self.setDirection(target.frame, limb_data[2])
self.toJ2000(et)
return self.data
[docs]class Limb_Point_Direction(Direction):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.dimension = 1 # force dimension to one; returned data is a list.
self.required_parameters_sets = [['target', 'latitude', 'max_incidence_angle']]
self.set_parameters(parameters)
[docs] def compute(self, et):
target = Target(self.parameters['target'])
lat0 = self.parameters['latitude']
max_incidence_angle = self.parameters['max_incidence_angle']
limb_data = Limb({'target':target.name}).compute(et) # [points, epochs, tangts]
points = limb_data[0]
epochs = limb_data[1]
tangts = limb_data[2]
# Find desired limb point, eg: could be the brightest point on the diurnal limb,
# or a specific latitude if visible. Implementation below is the latter.
limb_lats = []
limb_points = []
limb_tangts = []
limb_epochs = []
# max_incidence_angle = 90.0
for point, tangt, epoch in zip(points, tangts, epochs):
r, lon, lat = spice.reclat(point)
trgepc, srfvec, phase, incdnc, emissn = spice.ilumin('ELLIPSOID', target.name, et, target.frame, ABCORR,
OBSRVR, point)
if incdnc * spice.dpr() < max_incidence_angle:
limb_lats.append(lat * spice.dpr())
limb_points.append(point)
limb_tangts.append(tangt)
limb_epochs.append(epoch)
n_points = len(limb_lats)
if n_points < 2:
return None
## convert to np arrays
limb_lats = np.array(limb_lats)
limb_points = np.array(limb_points)
limb_tangts = np.array(limb_tangts)
limb_epochs = np.array(limb_epochs)
## sort arrays based on latitudes array sorted in ascending order
sort_idxs = limb_lats.argsort()
limb_lats = limb_lats[sort_idxs]
limb_points = limb_points[sort_idxs]
limb_tangts = limb_tangts[sort_idxs]
limb_epochs = limb_epochs[sort_idxs]
# search and interpolate
found = False
for i in range(n_points - 1):
if limb_lats[i] <= lat0 and limb_lats[i + 1] > lat0:
r = (lat0 - limb_lats[i]) / (limb_lats[i + 1] - limb_lats[i])
limb_tangt0 = (1-r) * limb_tangts[i] + r*limb_tangts[i+1]
epoch0 = (1-r) * limb_epochs[i] + r*limb_epochs[i+1]
found = True
#limb_r0, limb_lon0, limb_lat0 = spice.reclat(limb_point0)
#print('r = {}, bracket latitudes:({},{}), limb_lat0={}'.format(r,limb_lats[i], limb_lats[i+1], limb_lat0*spice.dpr()))
if not found:
return None
# set and transform direction vectors in J2000 frame
self.setDirection(target.frame, limb_tangt0)
self.toJ2000(et)
return self.data
[docs]class Limb_Min_Latitude(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['target', 'max_incidence_angle']]
self.set_parameters(parameters)
[docs] def compute(self, et):
self.data = None
limb_lats = Limb_Latitudes(self.parameters).compute(et)
if limb_lats:
self.data = min(limb_lats)
return self.data
[docs]class Limb_Max_Latitude(Scalar): # DOING...
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['target', 'max_incidence_angle']]
self.set_parameters(parameters)
[docs] def compute(self, et):
self.data = None
limb_lats = Limb_Latitudes(self.parameters).compute(et)
if limb_lats:
self.data = max(limb_lats)
return self.data
[docs]class Limb_Latitudes(Vector):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.data = [None, None, None, None]
self.units = 'deg'
self.required_parameters_sets = [['target', 'max_incidence_angle']]
self.set_parameters(parameters)
[docs] def compute(self, et):
target = Target(self.parameters['target'])
max_incidence_angle = self.parameters['max_incidence_angle']
# set spice.limbpt input parameters
method = 'TANGENT/ELLIPSOID'
corloc = 'CENTER'
refvec = [0., 0., 1.0]
ncuts = 100
rolstp = spice.twopi() / ncuts
schstp = 0
soltol = 0
maxn = 10000
npts, points, epochs, tangts = spice.limbpt(method, target.name, et, target.frame, ABCORR, corloc, OBSRVR,
refvec, rolstp, ncuts, schstp, soltol, maxn)
limb_lats = []
for point in points:
r, lon, lat = spice.reclat(point)
trgepc, srfvec, phase, incdnc, emissn = spice.ilumin('ELLIPSOID', target.name, et, target.frame, ABCORR,
OBSRVR, point)
if incdnc * spice.dpr() < max_incidence_angle: # diurnal limb
limb_lats.append(lat * spice.dpr())
self.data = limb_lats
return self.data
[docs]class Limb_Point_Tangent_Direction(Direction):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.dimension = 1 # force dimension to one; returned data is a list.
self.required_parameters_sets = [['target', 'latitude', 'max_incidence_angle', 'power_optimised']]
self.set_parameters(parameters)
[docs] def compute(self, et):
target = Target(self.parameters['target'])
lat0 = self.parameters['latitude']
max_incidence_angle = self.parameters['max_incidence_angle']
power_optimised = self.parameters['power_optimised']
limb_data = Limb({'target': target.name}).compute(et) # [points, epochs, tangts]
points = limb_data[0]
epochs = limb_data[1]
tangts = limb_data[2]
# Find desired limb point, eg: could be the brightest point on the diurnal limb, or a specific latitude if visible. Herebelow the latter.
limb_lats = []
limb_points = []
limb_tangts = []
limb_epochs = []
# max_incidence_angle = 90.0
for point, tangt, epoch in zip(points, tangts, epochs):
r, lon, lat = spice.reclat(point)
trgepc, srfvec, phase, incdnc, emissn = spice.ilumin('ELLIPSOID', target.name, et, target.frame, ABCORR,
OBSRVR, point)
if incdnc * spice.dpr() < max_incidence_angle:
limb_lats.append(lat * spice.dpr())
limb_points.append(point)
limb_tangts.append(tangt)
limb_epochs.append(epoch)
n_points = len(limb_lats)
if n_points < 2:
return None
## convert to np arrays
limb_lats = np.array(limb_lats)
limb_points = np.array(limb_points)
limb_tangts = np.array(limb_tangts)
limb_epochs = np.array(limb_epochs)
## sort arrays based on latitudes array sorted in ascending order
sort_idxs = limb_lats.argsort()
limb_lats = limb_lats[sort_idxs]
limb_points = limb_points[sort_idxs]
limb_tangts = limb_tangts[sort_idxs]
limb_epochs = limb_epochs[sort_idxs]
# search and interpolate
found = False
for i in range(n_points - 1):
if limb_lats[i] <= lat0 and limb_lats[i + 1] > lat0:
r = (lat0 - limb_lats[i]) / (limb_lats[i + 1] - limb_lats[i])
limb_point0 = (1-r) * limb_points[i] + r * limb_points[i+1]
limb_tangt0 = (1-r) * limb_tangts[i] + r * limb_tangts[i+1]
epoch0 = (1-r) * limb_epochs[i] + r * limb_epochs[i+1]
found = True
if not found:
return None
tangent_dir = spice.vcrss(limb_point0, limb_tangt0)
sunpos, ltime = spice.spkpos('SUN', et, target.frame, ABCORR, OBSRVR)
if power_optimised: # Sun direction angle must be minimal
if spice.vsep(tangent_dir, sunpos) > spice.halfpi():
tangent_dir = -tangent_dir
else: # Sun direction angle must be maximal
if spice.vsep(tangent_dir, sunpos) < spice.halfpi():
tangent_dir = -tangent_dir
# set and transform direction vectors in J2000 frame
self.setDirection(target.frame, tangent_dir)
self.toJ2000(et)
return self.data
[docs]class SC_In_Equatorial_Plane(Scalar):
"""SC_In_Equatorial_Plane.
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.type = 'boolean'
self.required_parameters_sets = [['target', 'margin']]
self.set_parameters(parameters)
[docs] def compute(self, et):
target = Target(self.parameters['target'])
margin = self.parameters['margin']
sc_pos, ltime = spice.spkpos('JUICE', et, target.frame, 'NONE', target.name)
r, lon, lat = spice.reclat(sc_pos)
lat = lat * spice.dpr()
sc_in_equatorial_plane = False
if (lat <= margin) and (lat >= -margin):
sc_in_equatorial_plane = True
self.data = sc_in_equatorial_plane
return self.data
from trianglesolver import solve
[docs]class Equatorial_Coverage(Scalar):
"""Geometry class representing the maximum MAJIS scanner latitudinal coverage on a given target, as seen from its equatorial plane.
Initially implemented to enable the definition of MAJIS Jupiter Disk Nadir Scan observation opportunity.
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.required_parameters_sets = [['target']]
self.set_parameters(parameters)
self.units = 'deg'
[docs] def compute(self, et):
target = Target(self.parameters['target'])
# retrieve target body mean radius
radius = np.max(target.radii)
# Compute distance from target centre to spacecraft
targpos, ltime = spice.spkpos(target.name, et, target.frame, ABCORR, OBSRVR)
distance = spice.vnorm(targpos)
# Compute latitudinal coverage A, using trianglesolver.
try:
a, b, c, A, B, C = solve(b=distance, c=radius, C=MAJIS_SCAN_ANGLE_MAX, ssa_flag='obtuse')
except:
if distance <= radius:
A = 0
else:
A = 0.5 * spice.pi() - np.arcsin(radius / distance)
self.data = A * spice.dpr()
return self.data
[docs]class ROI_In_FOV(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.type = 'boolean'
self.required_parameters_sets = [['roi', 'detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
roi = ROI(self.parameters['roi'])
detector = Detector(self.parameters['detector'])
# Determine ray pointing from ROI center location on the target body surface to the spacecraft/detector
roi_spoint = roi.getSurfacePoint()
target_pos, ltime = spice.spkpos(roi.target.name, et, roi.target.frame, ABCORR, OBSRVR)
raydir = -(target_pos + roi_spoint)
# Determine of ROI is occulted by the target (found = True)
spoint, trgepc, srfvec, found = spice.sincpt('ELLIPSOID', roi.target.name, et, roi.target.frame, ABCORR, OBSRVR,
roi.target.frame, raydir)
# Determine if a ray is within the FOV of a specified detector
if found:
self.data = False
else:
self.data = spice.fovray(detector.name, -raydir, roi.target.frame, 'NONE', OBSRVR, et)
return self.data
[docs]class ROI_Boresight_Angle(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['roi', 'detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
roi = ROI(self.parameters['roi'])
detector = Detector(self.parameters['detector'])
# Determine ROI position/direction in detector frame at epoch et
roi_spoint = roi.getSurfacePoint()
target_pos, ltime = spice.spkpos(roi.target.name, et, roi.target.frame, ABCORR, OBSRVR)
roi_dir = target_pos + roi_spoint
rotmat = spice.pxform(roi.target.frame, detector.frame, et) # from body-fixed to detector frame at epoch et
roi_dir = spice.mxv(rotmat, roi_dir) # in detector frame
# Determine separation angle between ROI position vector and detector boresight
roi_boresight_angle = spice.vsep([0., 0., 1.], roi_dir) * spice.dpr()
self.data = roi_boresight_angle
return self.data
from gfinder.roi import ROI
[docs]class ROI_X_Depointing_Angle(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['target', 'detector', 'roi']]
self.set_parameters(parameters)
[docs] def compute(self, et):
roi_def = {
"id": "NA",
"target": self.parameters['target'],
"center_lon": self.parameters['roi']['parameters']['lon'],
"center_lat": self.parameters['roi']['parameters']['lat'],
"diameter": 1.0
}
roi = ROI(roi_def)
detector = Detector(self.parameters['detector'])
# Determine ROI position in detector frame at epoch et
roi_spoint = roi.getSurfacePoint()
target_pos, ltime = spice.spkpos(roi.target.name, et, roi.target.frame, ABCORR, OBSRVR)
roi_dir = target_pos + roi_spoint
rotmat = spice.pxform(roi.target.frame, detector.frame, et) # from body-fixed to detector frame at epoch et
roi_dir = spice.mxv(rotmat, roi_dir) # in detector frame
# Determine separation angle between ROI position vector and XZ plane of detector frame
plane = spice.nvp2pl([0., 1., 0.], [0., 0., 0.])
roi_x_depointing_angle = spice.vsep(spice.vprjp(roi_dir, plane), roi_dir) * spice.dpr()
self.data = roi_x_depointing_angle
return self.data
[docs]class ROI_Y_Depointing_Angle(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['roi', 'detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
roi = ROI(self.parameters['roi'])
detector = Detector(self.parameters['detector'])
# Determine ROI position in detector frame at epoch et
roi_spoint = roi.getSurfacePoint()
target_pos, ltime = spice.spkpos(roi.target.name, et, roi.target.frame, ABCORR, OBSRVR)
roi_dir = target_pos + roi_spoint
rotmat = spice.pxform(roi.target.frame, detector.frame, et) # from body-fixed to detector frame at epoch et
roi_dir = spice.mxv(rotmat, roi_dir) # in detector frame
# Determine separation angle between ROI position vector and YZ plane of detector frame
plane = spice.nvp2pl([1., 0., 0.], [0., 0., 0.])
roi_y_depointing_angle = spice.vsep(spice.vprjp(roi_dir, plane), roi_dir) * spice.dpr()
self.data = roi_y_depointing_angle
return self.data
[docs]class Angular_Separation(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['target1', 'target2']]
self.set_parameters(parameters)
[docs] def compute(self, et):
target1 = Target(self.parameters['target1'])
target2 = Target(self.parameters['target2'])
ref_frame = 'J2000'
# retrieve position vector of target1 as seen from observer in ref frame
target1_pos, ltime = spice.spkpos(target1.name, et, ref_frame, ABCORR, OBSRVR)
# retrieve position vector of target2 as seen from observer in ref frame
target2_pos, ltime = spice.spkpos(target2.name, et, ref_frame, ABCORR, OBSRVR)
# compute separation with taking account for target angular size
angular_sep = spice.vsep(target1_pos, target2_pos) * spice.dpr()
target1_ang_size = Target_Angular_Size({'target': target1.name}).compute(et)
target2_ang_size = Target_Angular_Size({'target': target2.name}).compute(et)
angular_sep = angular_sep - (target1_ang_size + target2_ang_size) / 2
self.data = angular_sep * spice.dpr()
return self.data
[docs]class Target_Solar_Eclipse(Scalar):
"""Target_Solar_Eclipse.
https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/occult_c.html
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.required_parameters_sets = [['target', 'occult_type']]
self.set_parameters(parameters)
self.type = 'boolean'
[docs] def compute(self, et):
target = Target(self.parameters['target'])
occult_type = self.parameters['occult_type']
occult_code = spice.occult('SUN', 'ELLIPSOID', 'IAU_SUN', target.name, 'ELLIPSOID', target.frame, ABCORR, OBSRVR, et)
occult_types = {
'TOTAL': -3,
'ANNULAR': -2,
'PARTIAL': -1
}
self.data = False
if occult_code < 0:
self.data = (occult_code == occult_types[occult_type])
return self.data
[docs]class No_Compensation_Required(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.required_parameters_sets = [['margin']]
self.set_parameters(parameters)
self.type = 'boolean'
[docs] def compute(self, et):
margin = self.parameters['margin']
# Retrieve dwell time from sibling Dwell_Time geometry data.
# or compute using Dwell_Time class (TODO)
dwell_time = self.geoevt.parent_event.get_sub_events_geometry_data('Dwell_Time')[self.geoevt.index]
time_step = self.geoevt.parent_event.time_step # observation measurements time step
# !!Geometries depending of parent GeometryEvent object (Geometry objects container) are NOT SEARCHABLE !!
# dwell_time = time_step = 1.0
# Compensation is required when dwell time is significantly different from the repetition time step.
# If time_step - margin (%) < Dwell Time < time_step + margin (%), then no compensation is required.
#
no_compensation_required = False
if (dwell_time <= time_step * (1 + margin)) and dwell_time >= time_step * (1 - margin):
no_compensation_required = True
self.data = no_compensation_required
return self.data
[docs]class Motion_Compensation_Angle(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'rad'
self.required_parameters_sets = [['target', 'detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
if et == self.geoevt.reference_time:
# TODO: check that parameters match Motion_Compensation parameters.
mc_data = self.geoevt.parent_event.get_sub_events_geometry_data('Motion_Compensation')[self.geoevt.index]
self.data = mc_data[0]
else:
self.data = None
# if self.data > MAJIS_SCAN_ANGLE_MAX*spice.rpd() or self.data < -MAJIS_SCAN_ANGLE_MAX*spice.rpd():
# print('Invalid scanner angle: line number = {}, angle = {}'.format(self.geoevt.index, self.data*spice.dpr()/2))
return self.data
[docs]class Motion_Compensation(Multidimentional):
"""Geometry class that computes the detector boresight rotation angle required for two consecutive detector FOV
(slit) footprints to be contiguous (avoiding overlapping or drifting). Motion compensation is required when dwell
time is greater or lower than the repetition time between two detector acquisitions.
It holds the coordinates of the rotated (compensated) forward- and backward-boresights surface intersect
points in the target body-fixed frame. These coordinates are used to derive the required rotation angle of the
next detector acquisition, and scan mirror position.
In our computation approach, we currently make the following hypothesis:
- Boresight points in the "nadir" direction (vector from spacecraft to sub-spacecraft point or target center)
- Slit (detector FOV) is perpendicular to orbital plane (velocity vector)
- Spacecraft is in forward motion (velocity vector co-aligned with +Ysc axis)
Note for later: Direction-class-based method, by which only the SC pointing would be simulated (disabling the MAJIS
Scan frame simulator). This would allow to determine the scanner motion required by any SC pointings (off-track
pointing for example).
simulator.turn_off(frame='JUICE_MAJIS_SCAN')
detector_bsight = Direction()
detector_bsight.setDirection(self.detector.frame, self.detector.get_forward_bsight())
#detector_bsight.toJ2000(et) Q: why would I transform to J2000 ?
forward_bsight_pvec = detector_bsight.data
simulator.turn_on(frame='JUICE_MAJIS_SCAN')
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.required_parameters_sets = [['target', 'detector', 'start_scan_mode']]
self.set_parameters(parameters)
[docs] def compute(self, et):
# Set parameters
self.target = Target(self.parameters['target'])
self.detector = Detector(self.parameters['detector'])
self.start_scan_mode = self.parameters['start_scan_mode']
# set start scanner angle depending of input start_scan_mode: "zero" or "max_range" at the moment.
start_angle = 0.0
if self.start_scan_mode == 'max_range':
start_angle = MAJIS_Scan_Angle({}, silent=True).scanner_angular_range * spice.rpd()
# set number of detector binning (used to derive forward- and backward-boresight directions)
binning = float(self.detector.binning)
# print(f'> binning = {binning}')
# Compute detector boresight / nadir direction, spacecraft position and scan mirror rotation axis in body-fixed frame
spoint, trgepc, srfvec = spice.subpnt('INTERCEPT/ELLIPSOID', self.target.name, et, self.target.frame, ABCORR, OBSRVR)
bsight = spice.vhat(srfvec) # ~ NADIR. unit vector from observer to sub-observer point at epoch et, in body-fixed frame.
target_state, ltime = spice.spkezr(self.target.name, et, self.target.frame, 'NONE', OBSRVR)
sc_position = -target_state[0:3]
sc_velocity = -target_state[3:]
rot_axis = spice.vcrss(sc_velocity, bsight) # hypothesis for now: scan mirror rotation axis.
# Adjust start angle and rot_axis depending on whether +Ysc (simulated or not) points forward or backward w.r.t
# spacecraft motion.
if SC_Y_Axis_Velocity_Angle({'target': self.target.name}, geoevt=self.geoevt).compute(et) < 90.0:
# forward motion
start_angle = -start_angle
# print(f'> SC forward motion')
else:
# backward motion
rot_axis = -rot_axis
# print(f'> SC backward motion')
# print(f'> sc_position = {sc_position}')
# print(f'> rot_axis = {rot_axis}')
# print(f'> start_angle = {start_angle}')
# Compute rotation matrix from target body-fixed to J2000 frame
rot_mat = spice.pxform(self.target.frame, 'J2000', et-ltime)
# Compute the first acquisition forward-boresight surface intersect point
#
if self.geoevt.index == 0:
# Note applicable anymore:
# # Adjust start angle depending of whether we are in the "too-fast" (dwell time < Trep) or "too-slow"
# # (dwell time < Trep) scenario.
# if self.start_scan_mode == 'max_range': # otherwise not needed because start scanner angle is zero.
# t_rep = self.geoevt.parent_event.time_step
# dwell_time = Dwell_Time({'target':self.target.name, 'detector': self.detector.name}).compute(et)
# if dwell_time > t_rep:
# start_angle = -start_angle
# Set initial required scan angle to start angle
req_scanner_angle = start_angle
# Simulate initial rotated boresight, taking into account for start scanner angle
# bsight = spice.vrotv(bsight, rot_axis, start_angle)
forward_bsight_spoint = [0.0, 0.0, 0.0]
backward_bsight_spoint = [0.0, 0.0, 0.0]
# Compute forward-boresight pointing vector at epoch et.
# forward_bsight_pvec = spice.vrotv(bsight, rot_axis, -0.5*IFOV*binning) # body-fixed
forward_bsight_pvec = spice.vrotv(bsight, rot_axis, -0.5*IFOV*binning + req_scanner_angle) # body-fixed
forward_bsight_pvec = spice.mxv(rot_mat, forward_bsight_pvec) # J2000
# compute forward-boresight surface intersect point
spoint, trgepc, srfvec, found = spice.sincpt('ELLIPSOID', self.target.name, et, self.target.frame, ABCORR,
OBSRVR, 'J2000', forward_bsight_pvec)
if found:
forward_bsight_spoint = spoint
else:
print('[WARNING] Motion_Compensation #1 : No surface intercept point found.')
# DEBUG
# # compute backward-boresight point vector at epoch et, in body-fixed frame.
# backward_bsight_pvec = spice.vrotv(bsight, rot_axis, 0.5 * IFOV * binning)
# backward_bsight_pvec = spice.mxv(rot_mat, backward_bsight_pvec) # J2000
# print(f'> start_angle = {start_angle}')
# print(f'> backward_bsight_pvec = {backward_bsight_pvec}')
# print(f'> req_scanner_angle = {req_scanner_angle}')
else:
# Retrieve previous forward-boresight surface intersect point
prev_motion_comp_data = self.geoevt.parent_event.get_sub_events_geometry_data('Motion_Compensation')[self.geoevt.index - 1]
prev_forward_bsight_spoint = prev_motion_comp_data[2]
# Compute the required backward-boresight pointing vector, at epoch et, to intersect with previous
# forward-boresight pointing vector surface intersection point (leading to acquiring contiguous iFOV
# footprints); defined by the difference between the latter and the detector/SC position at epoch et.
#
sp_ltime = ltime * spice.vnorm(prev_forward_bsight_spoint - sc_position) / spice.vnorm(sc_position)
target_pos, ltime = spice.spkpos(self.target.name, et-sp_ltime, self.target.frame, 'NONE', OBSRVR)
#rot_mat = spice.pxfrm2(self.target.frame, 'J2000', et-sp_ltime, et) ...not needed at the end.
rot_mat = spice.pxform(self.target.frame, 'J2000', et)
req_backward_bsight_pvec = spice.vhat(prev_forward_bsight_spoint + target_pos) # body-fixed
req_backward_bsight_pvec = spice.mxv(rot_mat, req_backward_bsight_pvec) # J2000
# Compute required scanner mirror rotation at epoch et, defined as the angular separation between the
# required backward-boresight pointing vector and the backward-boresight point vector.
#
# compute backward-boresight point vector at epoch et, in body-fixed frame.
# backward_bsight_pvec = spice.vrotv(bsight, rot_axis, 0.5*IFOV*binning)
backward_bsight_pvec = spice.vrotv(bsight, rot_axis, 0.5 * IFOV * binning)
backward_bsight_pvec = spice.mxv(rot_mat, backward_bsight_pvec) # J2000
# print(f'> backward_bsight_pvec = {backward_bsight_pvec}')
# print(f'> req_backward_bsight_pvec = {req_backward_bsight_pvec}')
req_scanner_angle = spice.vsep(backward_bsight_pvec, req_backward_bsight_pvec)
# Get oriented angle (invert angle if backward to required-backward pointing vector rotation axis is
# opposite to scanner rotation axis).
# print(f'> vdot_rotaxis = {spice.vdot(spice.vcrss(backward_bsight_pvec, req_backward_bsight_pvec), rot_axis)}')
if spice.vdot(spice.vcrss(backward_bsight_pvec, req_backward_bsight_pvec), rot_axis) < 0:
req_scanner_angle = -req_scanner_angle
# print('> inverted')
# print(f'> req_scanner_angle = {req_scanner_angle}')
# Compute corresponding (required) forward-boresight surface intercept point, from forward-boresight
# pointing vector and required scanner rotation angle.
# req_forward_bsight_pvec = spice.vrotv(bsight, rot_axis, -0.5*IFOV*binning + req_scanner_angle)
req_forward_bsight_pvec = spice.vrotv(bsight, rot_axis, -0.5 * IFOV * binning + req_scanner_angle)
req_forward_bsight_pvec = spice.mxv(rot_mat, req_forward_bsight_pvec) # J2000
# compute forward-boresight surface intersect point
spoint, trgepc, srfvec, found = spice.sincpt('ELLIPSOID', self.target.name, et, self.target.frame, ABCORR,
OBSRVR, 'J2000', req_forward_bsight_pvec)
if found:
forward_bsight_spoint = spoint
else:
forward_bsight_spoint = [0.0, 0.0, 0.0]
print('[WARNING] Motion_Compensation #2 : No surface intercept point found.')
backward_bsight_spoint = prev_forward_bsight_spoint
self.data = [req_scanner_angle, list(backward_bsight_spoint), list(forward_bsight_spoint)]
# print('>')
return self.data
[docs]class Target_Angular_Size(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['target']]
self.set_parameters(parameters)
[docs] def compute(self, et):
target = self.parameters['target']
fixref = 'IAU_' + target
# retrieve target body mean radius
dim, radii = spice.bodvrd(target, 'RADII', 3)
radius = np.max(radii)
# Compute distance from target centre to spacecraft
targpos, ltime = spice.spkpos(target, et, fixref, ABCORR, OBSRVR)
distance = spice.vnorm(targpos)
self.data = 2.0 * np.arcsin(radius / distance) * spice.dpr()
return self.data
[docs]class SC_Pointing_Offsets(Vector):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.prefix = ['x_sc_angle','y_sc_angle', 'z_sc_angle']
self.data = [None]*3
self.units = 'deg'
self.required_parameters_sets = [['reference_frame']]
self.set_parameters(parameters)
[docs] def compute(self, et):
reference_frame = self.parameters['reference_frame']
# Compute SC +Z axis/vector in reference frame
rotmat = spice.pxform('JUICE_SPACECRAFT', reference_frame, et)
# could the reference frame be the Simulated_SC_Frame ?
# It problably should!
# z_axis = spice.mxv(rotmat, [0., 0., 1.])
# self.data = [z_axis[0], z_axis[1]]
rotmat = spice.pxform('JUICE_SPACECRAFT', reference_frame, et)
x_sc = spice.vsep([1., 0., 0.], spice.mxv(rotmat, [1., 0., 0.])) * spice.dpr()
y_sc = spice.vsep([0., 1., 0.], spice.mxv(rotmat, [0., 1., 0.])) * spice.dpr()
z_sc = spice.vsep([0., 0., 1.], spice.mxv(rotmat, [0., 0., 1.])) * spice.dpr()
self.data = [x_sc, y_sc, z_sc]
return self.data
# Direction vector in J2000 specified by a vector in the spacecraft reference frame (JUICE_SPACECRAFT)
#
[docs]class SC_Vector(Direction):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.dimension = 1 # force dimension to one; returned data is a list.
self.required_parameters_sets = [['vector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
vector = self.parameters['vector']
self.setDirection('JUICE_SPACECRAFT', vector)
self.toJ2000(et)
return self.data
# Boresight direction vector in J2000 for a given input detector name.
#
[docs]class Detector_Boresight(Direction):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.dimension = 1 # force dimension to one; returned data is a list.
self.required_parameters_sets = [['detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
# set required parameters
detector = Detector(self.parameters['detector'])
# set and transform direction vector in J2000 frame
self.setDirection(detector.frame, detector.get_bsight())
self.toJ2000(et)
return self.data
[docs]class Detector_Forward_Boresight(Direction):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.dimension = 1 # force dimension to one; returned data is a list.
self.required_parameters_sets = [['detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
# set required parameters
detector = Detector(self.parameters['detector'])
# set and transform direction vector in J2000 frame
self.setDirection(detector.frame, detector.get_forward_bsight())
self.toJ2000(et)
return self.data
[docs]class Detector_Backward_Boresight(Direction):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.dimension = 1 # force dimension to one; returned data is a list.
self.required_parameters_sets = [['detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
# set required parameters
detector = Detector(self.parameters['detector'])
# set and transform direction vector in J2000 frame
self.setDirection(detector.frame, detector.get_backward_bsight())
self.toJ2000(et)
return self.data
# FOV boundary direction vectors in J2000 for a given input detector name.
#
[docs]class Detector_FOV(Direction):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.required_parameters_sets = [['detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
# set required parameters
detector = Detector(self.parameters['detector'])
# set and transform direction vectors in J2000 frame
self.setDirection(detector.frame, detector.get_bounds())
self.toJ2000(et)
# print('Detector_FOV.compute(et) direction vectors:')
# print(self.data)
return self.data
[docs]class Detector_Boresight_Altitude(Scalar):
"""Altitude of the detector boresight tangent point above the target surface."""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'km'
self.required_parameters_sets = [['target', 'detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
# set required parameters
#self.detector = Detector(self.parameters['detector'])
bsight_target_point_geometry = \
self.geoevt.parent_event.get_sub_events_geometries('Detector_Boresight_Target_Point')[self.geoevt.index]
# print(bsight_target_point_geometry.data)
self.data = bsight_target_point_geometry.data[0][2] # altitude
# print(self.data)
# print()
return self.data
[docs]class Detector_Boresight_Ground_Drift(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'IFOV'
self.required_parameters_sets = [['target', 'detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
# set required parameters
self.detector = Detector(self.parameters['detector'])
# # Compute detector boresight ground angular drift between two computation time step
# time_step = self.get_time_step()
# geoloc1 = Detector_Boresight_Target_Point(self.parameters, geoevt=self.geoevt)
# geoloc1.compute(et)
# geoloc2 = Detector_Boresight_Target_Point(self.parameters, geoevt=self.geoevt)
# geoloc2.compute(et+time_step)
# bsight_drift = geoloc1.angular_sep(geoloc2)
# # print('> geoloc1.data={}, geoloc2.data={}'.format(geoloc1.data, geoloc2.data))
if self.geoevt.index == 0:
self.data = None
return self.data
else:
# Retrieve current and previous Detector_Boresight_Target_Point.
# Rule: Detector_Boresight_Target_Point must be placed before Detector_Boresight_Ground_Drift in ODF file.
bsight_target_point_geometry = \
self.geoevt.parent_event.get_sub_events_geometries('Detector_Boresight_Target_Point')[self.geoevt.index]
prev_bsight_target_point_geometry = \
self.geoevt.parent_event.get_sub_events_geometries('Detector_Boresight_Target_Point')[self.geoevt.index - 1]
# Compute angular separation between both geolocations
bsight_drift = bsight_target_point_geometry.angular_sep(prev_bsight_target_point_geometry)
# Compute
dir_def = {'name': 'Detector_Forward_Boresight', 'parameters': {'detector': self.parameters['detector'] }}
forward_bsight_geoloc = Observed_Target_Points({'target': self.parameters['target'], 'direction': dir_def},
geoevt=self.geoevt)
forward_bsight_geoloc.compute(et)
dir_def = {'name': 'Detector_Backward_Boresight', 'parameters': {'detector': self.parameters['detector']}}
backward_bsight_geoloc = Observed_Target_Points({'target': self.parameters['target'], 'direction': dir_def},
geoevt=self.geoevt)
backward_bsight_geoloc.compute(et)
ifov_ground_res = forward_bsight_geoloc.angular_sep(backward_bsight_geoloc)
#print('> time_step={}, bsight_drift={}, ifov_ground_res={}, ratio={}'.format(time_step, bsight_drift, ifov_ground_res, bsight_drift / ifov_ground_res))
if ifov_ground_res != 0.0:
self.data = bsight_drift / ifov_ground_res
else:
self.data = 0.0
#print('This?! > ', bsight_drift, ifov_ground_res)
return self.data
[docs]class Directions_Angular_Separation(Scalar):
"""Compute separation angle between two direction vectors.
Example of usage in ODF file::
"name": "Directions_Angular_Separation",
"parameters": {
"direction1": {
"name": "SC_Vector",
"parameters": {
"vector": [0.0, 1.0, 0.0]
}
},
"direction2": {
"name": "Observer_Target_Velocity",
"parameters": {
"target": "JUPITER"
}
},
}
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['direction1', 'direction2']]
self.set_parameters(parameters)
[docs] def compute(self, et):
dir_vec1 = GeometryFactory().create(self.parameters['direction1'], geoevt=self.geoevt).compute(et)[0]
dir_vec2 = GeometryFactory().create(self.parameters['direction2'], geoevt=self.geoevt).compute(et)[0]
self.data = spice.vsep(dir_vec1, dir_vec2) * spice.dpr()
return self.data
[docs]class SC_X_Axis_Sun_Angle(Scalar):
# Useful to check solar array and radiator illumination.
# Ideally, one could use it to through warnings, eg: "Solar Array not in light!", "Radiator not in shadow!"
#
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
#self.required_parameters_sets = [['target']]
self.units = 'deg'
self.set_parameters(parameters)
[docs] def compute(self, et):
dir1_def = {'name': 'SC_Vector', 'parameters': {'vector': [1.0, 0.0, 0.0]}}
dir2_def = {'name': 'Target_Position', 'parameters': {'target': 'SUN'}}
self.data = Directions_Angular_Separation({'direction1': dir1_def, 'direction2': dir2_def},
geoevt=self.geoevt).compute(et)
return self.data
[docs]class SC_Solar_Array_Sun_Angle(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.set_parameters(parameters)
[docs] def compute(self, et):
dir1_def = {'name': 'SC_Vector', 'parameters': {'vector': [-1.0, 0.0, 0.0]}}
dir2_def = {'name': 'Target_Position', 'parameters': {'target': 'SUN'}}
self.data = Directions_Angular_Separation({'direction1': dir1_def, 'direction2': dir2_def},
geoevt=self.geoevt).compute(et)
return self.data
[docs]class SC_Y_Axis_Velocity_Angle(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['target']]
self.set_parameters(parameters)
[docs] def compute(self, et):
dir1_def = {'name': 'SC_Vector', 'parameters': {'vector': [0.0, 1.0, 0.0]}}
dir2_def = {'name': 'SC_Velocity', 'parameters': {'target': self.parameters['target']}}
self.data = Directions_Angular_Separation({'direction1': dir1_def, 'direction2': dir2_def},
geoevt=self.geoevt).compute(et)
return self.data
[docs]class SC_Z_Axis_Sun_Angle(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.set_parameters(parameters)
[docs] def compute(self, et):
dir1_def = {'name': 'SC_Vector', 'parameters': {'vector': [0.0, 0.0, 1.0]}}
dir2_def = {'name': 'Target_Position', 'parameters': {'target': 'SUN'}}
self.data = Directions_Angular_Separation({'direction1': dir1_def, 'direction2': dir2_def},
geoevt=self.geoevt).compute(et)
return self.data
[docs]class SC_Z_Axis_Target_Angle(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['target']]
self.set_parameters(parameters)
[docs] def compute(self, et):
dir1_def = {'name': 'SC_Vector', 'parameters': {'vector': [0.0, 0.0, 1.0]}}
dir2_def = {'name': 'Target_Position', 'parameters': {'target': self.parameters['target']}}
self.data = Directions_Angular_Separation({'direction1': dir1_def, 'direction2': dir2_def},
geoevt=self.geoevt).compute(et)
return self.data
[docs]class Boresight_Target_Angle(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['detector', 'target']]
self.set_parameters(parameters)
[docs] def compute(self, et):
dir1_def = {'name': 'Detector_Boresight', 'parameters': {'detector': self.parameters['detector']}}
dir2_def = {'name': 'Target_Position', 'parameters': {'target': self.parameters['target']}}
self.data = Directions_Angular_Separation({'direction1': dir1_def, 'direction2': dir2_def},
geoevt=self.geoevt).compute(et)
return self.data
[docs]class Observed_Target_Points(Geolocation):
# Compute the geodetic coordinates of the point(s) observed by one or several direction vectors, from observer (SC).
# Used by Detector_FOV_Footprint, Detector_Boresight_Target_Point.
#
# "name": "Observed_Target_Points",
# "parameters": {
# "target": "CALLISTO",
# "direction": {
# "name": "Detector_Boresight",
# "parameters": {
# "detector": "JUICE_MAJIS_VISNIR"
# }
# },
# }
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.required_parameters_sets = [['target', 'direction']]
self.set_parameters(parameters)
self.data = []
[docs] def compute(self, et):
# Set/initialise required parameters
self.target = Target(self.parameters['target'])
obs_points = [] # geodetic coordinates of observed target points wrt the target's surface [lon, lat, altitude]
# set parameters
dref = 'J2000'
method = 'ELLIPSOID' # should later become a parameter taking into account for digital shape model (if required)
# Retrieve or compute direction vectors in J2000 frame
if self.geoevt:
#print('>>> retrieving <{}> direction vectors.'.format(self.parameters['direction']['name']))
dir_vectors = self.geoevt.get_geometry_data(self.parameters['direction']['name'])
if dir_vectors is None:
#print('>>>>> Direction object not pre-computed (because not in ODF)')
#print('>>>>> computing <{}> direction vectors.'.format(self.parameters['direction']['name']))
dir_vectors = GeometryFactory().create(self.parameters['direction'], geoevt=self.geoevt).compute(et)
else:
# print('>>> computing <{}> direction vectors.'.format(self.parameters['direction']['name']))
dir_vectors = GeometryFactory().create(self.parameters['direction'], geoevt=self.geoevt).compute(et)
# compute observed target point for each direction vector (whether surface intersect or nearest point)
for dir_vector in dir_vectors:
dvec = np.array(dir_vector)
spoint, trgepc, srfvec, found = spice.sincpt(method, self.target.name, et, self.target.frame, ABCORR, OBSRVR, dref, dvec)
obspt = spoint
if not found:
# compute tangent point, using the new SPICE `tangpt` routine, CSPICE Version 1.0.0, 28-OCT-2021 (MCS)
tanpt, alt1, range, srfpt, trgepc, srfvec = spice.tangpt(method, self.target.name, et, self.target.frame, ABCORR, 'TANGENT POINT', OBSRVR, dref, dvec)
obspt = tanpt
obs_points.append(obspt.tolist())
# Set geolocation points from observation points
self.setGeolocation(obs_points)
return self.data
[docs]class Detector_Boresight_Target_Point(Geolocation):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.dimension = 1 # force dimension to one; returned data is a list.
self.required_parameters_sets = [['target', 'detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
# Set/initialise required parameters
self.target = Target(self.parameters['target'])
dir_def = {'name': 'Detector_Boresight', 'parameters': {'detector': self.parameters['detector'] }}
self.data = Observed_Target_Points({'target': self.parameters['target'], 'direction': dir_def}, geoevt=self.geoevt).compute(et)
self.n_points = len(self.data) # update Geolocation number of points
return self.data
[docs]class Vertical_Resolution(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'km'
self.required_parameters_sets = [['target', 'detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
# compute distance to limb
limb_distance = Limb_Distance({'target': self.parameters['target']}).compute(et)
# compute corresponding spatial resolution
self.data = limb_distance * IFOV * Detector(self.parameters['detector']).binning # .ifov attribute could be used instead
return self.data
[docs]class Limb_Distance(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'km'
self.required_parameters_sets = [['target']]
self.set_parameters(parameters)
[docs] def compute(self, et):
self.target = Target(self.parameters['target'])
# compute spacecraft altitude
sc_altitude = SC_Altitude({'target': self.target.name}).compute(et)
# retrieve target radius
radius = self.target.re
# compute to distance to limb: radius^2 + limb_distance^2 = (radius+ sc_altitude)^2
self.data = np.sqrt(np.square(radius+sc_altitude)-np.square(radius))
return self.data
[docs]class Spatial_Resolution(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'km'
self.required_parameters_sets = [['target', 'detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
# compute spacecraft altitude
sc_altitude = SC_Altitude({'target': self.parameters['target']}).compute(et)
# compute corresponding spatial resolution
self.data = sc_altitude * IFOV * Detector(self.parameters['detector']).binning # .ifov attribute could be used instead
return self.data
[docs]class Ground_Speed(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'km/s'
self.required_parameters_sets = [['target']]
self.set_parameters(parameters)
[docs] def compute(self, et):
self.target = Target(self.parameters['target'])
# Simple method
spoint1, trgepc, srfvec = spice.subpnt('NEAR POINT/ELLIPSOID', self.target.name, et, self.target.frame, ABCORR, OBSRVR)
spoint2, trgepc, srfvec = spice.subpnt('NEAR POINT/ELLIPSOID', self.target.name, et + 1.0, self.target.frame, ABCORR, OBSRVR)
self.data = spice.vnorm(spoint1 - spoint2)
return self.data
[docs]class Dwell_Time(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 's'
self.required_parameters_sets = [['target', 'detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
ifov_ground_resolution = Spatial_Resolution(self.parameters).compute(et)
ground_speed = Ground_Speed({'target': self.parameters['target']}).compute(et)
self.data = (ifov_ground_resolution / ground_speed) #/ self.get_time_step()
return self.data
[docs]class SC_Altitude(Scalar): # apparent sub-observer point altitude
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'km'
self.required_parameters_sets = [['target']]
self.set_parameters(parameters)
[docs] def compute(self, et):
self.target = Target(self.parameters['target'])
spoint, trgepc, srfvec = spice.subpnt('NEAR POINT/ELLIPSOID', self.target.name, et, self.target.frame, ABCORR, OBSRVR)
self.data = spice.vnorm(srfvec)
return self.data
[docs]class Obs_Min_SC_Altitude(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'km'
self.set_parameters(parameters)
[docs] def compute(self, et):
if self.geoevt is not None:
sc_altitudes = self.geoevt.get_sub_events_geometry_data('SC_Altitude')
self.data = min(sc_altitudes)
else:
self.data = None
self.valid = False
print('Missing GeometryEvent object containing this geometry.')
return self.data
[docs]class Obs_Max_SC_Altitude(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'km'
self.set_parameters(parameters)
[docs] def compute(self, et):
if self.geoevt is not None:
sc_altitudes = self.geoevt.get_sub_events_geometry_data('SC_Altitude')
self.data = max(sc_altitudes)
else:
self.data = None
self.valid = False
print('Missing GeometryEvent object containing this geometry.')
return self.data
[docs]class Distance_To_Target(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'km'
self.required_parameters_sets = [['target']]
self.set_parameters(parameters)
[docs] def compute(self, et):
target = self.parameters['target']
fixref = 'IAU_' + target
# Compute vector from target centre to spacecraft
targpos, ltime = spice.spkpos(target, et, fixref, ABCORR, OBSRVR)
self.data = spice.vnorm(targpos)
return self.data
[docs]class Obs_Min_Spatial_Resolution(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'km'
self.set_parameters(parameters)
[docs] def compute(self, et):
if self.geoevt is not None:
spatial_resolutions = self.geoevt.get_sub_events_geometry_data('Spatial_Resolution')
self.data = min(spatial_resolutions)
else:
self.data = None
self.valid = False
print('Missing GeometryEvent object containing this geometry.')
return self.data
[docs]class Obs_Max_Spatial_Resolution(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'km'
self.set_parameters(parameters)
[docs] def compute(self, et):
if self.geoevt is not None:
spatial_resolutions = self.geoevt.get_sub_events_geometry_data('Spatial_Resolution')
self.data = max(spatial_resolutions)
else:
self.data = None
self.valid = False
print('Missing GeometryEvent object containing this geometry.')
return self.data
[docs]class Phase_Angle(Scalar): # parameters = {'target':'GANYMEDE', 'detector': 'JUICE_MAJIS_VISNIR'}
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.required_parameters_sets = [['target']]
self.set_parameters(parameters)
[docs] def compute(self, et):
target = self.parameters['target']
fixref = 'IAU_' + target
spoint, trgepc, srfvec = spice.subpnt('NEAR POINT/ELLIPSOID', target, et, fixref, ABCORR, OBSRVR)
trgepc, srfvec, phase, incdnc, emissn = spice.ilumin('ELLIPSOID', target, et, fixref, ABCORR, OBSRVR, spoint)
self.data = phase * spice.dpr()
return self.data
[docs]class SubSC_Point(Geolocation):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.dimension = 1
self.required_parameters_sets = [['target']]
self.set_parameters(parameters)
# def __format__(self, geometry_fmt):
# if self.data is not None:
# str = '[{:f},{:f}]'.format(self.data[0],self.data[1])
# else:
# str = '[{:f},{:f}]'.format(0,0)
# return str
[docs] def compute(self, et):
# set target
self.target = Target(self.parameters['target'])
# compute sub-spacecraft point rectangular coordinates
spoint, trgepc, srfvec = spice.subpnt('NEAR POINT/ELLIPSOID', self.target.name, et, self.target.frame, ABCORR, OBSRVR)
self.setGeolocation(spoint)
return self.data
[docs]class Obs_SC_Slew_Parameters(Vector):
"""Obs_SC_Slew_Parameters.
Requires Target_Angular_Size data.
"""
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.data = [None] * 5
self.prefix = ['x_start', 'x_stop', 'x_rate', 'slew_duration', 'n_lines']
self.units = ['deg', 'deg', 'deg/s', 's', '']
self.required_parameters_sets = [['detector']]
self.set_parameters(parameters)
[docs] def compute(self, et):
target_angular_sizes = self.geoevt.get_measurements_geometry_data('Target_Angular_Size')
# x_start = -0.5*target_angular_sizes[0]
# x_stop = 0.5*target_angular_sizes[-1]
# angular_range = x_stop-x_start
# detector_ifov = Detector(self.parameters['detector']).ifov[1] * spice.dpr()
# slew_duration = self.geoevt.time_step * (angular_range/detector_ifov)
# x_rate = angular_range/slew_duration
# n_lines = slew_duration/self.geoevt.time_step
x_start = -0.5 * target_angular_sizes[0]
x_stop = 0.5 *target_angular_sizes[-1]
angular_range = x_stop - x_start
detector_ifov = Detector(self.parameters['detector']).ifov[1] * spice.dpr()
time_step = self.geoevt.time_step
x_rate = detector_ifov / time_step
n_lines = int(angular_range / detector_ifov)
slew_duration = n_lines * time_step
self.data = [x_start, x_stop, x_rate, slew_duration, n_lines]
return self.data
from numpy.polynomial import Polynomial
[docs]class Obs_Scanner_Motion_Parameters(Vector):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.data = [None] * 5
self.prefix = ['start_angle', 'start_scan_rate', 'stop_scan_rate', 'scan_duration', 'scan_lines']
self.units = ['deg', 'deg/s', 'deg/s', 's', '']
self.set_parameters(parameters)
[docs] def compute(self, et):
if self.geoevt is not None:
# Retrieve pre-computed scanner motion compensation "optimised" optical angles and acquisition times
mc_angles = np.array(self.geoevt.get_measurements_geometry_data('Motion_Compensation_Angle')) * spice.dpr() / 2
etimes = np.array(self.geoevt.get_measurements_times(format='et'))
times = etimes - etimes[0]
# Determine if and when MC angle goes beyond scanner angular range (+/- 2.0 degrees)
majis_scan_angle_obj = MAJIS_Scan_Angle({}, silent=True)
# check that start scanner angle is +/- 2.0 degrees
start_angle = mc_angles[0]
if abs(start_angle) - majis_scan_angle_obj.scanner_angular_range / 2 != 0.0:
print(f'[WARNING] Invalid motion compensation start angle: {start_angle}. It must correspond to +/-2.0 degrees.')
return self.data
# Symmetrical scan
stop_angle = -start_angle
last_valid_line = -1
for i in range(len(mc_angles) - 1):
if stop_angle < 0:
if (mc_angles[i] > stop_angle) and (mc_angles[i+1] < stop_angle):
last_valid_line = i
else:
if (mc_angles[i] < stop_angle) and (mc_angles[i+1] > stop_angle):
last_valid_line = i
# No symmetrical scan is possible
if last_valid_line < 0:
return self.data
# Compute valid scanner angles and positions corresponding to input optimal motion compensation (MC) angles
valid_mc_scanner_angles = []
valid_mc_scanner_positions = []
for mc_angle in mc_angles[0:last_valid_line + 2]:
valid_angle, valid_position = majis_scan_angle_obj.get_valid_scanner_angle(mc_angle)
valid_mc_scanner_angles.append(valid_angle)
valid_mc_scanner_positions.append(valid_position)
valid_mc_scanner_angles = np.array(valid_mc_scanner_angles)
valid_mc_scanner_positions = np.array(valid_mc_scanner_positions)
# Estimate optimal scan duration and update last acquisition time
scan_duration = np.interp(stop_angle, mc_angles, times, period=360)
mc_scanner_times = times[0:last_valid_line + 2].copy()
mc_scanner_times[-1] = scan_duration
# Compute duration of each optimal acquisition
delta_mc_acq_times = mc_scanner_times[1:] - mc_scanner_times[0:-1]
# Compute angular variation between each acquisition
valid_delta_mc_scanner_angles = valid_mc_scanner_angles[1:] - valid_mc_scanner_angles[0:-1]
# Compute optimal scanner angular velocity
valid_mc_scanner_speeds = valid_delta_mc_scanner_angles / delta_mc_acq_times
# Compute start/stop scanner angular velocities based on a first order polynomial fit
coefs = Polynomial.fit(valid_mc_scanner_positions[0:-1], valid_mc_scanner_speeds, 1).convert().coef
fit_scan_rates = coefs[0] + coefs[1] * valid_mc_scanner_positions
start_scan_rate = fit_scan_rates[0]
stop_scan_rate = fit_scan_rates[-1]
# Compute number of acquisitions fitting in the scan duration
t_int = mc_scanner_times[1]-mc_scanner_times[0]
scan_lines = int(np.floor(scan_duration / t_int))
self.data = [start_angle, start_scan_rate, stop_scan_rate, scan_duration, scan_lines]
else:
self.data = None
self.valid = False
print('Missing GeometryEvent object containing this geometry.')
return self.data
[docs]class Obs_Limb_Latitudinal_Range(Vector):
"""Observable limb latitudinal range during an observation opportunity.
"""
def __init__(self, parameters, geoevt=None, silent=None):
"""Constructor.
Args:
parameters:
geoevt:
silent:
"""
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.prefix = ['min', 'max']
self.data = [None, None]
self.units = 'deg'
self.set_parameters(parameters)
[docs] def compute(self, et):
min_lat = max(self.geoevt.get_sub_events_geometry_data('Limb_Min_Latitude'))
max_lat = min(self.geoevt.get_sub_events_geometry_data('Limb_Max_Latitude'))
if min_lat < max_lat:
self.data = [min_lat, max_lat]
return self.data
[docs]class Obs_Min_Phase_Angle(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.set_parameters(parameters)
[docs] def compute(self, et):
if self.geoevt is not None:
phase_angles = self.geoevt.get_sub_events_geometry_data('Phase_Angle')
self.data = min(phase_angles)
else:
self.data = None
self.valid = False
print('Missing GeometryEvent object containing this geometry.')
return self.data
[docs]class Obs_Max_Phase_Angle(Scalar):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.units = 'deg'
self.set_parameters(parameters)
[docs] def compute(self, et):
if self.geoevt is not None:
phase_angles = self.geoevt.get_sub_events_geometry_data('Phase_Angle')
self.data = max(phase_angles)
else:
self.data = None
self.valid = False
print('Missing GeometryEvent object containing this geometry.')
return self.data
[docs]class Obs_Ref_SubSC_Position(Vector):
def __init__(self, parameters, geoevt=None, silent=None):
super().__init__(parameters, geoevt=geoevt, silent=silent)
self.prefix = ['x', 'y', 'z']
self.data = [None] * 3
self.units = 'km'
self.set_parameters(parameters)
[docs] def compute(self, et):
ref = 0 # temporary
sub_sub_sc_positions = []
if self.geoevt != None:
sub_sub_sc_positions = self.geoevt.get_sub_events_geometry_data('SubSC_Point')
self.data = sub_sub_sc_positions[ref]
else:
self.data = None
self.valid = False
print('Missing GeometryEvent object containing this geometry.')
return self.data
