//
//++++++++++++++++
//.IDENTIFICATION $Id: dmstest.java,v 2.0 2002/01/09 17:04:28 arenou Exp $
//.TYPE		class
//.LANGUAGE	java
//.AUTHOR	F.Arenou
//.ENVIRONMENT	tested on MacOS X, should run anywhere ;-)
//.KEYWORDS	orbits, Gaia
//.PURPOSE	this is a test program of the classes dms_simu & dms_orbit
//.INPUT	choose 1, 2, 3, 4 or  5 on command line to run one of the tests
//----------------
//
import java.text.* ;
import java.io.*;

public class dmstest {

    //
    // choose which test
    //
    public static void main (String args[]) {
        String myarg="0" ; if (args.length != 0) myarg=args[0] ;
        if (myarg.startsWith("1")) creates_one_random_system() ;
        else if (myarg.startsWith("2")) creates_a_random_sample() ;
        else if (myarg.startsWith("3")) test_a_true_star() ;        
        else if (myarg.startsWith("4")) test_a_multiple_planet() ;
        else if (myarg.startsWith("5")) test_a_LTT_eclipsing() ;
        else System.out.println("Usage: dmstest #, with 1<=#<=5") ;
    }


    //
    // this test creates a random system (if the star is not single)
    // then prints the information about this system
    // and the position offset with respect to barycenter at some date
    //
    public static void creates_one_random_system () {
        // creates a random new system given the absolute mag of primary 
        dms_simu my_system = new dms_simu(4.77) ; // abs magnitude of a G star
        int number_comp = my_system.getNumberOfCompanions() ;
        System.out.println("A system randomly chosen:") ;
        System.out.println("Mass of primary = " + my_system.getMass(0)) ;
        System.out.println("V-I of primary = " + my_system.getVmICol(0)) ;
        if (number_comp == 0) {
            System.out.println("This is a single star") ;
        } else {
            System.out.println("This is a system with a primary and " + 
            number_comp + " companion(s)") ;
            for (int i=1; i<= number_comp; i++) {
                System.out.println("Companion " + i) ;
                System.out.println("\tMass = " + my_system.getMass(i)) ;
                System.out.println("\tRadius = " + my_system.getRadius(i)) ;
                System.out.println("\tV absolute magnitude = " + 
                my_system.getAbsMag(i)) ;
                System.out.println("\tV-I color = " + my_system.getVmICol(i)) ;
                my_system.printCompanionInfo(i) ;
            }
            System.out.println("Position wrt barycenter at day "+
            "(e.g.) -100 will be (in A.U.):") ;
            my_system.OrbitalPosition(-100.) ; // computes position of all comp. 
            for (int i=0; i<= number_comp; i++) {
                System.out.println("Component " + i + ", ksi=" + 
                my_system.getKsi(i) + ", eta=" + my_system.getEta(i)) ;
            }
        }
    }

    //
    // statistics for a large sample
    //
    public static void creates_a_random_sample () {
        DecimalFormat fmt = new DecimalFormat("##.####");
        int SampleSize = 1000 ;			// of 1000 stars
        dms_simu [] sample = new dms_simu[SampleSize] ; 		
        try {
            String FileName = "sample.system";	// name to be changed...
            PrintWriter OutputFile=new PrintWriter(new 
                    BufferedWriter(new FileWriter(FileName)));
            System.out.println("Writing parameters of a simulated sample "+
            "in file "+FileName) ;
            System.out.println("Format: number of component, then "+
            "for each component: mass, absmag, V-I, radius, period");
            for (int star=0; star<SampleSize; star++) {
                sample[star] = new dms_simu(4.77) ; 		// G-type primary
                int number_comp = sample[star].getNumberOfCompanions() ;
                OutputFile.print(number_comp) ;
                for (int i=0; i<= number_comp; i++) { 
                    OutputFile.print("\t" + 
                    fmt.format(sample[star].getMass(i)) + "\t" +
                    fmt.format(sample[star].getAbsMag(i)) + "\t" + 
                    fmt.format(sample[star].getVmICol(i)) + "\t" +
                    fmt.format(sample[star].getRadius(i)) + "\t" +
                    fmt.format(sample[star].getCompanionInfo(i)[0])) ;	// e.g. 
                    // abs mag and V-I should be ignored for BD and EP !!!
                }
                OutputFile.println() ;
            }
            OutputFile.close();
        } 
        catch (Exception e) { e.printStackTrace(); } // in case of file IO error
    }


