/*
TroisCorps - algorithme de simulation pour SimuApplet

Copyright (C) 2002 Observatoire de Paris-Meudon

Ce programme est un logiciel libre ; vous pouvez le redistribuer et/ou le modifier conformément aux dispositions de la Licence Publique Générale GNU, telle que publiée par la Free Software Foundation ; version 2 de la licence, ou encore (à votre choix) toute version ultérieure.

Ce programme est distribué dans l'espoir qu'il sera utile, mais SANS AUCUNE GARANTIE ; sans même la garantie implicite de COMMERCIALISATION ou D'ADAPTATION A UN OBJET PARTICULIER. Pour plus de détail, voir la Licence Publique Générale GNU .

Vous devez avoir reçu un exemplaire de la Licence Publique Générale GNU en même temps que ce programme ; si ce n'est pas le cas, écrivez à la Free Software Foundation Inc., 675 Mass Ave, Cambridge, MA 02139, Etats-Unis.
*/

import java.util.Vector;

import simu.*;

/**
 * Algorithme de simulation pour SimuApplet.
 * Problème à trois corps.
 * compilation: javac -source 1.2 -target 1.1 -encoding ISO-8859-1 -classpath SimuLab.jar TroisCorps.java
 */
public class TroisCorps extends SimuApplet {

    static final double G = 6.67e-11;
    double x2, y2, z2, x3, y3, z3, x20, y20, z20, x30, y30, z30;
    double xp2, yp2, zp2, xp3, yp3, zp3, xp20, yp20, zp20, xp30, yp30, zp30;
    double xs2, ys2, zs2, xs3, ys3, zs3;
    double m1, m2, m3;
    double d12, d23, d31, g12, g31, A, B, C, D, T, Tps, T0=0, h=600;
    double k11, k12, k13, k14, k15, k16, l11, l12, l13, l14, l15, l16; 
    double k21, k22, k23, k24, k25, k26, l21, l22, l23, l24, l25, l26; 
    double k31, k32, k33, k34, k35, k36, l31, l32, l33, l34, l35, l36; 
    double k41, k42, k43, k44, k45, k46, l41, l42, l43, l44, l45, l46; 
    
