

/**
Represnts Universe, stores main characteristics:
		<OL>
		<LI>Light speed
		<LI>Reference to central body
		</OL>
*/
class Universe {	
	double c;				//light speed of universe
	CentrBody cb;		//central body of universe

	/**
	Constructs new Universe with given central body and light speed set to 5;
	@param cb  Reference to central body
	*/
	public Universe ( CentrBody cb ) {
		this.cb = cb;
		c = 5;
	}

	/**
	Sets light speed of Universe
	@param c  New light speed
	*/
	public void set_c( double c ) {
		this.c = c;
	}
	
	/**
	Computes Schwartzield radius of central body
	@return Schwartzield radius
	*/
	double getSchwRadius() {
		return 2 * cb.mass / c / c;
	}
}

/**
Stores data about Central Body of Universe
Theese are Mass, Size
*/
class CentrBody {
	double mass; 		//mass of central body
	double size;
	
	/**
	Constructs new Central Body
	@param m  Mass
	@param s  Size
	*/
	public CentrBody( double m, double s ) {
		mass = m;
		size = s;
	}
}

/**
Dependent body
Stores
<OL>
	<LI> position vector in polar
	<LI> radial and angular dynamics in another polar vector
	<LI> reference to Universe
</OL>
*/		
class DependBody {
	Vector pos;			//position vector of this body in POLAR
	Vector dynam;		//dynamics vector
									//mass = 1
	Universe u;			//universe, where body exists
	static double cor_rd = 0.1;					//correction for input of radial dynamics
	static double r_k_correct = 2.80;		//koeficient for Runge-Kutta method (empirical)
	
	/**
	@param u  Belonging universe
	*/
	public DependBody( Universe u ) {
		this.u = u;
	}
	
	/**
	Sets the initial position of body at start of simulation
	@param cp  Cartesian position vector
	@return  Distance from central body
	*/
	public double initPos( Vector cp ) {
		pos = Vector.cart_2_pol( cp );
		//DEBUG System.out.println("r = " + pos.x1 + " fi = " + pos.x2);
		return pos.x1;
	}
	
	/**
	Sets the initial dynamics
	This method is called at start of simulation to set-up dynamics of body
	@param nw  Cartesian vector
	@return  Vector of dynamics in polar. It is computed as vector substraction of <B>nw</B> and current position and some corrections are made.
	*/
	public Vector initDynam( Vector nw ) {
		nw = Vector.cart_2_pol( Vector.minus( nw, Vector.pol_2_cart( pos ) ) );
		double sub_r = nw.x2 - pos.x2;
		dynam = new Vector( nw.x1 * Math.cos(sub_r) * cor_rd, nw.x1 * Math.sin(sub_r) );
		return dynam;
	}
	
	/**
	Updates position of body
	This is main method for simulation It computes new position by using Kanonic form of Hamiltonian equotations and Runge-Kutta method for solving theese diferential equations
	@param time  Time period
	@param relativ  If computing is made ralativistacaly or not
	@return  new position vector in polar
	*/
	public Vector move( double time, boolean relativ ) {
		double k_r[] = new double[4];	
		double k_fi[] = new double[4];
		double k_pr[] = new double[4];		//fields used for Runge-Kutta
		double c, r, pfi, pr, fi;					//aliases for needed variables
		double r2, c2, pfi2, mc, sq;			//frequently used sub-expressions
		
		c = u.c;
		r = pos.x1;
		fi = pos.x2;
		pr = dynam.x1;
		pfi = dynam.x2;
		time /= r_k_correct;
		for( int i = 0; i<4; i++ ) {
			r2 = 1 /r / r;		
			if( relativ ) {			//relativistic 
				c2 = 1 / c / c;
				pfi2 = pfi * pfi;
				mc = 2 * u.cb.mass * c2;
				sq = c / 2 / Math.sqrt( pr * pr + c * c * (1 - mc/r ) * (1 + pfi2 * c2 * r2) );
				k_fi[i] = sq * c * c * ( 1 - mc/r ) * 2 * pfi * c2 * r2;
				k_r[i] = sq * 2 * pr;
				k_pr[i] = -sq * ( mc * r2 * ( c * c + pfi2 * r2 ) - 2 * pfi2 * r2 / r * ( 1 - mc / r ) ); 
			} else {						//nonrelativistic
				k_fi[i] = pfi * r2;
				k_pr[i] = pfi * pfi * r2 / r - u.cb.mass * r2;
				k_r[i] = pr;
				//DEBUG System.out.println("k " + k_r[i] + " " + k_fi[i] + " " + k_pr[i]);
			}	
			if( i==2 )
				time *= r_k_correct;
			if( i!=2 )	{
				r += k_r[i] * time;
				fi += k_fi[i] * time;
				pr += k_pr[i] * time;
				//DEBUG System.out.println("r " + r + " fi " + fi + " pr " + pr );
			}
		}
		//updating of variables
		pos.x1 += ( k_r[0] + r_k_correct * k_r[1] + r_k_correct * k_r[2] + k_r[3]) * time / (2 + 2*r_k_correct);
		pos.x2 += ( k_fi[0] + r_k_correct * k_fi[1] + r_k_correct * k_fi[2] + k_fi[3]) * time / (2 + 2*r_k_correct);
		dynam.x1 += ( k_pr[0] + r_k_correct * k_pr[1] + r_k_correct * k_pr[2] + k_pr[3]) * time / (2 + 2*r_k_correct);		
		
		return pos;
	}
}