public abstract class Tokamak {
	Tokamak() {}
		
	// Implented by sub-class: inductive or reversed
	public abstract int Set_HH_bN_nnGW( double HH,  double bN,  double nnGW);

		
	public int SetConstantAspectRatio( double Rp,  double Bmaxp)
	{
	R = Rp;
	a = a0R0 * Rp;
	Bmax = Bmaxp;

	B = Bmax * (1 - a0R0 - Shield0 / R);

	R2 = R * R;
	a2 = a * a;
	return(SUCCESS);
	}

	public int SetConstantInductiveFlux( double Rp,  double Bmaxp)
	{ 
	R = Rp;
	a = 0.632 * R - 2.34;
	Bmax = Bmaxp;

	B = Bmax * (1 - a0R0 - Shield0 / R);

	R2 = R * R;
	a2 = a * a;

	return(SUCCESS);
	}

	public int Set_HH_Paux_nnGW( double HHp,  double Prfp,  double nnGWp)
	{
	//*******************************************************************************
	//*   HH Paux nnGW								
	//*										
	//*   Finds bN such that: (HH bN nnGW) => Paux					
	//*										
	//*******************************************************************************
	double bNtry = 2.0;

	double step  = 0.1 * bNtry;
	Set_HH_bN_nnGW(HHp, bNtry, nnGWp);
	if(Paux() > Prfp)
		{// Decrease bN
		do 
			{
			if(Paux() > Prfp)
				bNtry -= step;
			else
				{bNtry += step;
				step *= 0.5;
				bNtry -= step;
				}
			Set_HH_bN_nnGW(HHp, bNtry, nnGWp);
			}
		while(step > 1e-6 && Math.abs(Paux() - Prfp) > 1e-4);
		}
	else
		{// Increase bN
		do 
			{
			if(Paux() < Prfp)
				bNtry += step;
			else
				{bNtry -= step;
				step *= 0.5;
				bNtry += step;
				}
			Set_HH_bN_nnGW(HHp, bNtry, nnGWp);
			}
		while(step > 1e-6 && Math.abs(Paux() - Prfp) > 1e-4);
		}
	return(SUCCESS);
	}	

	public int Set_HH_Paux_Palpha( double HHp,  double Prfp,  double Pap)
	{
	//*******************************************************************************
	//   HH Paux Palpha								
	//										
	//   Finds nnGW such that: (HH Paux nnGW) => Palpha				
	//										
	//******************************************************************************
	double nGWtry = 1.0;
	double step = 0.1 * nGWtry;
	Set_HH_Paux_nnGW(HHp, Prfp, nGWtry);
	if (Palpha() > Pap) 
		{// Decrease nGW
		do
			{if(Palpha() > Pap)
				nGWtry -= step;
			else
				{nGWtry += step;
				step *= 0.5;
				nGWtry -= step;
				}
 			Set_HH_Paux_nnGW(HHp, Prfp, nGWtry);
			}
		while(step > 1e-6 && Math.abs(Palpha() - Pap) > 1e-4);
		}
	else
		{// Increase nGW
		double check;
		do
			{if(Palpha() < Pap)
				{nGWtry += step;
				check = Palpha();
				}
			else
				{nGWtry -= step;
				step *= 0.5;
				nGWtry += step;
				check = -1.;
				}
 			Set_HH_Paux_nnGW(HHp, Prfp, nGWtry);
	 		if(check > 0. && check >= Palpha()) return(FAILURE);
			}
		while(step > 1e-6 && Math.abs(Palpha() - Pap) > 1e-4);
		}
	
	return(SUCCESS);
	}

	public int Set_HH_Palpha( double HHp,  double Pap)
	{
	//******************************************************************************
	//   HH Palpha									
	//										
	//   Finds Paux such that: (HH Paux Palpha) succeeds				
	//										
	//******************************************************************************
	
	double PauxTry = 100.;
	double Px = PauxTry;
	
	// Increases Paux until solution is found
	
	while(Set_HH_Paux_Palpha(HHp, Px, Pap) == FAILURE)
		Px += 0.1 * PauxTry;
		
	// Decrease to within 1MW
	
	double step = 0.1 * PauxTry; 
	do
		{if(Set_HH_Paux_Palpha(HHp, Px, Pap) == SUCCESS)
			{if(Px <= 0.) 
				{Px = 0.;
				Set_HH_Paux_Palpha(HHp, Px, Pap);
				break;
				}
			Px -= step;
			}
		else
			{Px += step;
			step *= 0.5;
			Px -= step;
			}
		}
	while(step > 1e-2 * PauxTry);
	
	if(Px > 0.)
		{// Ended on a failure
		Px += step;
		Set_HH_Paux_Palpha(HHp, Px, Pap);
		}
	return(SUCCESS);	
	}
	
