#include "TROOT.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TF1.h"
#include "TH1.h"
#include "TAxis.h"
#include "TLatex.h"
#include "TString.h"
#include "TGraph.h"
#include "TGaxis.h"
#include "TLine.h"
#include "TEllipse.h"
#include "TArrow.h"
#include "Math/SpecFunc.h"
#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <string>
#include <iomanip>

Double_t real_sph_hankel_1st(Int_t n, Double_t x) {
  /*same as the output of the AHNKEV subroutine in the Hamamoto's code*/
  return sqrt(2.0/(3.14159265358979323846*x))
    *ROOT::Math::cyl_bessel_k((Double_t)n+0.5,x);
}
Double_t x_sph_bessel(Int_t n, Double_t x) {
  return x*ROOT::Math::sph_bessel(n,x);
}
Double_t x_sph_neumann(Int_t n, Double_t x) {
  return x*ROOT::Math::sph_neumann(n,x);
}
Double_t deriv_sph_bessel(Int_t n, Double_t x) {
  return -ROOT::Math::sph_bessel(n+1,x)
    +(n/x)*ROOT::Math::sph_bessel(n,x);
}
Double_t deriv_sph_neumann(Int_t n, Double_t x) {
  return -ROOT::Math::sph_neumann(n+1,x)
    +(n/x)*ROOT::Math::sph_neumann(n,x);
}
Double_t deriv_x_sph_bessel(Int_t n, Double_t x) {
  return ROOT::Math::sph_bessel(n,x)
    +x*deriv_sph_bessel(n,x);
}
Double_t deriv_x_sph_neumann(Int_t n, Double_t x) {
  return ROOT::Math::sph_neumann(n,x)
    +x*deriv_sph_neumann(n,x);
}
Double_t sph_bessel_neumann(Int_t n, Double_t x, Double_t ps) {
  return cos(ps)*ROOT::Math::sph_bessel(n,x)
    -sin(ps)*ROOT::Math::sph_neumann(n,x);
}
Double_t x_sph_bessel_neumann(Int_t n, Double_t x, Double_t ps) {
  return cos(ps)*x_sph_bessel(n,x)
    -sin(ps)*x_sph_neumann(n,x);
}
Double_t deriv_x_sph_bessel_neumann(Int_t n, Double_t x, Double_t ps) {
  return cos(ps)*deriv_x_sph_bessel(n,x)
    -sin(ps)*deriv_x_sph_neumann(n,x);
}

class WaveFunction {
public:
  WaveFunction(Double_t rmassin, Double_t vin, Double_t ain, Double_t Rin)
    : rmass(rmassin), v(vin), a(ain),  R(Rin), 
      ndiv(300), rmin(-0.1), rmax(59.9), dr((rmax-rmin)/ndiv),
      hbarc(197.32891),/*  hbarc value in Hamamoto-san's code (main-eg.f)*/
      /*hbarc(197.11),  hbarc value in Hamamoto-san's code (cc-nilson-g.f)*/
      hLambda(hbarc/(2.*rmass)), cfunit_inv(2.*rmass/(hbarc*hbarc))
  {
    calc_faclog();
    wf          = new double[ndiv+4];
    r           = new double[ndiv+4];
    Vcenter     = new double[ndiv+1];
    Vcenter_cfb = new double[ndiv+1];
    wfasymp     = new double[ndiv+1];
    wfasymp_sin = new double[ndiv+1];
  }
  
  Double_t* wf_c(Int_t l2in, Int_t j2in,
		 Double_t V0in, Double_t Eljin) {
    l1 = l2in/2;
    l2 = l2in;
    j1 = j2in/2.;
    j2 = j2in;
    ls = (j2*(j2+2) - l2*(l2+2) - 1*(1+2)) / 8.;

    V0 = V0in;
    Elj = Eljin;
    
    for (Int_t i = 0; i <= ndiv; i++) {
      r[i] = dr * (i+0.5)+rmin;
    }
    stormer(1,ndiv+2);
    ps      = FindPS(ndiv);

    Double_t kouter  = sqrt(2.0*rmass*Elj)/hbarc;
    Double_t normfac = sqrt(cfunit_inv/(3.141592653589793238*kouter));
    Double_t uouter  = normfac*x_sph_bessel_neumann(l1, kouter*r[ndiv], ps);
    Double_t uratio  = uouter/wf[ndiv];
    if (uouter*wf[ndiv] < 0) {
      uratio = -uratio;
    }
    for (Int_t i = 1; i <=ndiv; i++) {
      wf[i] = uratio*wf[i];
    }
    for (Int_t i = ndiv+1; i > 1; i--) {
      wf[i] = wf[i-1];
      r[i] = r[i-1];
    }
    wf[1] = 0.;
    r[1] = 0.;
    return wf;
  }
  
