Access to ephemeris of celestial bodies

Spacecraft, Thruster, Tank, Propagator, CoordinateSystem, etc.

Access to ephemeris of celestial bodies

Postby prince » Thu Jul 12, 2012 11:57 pm

Is it possible to access the orbital parameters (cartesian or keplerian) of solar system bodies (Sun, Earth, Luna, etc) in the GMAT script just like we can do for a spacecraft (like DefaultSC.SMA)?
prince
 
Posts: 16
Joined: Mon Jul 02, 2012 8:02 pm

Re: Access to ephemeris of celestial bodies

Postby jjkparker » Fri Jul 13, 2012 3:51 pm

I don't believe this is possible. You can report parameters, and these are all listed when you click the Edit button on the ReportFile properties window (after adding a ReportFile to the Output folder in the GUI).

One option here is to write a MATLAB script to use SPICE to read the same DE file that GMAT is using, then pass back the state to GMAT. This forum has some extensive topics describing how to use the MATLAB interface (just search for them).
Joel J. K. Parker
Flight dynamics engineer, GMAT team
http://gmatcentral.org
jjkparker
 
Posts: 617
Joined: Thu Jan 07, 2010 9:48 pm
Location: NASA Goddard Space Flight Center, Greenbelt, MD

Re: Access to ephemeris of celestial bodies

Postby shughes » Mon Jul 16, 2012 11:47 am

There is a way to do this, although it is not as clean as I would like. We're working on defining features for R2013b and the ability to get access to lower level data, in a more elegant manner, is currently on the list.

Say you want body1 w/r/t to body2. Put a spacecraft at the center of a coordinate system that has body1 as the origin (i.e. coordinate [0,0,0]), and ask for the spacecraft state in a coordinate system centered at body2. Assume you need earth w/r/t sun:

Code: Select all
Create Spacecraft earthSat
earthSat.Epoch.UTCGregorian = '04 Jul 2012 12:00:00.000'
earthSat.CoordinateSystem   = EarthMJ2000Eq;
earthSat.X = 0;
earthSat.Y = 0;
earthSat.Z = 0;


Create CoordinateSystem SunMJ2000Eq
SunMJ2000Eq.Axes    = MJ2000Eq
SunMJ2000Eq.Origin  = Sun;

Create ReportFile rf

BeginMissionSequence

Report rf earthSat.SunMJ2000Eq.X earthSat.SunMJ2000Eq.Y earthSat.SunMJ2000Eq.Z
shughes
 
Posts: 443
Joined: Mon Jun 09, 2008 6:27 pm

Re: Access to ephemeris of celestial bodies

Postby prince » Mon Jul 16, 2012 6:29 pm

thanks that's a nice trick! :)
prince
 
Posts: 16
Joined: Mon Jul 02, 2012 8:02 pm

Re: Access to ephemeris of celestial bodies

Postby prince » Sun Aug 05, 2012 3:07 am

shughes wrote:There is a way to do this, although it is not as clean as I would like. We're working on defining features for R2013b and the ability to get access to lower level data, in a more elegant manner, is currently on the list.

Say you want body1 w/r/t to body2. Put a spacecraft at the center of a coordinate system that has body1 as the origin (i.e. coordinate [0,0,0]), and ask for the spacecraft state in a coordinate system centered at body2. Assume you need earth w/r/t sun:

Code: Select all
Create Spacecraft earthSat
earthSat.Epoch.UTCGregorian = '04 Jul 2012 12:00:00.000'
earthSat.CoordinateSystem   = EarthMJ2000Eq;
earthSat.X = 0;
earthSat.Y = 0;
earthSat.Z = 0;


Create CoordinateSystem SunMJ2000Eq
SunMJ2000Eq.Axes    = MJ2000Eq
SunMJ2000Eq.Origin  = Sun;

Create ReportFile rf

BeginMissionSequence

Report rf earthSat.SunMJ2000Eq.X earthSat.SunMJ2000Eq.Y earthSat.SunMJ2000Eq.Z



