#include <math.h>
#include <iostream>
void EField( int model, double VB, double VD, double x, double y, double &Ex, 
double &Ey );
double mud_e(double T, double E) ;
double r_e(double T) ;
double mud_h(double T, double E) ;
double r_h(double T) ;
double EyValue150(int& ix, int& iy), ExValue150(int& ix, int& iy);

bool electron(int model, double VB, double VD, double B, double T, double x_e, double y_e, double &vx_e, double &vy_e, double &D_e) 
{
double kB= 1.38E-23;       // [m^2*kg/s^2/K]
double e= 1.602E-19;       // [Coulomb]
double E, Ex, Ey, mu_e, v_e, tanLA_e, secLA, cosLA, sinLA;
EField(model, VB, VD, x_e, y_e, Ex, Ey);    // [V/cm]
if( Ey > 0.) {
   double REx = -Ex; // because electron has negative charge 
   double REy = -Ey; // because electron has negative charge
   E = sqrt(Ex*Ex+Ey*Ey);
   mu_e = mud_e(T, E);
   v_e = mu_e * E;
   tanLA_e = r_e(T) * mu_e * (-B) * 1.E-4;  // because e has negative charge
   secLA = sqrt(1.+tanLA_e*tanLA_e);
   cosLA=1./secLA;
   sinLA = tanLA_e / secLA;
   vy_e = v_e * (REy*cosLA - REx*sinLA)/E;
   vx_e = v_e * (REx*cosLA + REy*sinLA)/E;
   D_e = kB * T * mu_e/ e;  
   return true;
}  else return false;
}

bool hole(int model, double VB, double VD, double B, double T, double x_h, double y_h, double &vx_h, double &vy_h, double &D_h)
{
double kB= 1.38E-23;       // [m^2*kg/s^2/K]
double e= 1.602E-19;       // [Coulomb]
double E, Ex, Ey, mu_h, v_h, tanLA_h, secLA, cosLA, sinLA;
EField(model, VB, VD, x_h, y_h, Ex, Ey);  // [V/cm]
if( Ey > 0.) {
   E = sqrt(Ex*Ex+Ey*Ey);
   mu_h = mud_h(T, E);
   v_h = mu_h * E;
   tanLA_h = r_h(T) * mu_h * B * 1.E-4;
   secLA = sqrt(1.+tanLA_h*tanLA_h);
   cosLA=1./secLA;
   sinLA = tanLA_h / secLA;
   vy_h = v_h * (Ey*cosLA - Ex*sinLA)/E;
   vx_h = v_h * (Ex*cosLA + Ey*sinLA)/E;
   D_h = kB * T * mu_h/ e;  
   return true;
}  else return false;
}

double mud_e(double T, double E) {
   double vs_e=1.53E9*pow(T,-0.87);
   double Ec_e = 1.01*pow( T, 1.55);
   double beta_e = 2.57E-2* pow(T,0.66);
   return vs_e/Ec_e/pow(1+pow(E/Ec_e,beta_e),1./beta_e);
}

double mud_h(double T, double E) {
   double vs_h = 1.62E8 * pow( T, -0.52);
   double Ec_h = 1.24 * pow( T, 1.68);
   double beta_h = 0.46 * pow(T,0.17);
   return vs_h/Ec_h/pow(1+pow(E/Ec_h,beta_h),1./beta_h);
}

double r_e(double T) {
   return 1.13+0.0008*(T-273.15);
}

double r_h(double T) {
   return 0.72 - 0.0005*(T-273.15);
}

void EField( int model, double VB, double VD, double x, double y, double &Ex, 
double &Ey ) { 
// model==0; uniform E-field model
// model==1; flat diode model
// model==2; FEM solusions
// x ==0 : at the center of strip

double d = 0.0285;            // [cm]
double deltay=0.00025, deltax=0.00025, xpitch = 0.008 ;   // [cm]  
double xpedestal = xpitch*1000., xhalfpitch=xpitch/2.;
double wd, Emax;
double fx, fy, xx, xxx;
int ix, iy;
Ex=0.;
wd = sqrt(VB/VD) * d;
if (VB>VD) wd = d;
if( model ==0 ) Ey = VB/wd ;   // flat E field
else {
       if(model==1) {        // flat diode model
          if(VB > VD) Ey = (VB+VD)/d - 2.*VD*(d-y)/(d*d);
          else {
             Emax = 2.* wd * VD /(d*d);
             Ey = Emax*(1-(d-y)/wd);
          }
       }
       else {     // FEM analysis
         iy = int (y/deltay);
         fy = (y-iy*deltay) / deltay;
         xx = fmod (x+xpedestal, xpitch);
         if (xx <= xhalfpitch ) xxx = xx;
         else xxx = xpitch - xx; 
         ix = int(xxx/deltax) ;
         fx = (xxx - ix*deltax) / deltax;         
         //std::cout <<"x,y,ix,iy,fx,fy,xx,xxx= "<<x<<" "<<y<<" "<<ix
         //<<" "<<iy<<" "<<fx<<" "<<fy<<" "<<xx<<" "<<xxx<<std::endl;
         int ix1=ix+1;
         int iy1=iy+1;
         double Ex00 = 0., Ex10 = 0., Ex01 = 0., Ex11 = 0.;
         double Ey00 = 0., Ey10 = 0., Ey01 = 0., Ey11 = 0.;
 //-------- pick up appropriate data file for VD-----------------------
         Ex00 = ExValue150(ix,iy);  Ex10 = ExValue150(ix1,iy);
         Ex01 = ExValue150(ix,iy1); Ex11 = ExValue150(ix1,iy1);
         Ey00 = EyValue150(ix,iy);  Ey10 = EyValue150(ix1,iy);
         Ey01 = EyValue150(ix,iy1); Ey11 = EyValue150(ix1,iy1);
//---------------- end of data bank search---

         Ex = Ex00*(1.-fx)*(1.-fy) + Ex10*fx*(1.-fy)
                + Ex01*(1.-fx)*fy + Ex11*fx*fy ;
         //std::cout <<"xx, xhalfpitch = "<< xx << " "<< xhalfpitch<<std::endl;
         if(xx > xhalfpitch ) Ex = -Ex;
         Ey = Ey00*(1.-fx)*(1.-fy) + Ey10*fx*(1.-fy)
                + Ey01*(1.-fx)*fy + Ey11*fx*fy ;
       }
}
return;
}