  Double_t* wf_b(Int_t nodein, Int_t l2in, Int_t j2in,
		 Int_t opt, Double_t energy)
  {
    node = nodein;
    l1   = l2in/2;
    l2   = l2in;
    j1   = j2in/2.;
    j2   = j2in;
    ls   = (j2*(j2+2) - l2*(l2+2) - 1*(1+2)) / 8.;
    
    for (Int_t i = 0; i <= ndiv+3; i++) {
      r[i] = dr * (i+0.5)+rmin;
    }
    
    Int_t n = (Int_t) (R/0.2);
    Double_t kmin;
    Double_t kmax;
    if (node == 0) {
      kmin = 0.00001;
      kmax = GetZeroPosOfj(l1,0) / R;
    }else{
      kmin = GetZeroPosOfj(l1,node-1) / R;
      kmax = GetZeroPosOfj(l1,node) / R;
    }
    Double_t EljminusV0_min = pow((kmin*hbarc),2)/(2*rmass);
    Double_t EljminusV0_max = pow((kmax*hbarc),2)/(2*rmass);
    if (opt == 0) {
      Elj = energy;
      Double_t V0_min = Elj - EljminusV0_max;
      Double_t V0_max = Elj - EljminusV0_min;
      V0 = FindEne(n, opt, V0_min, V0_max);
    }else{
      V0 = energy;
      Double_t Elj_min = 0.95*(1.05*EljminusV0_min + V0);
      Double_t Elj_max = 0.95*(1.05*EljminusV0_max + V0);
      std::cout << "Elj_min = " << Elj_min << ", Elj_max = " << Elj_max << std::endl;
      if (Elj_min > 0) {
    	exit (0);
    	Elj_min = 0.;
      }
      if (Elj_max > 0) {Elj_max = 0.;}
      Elj = FindEne(n, opt, Elj_min, Elj_max);
    }
    std::cout << "N,l2,j2: " <<  node << "," << l2 << "," << j2 <<  std::endl;
    //std::cout << "V0 = " << V0 << " MeV" << std::endl;
    std::cout << "Elj = " << Elj << " MeV" << std::endl;
    
    stormer(1,n);
    Double_t uinner   = wf[n];
    stormer(ndiv+3,n);
    Double_t uouter   = wf[n];
    Double_t uratio = uinner/uouter;
    
    for (Int_t i = n; i <= ndiv; i++) {
      wf[i] = uratio*wf[i];
    }
    Double_t sum = 0;
    for (Int_t i = ndiv; i >= 1; i--) {
      sum = sum + pow(wf[i],2);
    }
    Double_t norm = 1./sqrt(sum*dr);
    for (Int_t i = 1; i <=ndiv; i++) {
      wf[i] = wf[i] * norm;
    }
    for (Int_t i = ndiv+1; i > 1; i--) {
      wf[i] = wf[i-1];
      r[i] = r[i-1];
    }
    wf[1] = 0.;
    r[1] = 0.;
    return wf;
  }
  ~WaveFunction(){
    delete[] wf;
    delete[] r;
    delete[] Vcenter;
    delete[] Vcenter_cfb;
    delete[] wfasymp_sin;
  }
  
private:
  Double_t FindPS(Int_t n) {
    Double_t a_outer = sqrt(2.0*rmass*Elj)/hbarc;
    Double_t x_outer = a_outer * r[n];
    Double_t slope = (8.0*(wf[n+1]-wf[n-1])-(wf[n+2]-wf[n-2]))/(12.*dr);
    Double_t Llj = r[n]*slope/wf[n];
    Double_t fcosps = Llj * ROOT::Math::sph_neumann(l1,x_outer)
      - deriv_x_sph_neumann(l1,x_outer);
    Double_t fsinps = Llj * ROOT::Math::sph_bessel(l1,x_outer)
      - deriv_x_sph_bessel(l1,x_outer);
    Double_t cosps = fcosps / sqrt(pow(fcosps,2.0)+pow(fsinps,2.0));
    Double_t pstmp = acos(cosps);
    if (fsinps < 0) {pstmp = -pstmp;}
    if (pstmp < 0)  {pstmp = pstmp + 3.141592653589793238;}
    return pstmp;
  }
  
