#include <iostream>
#include <random>
#include <sys/time.h>

using namespace std;

struct vec{
  double x;
  double y;
  double z;
  void set(double fx, double fy, double fz);
  vec(double fx=0, double fy=0, double fz=0){set(fx,fy,fz);}
  vec& operator=(const vec& rhs){
    set(rhs.x,rhs.y,rhs.z);
    return *this;
  }
};

void vec::set(double fx, double fy, double fz){
  x=fx;
  y=fy;
  z=fz;
}

vec operator+(vec v1, vec v2){
  vec v;
  v.x=v1.x+v2.x;
  v.y=v1.y+v2.y;
  v.z=v1.z+v2.z;
  return v;
}

vec operator*(double a, vec v1){
  vec v;
  v.x=a*v1.x;
  v.y=a*v1.y;
  v.z=a*v1.z;
  return v;
}

vec cMult(vec a, vec b){
  vec v;
  v.x=a.x*b.x;
  v.y=a.y*b.y;
  v.z=a.z*b.z;
  return v;
}

double mag(vec a){
  return a.x+a.y+a.z;
}

double mag2(vec a){
  vec v=cMult(a,a);
  return mag(v);
}

double dot(vec a, vec b){
  vec v=cMult(a,b);
  return mag(v);
}

std::ostream& operator <<(std::ostream& os, vec& v){
  return os<<v.x<<' '<<v.y<<' '<<v.z;
}

double cThetaFromE(double e, double mu){
  return 1-(1/e-1)/mu;
}

double rnd(std::mt19937& gen){
  return (gen()-mt19937::min())/(mt19937::max()-mt19937::min()+0.);

}

void sampleState(std::mt19937& gen, double mu, double& e){

    //compton edge
    double e0=1./(1+2*mu);

    //integrals for f1 and f2
    double f1=-log(e0);
    double f2=0.5*(1-e0*e0);
    double q1=f1/(f1+f2);
    while (1){
      double x1=rnd(gen);
      double x2=rnd(gen);
      if (x1<q1)
        //generate according to 1/e
        e=pow(e0,x2);

      else
        //generate according to e
        e=sqrt(e0*e0*(1-x2)+x2);

      //apply g
      double x3=rnd(gen);
      double cTh=cThetaFromE(e,mu);
      double sTh2=1-cTh*cTh;
      double g=1-e*sTh2/(1+e*e);
      if (x3<g) return;
    }

}

vec randomDirection(std::mt19937& rnd){



}

vec getR1(std::mt19937& gen, vec rc, vec fsize, vec& r0, vec& ri, double& ctheta){
  double pi=3.14159;
  double phi=2*pi*rnd(gen);
  double sphi=sin(phi);
  double cphi=cos(phi);
  ctheta=rnd(gen);//only fwd
  double stheta=sqrt(1-ctheta*ctheta);
  vec v(cphi*stheta,sphi*stheta,ctheta);
  vec eps(2*rnd(gen)-1,2*rnd(gen)-1,2*rnd(gen)-1);
  vec feps=cMult(fsize,eps);
  r0=rc+feps;
  double len1=-r0.z/ctheta;
  ri=r0+len1*v;
  double len=-log(rnd(gen))+len1;
  vec r1=r0+len*v;
  return r1;
}

vec getR1Cube(std::mt19937& gen, vec rc, vec fsize, vec& r0, vec& ri, double& ctheta){
    vec eps(2*rnd(gen)-1,2*rnd(gen)-1,2*rnd(gen)-1);
    vec feps=cMult(fsize,eps);
    r0=rc+feps;


}

int main(int argc, char **argv){
  std::mt19937 mtGen;
  timeval t;
  gettimeofday(&t,0);
  mtGen.seed(t.tv_usec);
  vec rc(0,0,-10);
  vec fsize(10,10,10);
  double ctheta;
  vec r0;
  vec ri;
  int n=1000;
  if (argc>1) n=atoi(argv[1]);

  double mu=0.5;
  if (argc>2) mu=atof(argv[2]);
	
  double s=0;
  for (int i=0;i<n;i++){
    double e;
    //vec r1=getR1(mtGen,rc,fsize,r0,ri,ctheta);
    sampleState(mtGen,mu,e);
    //cout<<r0<<' '<<ri<<' '<<r1<<' '<<ctheta<<' '<<mag2(r0+(-1)*r1)<<endl;
    cout<<e<<endl;
    s+=e;
  }
  cout<<"Avg: "<<s/n<<endl;
  return 0;



}