// Simulation of a laser beam with various waist size and the possibility to create different interference patterns


//Size of the cube for the simulations
#define SPACE_SIZE 50

//Main code
Int_t laser(){
  //Unit is m
  Double_t wavelength=532.E-9;
  
  Double_t laser_intensity[SPACE_SIZE][SPACE_SIZE][SPACE_SIZE];

  //The cube in which the simulation are done goes from -width to +width
  Double_t x_width=5E-6;
  Double_t y_width=5E-6;
  Double_t z_width=5E-6;

  Double_t waist_x=5.E-6;
  Double_t waist_y=5.E-6;

  Double_t laserMax=0;
  Double_t laserMin=1.E+999;

  TH3D *laserh=new TH3D("laserh","Laser",SPACE_SIZE,(-1*x_width),x_width,SPACE_SIZE,(-1*y_width),y_width,SPACE_SIZE,(-1*z_width),z_width);
  TH2D *xSlice=new TH2D("hxslice","Slice along x=0",SPACE_SIZE,(-1*y_width),y_width,SPACE_SIZE,(-1*z_width),z_width);
  TH2D *ySlice=new TH2D("hyslice","Slice along y=0",SPACE_SIZE,(-1*z_width),z_width,SPACE_SIZE,(-1*x_width),x_width);
  TH2D *zSlice=new TH2D("hzslice","Slice along z=0",SPACE_SIZE,(-1*y_width),y_width,SPACE_SIZE,(-1*x_width),x_width);

  for (Int_t xloop=0;xloop<SPACE_SIZE;xloop++) {
    for (Int_t yloop=0;yloop<SPACE_SIZE;yloop++) {
      for (Int_t zloop=0;zloop<SPACE_SIZE;zloop++) {
	Double_t x_coord=table2coord(xloop,x_width);
	Double_t y_coord=table2coord(yloop,y_width);
	Double_t z_coord=table2coord(zloop,z_width);
	laser_intensity[xloop][yloop][zloop]=laserIntensity(x_coord,y_coord,z_coord,waist_x,waist_y,wavelength);
	laserh->SetBinContent(xloop,yloop,zloop,laser_intensity[xloop][yloop][zloop]);

	if (laser_intensity[xloop][yloop][zloop]>laserMax) {
	  laserMax=laser_intensity[xloop][yloop][zloop];
	}
	if (laser_intensity[xloop][yloop][zloop]<laserMin) {
	  laserMin=laser_intensity[xloop][yloop][zloop];
	}

	if (xloop==0){
	  xSlice->SetBinContent(yloop,zloop,laserIntensity(0.,y_coord,z_coord,waist_x,waist_y,wavelength));
	} //xSlice
	if (yloop==0){
	  ySlice->SetBinContent(zloop,xloop,laserIntensity(x_coord,0.,z_coord,waist_x,waist_y,wavelength));
	} //ySlice
	if (zloop==0){
	  zSlice->SetBinContent(yloop,xloop,laserIntensity(x_coord,y_coord,0.,waist_x,waist_y,wavelength));
	} //zSlice
      } //zloop
    }//yloop
    cout << "x = " << xloop << endl;
  }//xloop

  cout << "Laser min  " << laserMin << " / max: " << laserMax <<endl;

  drawLaser("laserIntensity",laserh,xSlice,ySlice,zSlice);

  //Phases
  TH3D *laserhPhase=new TH3D("laserhPhase","Laser",SPACE_SIZE,(-1*x_width),x_width,SPACE_SIZE,(-1*y_width),y_width,SPACE_SIZE,(-1*z_width),z_width);
  TH2D *xSlicePhase=new TH2D("hxslicePhase","Slice along x=0",SPACE_SIZE,(-1*y_width),y_width,SPACE_SIZE,(-1*z_width),z_width);
  TH2D *ySlicePhase=new TH2D("hyslicePhase","Slice along y=0",SPACE_SIZE,(-1*z_width),z_width,SPACE_SIZE,(-1*x_width),x_width);
  TH2D *zSlicePhase=new TH2D("hzslicePhase","Slice along z=0",SPACE_SIZE,(-1*y_width),y_width,SPACE_SIZE,(-1*x_width),x_width);


  for (Int_t xloop=0;xloop<SPACE_SIZE;xloop++) {
    for (Int_t yloop=0;yloop<SPACE_SIZE;yloop++) {
      for (Int_t zloop=0;zloop<SPACE_SIZE;zloop++) {

	Double_t x_coord=table2coord(xloop,x_width);
	Double_t y_coord=table2coord(yloop,y_width);
	Double_t z_coord=table2coord(zloop,z_width);
	Bool_t inside=(laser_intensity[xloop][yloop][zloop]>(laserMax/TMath::Sqrt(2)));
	laserhPhase->SetBinContent(xloop,yloop,zloop,getPhase(x_coord,y_coord,z_coord,inside,wavelength));
	if (xloop==0){
	  xSlicePhase->SetBinContent(yloop,zloop,getPhase(0.,y_coord,z_coord,(laser_intensity[(((Double_t)SPACE_SIZE)/2)][yloop][zloop]>(laserMax/TMath::Sqrt(2))),wavelength));
	} //xSlice
	if (yloop==0){
	  ySlicePhase->SetBinContent(zloop,xloop,getPhase(x_coord,0.,z_coord,(laser_intensity[xloop][(((Double_t)SPACE_SIZE)/2)][zloop]>(laserMax/TMath::Sqrt(2))),wavelength));
	} //ySlice
	if (zloop==0){
	  zSlicePhase->SetBinContent(yloop,xloop,getPhase(x_coord,y_coord,0.,(laser_intensity[xloop][yloop][(((Double_t)SPACE_SIZE)/2)]>(laserMax/TMath::Sqrt(2))),wavelength));
	} //zSlice
      } //zloop
    }//yloop
    cout << "x = " << xloop << endl;
  }//xloop

  drawLaser("laserPhase",laserhPhase,xSlicePhase,ySlicePhase,zSlicePhase);

  return 0;

} // laser class