  Double_t FindEne(Int_t n, Int_t opt, Double_t ene_min, Double_t ene_max) {
    if (opt == 0) {V0  = ene_min;}
    if (opt == 1) {Elj = ene_min;}
    
    stormer(ndiv+3,n-2);
    Double_t numerator_outer   = wf[n];
    std::cout << "numerator_outer = " << numerator_outer << std::endl;
    Double_t denominator_outer = (8.0*(wf[n+1]-wf[n-1])-(wf[n+2]-wf[n-2]))/(12.*dr);
    stormer(1,n+2);
    Double_t numerator_inner   = wf[n];
    Double_t denominator_inner = (8.0*(wf[n+1]-wf[n-1])-(wf[n+2]-wf[n-2]))/(12.*dr);
    Double_t delta0 = numerator_outer*denominator_inner
      - numerator_inner*denominator_outer;
    for (Int_t k=0; k < 50; k++){
      if (opt == 0) {V0  = 0.5*(ene_min+ene_max);}
      if (opt == 1) {Elj = 0.5*(ene_min+ene_max);}
      stormer(ndiv+3,n-2);
      numerator_outer   = wf[n];
      denominator_outer = (8.0*(wf[n+1]-wf[n-1])-(wf[n+2]-wf[n-2]))/(12.*dr);
      stormer(1,n+2);
      numerator_inner   = wf[n];
      denominator_inner = (8.0*(wf[n+1]-wf[n-1])-(wf[n+2]-wf[n-2]))/(12.*dr);
      Double_t delta = numerator_outer*denominator_inner
	- numerator_inner*denominator_outer;
      if (delta0*delta > 0) {
	if (opt == 0) {ene_min = V0;}
	if (opt == 1) {ene_min = Elj;}
      }else{
	if (opt == 0) {ene_max = V0;}
	if (opt == 1) {ene_max = Elj;}
      }
    }
    if (opt == 0) {return V0;}
    if (opt == 1) {return Elj;}
    return 0;
  }

