#include "TROOT.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TH1.h"
#include "TH2.h"
#include <iostream>

Double_t simps(Double_t *f, Double_t a, Double_t b, Int_t n2) {
  /*Based on cernlib simps(D101) subroutine (simps.F & simpscod.inc)*/
  Double_t h=0;
  if((n2 <= 0)||(2*(n2/2)!=n2)) {
    h=0;
    std::cout << "NON-POSITIVE OR EVEN NUMBER OF FUNCTION VALUES =" << n2 <<std::endl;
  }else{
    Double_t s1 = f[n2-1];
    Double_t s2 = 0.;
    for (Int_t n = 1;n<=n2-3;n+=2){
      s1=s1+f[n];
      s2=s2+f[n+1];
    }
    s1=s1+s1+s2;
    h=(f[0]+f[n2]+s1+s1)*(b-a)/(3*n2);
  }
  return h;
}

Double_t V_real(Double_t *z, Double_t *b) {
  Double_t V0 = -40.;
  Double_t R  = 5.0;
  Double_t r  = sqrt(pow(z[0],2)+pow(b[0],2));
  Double_t a  = 0.65;
  return V0/(1+exp((r-R)/a));
}

Double_t V_comp(Double_t *z, Double_t *b) {
  Double_t W0 = -10.;
  Double_t R  = 5.;
  Double_t r  = sqrt(pow(z[0],2)+pow(b[0],2));
  Double_t a  = 0.65;
  return W0/(1+exp((r-R)/a));
}

Double_t ps(Int_t i, TH2D *pot) {
  Double_t hbar   = 197.;
  Double_t m      = 931.;
  Double_t E      = 100.;
  Double_t k      = sqrt(2.*m*E)/hbar;
  Double_t v      = hbar*k/m;
  Int_t    xnbins = pot->GetNbinsX();
  Double_t *arr   = new Double_t[xnbins];
  for (Int_t j=1; j<=xnbins; j++){
    arr[j-1] = pot->GetBinContent(i,j);
  }
  return -(1/(hbar*v))*simps(arr,
			     pot->GetXaxis()->GetBinCenter(1),
			     pot->GetXaxis()->GetBinCenter(xnbins),
			     xnbins-1);
  delete [] arr;
}

Double_t sa_real(Double_t theta, TH1D *sa_real_h, TH1D *sa_comp_h) {
  Double_t hbar   = 197.;
  Double_t m      = 931.;
  Double_t E      = 100.;
  Double_t k      = sqrt(2.*m*E)/hbar;
  Int_t    xnbins = sa_real_h->GetNbinsX();
  Double_t *arr   = new Double_t[xnbins];
  for (Int_t j=1; j<=xnbins; j++){
    Double_t b = sa_real_h->GetXaxis()->GetBinCenter(j);
    arr[j-1] = k*b*TMath::BesselJ0(2*k*sin(theta/2.)*b)*
      sin(sa_real_h->GetBinContent(j))/exp(sa_comp_h->GetBinContent(j));
  }
  return simps(arr,
	       sa_real_h->GetXaxis()->GetBinCenter(1),
	       sa_real_h->GetXaxis()->GetBinCenter(xnbins),
	       xnbins-1);
  delete [] arr;
}

Double_t sa_comp(Double_t theta, TH1D *sa_real_h, TH1D *sa_comp_h) {
  Double_t hbar   = 197.;
  Double_t m      = 931.;
  Double_t E      = 100.;
  Double_t k      = sqrt(2.*m*E)/hbar;
  Int_t    xnbins = sa_real_h->GetNbinsX();
  Double_t *arr   = new Double_t[xnbins];
  for (Int_t j=1; j<=xnbins; j++){
    Double_t b = sa_real_h->GetXaxis()->GetBinCenter(j);
    arr[j-1] = k*b*TMath::BesselJ0(2*k*sin(theta/2.)*b)*
      (1-cos(sa_real_h->GetBinContent(j))/exp(sa_comp_h->GetBinContent(j)));
  }
  return simps(arr,
	       sa_real_h->GetXaxis()->GetBinCenter(1),
	       sa_real_h->GetXaxis()->GetBinCenter(xnbins),
	       xnbins-1);
  delete [] arr;
}

