#include "SiRDPTrap.h"

ClassImp(SiRDPTrap)

 void SiRDPTrap::SetDefName(Char_t *x) 
{
  if(strlen(x)<10) name=x;
  else 
    printf("Name must be shorter than 10 chars!\n");
}

 void SiRDPTrap::SetNoSteps(Int_t x) 
{
  if(x<1 || x>100000) 
    {
    printf("Number of step is out of bounds!!\n");
    return;
    }
  Na.Set(x);
  Nb.Set(x);
  Nc.Set(x);
}


 void SiRDPTrap::SetModel(Int_t x) 
{
  if(x<1 || x>2) 
    {
    printf("Only models of 1. and 2. order are allowed!\n");
    return;
    }
  else Model=x;
}

 void SiRDPTrap::SetGenRate(Float_t *x) 
{
  for(Int_t i=0;i<Model;i++)
  if(x[i]<0) 
    {
    printf("Generation rate of defects must be a positive number!\n");
    return;
    }
  else gr[i]=x[i];
}

 void SiRDPTrap::SetTref(Float_t x) 
{
  if(x<0 || x> 1000) 
    {
    printf("Reference temperature must be ]0,1000[ !\n");
    return;
    }
  else T_ref=x;
}

 void SiRDPTrap::SetEa(Float_t x) 
{
  if(x<0 || x>5) 
    {
    printf("Activation energy must be ]0,5[ !\n");
    return;
    }
  else Ea=x;
}

 void SiRDPTrap::SetReacConst(Float_t x) 
{
  if(x<0) 
    {
    printf("Reaction constant must be a positive number\n");
    return;
    }
  else kr_Tref=x;
}



 void SiRDPTrap::rk4(float y[], float dydx[], int n, float x, float h, float yout[])
{
	int i;
	float xh,hh,h6,*dym,*dyt,*yt;

	dym=new float[n];
	dyt=new float[n];
	yt=new float[n];
	hh=h*0.5;
	h6=h/6.0;
	xh=x+hh;
	for (i=0;i<n;i++) yt[i]=y[i]+hh*dydx[i];
	Dynamics(xh,yt,dyt);
	for (i=0;i<n;i++) yt[i]=y[i]+hh*dyt[i];
	Dynamics(xh,yt,dym);
	for (i=0;i<n;i++) {
		yt[i]=y[i]+h*dym[i];
		dym[i] += dyt[i];
	}
	Dynamics(x+h,yt,dyt);
	for (i=0;i<n;i++)
		yout[i]=y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]);
        delete [] dym;
	delete [] dyt;
	delete [] yt;
}

 void SiRDPTrap::Dynamics(Float_t x,Float_t y[],Float_t dydx[])
{
  switch(Model)
    {
    case 1:
    dydx[0]= -y[0]*kr+gr[0]*flux;
    dydx[1]=  y[0]*kr; 
    break;
    case 2:
    dydx[0]= -kr*y[0]*y[1]+gr[0]*flux;
    dydx[1]= -kr*y[0]*y[1]+gr[1]*flux;
    dydx[2]=  kr*y[0]*y[1];
    case 3:
    break;
    }
} 

 void SiRDPTrap::Reset(Int_t opt)
{
  switch(opt)
    {
    case 1:
      Na.Set(1000);
      Nb.Set(1000);
      if(Model==2) Nc.Set(1000);
    case 2:
      Model=1;
      gr[0]=0; gr[1]=0;
      T_ref=293;
      kr_Tref=0;
      flux=0; U=0; T=293; ts=1; kr=0;
      Na.Set(1000);
      Nb.Set(1000);
      if(Model==2) Nc.Set(1000); 
      else  Nc.Set(0); 
      break;
    default: 
      cs=0;
      Na.Reset(-1); Na[0]=0;
      Nb.Reset(-1); Nb[0]=0; 
      if(Model==2) {Nc.Reset(-1); Nc[0]=0;}
    }
  
}


 void SiRDPTrap::ScaleReacConst()
{
  Double_t kB=8.617385e-5;  // eV/K
  kr=kr_Tref*TMath::Exp(Ea/kB*(1/T_ref-1/T)); 
}



 void SiRDPTrap::DoStep(Float_t time_step)
{
  Float_t y[3],x=0,dydx[3],yout[3];
  if(time_step!=0) ts=time_step;
  //   if(cs<10) printf("Let's see: %d %f %fn",cs,Na[cs],Nb[cs]);
   if(Na[cs+1]>0) Na[cs]=Na[cs+1]; 
   if(Nb[cs+1]>0) Nb[cs]=Nb[cs+1];
   if(Model>1 && Nc[cs+1]>0) Nc[cs]=Nc[cs+1];

   
  y[0]=Na[cs]; 
  y[1]=Nb[cs];
  if(Model>1) y[2]=Nc[cs];
  Dynamics(x,y,dydx);
  rk4(y,dydx,Model+1,x,ts,yout);
  cs++;
  
    Na[cs]=yout[0];
    Nb[cs]=yout[1];
    if(Model>1) Nc[cs]=yout[2];
  
  
}

 void SiRDPTrap::Info()
{
  printf("Properties of the defect     : %s n",(const char *)name);
  printf("Model dynamcis of order      : %d n",Model);
  for(Int_t i=0;i<Model;i++)
  printf("Defect (%d) generation rate  : %6.3f n",i,gr[i]);
  printf("Reaction constant            : %6.3e n",kr_Tref);
  printf("Activation Energy            : %6.3f n",Ea);
  printf("Reference temperature        : %4.2f n",T_ref);
  printf("Contribution to space charge : A(%4.2f) B(%4.2f) C(%4.2f) n",NaVfd,NbVfd,NcVfd);
  printf("Defect sensitive to bias     : %d n",BiasSens);
  printf("Current settings applied to detector:n");
  printf("t time step  : %6.3f n",ts);
  printf("t flux       : %5.1f n",flux);	 
  printf("t bias       : %4.2f n",U);	 
  printf("t temperature: %4.2f n",T);	 
  printf("t reaction c.: %6.3e n",kr);	 

}

 SiRDPTrap::SiRDPTrap(Int_t x0,Float_t x1a,Float_t x1b,Float_t x2,Float_t x3,Float_t x4,Char_t *dname)
{
  Float_t g[2];
  g[0]=x1a;
  g[1]=x1b;
  SetModel(x0);
  SetGenRate(g); 
  SetReacConst(x2);
  SetEa(x3);
  SetTref(x4);
  Na.Set(1000); Na.Reset(-1); Na[0]=0;
  Nb.Set(1000); Nb.Reset(-1); Nb[0]=0;
  if(Model==2) 
    { Nc.Set(1000); Nc.Reset(-1); Nc[0]=0;}
  SetDefName(dname);
  flux=100;
  U=0;
  SetT(293);
  ts=1;
  BiasSens=kFALSE;
  cs=0;  
}

 TGraph *SiRDPTrap::TimePlot(Option_t *what,Float_t *tim)
{
  Int_t i=0;
  TString opt=what;
  opt.ToLower();
  Float_t *time;
  Float_t *plot;
  if(tim==NULL)
    {
      time=new Float_t [cs];
      for(i=0;i<cs;i++) time[i]=i*ts;
    }
  else
  time=tim;

  if (opt.Contains("na")) plot=Na.GetArray();
  if (opt.Contains("nb")) plot=Nb.GetArray();
  if (opt.Contains("nc")) plot=Nc.GetArray();
  

  TGraph *gr=new TGraph(cs,time,plot);
  gr->SetMarkerStyle(21);
  
  delete [] time;
return gr;
}


 SiRDPTrap::~SiRDPTrap()
{

}

 void SiRDPTrap::Copy(SiRDPTrap *dest)
{
 Int_t i;
 dest->name=name;  
 dest->Model=Model;
 for(i=0;i<Model;i++) dest->gr[i]=gr[i];
dest->T_ref=T_ref;
dest->kr_Tref=kr_Tref;
dest->Ea=Ea;
dest->BiasSens=BiasSens;
dest->cs=cs;
dest->flux=flux;
dest->U=U;
dest->T=T;
dest->ts=ts;
dest->kr=kr;
dest->NaVfd=NaVfd;
dest->NbVfd=NbVfd;
dest->NcVfd=NcVfd;

dest->Na.Set(Na.GetSize());
dest->Nb.Set(Nb.GetSize());
dest->Nc.Set(Nc.GetSize());

 for(i=0;i<Na.GetSize();i++) dest->Na[i]=Na[i];
 for(i=0;i<Nb.GetSize();i++) dest->Nb[i]=Nb[i];
 for(i=0;i<Nc.GetSize();i++) dest->Nc[i]=Nc[i];

}


 Int_t SiRDPTrap::ReadTrapFile(FILE *in)
{
  Char_t Tname[10];
  Int_t x1,x9,ok;
  Float_t x2[2],x3,x4,x5,x6,x7,x8;
  ok=fscanf(in,"%s %d %f %f %e %f %f %f %f %f %d",Tname,&x1,&x2[0],&x2[1],&x3,&x4,&x5,&x6,&x7,&x8,&x9);
  if(ok==11) 
    {
  SetDefName(Tname);
  SetModel(x1);
  SetGenRate(x2); 
  SetReacConst(x3);
  SetEa(x4);
  SetTref(x5); 
  NaVfd=x6;
  NbVfd=x7;
  NcVfd=x8;
  SetBiasSens((Bool_t) x9);
  return 0;
   }
  else
    return 1;
}