content: /*
   * Copyright (C) 2010-2014 Rotimi X Ojo, Andreas Maier
   * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
   */
  package edu.stanford.rsl.conrad.physics;
  
  import edu.stanford.rsl.conrad.utils.CONRAD;
  
  /**
   * <p> This class creates a model of a polychromatic X Ray spectrum.</p>
   * We model a spectrum as a set of bins.<br>
   *
   * Note that we use the term photon flux here for a quantity that reports the number of photons per area.<br>
   *
   * <p>The default parameters are:<br/>
   * minimum energy = 10keV
   * maximum energy = 150kev
   * resolution delta= 0.5
   * peak voltage = 90kVp
   * time current product = 2.5;</p>
   * 
   * @author Rotimi X Ojo
   * @author Andreas Maier
   */
  public class PolychromaticXRaySpectrum {
  
  	//private NumberInterpolatingTreeMap map = new NumberInterpolatingTreeMap();
  	private double min = 10.0; // Minimum [keV] 
  	private double max = 150.0; // Maximum [keV]
  	private double delta = 0.5; // Delta [keV]
  	private double peakVoltage = 90.0; //125.0; // Peak Voltage [kVp]			
  	private double timeCurrentProduct = 2.5;// 1.0; // Time Current Product [mAs]
  	// (X-Ray Tube current 447 mA * Average Pulse Width 5.6 ms )/1000 gives mAs = 2.5
  	private double totalPhotonFlux;
  	private double[] energies;
  	private double[] photonFlux;
  	private double avgPhotonEnergy;
  
  	/**
  	 * Creates a new polychromatic X-Ray spectrum satisfying default parameters.
  	 */
  	public PolychromaticXRaySpectrum(){	
  		generateSpectrum(min, max, delta, peakVoltage, "W", timeCurrentProduct, 1, 10, 2.38, 3.06, 2.66, 10.5);
  	}
  
  	/**
  	 * Creates a new polychromatic X-Ray spectrum with successive energies having a difference of delta
  	 * @param delta is the difference between successive energies starting at min
  	 */
  	public PolychromaticXRaySpectrum(double delta){	
  		this.delta = delta;
  		generateSpectrum(min, max, delta, peakVoltage, "W", timeCurrentProduct, 1, 10, 2.38, 3.06, 2.66, 10.5);
  	}
  
  	/**
  	 * Creates a new polychromatic X-Ray spectrum satisfying the parameters below.
  	 * @param min is minimum energy in keV
  	 * @param max is maximum energy in keV
  	 * @param delta is resolution
  	 * @param peakVoltage is peak voltage 
  	 * @param timeCurrentProduct is time current product
  	 */
  	public PolychromaticXRaySpectrum(double min, double max, double delta, double peakVoltage, double timeCurrentProduct){
  		this.min = min;
  		this.max = max;
  		this.delta = delta;
  		this.peakVoltage = peakVoltage;
  		this.timeCurrentProduct = timeCurrentProduct;
  		generateSpectrum(min, max, delta, peakVoltage, "W", timeCurrentProduct, 1, 10, 2.38, 3.06, 2.66, 10.5);
  	}
  
  	/**
  	 * Constructor to create a fully configured spectrum
  	 * @param min is minimum energy in keV
  	 * @param max is maximum energy in keV
  	 * @param delta is resolution
  	 * @param peakVoltage is peak voltage 
  	 * @param anodeMaterial the target material ("W" or "Mo")
  	 * @param timeCurrentProduct the acceleration time times the current in [mAs]
  	 * @param mdis amount of air in [m] (i.e. the source detector distance)
  	 * @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]
  	 */
  	public PolychromaticXRaySpectrum(double min, double max, double delta, double peakVoltage, String anodeMaterial, double timeCurrentProduct, double mdis, double degrees, double mmpyrex, double mmoil, double mmlexan, double mmAl){
  		generateSpectrum(min, max, delta, peakVoltage, anodeMaterial, timeCurrentProduct, mdis, degrees, mmpyrex, mmoil, mmlexan, mmAl);
  		this.min = min;
  		this.max = max;
  		this.delta = delta;
  		this.peakVoltage = peakVoltage;
  		this.timeCurrentProduct = timeCurrentProduct;
  	}
  
  	/**
  	 * Constructor to create an object of an already known spectrum
  	 * @param min minimum energy in keV
  	 * @param max maximum energy in keV
  	 * @param delta resolution
  	 * @param values known values for each bin
  	 * @param timeCurrentProduct the acceleration time times the current in [mAs]
  	 */
  	public PolychromaticXRaySpectrum(double min, double max, double delta,
  			double[] values, double timeCurrentProduct) {
  		this.min = min;
  		this.max = max;
  		this.delta = delta;
  		this.timeCurrentProduct = timeCurrentProduct;
  		int numberOfBins = (int) ((max - min) / delta);
  		energies = new double[numberOfBins];
  		photonFlux = new double[numberOfBins];
  		double weightedIntegral = 0.0;
  		double peakVoltage = 0.0;
  		for (int i = 0; i < numberOfBins; i++) {
  			energies[i] = min + delta * i;
  			photonFlux[i] = values[i];
  			totalPhotonFlux += photonFlux[i];
  			weightedIntegral += energies[i] * photonFlux[i];
  			if (values[i] > CONRAD.DOUBLE_EPSILON) {
  				peakVoltage = energies[i];
  			}
  		}
  		this.peakVoltage = peakVoltage;
  		avgPhotonEnergy = weightedIntegral / totalPhotonFlux;
  	}
  
  	/**
  	 * Generate a X-Ray spectrum satisfying the parameters below.
  	 * @param min is minimum energy in keV
  	 * @param max is maximum energy in keV
  	 * @param delta is resolution
  	 * @param peakVoltage is peak voltage 
  	 * @param timeCurrentProduct is time current product
  	 */
  	private void generateSpectrum(double min, double max, double delta,
  			double peakVoltage, String anodeMaterial, double timeCurrentProduct, double mdis, double degrees, double mmpyrex, double mmoil, double mmlexan, double mmAl) {		
  		energies = generateEnergies(min, max, delta);
  		try {
  			photonFlux = XRaySpectrum.generateXRaySpectrum(energies, peakVoltage, anodeMaterial, timeCurrentProduct, mdis, degrees, mmpyrex, mmoil, mmlexan, mmAl);
  			totalPhotonFlux = 0;
  			double weightedIntegral = 0;
  			for (int i = 0; i < energies.length; i++) {
  				totalPhotonFlux += photonFlux[i];
  				weightedIntegral += energies[i] * photonFlux[i];
  			}
  			avgPhotonEnergy = weightedIntegral / totalPhotonFlux;
  		} catch (Exception e) {
  			e.printStackTrace();
  		}
  	}
  
  	/**
  	 * Photon intensity of spectral bin given at a certain energy (kev) In this
  	 * implementation, we round the keV to the next bin of the spectrum.
  	 * 
  	 * @param energy is energy in [keV]
  	 * @return photon intensity in [keV/mm2] 
  	 */
  	public double getIntensity(double energy) {
  		int bin = (int) Math.round((energy - min) / delta);
  		if (bin < energies.length)
  			return photonFlux[bin] * energies[bin];
  		else
  			return 0;
  	}
  
  	/**
  	 * Photon flux of XRay at given energy (kev) In this implementation, we
  	 * round the keV to the next bin of the spectrum.
  	 * 
  	 * @param energy is energy in [keV]
  	 * @return x-ray photon flux in [photons/mm2] 
  	 */
  	public double getPhotonFlux(double energy) {
  		int bin = (int) Math.round((energy - min) / delta);
  		if (bin < energies.length)
  			return photonFlux[bin];
  		else
  			return 0;
  	}
  
  	/**
  	 * Returns the total energy of the spectrum per square millimeter, i.e. the integral over all energy bins.
  	 * 
  	 * @return the total energy in [keV/mm^2]
  	 */
  	public double getTotalIntensity() {
  		return totalPhotonFlux * avgPhotonEnergy;
  	}
  
  
  	/**
  	 * Returns the total photon flux of the spectrum, i.e. the integral over all energy bins.
  	 * 
  	 * @return the total photon flux in [photons/mm^2]
  	 */
  	public double getTotalPhotonFlux() {
  		return totalPhotonFlux;
  	}
  
  	/**
  	 * Computes the average phonton energy of this spectrum given as the weighted sum of energies multiplied with photon flux normalized by the photon flux.
  	 * 
  	 * @return average photon enegy in [keV]
  	 */
  	public double getAveragePhotonEnergy() {
  		return avgPhotonEnergy;
  	}
  
  
  	/**
  	 * @return photon energies in kev
  	 */
  	public double[] getPhotonEnergies() {
  		return energies;
  	}
  	
  	/**
  	 * set photon energies in kev
  	 */
  	public void setPhotonEnergies(double[] e) {
  		energies=e;
  	}
  
  	/**
  	 * Determine the number of discrete energies used to describe the spectrum
  	 * @return 0 if there are energies in the spectrum.
  	 */
  	public int size() {
  		return energies.length;
  	}
  
  	/**
  	 * returns the photon flux array containing for a spectral bin.
  	 * @return the photon flux the array of fluxes in [photons/mm^2]
  	 */
  	public double[] getPhotonFlux() {
  		return photonFlux;
  	}
  	
  	/**
  	 *
  	 * @return the timeCurrentProduct of the spectrum
  	 */
  	public double getTimeCurrentProduct(){
  		return timeCurrentProduct;
  	}
  
  	/**
  	 * @param photonFlux the photon flux to set
  	 */
  	public void setPhotonFlux(double[] photonFlux) {
  		this.photonFlux = photonFlux;
  	}
  
  	private double[] generateEnergies(double min, double max, double delta) {
  		int steps = (int) ((max - min) / delta);
  		double[] pen = new double[steps];
  		for (int i = 0; i < steps; i++) {
  			pen[i] = min + (i * delta);
  		}
  		return pen;
  	}
  
  	/**
  	 * @return the min
  	 */
  	public double getMin() {
  		return min;
  	}
  
  	/**
  	 * @param min the min to set
  	 */
  	public void setMin(double min) {
  		this.min = min;
  		generateSpectrum(min, max, delta, peakVoltage, "W", timeCurrentProduct, 1, 10, 2.38, 3.06, 2.66, 10.5);
  	}
  
  	/**
  	 * @return the max
  	 */
  	public double getMax() {
  		return max;
  	}
  
  	/**
  	 * @param max the max to set
  	 */
  	public void setMax(double max) {
  		this.max = max;
  		generateSpectrum(min, max, delta, peakVoltage, "W", timeCurrentProduct, 1, 10, 2.38, 3.06, 2.66, 10.5);
  	}
  
  	/**
  	 * @return the delta
  	 */
  	public double getDelta() {
  		return delta;
  	}
  
  	/**
  	 * @param delta the delta to set
  	 */
  	public void setDelta(double delta) {
  		this.delta = delta;
  		generateSpectrum(min, max, delta, peakVoltage, "W", timeCurrentProduct, 1, 10, 2.38, 3.06, 2.66, 10.5);
  	}
  
  	/**
  	 * @return the peakVoltage
  	 */
  	public double getPeakVoltage() {
  		return peakVoltage;
  	}
  
  	/**
  	 * @param peakVoltage the peakVoltage to set
  	 */
  	public void setPeakVoltage(double peakVoltage) {
  		this.peakVoltage = peakVoltage;
  		generateSpectrum(min, max, delta, peakVoltage, "W", timeCurrentProduct, 1, 10, 2.38, 3.06, 2.66, 10.5);
  	}
  
  	/**
  	 * @param timeCurrentProduct the timeCurrentProduct to set
  	 */
  	public void setTimeCurrentProduct(double timeCurrentProduct) {
  		this.timeCurrentProduct = timeCurrentProduct;
  		generateSpectrum(min, max, delta, peakVoltage, "W", timeCurrentProduct, 1, 10, 2.38, 3.06, 2.66, 10.5);
  	}
  
  }
  