    public TroisCorps() {
        appletDoc("troiscorps_doc.xml");
        /* nouvelle version de SimuLab : les paramètres sont dans le fichier XML
        // paramètres d'entrée
        ParamIn x20 = new ParamIn();
        x20.type = "nombre";
        x20.label = "x20";
        x20.titre = "abscisse initiale du corps 2";
        x20.unit = "km";
        x20.acquisition = "champ";
        x20.defaut = "384600000";
        addParamIn(x20);
        ParamIn y20 = new ParamIn();
        y20.type = "nombre";
        y20.label = "y20";
        y20.titre = "ordonnée initiale du corps 2";
        y20.unit = "km";
        y20.acquisition = "champ";
        y20.defaut = "0";
        addParamIn(y20);
        ParamIn z20 = new ParamIn();
        z20.type = "nombre";
        z20.label = "z20";
        z20.titre = "cote initiale du corps 2";
        z20.unit = "km";
        z20.acquisition = "champ";
        z20.defaut = "0";
        addParamIn(z20);
        
        ParamIn x30 = new ParamIn();
        x30.type = "nombre";
        x30.label = "x30";
        x30.titre = "abscisse initiale du corps 3";
        x30.unit = "km";
        x30.acquisition = "champ";
        x30.defaut = "1197000";
        addParamIn(x30);
        ParamIn y30 = new ParamIn();
        y30.type = "nombre";
        y30.label = "y30";
        y30.titre = "ordonnée initiale du corps 3";
        y30.unit = "km";
        y30.acquisition = "champ";
        y30.defaut = "-9928000";
        addParamIn(y30);
        ParamIn z30 = new ParamIn();
        z30.type = "nombre";
        z30.label = "z30";
        z30.titre = "cote initiale du corps 3";
        z30.unit = "km";
        z30.acquisition = "champ";
        z30.defaut = "0";
        addParamIn(z30);

        ParamIn xp20 = new ParamIn();
        xp20.type = "nombre";
        xp20.label = "xp20";
        xp20.titre = "vitesse en x initiale du corps 2";
        xp20.unit = "km/s";
        xp20.acquisition = "champ";
        xp20.defaut = "0";
        addParamIn(xp20);
        ParamIn yp20 = new ParamIn();
        yp20.type = "nombre";
        yp20.label = "yp20";
        yp20.titre = "vitesse en y initiale du corps 2";
        yp20.unit = "km/s";
        yp20.acquisition = "champ";
        yp20.defaut = "1025";
        addParamIn(yp20);
        ParamIn zp20 = new ParamIn();
        zp20.type = "nombre";
        zp20.label = "zp20";
        zp20.titre = "vitesse en z initiale du corps 2";
        zp20.unit = "km/s";
        zp20.acquisition = "champ";
        zp20.defaut = "0";
        addParamIn(zp20);

        ParamIn xp30 = new ParamIn();
        xp30.type = "nombre";
        xp30.label = "xp30";
        xp30.titre = "vitesse en x initiale du corps 3";
        xp30.unit = "km/s";
        xp30.acquisition = "champ";
        xp30.defaut = "8490";
        addParamIn(xp30);
        ParamIn yp30 = new ParamIn();
        yp30.type = "nombre";
        yp30.label = "yp30";
        yp30.titre = "vitesse en y initiale du corps 3";
        yp30.unit = "km/s";
        yp30.acquisition = "champ";
        yp30.defaut = "-2455";
        addParamIn(yp30);
        ParamIn zp30 = new ParamIn();
        zp30.type = "nombre";
        zp30.label = "zp30";
        zp30.titre = "vitesse en z initiale du corps 3";
        zp30.unit = "km/s";
        zp30.acquisition = "champ";
        zp30.defaut = "100";
        addParamIn(zp30);

        ParamIn m10 = new ParamIn();
        m10.type = "nombre";
        m10.label = "m10";
        m10.titre = "masse du corps 1";
        m10.unit = "kg";
        m10.acquisition = "champ";
        m10.defaut = "5.974e24";
        addParamIn(m10);
        ParamIn m20 = new ParamIn();
        m20.type = "nombre";
        m20.label = "m20";
        m20.titre = "masse du corps 2";
        m20.unit = "kg";
        m20.acquisition = "champ";
        m20.defaut = "7.348e22";
        addParamIn(m20);
        ParamIn m30 = new ParamIn();
        m30.type = "nombre";
        m30.label = "m30";
        m30.titre = "masse du corps 3";
        m30.unit = "kg";
        m30.acquisition = "champ";
        m30.defaut = "10000";
        addParamIn(m30);

        ParamIn tps = new ParamIn();
        tps.type = "nombre";
        tps.label = "tps";
        tps.titre = "temps de simulation (pour plot uniquement)";
        tps.unit = "min";
        tps.acquisition = "champ";
        tps.defaut = "3000";
        addParamIn(tps);

        // paramètres de sortie
        ParamOut x2 = new ParamOut();
        x2.type = "nombre";
        x2.label = "x";
        addParamOut(x2);
        ParamOut y2 = new ParamOut();
        y2.type = "nombre";
        y2.label = "y";
        addParamOut(y2);
        ParamOut x3 = new ParamOut();
        x3.type = "nombre";
        x3.label = "x3";
        addParamOut(x3);
        ParamOut y3 = new ParamOut();
        y3.type = "nombre";
        y3.label = "y3";
        addParamOut(y3);
        
        // options d'affichage
        Affichage p = new Affichage();
        p.type = "plotlive";
        p.titre = "mouvement des corps 2 (rouge) et 3 (bleu) live";
        String[][] pparams = {{"x", "y"},{"x3","y3"}};
        p.params = pparams;
        p.fond = "terre.gif";
        p.xdebut="-4.5e8";
        p.xfin="4.5e8";
        p.ydebut="-4.5e8";
        p.yfin="4.5e8";
        p.zooming=true;
        addAffichage(p);
        
        p = new Affichage();
        p.type = "plot";
        p.titre = "mouvement des corps 2 (rouge) et 3 (bleu)";
        p.params = pparams;
        p.fond = "terre.gif";
        p.xdebut="-4.5e8";
        p.xfin="4.5e8";
        p.ydebut="-4.5e8";
        p.yfin="4.5e8";
        p.zooming=true;
        addAffichage(p);

        p = new Affichage();
        p.type = "tableau";
        p.titre = "mouvement des corps 2 et 3";
        addAffichage(p);
        */
    }
    