  void stormer(Int_t n1, Int_t n2) {
    Double_t pre41=19./240.;
    Double_t pre42=-2./5.;
    Double_t pre43=97./120.;
    Double_t pre44=-11./15.;
    Double_t pre45=299./240.;
    Double_t fprer[6];
    Double_t dxx_var[5] = {0.000, 0.025, 0.050, 0.1, 0.2};
    Double_t xx_var[315] = 
      {0.000,
       0.025,0.050,0.075,0.100,
       0.125,0.150,0.175,0.200,
       0.250,0.300,0.350,0.400,
       0.500,0.600,0.700,0.800,
      };
    for (Int_t i = 17;i<315;i++){
      xx_var[i] = (i-17)*0.2+1.0;
    }
    if (n1 < n2) {
      Int_t    n_start[6] = {0, 5, 9, 13, 17};
      n_start[5] = n2+16;
      Double_t fprer_stored[315];
      Double_t wf_stored[315];
      
      Double_t k_inner = sqrt(2.0*rmass*(Elj-V0))/hbarc;
      wf_stored[0] = 0.;
      //wf_stored[1] = (0.5/k_inner)*pow(2.*dxx_var[1]*k_inner,l1+1)*exp(faclog[l1+1]-faclog[2*(l1+1)]);
      wf_stored[1] = (0.5/k_inner)*pow(2.*dxx_var[1]*k_inner,l1+1)*exp(0.0);
      for (Int_t i=2;i<=5;i++) {
	wf_stored[i] = 2.*wf_stored[i-1]-wf_stored[i-2]+pow(dxx_var[1],2)*pot(xx_var[i-1])*wf_stored[i-1];
      }
      fprer[1]        = 0.0;
      fprer_stored[0] = fprer[1];
      if(l1==1) {
	fprer[1] = 0.6666666666*pow(dxx_var[1],2)*k_inner;
	fprer_stored[0] = fprer[1];
      }
      for (Int_t i=2;i<=5;i++) {
	fprer[i] = wf_stored[n1+i-2]*pow(dxx_var[1],2)*pot(xx_var[n1+i-2]);
	fprer_stored[i-1] = fprer[i];
      }
      for (Int_t inex=1; inex<=4; inex++){
	for (Int_t k=n_start[inex]; k<n_start[inex+1]; k++) {
	  Double_t termr = pre41*fprer[1]+pre42*fprer[2]+pre43*fprer[3]+pre44*fprer[4]+pre45*fprer[5];
	  if (inex!=1 && k==n_start[inex]) {
	    wf_stored[k] = 2.*wf_stored[k-1]-wf_stored[k-3]+termr;
	  }else{
	    wf_stored[k] = 2.*wf_stored[k-1]-wf_stored[k-2]+termr;
	  }
	  for (Int_t i=1;i<=4;i++) {
	    fprer[i]=fprer[i+1];
	  }
	  fprer[5]=wf_stored[k]*pow(dxx_var[inex],2)*pot(xx_var[k]);
	  fprer_stored[k]=fprer[5];
	}
	if (inex == 1) {
	  for (Int_t i=0;i<=4;i++) {
	    fprer[i+1] = fprer_stored[i*2]*4.0;
	  }
	} else if (inex == 2) {
	  fprer[5] = fprer_stored[12]*4.0; /*n=12,x=0.4*/
	  fprer[4] = fprer_stored[10]*4.0; /*n=10,x=0.3*/
	  fprer[3] = fprer_stored[8]*16.0; /*n=8,x=0.2*/
	  fprer[2] = fprer_stored[4]*16.0; /*n=4,x=0.1*/
	  fprer[1] = fprer_stored[0]*16.0; /*n=0,x=0.0*/
	} else if (inex == 3) {
	  fprer[5] = fprer_stored[16]*4.0;  /*n=16,x=0.8*/
	  fprer[4] = fprer_stored[14]*4.0;  /*n=14,x=0.6*/
	  fprer[3] = fprer_stored[12]*16.0; /*n=12,x=0.4*/
	  fprer[2] = fprer_stored[8]*64.0;  /*n=8,x=0.2*/
	  fprer[1] = fprer_stored[0]*64.0;  /*n=0,x=0.0*/
	}
      }	
      wf[1] = wf_stored[8];
      wf[2] = wf_stored[12];
      wf[3] = wf_stored[14];
      wf[4] = wf_stored[16];
      for (Int_t i=5;i<=n2;i++) {
	wf[i] = wf_stored[i+12];
      }
    }else{
      Double_t k_outer = sqrt(2.0*rmass*fabs(Elj))/hbarc;
      for (Int_t i=0;i<5;i++) {
      	wf[n1-i] = r[n1-i]*real_sph_hankel_1st(l1,k_outer*r[n1-i]);
      }
      for (Int_t i=0;i<5;i++) {
      	fprer[i+1] = wf[n1-i]*pow(dr,2)*pot(r[n1-i]);
      }
      for (Int_t k=n1-5;k>=n2;k--) {
	Double_t termr = 
	  pre41*fprer[1]+
	  pre42*fprer[2]+
	  pre43*fprer[3]+
	  pre44*fprer[4]+
	  pre45*fprer[5];
	wf[k] = 2.*wf[k+1]-wf[k+2]+termr;

	for (Int_t i=1;i<=4;i++) {
	  fprer[i]=fprer[i+1];
	}
	fprer[5]=wf[k]*pow(dr,2)*pot(r[k]);
      }
    }
  }

  Double_t pot(Double_t x)
  {
    Double_t cfb = l1*(l1+1)/pow(x,2);
    Double_t tmp = (x-R)/a;
    Double_t ex, f, V, Vso,Vsum;
    ex   = exp(tmp);
    f    = 1./(1.+ex);
    V    = V0 * f;
    Vso  = -V0*v*pow(hLambda,2)/x*(-ex*pow(f,2)/a)*2.*ls;
    //Vso  = -V0*0.44*pow(1.27,2)/x*(-ex*pow(f,2)/a)*ls;
    Vsum = cfb-cfunit_inv*(Elj-V-Vso);
    return Vsum;
  }

