Computing Sun position from Mars

# Computing Sun position from Mars

To compute planetary ephemerides from a given point located on the surface of another planet is something initially beyond the scope of the JPARSEC library. However, this question has recently attracted my attention since it is something interesting in certain projects like the Mars exploration rovers. Spirit and Opportunity were a real success, so the design of new rovers to explore Mars with better technology is something that can easy earn the attention (and the money) from USA goverment. The new rovers will be Curiosity and Phoenix, and are now in the final phase of construction. Curiosity will be launched in december, 2011, and will reach Mars on August, 2012.

In terms of software development, the most visible program produced by JPL during the Spirit and Opportunity adventure was Maestro. This software was designed to schedule the activities of the rovers and to analyze the data taken with their instruments, to obtain the maximum science from them. A public version of Maestro allowed everyone to download very interesting images and data just after they were taken and analyzed by scientists.

Mark W. Powell is one of the engineers working at JPL in the development team of the next version of Maestro. He recently contacted me asking for a program to obtain the ephemerides of the Sun from any given point on Mars. I always thought there should be a number of experts in ephemerides calculations at JPL (that's sure), but it seems that the different tools they have for ephemerides (let's say Horizons, NOVAS, among others) are not adequate since Maestro is developed in Java. I never thought to collaborate with JPL in the field of ephemerides calculation, but the characteristics of the JPARSEC library (open source, well documented and organized in Mark's words, combined with some lack of good Java software in this field) are appreciated in new projects.

The point in this case is to calculate if a given activity (like to capture an image of a given rock from a given direction and at a given time) is going to be made in shadow or not. The new rovers will be capable or shooting images at night and taking some decisions by themselves, but the presence of the Sun changes the lighting conditions and therefore the steps to be taken by the science investigation team, or the schedule to give priority to some scientific proposals. Once the position of the Sun is known, the lighting conditions on Mars can be simulated using 3d models for both the rovers and the terrain, giving a direct look of the shadows on the terrain, coming both from the rovers or from other rocks.

So I helped Mark (of course) in this interesting real case, providing him with a program that uses a reduced version of the JPARSEC library. The position of an object on Mars can be calculated in several ways. The most elaborated one is to use Mars as reference instead of the Earth to calculate the 'geocentric' coordinates of the objects, and then to use the Mars geoid for 'topocentric' positions. I discarded this since I realized I needed more than a few days to do it properly, with all necessary tests and developments to give him a clean code with updated javadoc. The second approach was to obtain the position of the Sun using the subsolar position from Earth that JPARSEC already calculates with excellent precision (using geometric coordinates instead of apparent ones), and doing a simple rotation. The third method was to implement the algorithm described in the program Mars24, done by a team at Goddard Space Flight Center (see http://www.giss.nasa.gov/tools/mars24/help/algorithm.html for a sketch of the algorithm, and http://pubs.giss.nasa.gov/docs/2000/2000_Allison_McEwen.pdf for the details). As reference for accuracy tests I used Horizons. After a complete check I concluded that the mean difference between JPARSEC and Horizons was in the 0.02 to 0.05 degrees range, and 0.15 to 0.25 relative to Mars24. Half of the difference with JPARSEC seems to be due to a different correction from UTC to TT (or maybe Horizons uses UT1?), and the other half to that simple rotation.

Both JPARSEC and Mars24 methods were suitable for Mark, but the better accuracy of JPARSEC combined with the organization, documentation, and simplicity of the distribution as a .jar file convinced him to use it for Maestro. I'm quite confident that the quality of my work is good enough, but now I feel somewhat responsible for such a large project involving so many engineers and scientists. Hopefully we will have a public version of Maestro even better than the one released for Spirit and Opportunity to allow people interested in this great project to participate in the exploration of the Mars surface.

Here is the code for the MarsEphem class I produced for Mark. The class implements a test over 5000 random locations on Mars and dates between 2010 and 2020 and returns maximum difference between both algorithms.

package jparsec.ephem;

import jparsec.ephem.planets.EphemElement;
import jparsec.math.Constant;
import jparsec.observer.LocationElement;
import jparsec.observer.ObserverElement;
import jparsec.time.AstroDate;
import jparsec.time.TimeElement;
import jparsec.time.TimeFormat;
import jparsec.time.TimeScale;

/**
* Ephemerides for the Sun from Mars.
* @author T. Alonso Albi - OAN (Spain)
* @version 1
*/
public class MarsEphem {

/**
* Test program.
* @param args Not used.
*/
public static void main(String args[]) {

try {
double maxDist = -1, maxDist_TT = 0.0;

for (int iii=0; iii< 5000; iii++) {

int year = 2010 + (int) (Math.random() * 10);
int month = 1 + (int) (Math.random() * 11);
int day = 1 + (int) (Math.random() * 27);
int hour = 0 + (int) (Math.random() * 23);
int minute = 0 + (int) (Math.random() * 59);
int second = 0 + (int) (Math.random() * 59);
double lon = Math.random() * 360, lat = 90 - Math.random() * 180;

// Example 1: 2000 Jan 6 00:00:00 UTC, from lon = lat = 0
// Example 2: 2004 Jan 3 13:46:31 UTC, from lon = 184.702, lat = -14.64
// Example 3: 2010 Jan 6 00:00:00 UTC, from lon = 184.702, lat = -14.64
/*
year = 2000;
month = 1;
day = 6;
hour = 0;
minute = 0;
second = 0;
lon = 0;
lat = 0;
*/
AstroDate astro = new AstroDate(year, month, day, hour, minute, second); // Example 1
TimeElement time = new TimeElement(astro.toGCalendar(), TimeElement.TERRESTRIAL_TIME);
EphemerisElement eph = new EphemerisElement(Target.MARS, EphemerisElement.EPHEM_GEOMETRIC,
EphemerisElement.EQUINOX_OF_DATE, EphemerisElement.GEOCENTRIC, EphemerisElement.APPLY_IAU2000,
EphemerisElement.FRAME_J2000, EphemerisElement.ALGORITHM_MOSHIER);
EphemElement ephem = Ephem.getEphemeris(time, observer, eph, false);

System.out.println("USING JPARSEC");

// Obtain Sun position in Mars knowing subsolar position
double dif = ephem.subsolarLatitude - (float) Math.atan(Math.tan(ephem.subsolarLatitude) / Math.pow(3376.2 / 3396.19, 2.0));
LocationElement subSolar = new LocationElement(ephem.subsolarLongitude, ephem.subsolarLatitude + dif, 1.0);
LocationElement refPos = new LocationElement(observer.longitude, observer.latitude, 1.0);
double elev = Math.PI * 0.5 - LocationElement.getAngularDistance(subSolar, refPos);
System.out.println("Sun azimuth     " + pa * Constant.RAD_TO_DEG);
System.out.println("Sun elevation   " + elev * Constant.RAD_TO_DEG);

// Obtain Earth position in Mars knowing sub-Earth position
/*				dif = ephem.positionAngleOfPole - (float) Math.atan(Math.tan(ephem.positionAngleOfPole) / Math.pow(3376.2 / 3396.19, 2.0));
LocationElement subEarth = new LocationElement(ephem.longitudeOfCentralMeridian, ephem.positionAngleOfPole + dif, 1.0);
elev = Math.PI * 0.5 - LocationElement.getAngularDistance(subEarth, refPos);
pa = LocationElement.getPositionAngle(refPos, subEarth);
*/
// Obtain Sun position in Mars using simple method by Goddard Space Flight Center
// See http://www.giss.nasa.gov/tools/mars24/help/algorithm.html and http://pubs.giss.nasa.gov/docs/2000/2000_Allison_McEwen.pdf
System.out.println("USING MARS24");
double JD_TT = TimeScale.getJD(time, observer, eph, TimeScale.JD_TT);
double dt = JD_TT - 2451545.0;
double M = Math.toRadians(19.387 + 0.52402075 * dt);
double afms = Math.toRadians(270.3863 + 0.5240384 * dt);
double pbs = 0.0;
double a[] = new double[] {0.0071, 0.0057, 0.0039, 0.0037, 0.0021, 0.0020, 0.0018};
double tau[] = new double[] {2.2353, 2.7543, 1.1177, 15.7866, 2.1354, 2.4694, 32.8493};
double phi[] = new double[] {49.409, 168.173, 191.837, 21.736, 15.704, 95.528, 49.095};
for (int i=0; i<7; i++) {
pbs += a[i] * Math.cos(Math.toRadians(phi[i] + 0.985626 * dt / tau[i]));
}
double vminusM = Math.toRadians((10.691 + 3.0E-7*dt) * Math.sin(M) + 0.623 * Math.sin(2.0 * M) +
0.050 * Math.sin(3.0 * M) + 0.005 * Math.sin(4.0 * M) + 0.0005 * Math.sin(5.0 * M) + pbs);
double ls = afms + vminusM;
double eot = Math.toRadians(2.861 * Math.sin(2.0 * ls) - 0.071 * Math.sin(4.0 * ls) + 0.002 * Math.sin(6.0 * ls)) - vminusM;
double mtc = 24.0 * ((JD_TT - 2451549.5) / 1.027491252 + 44796.0 - 0.00096);
mtc = (mtc - (int) (mtc / 24.0) * 24.0) * (2.0 * Math.PI) / 24.0;
double deltas = Math.asin(0.42565 * Math.sin(ls)) + Math.toRadians(0.25 * Math.sin(ls));
double lons = Functions.normalizeRadians((mtc + eot) + Math.PI);
double H = observer.longitude - lons;
double z = Math.acos(Math.sin(deltas) * Math.sin(observer.latitude) + Math.cos(deltas) * Math.cos(observer.latitude) * Math.cos(H));
// Azimuth = angle from north direction, measured clockwise
double azimuth = Functions.normalizeRadians(Math.atan2(Math.sin(H), (Math.cos(observer.latitude) * Math.tan(deltas) - Math.sin(observer.latitude) * Math.cos(H))));
double elevation = Math.PI * 0.5 - z;

// COMPARISON TEST
//                       JPL HORIZONS             JPARSEC            MARS24
// 2000-Jan-06 00:00     191.1564 -64.5049    191.1722  -64.5053     191.0398  -64.2616
// 2004-Jan-03 13:46:31  179.9952 -62.0741    180.0531  -62.0770     179.9890  -61.9392
// 2010-Jan-06 00:00     284.5225   1.1738    284.5392    1.2354     284.6525    1.1303
// RMS RELATIVE TO JPL HORIZONS:                     0.05                   0.15
// RMS RELATIVE TO JPL HORIZONS (with TT, no UT)     0.02                   0.15

/*
System.out.println(JD_TT);
System.out.println(dt);
System.out.println(Math.toDegrees(M));
System.out.println(Math.toDegrees(afms));
System.out.println("pbs "+pbs);
System.out.println(Math.toDegrees(vminusM));
System.out.println(Math.toDegrees(ls));
System.out.println(eot*24.0/(2.0*Math.PI));
System.out.println("mtc "+mtc*24.0/(2.0*Math.PI));
System.out.println(Math.toDegrees(lons)/15.0);
System.out.println(Math.toDegrees(deltas));
*/				System.out.println("Sun azimuth     " + Math.toDegrees(azimuth));
System.out.println("Sun elevation   " + Math.toDegrees(elevation));

LocationElement loc1 = new LocationElement(pa, elev, 1);
LocationElement loc2 = new LocationElement(azimuth, elevation, 1);
double dist = LocationElement.getAngularDistance(loc1, loc2) * Constant.RAD_TO_DEG;
if (dist > maxDist || maxDist == -1) {
maxDist = dist;
maxDist_TT = JD_TT;
}
}
System.out.println(maxDist);
System.out.println(maxDist_TT);
System.out.println(TimeFormat.formatJulianDayAsDateAndTime(maxDist_TT));
} catch (Exception exc) {
exc.printStackTrace();
}
}
}

NOTE: this is al old post written before the implementation of some improvements for computing planetary positions in JPARSEC for an observer located on the surface of any other body, and also many improvements in accuracy. Currently the position of the Sun from Mars is computed using an ObserverElement object for Mars, and the azimuth and elevation returned agree with JPL Horizons up to the fourth decimal place presented in this post.

## Discussion

blog/computing_sun_position_from_mars.txt · created: 2010/10/27 12:18 (Last modified 2018/05/14 01:10) by Tomás Alonso Albi