Exoplanet transits from Solar System barycenter

Exoplanet transits from Solar System barycenter

To divulgate or just answer questions is something ussual (maybe much less than in other institutes) at OAN, and some astronomers here do that. In my case I'm open to questions at the OAN ephemerides server (there's an email for that in the latest version upploaded a few days ago, currently no questions at all) and also here for the JPARSEC library. I receive comments quite often and, from time to time, some interesting questions not taken into account in JPARSEC. The latest one was recently asked by Karen Collins at University of Louisville. She is working in an astronomical version of the ImageJ app (called AstroImageJ) to use it, among other aims, to reduce observations of transits of exoplanets. But before going in depth to this I would like to write a little summary of all questions I've answered so far, some were posted here but not all of them. Surprisingly I have almost one per month, and it is also nice to see that almost 50% of them came from women.

  • Since 2009 I have colaborated from time to time with Jan Kotek, who is developing a planetarium program called Asterope. I would like to obtain more free time to use some of his excellent libraries (JDBM3 combined with Healpix could be a good way to index the data of the Martian surface in a 3d app I'm working on, although I will probably use Ardor3d).
  • In October, 2010, I helped Mark Powell at JPL to obtain ephemerides for the Sun as seen for a given point on the surface of Mars. It was a really nice (but brief!) collaboration.
  • In April and September, 2011, I helped Aina Palau from Universidad Autónoma de Barcelona to use some of my models for research. She in the only astronomer that has shown interest in my technical work for research.
  • Joe Novak sent some comments about the library in different times. In May, 2011, he sent a little but nice project done using JPARSEC (AstroClock).
  • In September 2011 I helped Asrita Chetan to obtain accurate ephemerides for planets. It seem's the object was to use JPARSEC for astrology, anyway I helped …
  • In October-November, 2011, Doug Brann asked about the mods for JFreeChart at this entry of my blog. He is working in some kind of plotings of elements of satellites, maybe in a commercial app, anyway that code is released under LGPL.
  • In November, 2011, I helped Francesca de Angeli to use the TextLabel class to draw complex labels in any kind of chart, including Latex formulae, for the Gaia project. The code could become integrated inside the GaiaTools package used by different teams across Europe. Francesca is working in Gaia at Cambridge University.
  • In December, 2011, Matteo Facchin asked about the possibility of using JPARSEC in a commercial app, something explained in the section about what GPL license means in the JPARSEC page.
  • In January, 2012, Achim Brunner, a bioinformatician scientist at Max-Planck Institute for Biochemistry (Munich) asked about the modifications I did to the JMathPlot library to improve the translation/rotation behavior with mouse.
  • This month it's time for Karen Collins from Lousville University and her AstroImageJ project.

AstroImageJ is a tool used in the observatories of some universities in USA to reduced different kind of observations, and there is interest in exoplanets transits, probably because they are easy to observe since current CCD and DSLR cameras allow to reach a photometric accuracy similar to that of the Hubble Space Telescope. Maybe the most popular is HD 189733, very close to M27. What Karen asked for is a Java method to calculate at what time a given event observed from Earth would be observed from the Solar System barycenter (SSB), which is very close to the Sun. This is useful to reduce observations to an uniform time scale (and then study them using Fourier analysis or whatever) since the Earth moves around the Sun and this produces a sinusoidal displacement of +/- 8 minutes during the year. The accurate calculation of the SSB is only possible using JPL ephemerides. Since JPARSEC also considers the proper motions of the stars and this affects the timing, I decided to take also this into account in my implementation.

Testing the new implementation

Karen sent a link to a webform page written by Jason Eastman to do the same calculation. I did a very quick test using as input dates 2451545.0 and 2451645.0 UTC, from the geocenter. Output was 2451545.002244838 and 2451645.005895009. Below I include the code of the class I wrote to do the test, but a few parameters must be modified to set the dates, geocentric calculations, static star, and so on. After that, the difference between the light-time from barycenter and Earth (the two parameters calculated) and the previous output was -7.42868938208103E-4 and -7.42887813743591E-4 days. The mean value is -7.4288E-4 +/- 1E-8 days = 64.1848 +/- 0.0009 seconds. Eastman calculates what he called BJD_TDB, while I obtain the true BJD as a coordinate time. The difference is 32.184 + N, where N is the number of leap seconds since 1961. N = 34 now, but in J2000 it was 32, so the expected difference is 64.184s, which is within the millisecond. The Eastman page has a precision of 20 ms from the geocenter.

This means that my implementation has at least the same precision as that of Eastman. Of course, I also did my own tests constructing a few charts and testing numerical stability during the years. The first difference is that my implementation takes into account the stellar proper motions, something required to reduce correctly the events observed from Earth (exoplanets transits or whatever). The second difference is the accuracy itself, which I estimate in +/- 5E-10 days = +/- 50 microseconds, which is the maximum possible level using 17-digit double precision values. To improve it, BigDecimal is required. The third difference is that I don't use a planar wave approximation, so with few modifications it can be used for Solar System objects (or none at all, using the StarElement object to store the position and movement of a planet between two instants). These differences, combined with the fact that my implementation was completed in just a few hours of work, show the power of JPARSEC after six years of development.




 

The variation of the light time from SSB seems to be a straight line, but this is only apparent. The real shape it describes (during millenia) is a huge hyperbola. Karen is interested in incorporating JPARSEC inside AstroImageJ, which is great since it would be the second app using my library.

Code listings

Here is the new function in StarEphem class that calculates the light time to a star as seen from the Solar System barycenter and from the Earth. This function will be included in the next release of JPARSEC, that will provide slightly better ephemerides for stars when using high precision mode (JPL ephemerides will use now the barycenter as reference instead of the center of the Sun). To obtain a decent level of performance the new method does not correct light time iteratively. Also, for consistency inside my library, I use astrometric positions for the stars, which means that I intenctionally neglect light deflection (I think this is the only approximation). The level of precision is around +/- 50 microseconds, which is very likely good enough unless someone wants to do some crazy thing, for instance exoplanets transiting a microsecond pulsar.

private static double[] getSSBPosition(TimeElement time, ObserverElement observer, EphemerisElement eph,
		double JD_TDB, double lightTime) throws JPARSECException {
	double pos_SSB[] = null;
	String a = "JPL DE406";
	if (Configuration.PREFER_PRECISION_IN_EPHEMERIDES) {
		try {
			JPLEphemeris jpl = new JPLEphemeris(EphemerisElement.ALGORITHM.JPL_DE406);
			if (eph.algorithm.name().indexOf("JPL") >= 0) {
				a = eph.algorithm.name();
				jpl = new JPLEphemeris(eph.algorithm, Configuration.JPL_EPHEMERIDES_FILES_EXTERNAL_PATH);
			}
			pos_SSB = jpl.getGeocentricPosition(JD_TDB, TARGET.Solar_System_Barycenter, lightTime, false);
		} catch (Exception exc) {
			JPARSECException.addWarning(a+" could not be found in the classpath, using Sun position from Moshier algorithms instead.");
			pos_SSB = Ephem.eclipticToEquatorial(PlanetEphem.getGeocentricPosition(JD_TDB, TARGET.SUN, lightTime, false),
					Constant.J2000, eph.ephemMethod);
		}
	} else {
		pos_SSB = Ephem.eclipticToEquatorial(PlanetEphem.getGeocentricPosition(JD_TDB, TARGET.SUN, lightTime, false),
				Constant.J2000, eph.ephemMethod);			
	}
	// Topocentric observer should not be corrected for nutation, although the effect is 
	// well below the microsecond level. A better way is to avoid aberration in stars,
	// so that deflection is also considered.
	return Functions.substract(pos_SSB, Ephem.topocentricObserver(time, observer, eph));
}
 
/**
  * Returns the light time to a star from the Solar System Barycenter (SSB) and the Earth.
  * Precision of this method should be at the level of +/- 0.05 ms if JPL DE406 (or the selected JPL
  * integration) is available. In case the files are not found (or the ephemerides are not configured as
  * high precision at {@linkplain Configuration#PREFER_PRECISION_IN_EPHEMERIDES}), Sun position (using Moshier algorihtms)
  * is used instead of the Solar System barycenter. In that case precision
  * will be +/- 4s or better (obviously these values are relative between computations,
  * not absolute since stellar distances are not known to such level of precision).<P>
  * The difference between the light-time from SSB and from Earth (calculated from the
  * distance of the star to the observer) allows to correct time events observed from Earth
  * and refer them to an uniform time scale defined by the movement of the barycenter
  * of the Solar System in the space.<P>
  * Note that the difference between JPL DE406 and DE422 is of a few us, so precision can never be
  * better than that.
  * @param time The time.
  * @param observer The observer at Earth (possibly) used to obtain the star position.
  * @param ephIn The ephemeris properties. Only geocentric/topocentric flag and ephemeris 
  * reduction method are considered, to rest is set to J2000 equinox, ICRF frame, and astrometric
  * coordinates. The algorithm is also taken into account if it corresponds to a JPL ephemeris
  * version supported and available, otherwise DE406 is used (if available). As a last chance, 
  * Moshier is used.
  * @param star The properties of the star.
  * @return The light time to the star from the Solar System Barycenter and the Earth (second
  * component of the output array), in days.
  * @throws JPARSECException If an error occurs.
  */
public static double[] getLightTimeFromSSBandEarth(TimeElement time, ObserverElement observer, EphemerisElement ephIn,
		StarElement star) throws JPARSECException {
 
	EphemerisElement eph = ephIn.clone();
	eph.equinox = EphemerisElement.EQUINOX_J2000;
	eph.frame = EphemerisElement.FRAME.ICRF;
	eph.targetBody = TARGET.NOT_A_PLANET;
	eph.algorithm = EphemerisElement.ALGORITHM.STAR;
	eph.ephemType = EphemerisElement.COORDINATES_TYPE.GEOMETRIC; // Neglect light deflection
 
	// Get star properties to J2000 equinox and FK5/ICRF frame
	StarElement newStar = star.clone();
	if (newStar.frame == EphemerisElement.FRAME.FK4) newStar = StarEphem.transform_FK4_B1950_to_FK5_J2000(newStar);
	if (newStar.equinox != Constant.J2000) {
		// Precess to J2000
		LocationElement loc = new LocationElement(newStar.rightAscension, newStar.declination, 1.0);
		loc = LocationElement.parseRectangularCoordinates(Precession.precessToJ2000(newStar.equinox, LocationElement.parseLocationElement(loc), eph.ephemMethod));
		newStar.rightAscension = loc.getLongitude();
		newStar.declination = loc.getLatitude();
		newStar.equinox = Constant.J2000;			
	}
 
	double pc2au = Constant.PARSEC / (1000.0 * Constant.AU);
	StarEphemElement ephem_star = StarEphem.starEphemeris(time, observer, eph, newStar, false);
	double lightTimeInDaysFromEarth = ephem_star.distance * pc2au * Constant.LIGHT_TIME_DAYS_PER_AU;
 
	double JD_TDB = TimeScale.getJD(time, observer, eph, SCALE.BARYCENTRIC_DYNAMICAL_TIME);
	double pos_star[] = ephem_star.getEquatorialLocation().getRectangularCoordinates();
	pos_star = Functions.scalarProduct(pos_star, Constant.PARSEC / (1000.0 * Constant.AU)); // to AU, same as planetary ephemerides
 
	// Get topocentric/geocentric position of the Solar System barycenter
	double lightTime = 0.0;
	double pos_SSB[] = getSSBPosition(time, observer, ephIn, JD_TDB, lightTime);
 
	// Get position from SSB
	double pos_star_SSB[] = Functions.substract(pos_star, pos_SSB);
	double lightTimeInDaysFromSSB = LocationElement.parseRectangularCoordinates(pos_star_SSB).getRadius() * Constant.LIGHT_TIME_DAYS_PER_AU;
 
	// Correct from different light-time SSB-Earth (The barycenter and the star have moved in those +/- 8 minutes)
	// This effect is generally below 1s, and should be done with a StarElement referred to an
	// equinox close to J2000. This is corrected at the beginning of this method.
	// This should be an iteration, but here I simply use a first order approximation.
	lightTime = -(lightTimeInDaysFromSSB - lightTimeInDaysFromEarth);
	pos_SSB = getSSBPosition(time, observer, ephIn, JD_TDB, lightTime);
	pos_star_SSB = Functions.substract(pos_star, pos_SSB);
	lightTimeInDaysFromSSB = LocationElement.parseRectangularCoordinates(pos_star_SSB).getRadius() * Constant.LIGHT_TIME_DAYS_PER_AU;
 
	double dif = (lightTimeInDaysFromSSB - lightTimeInDaysFromEarth) * newStar.properMotionRadialV * 1000.0 / Constant.SPEED_OF_LIGHT;
	lightTimeInDaysFromSSB += dif;
 
	return new double[] {lightTimeInDaysFromSSB, lightTimeInDaysFromEarth};
}

The code of the class I produced to test the calculations is listed next.

LightTimeBarycenter.java
package research.tests;
 
import java.awt.Color;
 
import jparsec.ephem.EphemerisElement;
import jparsec.ephem.Functions;
import jparsec.ephem.Target.TARGET;
import jparsec.ephem.stars.StarElement;
import jparsec.ephem.stars.StarEphem;
import jparsec.ephem.stars.StarEphemElement;
import jparsec.graph.ChartElement;
import jparsec.graph.ChartSeriesElement;
import jparsec.graph.CreateChart;
import jparsec.graph.DataSet;
import jparsec.math.Constant;
import jparsec.math.GenericFit;
import jparsec.math.LinearFit;
import jparsec.observer.City;
import jparsec.observer.CityElement;
import jparsec.observer.ObserverElement;
import jparsec.time.AstroDate;
import jparsec.time.TimeElement;
import jparsec.time.TimeFormat;
import jparsec.time.TimeElement.SCALE;
import jparsec.util.Configuration;
import jparsec.util.JPARSECException;
 
/**
 * A program to obtain the light time to a star as seen from the Solar System Barycenter (Feb 13, 2012).
 * @author T. Alonso Albi - OAN (Spain)
 */
public class LightTimeBarycenter {
 
	/**
	 * For unit testing only.
	 */
	public static void main(String args[])
	{
		System.out.println("LightTimeBarycenter test");
 
		// Uncomment this to use Sun instead of Solar System Barycenter. Slightly lower precision (+/- 4s at most, but
		// JPL files are not required.
		// Configuration.PREFER_PRECISION_IN_EPHEMERIDES = false;
 
		try {
			// Select a star with great proper motions. Note radial velocity is -107.8 km/s, which means the
			// source moves toward us. BE CAREFUL HERE WHEN SETTING THE CORRECT PROPER MOTION IN RA, SEE
			// JAVADOC AT StarElement.properMotionRA.
			StarElement star = new StarElement("Barnard", Functions.parseRightAscension(17, 55, 23.0), Functions
					.parseDeclination("+4", 33, 18.00), 1.0 / 0.548, 9.54f, (float) (-5.0 * 15.0 * Constant.ARCSEC_TO_RAD / 100.0),
					(float) (1031.0 * Constant.ARCSEC_TO_RAD / 100.0), -107.8f, Constant.B1950, EphemerisElement.FRAME.FK4);
 
			// An example of a more realistic star at 1000 pc.  Note I use here a VERY OLD FK4 catalog, it is better to use
			// Bright Star Catalog (FK5), or Sky Master 2000 (ICRF), both referred to J2000 and already integrated into JPARSEC
			// (JPL ephems uses ICRF and J2000).
/*			StarElement star = new StarElement("HD 119288", Functions.parseRightAscension(13, 39, 44.526), Functions
					.parseDeclination("8", 38, 28.63), 1000, 2.06f, (float) (-0.0259 * 15.0 * Constant.ARCSEC_TO_RAD), (float) (-0.093 * Constant.ARCSEC_TO_RAD), 0, 
					Constant.B1950, EphemerisElement.FRAME.FK4);
*/
			// Static testing star
/*			StarElement star = new StarElement("ficticious static star", Functions.parseRightAscension(13, 39, 44.526), Functions
					.parseDeclination("8", 38, 28.63), 1000, 2.06f, (float) (-0.0259 * 15.0 * Constant.ARCSEC_TO_RAD*0), (float) (-0.093 * Constant.ARCSEC_TO_RAD*0), 0, 
					Constant.J2000, EphemerisElement.FRAME.ICRF);
*/
			// Select a date, city, and ephemeris properties.
			AstroDate astro = new AstroDate(2012, AstroDate.JANUARY, 1, 0, 0, 0);
			TimeElement time = new TimeElement(astro.jd(), SCALE.BARYCENTRIC_DYNAMICAL_TIME); // Time scale can be any other
			CityElement city = City.findCity("Madrid");
			ObserverElement observer = ObserverElement.parseCity(city);
			EphemerisElement eph = new EphemerisElement(TARGET.NOT_A_PLANET, EphemerisElement.COORDINATES_TYPE.APPARENT,
					EphemerisElement.EQUINOX_J2000, EphemerisElement.TOPOCENTRIC, EphemerisElement.REDUCTION_METHOD.IAU_2006,
					EphemerisElement.FRAME.ICRF, EphemerisElement.ALGORITHM.JPL_DE406);
 
			// Obtain the light time from barycenter
			double lightTimes[] = StarEphem.getLightTimeFromSSBandEarth(time, observer, eph, star);
			double lightTimeInDaysFromSSB = lightTimes[0];
			double lightTimeInDaysFromEarth = lightTimes[1];
 
			System.out.println("JD "+astro.jd());
			System.out.println("Light time in days from Earth: "+lightTimeInDaysFromEarth);
			System.out.println("Light time in days from SSB: "+lightTimeInDaysFromSSB);
 
			// Now repeat calculations 100 days after
			astro = new AstroDate(astro.jd() + 100);
			time.astroDate = astro;
			lightTimes = StarEphem.getLightTimeFromSSBandEarth(time, observer, eph, star);
			lightTimeInDaysFromSSB = lightTimes[0];
			lightTimeInDaysFromEarth = lightTimes[1];
 
			System.out.println("JD "+astro.jd());
			System.out.println("Light time in days from Earth: "+lightTimeInDaysFromEarth);
			System.out.println("Light time in days from SSB: "+lightTimeInDaysFromSSB);
 
			// As expected, the source is closer to us 100 days after, but light time from SSB is greater due to Earth's orbit
			// around the Sun. To get a sequence of events referred to an uniform time scale it is necessary to take as
			// reference the distance to the source at a given time and how that distance changes with time, not only the
			// difference lightTimeInDaysFromEarth-lightTimeInDaysFromSSB. The drift in case of the Barnard's star is
			// 0.042 days = 1 hour of (light) time in 100 days.
 
			// Now lets see a little chart with the variation during 10 years
			astro = new AstroDate(2012, AstroDate.JANUARY, 1, 0, 0, 0);
			int np = 365, stepDays = 10;
			double x[] = new double[np];
			double y0[] = new double[np];
			double y1[] = new double[np];
			double referenceJD = astro.jd() + stepDays;
			for (int i=0; i<np; i++) {
				System.out.println((i+1)+"/"+np);
 
				astro = new AstroDate(referenceJD + i * stepDays);
				time.astroDate = astro;
 
				lightTimes = StarEphem.getLightTimeFromSSBandEarth(time, observer, eph, star);
				lightTimeInDaysFromSSB = lightTimes[0];
				lightTimeInDaysFromEarth = lightTimes[1];
 
				x[i] = astro.jd();
				y0[i] = lightTimeInDaysFromEarth;
				y1[i] = lightTimeInDaysFromSSB;
			}
			ChartSeriesElement chartSeries1 = new ChartSeriesElement(x, y0, null, null,
					"Light-time from Earth", true, Color.BLUE, ChartSeriesElement.SHAPE_EMPTY,
					ChartSeriesElement.REGRESSION.NONE);
			ChartSeriesElement chartSeries2 = new ChartSeriesElement(x, y1, null, null,
					"Light time from SSB", true, Color.RED, ChartSeriesElement.SHAPE_EMPTY,
					ChartSeriesElement.REGRESSION.NONE);
			chartSeries1.showLines = true;
			chartSeries1.showShapes = false;
			chartSeries2.showLines = true;
			chartSeries2.showShapes = false;
			ChartElement chart = new ChartElement(new ChartSeriesElement[] {chartSeries1, chartSeries2}, 
					ChartElement.TYPE.XY_CHART, ChartElement.SUBTYPE.XY_TIME,
					"Light-time to "+star.name, 
					"Date ("+TimeFormat.getTimeScale(time)+")", "Light-time (days)", false, 800, 600);
			CreateChart c = new CreateChart(chart);
			c.showChartInJFreeChartPanel();
 
			// The blue line (drag with the mouse to zoom in/out) shows consecutive minima during June, which is 
			// correct since at that time we have the alineation Sun - Earth - Barnard. The red line from SSB shows 
			// the constant drift just in the middle, as it should be.
 
			// Now let's reduce everything to an uniform time scale, using the first point as reference
			boolean reduce = true; // In this case expected slope is 0 from SSB
			if (reduce) {
				for (int i=0; i<np; i++) {
					// Correction at first order, supposing the star is very far or static.
					double drift = (x[i] - referenceJD) * Constant.SECONDS_PER_DAY * star.properMotionRadialV * Constant.LIGHT_TIME_DAYS_PER_AU / Constant.AU;
					//double drift = y1[i]; // This would be the 'exact' correction, would yield not only slope 0, but also y1 = 0 always
					y0[i] -= drift;
					y1[i] -= drift;			
					System.out.println(y1[i]-y1[0]);
				}
			}
 
			// Now let's fit the results to a straight line
			LinearFit lf = new LinearFit(x, y1);
			lf.linearFit();
			double maxError = -1, maxErrorX = -1;
			for (int i=0; i<np; i++) {
				double predictedY = lf.evaluateFittingFunction(x[i]);
				double predictedYError = lf.evaluateFittingFunctionError(x[i]);
				double error = Math.abs(predictedY - y1[i]);
				if (error > maxError || maxError == -1) {
					System.out.println("Maximum error is "+error+" at "+x[i]);
					maxError = error;
					maxErrorX = x[i];
				}
				if (error > predictedYError) System.out.println("The point at x = "+x[i]+" is wrongly calculated ?!");
			}
			System.out.println("Maximum error found from LinearFit: "+maxError+" days at x (JD) = "+maxErrorX+" "+TimeFormat.getTimeScale(time)); 
			if (reduce) System.out.println("Maximum value of reduced light-time from SSB (should be 0): "+(DataSet.getMaximumValue(y1)-y1[0]));
			if (reduce) System.out.println("Minimum value of reduced light-time from SSB (should be 0): "+(DataSet.getMinimumValue(y1)-y1[0]));
 
			// Conclusions
			// There's always a 'maximum error found' equal or above 3.4E-5 d, seems due to a limitation in the linear fitting class.
			// With a modified HD 119288 (J2000, ICRF, and proper motions to 0) the maximum value of reduced light-time in 1st order approx is only +/- 5E-10 d in 10 yr ! 
			//    This can be considered to be the maximum internal precision. In fact, for J2000 epoch that is the precision of 17-digit double values.
			// With true HD 119288 maximum reduced light-time is 0.003 d, and with Barnard 0.004 d. These values are probably the 2nd order effect because the proper motions change with time, and amounts to 2 min per year = 1s per day.
			//    (IMPORTANT) => To avoid this issue the 'exact' correction should be used by subtracting the absolute light-time distance and not its change respect to the initial calculation time, which is not constant with time.
 
			JPARSECException.showWarnings(); // To see if Moshier algorithms were used instead of JPL ephemerides
		} catch (Exception exc) {
			exc.printStackTrace();
		}
	}
}

Discussion

Enter your comment
   ____     __ __  __   __ __   ___ 
  / __/ __ / / \ \/ /  / //_/  / _ \
 _\ \  / // /   \  /  / ,<    / // /
/___/  \___/    /_/  /_/|_|  /____/
 
 
blog/exoplanet_transits.txt · created: 2012/02/17 10:15 (Last modified 2019/04/14 16:36) by Tomás Alonso Albi
 
Recent changes RSS feed Creative Commons License Donate Powered by PHP Valid XHTML 1.0 Valid CSS Driven by DokuWiki