#ifndef __JPHYSICS_JPHYSICSSUPPORTKIT__
#define __JPHYSICS_JPHYSICSSUPPORTKIT__

#include <cmath>

#include "JTools/JFunction1D_t.hh"

#include "JPhysics/JConstants.hh"


/**
 * \author mdejong
 */

namespace JPHYSICS {}
namespace JPP { using namespace JPHYSICS; }

/**
 * Auxiliary methods for light properties of deep-sea water.
 */
namespace JPHYSICS {

  using JTOOLS::JPolint1Function1D_t;
  using JTOOLS::JGridSplineFunction1D_t;

  
  /**
   * Get multiple Coulomb scattering angle.
   * 
   * The formula is taken from reference:
   * Particle Data Book, formula 27.14.
   *
   * \param  E          Energy           [GeV]
   * \param  x          distance         [m]
   * \param  X0         radiation length [m]
   * \param  M          mass             [GeV]
   * \param  Q          charge           [unit]
   * \return            angle            [rad]
   */
  inline double getThetaMCS(const double E,
			    const double x,
			    const double X0,
			    const double M,
			    const double Q)
  {
    if (E > M) {

      const double p    = sqrt((E + M) * (E - M));
      const double beta = p / E;
      
      return THETA_MCS * Q * sqrt(x/X0) * (1.0 + 0.038*log(x/X0)) / (beta*p);
    }

    return 0.0;
  }

  
  /**
   * Get multiple Coulomb scattering angle for muon.
   *
   * \param  E          Energy [GeV]
   * \param  x          distance [m]
   * \return            angle [rad]
   */
  inline double getThetaMCS(const double E, const double x)
  {
    return getThetaMCS(E, x, X0_WATER_M, MASS_MUON, 1.0);
  }

  
  /**
   * Auxiliary method to describe light scattering in water (Henyey-Greenstein).
   *
   * \param  g      angular dependence parameter
   * \param  x      cosine scattering angle
   * \return        probability
   */
  inline double henyey_greenstein(const double g,
				  const double x)
    
  {
    const double a0 = (1.0 - g*g) / (4*PI);
    const double y  =  1.0 + g*g - 2.0*g*x;
    
    return a0 / (y*sqrt(y));
  }
  

  /**
   * Auxiliary method to describe light scattering in water (Heneyey-Greenstein).
   *
   * \param  x      cosine scattering angle
   * \return        probability
   */
  inline double henyey_greenstein(const double x)
  {
    static const double g = 0.924;
    
    return henyey_greenstein(g, x);
  }


  /**
   * Auxiliary method to describe light scattering in water (Rayleigh).
   *
   * \param  a      angular dependence parameter
   * \param  x      cosine scattering angle
   * \return        probability
   */
  inline double rayleigh(const double a,
			 const double x)

  {
    const double a0 = 1.0 / (1.0 + a/3.0) / (4*PI);

    return a0 * (1.0 + a*x*x); 
  }


  /**
   * Auxiliary method to describe light scattering in water (Rayleigh).
   *
   * \param  x      cosine scattering angle
   * \return        probability
   */
  inline double rayleigh(const double x)
  {
    return rayleigh(0.835, x);
  }


  /**
   * Model specific function to describe light scattering in water (f4).
   *
   * \param  x      cosine scattering angle
   * \return        probability
   */
  inline double f4(const double x)
  {
    static const double g1 = 0.77;
    static const double g2 = 0.00;
    static const double f  = 1.00;

    return f * henyey_greenstein(g1,x)  +  (1.0 - f) * henyey_greenstein(g2,x);
  }


  /**
   * Model specific function to describe light scattering in water (p00075).
   *
   * \param  x      cosine scattering angle
   * \return        probability
   */
  inline double p00075(const double x)
  {
    static const double g = 0.924;
    static const double f = 0.17;

    return f * rayleigh(x)  +  (1.0 - f) * henyey_greenstein(g,x);
  }


