<p>It looks like you're trying to find the Moon's ascending node from position and velocity vectors as in <a href="https://en.wikipedia.org/wiki/Longitude_of_the_ascending_node#Calculation_from_state_vectors" rel="nofollow noreferrer">this Wikipedia article</a>.
<strong>h</strong> = <strong>r</strong> × <strong>v</strong> is perpendicular to the Moon's orbital plane, and its projection onto the ecliptic plane (ignoring <em>h<sub>z</sub></em>) is 90° behind the ascending node.</p>
<p>These lines of your HORIZONS output show a coordinate center in southern New Mexico, possibly from looking up 'earth' instead of '@earth'.</p>
<pre><code>Center-site name: Heaven on Earth Observatory, Mayhill
⋮
Center geodetic : 254.471, 32.9035384, 2.23607    {E-lon(deg),Lat(deg),Alt(km)}
</code></pre>
<p>This puts <strong>r</strong> off by 6000 km, <strong>v</strong> off by 0.4 km/s, and the computed longitude off by several degrees.
With coordinate center Geocentric [code: 500], HORIZONS gives this output:</p>
<pre><code>2460215.500000000 = A.D. 2023-Sep-28 00:00:00.0000 TDB 
 X = 3.461836628403760E+05 Y =-9.613409075264385E+04 Z =-2.126911943880732E+04
 VX= 2.954176491446214E-01 VY= 1.050043302428879E+00 VZ= 7.630639995108085E-02
</code></pre>
<p>Your cross product implementation looks correct,
but you don't need the magnitude just to compute an angle.
I recommend <a href="https://numpy.org/" rel="nofollow noreferrer">numpy</a> for vector math.</p>
<pre><code>import math
import numpy

# HORIZONS geocentric state vectors for the Moon at 2023-09-28 0:00 TDB
r = numpy.array([3.461837e+5, -9.613409e+4, -2.126912e+4])
v = numpy.array([2.954176e-1, 1.050043, 7.630640e-2])
h = numpy.cross(r, v)

# J2000 longitude of ascending node
omega = math.degrees(math.atan2(h[1], h[0])) + 90
if omega &lt; 0:
    omega += 360

print('omega = {0:.2f}'.format(omega))
</code></pre>
<p>gives this output:</p>
<pre><code>omega = 24.64
</code></pre>
<p>which agrees with OM in HORIZONS's osculating orbital elements:</p>
<pre><code>2460215.500000000 = A.D. 2023-Sep-28 00:00:00.0000 TDB 
 EC= 6.651987996774227E-02 QR= 3.599111933878179E+05 IN= 5.244703685174947E+00
 OM= 2.463901296089877E+01 W = 3.204209825950677E+02 Tp=  2460215.546430446208
 N = 1.520234985695220E-04 MA= 3.593901439672831E+02 TA= 3.593016826209310E+02
</code></pre>
<p>As PM 2Ring commented, that's relative to the J2000 equinox.
For 2023-09-28, add 0.33° for <a href="https://en.wikipedia.org/wiki/Axial_precession#Values" rel="nofollow noreferrer">precession</a>,
making Ω = 24.97° relative to the equinox of date.</p>

CC BY-SA 4.0
2023-10-02T12:30:43.393Z
2023-10-02T12:30:43.393Z
Mike G
<p>I have to be honest. Been learning python for few months  but when I saw this code I was blown away, mainly by its simplicity  and secondly by not understanding it .. can someone help me understand how to use &quot;h_x &amp;&amp; h_y&quot; ?  I have no clue where to start looking for this info. Thank you &lt;3</p>
<pre><code>import math
 
def calculate_longitude_of_ascending_node(h_x, h_y, h_z):
    # Calculate the magnitude of the specific relative angular momentum vector (|h|)
    h_magnitude = math.sqrt(h_x**2 + h_y**2 + h_z**2)
 
    # Calculate the unit vector n using the cross product of r and v
    n_x = h_y * h_z
    n_y = -h_x * h_z
    n_z = h_x**2 + h_y**2
 
    # Calculate the longitude of the ascending node (Ω) in radians
    omega = math.atan2(n_y, n_x)
    if omega &lt; 0:
        omega += 2 * math.pi
 
    # Convert Ω from radians to degrees, minutes, and seconds
    omega_deg = math.degrees(omega)
    omega_deg_int = int(omega_deg)
    omega_min = (omega_deg - omega_deg_int) * 60
    omega_min_int = int(omega_min)
    omega_sec = (omega_min - omega_min_int) * 60
 
    return omega_deg_int, omega_min_int, omega_sec
 
