/**
 * @file nubeam.h 
 * @brief Nubeam C++ class
 * @version $Id: nubeam.h 863 2011-02-24 19:46:44Z pletzer $
 */

// std includes
#include <string>
#include <vector>
#include <map>
#include <cstdlib>

// nubeam includes
#include "cnubeam_comp.h"

#ifndef NUBEAM_H
#define NUBEAM_H

class NubeamError {
 public:
  NubeamError(const std::string& fN, int lN, int ier) :
    fileName(fN), lineNumber(lN), errorCode(ier) {};
  ~NubeamError() {};
 private:
  std::string fileName;
  int lineNumber;
  int errorCode;
};

class Nubeam {

 public:
  /** 
   * Constructor, no argument
   */
  Nubeam();

  /** 
   * Destructor, virtual
   */
  virtual ~Nubeam();

  /**
   * Set Fortran MPI communicator 
   * @param comm communicator
   */
  void setFortranComm(int commF);

  /**
   * Set PREACTDIR environment variable to where atomic reaction tables reside
   * @param path
   */
  void setPreactdir(const std::string& path);

  /**
   * Set ADASDIR environment variable to where atomic reaction tables reside
   * @param path
   */
  void setAdasdir(const std::string& path);

  /**
   * Set plasma state file
   * @param file file name
   */
  void setPlasmaStateFile(const std::string& file);

  /**
   * Set run ID
   * @param runid string
   */
  void setRunid(const std::string& runid);

  /**
   * Set time interval
   * @param t0 start time
   * @param t1 end time
   */  
  void setTimeInterval(const double t0, const double t1);

  /**
   * Set random number generator seed
   * @param n seed
   */
  void setRandomNumberGeneratorSeed(const int n);

  /**
   * Set g-eqdsk file
   * @param geqdsk file name
   * @param minRadius minimal radius at X-point
   */
  void setGeqdskFile(const std::string& geqdsk, const double minRadius=0.75);

  /**
   * Load machine description from file
   * @param file machine description file 
   * @param geqdsk g-eqdsk file to extract limiter geometry
   */
  void setMachineDescriptionFile(const std::string& file);

  /**
   * Set grid by name
   * @param name grid name (e.g "rho", "rho_rw", ...)
   * @param values nodal values within [0,1]
   */
  void setGrid(const std::string& name, const std::vector<double>& values);

  /**
   * Set thermal species
   * @param z number of protons
   * @param c ionization (<=z)
   * @param a atomic mass 
   */
  void setThermalSpecies(const int z, const int c, const int a);

  /** 
   * Set neutral beam species
   * @param z number of protons
   * @param c ionization (<=z)
   * @param a atomic mass 
   */
  void setNeutralBeamSpecies(const int z, const int c, const int a);

  /** 
   * Set impurity species
   * @param z number of protons
   * @param c ionization (<=z)
   * @param a atomic mass 
   */
  void setImpuritySpecies(const int z, const int c, const int a);

  /**
   * Set the file name for initializing the nubeam control struct
   * @param file file name
   */
  void setNbCtrlInitFile(const std::string& file);

  /**
   * Initialize, allocate data, etc
   */
  void initialize();

  /**
   * Set number of (normal) beam particles
   * @param number of particles
   */
  void setNumberOfBeamParticles(const int n);

  /**
   * Set number of fusion particles
   * @param number of particles
   */
  void setNumberOfFusionParticles(const int n);

  /**
   * Set fractions of beam currents 
   * @param atHalfVoltage at half voltage
   * @param atFullVoltage at full voltage
   */
  void setFractionsOfBeamCurrents(const std::vector<double>& atHalfVoltage, 
				  const std::vector<double>& atFullVoltage);

  /**
   * Set Z effective (=sum[nj*Zj**2]/ne)
   * @param zeff values
   */
  void setZEffective(const std::vector<double>& zeff);  

  /**
   * Set cell centered temperatures
   * @param te electron temperature in eV
   * @param ti ion temperature in eV
   */
  void setTemperatures(const std::vector<double>& te,
		       const std::vector<double>& ti);

  /**
   * Set cell centered densities
   * @param ne electron density in m^-3
   * @param ni ion density in m^-3
   */
  void setDensities(const std::vector<double>& ne,
		       const std::vector<double>& ni);

  /**
   * Get the number of beams
   * @return number
   */
  int getNumBeams();

  /**
   * Set beam energies
   * @param evolt_nbi energies in eV, one per beam 
   */
  void setBeamEnergies(const std::vector<double>& evolt_nbi);

