package newExamples;
import gaia.cu1.params.GaiaParam;
import gaia.cu1.params.MissingParams;
import gaia.cu1.tools.exception.GaiaException;
import gaia.cu1.tools.numeric.algebra.GVector2d;
import gaia.cu1.tools.numeric.algebra.GVector3d;
import gaiasimu.GaiaSimuTime;
import gaiasimu.SimuException;
import gaiasimu.GaiaSimuEnums.FoV;
import gaiasimu.gaia.Gaia;
import gaiasimu.gaia.spacecraft.Transit;
import gaiasimu.universe.source.AstrometricParam;
import gaiasimu.universe.source.stellar.ComponentAstrometry;
import gaiasimu.universe.source.stellar.ComponentPhotometry;
import gaiasimu.universe.source.stellar.ExoPlanet;
import gaiasimu.universe.source.stellar.ExoPlanetPhysicalParameters;
import gaiasimu.universe.source.stellar.OrbitalParams;
import gaiasimu.universe.source.stellar.SpectralType;
import gaiasimu.universe.source.stellar.Star;
import gaiasimu.universe.source.stellar.StarPhysicalParameters;
import gaiasimu.universe.source.stellar.StarSystem;
import gaiasimu.universe.source.stellar.StarSystemSimuException;
import gaiasimu.universe.source.stellar.StellarAstrometry;
import gaiasimu.universe.source.stellar.StellarPhysicalParameters;
import gaiasimu.universe.source.stellar.StellarSource;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Date;
import java.util.Random;
/**
* This class Model represents the model of a particular system. It takes care of all GaiaSimu
* related simulation procedures and provides an interface to obtain a simulated orbit from GaiaSimu.
*
* @author ebopp
*/
public class Model {
/** An array of the planets in the system */
private ArrayList<ExoPlanet> planets;
/** The star of the system */
private Star star;
/** GaiaSimu star system */
private StarSystem system;
/** System astrometry */
private StellarAstrometry astrometry;
/** Reference astrometry for coordinate transformation */
private StellarAstrometry refAstro;
/** Simulation to which the model belongs; used for Gaia interaction */
private Simulation sim;
/** Array list of transits */
private ArrayList<Transit> transits;
/** Array list of filtered transits */
private ArrayList<Transit> filteredTransits;
/** Limit for across scan information */
public final double ACLIMIT = 12.0;
/**
* Default constructor.
*
* @param sim The simulation
*/
public Model ( Simulation sim ) {
planets = new ArrayList<ExoPlanet>();
this.sim = sim;
transits = null;
}
// Interface to simulation
/**
* Calculates field angles of star for given transit
*
* @param tr The transit
*
* @return The 2D-vector containing the field angles (AL, AC) in mas
*
* @throws SimuException
*/
public synchronized GVector2d getFieldAngles(Transit tr) throws SimuException {
// do NOT take into account relativistic effects, GaiaSimu algorithm does NOT work as it is supposed to
boolean relativity = false;
// Get CoMRS position from transit
GaiaSimuTime time = tr.getGSTime();
GVector3d comrs = star.getAstrometry().getCoMRSPosition(time, relativity);
// Calculate local scan coordinates
double scanAngle;
Gaia gaia = sim.getGaia();
synchronized(gaia) {
scanAngle = gaia.getAttitude().getScanAngle(comrs, time);
}
StellarAstrometry astro = (refAstro != null)? refAstro: astrometry;
GVector2d lsc = Conversion.fromCoMRStoLSC(comrs, astro.getAlpha(), astro.getDelta(), scanAngle);
// Return field angles
return lsc;
}
/**
* Gets transits of the system as observed by Gaia. This implementation filters the transits
* for membership in the sub-catalogue specified by the parent simulation.
*
* @return The filtered array of transits.
*/
public synchronized ArrayList<Transit> getFilteredTransits() {
// Filter the transit list if necessary
if(filteredTransits == null) {
filteredTransits = new ArrayList<Transit>();
for(Transit tr: getTransits())
if(isPartOfSubCatalogue(tr))
filteredTransits.add(tr);
}
return filteredTransits;
}
/**
* Checks if a given transit is part of the sub-catalogue of the parent simulation.
*
* @param tr The transit to be checked
*
* @return True, iff the transit is part of the subcatalogue
*/
private boolean isPartOfSubCatalogue(Transit tr) {
// Compute subcatalogue end time (ns)
final double YR2NS = GaiaParam.Nature.JULIANYEAR_DAY * GaiaParam.Nature.DAY_SECOND * 1e9;
double missionStart = MissingParams.MISSION_START * YR2NS;
double timeSpan = sim.getSubCatalogue() * YR2NS;
double subCatalogueTimeEnd = missionStart + timeSpan;
// Return if transit is in the subcatalogue
double timeNanoSecs = (double) tr.getGSTime().getElapsedNanoSecs();
return ((timeNanoSecs >= missionStart) && (timeNanoSecs <= subCatalogueTimeEnd));
}
/**
* Gets all transits of the system.
*
* @return Array of transits
*/
public synchronized ArrayList<Transit> getTransits() {
Random rand = new Random();
int id = rand.nextInt(50000);
// Calculate transit array from inverse scanning law
if (transits != null)
return transits;
Logger.getInstance().debug("Start of getTransits with ID "+id);
transits = new ArrayList<Transit>();
try {
// Set up parameters
GVector3d vICRS = new GVector3d(astrometry.getAlpha(), astrometry.getDelta());
double halfAperture = GaiaParam.Satellite.FOV_AC * GaiaParam.Nature.DEGREE_RADIAN/2;
// Get transits for both fields of view
ArrayList<Transit> fov1Transits, fov2Transits;
Gaia gaia = sim.getGaia();
synchronized(gaia) {
fov1Transits = gaia.attitude.inverseScan(vICRS, halfAperture, gaia.getTelescope(FoV.FOV1) );
fov2Transits = gaia.attitude.inverseScan(vICRS, halfAperture, gaia.getTelescope(FoV.FOV2) );
}
// Merge in one array
transits.addAll(fov1Transits);
transits.addAll(fov2Transits);
} catch (SimuException e) {
System.err.println("SimuException occurred while calculating transits!");
e.printStackTrace();
}
Logger.getInstance().debug("End of getTransits with ID "+id+" -- #Tr. = "+transits.size());
return transits;
}
/**
* Sets the transits to a fixed list. (Intended for fit models, that should comply with the observational model.)
*
* @param transits
*/
public synchronized void setTransits(ArrayList<Transit> transits) {
this.transits = transits;
}
// Initialization methods
/**
* Adds a new planet to the system
*
* @param star
* @param massPlanet Mass of planet in Jupiter masses
* @param period Period of planetary orbit in days
* @param timePeriastron Time of periastron in JD
* @param eccentricity
* @param inclination
* @param omega2
* @param nodeAngle
* @throws SimuException
* @throws GaiaException
* @throws IOException
*/
public ExoPlanet createPlanet (
double massPlanet,
double period,
double timePeriastron,
double eccentricity,
double inclination,
double omega2,
double nodeAngle ) {
// set orbital parameters
OrbitalParams orbParam = new OrbitalParams();
double massM1 = ((StellarPhysicalParameters) star.getPhysParam()).getMass();
double massM2 = massPlanet * GaiaParam.Nature.JUPITER_MASS / GaiaParam.Nature.SUN_MASS;
double systemSeparation = Math.pow(massM1 + massM2, 1.0/3.0) * Math.pow(period/GaiaParam.Nature.JULIANYEAR_DAY, 2.0/3.0);
orbParam.semiMajorAxis = massM1 / (massM1 + massM2) * systemSeparation;
orbParam.period = period;
orbParam.timeperiastron = timePeriastron;
orbParam.eccentricity = eccentricity;
orbParam.omega2 = omega2;
if (orbParam.eccentricity == 0)
orbParam.omega2 = 0;
orbParam.nodeangle = nodeAngle;
orbParam.inclination = inclination;
// Dummy values
double tEq = 400.0; // dummy value
double plRadius = 0.1; // dummy value
// Obtain distance via parallax from astrometry
double distance = 1e3 / astrometry.getVarpi();
// TODO put a correct value for absolute visual and bolometric magnitude, vsini, V-I
double magMV2 = 32.58, MBol2 = magMV2-4; //Mag MV from Baraffe models.
double vSini = 0.0, vmi2 = 9.9999;
// Initialize ComponentAstrometry pointing to the primary and then correct with setCompanion()
ComponentAstrometry compAstrom = new ComponentAstrometry(planets.size(), orbParam);
ExoPlanetPhysicalParameters compPhys =
new ExoPlanetPhysicalParameters(tEq, plRadius, vSini, massM2, MBol2, star.getStarPhysicalParams());
ComponentPhotometry compPhot =
new ComponentPhotometry(compAstrom, magMV2, distance, vmi2, star.getPhotometry().getAbsV(), star.getPhotometry().getExtinctionCurve());
ExoPlanet planet = new ExoPlanet(star.getRootId(), compPhot, compAstrom, compPhys);
compPhot.setTargetAstroSource(planet);
compAstrom.setTargetAstroSource(planet);
compPhys.setTargetAstroSource(planet);
return this.addPlanet(planet);
}
/**
* Setter for model star.
*
* @param star The star
*/
public void setStar ( Star star ) {
this.star = star;
}
/**
* Finalizes the system. To be called after star has been set and all planets have been added.
* NOTE: Needs to be called every time, when parameters are altered!
*
* @throws StarSystemSimuException
* @throws SimuException
*/
public void finalizeSystem()
throws StarSystemSimuException, SimuException {
StellarSource[] components = new StellarSource[planets.size()+1];
components[0] = star;
int i=1;
for (ExoPlanet planet: planets) {
components[i] = planet;
i++;
}
system = new StarSystem(star.getRootId(), components, astrometry);
}
// Setters, getters and adders
/**
* Setter for astrometry given position and velocity parameters
*
* @param ra Right Ascension ICRS [deg]
* @param dec Declination ICRS [deg]
* @param parallax Parallax [mas]
* @param muRa Proper motion RA [mas/yr]
* @param muDec Proper motion DEC [mas/yr]
* @param vRad Radial velocity [km/s]
*/
public void setAstrometry ( double ra, double dec, double parallax, double muRa, double muDec, double vRad ) {
setAstrometry ( new StellarAstrometry ( new AstrometricParam ( ra, dec, parallax, muRa, muDec, vRad ) ) );
}
public boolean useAcrossScan() {
return (star.getPhotometry().getMeanV0() < ACLIMIT);
}
/**
* Setter for astrometry given StellarAstrometry object
*
* @param astrometry
*/
public void setAstrometry ( StellarAstrometry astrometry ) {
this.astrometry = astrometry;
}
/**
* Getter for astrometry
*
* @return Astrometry
*/
public StellarAstrometry getAstrometry() {
return astrometry;
}
/**
* Setter for array of planets
*
* @param planets
*/
public void setPlanets ( ArrayList<ExoPlanet> planets ) {
this.planets = planets;
}
/**
* Getter for array of planets
*
* @return Array of planets in the system
*/
public ArrayList<ExoPlanet> getPlanets() {
return planets;
}
/**
* Gets a planet with a certain index.
*
* @param i Number of requested planet
* @return The i-th Planet
*/
public ExoPlanet getPlanet ( int i ) {
return planets.get(i);
}
/**
* Removes planet from system
*
* @param planet Planet to be removed
* @return true, if planet existed in system and was removed accordingly.
*/
public boolean removePlanet(ExoPlanet planet) {
return planets.remove(planet);
}
/**
* Adds a new planet to the system
*
* @param planet GaiaSimu exoplanet object
*/
public ExoPlanet addPlanet (
ExoPlanet planet ) {
planets.add(planet);
return planet;
}
/**
* Getter for star.
*
* @return Star
*/
public Star getStar() {
return star;
}
public void setRefAstro ( double ra, double dec, double parallax, double muRa, double muDec, double vRad ) {
setRefAstro ( new StellarAstrometry ( new AstrometricParam ( ra, dec, parallax, muRa, muDec, vRad ) ) );
}
public void setRefAstro(StellarAstrometry refAstro) {
this.refAstro = refAstro;
}
public StellarAstrometry getRefAstro() {
if ( refAstro == null )
return astrometry;
return refAstro;
}
/**
* Gets biggest planetary astrometric signature.
*
* @return Astrometric signature in muas
*/
public double getAstroSig() {
// Initialize array list of astrometric signatures
ArrayList<Double> sigs = new ArrayList<Double>();
// Mass of star in solar masses
double massStar = ((StarPhysicalParameters) star.getPhysParam()).getMass();
// Parallax in mas
double parallax = astrometry.getVarpi();
// Loop through planets
for(ExoPlanet pl: planets) {
// Mass of the planet in solar masses
double massPlanet = ((ExoPlanetPhysicalParameters) pl.getPhysParam()).getMass();
try {
// Period of the planetary orbit in years
double period = ((ComponentAstrometry) pl.getAstrometry()).getOrbitalParams().period / GaiaParam.Nature.JULIANYEAR_DAY;
// Calculate astrometric signature in microarcseconds
double sig = massPlanet * Math.pow(period / massStar, 0.666666666e0) * parallax * 1e3;
sigs.add(sig);
} catch (SimuException e) {
// Just skip it, no handling possible
}
}
// Find maximum and return it
double maxSig = 0.0;
for(Double sig: sigs)
maxSig = (sig > maxSig)? sig: maxSig;
return maxSig;
}
}