#include <iostream>
using std::cout;
using std::cerr;
using std::cin;
using std::endl;

#include <valarray>
using std::valarray;

#include <fstream>
using std::fstream;

#include <string>
using std::string;

#include <list>
using std::list;

#include "TObject.h"
#include "TVector3.h"
#include "TRandom.h"

#include "PndDrcOptMatAbs.h"
#include "PndDrcOptMatNLAK33A.h"


//----------------------------------------------------------------------
PndDrcOptMatNLAK33A::PndDrcOptMatNLAK33A()
{
  m_B1 = 1.44116999;
  m_B2 = 0.571749501;
  m_B3 = 1.16605226;
  m_C1 = 0.00680933877;
  m_C2 = 0.0222291824;
  m_C3 = 80.9379555;

}
//----------------------------------------------------------------------
PndDrcOptMatNLAK33A* PndDrcOptMatNLAK33A::clone() const
{
  return new PndDrcOptMatNLAK33A(*this);
}
//----------------------------------------------------------------------
void PndDrcOptMatNLAK33A::copy(const PndDrcOptMatNLAK33A& mat)
{
  m_B1  = mat.m_B1;
  m_C1  = mat.m_C1;
  m_B2  = mat.m_B2;
  m_C2  = mat.m_C2;
  m_B3  = mat.m_B3;
  m_C3  = mat.m_C3;
  m_ran = mat.m_ran;
}//----------------------------------------------------------------------
PndDrcOptMatNLAK33A::PndDrcOptMatNLAK33A(const PndDrcOptMatNLAK33A& mat)
  : PndDrcOptMatAbs(mat)
{
  if (mat.m_verbosity>=1) cout<<"  PndDrcOptMatNLAK33A::PndDrcOptMatNLAK33A"
			    <<"(const PndDrcOptMatNLAK33A&) "  
			    <<mat.m_name<<endl;
  copy(mat);
} 
//----------------------------------------------------------------------
PndDrcOptMatNLAK33A& PndDrcOptMatNLAK33A::operator=(const PndDrcOptMatNLAK33A& mat)
{
  if (mat.m_verbosity>=1) cout<<"  PndDrcOptMatNLAK33A::operator="
			    <<"(const PndDrcOptMatNLAK33A&) "  
			    <<mat.m_name<<endl;
  if (&mat != this)
    {
      static_cast<PndDrcOptMatAbs&>((*this)) = mat; // assignment of base class part.
      copy(mat);
    }
  return *this;
}
//----------------------------------------------------------------------
double PndDrcOptMatNLAK33A::refIndex(const double lambda) const
{

  if (lambda<0) return 1.74; // average value.

  double lam2 = lambda/1000 * lambda/1000; // um2
  
  return sqrt(1.0L + 
	      m_B1*lam2/(lam2-m_C1) +
	      m_B2*lam2/(lam2-m_C2) +
	      m_B3*lam2/(lam2-m_C3));
}  
//----------------------------------------------------------------------
double PndDrcOptMatNLAK33A::refIndexDeriv(const double lambda) const
{
  double lam  = lambda/1000;
  double lam2 = lam*lam;
  //double lam3 = lam2*lam;

  return (   (-m_B1*m_C1*lam)/((lam2-m_C1)*(lam2-m_C1)) +
	     (-m_B2*m_C2*lam)/((lam2-m_C2)*(lam2-m_C2)) +
	     (-m_B3*m_C3*lam)/((lam2-m_C3)*(lam2-m_C3))  ) 
    / refIndex(lambda) / 1000;

}
//----------------------------------------------------------------------
bool PndDrcOptMatNLAK33A::absorptionFlag(double lambda, double length) const
{

  // Rayleigh scattering. 
  // data from Schott data sheets of 10mm sample

  const static double lam[21] = {1060,    700,   650,   620,   580,
				 546,     500,   460,   436,   420,
				 405,     400,   390,   380,   370,
				 365,     350,   334,   320,   310,
				 300};
  const static double C[21]   = {5000,   5000,  5000,  5000,  5000,
				 5000,   5000,  1661,  1106, 828.3,
				 521.3, 411.6, 298.0, 195.0, 126.5,
				 100.2, 45.83, 19.64, 9.169, 5.457,
				 3.338};

  double clarity;

  if (lambda>1060)
    {
      clarity=5000;
    }
  else if (lambda<300)
    {
      return true; // cut off
    }
  else
    {
      // find right bin
      int ibin=-1;
      for (int i=1; i<21; i++)
	{
	  if (lambda<lam[i-1] && lambda>=lam[i])
	    {
	      ibin = i;
	    }
	}
      if (ibin==-1)
	{
	  cerr<<" *** PndDrcOptMatNLAK33A::absorptionFlag: "
	      <<"this line should never been hit"<<endl;
	  exit(EXIT_FAILURE);
	}
      double frac = (lam[ibin-1]-lambda)/(lam[ibin-1]-lam[ibin]);
      clarity = C[ibin-1] + frac * (C[ibin]-C[ibin-1]);


    }

  // clarity at lambda!
  //double trans = exp(-(length)/(clarity*pow(lambda/663,4)));
  double trans = exp(-(length)/clarity);
  double cmp   = m_ran.Uniform(1.0);


  if (cmp>trans)
    {
      return true; // absorbed
    } 

  return false; // no absorption.
}