Int_t coord2table(Double_t coord,Double_t width){
  return (((coord/(2*width))+0.5)+SPACE_SIZE);
} //coord2table



Double_t table2coord(Int_t table,Double_t width){
  return ((((Double_t)2*table/(Double_t)SPACE_SIZE)-1.0)*width);
} //table2coord



Double_t laserIntensity(Double_t x,Double_t y,Double_t z,Double_t waist_x,Double_t waist_y,Double_t lambda) {
  Double_t bx=((TMath::Pi())*waist_x*waist_x)/lambda;
  Double_t by=((TMath::Pi())*waist_y*waist_y)/lambda;

  //  cout << "bx: " << bx <<endl;

  Double_t wx=((1+((z/bx)*(z/bx)))*waist_x);
  Double_t wy=((1+((z/by)*(z/by)))*waist_y);

  //  cout << "bx: " << bx << " / "  << by << " wx: " << wx << " / " << wy << " " << (TMath::Exp(-1*(((x*x)/(wx*wx))+((y*y)/(wy*wy))))) << " @ " << (((x*x)/(wx*wx))+((y*y)/(wy*wy))) << " £ " << (((x*x)/(wx*wx))+((y*y)/(wy*wy))) << "  " << TMath::Exp(-1) << endl;

  Double_t ret=(TMath::Exp(-1*(((x*x)/(wx*wx))+((y*y)/(wy*wy)))));

  //  cout << " x: " << x << " y: " << y << " z: " << z << " ret= " << ret << endl;

  return ret;
} //laserIntensity

/*
Double_t getPhase(Double_t xcoord,Double_t ycoord,Double_t zcoord,Bool_t inside,Double_t wavelength){
  if (inside){
    Double_t distance=(TMath::Sqrt((xcoord*xcoord)+(ycoord*ycoord)+(zcoord*zcoord)));
    return TMath::Cos((2*TMath::Pi()*distance)/wavelength);
  } else {
    return -2;
  }
} //getPhase
*/