  Double_t GetZeroPosOfj(Int_t n, Int_t i) {
    Double_t zeropos[9][6] = {
      {3.14159, 6.28319, 9.42478, 12.5664, 15.7080, 18.8496},
      {4.49341, 7.72525, 10.9041, 14.0662, 17.2208, 20.3713},
      {5.76346, 9.09501, 12.3229, 15.5146, 18.6890, 21.8539},
      {6.98793, 10.4171, 13.6980, 16.9236, 20.1218, 23.3042},
      {8.18256, 11.7049, 15.0397, 18.3013, 21.5254, 24.7276},
      {9.35581, 12.9665, 16.3547, 19.6532, 22.9046, 26.1278},
      {10.5128, 14.2074, 17.6480, 20.9835, 24.2628, 27.5079},
      {11.6570, 15.4313, 18.9230, 22.2953, 25.6029, 28.8704},
      {12.7908, 16.6410, 20.1825, 23.5913, 26.9270, 30.2173}};
    if ((n < 0)||(n > 8)||(i < 0)||(i > 5)) {
      std::cout << "error in GetZerosPosOfj" << std::endl;
      return -1;
    }
    return zeropos[n][i];
  }
  
  void calc_faclog() {
    faclog[1]=(Float_t)0.0;
    faclog[2]=(Float_t)0.0;
    Double_t FN=(Float_t)1.0;
    for (Int_t  N=3; N<=500; N++) {
      FN = FN + (Float_t)1.0;
      faclog[N]=faclog[N-1]+(Float_t)log(FN);
    }
  }
  
public:
  Double_t* Getwf(){return wf;}
  Double_t* GetVcenter(){
    for (Int_t i=1;i<=ndiv;i++) {
      Double_t tmp = (r[i]-R)/a;
      Double_t ex, f, V;
      ex   = exp(tmp);
      f    = 1./(1.+ex);
      V    = V0 * f;
      Vcenter[i] = V;
    }
    return Vcenter;
  }
  Double_t* GetVcenter_cfb(){
    for (Int_t i=1;i<=ndiv;i++) {
      Double_t cfb = l1*(l1+1)/pow(r[i],2);
      Double_t tmp = (r[i]-R)/a;
      Double_t ex, f, V, Vso;
      ex   = exp(tmp);
      f    = 1./(1.+ex);
      V    = V0 * f;
      Vso  = -V0*v*pow(hLambda,2)/r[i]*(-ex*pow(f,2)/a)*2.*ls;
      Vcenter_cfb[i] = V+Vso+cfb/cfunit_inv;
      //std::cout <<  Vcenter_cfb[i] << std::endl;
    }
    return Vcenter_cfb;
  }
  Double_t* GetAsymp(){
    for (Int_t i=2;i<=ndiv;i++) {
      Double_t kouter  = sqrt(2.0*rmass*Elj)/hbarc;
      Double_t normfac = sqrt(cfunit_inv/(3.141592653589793238*kouter));
      Double_t uouter  = normfac*x_sph_bessel_neumann(l1, kouter*r[i], ps);
      wfasymp[i] = uouter;
    }
    return wfasymp;
  }
  Double_t* GetAsympSin(){
    for (Int_t i=1;i<=ndiv;i++) {
      Double_t kouter  = sqrt(2.0*rmass*Elj)/hbarc;
      //Double_t normfac = sqrt(cfunit_inv/(3.141592653589793238*kouter));
      Double_t uouter  = sin(kouter*r[i]+ps-l1*3.1415926535/2.0);
      wfasymp_sin[i] = uouter;
    }
    return wfasymp_sin;
  }

  Double_t GetV(){return V0;}
  Double_t GetE(){return Elj;}
  Double_t GetPS(){return ps;}
  
private:
  Int_t    node, l1, l2;
  Double_t j1;
  Int_t    j2;
  Double_t rmass, v, a, R, V0, Elj;
  Int_t    ndiv;
  Double_t rmin, rmax, dr;
  Double_t hbarc, ls, hLambda, cfunit_inv;
  Double_t ps;
  Double_t faclog[501];
  Double_t *wf;
  Double_t *r;
  Double_t *Vcenter;
  Double_t *Vcenter_cfb;
  Double_t *wfasymp;
  Double_t *wfasymp_sin;
};

std::vector<Double_t> widen(Double_t width, std::vector<Double_t> vecin) {
  std::vector<Double_t> vecout;
  UInt_t i = 0;
  while (i < vecin.size()) {
    UInt_t j;
    for (j = vecin.size() - 1; j > i; j--) {
      if ((vecin[j] - vecin[i]) < width * (j-i)) {
        Double_t center = (vecin[j]+vecin[i])/2.;
        Double_t lowest = center - width*(j-i)/2.;
        for (UInt_t k = i; k <= j; k++){
          vecout.push_back(lowest + width*(k-i));
        }
        break;
      }
    }
    if (i==j){
      vecout.push_back(vecin[i]);
    }
    i = j + 1;
    Int_t ground_state_fixed = 0;
    if ((ground_state_fixed == 1) && (vecout[0] != vecin[0])) {
      vecout[0] = vecin[0];
      for (UInt_t k = 1; k < vecin.size(); k++) {
        if ((vecout[k]-vecout[k-1]) < width) {
          vecout[k] = vecout[k-1] + width;
        }
      }
    }
  }
  return vecout;
}

