/*
  2003/10/25  小松 研吾
*/

import java.io.*;

public class TrappedParticlesM {
    public static double pi = Math.PI;  // 円周率 π
    public static double dt;            // 時間刻み幅
    public static double r_eq;          // 赤道面上での粒子の距離
    public static double l;             // L 値
    public static double alpha_eq;      // 赤道面でのピッチ角
    public static double sinA_eq;       // sin(alpha_eq)^2
    public static double bp;            // 赤道面上惑星表面での磁場の強さ
    public static double v;             // 全速度

    public static void main(String args[]) {
        // 変数宣言
        double x;             // x 座標
        double y;             // y 座標
        double t;             // 経過時間
        double alpha;         // ピッチ角
        double lamda;         // 粒子位置の緯度
        double v_pala;        // 磁力線方向の速度
        double r;             // 粒子の距離
        double sinA;          // sin(alpha)^2
        double b;             // 磁場の強さ
        double bx;            // x 方向の磁場の強さ
        double by;            // y 方向の磁場の強さ
        double bm;            // ミラーポイントでの磁場の強さ
        double cosLamda;      // cos(lamda)
        double sinLamda;      // sin(lamda)

        // 初期化
        dt       = 0.0;
        r_eq     = 0.0;
        l        = 0.0;
        alpha_eq = 0.0;
        sinA_eq  = 0.0;
        bp       = 0.0;
        v        = 0.0;

        x        = 0.0;
        y        = 0.0;
        t        = 0.0;
        alpha    = 0.0;
        lamda    = 0.0;
        v_pala   = 0.0;
        r        = 0.0;
        sinA     = 0.0;
        b        = 0.0;
        bx       = 0.0;
        by       = 0.0;
        bm       = 0.0;
        cosLamda = 0.0;
        sinLamda = 0.0;


        // 刻み幅
        dt = 0.001;

        // 初期状態
        r_eq = 3.0;           // l=r_eq ならば惑星半径が 1 であるということ
        l = 3.0;              // 惑星半径の 3 倍
        alpha_eq = pi/6.0;    // ピッチ角 30°
        sinA_eq = sin(alpha_eq) * sin(alpha_eq);   // (sin(alpha_eq))^2 の値 sinA_eq
        bp = -1.0;            // 地球の場合 3.11x10^-5 T 。
                              // 負号は双極子モーメントの向きが南向きであることを示す。
        v = 10.0;             // 1 keV の電子の速度は 1.87x10^7 m/s
        x = l;                // 赤道面での x
        y = 0.0;              // 赤道面での y
        alpha = alpha_eq;     // 赤道面でのピッチ角
        lamda = 0.0;          // 赤道面での λ
        r = r_eq;             // 赤道面での r
        sinA = sinA_eq;       // 赤道面での (sin(alpha))^2
        cosLamda = 1.0;
        sinLamda = 0.0;


        // ミラーポイントでの磁場の強さ bm
        //        bm = b(x, y, r);

        // ファイル出力
        String fn1 = new String();
        // ファイル名
        fn1 = "TrappedParticlesM.txt";
        for(int i=0; i<1000; i++) {
            // 現在の位置の距離
            r = r(x, y);
            // 現在の位置での磁場の強さ
            b = b(x, y, r);
            bx = bx(x, y, r);
            by = by(x, y, r);
            // (sin(alpha))^2 の値 sinA_eq
            sinA = sinA(b);
            // v_pala の値
            v_pala = v_pala(sinA) ;
            // x, y の更新
            x += v_pala * dt * bx/b;
            y += v_pala * dt * by/b;
            // 時間 t の更新
            t += dt;
            // λを求める
            cosLamda = x / r;
            sinLamda = y / r;
            // 座標変換
            //            x = r * cos(lamda);
            //            y = r * sin(lamda);

            // 書き込み
            try {
                PrintWriter pw = new PrintWriter
                    (new BufferedWriter(new FileWriter(fn1,true)));
                if((i%20)==0) {   // 20 ステップごとに記録
                    pw.println(t + " " + x + " " + y + " " + v_pala
                               + " " + bx/b + " " + by/b + " " + cosLamda);
                    pw.close();
                }
            }
            catch(IOException e) {
                System.out.println("入出力エラー");
            }
        }

        // ファイル出力
        String fn2 = new String();
        // ファイル名
        fn2 = "TrappedParticlesTrue.txt";
        lamda = 0.0;
        for(int i=0; i<150; i++) {
            r = r_eq * cos(lamda) * cos(lamda);
            // 座標変換
            x = r * cos(lamda);
            y = r * sin(lamda);
            // lamda の更新
            lamda = lamda + 0.01;

            // 書き込み
            try {
                PrintWriter pw = new PrintWriter
                    (new BufferedWriter(new FileWriter(fn2,true)));
                pw.println(x + " " + y);
                pw.close();
            }
            catch(IOException e) {
                System.out.println("入出力エラー");
            }
        }
    }

    // sin(theta)
    private static double sin(double theta) {
        double value = 0;
        value = Math.sin(theta);
        return value;
    }

    // cos(theta)
    private static double cos(double theta) {
        double value = 0;
        value = Math.cos(theta);
        return value;
    }

    // (sin(alpha))^2 の値 sinA
    private static double sinA(double b) {
        double value = 0;
        value = sinA_eq * l*l*l * b / Math.abs(bp);
        return value;
    }

    // (x*x+y*y)^(1/2) の値 r
    private static double r(double x, double y) {
        double value = 0;
        value = Math.pow((x*x + y*y), 0.5);
        return value;
    }

    // 磁場の強さ b
    private static double b(double x, double y, double r) {
        double value = 0;
        value = Math.abs(bp) * Math.pow(r_eq, 3) * Math.pow((4.0*y*y + x*x), 0.5)
            / (l*l*l * Math.pow(r, 4));
        return value;
    }

    // 磁場の x 方向の強さ bx
    private static double bx(double x, double y, double r) {
        double value = 0;
        value = bp * Math.pow(r_eq, 3) * 3.0*x*y / (l*l*l * Math.pow(r, 5));
        return value;
    }

    // 磁場の y 方向の強さ by
    private static double by(double x, double y, double r) {
        double value = 0;
        value = bp * Math.pow(r_eq, 3) * (2.0*y*y - x*x) / (l*l*l * Math.pow(r, 5));
        return value;
    }

    // 磁場に平行な速度 v_pala
    private static double v_pala(double sinA) {
        double value = 0;
        value = Math.pow(1.0-sinA, 0.5) * v;
        return value;
    }
}