class BlackholeCalculator{

double dr;                                                                                //変数、定数の宣言
double x,y;


double EPS_TURN=0.0000001;
double EPS_BH=0.00001;
double t=0;

double x0=10.0;                                                                         //初期状態
double xold=x0;                                                                          //一つ前のx1
double x1=x0;                                                                               //x1の初期状態
double x2=Math.asin(b/x0);                                                            //x2の初期状態
  

double fugo=-1.0;
double b=5.2;                                                                                    //衝突係数


public BlackholeCalculator(){
    


}




private double diff_eq1(double t, double x1, double x2,double fugo,double b){
    return fugo*(1.0-(1.0/this.x1))*Math.sqrt(1.0-((this.b*this.b)/(this.x1*this.x1)) +((this.b*this.b)/(this.x1*this.x1*this.x1)));}　　　　　//微分方程式1

private double diff_eq2(double t, double x1, double x2,double b){                                                                                                     //微分方程式2
    return (1.0-(1.0/this.x1))*((this.b/(this.x1*this.x1)));    }



public void step(){

    double dt=0.1*Math.sqrt(((x1*x1*x1)-(b*b)*(x1-1.0))/(x1*x1*x1));                                  //刻み幅dt

    
    double k1=diff_eq1(t,x1,x2,fugo,b);                                                                                           //ルンゲ・クッタ法
    double l1=diff_eq2(t,x1,x2,b);
    double k2=diff_eq1(t+0.5*dt,x1+0.5*k1*dt,x2+0.5*l1*dt,fugo,b);
    double l2=diff_eq2(t+0.5*dt,x1+0.5*k1*dt,x2+0.5*l1*dt,b);
    double k3=diff_eq1(t+0.5*dt,x1+0.5*k2*dt,x2+0.5*l2*dt,fugo,b);
    double l3=diff_eq2(t+0.5*dt,x1+0.5*k2*dt,x2+0.5*l2*dt,b);
    double k4=diff_eq1(t+dt,x1+k3*dt,x2+l3*dt,fugo,b);
    double l4=diff_eq2(t+dt,x1+k3*dt,x2+l3*dt,b);

    x1=x1+(1/6.0)*(k1+2.0*k2+2.0*k3+k4)*dt;
    x2=x2+(1/6.0)*(l1+2.0*l2+2.0*l3+l4)*dt;



    this.dr=this.x1-this.xold;　　　　　　　　　　　　　　　　　　　　　　　　　　//現在のｘ１と一つ前のｘ１の差分
    if(Math.abs(this.dr)<=this.EPS_TURN){this.fugo=1.0;}　　　　　　　　　//もしdrの絶対値がある値以下であればfugoを変える
    
    this.xold=this.x1;                                                                                        //xoldの更新
    
    this.x=x1*Math.cos(x2);　　　　　　　　　　　　　　　　　　　　　　　　　　　　   //直角座標に変換
    this.y=x1*Math.sin(x2);       

}

public double gett(){　　　　　　　　　　　　　　//各変数を取得するメソッド
    return this.t;
}


public double getx(){
    return this.x;
}

public double gety(){
    return this.y;
}

public double getx1(){
    return this.x1;

}

}