<p>Answer provided by @PM2Ring, who linked to <a href="https://iopscience.iop.org/article/10.3847/1538-3881/abd414" rel="nofollow noreferrer">The JPL Planetary and Lunar Ephemerides DE440 and DE441</a> which contained the necessary calculation. The main issue was the was the angles were being rotated with <em>eul2rotm</em>. Using their calculations, I was able to get reasonable values and calculate the sub-solar point, the sub-earth point, the sub-anything really! Thank you again!</p>
<p>Here is the finished code for anyone stumbling across this in the future.</p>
<p>I made it a little more complex to be able to show the progression of the subsolar latitude (important for determining sun conditions for a specific spot)</p>
<pre><code>mission_start = datetime(2022, 1, 1);
mission_end = datetime(2023, 1, 1);
mission_time = linspace(mission_start, mission_end, 100)';

jdt = juliandate(mission_time);
sun_pos = planetEphemeris(jdt, 'Moon', 'Sun');
moon_rot = moonLibration(jdt);

ss_long = zeros(size(jdt));
ss_lat = zeros(size(jdt));
for i = 1:length(jdt)
    phi = moon_rot(i, 1); tet = moon_rot(i, 2); psi = moon_rot(i, 3);
    rotm_z1 = [ cos(-phi)    sin(-phi)    0;
                -sin(-phi)   cos(-phi)    0;
                0           0           1];
    
    rotm_x = [  1   0           0;
                0   cos(-tet)    sin(-tet);
                0   -sin(-tet)   cos(-tet)];
    
    rotm_z2 = [ cos(-psi)    sin(-psi)    0;
                -sin(-psi)   cos(-psi)    0;
                0           0           1];
    
    rotm = rotm_z1 * rotm_x * rotm_z2;
    sun_vec = sun_pos(i, :) / norm(sun_pos(i, :));
    sun_vec = sun_vec * rotm;
    
    [long, lat] = cart2sph(sun_vec(1), sun_vec(2), sun_vec(3));
    ss_long(i) = rad2deg(long); ss_lat(i) = rad2deg(lat);
end

plot(mission_time, ss_lat, '-r')
xtickangle(45)
datetick('x', 'mmm-dd-yy', 'keepticks')
xlabel('Mission Time')
ylabel('Subsolar Latitude [°]')
</code></pre>
<p><a href="https://i.sstatic.net/m65hF.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/m65hF.jpg" alt="Example image of the subsolar latitude over 1 Earth year, showing the near sinusoidal relationship with tiny wobbles." /></a></p>

