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(); } } }
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
Execuse me, I am a student who is interested in the subsolar point of the Mars. I recently used the way provided by this page to calculate the subsolar longitude of the Mars. But I find the subsolar longitude is different from the result of Mars Ephemeris Generator 2.6, which is derived from MAR097 + DE430 (https://pdsrings.seti.org/tools/ephem2_mar.html). The difference is about 3 degrees. May I have some advice about the difference?
Hi Liu,
An accurate response is difficult because I don't know in detail how these values are computed in the web page you mention. I can only try to reproduce them. First, notice the author mentions the result are referred to J2000 equinox, so the precesion of the ecliptic is not considered. I have also checked that the right ascension and declination are astrometric, not apparent. Using JPARSEC I can suggest the following code:
The output is:
The same example with the page you provide yields:
By the way, The horizons server at JPL gives this output:
As you see, although I'm not using JPL DE ephemerides the RA and DEC are exactly the same for those decimal places. The difference is 0.1 deg in the subsolar latitude and 0.25 deg in the subsolar longitude, but as I said at the beginning I can't provide more details to explain this little discrepancy. Probably related to a little different TTUT1 correction and Mars rotation/axis model. The results from Horizon are closer to the page you mention for the subsolar longitude, and closer to my results for the latitude.
Regards,
Tomás.
Thank you very much for your prompt reply and useful advice. I thought that the subsolar longitude is derived from the equation (C5) in (https://www.giss.nasa.gov/tools/mars24/help/algorithm.html). It seems that the subsolar longitude is derived from the package jparsec.ephem. I will try your method. Thank you very much.
The method on that page is also correct, but light time is not considered. The example uses Earth's time Jan 6, 2000, 00h UT, and computes the angles for that specific time on Earth, without considering that Mars was 21s of light time away that moment. This is adequate when you have a rover there but you prefer to use the time on Earth to compute what's going on there. To use Earth's time you only have to change EphemerisElement.COORDINATES_TYPE.ASTROMETRIC to EphemerisElement.COORDINATES_TYPE.GEOMETRIC in the previous code to obtain with JPARSEC:
Which agrees with that page within 0.02 degrees. The rest of the stuff including the Sun position from Mars can be obtained using an observer on Mars, for lon = lat = 0 deg would be something like:
Regards
Hello Tomás
Thank you very much for answering my question! I used the Equation (C5) to calculate the subsolar longitude. The time is 2019Jan01 00:00. The Julian date is 2458484.5. The calculated subsolar longitude is 329.090198. It is different from the 326.811 by about 2 degrees. The time of light, 21 seconds, may not account for this error. The distance between the Mars and the Sun are in good aggrement. My code is written by Matlab. It is easy to read. I have checked the parameters with the website and the paper (Allion, 1997 and Allion et al. 2000). I am not sure where my mistake is. The code is as follows. May be I need to check it more carefully. Thank you very much!
Best Regards
Liu, N
I suggest to apply the TTUT1 correction in section A, since you are using the UT date (JD) in equations that should use the date in TT scale. The rest is just testing the result of each line of code until it is correct.
Thank you for your reply. I have checked my code When the time on Earth was 00:00:00 on Jan. 6, 2000 (UTC). The result is also provided in the part II of the website(https://www.giss.nasa.gov/tools/mars24/help/algorithm.html). My subsolar Longitude is the same with the website, 174.72600 degree. But the result in another website (https://pdsrings.seti.org/tools/ephem2_mar.html) is 170.938. The difference is about 4 degrees. The result really confuse me. Maybe there is some difference between the formula and the ephemeris. Thank you again for your useful advice.
Hello Tomás
Execuse me, sir. I want to use the code provided above to obtain the subsolar latitude and longitude of the Mars because it matches the DE430. I have some problems when installing the code. There are so many packages provided at (http://conga.oan.es/~alonso/jparsec/lib/). Error occurs that package org.stathissideris.ascii2image.graphics does not exist when I download the package org already. In addition, I can not find the package 'gov.noaa.pmel.sgt' in the website above. Symbol JPlotLayout can not be found at jparsec.garph.CreateGridChart although I have download the package jparsec.jar and added it to the directory. I am new in Java. May I have some advice from you.
Regards Liu N.
Hello Liu. For (basic) ephemerides calculations you don't need anything from the lib directory, jparsec.jar is enough. You could need other files like eop.jar for Earth Orientation Parameters, or jpl_ephem.jar for some DE4xx integrations (or vsop.jar, elp2000 for other theories), but they are not required when the greatest accuracy is not needed. In case you install the code instead of simply adding jparsec.jar to the classpath of your program, you will see a lot of errors because many libraries from the lib directory need to be downloaded and added to the classpath of the project, but none of them are required for ephemerides, you can just ignore the compilation errors. In case of problems just add jparsec.jar to the classpath instead of creating a project with the source code of the library.
Here is a basic program that only uses jparsec.jar. The algorithm is set to Moshier which uses data already integrated in jparsec.jar, is a low precision fit to DE405. The optimization for speed in the eph object will not use eop.jar, so that only jparsec.jar is needed.
Hello Tomás Your advice is really useful. The problem has been solved. Thank you very much.
Regards Liu N