public class Inductive extends Tokamak {

	private double sig0;
	SigT2	sigT2;

	Inductive() 
	{sigT2 = new SigT2();
	sig0 = sigT2.sigv(T0);

	SetConstantInductiveFlux(R0,Bmax0);
	Set_HH_bN_nnGW(HH0, bN0, nnGW0);
	Burntime();
	}
		
	public int Set_HH_bN_nnGW(double HHp, double bNp, double nnGWp)
	{
	//******************************************************************************
	//   HH betaN nnGW solver							
	//******************************************************************************
	bN   = bNp;
	HH   = HHp;
	nnGW = nnGWp;
	
	//-- Normalized inputs --------
	
	double bNc   = bN / bN0;
	double HHc   = HH / HH0;
	double nnGWc = nnGW / nnGW0;
	double Rc    = R/R0;
	double R2c   = R2/R02;
	double ac    = a/a0;
	double a2c   = a2/a02;
	double Bc    = B/B0;
	
	// In normalized units:
	// 1) I = a^2 B / R
	// 2) n = nnGW I/a^2
	// 3) Pa = W^2/(R^2 a)
	// 4) P = (Pa0 Pa + Prf) / (Pa0)
	// 5) W = HH I^0.9 R^1.8 a^0.2 P^0.34 B^0.2 n^0.4
	// 6) W = bN R a I B
	// 7) T = W / (n R2 a)
	
	I = a2c * Bc / Rc;
	n = nnGWc * I / a2c;
	W = bNc * ac * Rc * I * Bc;
	T = W / (n * R2c * ac);
	P = Math.pow(W / (HHc * Math.pow(I, 0.9) * Math.pow(n, 0.4) * Math.pow(Rc, 1.8) * Math.pow(ac, 0.2) * Math.pow(Bc, 0.2)), 2.9412);
	Pbr  = n * n * Math.sqrt(T) * R2c * ac;
	Pece = Math.sqrt(n) * Math.pow(T, 2.5) * Math.pow(Bc, 2.5) * R2c * ac;
	
	double Pa1, Pa2, Pa3;
	
	Pa1 = W * W / (R2c * ac) * sigT2.sigv(T0 * T) / (sig0 * (1. - 2. * fHe0) * (1. - 2. * fHe0));
	Pa2 = 2. * fHe0 * W * ac / (nnGWc * I * R2c * P);
	Pa3 = 2. * Pa1 * Pa2;
	Pa  = Pa1 * 2. * (1. + Pa3 - Math.sqrt(1. + 2. * Pa3))/(Pa3 * Pa3);
	fHe = W * ac / (nnGWc * I * R2c * P) * Pa;
	
	Prf = (P * (Pa0 - Pbr0 - Pece0) - Pa0 * Pa + Pbr0 * Pbr + Pece0 * Pece) / Prf0;
	
	//--- Dimension outputs ----
	Prf *= Prf0;
	I   *= I0;
	n   *= n0;
	T   *= T0;
	fHe *= fHe0;
	Pa  *= Pa0;
	P   *= P0;
	Pbr *= Pbr0;
	Pece*= Pece0;
	W   *= W0;
	
	return(SUCCESS);
	}

	public double RelativeCOE()
	{
	//******************************************************************************
	//   Relative Cost of Electricity							
	//******************************************************************************
	//--- Normalize inputs
	double Rc  = R / R0;
	double ac  = a / a0;
	double Pac = Pa / Pa0;
	double Prfc= Prf / Prf0;
	double Wlc = Pac / (Rc * ac);
	
	double coe = 
		(0.1 * RelativeCapitalCost() +   // Annual capital cost
			0.07 * Rc * ac * WL0 * Wlc / 5.	 // Annual blanket replacement cost
        					// Blanket is 7% of original cost
	        				
			) / (0.1 + 0.07 * WL0 / 5.) 	// Normalization factor
	        
			/ (Pac - 3./5. * Prfc * Prf0 / Pa0) 	    // Electricity production
		
			* (1. + WL0 * Wlc / 15.) / (1. + WL0 / 15.) // Availibility restriction due to blanket change
        							// Replacement time 3 months
	
			* (1. + 500. / burn) / (1. + 500. / burn0); // Availibility restriction due to pulsed operation
        							// Dwell time 500s
	        			
	return(coe);
	}

	public double Burn()    
	{return(burn);}
	
	public double q95() 
	{return(q0 * (R0 * I0) / (a02 * B0) * a2 * B / (R * I));}
	
	private void Burntime()
	{// Normalize
	double fHec = fHe / fHe0;
	double Tc   = T / T0;
	double a2c  = a2 / a02;
	double Ic   = I / I0;
	double Rc   = R / R0;
	
	/*---------------------------
	Flux = (R - a - Shield0 - TF0 - CS0)^2 
	Vl = Zeff I R / (T^3/2 a^2)
	burn = Flux / Vl = T^3/2 a^2 / Zeff
	----------------------------*/
	
	double Fluxc = (R - a - 3.5) / (R0 - a0 - 3.5);
	if(Fluxc < 0.) Fluxc = 0.;
	double Zeffc = (1.3 + 2. * fHe0 * fHec) / (1.3 + 2. * fHe0);
	double Vlc = Zeffc * Ic * Rc / (Tc * Math.sqrt(Tc) * a2c);
	
	burn = burn0 * Fluxc / Vlc;
	}

		
	// ITER Reference Inductive parameters
		
	private static final double HH0  = 1.0;
	private static final double nnGW0 = 1.15;
	private static final double bN0  = 2.2;
		
	private static final double I0   = 21.;
	private static final double B0	  = 5.68;
	private static final double Pa0  = 300.;
	private static final double Prf0 = 100.;
	private static final double q0   = 3.0;
	private static final double WL0  = 1.0;
	private static final double T0   = 12.;
	private static final double Pbr0 = 62.5;
	private static final double Pece0= 42.0;
	private static final double fHe0 = 0.09;
	private static final double P0   = Pa0 + Prf0 - Pbr0 - Pece0;
	private static final double W0   = 1066.;	
	private static final double n0   = 9.80;
	private static final double burn0= 1200.;
		
	// Additional plasma state
		
	private double burn; 
}