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


//----------------------------------------------------------------------
PndDrcOptMatLiF::PndDrcOptMatLiF()
{
  double temp=300; // Kelvin
   
  calcSellmeier(m_B1,m_B2,m_B3,m_C1,m_C2,m_C3,temp);
  
  //cout<<" PndDrcOptMatLiF::PndDrcOptMatLiF: check: n(400)="<<refIndex(400)<<endl;
  
}
//----------------------------------------------------------------------
PndDrcOptMatLiF* PndDrcOptMatLiF::clone() const
{
  return new PndDrcOptMatLiF(*this);
}
//----------------------------------------------------------------------
void PndDrcOptMatLiF::copy(const PndDrcOptMatLiF& 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;
}//----------------------------------------------------------------------
PndDrcOptMatLiF::PndDrcOptMatLiF(const PndDrcOptMatLiF& mat)
  : PndDrcOptMatAbs(mat)
{
  if (mat.m_verbosity>=1) cout<<"  PndDrcOptMatLiF::PndDrcOptMatLiF"
			    <<"(const PndDrcOptMatLiF&) "  
			    <<mat.m_name<<endl;
  copy(mat);
} 
//----------------------------------------------------------------------
PndDrcOptMatLiF& PndDrcOptMatLiF::operator=(const PndDrcOptMatLiF& mat)
{
  if (mat.m_verbosity>=1) cout<<"  PndDrcOptMatLiF::operator="
			    <<"(const PndDrcOptMatLiF&) "  
			    <<mat.m_name<<endl;
  if (&mat != this)
    {
      static_cast<PndDrcOptMatAbs&>((*this)) = mat; // assignment of base class part.
      copy(mat);
    }
  return *this;
}
//----------------------------------------------------------------------
double PndDrcOptMatLiF::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 PndDrcOptMatLiF::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 PndDrcOptMatLiF::absorptionFlag(double lambda, double length) const
{
  // not yet implemented
  return false; // no absorption.
}
//----------------------------------------------------------------------
void PndDrcOptMatLiF::calcSellmeier(double& b1, double& b2, double& b3,
				 double& c1, double& c2, double& c3,
				 double temp)
{
  const double s1[5] = {9.33385e-1,1.84043e-5,-1.68868e-7,4.08865e-12,3.36861e-13};
  const double s2[5] = {2.91928,6.91076e-3,8.47693e-6,-3.20709e-8,1.36856e-10};
  const double s3[5] = {2.73548,6.71841e-4,-1.24327e-5,-2.97405e-8,1.21434e-10};
  
  const double l1[5] = {7.15768e-2,8.52651e-6,-5.56075e-8,2.37456e-10,-3.54594e-13};
  const double l2[5] = {2.93086e1,3.64286e-2,3.22167e-5,7.85825e-8,-2.81432e-10};
  const double l3[5] = {2.94083e1,1.04614e-4,-8.56271e-5,-1.18786e-7,7.82576e-10};// nm
  
  b1=0;
  b2=0;
  b3=0;
  c1=0;
  c2=0;
  c3=0;
  for (int i=0; i<5; i++)
    {
      b1 += s1[i]*pow(temp,i);
      b2 += s2[i]*pow(temp,i);
      b3 += s3[i]*pow(temp,i);
      c1 += l1[i]*pow(temp,i);
      c2 += l2[i]*pow(temp,i);
      c3 += l3[i]*pow(temp,i);
    }
  c1 = c1*c1;
  c2 = c2*c2;
  c3 = c3*c3;
}