    // calcul live (utilisé avec affichage d'un plot live)
    public void initCalculLive(ListeValeurs in) throws SimuException {
        initParamIn(in);        
    }
    
    public ListeValeurs calculLive() throws SimuException {
        ListeValeurs out = new ListeValeurs();
        // comme on ne connaît pas les rayons, on prend r(km)=masse(kg)/10^21,
        double r1 = m1/1e21;
        double r2 = m2/1e21;
        double r3 = m3/1e21;
        if((Math.sqrt((x2-x3)*(x2-x3) + (y2-y3)*(y2-y3) + (z2-z3)*(z2-z3)) < r2 + r3) ||
            (Math.sqrt(x2*x2 + y2*y2 + z2*z2) < r1+r2) ||
            (Math.sqrt(x3*x3 + y3*y3 + z3*z3) < r1+r3)) {
                System.out.println("Collision entre deux corps");
                return null;
            }

        rk4();

        out.ajouterDouble("x", x20);
        out.ajouterDouble("y", y20);
        out.ajouterDouble("x3", x30);
        out.ajouterDouble("y3", y30);
        //try { // pour bloquer le programme
            //Thread.currentThread().sleep(Math.round(5));
        //} catch (InterruptedException e) {}
        return out;
    }
    
    // calcul complet (pas utilisé avec affichage d'un plot live)
    public ListeValeurs calcul(ListeValeurs in) throws SimuException {
        ListeValeurs out = new ListeValeurs();
        
        T=0;
        T0=0;
        initParamIn(in);
        Tps*=60;
        
        double r1 = m1/1e21;
        double r2 = m2/1e21;
        double r3 = m3/1e21;
        while(T<Tps) {
        // comme on ne connaît pas les rayons, on prend r(km)=masse(kg)/10^21,
            if((Math.sqrt((x2-x3)*(x2-x3) + (y2-y3)*(y2-y3) + (z2-z3)*(z2-z3)) < r2 + r3) ||
               (Math.sqrt(x2*x2 + y2*y2 + z2*z2) < r1+r2) ||
               (Math.sqrt(x3*x3 + y3*y3 + z3*z3) < r1+r3)) {
                System.out.println("Collision entre deux corps");
                break;
            }
            T=T0;
            rk4();
            T0=T;
            
            out.ajouterDouble("x", x20);
            out.ajouterDouble("y", y20);
            out.ajouterDouble("x3", x30);
            out.ajouterDouble("y3", y30);
        }
        
        return out;
    }
    
    //calcul des dérivées secondes
    public void calcDerSec() {
        d12=Math.pow((x2*x2+y2*y2+z2*z2), 1.5);
        d23=Math.pow(((x3-x2)*(x3-x2)+(y3-y2)*(y3-y2)+(z3-z2)*(z3-z2)), 1.5);
        d31=Math.pow((x3*x3+y3*y3+z3*z3), 1.5);
        g12=C/d12;
        g31=D/d31;
        xs2=B*((x3-x2)/d23-x3/d31)-g12*x2;
        ys2=B*((y3-y2)/d23-y3/d31)-g12*y2;
        zs2=B*((z3-z2)/d23-z3/d31)-g12*z2;
        xs3=A*((x2-x3)/d23-x2/d12)-g31*x3;
        ys3=A*((y2-y3)/d23-y2/d12)-g31*y3;
        zs3=A*((z2-z3)/d23-z2/d12)-g31*z3;
        }
    
    //initialisation des paramètres d'entrée
    public void initParamIn(ListeValeurs in) throws SimuException{
        
        x2 = in.lireDouble("x20");
        x20=x2;
        y2 = in.lireDouble("y20");
        y20=y2;
        z2 = in.lireDouble("z20");
        z20=z2;
        
        x3 = in.lireDouble("x30");
        x30=x3;
        y3 = in.lireDouble("y30");
        y30=y3;
        z3 = in.lireDouble("z30");
        z30=z3;
        
        xp2 = in.lireDouble("xp20");
        xp20=xp2;
        yp2 = in.lireDouble("yp20");
        yp20=yp2;
        zp2 = in.lireDouble("zp20");
        zp20=zp2;
        
        xp3 = in.lireDouble("xp30");
        xp30=xp3;
        yp3 = in.lireDouble("yp30");
        yp30=yp3;
        zp3 = in.lireDouble("zp30");
        zp30=zp3;
        
        m1 = in.lireDouble("m10");
        m2 = in.lireDouble("m20");
        m3 = in.lireDouble("m30");
        
        Tps = in.lireDouble("tps");
        
        A=G*m2;
        B=G*m3;
        C=G*(m1+m2);
        D=G*(m1+m3);
    }
    
