#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 "PndDrcOptMatLLF1.h"


//----------------------------------------------------------------------
PndDrcOptMatLLF1::PndDrcOptMatLLF1()
{
  m_B1 = 1.21640125;
  m_B2 = 0.13366454;
  m_B3 = 0.883399468;
  m_C1 = 0.00857807248;
  m_C2 = 0.0420143003;
  m_C3 = 107.59306;

}
//----------------------------------------------------------------------
PndDrcOptMatLLF1* PndDrcOptMatLLF1::clone() const
{
  return new PndDrcOptMatLLF1(*this);
}
//----------------------------------------------------------------------
void PndDrcOptMatLLF1::copy(const PndDrcOptMatLLF1& 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;
}//----------------------------------------------------------------------
PndDrcOptMatLLF1::PndDrcOptMatLLF1(const PndDrcOptMatLLF1& mat)
  : PndDrcOptMatAbs(mat)
{
  if (mat.m_verbosity>=1) cout<<"  PndDrcOptMatLLF1::PndDrcOptMatLLF1"
			    <<"(const PndDrcOptMatLLF1&) "  
			    <<mat.m_name<<endl;
  copy(mat);
} 
//----------------------------------------------------------------------
PndDrcOptMatLLF1& PndDrcOptMatLLF1::operator=(const PndDrcOptMatLLF1& mat)
{
  if (mat.m_verbosity>=1) cout<<"  PndDrcOptMatLLF1::operator="
			    <<"(const PndDrcOptMatLLF1&) "  
			    <<mat.m_name<<endl;
  if (&mat != this)
    {
      static_cast<PndDrcOptMatAbs&>((*this)) = mat; // assignment of base class part.
      copy(mat);
    }
  return *this;
}
//----------------------------------------------------------------------
double PndDrcOptMatLLF1::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 PndDrcOptMatLLF1::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 PndDrcOptMatLLF1::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,  5000,  5000,  5000,
				 5000,   3328,  3328,  1995,  1662,
				 1245,  550.5, 118.4, 20.78, 7.007,
				2.556};
  // C = -10mm /ln (t_i)

  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<<" *** PndDrcOptMatLLF1::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.
}