  /**
   * Measurement of light scattering in water.
   *
   * Values are taken from reference C.D. Mobley "Light and Water", ISBN 0-12-502750-8, pag. 111, Table 3.10 (right column).
   */
  class JPetzhold :
    public JPolint1Function1D_t
  {
  public:
    /**
     * Default constructor.
     */
    JPetzhold()
    {
      //       angle	probability
      //       [deg]
      (*this)[  0.100] = 1.767e+03;
      (*this)[  0.126] = 1.296e+03;
      (*this)[  0.158] = 9.502e+02;
      (*this)[  0.200] = 6.991e+02;
      (*this)[  0.251] = 5.140e+02;
      (*this)[  0.316] = 3.764e+02;
      (*this)[  0.398] = 2.763e+02;
      (*this)[  0.501] = 2.188e+02;
      (*this)[  0.631] = 1.444e+02;
      (*this)[  0.794] = 1.022e+02;
      (*this)[  1.000] = 7.161e+01;
      (*this)[  1.259] = 4.958e+01;
      (*this)[  1.585] = 3.395e+01;
      (*this)[  1.995] = 2.281e+01;
      (*this)[  2.512] = 1.516e+01;
      (*this)[  3.162] = 1.002e+01;
      (*this)[  3.981] = 6.580e+00;
      (*this)[  5.012] = 4.295e+00;
      (*this)[  6.310] = 2.807e+00;
      (*this)[  7.943] = 1.819e+00;
      (*this)[ 10.000] = 1.153e+00;
      (*this)[ 15.000] = 4.893e-01;
      (*this)[ 20.000] = 2.444e-01;
      (*this)[ 25.000] = 1.472e-01;
      (*this)[ 30.000] = 8.609e-02;
      (*this)[ 35.000] = 5.931e-02;
      (*this)[ 40.000] = 4.210e-02;
      (*this)[ 45.000] = 3.067e-02;
      (*this)[ 50.000] = 2.275e-02;
      (*this)[ 55.000] = 1.699e-02;
      (*this)[ 60.000] = 1.313e-02;
      (*this)[ 65.000] = 1.046e-02;
      (*this)[ 70.000] = 8.488e-03;
      (*this)[ 75.000] = 6.976e-03;
      (*this)[ 80.000] = 5.842e-03;
      (*this)[ 85.000] = 4.953e-03;
      (*this)[ 90.000] = 4.292e-03;
      (*this)[ 95.000] = 3.782e-03;
      (*this)[100.000] = 3.404e-03;
      (*this)[105.000] = 3.116e-03;
      (*this)[110.000] = 2.912e-03;
      (*this)[115.000] = 2.797e-03;
      (*this)[120.000] = 2.686e-03;
      (*this)[125.000] = 2.571e-03;
      (*this)[130.000] = 2.476e-03;
      (*this)[135.000] = 2.377e-03;
      (*this)[140.000] = 2.329e-03;
      (*this)[145.000] = 2.313e-03;
      (*this)[150.000] = 2.365e-03;
      (*this)[155.000] = 2.506e-03;
      (*this)[160.000] = 2.662e-03;
      (*this)[165.000] = 2.835e-03;
      (*this)[170.000] = 3.031e-03;
      (*this)[175.000] = 3.092e-03;
      (*this)[180.000] = 3.154e-03;

      this->mul(1.166);  // normalisation
      this->compile();
    }


    /**
     * Measurement of light scattering in water.
     *
     * \param  x      cosine scattering angle
     * \return        probability
     */
    double operator()(const double x) const
    {
      double a;

      if      (x >= +1.0)
	a =   0.0;
      else if (x <= -1.0)
	a = 180.0;
      else
	a = acos(x) * 180.0 / PI;

      if      (a >= this->getXmax()) 
	return this->rbegin()->getY();
      else if (a <= this->getXmin()) 
	return this-> begin()->getY();
      else
	return JPolint1Function1D_t::operator()(a);
    }
  };
  

  /**
   * Function object for measurement of light scattering in water.
   */
  static const JPetzhold petzhold;


