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

import java.io.*;

public class TrappedParticles {
    public static void main(String args[]) {
        // 変数宣言
        double x;             // x 座標
        double y;             // y 座標
        double dt;            // 時間刻み幅
        double t;             // 経過時間
        double alpha;         // ピッチ角
        double alpha_eq;      // 赤道面でのピッチ角
        double lamda;         // 粒子位置の緯度
        double v;             // 全速度
        double v_pala;        // 磁力線方向の速度
        double bm;            // ミラーポイントでの磁場の強さ
        double bp;            // 赤道面上惑星表面での磁場の強さ
        double l;             // L 値
        double r_eq;          // 赤道面上での粒子の距離
        double r;             // 粒子の距離
        double sinA_eq;       // sin(alpha_eq)^2
        double sinA;          // sin(alpha)^2
        double pi = Math.PI;

        // 刻み幅
        dt = 0.001;

        // 初期条件
        x = 0.0;
        y = 0.0;
        t = 0.0;
        alpha_eq = pi/9.0; // ピッチ角 30°
        alpha = alpha_eq;
        lamda = 0.0;
        v = 10.0;          // 1 keV の電子の速度は 1.87x10^7 m/s
        bp = 1.0;          // 地球の場合 3.11x10^-5 T
        l = 3.0;           // 惑星半径の 3 倍
        r_eq = 3.0;        // l=r_eq ならば惑星半径が 1 であるということ
        r = r_eq;
        sinA_eq = 0.0;
        sinA = 0.0;

        // (sin(alpha_eq))^2 の値 sinA_eq
        sinA_eq = sin(alpha_eq) * sin(alpha_eq);

        // ミラーポイントでの磁場の強さ bm
        bm = b(bp, l, 1.0, sinA_eq);

        // ファイル出力
        String fn1 = new String();
        // ファイル名
        fn1 = "TrappedParticles.txt";
        for(int i=0; i<300; i++) {
            // (sin(alpha))^2 の値 sinA_eq
            sinA = sinA(sinA_eq, lamda);
            // v_pala の値
            v_pala = v_pala(sinA, v) ;
            // lamda の更新
            lamda = lamda + v_pala*dt/r;
            // r の更新
            r = r_eq * cos(lamda) * cos(lamda);
            // 時間 t の更新
            t = t + dt;
            // 座標変換
            x = r * cos(lamda);
            y = r * sin(lamda);

            // 書き込み
            try {
                PrintWriter pw = new PrintWriter
                    (new BufferedWriter(new FileWriter(fn1,true)));
                if((i%10)==0) {
                    pw.println(t + " " + x + " " + y + " " + v_pala + " " + lamda);
                    pw.close();
                }
            }
            catch(IOException e) {
                System.out.println("入出力エラー");
            }
        }

        // ファイル出力
        String fn2 = new String();
        // ファイル名
        fn2 = "TrappedParticlesTrue.txt";
        lamda = 0.0;
        for(int i=0; i<300; 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 sinA_eq, double lamda) {
        double value = 0;
        value = sinA_eq*sinA_eq
            * Math.pow(1.0+3.0*sin(lamda)*sin(lamda), 0.5)
            / Math.pow(cos(lamda), 6);
        return value;
    }

    // 磁場の強さ b
    private static double b(double bp, double l, double sinA, double sinA_eq) {
        double value = 0;
        value = (bp * sinA * sinA) / (l*l*l * sinA_eq * sinA_eq);
        return value;
    }

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