I tried but it seems like using this way the state of earth is not propagated with time. See the following script
Code: Select all

Create Spacecraft earthSat
earthSat.Epoch.UTCGregorian = '04 Jul 2012 12:00:00.000'
earthSat.CoordinateSystem   = EarthMJ2000Eq;
earthSat.X = 0;
earthSat.Y = 0;
earthSat.Z = 0;

Create Spacecraft earthSat1
earthSat1.Epoch.UTCGregorian = '04 Jul 2012 12:00:00.000'
earthSat1.CoordinateSystem   = EarthMJ2000Eq;
earthSat1.X = 6678;
earthSat1.Y = 0;
earthSat1.Z = 0;

Create String newrow;
GMAT newrow = ' ';

Create ForceModel DefaultProp_ForceModel;
GMAT DefaultProp_ForceModel.CentralBody = Earth;
GMAT DefaultProp_ForceModel.PrimaryBodies = {Earth};
GMAT DefaultProp_ForceModel.Drag = None;
GMAT DefaultProp_ForceModel.SRP = Off;
GMAT DefaultProp_ForceModel.RelativisticCorrection = Off;
GMAT DefaultProp_ForceModel.ErrorControl = RSSStep;
GMAT DefaultProp_ForceModel.GravityField.Earth.Degree = 4;
GMAT DefaultProp_ForceModel.GravityField.Earth.Order = 4;
GMAT DefaultProp_ForceModel.GravityField.Earth.PotentialFile = 'JGM2.cof';
GMAT DefaultProp_ForceModel.GravityField.Earth.EarthTideModel = 'None';

Create Propagator DefaultProp;
GMAT DefaultProp.FM = DefaultProp_ForceModel;
GMAT DefaultProp.Type = RungeKutta89;
GMAT DefaultProp.InitialStepSize = 60;
GMAT DefaultProp.Accuracy = 9.999999999999999e-012;
GMAT DefaultProp.MinStep = 0.001;
GMAT DefaultProp.MaxStep = 2700;
GMAT DefaultProp.MaxStepAttempts = 50;
GMAT DefaultProp.StopIfAccuracyIsViolated = true;


Create CoordinateSystem SunMJ2000Eq
SunMJ2000Eq.Axes    = MJ2000Eq
SunMJ2000Eq.Origin  = Sun;

Create ReportFile rf

BeginMissionSequence

Report rf earthSat.SunMJ2000Eq.X earthSat.SunMJ2000Eq.Y earthSat.SunMJ2000Eq.Z

Report rf earthSat1.SunMJ2000Eq.X earthSat1.SunMJ2000Eq.Y earthSat1.SunMJ2000Eq.Z
Report rf earthSat1.EarthMJ2000Eq.X earthSat1.EarthMJ2000Eq.Y earthSat1.EarthMJ2000Eq.Z

Propagate DefaultProp(earthSat1) {earthSat1.ElapsedSecs = 120000.0};

Report rf newrow

Report rf earthSat.SunMJ2000Eq.X earthSat.SunMJ2000Eq.Y earthSat.SunMJ2000Eq.Z

Report rf earthSat1.SunMJ2000Eq.X earthSat1.SunMJ2000Eq.Y earthSat1.SunMJ2000Eq.Z
Report rf earthSat1.EarthMJ2000Eq.X earthSat1.EarthMJ2000Eq.Y earthSat1.EarthMJ2000Eq.Z


Results
Code: Select all

earthSat.SunMJ2000Eq.X   earthSat.SunMJ2000Eq.Y   earthSat.SunMJ2000Eq.Z   
33526136.91647737        -136111176.19852         -59006945.77960628       
earthSat1.SunMJ2000Eq.X   earthSat1.SunMJ2000Eq.Y   earthSat1.SunMJ2000Eq.Z   
33532814.91647737         -136111176.19852          -59006945.77960628       
earthSat1.EarthMJ2000Eq.X   earthSat1.EarthMJ2000Eq.Y   earthSat1.EarthMJ2000Eq.Z   
6678                        0                           0                           
newrow                 
                       