void DrawCircles(Double_t xpos, std::vector<Double_t> ypos, Double_t r1, Double_t r2, Double_t distance,
		 std::vector<Int_t> nClosedCircles, std::vector<Int_t> nOpenCircles) {
  TEllipse el1,el2;
  el1.SetFillColor(kAzure-2);
  //el2.SetFillStyle(4000);
  el1.SetLineColor(1);
  el2.SetFillColor(0);
  el2.SetLineStyle(1);
  
  for (UInt_t i=0; i < ypos.size(); i++){
    Double_t xleft  = xpos - distance*(nClosedCircles[i]+nOpenCircles[i]-1)/2.0;
    Double_t xright = xpos + distance*(nClosedCircles[i]+nOpenCircles[i]-1)/2.0;
    for (Int_t j = 0; j < nClosedCircles[i];j++) {
      el1.DrawEllipse(xleft+distance*j,ypos[i],r1,r2,0.,360.,0);
    }
    for (Int_t j = 0; j < nOpenCircles[i];j++) {
      el2.DrawEllipse(xright-distance*j,ypos[i],r1,r2,0.,360.,0);
    }
  }   
}

void init() {
  gROOT->SetStyle("Plain");
  gStyle->SetPadLeftMargin(0.15);
  //  gStyle->SetAxisColor(1,"X");
  gStyle->SetTickLength(-0.01,"Y");
  gStyle->SetTickLength(-0.01,"X");
  gStyle->SetLabelOffset(0.03,"Y");
  gStyle->SetLabelOffset(0.01,"X");
  gStyle->SetPadTickY(1);
  gStyle->SetPadTickX(1);
  gStyle->SetTitleYOffset(1.6);
  gStyle->SetTitleXOffset(1.2);
  gStyle->SetLabelFont(132,"xyz");
  gStyle->SetTitleFont(132,"xyz");
  gStyle->SetTitleFont(132,"");
  gStyle->SetTextFont(132);
  gStyle->SetLabelSize(0.04,"xyz");
  gStyle->SetTitleSize(0.04,"xyz");
  gStyle->SetTextSize(0.04);
  gStyle->SetPadRightMargin(0.05);
  gStyle->SetPadLeftMargin(0.15);
  gStyle->SetPadBottomMargin(0.12);
  gStyle->SetPadTopMargin(0.03);
}