    //
    // shows the astrometric motion of a 900 day period binary (HIP 39903)
    //
    public static void test_a_true_star () {
        // P (days), T (HJD-2440000), e, omega_2 (deg), i (deg), Omega (deg)
        double [] param39903 = { 899.2761,1845.6018,0.1117,136.49+180.,21.81,186.85 } ;
        double massprimary = 1.26, masssecond = 0.434 ;	// masses primary& secondary
        double plx39903 = 48.941 ; 		// parallax in mas
        double delta_alf0, delta_del0 ;		// reflex motion on the sky
        double delta_alf1, delta_del1 ;		// secondary motion on the sky
        double rho ;				// angular separation
        DecimalFormat fmt = new DecimalFormat("######.###");

        System.out.println("This is the astrometric binary HIP 39903:") ;
        // We could have constructed in one instruction the system
        // but we show how it can be done in 3 steps
        dms_orbit double_star = new dms_orbit() ;	// empty construction
        double_star.primary(1, massprimary) ; 		// primary with 1 companion
        double_star.companion(1, masssecond, param39903) ; // then init secondary 

        // foreach time (interval of Hipparcos observation in HJD-2440000),
        // the reflex motion of the primary in milliarcsec is printed        
        try {
            String FileName = "HIP39903.pos";	// to be changed...
            PrintWriter OutputFile=new PrintWriter(new 
                    BufferedWriter(new FileWriter(FileName)));
            System.out.println("Writing the reflex motion alfa, delta and rho (mas) "+
            "versus time (HJD-2440000) in file "+FileName) ;
            for (double t=7800; t<9100; t+=10) {
                // compute position on tangent plane of all components at time t
                double_star.OrbitalPosition(t) ; 
                // get position of primary and convert it to angular coordinates
                delta_alf0 = plx39903*double_star.getKsi(0) ;	// in alpha
                delta_del0 = plx39903*double_star.getEta(0);	// and delta
                rho = plx39903*double_star.getRho(1);		// separation 0-1 (AB)
                OutputFile.println(t + "\t" + fmt.format(delta_alf0) + "\t" + 
                fmt.format(delta_del0) + "\t" + rho) ;
            }
            OutputFile.close();
        } 
        catch (Exception e) { e.printStackTrace(); } // in case of file IO error
    }
    

    //
    // orbital parameters of star Ups And, assuming i=75, O=0 for all
    //
    public static void test_a_multiple_planet () {
        int number_comp = 3 ;			// 3 planets
        // orbital parameters of star Ups And, assuming i=75, O=0 for all
        // P (days), T (HJD-2450000), e, omega_2 (deg), i (deg), Omega (deg)
        double [][] paramsUpsAnd = { {}, 	// nothing for the star then
        {    4.617, 2.24, 0.034, 83  +180,75,0 },// orbital param for the 3 planets
        {  241.2,  154.9, 0.18, 243.6+180,75,0 }, 
        { 1308.5,     44, 0.31, 255  +180,75,0 } } ;
        double [] massesUpsAnd = {1.3,  0.000678, 0.00196, 0.004097} ; // solar masses

        // creates the system: number of companions, masses, orbital parameters 
        dms_orbit planet_star=new dms_orbit(number_comp,massesUpsAnd,paramsUpsAnd) ; 
        
        System.out.println("This is the Ups And system with a primary and " + 
            number_comp + " planets") ;
        for (int i=1; i<= number_comp; i++) {
            System.out.println("Companion " + i) ;
            System.out.println("\tMass = " + planet_star.getMass(i)) ;
            planet_star.printCompanionInfo(i) ;
        }
        
        try {
            String FileName = "upsand.rv";	// to be changed...
            PrintWriter OutputFile=new PrintWriter(new 
                    BufferedWriter(new FileWriter(FileName)));
            // the reflex motion of the primary in m/s is printed
            System.out.println("Writing the radial velocity (m/s) of primary "+
            "versus time (HJD-2450000) in file "+FileName) ;
            for (double t=0; t<3500; t+=0.5) {
                // compute position of all components at time t
                planet_star.OrbitalPosition(t) ; 
                // get radial velocity
                double delta_vr = 1000*planet_star.getVra(0) ;	// in m/s
                OutputFile.println(t + "\t" + delta_vr) ;
            }
            OutputFile.close();
        } 
        catch (Exception e) { e.printStackTrace(); } // in case of file IO error
    }