earthSat.SunMJ2000Eq.X   earthSat.SunMJ2000Eq.Y   earthSat.SunMJ2000Eq.Z   
33526136.91647737        -136111176.19852         -59006945.77960628       
earthSat1.SunMJ2000Eq.X   earthSat1.SunMJ2000Eq.Y   earthSat1.SunMJ2000Eq.Z   
36947213.61357012         -135369604.859301         -58683818.51481943       
earthSat1.EarthMJ2000Eq.X   earthSat1.EarthMJ2000Eq.Y   earthSat1.EarthMJ2000Eq.Z   
3270.819287809069           -5335.765853022342          -577.4975532559364         
prince
 
Posts: 16
Joined: Mon Jul 02, 2012 8:02 pm

Re: Access to ephemeris of celestial bodies

Postby jjkparker » Mon Aug 06, 2012 6:09 pm

This is because time isn't universal to the script; it's local to each Spacecraft.

At the beginning of the script you have earthSat and earthSat1, both of which are at 04 Jul 2012 12:00:00.000. Then you're propagating earthSat1 forward by 120000 seconds, so now it's at 05 Jul 2012 21:20:00.000. But when you run this line the second time:
Code: Select all
Report rf earthSat.SunMJ2000Eq.X earthSat.SunMJ2000Eq.Y earthSat.SunMJ2000Eq.Z

earthSat is still at 04 Jul 2012 12:00:00.000.

Try something like this (mission sequence only). Notice the new line after Report rf newrow:
Code: Select all
Report rf earthSat.SunMJ2000Eq.X earthSat.SunMJ2000Eq.Y earthSat.SunMJ2000Eq.Z

Report rf earthSat1.SunMJ2000Eq.X earthSat1.SunMJ2000Eq.Y earthSat1.SunMJ2000Eq.Z
Report rf earthSat1.EarthMJ2000Eq.X earthSat1.EarthMJ2000Eq.Y earthSat1.EarthMJ2000Eq.Z

Propagate DefaultProp(earthSat1) {earthSat1.ElapsedSecs = 120000.0};

Report rf newrow
earthSat.Epoch = earthSat1.UTCGregorian
Report rf earthSat.SunMJ2000Eq.X earthSat.SunMJ2000Eq.Y earthSat.SunMJ2000Eq.Z

Report rf earthSat1.SunMJ2000Eq.X earthSat1.SunMJ2000Eq.Y earthSat1.SunMJ2000Eq.Z
Report rf earthSat1.EarthMJ2000Eq.X earthSat1.EarthMJ2000Eq.Y earthSat1.EarthMJ2000Eq.Z
Joel J. K. Parker
Flight dynamics engineer, GMAT team
http://gmatcentral.org
jjkparker
 
Posts: 617
Joined: Thu Jan 07, 2010 9:48 pm
Location: NASA Goddard Space Flight Center, Greenbelt, MD

Re: Access to ephemeris of celestial bodies

Postby prince » Tue Aug 07, 2012 5:00 am

thanks it works.

Also I saw a minor issue with STM propagation that if we propagate the same states with and without STM, the results are not exactly same. I have attached the script. The error is of the order of 1e-4 when I propagated for 120000 secs.
Attachments
test.script
(6.7 KiB) Downloaded 364 times
prince
 
Posts: 16
Joined: Mon Jul 02, 2012 8:02 pm

Re: Access to ephemeris of celestial bodies

Postby shughes » Thu Aug 09, 2012 4:46 pm

The difference is because the adaptive step algorithm has 36 additional states that can affect the step size and so the integrator is taking slightly different steps and arriving at a slightly different answer that still meets the integration tolerance of 1e-11. I tried two things. (1) I reduced the Acurracy setting to 1e-14 and get sub millimeter agreement as shown below. I set to fixed step (DefaultProp_ForceModel.ErrorControl = None;) and see identical results.