CC BY-SA 4.0
2022-09-01T14:51:15.400Z
2022-09-01T14:51:15.400Z
Timothy Barrett
<p>So I have been trying to write a program to generate the latitude and longitude of the lunar subsolar point in Matlab.</p>
<p>I have the ephemeris data from aero toolbox, but I can't seem to get values that make sense.</p>
<p>The method I'm currently trying (without success) is...</p>
<ul>
<li>Generate a juliandate for a specific time</li>
<li>Find the position of the sun relative to the moon from the ephemeris data for the specified date</li>
<li>Find the (φ, θ, ψ) lunar attitude from the ephemeris data</li>
<li>Create a rotational matrix from the lunar attitude values</li>
<li>Transpose the matrix to get the inverse rotation (to convert the vector from the ICRF frame to the moon's coordinate system)</li>
<li>Apply that rotation to the direction of the sun vector</li>
<li>Convert to spherical coordinates (longitude, latitude)</li>
</ul>
<pre><code>mission_time = juliandate(2022, 1, 1);
sun_pos = planetEphemeris(mission_time, 'Moon', 'Sun');
moon_rot = moonLibration(mission_time);

rotm = eul2rotm(moon_rot);
rotm = transpose(rotm);

sun_vec = sun_pos / norm(sun_pos);
sun_vec = sun_vec * rotm;

[ss_long, ss_lat] = cart2sph(sun_vec(1), sun_vec(2), sun_vec(3))
fprintf(&quot;Subsolar Lat: %2.2f°\tSubsolar Long: %2.2f°\n&quot;, rad2deg(ss_long), rad2deg(ss_lat))
</code></pre>
<p>The subsolar latitude should be like ±1.57° but this calculation goes all over the place. What am I missing?</p>

CC BY-SA 4.0
2022-09-01T05:19:27.193Z
You can check your calculations against Horizons. Eg, https://ssd.jpl.nasa.gov/api/horizons.api?format=text&MAKE_EPHEM=YES&EPHEM_TYPE=O&OBJ_DATA=NO&COMMAND=301&CENTER=%4010&QUANTITIES=14&START_TIME=2020-1-1&STOP_TIME=2020-2-1&STEP_SIZE=1d Note that the Sun (10) is specified as the center and the Moon (301) is the command target.
PM 2Ring
CC BY-SA 4.0
2022-09-01T12:34:10.337Z
@PM2Ring thanks for the tip! I've gotten the lat, long of the subsolar and sub-earth points and overlayed them (even did a time history to see if they are close) and my code is making a really big sweep, when the real subsolar point stays about ±1.57° about the lunar equator.
Timothy Barrett
CC BY-SA 4.0
2022-09-01T12:36:33.377Z
@CarlWitthoft I've checked the documentation pages and they just say that it returns "The columns contain the Euler angles (φ θ ψ) for Moon attitude, in radians." I'm guessing these angles relate the selenographic coordinate system with the ICRF, but if that's true, I'm not sure where my calculations have gone wrong...

Also long, lat are spherical coordinates? Unless you're one of those flat Earthers haha
Timothy Barrett
CC BY-SA 4.0
2022-09-01T13:05:42.223Z
See page 4 of [The JPL Planetary and Lunar Ephemerides DE440 and DE441](https://iopscience.iop.org/article/10.3847/1538-3881/abd414) for details on the lunar Euler angles & the ICRF3 frame.
PM 2Ring
CC BY-SA 4.0
2022-09-01T13:15:11.560Z
@Carl On a sphere, latitude & longitude are certainly spherical coordinates (although sometimes colatitude is preferred). On an ellipse, there are several useful types of [latitude](https://en.wikipedia.org/wiki/Latitude). To make a flat map, we can treat latitude & longitude as Cartesian coords. From a sphere, that gives us an [equirectangular projection](https://en.wikipedia.org/wiki/Equirectangular_projection).
PM 2Ring
CC BY-SA 4.0
2022-09-01T13:51:36.443Z
@PM2Ring That did it!! Thank you so much! I was using the wrong Euler sequence!! Page 4 was really helpful. I looked at a similar document for DE421, but they didn't go into the details of the calculation like this did. I'll write up a fix of the code as an answer in-case anyone else is interested. Thank you again!
Timothy Barrett
CC BY-SA 4.0
2022-09-01T16:44:36.153Z
@PM2Ring fair enough
Carl Witthoft
<the-moon><software><navigation><lunar-libration>
Lunar Librations and Subsolar Point in Matlab
Timothy Barrett
# Lunar Librations and Subsolar Point in Matlab

(Asked by: username_0 on 2022-08-31)

So I have been trying to write a program to generate the latitude and longitude of the lunar subsolar point in Matlab.
I have the ephemeris data from aero toolbox, but I can't seem to get values that make sense.
The method I'm currently trying (without success) is...

Generate a juliandate for a specific time
Find the position of the sun relative to the moon from the ephemeris data for the specified date
Find the (φ, θ, ψ) lunar attitude from the ephemeris data
Create a rotational matrix from the lunar attitude values
Transpose the matrix to get the inverse rotation (to convert the vector from the ICRF frame to the moon's coordinate system)
Apply that rotation to the direction of the sun vector
Convert to spherical coordinates (longitude, latitude)

mission_time = juliandate(2022, 1, 1);
sun_pos = planetEphemeris(mission_time, 'Moon', 'Sun');
moon_rot = moonLibration(mission_time);

rotm = eul2rotm(moon_rot);
rotm = transpose(rotm);

sun_vec = sun_pos / norm(sun_pos);
sun_vec = sun_vec * rotm;

[ss_long, ss_lat] = cart2sph(sun_vec(1), sun_vec(2), sun_vec(3))
fprintf("Subsolar Lat: %2.2f°\tSubsolar Long: %2.2f°\n", rad2deg(ss_long), rad2deg(ss_lat))

The subsolar latitude should be like ±1.57° but this calculation goes all over the place. What am I missing?

**Tags:** the-moon, software, navigation, lunar-libration
> **Comments:**
> **username_1**: You can check your calculations against Horizons. Eg, https://ssd.jpl.nasa.gov/api/horizons.api?format=text&MAKE_EPHEM=YES&EPHEM_TYPE=O&OBJ_DATA=NO&COMMAND=301&CENTER=%4010&QUANTITIES=14&START_TIME=2020-1-1&STOP_TIME=2020-2-1&STEP_SIZE=1d Note that the Sun (10) is specified as the center and the Moon (301) is the command target.  
> **username_0**: @PM2Ring thanks for the tip! I've gotten the lat, long of the subsolar and sub-earth points and overlayed them (even did a time history to see if they are close) and my code is making a really big sweep, when the real subsolar point stays about ±1.57° about the lunar equator.  
> **username_0**: @CarlWitthoft I've checked the documentation pages and they just say that it returns "The columns contain the Euler angles (φ θ ψ) for Moon attitude, in radians." I'm guessing these angles relate the selenographic coordinate system with the ICRF, but if that's true, I'm not sure where my calculations have gone wrong...  Also long, lat are spherical coordinates? Unless you're one of those flat Earthers haha  
> **username_1**: See page 4 of [The JPL Planetary and Lunar Ephemerides DE440 and DE441](https://iopscience.iop.org/article/10.3847/1538-3881/abd414) for details on the lunar Euler angles & the ICRF3 frame.  
> **username_1**: @Carl On a sphere, latitude & longitude are certainly spherical coordinates (although sometimes colatitude is preferred). On an ellipse, there are several useful types of [latitude](https://en.wikipedia.org/wiki/Latitude). To make a flat map, we can treat latitude & longitude as Cartesian coords. From a sphere, that gives us an [equirectangular projection](https://en.wikipedia.org/wiki/Equirectangular_projection).  
> **username_0**: @PM2Ring That did it!! Thank you so much! I was using the wrong Euler sequence!! Page 4 was really helpful. I looked at a similar document for DE421, but they didn't go into the details of the calculation like this did. I'll write up a fix of the code as an answer in-case anyone else is interested. Thank you again!  
> **username_2**: @PM2Ring fair enough  

---

## 1 answers

### Answer by username_0 on 2022-09-01. Score: 3

Answer provided by @PM2Ring, who linked to The JPL Planetary and Lunar Ephemerides DE440 and DE441 which contained the necessary calculation. The main issue was the was the angles were being rotated with eul2rotm. Using their calculations, I was able to get reasonable values and calculate the sub-solar point, the sub-earth point, the sub-anything really! Thank you again!
Here is the finished code for anyone stumbling across this in the future.
I made it a little more complex to be able to show the progression of the subsolar latitude (important for determining sun conditions for a specific spot)
mission_start = datetime(2022, 1, 1);
mission_end = datetime(2023, 1, 1);
mission_time = linspace(mission_start, mission_end, 100)';

jdt = juliandate(mission_time);
sun_pos = planetEphemeris(jdt, 'Moon', 'Sun');
moon_rot = moonLibration(jdt);

ss_long = zeros(size(jdt));
ss_lat = zeros(size(jdt));
for i = 1:length(jdt)
    phi = moon_rot(i, 1); tet = moon_rot(i, 2); psi = moon_rot(i, 3);
    rotm_z1 = [ cos(-phi)    sin(-phi)    0;
                -sin(-phi)   cos(-phi)    0;
                0           0           1];
    
    rotm_x = [  1   0           0;
                0   cos(-tet)    sin(-tet);
                0   -sin(-tet)   cos(-tet)];
    
    rotm_z2 = [ cos(-psi)    sin(-psi)    0;
                -sin(-psi)   cos(-psi)    0;
                0           0           1];
    
    rotm = rotm_z1 * rotm_x * rotm_z2;
    sun_vec = sun_pos(i, :) / norm(sun_pos(i, :));
    sun_vec = sun_vec * rotm;
    
    [long, lat] = cart2sph(sun_vec(1), sun_vec(2), sun_vec(3));
    ss_long(i) = rad2deg(long); ss_lat(i) = rad2deg(lat);
end

plot(mission_time, ss_lat, '-r')
xtickangle(45)
datetick('x', 'mmm-dd-yy', 'keepticks')
xlabel('Mission Time')
ylabel('Subsolar Latitude [°]')


---

