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 ObserverElement observer = new ObserverElement("", Math.toRadians(lon), Math.toRadians(lat), 0); 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); double pa = Functions.normalizeRadians(LocationElement.getPositionAngle(refPos, subSolar)); 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 subEarth 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); System.out.println("Earth azimuth "+Functions.normalizeRadians(pa) * Constant.RAD_TO_DEG); System.out.println("Earth elevation "+elev * Constant.RAD_TO_DEG); */ // 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.0E7*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 // 2000Jan06 00:00 191.1564 64.5049 191.1722 64.5053 191.0398 64.2616 // 2004Jan03 13:46:31 179.9952 62.0741 180.0531 62.0770 179.9890 61.9392 // 2010Jan06 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(); } } }
Discussion