# Example usage:
# Use the provided position data
x = 3.470705658033574E+05
y = -9.265520886744432E+04
z = -2.653599067380910E+04
vx = -9.011403192706635E-02
vy = 1.110259691085033E+00
vz = 5.115987924302562E-02
 
# Calculate h_x, h_y, h_z by taking the cross product of r and v
h_x = y * vz - z * vy
h_y = z * vx - x * vz
h_z = x * vy - y * vx
 
# Calculate the longitude of the ascending node in DMS
longitude_deg, longitude_min, longitude_sec = calculate_longitude_of_ascending_node(h_x, h_y, h_z)
print(f&quot;Longitude of Ascending Node (Ω): {longitude_deg}° {longitude_min}' {longitude_sec:.2f}\&quot;&quot;)
</code></pre>
<p>edit***  so I was able to get the x, y and z from Horizon</p>
<pre><code>*******************************************************************************
 Revised: July 31, 2013             Moon / (Earth)                          301
 
 GEOPHYSICAL DATA (updated 2018-Aug-15):
  Vol. mean radius, km  = 1737.53+-0.03    Mass, x10^22 kg       =    7.349
  Radius (gravity), km  = 1738.0           Surface emissivity    =    0.92
  Radius (IAU), km      = 1737.4           GM, km^3/s^2          = 4902.800066
  Density, g/cm^3       =    3.3437        GM 1-sigma, km^3/s^2  =  +-0.0001  
  V(1,0)                =   +0.21          Surface accel., m/s^2 =    1.62
  Earth/Moon mass ratio = 81.3005690769    Farside crust. thick. = ~80 - 90 km
  Mean crustal density  = 2.97+-.07 g/cm^3 Nearside crust. thick.= 58+-8 km 
  Heat flow, Apollo 15  = 3.1+-.6 mW/m^2   Mean angular diameter = 31'05.2&quot;
  Heat flow, Apollo 17  = 2.2+-.5 mW/m^2   Sid. rot. rate, rad/s = 0.0000026617
  Geometric Albedo      = 0.12             Mean solar day        = 29.5306 d
  Obliquity to orbit    = 6.67 deg         Orbit period          = 27.321582 d
  Semi-major axis, a    = 384400 km        Eccentricity          = 0.05490
  Mean motion, rad/s    = 2.6616995x10^-6  Inclination           = 5.145 deg
  Apsidal period        = 3231.50 d        Nodal period          = 6798.38 d
                                 Perihelion  Aphelion    Mean
  Solar Constant (W/m^2)         1414+-7     1323+-7     1368+-7
  Maximum Planetary IR (W/m^2)   1314        1226        1268
  Minimum Planetary IR (W/m^2)      5.2         5.2         5.2
********************************************************************************
 
 
*******************************************************************************
Ephemeris / WWW_USER Fri Sep 29 12:13:59 2023 Pasadena, USA      / Horizons
*******************************************************************************
Target body name: Moon (301)                      {source: DE441}
Center body name: Earth (399)                     {source: DE441}
Center-site name: Heaven on Earth Observatory, Mayhill
*******************************************************************************
Start time      : A.D. 2023-Sep-28 00:00:00.0000 TDB
Stop  time      : A.D. 2023-Oct-01 00:00:00.0000 TDB
Step-size       : 1440 minutes
*******************************************************************************
Center geodetic : 254.471, 32.9035384, 2.23607    {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 254.471, 5362.17112, 3446.19639 {E-lon(deg),Dxy(km),Dz(km)}
Center pole/equ : ITRF93                          {East-longitude positive}
Center radii    : 6378.137, 6378.137, 6356.752 km {Equator_a, b, pole_c}
Output units    : KM-S
Calendar mode   : Mixed Julian/Gregorian
Output type     : GEOMETRIC cartesian states
Output format   : 3 (position, velocity, LT, range, range-rate)
EOP file        : eop.230927.p231221
EOP coverage    : DATA-BASED 1962-JAN-20 TO 2023-SEP-27. PREDICTS-&gt; 2023-DEC-20
Reference frame : Ecliptic of J2000.0
*******************************************************************************
JDTDB
   X     Y     Z
   VX    VY    VZ
   LT    RG    RR
*******************************************************************************
<span class="math-container">$$SOE
2460215.500000000 = A.D. 2023-Sep-28 00:00:00.0000 TDB 
 X = 3.470705658033574E+05 Y =-9.265520886744432E+04 Z =-2.653599067380910E+04
 VX=-9.011403192706635E-02 VY= 1.110259691085033E+00 VZ= 5.115987924302562E-02
 LT= 1.201512151373686E+00 RG= 3.602042811771854E+05 RR=-3.761888927526311E-01
2460216.500000000 = A.D. 2023-Sep-29 00:00:00.0000 TDB 
 X = 3.612178246289230E+05 Y = 2.175102251481885E+02 Z =-1.933391485093577E+04
 VX=-3.538761876322111E-01 VY= 1.141964183533681E+00 VZ= 6.712048003763832E-02
 LT= 1.206617868501206E+00 RG= 3.617349366646974E+05 RR=-3.562710879736578E-01
2460217.500000000 = A.D. 2023-Sep-30 00:00:00.0000 TDB 
 X = 3.525900479296140E+05 Y = 9.334101334405206E+04 Z =-1.122942013890661E+04
 VX=-6.159050044487818E-01 VY= 1.104174316562248E+00 VZ= 7.721989291731685E-02
 LT= 1.217204527827483E+00 RG= 3.649087372861304E+05 RR=-3.150496228123838E-01
2460218.500000000 = A.D. 2023-Oct-01 00:00:00.0000 TDB 
 X = 3.221837310026606E+05 Y = 1.808849216335789E+05 Z =-2.748397019656826E+03
 VX=-8.568131300802702E-01 VY= 1.002000203732668E+00 VZ= 8.112992260778490E-02
 LT= 1.232514886108716E+00 RG= 3.694986672281219E+05 RR=-2.571795472029340E-01
$$</span>EOE
*******************************************************************************
 
TIME
 
  Barycentric Dynamical Time (&quot;TDB&quot; or T_eph) output was requested. This
continuous coordinate time is equivalent to the relativistic proper time
of a clock at rest in a reference frame co-moving with the solar system
barycenter but outside the system's gravity well. It is the independent
variable in the solar system relativistic equations of motion.
 
  TDB runs at a uniform rate of one SI second per second and is independent
of irregularities in Earth's rotation.
 
CALENDAR SYSTEM
 
  Mixed calendar mode was active such that calendar dates after AD 1582-Oct-15
(if any) are in the modern Gregorian system. Dates prior to 1582-Oct-5 (if any)
are in the Julian calendar system, which is automatically extended for dates
prior to its adoption on 45-Jan-1 BC.  The Julian calendar is useful for
matching historical dates. The Gregorian calendar more accurately corresponds
to the Earth's orbital motion and seasons. A &quot;Gregorian-only&quot; calendar mode is
available if such physical events are the primary interest.
 
REFERENCE FRAME AND COORDINATES
 
  Ecliptic at the standard reference epoch
 
    Reference epoch: J2000.0
    X-Y plane: adopted Earth orbital plane at the reference epoch
               Note: IAU76 obliquity of 84381.448 arcseconds wrt ICRF X-Y plane
    X-axis   : ICRF
    Z-axis   : perpendicular to the X-Y plane in the directional (+ or -) sense
               of Earth's north pole at the reference epoch.
 
  Symbol meaning:
 
    JDTDB    Julian Day Number, Barycentric Dynamical Time
      X      X-component of position vector (km)
      Y      Y-component of position vector (km)
      Z      Z-component of position vector (km)
      VX     X-component of velocity vector (km/sec)                           
      VY     Y-component of velocity vector (km/sec)                           
      VZ     Z-component of velocity vector (km/sec)                           
      LT     One-way down-leg Newtonian light-time (sec)
      RG     Range; distance from coordinate center (km)
      RR     Range-rate; radial velocity wrt coord. center (km/sec)
 
ABERRATIONS AND CORRECTIONS
</code></pre>
<p>and came up with: &quot;python nodemath.py
Longitude of Ascending Node (Ω): 4.156300512180232&quot; .. now I am not quite sure its correct or I got it right. This is not matching with the current NN coordinates so the code is off</p>

CC BY-SA 4.0
2023-09-29T15:49:56.283Z
Would need to see where you got the code to make judgments on how it's intended to work, and if `h_x` and `h_y` are the x- and y- components of the specific relative angular momentum vector, it looks like the calculation to get its magnitude is wrong, unless the z- component of said vector is zero.
notovny
CC BY-SA 4.0
2023-09-29T16:26:51.850Z
As notovny said, that code is a bit weird. The [specific angular momentum](https://en.wikipedia.org/wiki/Specific_angular_momentum) vector, $\mathbf{ h=r\times v}$ is normal to the orbit plane. And if your XY plane is the ecliptic, then the z component of the normal to the Moon's orbit plane isn't zero.
PM 2Ring
CC BY-SA 4.0
2023-09-29T16:34:55.973Z
FWIW, if the h vector is (a, b, c), the orbit plane intersects the XY plane in the line $ax+by=0$, which is the line of nodes. Which means $\tan(\Omega)$ is $(-a/b)$ or $(a/-b)$.
PM 2Ring
CC BY-SA 4.0
2023-09-29T18:09:16.307Z
Ah, that's the issue... Looking at Wikipedia's [Longitude of the Ascending Node](https://en.wikipedia.org/wiki/Longitude_of_the_ascending_node) article, the code has misnamed and incorrectly described a variable. `h_magnitude` should be called `n_magnitude`, because it's the magnitude of the vector $\textbf{n}$  (I don't know if it has a customary name, I usually call it the node vector), which points from the central body towards the Ascending Node.  Also, in the code as above, as in the wikipedia article, $\textbf{n}$ is _not_ a unit vector.
notovny
CC BY-SA 4.0
2023-09-29T19:21:13.357Z
lets see if it allows me to post a link https://pastebin.com/7rwzgNLD
dimitri33
CC BY-SA 4.0
2023-09-29T21:46:45.993Z
The original code was functionally correct, it's just that one of the variables was mislabeled, (the value it contained, and the equations it was being used in were fine), and the comments were wrong about the properties they described. No comment about the current version.
notovny
CC BY-SA 4.0
2023-09-30T00:35:09.177Z
On the Horizons Web App, just above the Generate Ephemeris button, there's a link where it says "the settings above can be viewed by using this link". You can paste those settings into your question (in a code block). I have Python code that sends Batch Data to Horizons.https://gist.github.com/PM2Ring/b1fec75e78cc08f6fc28c6f6c43529c3 Or just click on the link in the next comment.
PM 2Ring
CC BY-SA 4.0
2023-09-30T00:35:17.167Z
https://sagecell.sagemath.org/?z=eJxdj0FLwzAUx--FfodHPLTRktUyQQpF8SB6EMSrmyNr0zYjTWLyqmO67266ehADIeT9815-vzN4ME4ejPaw5Vj3DJ6foHiRuoMiLwp2WbD8Ko7kYI1DcCIL-30UHn0cjU5BBaRHtL5cLLxv2M4qprnnrDMfC27lov-dvmmlEixUSBxZjqHPCVabwYZy6kh6U76-JSu__k7CuT5P6AWhcXQrNQrHa4yjRrTQmfQEWT1qO-Kd2aefssG-WuYZ9EJ2PVZFTjPgI5rNaBuOorrnygtaxhGEJVvQBmfTMhDg6PScBKsT0-zGrPGYBr8MwhBefZHWuIEjKYGg2CM5ZjAJ-ZDIiSUEKamHhmRAVpqwnZE6DZ6slbrhSs3clNIjnf_7YzZdJ7tNevgvdv1HbJlNr_iosAqYbMLIQPGtUFWSUFqC5d7_ABQnk8E=&lang=sage
PM 2Ring
CC BY-SA 4.0
2023-09-30T00:59:52.017Z
The section of your code just under the comment "#Calculate the unit vector" is still wacky. But it doesn't matter, because the ratio `n_y/n_x` is correct, so you should get the correct omega. Note that the Horizons vector data is relative to the J2000 equinox, not the current equinox (aka the equinox of date). The XY plane should be ok, but the direction of the X axis (the First Point of Aries) in equinox of date coordinates has shifted slightly since the year 2000, due to precession.
PM 2Ring
CC BY-SA 4.0
2023-09-30T01:08:33.763Z
Also, you should be using the centre of the Earth as your origin, `@399` or `500@399`, not some observatory.
PM 2Ring
CC BY-SA 4.0
2023-09-30T06:10:47.867Z
@PM2Ring When a user looks up "earth" without a leading "@", Horizons unfortunately chooses Heaven on Earth Observatory of Mayhill NM instead of the geocenter.
Mike G
<the-moon><natural-satellites>
Determining the position of the North Node using math?
dimitri33
# Determining the position of the North Node using math?

(Asked by: username_0 on 2023-09-29)

I have to be honest. Been learning python for few months  but when I saw this code I was blown away, mainly by its simplicity  and secondly by not understanding it .. can someone help me understand how to use "h_x && h_y" ?  I have no clue where to start looking for this info. Thank you <3
import math
 
def calculate_longitude_of_ascending_node(h_x, h_y, h_z):
    # Calculate the magnitude of the specific relative angular momentum vector (|h|)
    h_magnitude = math.sqrt(h_x**2 + h_y**2 + h_z**2)
 
    # Calculate the unit vector n using the cross product of r and v
    n_x = h_y * h_z
    n_y = -h_x * h_z
    n_z = h_x**2 + h_y**2
 
    # Calculate the longitude of the ascending node (Ω) in radians
    omega = math.atan2(n_y, n_x)
    if omega < 0:
        omega += 2 * math.pi
 
    # Convert Ω from radians to degrees, minutes, and seconds
    omega_deg = math.degrees(omega)
    omega_deg_int = int(omega_deg)
    omega_min = (omega_deg - omega_deg_int) * 60
    omega_min_int = int(omega_min)
    omega_sec = (omega_min - omega_min_int) * 60
 
    return omega_deg_int, omega_min_int, omega_sec
 
# Example usage:
# Use the provided position data
x = 3.470705658033574E+05
y = -9.265520886744432E+04
z = -2.653599067380910E+04
vx = -9.011403192706635E-02
vy = 1.110259691085033E+00
vz = 5.115987924302562E-02
 
# Calculate h_x, h_y, h_z by taking the cross product of r and v
h_x = y * vz - z * vy
h_y = z * vx - x * vz
h_z = x * vy - y * vx
 
# Calculate the longitude of the ascending node in DMS
longitude_deg, longitude_min, longitude_sec = calculate_longitude_of_ascending_node(h_x, h_y, h_z)
print(f"Longitude of Ascending Node (Ω): {longitude_deg}° {longitude_min}' {longitude_sec:.2f}\"")

edit***  so I was able to get the x, y and z from Horizon
*******************************************************************************
 Revised: July 31, 2013             Moon / (Earth)                          301
 
 GEOPHYSICAL DATA (updated 2018-Aug-15):
  Vol. mean radius, km  = 1737.53+-0.03    Mass, x10^22 kg       =    7.349
  Radius (gravity), km  = 1738.0           Surface emissivity    =    0.92
  Radius (IAU), km      = 1737.4           GM, km^3/s^2          = 4902.800066
  Density, g/cm^3       =    3.3437        GM 1-sigma, km^3/s^2  =  +-0.0001  
  V(1,0)                =   +0.21          Surface accel., m/s^2 =    1.62
  Earth/Moon mass ratio = 81.3005690769    Farside crust. thick. = ~80 - 90 km
  Mean crustal density  = 2.97+-.07 g/cm^3 Nearside crust. thick.= 58+-8 km 
  Heat flow, Apollo 15  = 3.1+-.6 mW/m^2   Mean angular diameter = 31'05.2"
  Heat flow, Apollo 17  = 2.2+-.5 mW/m^2   Sid. rot. rate, rad/s = 0.0000026617
  Geometric Albedo      = 0.12             Mean solar day        = 29.5306 d
  Obliquity to orbit    = 6.67 deg         Orbit period          = 27.321582 d
  Semi-major axis, a    = 384400 km        Eccentricity          = 0.05490
  Mean motion, rad/s    = 2.6616995x10^-6  Inclination           = 5.145 deg
  Apsidal period        = 3231.50 d        Nodal period          = 6798.38 d
                                 Perihelion  Aphelion    Mean
  Solar Constant (W/m^2)         1414+-7     1323+-7     1368+-7
  Maximum Planetary IR (W/m^2)   1314        1226        1268
  Minimum Planetary IR (W/m^2)      5.2         5.2         5.2
********************************************************************************
 
 
*******************************************************************************
Ephemeris / WWW_USER Fri Sep 29 12:13:59 2023 Pasadena, USA      / Horizons
*******************************************************************************
Target body name: Moon (301)                      {source: DE441}
Center body name: Earth (399)                     {source: DE441}
Center-site name: Heaven on Earth Observatory, Mayhill
*******************************************************************************
Start time      : A.D. 2023-Sep-28 00:00:00.0000 TDB
Stop  time      : A.D. 2023-Oct-01 00:00:00.0000 TDB
Step-size       : 1440 minutes
*******************************************************************************
Center geodetic : 254.471, 32.9035384, 2.23607    {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 254.471, 5362.17112, 3446.19639 {E-lon(deg),Dxy(km),Dz(km)}
Center pole/equ : ITRF93                          {East-longitude positive}
Center radii    : 6378.137, 6378.137, 6356.752 km {Equator_a, b, pole_c}
Output units    : KM-S
Calendar mode   : Mixed Julian/Gregorian
Output type     : GEOMETRIC cartesian states
Output format   : 3 (position, velocity, LT, range, range-rate)
EOP file        : eop.230927.p231221
EOP coverage    : DATA-BASED 1962-JAN-20 TO 2023-SEP-27. PREDICTS-> 2023-DEC-20
Reference frame : Ecliptic of J2000.0
*******************************************************************************
JDTDB
   X     Y     Z
   VX    VY    VZ
   LT    RG    RR
*******************************************************************************
$$SOE
2460215.500000000 = A.D. 2023-Sep-28 00:00:00.0000 TDB 
 X = 3.470705658033574E+05 Y =-9.265520886744432E+04 Z =-2.653599067380910E+04
 VX=-9.011403192706635E-02 VY= 1.110259691085033E+00 VZ= 5.115987924302562E-02
 LT= 1.201512151373686E+00 RG= 3.602042811771854E+05 RR=-3.761888927526311E-01
2460216.500000000 = A.D. 2023-Sep-29 00:00:00.0000 TDB 
 X = 3.612178246289230E+05 Y = 2.175102251481885E+02 Z =-1.933391485093577E+04
 VX=-3.538761876322111E-01 VY= 1.141964183533681E+00 VZ= 6.712048003763832E-02
 LT= 1.206617868501206E+00 RG= 3.617349366646974E+05 RR=-3.562710879736578E-01
2460217.500000000 = A.D. 2023-Sep-30 00:00:00.0000 TDB 
 X = 3.525900479296140E+05 Y = 9.334101334405206E+04 Z =-1.122942013890661E+04
 VX=-6.159050044487818E-01 VY= 1.104174316562248E+00 VZ= 7.721989291731685E-02
 LT= 1.217204527827483E+00 RG= 3.649087372861304E+05 RR=-3.150496228123838E-01
2460218.500000000 = A.D. 2023-Oct-01 00:00:00.0000 TDB 
 X = 3.221837310026606E+05 Y = 1.808849216335789E+05 Z =-2.748397019656826E+03
 VX=-8.568131300802702E-01 VY= 1.002000203732668E+00 VZ= 8.112992260778490E-02
 LT= 1.232514886108716E+00 RG= 3.694986672281219E+05 RR=-2.571795472029340E-01
$$EOE
*******************************************************************************
 
TIME
 
  Barycentric Dynamical Time ("TDB" or T_eph) output was requested. This
continuous coordinate time is equivalent to the relativistic proper time
of a clock at rest in a reference frame co-moving with the solar system
barycenter but outside the system's gravity well. It is the independent
variable in the solar system relativistic equations of motion.
 
  TDB runs at a uniform rate of one SI second per second and is independent
of irregularities in Earth's rotation.
 
CALENDAR SYSTEM
 
  Mixed calendar mode was active such that calendar dates after AD 1582-Oct-15
(if any) are in the modern Gregorian system. Dates prior to 1582-Oct-5 (if any)
are in the Julian calendar system, which is automatically extended for dates
prior to its adoption on 45-Jan-1 BC.  The Julian calendar is useful for
matching historical dates. The Gregorian calendar more accurately corresponds
to the Earth's orbital motion and seasons. A "Gregorian-only" calendar mode is
available if such physical events are the primary interest.
 
REFERENCE FRAME AND COORDINATES
 
  Ecliptic at the standard reference epoch
 
    Reference epoch: J2000.0
    X-Y plane: adopted Earth orbital plane at the reference epoch
               Note: IAU76 obliquity of 84381.448 arcseconds wrt ICRF X-Y plane
    X-axis   : ICRF
    Z-axis   : perpendicular to the X-Y plane in the directional (+ or -) sense
               of Earth's north pole at the reference epoch.
 
  Symbol meaning:
 
    JDTDB    Julian Day Number, Barycentric Dynamical Time
      X      X-component of position vector (km)
      Y      Y-component of position vector (km)
      Z      Z-component of position vector (km)
      VX     X-component of velocity vector (km/sec)                           
      VY     Y-component of velocity vector (km/sec)                           
      VZ     Z-component of velocity vector (km/sec)                           
      LT     One-way down-leg Newtonian light-time (sec)
      RG     Range; distance from coordinate center (km)
      RR     Range-rate; radial velocity wrt coord. center (km/sec)
 
ABERRATIONS AND CORRECTIONS

and came up with: "python nodemath.py
Longitude of Ascending Node (Ω): 4.156300512180232" .. now I am not quite sure its correct or I got it right. This is not matching with the current NN coordinates so the code is off

**Tags:** the-moon, natural-satellites
> **Comments:**
> **username_1**: Would need to see where you got the code to make judgments on how it's intended to work, and if `h_x` and `h_y` are the x- and y- components of the specific relative angular momentum vector, it looks like the calculation to get its magnitude is wrong, unless the z- component of said vector is zero.  
> **username_2**: As notovny said, that code is a bit weird. The [specific angular momentum](https://en.wikipedia.org/wiki/Specific_angular_momentum) vector, $\mathbf{ h=r\times v}$ is normal to the orbit plane. And if your XY plane is the ecliptic, then the z component of the normal to the Moon's orbit plane isn't zero.  
> **username_2**: FWIW, if the h vector is (a, b, c), the orbit plane intersects the XY plane in the line $ax+by=0$, which is the line of nodes. Which means $\tan(\Omega)$ is $(-a/b)$ or $(a/-b)$.  
> **username_1**: Ah, that's the issue... Looking at Wikipedia's [Longitude of the Ascending Node](https://en.wikipedia.org/wiki/Longitude_of_the_ascending_node) article, the code has misnamed and incorrectly described a variable. `h_magnitude` should be called `n_magnitude`, because it's the magnitude of the vector $\textbf{n}$  (I don't know if it has a customary name, I usually call it the node vector), which points from the central body towards the Ascending Node.  Also, in the code as above, as in the wikipedia article, $\textbf{n}$ is _not_ a unit vector.  
> **username_0**: lets see if it allows me to post a link https://pastebin.com/7rwzgNLD  
> **username_1**: The original code was functionally correct, it's just that one of the variables was mislabeled, (the value it contained, and the equations it was being used in were fine), and the comments were wrong about the properties they described. No comment about the current version.  
> **username_2**: On the Horizons Web App, just above the Generate Ephemeris button, there's a link where it says "the settings above can be viewed by using this link". You can paste those settings into your question (in a code block). I have Python code that sends Batch Data to Horizons.https://gist.github.com/PM2Ring/b1fec75e78cc08f6fc28c6f6c43529c3 Or just click on the link in the next comment.  
> **username_2**: https://sagecell.sagemath.org/?z=eJxdj0FLwzAUx--FfodHPLTRktUyQQpF8SB6EMSrmyNr0zYjTWLyqmO67266ehADIeT9815-vzN4ME4ejPaw5Vj3DJ6foHiRuoMiLwp2WbD8Ko7kYI1DcCIL-30UHn0cjU5BBaRHtL5cLLxv2M4qprnnrDMfC27lov-dvmmlEixUSBxZjqHPCVabwYZy6kh6U76-JSu__k7CuT5P6AWhcXQrNQrHa4yjRrTQmfQEWT1qO-Kd2aefssG-WuYZ9EJ2PVZFTjPgI5rNaBuOorrnygtaxhGEJVvQBmfTMhDg6PScBKsT0-zGrPGYBr8MwhBefZHWuIEjKYGg2CM5ZjAJ-ZDIiSUEKamHhmRAVpqwnZE6DZ6slbrhSs3clNIjnf_7YzZdJ7tNevgvdv1HbJlNr_iosAqYbMLIQPGtUFWSUFqC5d7_ABQnk8E=&lang=sage  
> **username_2**: The section of your code just under the comment "#Calculate the unit vector" is still wacky. But it doesn't matter, because the ratio `n_y/n_x` is correct, so you should get the correct omega. Note that the Horizons vector data is relative to the J2000 equinox, not the current equinox (aka the equinox of date). The XY plane should be ok, but the direction of the X axis (the First Point of Aries) in equinox of date coordinates has shifted slightly since the year 2000, due to precession.  
> **username_2**: Also, you should be using the centre of the Earth as your origin, `@399` or `500@399`, not some observatory.  
> **username_3**: @PM2Ring When a user looks up "earth" without a leading "@", Horizons unfortunately chooses Heaven on Earth Observatory of Mayhill NM instead of the geocenter.  

---

## 1 answers

### Answer by username_3 on 2023-10-02. Score: 1

It looks like you're trying to find the Moon's ascending node from position and velocity vectors as in this Wikipedia article.
h = r × v is perpendicular to the Moon's orbital plane, and its projection onto the ecliptic plane (ignoring hz) is 90° behind the ascending node.
These lines of your HORIZONS output show a coordinate center in southern New Mexico, possibly from looking up 'earth' instead of '@earth'.
Center-site name: Heaven on Earth Observatory, Mayhill
⋮
Center geodetic : 254.471, 32.9035384, 2.23607    {E-lon(deg),Lat(deg),Alt(km)}

This puts r off by 6000 km, v off by 0.4 km/s, and the computed longitude off by several degrees.
With coordinate center Geocentric [code: 500], HORIZONS gives this output:
2460215.500000000 = A.D. 2023-Sep-28 00:00:00.0000 TDB 
 X = 3.461836628403760E+05 Y =-9.613409075264385E+04 Z =-2.126911943880732E+04
 VX= 2.954176491446214E-01 VY= 1.050043302428879E+00 VZ= 7.630639995108085E-02

Your cross product implementation looks correct,
but you don't need the magnitude just to compute an angle.
I recommend numpy for vector math.
import math
import numpy

# HORIZONS geocentric state vectors for the Moon at 2023-09-28 0:00 TDB
r = numpy.array([3.461837e+5, -9.613409e+4, -2.126912e+4])
v = numpy.array([2.954176e-1, 1.050043, 7.630640e-2])
h = numpy.cross(r, v)

# J2000 longitude of ascending node
omega = math.degrees(math.atan2(h[1], h[0])) + 90
if omega < 0:
    omega += 360

print('omega = {0:.2f}'.format(omega))

gives this output:
omega = 24.64

which agrees with OM in HORIZONS's osculating orbital elements:
2460215.500000000 = A.D. 2023-Sep-28 00:00:00.0000 TDB 
 EC= 6.651987996774227E-02 QR= 3.599111933878179E+05 IN= 5.244703685174947E+00
 OM= 2.463901296089877E+01 W = 3.204209825950677E+02 Tp=  2460215.546430446208
 N = 1.520234985695220E-04 MA= 3.593901439672831E+02 TA= 3.593016826209310E+02

As PM 2Ring commented, that's relative to the J2000 equinox.
For 2023-09-28, add 0.33° for precession,
making Ω = 24.97° relative to the equinox of date.


---

