import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.MathContext;
import java.math.RoundingMode;
import java.util.ArrayList;

public class Variables {
	public Variables(){
		/*initial conditions*/
		double R=0; //rotor radius
		double theta=0;// Initial pitch angle
		double delta_3=0;// pitch angle flap coupling
		double beta=0;// flap angle (degrees) 
		double beta_p=0;//precone angle (degrees)
		double time=0;
		double timeStep=0;
		double alpha=-90;//initial angle of attack (degrees)
		double RPM=0;
		double theta_tw=0;//Twist Rate (degrees)
		double c =0;//bladeChord
		double dy=	0;//blade element width
		double N_b=0;//number of blades
		double m_total=0;//mass of the model
		double m_seg=0;//mass of a blade element
		double m_blade=0;//mass of the blade 
		double e=0;// the hinge offset
		double k_beta=0;// the flapping stiffness
		double y_n=0;// ?? blade width ??
		double y=0;// ?? blade width ??
		
		/*calculated*/
		double A =0;//Area of the rotor 
		double deltaTheta;// change in pitch generated by pitch-flap coupling
		double r_n;//non-dimensional spanwise location of the nth blade element
		double V_d;//vertical velocity (m/s)
		double gravity = 9.81;//(m/s)
		double rho = 1.2041;//Atmospheric pressure (kg/m^3)
		double U_pn;//out-of-plane velocity for blade element n (m/s)
		double W;//model weight
				
		
		/*calculating initial conditions*/
		r_n = y_n/R;
		/*deltaTheta = -Math.tan(delta_3)*0;//no change in beta*/
		theta = theta + r_n * theta_tw;// + deltaTheta;
		U_pn = V_d = gravity*timeStep;
		W = gravity * m_total;
		
		/*calculating aerodynamic forces*/
		double U_tn;//in-plane velocity for blade element n (m/s)
		double v_i;// Induced velocity
		double sigma;//rotor solidity
		double phi;//inflow angle
		double Usqrd_n;//resultant velocity
		double CL_n;//Coefficient of lift
		double CD_n;//Coefficient of drag 
		double a;//lift curve slope (for 2-D blade section)
		double dL_n;//Elemental lift on segment n
		double dD_n;//Elemental drag on segment n
		double dF_xn;//in-plane component (double check)
		double dF_zn;//out-of-plane component  (double check)
		double dT_n;// Elemental Thrust on an annulus
		double dQ_n;// Elemental Torque on an annulus
		double A_annul;//area of the annulus
		double v_hn;// induced velocity in hover 
		double m; //mass flow through the whole rotor disk
		double T =1;//rotor thrust (initial value irrelevant: used to ensure loop) 
		double w;//slipstream velocity
		double v_in;//induced velocity over one annul
		double M_af;//Moment created by the AF (Nm)
		double M_cf;//Moment created by the CF (Nm)
		double Q;// rotor torque
		double I_rotor;// rotor moment of inertia
		double a_body;//vertical acceleration 
		double deltaBeta;//change in beta over one step 

		
		while(T>=0.001){
			U_tn = RPM *r_n;
			sigma = (N_b *c)/(Math.PI*R);
			a=0;//hopefully in 2-d airfoil lookup tables 
			v_i = (sigma*a/16)*(Math.sqrt(1+((32/sigma*a)*theta*r_n))-1);
			U_pn = v_i - V_d;
			phi = Math.atan(U_pn/U_tn);
			alpha = delta_3 - phi;//needed twice?
			CL_n=0;//look up from 2-d airfoil lookup tables
			CD_n=0;//look up from 2-d airfoil lookup tables
			Usqrd_n = Math.pow(U_pn,2)+Math.pow(U_tn,2);
			dL_n = .5*rho*Usqrd_n*c*CL_n*dy;
			dD_n = .5*rho*Usqrd_n*c*CD_n*dy;
			dF_xn = dL_n*Math.sin(phi)+dD_n*Math.cos(phi);
			dF_zn = dL_n*Math.cos(phi)-dD_n*Math.sin(phi);
			dT_n = N_b*dF_zn;
			dQ_n = N_b*dF_xn*y_n;
			A_annul = 2*Math.PI*y_n*dy*Math.cos(beta);
			v_hn = Math.sqrt(dT_n/(2*rho*A_annul));
			//m = rho*A*(v_i-V_d);//needed?
			//w = 2*v_i;//needed?
			//T = -m*w;//needed?
			v_in = (V_d/2) - Math.sqrt(Math.pow(V_d/2,2)-Math.pow(v_hn,2));
			alpha = delta_3 - phi;//needed twice?
			M_af = N_b *(dF_zn*(y-e)*Math.cos(beta));//Summation of all blade elements
			M_cf = N_b * (m_seg*Math.pow(RPM,2)*(y_n - e)*Math.sin(beta)*(e+(y_n - e)*Math.cos(beta)));//Summation of all blade elements
			deltaBeta = beta;//store value for beta at current time step
			if(beta >= 0){ 
				beta = ((M_af - M_cf)/k_beta) + beta_p;
			}else{
				beta = ((M_af + M_cf)/k_beta) + beta_p;
			}
			deltaBeta = Math.abs(deltaBeta-beta);//calculate change in beta after update
			//update in & out of plane loads form blade pitch and alpha ??
			T = N_b * dT_n;
			Q = N_b * dQ_n;
			I_rotor = (N_b/3)*m_blade*Math.pow(R-e,2);
			RPM = N_b * Q * (timeStep/I_rotor) + RPM;
			a_body = (T-W)/m_total;
			V_d = V_d + timeStep * a_body;
			deltaTheta = -Math.tan(delta_3)*deltaBeta;//beta = change in beta ??
			theta = theta + r_n * theta_tw + deltaTheta;
			
		}
		//stable auto rotation is achieved at this point 
	}
	
	
	/*updates conditions
	public Variables(Variables v){
		
	}
	public double mu(){
		return (0.000001458 * Math.sqrt(Temp))/(1+(110.4/Temp));
	}
	public double Temperature(){
		if(A < 11){
			return Temp_0 + (-6.5)*(A - 0);
		}else if(11 <= A && A < 20.1){
			return Temp_0 + (0)*(A - 11);
		}else if(20.1 <= A && A < 32.2){
			return Temp_0 + (1)*(A - 20.1);
		}else if(32.2 <= A && A < 47.3){
			return Temp_0 + (2.8)*(A - 32.2);
		}else if(47.3 <= A && A < 52.4){
			return Temp_0 + (0)*(A - 47.3);
		}else if(52.4 <= A && A < 61.6){
			return Temp_0 + (-2)*(A - 52.4);
		}else if(61.6 <= A && A < 80){
			return Temp_0 + (-4)*(A - 61.6);
		}else{
			return Temp_0 + (0)*(A - 80);
		}
	}
	public double Airflow(){
		return omega* r;
	}
	public double RN(){
		return (vAirflow*rho*l)/mu;
	}
	public double PitchAttitudeAngle(){
		return Math.atan(RD/(omega*r))-alpha;
	}
	public double Weight(){
		return m*g;
	}
	public double trueVelocity(){
		return Math.sqrt((2*W)/(rho*CL*S));
	}
	public double terminalVelocity(){
		return Math.sqrt((2*W)/(rho*CD*S));
	}
	public double RD(){
		return ((vTrue*CD)/CL);
	}
	public double b(){
		return (.5*rho*CD*S);
	}
	public double gravity(){
		return ((-0.000003)*A)+ 9.8064;
	}
	public double airDensity(){
		return 1.4213*Math.exp((-0.0001)*A);
	}
	public double getAltitude(){
		return A;
	}
	public String Headers(){
		return "Altitude,Time,Terminal Velocity,Vertical Velocity,True Velocity,Pitch Attitude Angle,RN";
	}
	public String toString(){
		return A+","+t+","+vTerminal+","+vVerticle+","+vTrue+","+theta+","+RN;
	}
	//hyperbolic cos
	/*from: http://www.shodor.org/~jmorrell/interactivate/org/shodor/math11/FormulaMath.java
	public static BigDecimal cosh(BigDecimal x) {
		return (exp(x).add(x.negate())).divide(new BigDecimal("2"), x.scale(),
				BigDecimal.ROUND_HALF_EVEN);
	}
	public static BigDecimal exp(BigDecimal bd) {
		double d = bd.doubleValue();
		BigDecimal ret = new BigDecimal(Math.exp(d));
		return ret.setScale(bd.scale(), BigDecimal.ROUND_HALF_EVEN);
	}
	public static BigDecimal sqrt(BigDecimal bd) {
		double d = bd.doubleValue();
		BigDecimal ret = new BigDecimal(Math.sqrt(d));
		return ret.setScale(bd.scale(), BigDecimal.ROUND_HALF_EVEN);
	}
	/**
	 * COMMON log.  (Log base 10)
	 * @param bd
	 * @return Common log of bd
	 
	public static BigDecimal log10(BigDecimal x) {
		BigDecimal ret = log(x);
		BigDecimal denom = log(new BigDecimal("10"));
		BigDecimal result = ret.divide(denom, x.scale(),
				BigDecimal.ROUND_HALF_EVEN);
		return result;
	}
	/**
	 * NATURAL log.  (Log base e)
	 * @param bd
	 * @return Natural log of bd
	 
	public static BigDecimal log(BigDecimal bd) {
		double d = bd.doubleValue();
		double result = 0.0;
		result = Math.log(d);
		if(Double.isNaN(result) || Double.isInfinite(result))
			return new BigDecimal("0");
		BigDecimal ret = new BigDecimal(Math.log(d));
		return ret.setScale(bd.scale(), BigDecimal.ROUND_HALF_EVEN);
	}*/

}