forked from astropy/astropy
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathplot_mars-coordinate-frame.py
More file actions
152 lines (110 loc) · 5.13 KB
/
plot_mars-coordinate-frame.py
File metadata and controls
152 lines (110 loc) · 5.13 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
r"""
============================================
Create a new coordinate frame class for Mars
============================================
This example describes how to subclass and define a custom coordinate frame for a
planetary body which can be described by a geodetic or bodycentric representation,
as discussed in :ref:`astropy:astropy-coordinates-design` and
:ref:`astropy-coordinates-create-geodetic`.
Note that we use the frame here only to store coordinates. To use it to determine, e.g.,
where to point a telescope on Earth to observe Olympus Mons, one would need to add the
frame to the transfer graph, which is beyond the scope of this example.
To do this, first we need to define a subclass of a
`~astropy.coordinates.BaseGeodeticRepresentation` and
`~astropy.coordinates.BaseBodycentricRepresentation`, then a subclass of
`~astropy.coordinates.BaseCoordinateFrame` using the previous defined
representations.
*By: Chiara Marmo, Marten van Kerkwijk*
*License: BSD*
"""
##############################################################################
# Set up numpy, matplotlib, and use a nicer set of plot parameters:
import matplotlib.pyplot as plt
import numpy as np
from astropy.visualization import astropy_mpl_style, quantity_support
plt.style.use(astropy_mpl_style)
quantity_support()
##############################################################################
# Import the packages necessary for coordinates
import astropy.units as u
from astropy.coordinates.baseframe import BaseCoordinateFrame
from astropy.coordinates.representation import CartesianRepresentation
from astropy.coordinates.representation.geodetic import (
BaseBodycentricRepresentation,
BaseGeodeticRepresentation,
)
##############################################################################
# The first step is to create a new class, and make it a subclass of
# `~astropy.coordinates.BaseGeodeticRepresentation`.
# Geodetic latitudes are used and longitudes span from 0 to 360 degrees east positive
# It represent a best fit of the Mars spheroid to the martian geoid (areoid):
class MarsBestFitAeroid(BaseGeodeticRepresentation):
"""
A Spheroidal representation of Mars that minimized deviations with respect to the
areoid following
Ardalan A. A, R. Karimi, and E. W. Grafarend (2010)
https://doi.org/10.1007/s11038-009-9342-7
"""
_equatorial_radius = 3395.4280 * u.km
_flattening = 0.5227617843759314 * u.percent
#####################################################################################
# Now let's define a new geodetic representation obtained from MarsBestFitAeroid but
# described by planetocentric latitudes.
class MarsBestFitOcentricAeroid(BaseBodycentricRepresentation):
"""
A Spheroidal planetocentric representation of Mars that minimized deviations with
respect to the areoid following
Ardalan A. A, R. Karimi, and E. W. Grafarend (2010)
https://doi.org/10.1007/s11038-009-9342-7
"""
_equatorial_radius = 3395.4280 * u.km
_flattening = 0.5227617843759314 * u.percent
#############################################################################
# As a comparison we define a new spherical frame representation, we could
# have based it on `~astropy.coordinates.BaseBodycentricRepresentation` too.
class MarsSphere(BaseGeodeticRepresentation):
"""
A Spherical representation of Mars
"""
_equatorial_radius = 3395.4280 * u.km
_flattening = 0.0 * u.percent
#############################################################################
# The new planetary body-fixed reference system will be described using the
# previous defined representations.
class MarsCoordinateFrame(BaseCoordinateFrame):
"""
A reference system for Mars.
"""
name = "Mars"
#############################################################################
# Now we plot the differences between each component of the cartesian
# representation with respect to the spherical model, assuming the point on the
# surface of the body (``height = 0``)
mars_sphere = MarsCoordinateFrame(
lon=np.linspace(0, 360, 128) * u.deg,
lat=np.linspace(-90, 90, 128) * u.deg,
representation_type=MarsSphere,
)
mars = MarsCoordinateFrame(
lon=np.linspace(0, 360, 128) * u.deg,
lat=np.linspace(-90, 90, 128) * u.deg,
representation_type=MarsBestFitAeroid,
)
mars_ocentric = MarsCoordinateFrame(
lon=np.linspace(0, 360, 128) * u.deg,
lat=np.linspace(-90, 90, 128) * u.deg,
representation_type=MarsBestFitOcentricAeroid,
)
xyz_sphere = mars_sphere.represent_as(CartesianRepresentation)
xyz = mars.represent_as(CartesianRepresentation)
xyz_ocentric = mars_ocentric.represent_as(CartesianRepresentation)
fig, ax = plt.subplots(2, subplot_kw={"projection": "3d"})
ax[0].scatter(*((xyz - xyz_sphere).xyz << u.km))
ax[0].tick_params(labelsize=8)
ax[0].set(xlabel="x [km]", ylabel="y [km]", zlabel="z [km]")
ax[0].set_title("Mars-odetic spheroid difference from sphere")
ax[1].scatter(*((xyz_ocentric - xyz_sphere).xyz << u.km))
ax[1].tick_params(labelsize=8)
ax[1].set(xlabel="x [km]", ylabel="y [km]", zlabel="z [km]")
ax[1].set_title("Mars-ocentric spheroid difference from sphere")
plt.show()