content: package edu.stanford.rsl.conrad.physics;
  
  import java.text.NumberFormat;
  
  import java.util.Arrays;
  
  import edu.stanford.rsl.conrad.utils.CONRAD;
  import edu.stanford.rsl.conrad.utils.DoubleArrayUtil;
  import edu.stanford.rsl.conrad.utils.LinearInterpolatingDoubleArray;
  import static edu.stanford.rsl.conrad.physics.Constants.*;
  import static edu.stanford.rsl.conrad.utils.DoubleArrayUtil.*;
  
  public class XRaySpectrum {
  
  	/**
  	 * Creates a spectrum that is similar to a C-arm spectrum with default parameters.
  	 * @param E the energies in [keV]
  	 * @param kVp the peak voltage in [kV]
  	 * @param target the target material ("W" or "Mo")
  	 * @param mAs the acceleration time times the current 
  	 * @return the photon flux in [photons/mm2/bin] in an array matching the energies E
  	 * @throws Exception
  	 */
  	public static double [] generateXRaySpectrum(double [] E, double kVp, String target, double mAs) throws Exception{
  		double degrees, mmpyrex, mmoil, mmlexan, mmAl, mdis;
  		if (target.equals("W")) {
  				degrees = 10;
  				mmpyrex = 2.38;
  				mmoil = 3.06;
  				mmlexan = 2.66;
  				mmAl = 1.50;
  				mdis = 1;
  		}else if (target.equals("Mo")) {
  				degrees = 13;
  				mmpyrex = 1.11;
  				mmoil = 3.06;
  				mmlexan = 2.42;
  				mmAl = 0;
  				mdis = 1;
  		} else {
  			throw new Exception("Undefined target material");
  		}
  		return generateXRaySpectrum(E, kVp, target, mAs, mdis, degrees, mmpyrex, mmoil, mmlexan, mmAl);
  	}
  	
  	/**
  	 * 
  	 * @param E the energies in [keV]
  	 * @param kVp the peak voltage in [kV]
  	 * @param target the target material ("W" or "Mo")
  	 * @param mAs the acceleration time times the current 
  	 * @param mdis amount of air in [m]
  	 * @param degrees tube angle in [deg]
  	 * @param mmpyrex amount of pyrex filtration in [mm]
  	 * @param mmoil amount of oil filtration in [mm]
  	 * @param mmlexan amount of lexan filtration in [mm]
  	 * @param mmAl amount of Al filtration in [mm]
  	 * @return the photon flux in [photons/mm2/bin] in an array matching the energies E
  	 * @throws Exception
  	 */
  	public static double [] generateXRaySpectrum(double [] E, double kVp, String target, double mAs, double mdis, double degrees, double mmpyrex, double mmoil, double mmlexan, double mmAl) throws Exception{
  		double Emin, Emax;
  		double Emin_W	=	10.0;		/* lower keV cutoff for W spectra */
  		double Emin_Mo	=	3.0;		/* lower keV cutoff for Mo spectra */
  		double Emax_W	=	301.;		/* upper keV cutoff for W spectra */
  		double Emax_Mo	=	101.;		/* upper keV cutoff for Mo spectra */
  		if (target.equals("W")) {
  			Emin = Emin_W;
  			Emax = Emax_W;
  		}else if (target.equals("Mo")) {
  			Emin = Emin_Mo;
  			Emax = Emax_Mo;
  		} else {
  			throw new Exception("Undefined target material");
  		}
  		double radians = degrees / 180.0 * Math.PI;
  		double dE = E[1]-E[0];
  		if (kVp+dE > Emax)
  			throw new Exception("Desired kVp greater than Emax");
  		if (E[0]<Emin)
  			throw new Exception("E contains element less than Emin for this target material");
  		double [] phi = new double [E.length];
  		NumberFormat nf = NumberFormat.getInstance();
  		nf.setMaximumFractionDigits(10);
  		nf.setMaximumIntegerDigits(10);
  		DoubleArrayUtil.add(phi, characteristic(E, target, kVp, radians));
  		
  		DoubleArrayUtil.add(phi, bremsstrahlung(E, target, kVp, radians));
  		
  		/* Calculated spectrum transmitted through the inherent filter materials */
  		double [] inherentFilter = DoubleArrayUtil.multiply(mu(E, "pyrex"), PYREX_mdensity*mmpyrex/10.0);
  		DoubleArrayUtil.add(inherentFilter, DoubleArrayUtil.multiply(mu(E, "lexan") , LEXAN_mdensity*mmlexan/10.0));
  		
  		DoubleArrayUtil.add(inherentFilter, DoubleArrayUtil.multiply(mu(E, "Al") , Al_mdensity*mmAl/10.0));
  		DoubleArrayUtil.add(inherentFilter, DoubleArrayUtil.multiply(mu(E, "air"), AIR_mdensity*mdis*100.0));
  		
  		DoubleArrayUtil.multiply(inherentFilter, -1);
  		DoubleArrayUtil.exp(inherentFilter);
  		DoubleArrayUtil.multiply(phi, inherentFilter);
  		
  		/* Now scale to photons/mm2/bin.*/
  		double f = mAs * EpMAS / (4.0 * Math.PI * mdis * mdis * 1.0E6);
  		DoubleArrayUtil.multiply(phi, f);
  		return phi;
  	}
  
  	/**
  	 * Returns the mass attenuation coef for the selected material at
  	 * energy E.  The coefs are fit to Howerton data by N. Cardinal,
  	 * and are thought to be significantly better than the coefs used by Tucker.
  	 * Except for Re, which gives nonsense results, so use Tucker for Re
  	 * @param E the array of energies
  	 * @param material the material ("W", "Re", "Mo", "Al", "oil", "lexan", "pyrex", or "air")
  	 * @return the array of absorption coefficients.
  	 * @throws Exception may occur
  	 */
  	public static double [] mu(double [] E, String material) throws Exception{
  		double [] uT = Arrays.copyOf(E, E.length);
  		DoubleArrayUtil.divide(uT, 100);
  		double [] mu = new double[E.length];
  		int nEbin = E.length;
  		double [] a = null;
  		for (int j=0; j<nEbin; j++){
  			mu[j]=0;
  			if (material.equals("Re")){
  				if (E[j]>= Re_K){
  					a = new double [] {-5.803e-1, 3.336e0, 3.292e0, 1.502e0, 0};
  				} else {
  					a= new double [] {-2.987e-1, 5.815e-1, 6.230e-1, 7.727e-2, -6.155e-3};
  				}
  				mu[j] += a[0];
  				mu[j] += a[1] / Math.pow(uT[j], 1.6);
  				mu[j] += a[2] / Math.pow(uT[j], 2.7);
  				mu[j] += a[3] / Math.pow(uT[j], 3.5);
  				mu[j] += a[4] / Math.pow(uT[j], 4.5);
  			} else {
  				if (material.equals("W")){
  					if (E[j]> W_K){
  						a = aWh;
  					} else {
  						a = aWl;
  					}
  				} else if (material.equals("Re")){
  					if (E[j]> Re_K){
  						a = aReh;
  					} else {
  						a = aRel;
  					}
  				} else if (material.equals("Mo")){
  					if (E[j]> Mo_K){
  						a = aMoh;
  					} else {
  						a = aMol;
  					}
  				} else if (material.equals("Al")){
  					a = aAl;
  				} else if (material.equals("oil")){
  					a = aOIL;
  				} else if (material.equals("lexan")){
  					a = aLEX;
  				} else if (material.equals("pyrex")){
  					a = aPYR;
  				} else if (material.equals("air")){
  					a = aAIR;
  				} else {
  					throw new Exception("Undefined material in mu()");
  				}
  				
  				double u = Math.pow(100.0/E[j], 0.5);
  				mu[j] = a[0] * u;
  				
  				for(int n = 1; n < a.length; n++){
  					mu[j] +=  (a[n]) * Math.pow(u, n+1);
  					
  				}
  				//print("a", a);
  				//System.out.println(j+ " " + mu[j]);
  			}
  		}
  		
  		return mu;
  	}
  
  	/**
  	 * Generates characteristic lines and adds them to the fluence spectrum.  The
  	 * characteristic peaks are added to the nearest two bins, weighted by the
  	 * position of the peak between the bins
  	 * @param E the array of energies
  	 * @param tubeTarget the target material
  	 * @param To Peak kilo volage
  	 * @param theta angle
  	 * @return the array of characteristics.
  	 * @throws Exception may occur
  	 */
  	public static double [] characteristic (double [] E, String tubeTarget, double To, double theta) throws Exception{
  		double [] phi = new double [E.length];
  		if (tubeTarget.equals("W")){
  			if (To > Re_K) { /* add Rhenium characteristic lines*/
  				DoubleArrayUtil.add(phi, characteristicLine(E, tubeTarget, To, Re_K, Re_Ka1_f, Re_Ka1, theta));
  				DoubleArrayUtil.add(phi, characteristicLine(E, tubeTarget, To, Re_K, Re_Ka2_f, Re_Ka2, theta));
  				DoubleArrayUtil.add(phi, characteristicLine(E, tubeTarget, To, Re_K, Re_Kb1_f, Re_Kb1, theta));
  				DoubleArrayUtil.add(phi, characteristicLine(E, tubeTarget, To, Re_K, Re_Kb2_f, Re_Kb2, theta));
  			}
  			if (To > W_K) { /* Tungsten characteristic lines */
  				DoubleArrayUtil.add(phi, characteristicLine(E, tubeTarget, To, W_K, W_Ka1_f, W_Ka1, theta));
  				DoubleArrayUtil.add(phi, characteristicLine(E, tubeTarget, To, W_K, W_Ka2_f, W_Ka2, theta));
  				DoubleArrayUtil.add(phi, characteristicLine(E, tubeTarget, To, W_K, W_Kb1_f, W_Kb1, theta));
  				DoubleArrayUtil.add(phi, characteristicLine(E, tubeTarget, To, W_K, W_Kb2_f, W_Kb2, theta));
  			}
  		} else if (tubeTarget.equals("Mo")){ /* Molybdenum characteristic line */
  			if (To > Mo_K){
  				DoubleArrayUtil.add(phi, characteristicLine(E, tubeTarget, To, Mo_K, Mo_Ka1_f, Mo_Ka1, theta));
  				DoubleArrayUtil.add(phi, characteristicLine(E, tubeTarget, To, Mo_K, Mo_Ka2_f, Mo_Ka2, theta));
  				DoubleArrayUtil.add(phi, characteristicLine(E, tubeTarget, To, Mo_K, Mo_Kb1_f, Mo_Kb1, theta));
  				DoubleArrayUtil.add(phi, characteristicLine(E, tubeTarget, To, Mo_K, Mo_Kb2_f, Mo_Kb2, theta));
  			}
  		} else {
  			throw new Exception("Illegal target material in characteristic()");
  		}
  		return phi;
  	}
  
  	/** 
  	 * Return the vector spectrum phi that contains only the specified 
  	 * characteristic line with the appropriate split between the
  	 * nearest two bins.
  	 * 
  	 * @param E the energies
  	 * @param tubeTarget the material
  	 * @param To Peak kilo volage
  	 * @param Ek K-edge
  	 * @param char_line_f the flurescent yield
  	 * @param char_line_E the energy of the line
  	 * @param theta the angle
  	 * @return the array which contains the characteristic line.
  	 * @throws Exception may happen
  	 */
  	public static double [] characteristicLine(double [] E, String tubeTarget, double To, double Ek, double char_line_f, double char_line_E, double theta) throws Exception{
  		double [] phi = new double [E.length];
  		double Emin = E[0];
  		double dE = E[1] - E[0];
  		double fbin = ((char_line_E - Emin) / dE);
  		int il = (int) Math.floor(fbin);
  		int ih = il +1;
  		double n = N(tubeTarget, To, Ek, char_line_f, char_line_E, theta);
  		if (il < phi.length) {
  			phi[il] = (ih - fbin) * n;
  		}
  		if (ih < phi.length) {
  			phi[ih] = (fbin - il) * n;
  		}
  		return phi;
  	}
  
  	/**
  	 * Ref 1 Eq 22
  	 * Returns number of characteristic x rays per electron from target.
  	 * The coefficients Ak and nk are empirical values to fit the spectra
  	 * exposures to experimental values (see Tucker et al).
  	 * Micheal Moreau showed that there is an analytic solution to the integral
  	 * in this equation.  His solution is used
  	 * @param tubeTarget 'W' or 'Mo'
  	 * @param To
  	 * @param Ek K absorption energy (keV)
  	 * @param fi fractional emission of char line
  	 * @param Ei energy of char line
  	 * @param theta
  	 * @return number of x-rays
  	 * @throws Exception
  	 */
  	public static double N(String tubeTarget, double To, double Ek, double fi, double Ei, double theta) throws Exception{
  		/**
  		 * Ak		empir photons/electron
  		 * nk		empir exponent governing intensity
  		 * R         depth (cm) at which average kinetic energy of electrons 
  		 *           equals Ek, from T-W relation
  		 *
  		 * Evaluate integral over distance into target.  Use Gaussian
  		 * integration.  The constants Ak, nk, uf and R are first determined
  		 * appropriate for the target used
  		 */
  		double N = 1;
  		double Ak, nk, R;
  		double [] uf;
  		if (tubeTarget.equals("W")){
  			Ak = Ak_W;
  			nk = nk_W;
  			uf = DoubleArrayUtil.divide(DoubleArrayUtil.multiply(mu(new double [] {Ei}, tubeTarget), W_mdensity), Math.sin(theta));
  			R = (Math.pow(To,2) - Math.pow(Ek,2)) / (W_mdensity * c(To));
  		} else if (tubeTarget.equals("Mo")){
  			Ak = Ak_Mo;
  			nk = nk_Mo;
  			uf = DoubleArrayUtil.divide(DoubleArrayUtil.multiply(mu(new double [] {Ei}, tubeTarget), Mo_mdensity), Math.sin(theta));
  			R = (Math.pow(To,2) - Math.pow(Ek,2)) / (Mo_mdensity * c(To));
  		} else {
  			throw new Exception("Undefined target material in N(n)");
  		}
  		/* Evaluate the integral in Eq (22) using Micheal Moreau's solution.*/
  		double a = -uf[0];
  		double eaR = Math.exp(a*R);
  		double R2 = R*R;
  		double a2 = a*a;
  		double integral = 1.5/(R*a) * (eaR - 1. - 1./(R2)*(eaR*(R2 - 2.*R/a + 2./a2) - 2./a2));
  		N = Ak*Math.pow((To/Ek-1.0),nk)*fi*integral;	/*x rays/electro*/	
  		return N;
  	}
  
  	/**
  	 * Returns the Thomson-Whiddington constant using a rational approximation
  	 * (keV2 cm2 / g).  This routine fits the tabulated values in Birch and
  	 * Marshall's paper to within +/- 0.0025 x 1E6 and minimizes the rms error
  	 * subject to this restriction.  Maximum error less than +/- 0.4%
  	 * @param To incident electron energy (kVp)
  	 * @return  Thomson-Whiddington constant
  	 */
  	public static double c (double To){
  		double t0 = 0.20;
  		double t1 = 0.47;
  		double t2 = 0.21;
  		double t3 = 1.08;
  		double s0 = 1.00;
  		double s1 = -1.19;
  		double s2 = 3.00;
  		double x = 0.01 * To;
  		double c = ((t0 + x*(t1 + x*(t2 + x*t3))) / (s0 + x*(s1 + x*s2)) * 1.E06);
  		return c;
  	}
  
  	/**
  	 * Generates bremsstrahlung x-ray photon spectrum (photons/electron/bin)
  	 * using energies in vector E
  	 * @param E the energy
  	 * @param target the target 
  	 * @param To incident electron energy (kVp)
  	 * @param theta the angle
  	 * @return bremsstrahlung
  	 * @throws Exception 
  	 */
  	public static double [] bremsstrahlung(double [] E, String target, double To, double theta) throws Exception{
  		/**
  		 * These coefficients are used in the gaussian integration routines.
  		 * Use of 8 terms results in an integration accuracy of typically 0.01%.
  		 * The greater the number of terms, the longer is the computation time.
  		 */
  		NumberFormat nf = NumberFormat.getInstance();
  		nf.setMaximumFractionDigits(10);
  		nf.setMaximumIntegerDigits(10);
  		double NGAUSS = 4;
  		double [] Gx = {0.183434642495650, 0.525532409916329, 0.796666477413627, 0.960289856497536};
  		double [] Gw = {0.362683783378362, 0.313706645877887, 0.222381034453374, 0.101228536290376};
  		double dE = E[1] - E[0];
  		double f = 0; /* in cm2/g */
  		if (target.equals("W")){
  			f = dE*1000*ALPHA*Math.pow(RE,2)*Av * (TF_W*(Math.pow(W_Z,2)/W_A) + TF_Re*(Math.pow(Re_Z,2)/Re_A));	
  		} else if (target.equals("Mo")){
  			f = dE*1000*ALPHA*Math.pow(RE,2)*Av * Math.pow(Mo_Z,2)/Mo_A;
  		} else {
  			throw new Exception("Undefined target material in bremsstrahlung()");
  		}
  		/**
  		 * Evaluate integral over electron energies T from E to To.  Use
  		 * Gaussian integration
  		 */
  		double [] integral = new double [E.length];
  		double [] Tlolim = DoubleArrayUtil.min(Arrays.copyOf(E, E.length),To);
  		double Tuplim = To;
  		for (int j = 0; j < NGAUSS; j++) {
  			double [] T = multiply(add(multiply(Arrays.copyOf(Tlolim, Tlolim.length),(1. - Gx[j])), Tuplim *(1. + Gx[j])), 0.5);
  			DoubleArrayUtil.add(integral, DoubleArrayUtil.multiply(Integ(E,target,T,To,theta) , Gw[j]));
  			//DoubleArrayUtil.print("Integ", Integ(E,target,T,To,theta), nf);
  			T = multiply(add(multiply(Arrays.copyOf(Tlolim, Tlolim.length),(1. + Gx[j])), Tuplim *(1. - Gx[j])),0.5);
  			DoubleArrayUtil.add(integral, DoubleArrayUtil.multiply(Integ(E,target,T,To,theta), Gw[j]));
  			
  		}
  		/* g/(keV cm2) */
  		DoubleArrayUtil.multiply(integral, multiply(add(multiply(Arrays.copyOf(Tlolim,Tlolim.length),-1),Tuplim),0.5));  
  		double [] phi = DoubleArrayUtil.multiply(integral, f);
  		DoubleArrayUtil.divide(phi, E);
  		return phi;
  	}
  
  	public static double [] Integ(double [] Ei, String target, double [] T, double To, double theta) throws Exception {
  		double [] B = B(Ei,target,T,To);
  		NumberFormat nf = NumberFormat.getInstance();
  		nf.setMaximumFractionDigits(10);
  		nf.setMaximumIntegerDigits(10);
  		
  		DoubleArrayUtil.multiply(B, divide(add(Arrays.copyOf(T, T.length),MoC2),T));
  		
  		DoubleArrayUtil.multiply(B, F(Ei,target,T,To,theta));
  		//print("B", F(Ei,target,T,To,theta), nf);
  		DoubleArrayUtil.divide(B, MSP(target,T));
  		return B;
  	}
  
  	/**
  	 * 
  	 * Ref 1 Eq 18
  	 * Proportional to x-ray photons/electron (unitless).
  	 * The coefficients A0 and A1 are empirical values to fit the spectra
  	 * exposures to experimental values (see Tucker et al).  The values
  	 * of B_W and B_Mo are spectral shape parameters, and are not changed
  	 * from Tucker's empirically derived values
  	 * @param E
  	 * @param target
  	 * @param T
  	 * @param To
  	 * @return B
  	 * @throws Exception
  	 */
  	public static double [] B(double [] E, String target, double [] T, double To) throws Exception{
  		double [] B_W =  {-5.049, 10.847, -10.516, 3.842};
  		double [] B_Mo = {-4.238, 7.799, -6.739, 2.313};
  		double [] C = Arrays.copyOf(E, E.length);
  		DoubleArrayUtil.divide(C, T);
  		double A0, A1, B1, B2, B3, B4;
  		if (target.equals("W")){
  			A0 = A0_W;
  			A1 = A1_W;
  			B1 = B_W[0];
  			B2 = B_W[1];
  			B3 = B_W[2];
  			B4 = B_W[3];
  		} else if (target.equals("Mo")){
  			A0 = A0_Mo;
  			A1 = A1_Mo;
  			B1 = B_Mo[0];
  			B2 = B_Mo[1];
  			B3 = B_Mo[2];
  			B4 = B_Mo[3];
  		} else {
  			throw new Exception("Undefined target material in B()");
  		}
  		double factor = (A0 + A1*To);
  		double [] br = new double [C.length];
  		for (int i = 0; i < br.length; i++){
  			br[i] = factor *(1.0 + (B1*C[i]) + (B2*Math.pow(C[i],2)) + (B3*Math.pow(C[i],3)) + (B4*Math.pow(C[i],4)));
  		}
  		return br;
  	}
  
  	/**
  	 * Ref 1 Eq 16
  	 * Returns mass stopping power term for specified material and 
  	 * electron energy T
  	 * @param material
  	 * @param T
  	 * @return MSP
  	 * @throws Exception
  	 */
  	public static double [] MSP(String material, double [] T) throws Exception{
  		double Amsp_W = 2024.1;	    /* keV cm2 / g */
  		double Bmsp_W = 10361.0;	/* keV cm2 / g */
  		double Cmsp_W = 0.04695;	/* 1/keV */
  		double Amsp_Re = 2014.4;	/* keV cm2 / g */
  		double Bmsp_Re = 10276.0;	/* keV cm2 / g */
  		double Cmsp_Re = 0.04688;	/* 1/keV */
  		double Amsp_Mo = 2458.08;	/* keV cm2 / g */
  		double Bmsp_Mo = 14155.1;	/* keV cm2 / g */
  		double Cmsp_Mo = 0.04983;	/* 1/ke */
  		double [] msp = null;
  		if (material.equals("W")){
  			double [] temp = Arrays.copyOf(T, T.length);
  			multiply(temp, -Cmsp_W);
  			exp(temp);
  			multiply(temp, Bmsp_W);
  			multiply(add(temp, Amsp_W), TF_W);
  			msp = temp;
  			temp = Arrays.copyOf(T, T.length);
  			multiply(temp, -Cmsp_Re);
  			exp(temp);
  			multiply(temp, Bmsp_Re);
  			add(temp, Amsp_Re);
  			multiply(temp, TF_Re);
  			add(msp, temp);
  		} else if (material.equals("Mo")){
  			msp = new double[T.length];
  			for (int i = 0; i < T.length; i++) {
  				msp[i] = Amsp_Mo + Bmsp_Mo * Math.exp(-T[i] * Cmsp_Mo);
  			}
  		} else {
  			throw new Exception ("Undefined target material in MSP()");
  		}
  		return msp;
  	}
  
  	/**
  	 * Ref 1 Eq 11
  	 * Returns fraction of x-ray photons of energy E escaping target material.
  	 * Note: mu() returns mass attenuation data, so we can drop the density 
  	 * term in Tucker's paper
  	 * @param E
  	 * @param target
  	 * @param T
  	 * @param To
  	 * @param theta
  	 * @return F
  	 * @throws Exception
  	 */
  	public static double [] F(double [] E, String target, double [] T, double To, double theta) throws Exception{
  		double [] u;
  		if (target.equals("W")){
  			u = DoubleArrayUtil.multiply(mu(E, "W"), TF_W);
  			DoubleArrayUtil.add(u, DoubleArrayUtil.multiply(mu(E,"Re"),TF_Re));
  		} else if (target.equals("Mo")){
  			u = mu(E, "Mo");
  		} else {
  			throw new Exception ("Undefined target material in F()");
  		}
  		double factor [] = divide(add( multiply(pow(Arrays.copyOf(T, T.length),2), -1),Math.pow(To,2)), (c(To)*Math.sin(theta)));
  		NumberFormat nf = NumberFormat.getInstance();
  		nf.setMaximumFractionDigits(10);
  		nf.setMaximumIntegerDigits(10);
  		//print("factor",mu(E, "W"), nf);
  		DoubleArrayUtil.multiply(u, -1);
  		DoubleArrayUtil.multiply(u, factor);
  		DoubleArrayUtil.exp(u);
  		return u;
  	}
  	
  	/**
  	 * Returns the number of quanta per R of exposure at the energies
  	 * specified by vector E (keV).
  	 * @param airAbsorption
  	 * @return number of quanta
  	 * @throws Exception
  	 */
  	public static double xrqpRe(double airAbsorption) throws Exception{
  		
  		double E_CHG = 1.6021892E-19;   // charge on an electron (C)
  		double C_kg_R = 2.580E-4;       // Coul/kg/R for air
  		double W_air = 33.97;           // work function for air (eV)
  
  		// Calculate quanta/R values		
  		double f = W_air * C_kg_R * 10.0 / (E_CHG * 1E9);
  		
  		return f / airAbsorption;		
  	}
  
  	public static void main(String [] argv){	
  		double [] energy = {10, 20, 30, 40, 50, 60, 70, 80};
  		try {
  			double [] F = generateXRaySpectrum(energy, 120, "W", 1);
  			DoubleArrayUtil.print("E", energy);
  			NumberFormat nf = NumberFormat.getInstance();
  			nf.setMaximumFractionDigits(10);
  			nf.setMaximumIntegerDigits(10);
  			DoubleArrayUtil.print("F", F, nf);
  						
  			/**
  			 * Half value layer (HVL)
  			 * 1. Generate x-ray spectrum
  			 * 2. Optimize x thickness for HVL using optimizer
  			 */
  			double min = 10.0; // Minimum [keV]
  			double max = 150.0; // Maximum [keV]
  			double delta = 0.5; // Delta [keV]
  			double peakVoltage = 125.0; // Peak Voltage [kVp]			
  			double timeCurrentProduct = 1.0; // Time Current Product [mAs]
  			
  			int steps = (int) ((max - min) / delta);
  			double [] energies = new double [steps];
  			for (int i = 0; i < steps; i ++){
  				energies[i] = min + (i*delta);
  			}
  			double [] spectrum = XRaySpectrum.generateXRaySpectrum(energies, peakVoltage, "W", timeCurrentProduct, 1, 12, 2.38, 3.06, 2.66, 1.2); // C-arm: 1.2mm, CT: 10.5mm			
  			//* PCXMC 1.5, Al 1.2mm, 125KVp 
  			//double [] spectrum = {2.339734153, 4.041388504, 27.38163849, 15.39283896, 37.18279457, 73.56337068, 126.0434512, 194.3001899, 276.7188396, 370.9872404, 474.525746, 584.7174788, 698.9984541, 814.8879307, 930.0170453, 1042.178469, 1149.392328, 1249.971553, 1342.570514, 1426.207848, 1500.26219, 1564.445105, 1618.758252, 1663.442134, 1698.922723, 1725.760589, 1744.605441, 1756.157569, 1761.136588, 1760.257185, 1754.21114, 1743.654715, 1729.20045, 1711.412455, 1690.804388, 1667.839433, 1642.931697, 1616.448585, 1588.713784, 1560.0106, 1530.585431, 1500.651249, 1470.390973, 1439.960685, 1409.492634, 1379.098007, 1348.869472, 1318.883473, 4124.783498, 6130.065636, 1230.943841, 1202.436086, 1174.375023, 1146.776223, 1119.649532, 1092.999939, 1066.828356, 2502.684011, 1015.906379, 1366.607548, 747.2849712, 683.2666223, 675.4686445, 667.3257692, 658.8505088, 650.0551081, 640.9515346, 631.5514713, 621.8663123, 611.907161, 601.6848301, 591.2098444, 580.4924441, 569.542591, 558.3699745, 546.9840201, 535.3938979, 523.6085318, 511.6366095, 499.4865931, 487.1667296, 474.6850614, 462.0494375, 449.2675242, 436.3468158, 423.2946452, 410.118194, 396.824503, 383.4204816, 369.9129176, 356.308486, 342.6137582, 328.8352099, 314.9792297, 301.0521261, 287.0601356, 273.0094287, 258.9061174, 244.7562604, 230.5658698, 216.340916, 202.0873332, 187.8110241, 173.5178642, 159.2137066, 144.9043855, 130.59572, 116.2935175, 102.0035769, 87.73169164, 73.48365225, 59.26524903, 45.08227425, 30.94052435, 16.84580188, 3.363752905, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
  			//* PCXMC 1.5, Al 1.2mm, 109KVp
  			//double [] spectrum = {2.955266803, 5.534583134, 35.53063201, 20.78504989, 50.10945855, 98.9266474, 169.1105885, 260.0284174, 369.2513097, 493.3143776, 628.2639988, 769.9886413, 914.4277354, 1057.738864, 1196.452481, 1327.605535, 1448.83427, 1558.413734, 1655.244323, 1738.795765, 1809.023517, 1866.272443, 1911.179806, 1944.585811, 1967.456388, 1980.820069, 1985.718866, 1983.171895, 1974.149924, 1959.558842, 1940.230127, 1916.916613, 1890.292118, 1860.95375, 1829.426, 1796.165903, 1761.568761, 1725.974056, 1689.671287, 1652.905576, 1615.882911, 1578.774999, 1541.723672, 1504.844866, 1468.232175, 1431.960007, 1396.086368, 1360.655307, 3472.502797, 4966.973123, 1257.291784, 1223.861716, 1190.950986, 1158.5562, 1126.670174, 1095.282699, 1064.381199, 2122.305934, 1003.977217, 1253.30148, 786.1018455, 731.0427146, 716.8576874, 702.2039203, 687.1004379, 671.5660379, 655.6192709, 639.2784244, 622.5615106, 605.4862593, 588.070113, 570.3302259, 552.283465, 533.9464137, 515.3353769, 496.4663881, 477.3552179, 458.0173827, 438.4681553, 418.7225753, 398.7954603, 378.7014171, 358.4548537, 338.0699904, 317.5608713, 296.9413762, 276.2252309, 255.4260187, 234.5571907, 213.6320761, 192.6638916, 171.6657515, 150.6506761, 129.6316007, 108.6213836, 87.63281374, 66.67861822, 45.77146919, 24.92399031, 4.977138706, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
  			//* PCXMC 1.5, Al 1.2mm, 81KVp
  			//double [] spectrum = {5.378719397, 11.49101655, 67.12856602, 42.13840172, 100.9364506, 197.8129205, 335.2339375, 510.0244283, 714.7811096, 939.6352271, 1173.886102, 1407.319085, 1631.128801, 1838.432695, 2024.422216, 2186.243436, 2322.712599, 2433.959903, 2521.069511, 2585.757072, 2630.104018, 2656.352883, 2666.75903, 2663.490021, 2648.562707, 2623.808682, 2590.860072, 2551.149209, 2505.917269, 2456.228218, 2402.9855, 2346.949702, 2288.756026, 2228.930879, 2167.907149, 2106.038006, 2043.609155, 1980.849596, 1917.940972, 1855.025634, 1792.213549, 1729.588181, 1667.211485, 1605.128117, 1543.36897, 1481.954142, 1420.895412, 1360.19829, 1843.377528, 2165.763967, 1180.271382, 1121.004084, 1062.08199, 1003.49979, 945.2529603, 887.338147, 829.7534654, 1040.71992, 715.5755988, 727.4210604, 573.9540227, 517.2982556, 468.0250995, 418.2334156, 368.0000507, 317.4008233, 266.5105365, 215.4029968, 164.1510415, 112.8265696, 61.50057796, 12.2895214, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
  			//* PCXMC 1.5, Al 1.2mm, 70KVp
  			//double [] spectrum = {7.481922341, 16.50277126, 92.86281442, 60.83853289, 144.9671387, 282.2596166, 474.4437357, 714.5596079, 989.5192823, 1283.358776, 1580.139999, 1865.998848, 2130.235548, 2365.596519, 2568.006225, 2736.006804, 2870.104526, 2972.149214, 3044.810436, 3091.171012, 3114.433025, 3117.719489, 3103.951195, 3075.779185, 3035.556262, 2985.334486, 2926.878875, 2861.690343, 2791.033077, 2715.963189, 2637.35668, 2555.935555, 2472.291509, 2386.906945, 2300.173352, 2212.407152, 2123.863256, 2034.746552, 1945.221583, 1855.420659, 1765.450604, 1675.398367, 1585.335643, 1495.322682, 1405.411399, 1315.647914, 1226.074601, 1136.731752, 1052.909937, 967.8248154, 870.4834816, 782.4645221, 694.8838466, 607.7889466, 521.2301659, 435.2608781, 349.9375973, 267.8868369, 181.4711198, 99.11088352, 19.59335481, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
  			//* PCXMC 1.5, Al 0.0mm, 70KVp, Anode 6 deg
  			//double [] spectrum = {7.481922341, 16.50277126, 92.86281442, 60.83853289, 144.9671387, 282.2596166, 474.4437357, 714.5596079, 989.5192823, 1283.358776, 1580.139999, 1865.998848, 2130.235548, 2365.596519, 2568.006225, 2736.006804, 2870.104526, 2972.149214, 3044.810436, 3091.171012, 3114.433025, 3117.719489, 3103.951195, 3075.779185, 3035.556262, 2985.334486, 2926.878875, 2861.690343, 2791.033077, 2715.963189, 2637.35668, 2555.935555, 2472.291509, 2386.906945, 2300.173352, 2212.407152, 2123.863256, 2034.746552, 1945.221583, 1855.420659, 1765.450604, 1675.398367, 1585.335643, 1495.322682, 1405.411399, 1315.647914, 1226.074601, 1136.731752, 1052.909937, 967.8248154, 870.4834816, 782.4645221, 694.8838466, 607.7889466, 521.2301659, 435.2608781, 349.9375973, 267.8868369, 181.4711198, 99.11088352, 19.59335481, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
  			
  			
  			//System.out.println(spectrum.length);
  			
  			HalfValueLayerFunction hvlf = new HalfValueLayerFunction(energies, spectrum, 
  					EnergyDependentCoefficients.getPhotonMassAttenuationLUT(EnergyDependentCoefficients.Material.Aluminum), 
  					EnergyDependentCoefficients.getMassEnergyAbsorptionLUT(EnergyDependentCoefficients.Material.CsI),
  					new LinearInterpolatingDoubleArray(EnergyDependentCoefficients.airEnergies, EnergyDependentCoefficients.airAbsoprtion));
  			hvlf.runOptimization();
  	
  			System.out.println(peakVoltage + "kVp\t" + hvlf.getOptimalX() * 10 + "[mm]\t" + hvlf.evaluate(hvlf.optimalX)+ "\t" + CONRAD.DOUBLE_EPSILON);
  			
  		} catch (Exception e) {
  			e.printStackTrace();
  		}
  	}
  	
  }
  