void init() {
  gROOT->SetStyle("Plain");
  gStyle->SetTickLength(0.025,"Y");
  gStyle->SetTickLength(0.015,"X");
  gStyle->SetLabelOffset(0.02,"Y");
  gStyle->SetLabelOffset(0.01,"X");
  gStyle->SetPadTickY(1);
  gStyle->SetPadTickX(1);
  gStyle->SetTitleYOffset(1.4);
  gStyle->SetTitleXOffset(1.0);
  gStyle->SetLabelFont(132,"xyz");
  gStyle->SetTitleFont(132,"xyz");
  gStyle->SetTitleFont(132,"");
  gStyle->SetTextFont(132);
  gStyle->SetLabelSize(0.06,"xyz");
  gStyle->SetTitleSize(0.06,"xyz");
  gStyle->SetTextSize(0.05);
  gStyle->SetPadRightMargin(0.05);
  gStyle->SetPadLeftMargin(0.2);
  gStyle->SetPadBottomMargin(0.15);
  gStyle->SetPadTopMargin(0.05);
  gStyle->SetOptStat(0);
}

void eikonal_neutron() {
  Double_t x1  = -20.;
  Double_t y1  = -20.;
  Double_t x2  = -x1;
  Double_t y2  = -y1;
  Int_t    n   = 400;
  Double_t dx  = (x2-x1)/(2*n);
  Double_t dy  = (y2-y1)/(2*n);
  Double_t dxh = dx/2.;
  Double_t dyh = dy/2.;
  Double_t x1h = x1-dxh;
  Double_t x2h = x2+dxh;
  Double_t y1h = y1-dyh;
  Double_t y2h = y2+dyh;
  Double_t dtheta   = 180./n;
  Double_t dthetah  = dtheta/2.;

  TH2D *v_real_h  = new TH2D("v_real_h", "v_real_h",   2*n+1, x1h, x2h, 2*n+1, y1h, y2h);
  TH2D *v_comp_h  = new TH2D("v_comp_h", "v_comp_h",   2*n+1, x1h, x2h, 2*n+1, y1h, y2h);
  TH1D *ps_real_h = new TH1D("ps_real_h", "ps_real_h", n+1, -dyh, y2h);
  TH1D *ps_comp_h = new TH1D("ps_comp_h", "ps_comp_h", n+1, -dyh, y2h);
  TH1D *sa_real_h = new TH1D("sa_real_h", "sa_real_h", n+1, -dthetah, 180.+dthetah);
  TH1D *sa_comp_h = new TH1D("sa_comp_h", "sa_comp_h", n+1, -dthetah, 180.+dthetah);
  TH1D *sa2_h     = new TH1D("sa2_h",     "",       n+1, -dthetah, 180.+dthetah);

  for (Int_t i=1; i<=2*n+1; i++){
    Double_t x = x1h + dx*(i-1);
    for (Int_t j=1; j<=2*n+1; j++){
      Double_t y = y1h + dy*(j-1);
      v_real_h->SetBinContent(i,j, V_real(&x,&y));
      v_comp_h->SetBinContent(i,j, V_comp(&x,&y));
    }
  }
  for (Int_t j=1; j<=n+1; j++){
    ps_real_h->SetBinContent(j, ps(j+n,v_real_h));
    ps_comp_h->SetBinContent(j, ps(j+n,v_comp_h));
  }

  for (Int_t j=1; j<=n+1; j++){
    Double_t theta = dtheta*(j-1)*3.141592/180.;
    Double_t integ_real = sa_real(theta,ps_real_h,ps_comp_h);
    Double_t integ_comp = sa_comp(theta,ps_real_h,ps_comp_h);
    sa_real_h->SetBinContent(j, integ_real);
    sa_comp_h->SetBinContent(j, integ_comp);
    sa2_h->SetBinContent(j, (pow(integ_real,2)+pow(integ_comp,2))/100.);
  }

  init();
  TCanvas *c1 = new TCanvas("c1","c1",400,400);
  c1->SetLogy(1);
  sa2_h->GetXaxis()->SetRangeUser(-1.,72.);
  sa2_h->SetLineWidth(2);
  sa2_h->SetLineColor(1);
  TH1* frame = c1->DrawFrame(-dthetah,0.0008, 80.+dthetah,100);
  frame->GetXaxis()->SetTitle("Scattering Angle (degree)"); 
  frame->GetYaxis()->SetTitle("Differential cross section #it{d}#sigma/#it{d}#Omega (b/sr)");
  frame->GetYaxis()->CenterTitle();
  frame->GetXaxis()->CenterTitle();
  sa2_h->Draw("lsame");
}