    //méthode de Runge-Kutta d'ordre 4
    public void rk4(){
        x2=x20;
        y2=y20;
        z2=z20;
        
        x3=x30;
        y3=y30;
        z3=z30;
        
        xp2=xp20;
        yp2=yp20;
        zp2=zp20;
        
        xp3=xp30;
        yp3=yp30;
        zp3=zp30;

        calcDerSec();
        k11=h*xp2;
        k12=h*yp2;
        k13=h*zp2;
        k14=h*xp3;
        k15=h*yp3;
        k16=h*zp3;
        
        l11=h*xs2;
        l12=h*ys2;
        l13=h*zs2;
        l14=h*xs3;
        l15=h*ys3;
        l16=h*zs3;
        
        T=T0+h/2;
        
        x2=x20+k11/2;
        y2=y20+k12/2;
        z2=z20+k13/2;
        x3=x30+k14/2;
        y3=y30+k15/2;
        z3=z30+k16/2;
        
        xp2=xp20+l11/2;
        yp2=yp20+l12/2;
        zp2=zp20+l13/2;
        xp3=xp30+l14/2;
        yp3=yp30+l15/2;
        zp3=zp30+l16/2;
        
        calcDerSec();
        
        k21=h*xp2;
        k22=h*yp2;
        k23=h*zp2;
        k24=h*xp3;
        k25=h*yp3;
        k26=h*zp3;
        
        l21=h*xs2;
        l22=h*ys2;
        l23=h*zs2;
        l24=h*xs3;
        l25=h*ys3;
        l26=h*zs3;
        
        x2=x20+k21/2;
        y2=y20+k22/2;
        z2=z20+k23/2;
        x3=x30+k24/2;
        y3=y30+k25/2;
        z3=z30+k26/2;
        
        xp2=xp20+l21/2;
        yp2=yp20+l22/2;
        zp2=zp20+l23/2;
        xp3=xp30+l24/2;
        yp3=yp30+l25/2;
        zp3=zp30+l26/2;
        
        calcDerSec();
        
        k31=h*xp2;
        k32=h*yp2;
        k33=h*zp2;
        k34=h*xp3;
        k35=h*yp3;
        k36=h*zp3;
        
        l31=h*xs2;
        l32=h*ys2;
        l33=h*zs2;
        l34=h*xs3;
        l35=h*ys3;
        l36=h*zs3;
        
        T=T0+h;
        
        x2=x20+k31;
        y2=y20+k32;
        z2=z20+k33;
        x3=x30+k34;
        y3=y30+k35;
        z3=z30+k36;
        
        xp2=xp20+l31;
        yp2=yp20+l32;
        zp2=zp20+l33;
        xp3=xp30+l34;
        yp3=yp30+l35;
        zp3=zp30+l36;
        
        calcDerSec();

        k41=h*xp2;
        k42=h*yp2;
        k43=h*zp2;
        k44=h*xp3;
        k45=h*yp3;
        k46=h*zp3;
        
        l41=h*xs2;
        l42=h*ys2;
        l43=h*zs2;
        l44=h*xs3;
        l45=h*ys3;
        l46=h*zs3;
        
        x20+=(k11+2*k21+2*k31+k41)/6;
        y20+=(k12+2*k22+2*k32+k42)/6;
        z20+=(k13+2*k23+2*k33+k43)/6;
        x30+=(k14+2*k24+2*k34+k44)/6;
        y30+=(k15+2*k25+2*k35+k45)/6;
        z30+=(k16+2*k26+2*k36+k46)/6;
        
        xp20+=(l11+2*l21+2*l31+l41)/6;
        yp20+=(l12+2*l22+2*l32+l42)/6;
        zp20+=(l13+2*l23+2*l33+l43)/6;
        xp30+=(l14+2*l24+2*l34+l44)/6;
        yp30+=(l15+2*l25+2*l35+l45)/6;
        zp30+=(l16+2*l26+2*l36+l46)/6;
    }
        
}