package newExamples.archive;
import gaia.cu1.tools.exception.GaiaException;
import gaiasimu.SimuException;
import java.io.IOException;
import newExamples.FitResults;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.optimization.OptimizationException;
import org.apache.commons.math.optimization.VectorialPointValuePair;
import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer;
public class Exoplanet {
/** exoplanet inclination function to be optimized */
private ExoplanetInclination function;
/** containing exoplanetary system */
private ExoplanetSystem system;
/**
* Constructor
*/
public Exoplanet(ExoplanetSystem system, ExoplanetInclination function) {
this.setSystem(system);
this.setFunction(function);
}
/**
* This methods fits an exoplanet inclination angle.
*
* @return a BamFitResults object
* @throws IllegalArgumentException
* @throws FunctionEvaluationException
* @throws OptimizationException
* @throws IOException
* @throws GaiaException
* @throws SimuException
*/
public FitResults fitExoplanetInclination ()
throws OptimizationException, FunctionEvaluationException, IllegalArgumentException, SimuException, GaiaException,
IOException {
// Constants
final int MAX_ITERATIONS = 1000; // Maximum interations for Gauss-Newton
final double CONVERGENCE = 1.0e-5; // Convergenge threshold to stop optimizer
final double[] FIRST_GUESS = {Math.PI / 4}; // First guess: 45 degrees
// Compute noisy values for given inclination angle
double[] errors = function.getEpochAstrometryErrors();
double[] target = function.getEpochAstrometry(true, errors);
double[] time = function.getEpochAstrometryTime();
// Weights: 1 / uncertainty ^ 2: maximum likelihood
double[] weight = new double[errors.length];
double weightSum = 0;
for (int i = 0; i < weight.length; i++){
weight[i] = 1.0 / Math.pow( errors[i], 2 );
weightSum += weight[i];
};
// Weights normalisation: <weight> = 1.0
for (int i = 0; i < weight.length; i++){
weight[i] *= weight.length / weightSum;
};
System.out.println("Time:");
for (double value : time) {
System.out.printf("%20.12e\n", value);
}
System.out.println("End time");
int nTransits;
nTransits = function.getNTransits();
System.out.printf("Number of effective transits: %4d\n", nTransits);
// Optimizer set-up
LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
optimizer.setMaxIterations( MAX_ITERATIONS );
optimizer.setCostRelativeTolerance( CONVERGENCE );
optimizer.setParRelativeTolerance( CONVERGENCE );
optimizer.setOrthoTolerance( CONVERGENCE );
// BamFunction optimization
VectorialPointValuePair optimum = optimizer.optimize(function, target, weight, FIRST_GUESS);
// Covariances and errors
double[] fitErrors = optimizer.guessParametersErrors();
double[][] fitCovariance = optimizer.getCovariances();
// Output values
FitResults fitResults = new FitResults(
optimum.getValue(), optimum.getPoint(), fitErrors, fitCovariance);
/*
System.out.printf("Inclination angle [rad]: %15.7e +- %15.7e\n",
fitResults.getVariables()[0], fitResults.getErrors()[0]);
System.out.printf("First Fits: %15.7e, %15.7e, %15.7e %15.7e\n",
fitResults.getFit()[0], fitResults.getFit()[1], fitResults.getFit()[2], fitResults.getFit()[3]);
*/
for (int i=0; i<target.length; i++) {
System.out.printf("%18.10e, %18.10e, %18.10e\n", time[i], target[i], errors[i]);
}
return fitResults;
/*
// Read output stuff
System.out.println("Values:");
double[] epochAstrometry = function.value( new double[] {0} );
for (double value : epochAstrometry){
System.out.printf("%20.12e\n", value);
}
System.out.println("End values");
System.out.println("Values2:");
double[] epochAstrometry2 = function.value( new double[] { Math.PI / 2d } );
for (double value : epochAstrometry2){
System.out.printf("%20.12e\n", value);
}
System.out.println("End values2");
/*
System.out.println("Jacobian:");
double[][] epochAstrometryJacobian = function.jacobian().value( new double[] {0.5} );
for (double[] value : epochAstrometryJacobian){
System.out.printf("%20.12e\n", value[0]);
}
System.out.println("End Jacobian");
System.out.println("Errors:");
double[] epochAstrometryErrors = function.getEpochAstrometryErrors();
for (double value : epochAstrometryErrors){
System.out.printf("%20.12e\n", value);
}
System.out.println("End errors");
*/
}
/**
* Setter for exoplanetary system
*
* @param system Exoplanetary system
*/
public void setSystem(ExoplanetSystem system) {
this.system = system;
}
/**
* Getter for exoplanetary system
*
* @return Exoplanetary system
*/
public ExoplanetSystem getSystem() {
return system;
}
/**
* Setter for exoplanet inclination function
*
* @param function Exoplanet inclination function
*/
public void setFunction(ExoplanetInclination function) {
this.function = function;
}
/**
* Getter for exoplanet inclination function
*
* @return Exoplanet inclination function
*/
public ExoplanetInclination getFunction() {
return function;
}
}