Planet Viewing Example¶
This example shows how to plot geometric quantities using the frames system.
First we import numpy and the plotting module.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
Then we import GODOT modules and create a universe from configuration.
The universe file contains an entry to insert a keplerian orbit point, named 'SC' into the frames system.
from godot.core import tempo, astro
from godot import cosmos
# optionally avoid verbose logging messages
import godot.core.util as util
util.suppressLogger()
# create the universe
uni_config = cosmos.util.load_yaml('universe.yml')
uni = cosmos.Universe(uni_config)
To see what is in the universe file:
print (open('universe.yml').read())
We can check that the SC, and some planets are contained in the frames system:
for p in uni.frames.listPointNames():
print(f"Frames contains point {p}")
for a in uni.frames.listAxesNames():
print(f"Frames contains axes {a}")
Occultation Times¶
We want to compute the times when Mars is in occultation due to the Earth, as viewed from the spacecraft.
We can extract from the universe a time evaluable function for this purpose. In the universe file there are geometry functions added for occultations:
uni.evaluables.contains("SC_Mars_Earth_Occultation_Geometric")
We can plot this function:
# specify a time grid
start = tempo.Epoch('2022-04-01T00:00:00 TDB')
end = tempo.Epoch('2022-05-01T00:00:00 TDB')
ran = tempo.EpochRange( start, end )
grid = ran.createGrid(60.0) # 60 seconds stepsize
# get the time evaluable for the occulation distance of Mars by Earth
model = uni.evaluables["SC_Mars_Earth_Occultation_Geometric"]
min_value = 200
fig = plt.figure(figsize=(9, 5))
plt.title('Mars Earth Occultation distance')
plt.xlabel('Time (MJD)')
plt.ylabel('Mars Earth Distance (km)')
times = [e.mjd() for e in grid]
values = [model.eval(e) for e in grid]
plt.plot(times, values )
plt.plot(times,[min_value]*len(times),'--')
plt.tight_layout()
It looks like this function drops below zero in this time period, which means that Mars is occulted by the Earth. We can find the exact times of these occultations as follows:
Create an event generator for when this function drops below a minimum value (zero here)
from godot.model import eventgen
model = uni.evaluables["SC_Mars_Earth_Occultation_Geometric"]
min_value = 200
func = model - min_value
eps = 1e-6
tol = 1e-6
generator = eventgen.EventGenerator(-func, eps, tol) # use negative to compute periods with no visibility
Then we can compute the time periods when mars is not visible due to occultations by the Earth:
event_grid = ran.contract(eps).createGrid(60.0)
mars_occult = generator.computeEventIntervals(event_grid)
for interval in mars_occult:
start = interval.start().value()
end = interval.end().value()
print( f"Occultation from {start} to {end}, duration {end - start} secs" )
Instrument Observations¶
Now we look at the instrument boresight, which is the frames Axes 'SC_LOF' oriented cross-radial,along track. We assume the boresight is in the along track direction. We use the existing time grid.
from godot.model import common
class AnglefromInstrumentBoresight(common.ScalarTimeEvaluable):
def __init__(self,body):
common.ScalarTimeEvaluable.__init__(self)
self.__body = body
def eval(self, epoch):
vec = uni.frames.vector3('SC',self.__body,'SC_LOF',epoch)
# compute the colatitude in the RCA frame (z-axis is along track)
return np.pi/2 - astro.sphericalFromCart(vec)[2]
model = AnglefromInstrumentBoresight('MarsBarycenter')
boresight_angle = 20.0
fig = plt.figure(figsize=(9, 5))
plt.title('Mars Angle From Instrument Boresight')
plt.xlabel('Time (MJD)')
plt.ylabel('Mars Angle From Boresight (deg)')
times = [e.mjd() for e in grid]
values = np.rad2deg([model.eval(e) for e in grid])
plt.plot(times, values)
plt.plot(times,[boresight_angle]*len(times),'--')
plt.tight_layout()
The boresight angle to Mars drops below 20 deg, let's find the times when this happens:
# func is positive when Mars is in the instrument boresight
func = np.deg2rad(boresight_angle )- model
eps = 1e-6
tol = 1e-6
generator = eventgen.EventGenerator(func, eps, tol)
event_grid = ran.contract(eps).createGrid(600.0)
mars_in_fov = generator.computeEventIntervals(event_grid)
for interval in mars_in_fov:
print(f"From {interval.start().value()} to {interval.end().value()} lasting {interval.width()} seconds")
Now plot the events on top of the function plot
# extract the start and end times of the incursion into the FOV in plottable form
starts = [entry.start().value().mjd() for entry in mars_in_fov if entry.start().isSolved()]
ends = [entry.end().value().mjd() for entry in mars_in_fov if entry.end().isSolved()]
# generate the plots
fig = plt.figure(figsize=(9, 5))
plt.title('Mars Angle From Instrument Boresight')
plt.xlabel('Time (MJD)')
plt.ylabel('Mars Angle From Boresight (deg)')
times = [e.mjd() for e in grid]
values = np.rad2deg([model.eval(e) for e in grid])
plt.plot(times, values)
plt.plot(times,[boresight_angle]*len(times),'--')
plt.plot(starts, [boresight_angle] * len(starts), 'o', label="entry")
plt.plot(ends, [boresight_angle] * len(ends), 'o', label="exit")
# plot a single entry, for clarity
plt.ylim(boresight_angle-2.5 ,boresight_angle+2.5)
plt.xlim(8126,8129)
plt.tight_layout()
We want to know when Mars is in the boresight, but not in occulation by the Earth
from godot.core import events
# we could compute mars visibility event interval sets (not in occultation)
mars_visib = mars_occult.complement()
# find all event interval set with mars in occultation and in instrument fov
mars_occulted_and_in_fov = events.overlap([mars_in_fov, mars_occult])
print('MARS OCCULTED AND IN FOV')
for interval in mars_occulted_and_in_fov:
print(f"From {interval.start().value()} to {interval.end().value()} lasting {interval.width()} seconds")