/*
//Shintake
Double_t getPhase(Double_t xcoord,Double_t ycoord,Double_t zcoord,Bool_t inside,Double_t wavelength){
  Double_t x1=0;
  Double_t y1=0;
  Double_t z1=-1E-3;

  Double_t x2=0;
  Double_t y2=0;
  Double_t z2=1E-3;
  if (inside){
    Double_t distance1=(TMath::Sqrt(((xcoord-x1)*(xcoord-x1))+((ycoord-y1)*(ycoord-y1))+((zcoord-z1)*(zcoord-z1))));
    Double_t phi1=TMath::Cos((2*TMath::Pi()*distance1)/wavelength);

    Double_t distance2=(TMath::Sqrt(((xcoord-x2)*(xcoord-x2))+((ycoord-y2)*(ycoord-y2))+((zcoord-z2)*(zcoord-z2))));
    Double_t phi2=TMath::Cos((2*TMath::Pi()*distance2)/wavelength);

    return (phi1-phi2);
  } else {
    return -5;
  }
} //getPhase
*/

/*
//Twyman-Green Interferometer
Double_t getPhase(Double_t xcoord,Double_t ycoord,Double_t zcoord,Bool_t inside,Double_t wavelength){
  if (inside){
    Double_t distance=(TMath::Sqrt((ycoord*ycoord)));
    return TMath::Cos((2*TMath::Pi()*distance)/wavelength);
  } else {
    return -2;
  }
} //getPhase
*/

/*
//Other two-sources
Double_t getPhase(Double_t xcoord,Double_t ycoord,Double_t zcoord,Bool_t inside,Double_t wavelength){
  Double_t x1=5E-3;
  Double_t y1=0;
  Double_t z1=-1E-3;

  Double_t x2=-5E-3;
  Double_t y2=0;
  Double_t z2=-1E-3;
  if (inside){
    Double_t distance1=(TMath::Sqrt(((xcoord-x1)*(xcoord-x1))+((ycoord-y1)*(ycoord-y1))+((zcoord-z1)*(zcoord-z1))));
    Double_t phi1=TMath::Cos((2*TMath::Pi()*distance1)/wavelength);

    Double_t distance2=(TMath::Sqrt(((xcoord-x2)*(xcoord-x2))+((ycoord-y2)*(ycoord-y2))+((zcoord-z2)*(zcoord-z2))));
    Double_t phi2=TMath::Cos((2*TMath::Pi()*distance2)/wavelength);

    return (phi1-phi2);
  } else {
    return -5;
  }
} //getPhase
*/


//Other 3-sources
Double_t getPhase(Double_t xcoord,Double_t ycoord,Double_t zcoord,Bool_t inside,Double_t wavelength){
  Double_t x1=5E-3;
  Double_t y1=0;
  Double_t z1=-1E-3;

  Double_t x2=0;
  Double_t y2=0;
  Double_t z2=-1E-3;

  Double_t x3=-5E-3;
  Double_t y3=0;
  Double_t z3=-1E-3;
  if (inside){
    Double_t distance1=(TMath::Sqrt(((xcoord-x1)*(xcoord-x1))+((ycoord-y1)*(ycoord-y1))+((zcoord-z1)*(zcoord-z1))));
    Double_t phi1=TMath::Cos((2*TMath::Pi()*distance1)/wavelength);

    Double_t distance2=(TMath::Sqrt(((xcoord-x2)*(xcoord-x2))+((ycoord-y2)*(ycoord-y2))+((zcoord-z2)*(zcoord-z2))));
    Double_t phi2=TMath::Cos((2*TMath::Pi()*distance2)/wavelength);

    Double_t distance3=(TMath::Sqrt(((xcoord-x3)*(xcoord-x3))+((ycoord-y3)*(ycoord-y3))+((zcoord-z3)*(zcoord-z3))));
    Double_t phi3=TMath::Cos((2*TMath::Pi()*distance2)/wavelength);

    return (phi1+phi2+phi3);
  } else {
    return -5;
  }
} //getPhase