void plot_ws_res() {
  init();
  TCanvas *c1 = new TCanvas("c1","c1",500,500);
  TH1D* vws         = new TH1D("vws","#it{V}_{WS}",300,-0.1,59.9);
  TH1D* vwscfb      = new TH1D("vwscfb","#it{V}_{WS}",300,-0.1,59.9);
  TH1D* hcon        = new TH1D("hcon","hcon",300,-0.1,59.9);
  WaveFunction* wfpot = new WaveFunction((30./31.)*939.577218339,32.,0.67,1.27*pow(30,1./3.));

  wfpot->wf_c(6,7,-39.,1.);
  vws->SetContent(wfpot->GetVcenter());
  vwscfb->SetContent(wfpot->GetVcenter_cfb());
  vws->SetAxisRange(-0.1,9.,"X");
  vwscfb->SetAxisRange(0.1,40.,"X");
  vwscfb->SetAxisRange(-40,40.,"Y");
  hcon->SetAxisRange(-0.1,60.,"X");
  
  std::vector<std::string>  orb_str;
  orb_str.push_back("1#kern[0.1]{#it{s}_{1/2}}");
  orb_str.push_back("1#kern[0.1]{#it{p}_{3/2}}");
  orb_str.push_back("1#kern[0.1]{#it{p}_{1/2}}");
  orb_str.push_back("1#kern[0.1]{#it{d}_{5/2}}");
  orb_str.push_back("2#kern[0.1]{#it{s}_{1/2}}");
  orb_str.push_back("1#kern[0.1]{#it{d}_{3/2}}");
  
  const Int_t n_orb = 6;
  TH1D* h[n_orb];
  for (Int_t i =0;i<n_orb;i++) {
    h[i] = new TH1D(Form("h%d",i),orb_str[i].c_str(),300,-0.1,59.9);
  }

  std::vector<Double_t>     ene1,ene2,enew;
  h[0]->SetContent(wfpot->wf_b(0,0,1,1,-39.)); ene1.push_back(wfpot->GetE());
  h[1]->SetContent(wfpot->wf_b(0,2,3,1,-39.)); ene1.push_back(wfpot->GetE());
  h[2]->SetContent(wfpot->wf_b(0,2,1,1,-39.)); ene1.push_back(wfpot->GetE());
  h[3]->SetContent(wfpot->wf_b(0,4,5,1,-39.)); ene1.push_back(wfpot->GetE());
  h[4]->SetContent(wfpot->wf_b(1,0,1,1,-39.)); ene1.push_back(wfpot->GetE());
  h[5]->SetContent(wfpot->wf_b(0,4,3,1,-39.)); ene1.push_back(wfpot->GetE());
  enew = widen(3,ene1);
  
  TH1D* h_offset[n_orb];
  for (Int_t i =0;i<n_orb;i++) {
    h_offset[i]=(TH1D*)h[i]->Clone();
    for (Int_t j=1;j<=h[i]->GetNbinsX();j++) {
      h_offset[i]->SetBinContent(j,h[i]->GetBinContent(j)*8+ene1[i]);
    }
    h_offset[i]->SetAxisRange(-0.1,15.,"X");
    ene2.push_back(h_offset[i]->GetBinContent(h_offset[i]->FindBin(15.)));
  }

  Double_t ene=0.0;
  Int_t j=0;
  //while (ene<0.2) {
  while (ene<20) {
    j=j+1;
    if (ene < 1.8) {
      ene+=0.05;
    } else if (ene < 2.89){
      ene+=0.01;
    } else {
      ene+=0.2;
    }
    TH1* frame = c1->DrawFrame(-0.001,-40,70,30);
    frame->GetYaxis()->SetTitle("Eigen energy #it{#varepsilon} (MeV)");
    frame->GetYaxis()->CenterTitle();
    frame->GetXaxis()->SetTitle("Radius #it{r} (fm)");
    frame->GetXaxis()->CenterTitle();
    vws->SetLineWidth(3);
    vws->Draw("same,l");
    vwscfb->SetLineStyle(2);
    vwscfb->Draw("same,l");

    TLine line;
    TLatex text;
    text.SetTextAlign(12);
    for (Int_t i =0;i<n_orb;i++) {
      h_offset[i]->Draw("same,l");
      line.DrawLine(15.,ene2[i],18,enew[i]);
      text.SetTextAlign(12);
      text.DrawLatex(20,enew[i]-0.3,orb_str[i].c_str());
    }

    Double_t* wf_c     = wfpot->wf_c(6,7,-39.,ene);
    for (Int_t i=1;i<=hcon->GetNbinsX();i++) {
      hcon->SetBinContent(i,wf_c[i]*8+ene);
    }
    hcon->Draw("same,l");
    Double_t hcon_end = hcon->GetBinContent(hcon->FindBin(59.8));
    line.DrawLine(60.,hcon_end,62,ene);
    text.DrawLatex(62.5,ene-0.3,"1#kern[0.1]{#it{f}_{#kern[-0.3]{7/2}}}");
    text.DrawLatex(45,25,Form(" #it{#varepsilon}_{#it{c}} = %4.2f MeV",ene));
    
    TLatex notes;
    notes.SetTextSize(0.04);
    notes.DrawLatex(30,-30+15,"{}^{31}_{10}Ne_{21}");
    notes.DrawLatex(30,-30+10.,"#it{a} = 0.67 fm");
    notes.DrawLatex(30,-30+5.,"#it{R} = 3.964 fm (for #it{A} = 30)");
    notes.DrawLatex(30,-30+0.,"#it{V}_{WS} = - 39 MeV");
    notes.SetTextColor(2);
    notes.SetTextAlign(22);
    
    c1->Update();
    
    std::ostringstream os;
    os << "ws_res_ne31_" << std::setw(4) << std::setfill('0') << j << "_" << ene <<".gif";
    c1->Print(os.str().c_str());
  }
}