	public double RelativeCapitalCost()
	{ //--- Normalized paramters
	double Rc, ac, Ic, Prfc, Pfc; 
	Rc   = R / R0;
	ac   = a / a0;
	Ic   = I / Iplasma0;
	Prfc = Prf / Paux0;
	if(Prfc < Rc * ac) Prfc = Rc * ac; // Minimum installed power
	Pfc  = 5. * Pa / Pfusion0;

	double shield    = 150. * Rc * ac;
	double firstWall = 80. * Rc * ac;
	double Blanket   = 300. * Rc * ac;
	double fuelHandling = 80. * Pfc;
	double divertor  = 180. * Rc * ac;
	double PFcoil    = 580. * Rc * Ic;
	double TFcoil    = 900. * Rc * (a0 * ac + Shield0) / (a0 + Shield0);
	double Padd      = 250. * Prfc;
	double Remote    = 250. * (R0 * Rc + a0 * ac + Shield0) / (R0 + a0 + Shield0);
	double HeatSink  = 400. * Pfc;
	double Power	  = 300. * (Prfc + Ic) / 2.;
	double Vacuum	  = 50. * Rc * ac;
	double Building  = 850. * (2. + Rc * (a0 * ac + Shield0) / (a0 + Shield0))/(2. + 1.);
	 
	double total = shield + firstWall + Blanket + fuelHandling + divertor + PFcoil + 
	 		TFcoil + Padd + Remote + HeatSink + Power + Vacuum + Building;
		 
	total /= 150.+80.+300.+80.+180.+580.+900.+250.+250.+400.+300.+50.+850.;
	// 4370
	return(total);
	}

	public double RelativeWallLoading(){
		 return(a0 * R0 * 5. * Pa / (a * R * Pfusion0));}
	public double RelativePaux(){
		 return(a0 * R0 * Prf / (a * R * Paux0));}
	abstract public  double q95();
	abstract public double RelativeCOE();
		
	public double Pfusion() {return(5. * Pa);}
	public double NnGW()    {return(nnGW);}
	public double Palpha()  {return(Pa);}
	public double Paux()    {return(Prf);}
	public double BetaN()   {return(bN);}
	public double Qth()     {return((Prf <= 0.) ? 1000. : 5. * Pa / Prf);}

	public static final int SUCCESS = 0;
	public static final int FAILURE = 1;
		
	
	// ITER constants
		
	// Radial build:
	// |---|CS|TF|SHIELD| |-- plasma
	// 0  1.8 | 4.0    5.2
	// |   | 2.8 |      | 5.34
	
	protected static final  double	a0R0 	= 0.344;// Aspect ratio
	protected static final  double 	Shield0 = 1.64;	// Shield thickness
	protected static final  double 	TF0     = 1.20;	// TF thickness
	protected static final  double 	CS0     = 1.00;	// CS thickness		static final  double 	R0    	= 8.14;	// Major radius	
	protected static final  double	R0		= 8.14; // Major radius
	protected static final  double 	a0	= 2.80; // Minor radius
	protected static final  double 	Bmax0 	= 12.5;	// Max. field of the TF coil
	protected static final  double 	R02 = R0 * R0;
	protected static final  double 	a02 = a0 * a0;
	protected static final  double	Pfusion0 = 1500.; // Nominal fusion power    (for costing)
	protected static final  double	Paux0	 = 100.;  // Nominal auxiliary power (for costing)
	protected static final  double	Iplasma0 = 21.;   // Nominal plasma current  (for costing)
	
	//NON normalized plasma state:
		
		//Geometry
	protected double R, a, R2, a2;
	protected double Bmax;
	protected double B;
		
	//Plasma
	protected double HH;
	protected double Prf;
	protected double nnGW;
	protected double I;
	protected double n;
	protected double T;
	protected double fHe;
	protected double Pa;
	protected double P;
	protected double Pbr;
	protected double Pece;
	protected double W;
	protected double bN;
	
}