  /**
   * Auxiliary data structure for scattering lengths of deep-sea water.
   *
   * Use the Kopelevich model for spectral volume scattering functions.\n
   * The model separates the contributions by:
   * - pure:   pure sea water;
   * - small: 'small' particles (size < 1 um);
   * - large: 'large' particles (size > 1 um).
   *
   * Values are taken from reference C.D. Mobley "Light and Water", ISBN 0-12-502750-8, pag. 119.
   */
  struct JMobley {
    /**
     * Constructor.
     *
     * \param  lambda     wavelength of light [nm]
     */
    JMobley(const double lambda)
    {
      static const double Vs = 0.0075;
      static const double Vl = 0.0075;
      static const double bw = 0.0017;
      static const double bs = 1.340;
      static const double bl = 0.312;
      
      const double x = 550.0/lambda;

      Lis_pure  = bw *      pow(x, 4.3);
      Lis_small = bs * Vs * pow(x, 1.7);
      Lis_large = bl * Vl * pow(x, 0.3);
    }

    double Lis_pure;       //!< inverse scattering length due to pure water      [m^-1]
    double Lis_small;      //!< inverse scattering length due to small particles [m^-1]
    double Lis_large;      //!< inverse scattering length due to large particles [m^-1]
  };


  /**
   * Absorption length of pure water.
   *
   * CITATION:
   * Jonasz M. 2007. Absorption coefficient of water: Data sources (www.tpdsci.com/Tpc/AbsCfOfWaterDat.php).
   * In: Top. Part. Disp. Sci. (www.tpdsci.com). 			
   *
   * DATA FROM:
   * Wozniak B., Wozniak S. B., Tyszka K., Ostrowska M., Majchrowski R., Ficek D., Dera J. 2005.
   * Modelling the light absorption properties of particulate matter forming organic particles suspended in seawater. Part 2.
   * Modelling results. Oceanologia 47, 621-662.
   * see also
   * Wozniak B., Dera J. 2007.
   * Light absorption in sea water. Springer, Berlin, 456 pp. (see p. 62)
   *
   * NOTES:
   * As stated by the data authors, the data are based on measurement results obtained by various authors 
   * (interpolated by a linear approximation where applicable):
   * Wavelength	Reference
   * - 200-335 nm   Smith R.C., Baker K.S. 1981. Optical properties of the clearest natural waters (200-800 nm). Appl. Opt. 20, 177-184.
   * - 340-370 nm   Sogandares F.M., Fry, E.S. 1997. Absorption spectrum (340 -640 nm) of pure water. I. Photothermal measurements Appl. Opt. 36, 8699-8709.
   * - 380-700 nm   Pope R.M., Fry E.S. 1997. Absorption spectrum (380 -700 nm) of pure water. II. Integrating cavity measurements. Appl. Opt. 36, 8710-8723
   */
  class JAbsorptionLengthOfPureWater :
    public JGridSplineFunction1D_t
  {
  public:
    /**
     * Default constructor.
     */
    JAbsorptionLengthOfPureWater() 
    {
      //    wave-	absorption
      //    length	coefficient
      //    [um]	[1/m]
      (*this)[0.200e3]  =  3.07;
      (*this)[0.205e3]  =  2.53;
      (*this)[0.210e3]  =  1.99;
      (*this)[0.215e3]  =  1.65;
      (*this)[0.220e3]  =  1.31;
      (*this)[0.225e3]  =  1.1185;
      (*this)[0.230e3]  =  0.927;
      (*this)[0.235e3]  =  0.8235;
      (*this)[0.240e3]  =  0.72;
      (*this)[0.245e3]  =  0.6395;
      (*this)[0.250e3]  =  0.559;
      (*this)[0.255e3]  =  0.508;
      (*this)[0.260e3]  =  0.457;
      (*this)[0.265e3]  =  0.415;
      (*this)[0.270e3]  =  0.373;
      (*this)[0.275e3]  =  0.3305;
      (*this)[0.280e3]  =  0.288;
      (*this)[0.285e3]  =  0.2515;
      (*this)[0.290e3]  =  0.215;
      (*this)[0.295e3]  =  0.178;
      (*this)[0.300e3]  =  0.141;
      (*this)[0.305e3]  =  0.123;
      (*this)[0.310e3]  =  0.105;
      (*this)[0.315e3]  =  0.0947;
      (*this)[0.320e3]  =  0.0844;
      (*this)[0.325e3]  =  0.0761;
      (*this)[0.330e3]  =  0.0678;
      (*this)[0.335e3]  =  0.06195;
      (*this)[0.340e3]  =  0.0325;
      (*this)[0.345e3]  =  0.02645;
      (*this)[0.350e3]  =  0.0204;
      (*this)[0.355e3]  =  0.018;
      (*this)[0.360e3]  =  0.0156;
      (*this)[0.365e3]  =  0.0135;
      (*this)[0.370e3]  =  0.0114;
      (*this)[0.375e3]  =  0.011385;
      (*this)[0.380e3]  =  0.01137;
      (*this)[0.385e3]  =  0.00941;
      (*this)[0.390e3]  =  0.00851;
      (*this)[0.395e3]  =  0.00813;
      (*this)[0.400e3]  =  0.00663;
      (*this)[0.405e3]  =  0.0053;
      (*this)[0.410e3]  =  0.00473;
      (*this)[0.415e3]  =  0.00444;
      (*this)[0.420e3]  =  0.00454;
      (*this)[0.425e3]  =  0.00478;
      (*this)[0.430e3]  =  0.00495;
      (*this)[0.435e3]  =  0.0053;
      (*this)[0.440e3]  =  0.00635;
      (*this)[0.445e3]  =  0.00751;
      (*this)[0.450e3]  =  0.00922;
      (*this)[0.455e3]  =  0.00962;
      (*this)[0.460e3]  =  0.00979;
      (*this)[0.465e3]  =  0.01011;
      (*this)[0.470e3]  =  0.0106;
      (*this)[0.475e3]  =  0.0114;
      (*this)[0.480e3]  =  0.0127;
      (*this)[0.485e3]  =  0.0136;
      (*this)[0.490e3]  =  0.015;
      (*this)[0.495e3]  =  0.0173;
      (*this)[0.500e3]  =  0.0204;
      (*this)[0.505e3]  =  0.0256;
      (*this)[0.510e3]  =  0.0325;
      (*this)[0.515e3]  =  0.0396;
      (*this)[0.520e3]  =  0.0409;
      (*this)[0.525e3]  =  0.0417;
      (*this)[0.530e3]  =  0.0434;
      (*this)[0.535e3]  =  0.0452;
      (*this)[0.540e3]  =  0.0474;
      (*this)[0.545e3]  =  0.0511;
      (*this)[0.550e3]  =  0.0565;
      (*this)[0.555e3]  =  0.0596;
      (*this)[0.560e3]  =  0.0619;
      (*this)[0.565e3]  =  0.0642;
      (*this)[0.570e3]  =  0.0695;
      (*this)[0.575e3]  =  0.0772;
      (*this)[0.580e3]  =  0.0896;
      (*this)[0.585e3]  =  0.11;
      (*this)[0.590e3]  =  0.1351;
      (*this)[0.595e3]  =  0.1672;
      (*this)[0.600e3]  =  0.2224;
      (*this)[0.605e3]  =  0.2577;
      (*this)[0.610e3]  =  0.2644;
      (*this)[0.615e3]  =  0.2678;
      (*this)[0.620e3]  =  0.2755;
      (*this)[0.625e3]  =  0.2834;
      (*this)[0.630e3]  =  0.2916;
      (*this)[0.635e3]  =  0.3012;
      (*this)[0.640e3]  =  0.3108;
      (*this)[0.645e3]  =  0.325;
      (*this)[0.650e3]  =  0.34;
      (*this)[0.655e3]  =  0.371;
      (*this)[0.660e3]  =  0.41;
      (*this)[0.665e3]  =  0.429;
      (*this)[0.670e3]  =  0.439;
      (*this)[0.675e3]  =  0.448;
      (*this)[0.680e3]  =  0.465;
      (*this)[0.685e3]  =  0.486;
      (*this)[0.690e3]  =  0.516;
      (*this)[0.695e3]  =  0.559;
      (*this)[0.700e3]  =  0.624;

      compile();
    }


    /**
     * Absorption length of pure water.
     *
     * \param  lambda     wavelength of light [nm]
     * \return            absorption length   [m]
     */
    double operator()(const double lambda) const
    {
      const double y = JGridSplineFunction1D_t::operator()(lambda);
      
      return 1.0 / y;
    }
  };
  

  /**
   * Function object for absorption length of pure water.
   */
  static const JAbsorptionLengthOfPureWater getAbsorptionLengthOfPureWater;
}

#endif