  /**
   * Set beam powers
   * @param power_nbi powers in W, one per beam 
   */
  void setBeamPowers(const std::vector<double>& power_nbi);

  /**
   * Turn on the halo thermal neutral recapture using built-in FRANTIC
   */
  void haloThermalNeutralRecaptureOn();

  /**
   * Turn off the halo thermal neutral recapture
   */
  void haloThermalNeutralRecaptureOff();

  /** 
   * Turn debug on, will generate input/output files for running in standalone
   */
  void debugOn();

  /** 
   * Turn debug off
   */
  void debugOff();

  /**
   * Update nubeam control struct from file
   * @param file file name
   */
  void updateNbCtrl(const std::string& file);

  /**
   * Check and fix plasma state
   */
  void checkAndFix(); 

  /**
   * Save plasma state into file
   * @param file Netcdf file name
   */
  void savePlasmaState(const std::string& file);

  /**
   * Enforce quasi-neutrality on electrons
   */
  void enforceQuasiNeutrality();

  /** 
   * Compute sources
   */
  void run();

  /**
   * Load energy sources
   * @param rho target grid 
   * @param pe electron energy density source in W/m^3 (out)
   * @param pi ion energy density source in W/m^3 (out)
   */
  void loadEnergySources(const std::vector<double>& rho,
		  std::vector<double>& pe,
		  std::vector<double>& pi) const;

  /**
   * Load particle sources
   * @param rho target grid 
   * @param pe electron energy density source in W/m^3 (out)
   * @param pi ion energy density source in W/m^3 (out)
   */
  void loadParticleSources(const std::vector<double>& rho,
		  std::vector<double>& se,
		  std::vector<double>& si) const;

  /** 
   * Get density profile from plasma state (m^-3)
   * @param specieIndex 0=electron, >0 for ions
   * @return profile
   */
  std::vector<double> getDensityFromPlasmaState(int specieIndex);

  /** 
   * Get temperature profile from plasma state (keV)
   * @param specieIndex 0=electron, >0 for ions
   * @return profile
   */
  std::vector<double> getTemperatureFromPlasmaState(int specieIndex);

 private:

  /** Forbid copy-constructor */
  Nubeam(const Nubeam&);

  /** Forbid assignment */
  Nubeam& operator=(const Nubeam&);

  /**
   * Interpolation of 1d profile
   * @param rhoOri original grid
   * @param datOri original data
   * @param rho target grid
   * @return interpolated data
   */
  std::vector<double> interp1d(const std::vector<double>& rhoOri, 
			       const std::vector<double>& datOri,
			       const std::vector<double>& rho) const;
  /** 
   * Sets perBeamSpecies according to number of protons 
   * and atomic mass 
   * @param zi number of protons
   * @param ai atomic mass 
   */
  void setSpeciesNameFromZA(const int zi, const int ai);

  /** Catch error 
   * @param ier error code (0=OK)
   * @param lineNumber line number where error was detected
   * @param fileName file name where error occurred
   */
  void catchError(const int ier, const int lineNumber, 
		  const char* fileName) const;

  /** Grid name to grid values */
  std::map<std::string, std::vector<double> > gridMap;

  /** geqdsk file name */
  std::string geqdsk;
  
  /** minimum radius near X-point */
  double minRadius;

  /** Plasma state objects */
  PlasmaState* ps;

  /** contains initial settings for nubeam */
  struct nbiCtrlInit nbCtrlInit;

  /** contains nubeam settings that are likely updated stepwise */
  struct nbiCtrlUpdate nbCtrlUpdt;

  /** nb ctrl init file name */
  std::string fileNbCtrlInit;

  /** the species name for neutral beam channels (assumed all to tbe same) */
  std::string perBeamSpecies;

  /** unique run ID string */
  std::string runid;

  /** plasma state file, if non-empty will load from this file */
  std::string plasmaStateFile;

  /** machine description file */
  std::string machDescr;

  /** list of thermal species {ionization (z), charge, atomic number (A)} */
  std::vector<std::vector<int> > thermalSpecies;

  /** list of impurities species {ionization (z), charge, atomic number (A)} */
  std::vector<std::vector<int> > impuritySpecies;

  /** list of neutral beam species {ionization (z), charge, atomic number (A)} */
  std::vector<std::vector<int> > neutralBeamSpecies;

  /** random number generator seed */
  int rngSeed;
  
  /** number of beam particles */
  int numberOfBeamParticles;
  
  /** number of fusion particles */
  int numberOfFusionParticles;

  /** Fortran communicator */
  int commF;

  /** my process */
  int myid;

};
#endif // NUBEAM_H