Code: Select all
earthSat1.SunMJ2000Eq.X   earthSat1.SunMJ2000Eq.Y   earthSat1.SunMJ2000Eq.Z   earthSat1.SunMJ2000Eq.VX 
36947213.61359383         -135369604.8592888        -58683818.51481693        35.44631554785417          10.04496492114358   
earthSat1.EarthMJ2000Eq.X   earthSat1.EarthMJ2000Eq.Y   earthSat1.EarthMJ2000Eq.Z   earthSat1.EarthMJ2000Eq.VX     
3270.819311521271           -5335.765840896569          -577.4975507577153          7.044415173766448 
earthSat2.SunMJ2000Eq.X   earthSat2.SunMJ2000Eq.Y   earthSat2.SunMJ2000Eq.Z   earthSat2.SunMJ2000Eq.VX   
36947213.61359354         -135369604.859289         -58683818.51481707        35.44631554807062          10.0449649207858
earthSat2.EarthMJ2000Eq.X   earthSat2.EarthMJ2000Eq.Y   earthSat2.EarthMJ2000Eq.Z   earthSat2.EarthMJ2000Eq.VX   
3270.819311230655           -5335.765841041648          -577.4975508946158          7.044415173982896                     
shughes
 
Posts: 443
Joined: Mon Jun 09, 2008 6:27 pm

Re: Access to ephemeris of celestial bodies

Postby prince » Fri Sep 07, 2012 7:43 pm

that makes sense. However I am seeing one more behavior that I am not able to understand fully. I propagate 2 satellites using exactly same propagators. The satellites initial conditions however are not exactly same. After the propagation, I assign state (X,Y,Z,VX,Vy,VZ) of 1 satellite to another and immediately subtract the states of 2 satellites in a different coordinate system. The difference is not zero. If I subtract in the same coordinate system then difference is 0 as expected. I am not sure what I am doing wrong. I have attached my test script.
Attachments
test.script
(9.66 KiB) Downloaded 375 times
prince
 
Posts: 16
Joined: Mon Jul 02, 2012 8:02 pm

Re: Access to ephemeris of celestial bodies

Postby shughes » Fri Sep 07, 2012 9:54 pm

Here are my first thoughts, but I am not convinced I fully understand this issue yet. For starters, if you use a single propagate command like the code snippet below, which ensures the epochs of the spacecraft are precisely synchronized, then this issue goes away:

Code: Select all
Propagate DefaultProp(DefaultSC,DefaultSC3) {DefaultSC.ElapsedSecs = 60};


Here are some thoughts on the behavior for separate propagations: GMAT models time using a double precision real number in a modified Julian date format (less than ideal and we want to change that eventually). I can't prove it yet, but since the synchronized propagation avoids this issue, I believe the epochs of the separately propagated spacecraft are very slightly different due to ~5 microsecond time discretization that occurs when modelling time as a modified Julian date using a double precision real. When you take into account that the Earth moves about 30 km/s and error in time is about 5 microseconds you get 0.5e-5 s * 30 km/s = 1.5 cm. The RSS of difference of the actual data shown below is 1.1 cm. The precision of a double may come into play here in more ways though because for these values there are 9 digits to the left of the decimal leaving about 6 useful digits to the right which is what we see in the difference.

Code: Select all
DefaultSC.TAIModJulian    DefaultSC.SunMJ2000Ec.X    DefaultSC.SunMJ2000Ec.Y    DefaultSC.SunMJ2000Ec.Z   
21545.00069444444        -26494694.3116229           144697736.8044981          457.9547021538019     
           
DefaultSC3.TAIModJulian   DefaultSC3.SunMJ2000Ec.X   DefaultSC3.SunMJ2000Ec.Y   DefaultSC3.SunMJ2000Ec.Z 
21545.00069444444         -26494694.31161353         144697736.8044998          457.9547021538019         
     
X2(1,1)                   X2(2,1)                    X2(3,1)                    X2(4,1)                             
-9.365379810333252e-06   -1.728534698486328e-06      0 
shughes
 
Posts: 443
Joined: Mon Jun 09, 2008 6:27 pm

Next

Return to GMAT Resources

Who is online

Users browsing this forum: No registered users and 1 guest