//2015 1/9
// A. Miyazaki
//

#include <unistd.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <libgen.h>
#include <time.h>
#include <math.h>
#include <iostream>
#include <fstream>
#include "TCanvas.h"
#include "TFrame.h"
#include "TPad.h"
#include "TApplication.h"
#include "TDirectory.h"
#include "TFile.h"
#include "TTree.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "TF2.h"
#include "TLine.h"
#include "fstream"
#include "TCut.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TString.h"
#include "TRandom3.h"
#include "TMath.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TSystem.h"
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>

using namespace std;

//pulse function of Pf
const double t1 = 0.; //[sec]
const double t2 = 15.; //[sec]
const double tmax = t2 * 2.; //[sec]

//constants
const double U_over_Eacc2 = 0.207; //J / (MV/m)^2
const double omega = 2*TMath::Pi() * 101.28e6;
const double t_pulse_res = 0.0000000001;

//global constants set by argv
double Q0_0 = 1.e9;
double Qr = 1.e9;
double Qt = 3.5e10;
double Eaccth = 5.;//
double Pf = 2.; //[W]
//omake
double Uth = U_over_Eacc2 * pow(Eaccth, 2);

//------------CAUTION-----------------------------//
//   pulse function causes nan derivative
//   --> smooth function?
//------------------------------------------------//

double func_Pf(double t) {
/*
  if(t>t1 && t<t2) {
    return Pf;
  } else {
    return 0.;
  }
*/

 //Woods-Saxon potential
 return Pf / (TMath::Exp((t-t2)/t_pulse_res) + 1.);
}

//non-linear Q0-function
/*
double func_Q0(double U) {
  return Q0_0;
}
*/

double func_Q0(double U) {
  return Q0_0 * exp(-TMath::Sqrt(U/Uth));
}

/*
double func_Q0(double U) {
  return Q0_0 / (U/Uth + 1.0); // Analytically solvable --> compare numerical trace with truth
}
*/

/*
double func_Q0(double U) {
  return Q0_0 * (1. - (U/Uth));
}
*/

int DiffFunc(double t, const double y[], double f[], void *params) {
  double *_params = (double *)params;
  double _Pf = func_Pf(t);

  double U = y[0];
  double dU_dt = TMath::Sqrt((4.*_Pf*omega*U)/Qr) - (1./func_Q0(U) + 1./Qr + 1./Qt) * omega * U;

  f[0] = dU_dt;

  return GSL_SUCCESS;
}