    //
    // star RCMa, assuming spherical stars, an eclipsing binary with a third comp.
    //
    public static void test_a_LTT_eclipsing () {
        int number_comp = 2 ; 
        double [][] params = { {}, 
        { 1.13594197,-9563.4193, 0,3.612+180,79.47,0 }, // Omega ?, T=primary eclipse
        { 33892.569,9349.145,0.494,10.55+180,91.79,262.52 } } ; // T-2440000
        params[1][1] = Primary2Periastron(params[1][1],params[1][0],
            params[1][2],params[1][3]) ; // correct: periastron instead of primary eclipse
        double [] masses = {1.07, 0.17, 0.342} ;	// in M_sun
        double [] radii = {1.57, 1.06, 0.5} ; 		// in R_sun, 
        double [] Vmag = {5.7, 9, 15} ; 		// magnitudes (invented)
        int [] RefPrimary = {0,0,0} ;	// anything else not yet implemented.
        DecimalFormat fmt = new DecimalFormat("######.####");

        // create the new system: number of companions, masses, orbital param 
        dms_orbit eclipsing_star =
            new dms_orbit(number_comp,masses,params,RefPrimary,radii) ; 
        System.out.println("This is (approx.) the RCMa system with an eclipsed "+
            "primary and " + number_comp + " companions, with LTT effect") ;
        for (int i=1; i<= number_comp; i++) {
            if (eclipsing_star.isEclipsing(i)) 
            System.out.println("Companion " + i + " (will eclipse primary)") ;
            else System.out.println("Companion " + i + " (will not eclipse primary)") ;
            System.out.println("\tMass = " + eclipsing_star.getMass(i) +
            "\tRadius = " + eclipsing_star.getRadius(i)) ;
            eclipsing_star.printCompanionInfo(i) ;
        }
        
        try {
            String FileName = "rcma.lc";	// to be changed...
            System.out.println("Writing the light curve (mag) versus "+
            "time (days) in file "+FileName) ;
            PrintWriter OutputFile=new PrintWriter(new 
                    BufferedWriter(new FileWriter(FileName)));
            double flux0 = Math.pow(10,-0.4*Vmag[0]) ;
            double flux1 = Math.pow(10,-0.4*Vmag[1]) ;
            for (double t=-9564.2; t<-9563; t+=0.0001) { 
                eclipsing_star.OrbitalPosition(t) ; // around the primary eclipse T0
                double mag = -2.5*Math.log( eclipsing_star.getFlx(0)*flux0 +
                                eclipsing_star.getFlx(1)*flux1 )/Math.log(10) ;
                OutputFile.println(fmt.format(t) + "\t" + mag) ;
            }
            OutputFile.close();
        } 
        catch (Exception e) { e.printStackTrace(); } // in case of file IO error
        
        try {
            String FileName = "rcma.ltt";	// to be changed...
            System.out.println("Writing light-Time Travel (s) for primary "+
            "vs time (year) in file "+FileName) ;
            PrintWriter OutputFile=new PrintWriter(new 
                    BufferedWriter(new FileWriter(FileName)));
            for (double t=-30000; t<12000; t+=100) {
                // compute position of all components at time t
                eclipsing_star.OrbitalPosition(t) ; 
                // get LTT
                double delta_ltt = 86400*eclipsing_star.getLtt(0) ;// in seconds
                double year = (t-8349.0625)/365.25+1991.25 ;	// in julian year
                OutputFile.println(fmt.format(year) + "\t" + fmt.format(delta_ltt));
            }
            OutputFile.close();
        } 
        catch (Exception e) { e.printStackTrace(); } // in case of file IO error

    }

    // periastron date from date of primary eclipse
    public static double Primary2Periastron(double Tp, double P, 
                        double e, double omega2deg) {
        double Pi = 3.141592653589793 ; 
        double nup = Pi/2 - Pi*(omega2deg-180)/180 ; 
        double Ep = 2*Math.atan2(Math.sin(nup/2)*Math.sqrt((1-e)/(1+e)),
                    Math.cos(nup/2)) ;
        return Tp-P*(Ep-e*Math.sin(Ep))/(2*Pi) ;
    }

}