The JPL planetary ephemerides are saved as files of Chebyshev polynomials fit to the Cartesian positions and velocities of the planets, Sun, and Moon, typically in 32-day intervals. The positions are integrated in astronomical units (au), but with polynomials stored in units of kilometers.
Folkner et al., JPL Interplanetary Network Progress Report 42-196, 2014.
Meysam Mahooti (2021). NASA JPL Development Ephemerides (DE430) (https://www.mathworks.com/matlabcentral/fileexchange/60504-nasa-jpl-development-ephemerides-de430), MATLAB Central File Exchange. Retrieved .
The input parameter should be the Modified julian date of TDB, however, in the test, the input parameter is the Modified julian date of UTC. The example will mislead the user.
Thank you so much for your sharing! It helps me a lot. I compared the result of it with that of the STK software, the bias is extremely small. If it adopt the coordinate transformation function，it will be perfect.
Thanks for the coding. I can confirm the input should be in TDB not UTC as in the test file. At least, I checked it with Sun ephemeris. The error using UTC instead of TDB is about 3arcsec
Hi, could you add some information about the frame it adopts in output? I compare the result of this euphemisms and JPL HORIZONS system, The results in SSB are quite different. While the relative vectors between planets and sun are almost the same. I guess it is because of the coordinate.
Hi, many thanks for this kind sharing. While, the code might need some refinement to output the velocity of planets.
Found the problem. User error. When I built the ephemeris generator I stuck Julian Day in the first column. Then when I plotted, I used the first three columns.
Downloaded and ran. Worked fine. Modified test script to create a year-long ephemeris of earth-to-sun vectors. Hoped you could help me understand the output reference frame.
I made a simple 3D plot as a check to verify it was working right. Found the axes as output don't correspond to my understanding of the equatorial celestial reference frame. The Z-axis ought to correspond to the direction of the north pole which would mean it would be the smallest. In fact, as output, the x-axis is the small one. Could you please describe this reference frame?
Also found the axis with the small value was decreasing a few miles per year. I'm assuming this is precession? Cool!
It seems a discrepancy: function in TDB while test code calls Mjd_UTC(),
which means UTC time. Could someone clarify? Thx in advance.
% JPL_Eph_DE430: Computes the sun, moon, and nine major planets' equatorial
% position using JPL Ephemerides
% Mjd_TDB Modified julian date of TDB
% Last modified: 2018/01/11 M. Mahooti
% test code
format long g
PC = DE430Coeff;
Mjd_UTC = Mjday(2017,8,21,17,28,11.284);
PC is a matrix of 1020 columns. The first two double precision words contain Julian date of earliest data in record and Julian date of latest data in record. The remainder 1018 columns include 15 triplets corresponding to the following items:
Earth Nutations in longitude and obliquity (IAU 1980 model)
Lunar mantle libration
Lunar mantle angular velocity
TT-TDB (at geocenter)
3 171 231 309 342 366 387 405 423 441 753 819 899
14 10 13 11 8 7 6 6 6 13 11 10 10
4 2 2 1 1 1 1 1 1 8 2 4 4
Starting location of the Chebyshev coefficients for Mercury is 3rd column in each data record (1020). The number of Chebyshev coefficients for every component of the Mercury is 14 and there are 4 complete sets of coefficients in each data record. These information is used to produce following codes in JPL_Eph_DE430.m:
temp = (3:14:45);
Cx_Mercury = PCtemp(temp(1):temp(2)-1);
Cy_Mercury = PCtemp(temp(2):temp(3)-1);
Cz_Mercury = PCtemp(temp(3):temp(4)-1);
temp = temp+42;
Cx = PCtemp(temp(1):temp(2)-1);
Cy = PCtemp(temp(2):temp(3)-1);
Cz = PCtemp(temp(3):temp(4)-1);
Cx_Mercury = [Cx_Mercury,Cx];
Cy_Mercury = [Cy_Mercury,Cy];
Cz_Mercury = [Cz_Mercury,Cz];
if (0<=dt && dt<=8)
Mjd0 = t1;
elseif(8<dt && dt<=16)
Mjd0 = t1+8*j;
elseif (16<dt && dt<=24)
Mjd0 = t1+8*j;
elseif(24<dt && dt<=32)
Mjd0 = t1+8*j;
r_Mercury = 1e3*Cheb3D(Mjd_TDB, 14, Mjd0, Mjd0+8, Cx_Mercury(14*j+1:14*j+14),...
For more information please look at description sub-folder.
Thanks for providing this source code.
In AST_Const.m, the line 'GM_Pluto = 977.0000009; % [m^3/s^2]; DE430' is supposed to be 'GM_Pluto = 977.00000e9; % [m^3/s^2]; DE430'.
In addition, it is more efficient to set PC in JPL_Eph_DE430.m as a persistent and use the following codes:
i = find(PC(:,1)<=JD & JD<=PC(:,2),1,'first');
PCtemp = PC(i,:);
to replace the block for higher efficiency:
if (PC(i,1)<=JD && JD<=PC(i,2))
PCtemp = PC(i,:);
By the way, would you like to provide more details and structure of PC so that I could do some modification to boost the efficiency (mainly removing the for loops)? Thanks a lot.
Find the treasures in MATLAB Central and discover how the community can help you!Start Hunting!