mirror of
https://github.com/cosinekitty/astronomy.git
synced 2025-12-26 17:19:41 -05:00
53 lines
2.0 KiB
Python
Executable File
53 lines
2.0 KiB
Python
Executable File
#!/usr/bin/env python3
|
|
import math
|
|
|
|
class Vector:
|
|
def __init__(self, x, y, z):
|
|
self.x = x
|
|
self.y = y
|
|
self.z = z
|
|
|
|
def __str__(self):
|
|
return '[{:19.16f}, {:19.16f}, {:19.16f}]'.format(self.x, self.y, self.z)
|
|
|
|
def Normal(p, v):
|
|
# Compute the cross product of the position and velocity vectors.
|
|
nx = p.y*v.z - p.z*v.y
|
|
ny = p.z*v.x - p.x*v.z
|
|
nz = p.x*v.y - p.y*v.x
|
|
# Convert the normal to a unit vector.
|
|
mag = math.sqrt(nx*nx + ny*ny + nz*nz)
|
|
return Vector(nx/mag, ny/mag, nz/mag)
|
|
|
|
def ArcminBetween(a, b):
|
|
# Take the dot product of the two vectors.
|
|
dot = a.x*b.x + a.y*b.y + a.z*b.z
|
|
# The inverse cosine tells us the angle between them.
|
|
radians = math.acos(dot)
|
|
# Convert radians to arcminutes.
|
|
return 60.0 * math.degrees(radians)
|
|
|
|
# EMB-centered Moon state vector from JPL Horizons:
|
|
#2458932.750000000 = A.D. 2020-Mar-24 06:00:00.0000 TDB
|
|
# X = 2.673075076694347E-03 Y = 1.870517968021612E-04 Z =-1.731476992383482E-04
|
|
# VX=-1.967647146277114E-05 VY= 5.085581013947419E-04 VZ= 2.191688329215286E-04
|
|
moon_pos = Vector( 2.673075076694347E-03, 1.870517968021612E-04, -1.731476992383482E-04)
|
|
moon_vel = Vector(-1.967647146277114E-05, 5.085581013947419E-04, 2.191688329215286E-04)
|
|
|
|
# EMB-centered L4 state vector from JPL Horizons:
|
|
#2458932.750000000 = A.D. 2020-Mar-24 06:00:00.0000 TDB
|
|
# X = 1.229291654700910E-03 Y = 2.252359994857843E-03 Z = 8.460455897080097E-04
|
|
# VX=-4.932901244670627E-04 VY= 2.183360032469220E-04 VZ= 1.409671573944061E-04
|
|
L4_pos = Vector( 1.229291654700910E-03, 2.252359994857843E-03, 8.460455897080097E-04)
|
|
L4_vel = Vector(-4.932901244670627E-04, 2.183360032469220E-04, 1.409671573944061E-04)
|
|
|
|
moon_norm = Normal(moon_pos, moon_vel)
|
|
print('moon_norm = ', moon_norm)
|
|
|
|
L4_norm = Normal(L4_pos, L4_vel)
|
|
print('L4_norm = ', L4_norm)
|
|
|
|
# Find the angle between the two normal vectors.
|
|
arcmin = ArcminBetween(moon_norm, L4_norm)
|
|
print('Angle between vectors = {:0.6f} arcmin'.format(arcmin))
|