int main(int argc, char **argv) {

  TApplication theApp("App", &argc, argv);
  
   if(theApp.Argc()<6) {
     cerr << theApp.Argv(0) << " Q0(0) Qr Qt Eaccth Pf" << endl;
     return 1;
  }
  Q0_0 = atof(theApp.Argv(1)); //1.e9;
  Qr = atof(theApp.Argv(2)); //1.e9;
  Qt = atof(theApp.Argv(3)); //3.5e10;
  Eaccth = atof(theApp.Argv(4)); //5.;// parameter of nonlinearity. Given by experiments
  Pf = atof(theApp.Argv(5)); //2.; // incident power
  Uth = U_over_Eacc2 * pow(Eaccth, 2);

  const char *filename = TString::Format("data/simrf_Q00%2.1fe9_Qr%2.1fe9_Qt%2.1fe10_Eaccth%2.1f_Pf%2.1f.dat", Q0_0/1e9, Qr/1e9, Qt/1e10, Eaccth, Pf).Data();

  gROOT->Macro("~/rootlogon.C");
  gStyle->SetPalette(1, 0);

  //double _params[1] = {Pf};
  void * params;// = _params;

  const gsl_odeiv_step_type * Ty = gsl_odeiv_step_rkf45;
  //const gsl_odeiv_step_type * Ty = gsl_odeiv_step_rk4;
  gsl_odeiv_step * St = gsl_odeiv_step_alloc (Ty, 1);
  gsl_odeiv_control *Co = gsl_odeiv_control_y_new (1.e-20, 0.0);
  gsl_odeiv_evolve * Ev = gsl_odeiv_evolve_alloc(1);
  gsl_odeiv_system sys = {DiffFunc, NULL, 1, params};

  TGraph * g_U = new TGraph();
  TGraph * g_Eacc = new TGraph();
  TGraph * g_Pinc = new TGraph();
  TGraph * g_Pref = new TGraph();
  TGraph * g_Ptrn = new TGraph();

  TGraph *g_Q_U0 = new TGraph();
  g_Q_U0->SetMarkerColor(2);
  g_Q_U0->SetMarkerStyle(29);
  g_Q_U0->SetMarkerSize(5);
  

  //min and max
  double Pinc_min = 1000;
  double Pinc_max = -1000;
  double Pref_min = 1000;
  double Pref_max = -1000;
  double Ptrn_min = 1000;
  double Ptrn_max = -1000;

  //the variables calculated
  double t = 0;
  double dt = 1.e-9; //

  //double Eacc0 = 0.2;//
  //double U0 = U_over_Eacc2 * pow(Eacc0, 2);
  double U0 = 0.00001; // should be more than "0" otherwize RKF never evolves
  double y[1] =  {U0};

  int num=0;

  ofstream oData(filename);
  oData << "time [sec]\t" << "U [J]\t" << "Eacc [MV/m]\t" << "Q0\t" << "Qr\t" << "Qt\t" << "Q0\t" << "Pinc [W]\t" << "Pref [W]\t" << "Ptrn [W]\t" << endl;
  while (t < tmax) {
    int status = gsl_odeiv_evolve_apply(Ev, Co, St, &sys, &t, tmax, &dt, y);
 
    if(status != GSL_SUCCESS) {
       cout << status << endl;
       break;
    }

    double U = y[0];
    //if (U<0.) U=0.;
    //if(dt > 1.e-3) dt=1.3-6; //

    double Q0 = func_Q0(U);
    double Eacc = TMath::Sqrt(U/U_over_Eacc2); 
    double Pinc = func_Pf(t);
    double Ptrn = omega * U / Qt;
    double Pr_fromCav = omega * U / Qr;
    double Pref = TMath::Power( TMath::Sqrt(Pr_fromCav) - TMath::Sqrt(Pinc), 2);

    if(Pinc < Pinc_min) Pinc_min = Pinc;
    if(Pref < Pref_min) Pref_min = Pref;
    if(Ptrn < Ptrn_min) Ptrn_min = Ptrn;

    if(Pinc > Pinc_max) Pinc_max = Pinc;
    if(Pref > Pref_max) Pref_max = Pref;
    if(Ptrn > Ptrn_max) Ptrn_max = Ptrn;

    g_U->SetPoint(num, t, U);
    g_Eacc->SetPoint(num, t, Eacc);
    g_Pinc->SetPoint(num, t, sqrt(Pinc)); 
    g_Ptrn->SetPoint(num, t, sqrt(Ptrn)); 
    g_Pref->SetPoint(num, t, sqrt(Pref)); 
    if(t < t2 && t>t2-dt) {
      g_Q_U0->SetPoint(0, Eacc, Q0);
    }
    oData << t <<"\t" << U << "\t" << Eacc << "\t" << Q0 << "\t" << Qr << "\t" << Qt << "\t" << Pinc <<"\t" << Pref << "\t" << Ptrn <<"\t" << endl;

    if(num>10000000000) {
       cout << "ahoaho" << endl;
       break;
    }
    num++;
  }
  oData.close();
  gsl_odeiv_evolve_free (Ev);
  gsl_odeiv_control_free (Co);
  gsl_odeiv_step_free (St);

  Pinc_max = sqrt(Pinc_max);
  Pinc_min = sqrt(Pinc_min);
  Pref_max = sqrt(Pref_max);
  Pref_min = sqrt(Pref_min);
  Ptrn_max = sqrt(Ptrn_max);
  Ptrn_min = sqrt(Ptrn_min);

  //plot Q(Eacc)
  TCanvas *c1 = new TCanvas("c1", "c1", 600, 400);
  gPad->SetGridx();
  gPad->SetGridy();
  TGraph * g_Q = new TGraph();
  g_Q->SetMarkerStyle(20);
  for(int j=0; j<1000; j++) {
    double Eacc = 7./1000. * (double)j;
    double U = U_over_Eacc2 * pow(Eacc,2);
    double val = func_Q0(U);
    g_Q->SetPoint(j, Eacc, val); 
  }
  g_Q->GetXaxis()->SetTitle("Eacc [MV/m]");
  g_Q->GetXaxis()->SetTitleSize(0.06);
  g_Q->GetXaxis()->SetLabelSize(0.06);
  g_Q->GetYaxis()->SetTitle("Q_{0}");
  g_Q->GetYaxis()->SetTitleSize(0.06);
  g_Q->GetYaxis()->SetTitleOffset(1.2);
  g_Q->GetYaxis()->SetLabelSize(0.06);
  g_Q->Draw("AP");
  g_Q_U0->Draw("P");
  gPad->SetLogy();

  //plot observed powers
  TCanvas *c2 = new TCanvas("c2", "c2", 600, 900);
  TPad *pad1 = new TPad("pad1","This is pad1",0.02,0.67,0.98,0.98,0);
  TPad *pad2 = new TPad("pad2","This is pad2",0.02,0.34,0.98,0.66,0);
  TPad *pad3 = new TPad("pad3","This is pad3",0.02,0.02,0.98,0.33,0);
  pad1->Draw();
  pad2->Draw();
  pad3->Draw();

  pad1->cd();
  pad1->SetGridx();
  pad1->SetGridy();
  pad1->GetFrame()->SetFillColor(15);
  pad1->SetFillColor(12);

  TH1F *frame1 = pad1->DrawFrame(t1, (Pinc_min - (Pinc_max-Pinc_min)*0.1) , tmax, (Pinc_max + (Pinc_max-Pinc_min)*0.1));
  pad1->GetFrame()->SetFillColor(12);
  pad1->GetFrame()->SetFillStyle(3000);
  frame1->SetMarkerStyle(20);
  frame1->GetXaxis()->SetAxisColor(3);
  frame1->GetXaxis()->SetTitle("t [sec]");
  frame1->GetXaxis()->SetTitleSize(0.06);
  frame1->GetXaxis()->SetTitleColor(3);
  frame1->GetXaxis()->SetLabelSize(0.06);
  frame1->GetXaxis()->SetLabelColor(3);
  frame1->GetYaxis()->SetTitleColor(3);
  frame1->GetYaxis()->SetAxisColor(3);
  //frame1->GetYaxis()->SetTitle("Pinc [W]");
  frame1->GetYaxis()->SetTitle("Vi [V]");
  frame1->GetYaxis()->SetTitleSize(0.06);
  frame1->GetYaxis()->SetTitleOffset(1.2);
  frame1->GetYaxis()->SetTitleColor(3);
  frame1->GetYaxis()->SetLabelSize(0.06);
  frame1->GetYaxis()->SetLabelColor(3);
  g_Pinc->SetLineWidth(3);
  g_Pinc->SetLineColor(3);
  g_Pinc->Draw("L");

  pad2->cd();
  pad2->SetGridx();
  pad2->SetGridy();
  pad2->GetFrame()->SetFillColor(15);
  pad2->SetFillColor(12);

  TH1F *frame2 = pad2->DrawFrame(t1, (Pref_min - (Pref_max-Pref_min)*0.1) , tmax, (Pref_max + (Pref_max-Pref_min)*0.1));
  pad2->GetFrame()->SetFillColor(12);
  pad2->GetFrame()->SetFillStyle(3000);
  frame2->SetMarkerStyle(20);
  frame2->GetXaxis()->SetAxisColor(3);
  frame2->GetXaxis()->SetTitle("t [sec]");
  frame2->GetXaxis()->SetTitleSize(0.06);
  frame2->GetXaxis()->SetTitleColor(3);
  frame2->GetXaxis()->SetLabelSize(0.06);
  frame2->GetXaxis()->SetLabelColor(3);
  frame2->GetYaxis()->SetAxisColor(3);
  //frame2->GetYaxis()->SetTitle("Pref [W]");
  frame2->GetYaxis()->SetTitle("Vr [V]");
  frame2->GetYaxis()->SetTitleSize(0.06);
  frame2->GetYaxis()->SetTitleOffset(1.2);
  frame2->GetYaxis()->SetTitleColor(3);
  frame2->GetYaxis()->SetLabelSize(0.06);
  frame2->GetYaxis()->SetLabelColor(3);
  g_Pref->SetLineWidth(3);
  g_Pref->SetLineColor(3);
  g_Pref->Draw("L");

  pad3->cd();
  pad3->SetGridx();
  pad3->SetGridy();
  pad3->GetFrame()->SetFillColor(15);
  pad3->SetFillColor(12);

  TH1F *frame3 = pad3->DrawFrame(t1, (Ptrn_min - (Ptrn_max-Ptrn_min)*0.1) , tmax, (Ptrn_max + (Ptrn_max-Ptrn_min)*0.1));
  pad3->GetFrame()->SetFillColor(12);
  pad3->GetFrame()->SetFillStyle(3000);
  frame3->SetMarkerStyle(20);
  frame3->GetXaxis()->SetAxisColor(3);
  frame3->GetXaxis()->SetTitle("t [sec]");
  frame3->GetXaxis()->SetTitleSize(0.06);
  frame3->GetXaxis()->SetTitleColor(3);
  frame3->GetXaxis()->SetLabelSize(0.06);
  frame3->GetXaxis()->SetLabelColor(3);
  frame3->GetYaxis()->SetAxisColor(3);
  //frame3->GetYaxis()->SetTitle("Ptrn [W]");
  frame3->GetYaxis()->SetTitle("Vt [V]");
  frame3->GetYaxis()->SetTitleSize(0.06);
  frame3->GetYaxis()->SetTitleOffset(1.2);
  frame3->GetYaxis()->SetTitleColor(3);
  frame3->GetYaxis()->SetLabelSize(0.06);
  frame3->GetYaxis()->SetLabelColor(3);
  g_Ptrn->SetLineWidth(3);
  g_Ptrn->SetLineColor(3);
  g_Ptrn->Draw("L");

  //////////////////////////////////
  // Get tau_L
  //////////////////////////////////

  // LabView method for tau_L
  TCanvas *c_lab = new TCanvas("c_lab", "c_lab", 800, 600);
  TGraph *g_lab = (TGraph*)g_Ptrn->Clone();
  g_lab->Draw("AL");
  g_lab->SetLineColor(1);
  double markerX1, markerX2;
  double markerY1, markerY2;

  for(int j=0; j<g_lab->GetN(); j++) {
    double x, y;
    g_lab->GetPoint(j, x, y);
    if(fabs(x-t2)<dt) {
      markerX1 = x;
      markerY1 = y;
    }
    if(x>t2 && fabs(y-markerY1/2.)<0.00001) {
      markerX2 = x;
      markerY2 = y;
    }
  }
  TGraph *g_mX1 = new TGraph();
  g_mX1->SetLineColor(2);
  g_mX1->SetLineWidth(2);
  g_mX1->SetPoint(0, markerX1, 0);
  g_mX1->SetPoint(1, markerX1, markerY1);
  g_mX1->Draw("L");

  TGraph *g_mY1 = new TGraph();
  g_mY1->SetLineColor(2);
  g_mY1->SetLineWidth(2);
  g_mY1->SetPoint(0, t1, markerY1);
  g_mY1->SetPoint(1, markerX1, markerY1);
  g_mY1->Draw("L");

  TGraph *g_mX2 = new TGraph();
  g_mX2->SetLineColor(2);
  g_mX2->SetLineWidth(2);
  g_mX2->SetPoint(0, markerX2, 0);
  g_mX2->SetPoint(1, markerX2, markerY2);
  g_mX2->Draw("L");

  TGraph *g_mY2 = new TGraph();
  g_mY2->SetLineColor(2);
  g_mY2->SetLineWidth(2);
  g_mY2->SetPoint(0, t1, markerY2);
  g_mY2->SetPoint(1, markerX2, markerY2);
  g_mY2->Draw("L");

  double tau_half_lab = (markerX2-markerX1);
  double tau_L_lab = tau_half_lab / (2.*log(2.)) * 1000.;

  TF1 *f_tau_L_lab = new TF1("f_tau_L_lab", TString::Format("%f", tau_L_lab), 0, 1000);
  f_tau_L_lab->SetLineColor(2);

  cout << "tau(1/2) in LabView = " << tau_half_lab << " [s]" << endl;
  cout << "tau_L in LabView = " << tau_L_lab << " [ms]" << endl;


/*
  // Fit
  TCanvas *c_Fit = new TCanvas("c_Fit", "c_Fit", 800, 600);
  TGraphErrors *g_Fit = new TGraphErrors();
  for(int j=0; j<g_lab->GetN(); j++) {
    double x, y;
    g_lab->GetPoint(j, x, y);
    g_Fit->SetPoint(j, x, y);
    g_Fit->SetPointError(j, 0, y*0.01);
  }
  g_Fit->Draw("AP");
  g_Fit->GetXaxis()->SetRangeUser(t2-2, 2.*t2);
  TF1 *f_Fit = new TF1("f_Fit", "[0]*exp(-x/[1])");
  f_Fit->SetParameters(1, 1);
  f_Fit->SetLineColor(2);
  g_Fit->Fit(f_Fit, "MQ", "same", t2+0.01, t2+2.5);
  f_Fit->Draw("same");
  cout << "tau_L in Fit = " << f_Fit->GetParameter(1)*1000./2. << " [ms]" << endl;
*/

  // LLRF
  TCanvas *c_LLRF = new TCanvas("c_LLRF", "c_LLRF", 800, 600);
  TGraph *g_LLRF = (TGraph*)g_lab->Clone();
  for(int j=0; j<g_LLRF->GetN(); j++) {
    double x, y;
    g_LLRF->GetPoint(j, x, y);
    if(x>t2) g_LLRF->SetPoint(j, x, log(y));
  }
  g_LLRF->Draw("AL");
  g_LLRF->GetXaxis()->SetRangeUser(t2, 2.*t2);
  const int n_fit_end = 500;
  double fit_all_end = t2*2-0.001;
  double fit_start = t2 + 0.001;
  TGraph * g_tau_L_LLRF = new TGraph();
  TF1 *f_LLRF = new TF1("f_LLRF", "[0] - x/[1]");
  for(int j=0; j<n_fit_end; j++) {
    double fit_end = fit_start + (fit_all_end - fit_start) * ((double)j)/((double)n_fit_end);
    f_LLRF->SetParameter(1, 1);
    f_LLRF->SetLineColor(2);
    g_LLRF->Fit(f_LLRF, "MQ", "", fit_start, fit_end);
    //cout << "tau_L in LLRF = " << f_LLRF->GetParameter(1)*1000./2. << " [ms]" << endl;
    double tau_L_LLRF = f_LLRF->GetParameter(1)*1000./2.;
    g_tau_L_LLRF->SetPoint(j, fit_end, tau_L_LLRF);
  }

  //truth
  double Q0_truth, aho;
  g_Q_U0->GetPoint(0, aho, Q0_truth);
  double tau_L_truth = 1./(omega) / (1./Q0_truth + 1./Qr + 1./Qt) * 1000.;
  TF1 *f_tau_L_truth = new TF1("f_tau_L_truth", TString::Format("%f", tau_L_truth), 0, 1000);
  f_tau_L_truth->SetLineColor(4);
  cout << "true tau_L = " << tau_L_truth << " [ms]" << endl;

  TCanvas * c_tau_L = new TCanvas("c_tau_L", "c_tau_L", 800, 600);
  g_tau_L_LLRF->Draw("APL");
  f_tau_L_truth->Draw("same");
  f_tau_L_lab->Draw("same");
 
  theApp.Run();

  return 0;
}