# Lunar Librations and Subsolar Point in Matlab

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))


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

• You can check your calculations against Horizons. Eg, ssd.jpl.nasa.gov/api/… Note that the Sun (10) is specified as the center and the Moon (301) is the command target. Commented Sep 1, 2022 at 5:19
• @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. Commented Sep 1, 2022 at 12:34
• @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 Commented Sep 1, 2022 at 12:36
• See page 4 of The JPL Planetary and Lunar Ephemerides DE440 and DE441 for details on the lunar Euler angles & the ICRF3 frame. Commented Sep 1, 2022 at 13:05
• @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. To make a flat map, we can treat latitude & longitude as Cartesian coords. From a sphere, that gives us an equirectangular projection. Commented Sep 1, 2022 at 13:15

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));