/* intégration du mouvement du pendule par Runge-Kutta 4
 * version plate avec une classe Vecteur
 * compiler/executer comme suit:
 *   javac Pendule_rk4_moy
 *   java Pendule_rk4_moy > pendule.data
 * 2008 Adrian Daerr, Univ Paris Diderot
 * domaine public
 */

import static java.lang.Math.*;
import java.text.MessageFormat;

/* Définition d'une classe vecteur à deux composantes qui nous servira
 * à stocker l'état du système.
 */
class Vecteur {
    private final double x, y; // les deux composantes du vecteur
    public Vecteur(double x, double y) { // créateur
        this.x = x;
        this.y = y;
    }

    /** Renvoie un vecteur qui est la somme de ce vecteur et du
     * vecteur 'autre' */
    public Vecteur plus(Vecteur autre) {
        return new Vecteur(x+autre.getX(), y+autre.getY());
    }

    /** Renvoie un vecteur qui est la somme de ce vecteur et d'un
     * scalaire, composante par composante. */
    public Vecteur plus(double scalaire) {
        return new Vecteur(x+scalaire, y+scalaire);
    }

    /** Renvoie un vecteur qui est le produit de ce vecteur et d'un
     * scalaire. */
    public Vecteur fois(double scalaire) {
        return new Vecteur(x*scalaire, y*scalaire);
    }

    /** Renvoie la première composante du Vecteur. */
    public double getX() { return x; }

    /** Renvoie la deuxième composante du Vecteur. */
    public double getY() { return y; }
}

public class Pendule_rk4_moy {
    /* equations du mouvement du pendule rigide */
    static Vecteur deriv(Vecteur etat) {
        return new Vecteur(etat.getY(), -sin(etat.getX()));
    }

    /* boucle principale */
    public static void main(String args[]) {
        final double tfinal=50.0, dt=0.05;
        double t=0.0;
        final double theta0=0.0, omega0=1.99;// conditions initiales
        /* l'état du système est stocké dans un Vecteur;
           convention sur l'ordre des variables d'état: [theta, omega] */
        Vecteur etat = new Vecteur(theta0, omega0);
        Vecteur etat_int, k1, k2, k3, k4;
        int n=0;
        while (t<tfinal) {
            t = n*dt;

            // afficher l'état du système
            System.out.printf("%g %g %g\n",t,etat.getX(),etat.getY());
            // alternative un peu plus propre que:
            //   System.out.println(""+t+" "+etat.getX()+" "+etat.getY());

            // integrer d'un pas de temps
            k1 = deriv(etat);
            etat_int = etat.plus(k1.fois(0.5*dt));
            k2 = deriv(etat_int);
            etat_int = etat.plus(k2.fois(0.5*dt));
            k3 = deriv(etat_int);
            etat_int = etat.plus(k3.fois(dt));
            k4 = deriv(etat_int);
            etat = etat.plus((k1.plus(k2.fois(2)).plus(k3.fois(2)).plus(k4)).fois(dt/6));

            n++;
        }
    }
}