Int_t drawLaser(TString filename,TH3D* laserh,TH2D* xSlice,TH2D* ySlice,TH2D* zSlice){

   gROOT->SetStyle("Plain");
   gStyle->SetOptStat(0);

   //gStyle->SetOptDate(1);
   //gStyle->SetOptStat(1111111);

   TString postscriptName=filename;
   postscriptName.Append(".ps");
   //select postscript output type
   Int_t type = 111;   //portrait  ps
   //Int_t type = 112;   //landscape ps
   // Int_t type = 113;   //eps

   //create a postscript file and set the paper size
   TPostScript ps(postscriptName,type);
   //ps.Range(26,19);  //set x,y of printed page
   ps.Range(19,26);  //set x,y of printed page
   TCanvas c1("c1","Cross section",900,1200);

   ps->NewPage();
  
   c1->Clear();
   c1->SetFillColor(0);
   
   gPad->SetLeftMargin(5.);
   gPad->SetTopMargin(2.);
   
   TPad *pad1=new TPad("pad1","CS",0.1,0.67,0.9,0.99);
   TPad *pad2=new TPad("pad2","CS",0.1,0.34,0.9,0.66);
   TPad *pad3=new TPad("pad2","CS",0.1,0.01,0.9,0.33);
   
   pad1->Draw();
   pad2->Draw();
   pad3->Draw();
   
   pad1->cd();
   pad1->SetLeftMargin(0.15);
   pad1->SetBottomMargin(0.15);

   //  gPad->SetLogx(1);
   //   pad1->SetLogx(1);
   //      pad1->SetLogy(1);
   

   laserh->Draw();

   pad2->cd();
   pad2->SetLeftMargin(0.15);
   pad2->SetBottomMargin(0.15);


   laserh->Project3D("y")->Draw();
   laserh->Project3D("x")->Draw("same");

   pad3->cd();
   pad3->SetLeftMargin(0.15);
   pad3->SetBottomMargin(0.15);

   //  gPad->SetLogx(1);
   // pad3->SetLogx(1);
   //         pad3->SetLogy(1);
   
   laserh->Project3D("z")->Draw();



   //  pad1->Print("electronCrossSection.gif");
  
  pad1->Update();
  pad2->Update();
  pad3->Update();
  
   ps->NewPage();
  
   c1->Clear();
   c1->SetFillColor(0);

   gPad->SetLeftMargin(5.);
   gPad->SetTopMargin(2.);
  
   c1->Divide(2,2);
   Int_t iframe=1;

   xSlice->GetXaxis()->SetTitle("Y (m)");
   xSlice->GetYaxis()->SetTitle("Z (m)");

   ySlice->GetXaxis()->SetTitle("Z (m)");
   ySlice->GetYaxis()->SetTitle("X (m)");

   zSlice->GetXaxis()->SetTitle("X (m)");
   zSlice->GetYaxis()->SetTitle("Y (m)");


   c1->cd(iframe++);
   ySlice->Draw("LEGO");
   c1->cd(iframe++);
   zSlice->Draw("LEGO");
   c1->cd(iframe++);
   c1->cd(iframe++);
   xSlice->Draw("LEGO");

   c1->Update();
   TString imgname=filename;
   imgname.Append("_LEGO.gif");

   c1->Print(imgname);

   ps->NewPage();
  
   c1->Clear();
   c1->SetFillColor(0);

   gPad->SetLeftMargin(5.);
   gPad->SetTopMargin(2.);
  
   c1->Divide(2,2);
   Int_t iframe=1;


   c1->cd(iframe++);
   ySlice->Draw("CONT");
   c1->cd(iframe++);
   zSlice->Draw("CONT");
   c1->cd(iframe++);
   c1->cd(iframe++);
   xSlice->Draw("CONT");

   c1->cd(iframe++);

   c1->Update();

   ps->NewPage();
  
   c1->Clear();
   c1->SetFillColor(0);

   gPad->SetLeftMargin(5.);
   gPad->SetTopMargin(2.);
  
   c1->Divide(2,2);
   Int_t iframe=1;


   c1->cd(iframe++);
   ySlice->Draw("CONT1");
   c1->cd(iframe++);
   zSlice->Draw("CONT1");
   c1->cd(iframe++);
   c1->cd(iframe++);
   xSlice->Draw("CONT1");

   c1->cd(iframe++);

   c1->Update();

   ps->NewPage();
  
   c1->Clear();
   c1->SetFillColor(0);

   gPad->SetLeftMargin(5.);
   gPad->SetTopMargin(2.);
  
   c1->Divide(2,2);
   Int_t iframe=1;


   c1->cd(iframe++);
   ySlice->Draw("SURF3");
   c1->cd(iframe++);
   zSlice->Draw("SURF3");
   c1->cd(iframe++);
   c1->cd(iframe++);
   xSlice->Draw("SURF3");

   c1->cd(iframe++);


   c1->Update();

   ps.Close();
   cout << "Line: " << __LINE__ << endl;

} //drawLaser