void critical_current(void)
{
// 
// plot the critical current in a NbTi SSC cable. Parameterization IEEE Trans on Applied SC, Vol 10, No 1, March 2000, p. 1054
// Jc = (C0/B) * b**alpha * (1-b)**beta * (1 - t**nexpon)**gamma.
// 
// see also Brindza e-mail Jul 14, 2010.
//

// #include <TMath.h>
#include <TRandom.h>

gROOT->Reset();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
gStyle->SetPalette(1,0);
gStyle->SetOptStat(kFALSE);
gStyle->SetOptFit(kTRUE);
gStyle->SetOptFit(1111);
gStyle->SetPadRightMargin(0.15);
gStyle->SetPadLeftMargin(0.15);
gStyle->SetPadBottomMargin(0.15);
gStyle->SetFillColor(0);
//

   char string[256];
   char filename[80];
   Int_t j,jj;
   
   #define npts 3;
   #define mpts 8;
   #define xpts 5;
   #define max 201;

   // input data set first point (actually -100) to 100 so that it is not plotted.

   // use offset on x-axis to make data visible and give negative y-values for reference (not plotted)

   Double_t dummyx[npts]={-1,0,10};
   Double_t dummyy[npts]={-1,-1,-1};

   // Data from Log Book 6/2/72 Alcorn notes

   Double_t Ic_B_Bperp[mpts]={2.55,2.2,1.9,-1,1.6,-1,1.4,-1};
   Double_t Ic_B_Bpar[mpts]={3.1,2.9,2.5,-1,2.1,-1,1.7,-1};
   Double_t Ic_B_coil[mpts]={-1,-1,2.25,1.95,-1,1.65,1.4,1.15};

   /*Double_t Ic_B_fit[xpts]={2.25,1.95,1.65,1.4,1.15};   
   Double_t Ic_Bfit[xpts]={5,5.3,5.7,6,6.5};*/

   Double_t Ic_B_fit[xpts]={3.1,2.9,2.5,2.1,1.7};   
   Double_t Ic_Bfit[xpts]={4,4.5,5,5.5,6};

   Double_t Ic_A_Bperp[mpts]={-1,3.0,2.75,-1,2.3,-1,1.95,-1};
   Double_t Ic_A_Bpar[mpts]={-1,-1,3.3,-1,2.75,-1,2.25,-1};
   Double_t Ic_Bfield[mpts]={4,4.5,5,5.3,5.5,5.7,6,6.5};

   //
   TCanvas *c1 = new TCanvas("c1","c1 critical_current",200,10,700,700);

   c1->SetGridx();
   c1->SetGridy();
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   // c1->SetLogx();
   c1->SetLogy();


   Double_t xmin=0.;
   Double_t xmax=10;
   // Double_t xmin=0.1;
   // Double_t xmax=0.5;
   Double_t ymin=0.1;
   Double_t ymax=50;

   // dummy to draw axes

   TGraph *dummy = new TGraph (npts,dummyx,dummyy);
   TLegend *leg = new TLegend(0.23,0.7,0.57,0.9);
   dummy->SetMarkerColor(1);
   dummy->SetMarkerStyle(21);
   t1 = new TLatex(0.20,0.91,"Shapes");
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw(); 
   dummy->SetTitle("");
   dummy->GetXaxis()->SetRangeUser(xmin,xmax);
   dummy->GetYaxis()->SetRangeUser(ymin,ymax);
   dummy->GetXaxis()->SetTitleSize(0.04);
   dummy->GetYaxis()->SetTitleSize(0.04);
   dummy->GetYaxis()->SetTitleOffset(1.5);
   dummy->GetXaxis()->SetTitle("Magnetic Field (T)");
   dummy->GetYaxis()->SetTitle("Critical Current (kA)");
   dummy->Draw("Ap");
   // dummy->GetXaxis()->SetNdivisions(505);


  // plot data

   /*TGraph *B_Bperp = new TGraph (mpts,Ic_Bfield,Ic_B_Bperp);
   B_Bperp->SetMarkerColor(2);
   B_Bperp->SetMarkerStyle(24);
   B_Bperp->Draw("Samep");
   sprintf (string,"Short Sample 'B' Bperp\n");
   leg->AddEntry(B_Bperp,string,"p");

   TGraph *B_Bpar = new TGraph (mpts,Ic_Bfield,Ic_B_Bpar);
   B_Bpar->SetMarkerColor(2);
   B_Bpar->SetMarkerStyle(25);
   B_Bpar->Draw("Samep");
   sprintf (string,"Short Sample 'B' Bpar\n");
   leg->AddEntry(B_Bpar,string,"p");

   TGraph *B_coil = new TGraph (mpts,Ic_Bfield,Ic_B_coil);
   B_coil->SetMarkerColor(2);
   B_coil->SetMarkerStyle(20);
   B_coil->Draw("Samep");
   sprintf (string,"Short Sample 'B' Test Coil\n");
   leg->AddEntry(B_coil,string,"p");*/

   /*TGraph *B_fit = new TGraph (xpts,Ic_Bfit,Ic_B_fit);
   B_fit->SetMarkerColor(2);
   B_fit->SetMarkerStyle(20);
   B_fit->Draw("Samep");
   sprintf (string,"Short Sample 'B' Test Coil\n");
   leg->AddEntry(B_fit,string,"p");*/

   TGraph *A_Bperp = new TGraph (mpts,Ic_Bfield,Ic_A_Bperp);
   A_Bperp->SetMarkerColor(4);
   A_Bperp->SetMarkerStyle(24);
   A_Bperp->Draw("Samep");
   sprintf (string,"Short Sample 'A' Bperp\n");
   leg->AddEntry(A_Bperp,string,"p");

   TGraph *A_Bpar = new TGraph (mpts,Ic_Bfield,Ic_A_Bpar);
   A_Bpar->SetMarkerColor(4);
   A_Bpar->SetMarkerStyle(25);
   A_Bpar->Draw("Samep");
   sprintf (string,"Short Sample 'A' Bpar\n");
   leg->AddEntry(A_Bpar,string,"p");

   // plot temperature function. Consider x as a parameter and plot time dependence.

   Double_t Tnbp=4.2;    // nominal boiling point units are K
   // Double_t Toperating=4.42;    // nominal operating point. This is the 2-phase temperature of LHE at 1.2 atm.
   Double_t Toperating=4.42;    // nominal operating point. This is the 2-phase temperature of LHE at 1.2 atm.
   Double_t C0=27;  // overall normalization in Tesla
   Double_t C0=75;  // overall normalization in Tesla Type A
   // Double_t C0=55;  // overall normalization in Tesla Type B
   Double_t Tc0=9.2;  // maximum critical temperature at B=0
   Double_t Bc20=14.5;    // maximum upper crititcal field at T=0 in Tesla
   Double_t nexpon=1.7;   // Temperature dependence of critical field, assumed constant in Bottura
   /*Double_t alpha=0.57;   // first exponent in parameterization of Pinning force
   Double_t beta=0.9;      // second exponent in parameterization of Pinning force
   Double_t gamma=2.32;    // Exponent for temperature scaling law.*/

   Double_t alpha=0.65;   // first exponent in parameterization of Pinning force
   Double_t beta=1.0;      // second exponent in parameterization of Pinning force
   Double_t gamma=1.8;    // Exponent for temperature scaling law.

   // function computed at a fixed temperature T

   Double_t T=Toperating;    // 
   TF1 *Ic = new TF1("Ic_func",Ic_func,xmin,xmax,9);
   Ic->SetParameters(T,Tnbp,C0,Tc0,Bc20,alpha,beta,gamma,nexpon); 

   Double_t T=Toperating;    // 
   sprintf (string,"T=%g K\n",T);
   TF1 *Ic2 = Ic->DrawCopy("Samel");
   Ic2->SetParameter(0,T);
   Ic2->SetLineColor(2);
   leg->AddEntry(Ic2,string,"l");

   Double_t T=1.2*Toperating;    // 
   sprintf (string,"T=%g K\n",1.1*T);
   TF1 *Ic3 = Ic->DrawCopy("Samel");
   Ic3->SetParameter(0,T);
   Ic3->SetLineColor(1);
   leg->AddEntry(Ic3,string,"l");

   Double_t T=1.4*Toperating;    // 
   sprintf (string,"T=%g K\n",T);
   TF1 *Ic4 = Ic->DrawCopy("Samel");
   Ic4->SetParameter(0,T);
   Ic4->SetLineColor(4);
   leg->AddEntry(Ic4,string,"l");

   Double_t T=1.6*Toperating;    // 
   sprintf (string,"T=%g K\n",T);
   TF1 *Ic5 = Ic->DrawCopy("Samel");
   Ic5->SetParameter(0,T);
   Ic5->SetLineColor(3);
   leg->AddEntry(Ic5,string,"l");

   Double_t T=1.9*Toperating;    // 
   sprintf (string,"T=%g K\n",T);
   TF1 *Ic5 = Ic->DrawCopy("Samel");
   Ic5->SetParameter(0,T);
   Ic5->SetLineColor(5);
   leg->AddEntry(Ic5,string,"l");

   leg->Draw();

   sprintf (string,"T_{nbp}=%g K\n",Tnbp);
   t1 = new TLatex(0.6,0.88,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf (string,"C_{0}=%g T kA\n",C0);
   t1 = new TLatex(0.6,0.84,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf (string,"T_{c0}=%g K\n",Tc0);
   t1 = new TLatex(0.6,0.80,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf (string,"B_{c20}=%g T\n",Bc20);
   t1 = new TLatex(0.60,0.76,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf (string,"#alpha=%g\n",alpha);
   t1 = new TLatex(0.60,0.72,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf (string,"#beta=%g\n",beta);
   t1 = new TLatex(0.60,0.68,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf (string,"#gamma=%g\n",gamma);
   t1 = new TLatex(0.60,0.64,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf (string,"n=%g\n",nexpon);
   t1 = new TLatex(0.60,0.60,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
  
   Ic->FixParameter(0,Toperating);
   Ic->FixParameter(1,Tnbp);
   // Ic->FixParameter(2,C0);
   Ic->FixParameter(3,Tc0);
   Ic->FixParameter(4,Bc20);
   Ic->FixParameter(5,alpha);
   // Ic->FixParameter(6,beta);
   // Ic->FixParameter(7,gamma);
   Ic->FixParameter(8,nexpon);
   // B_fit->Fit(Ic,"","",3.5,6.5);

   // Take normalization at 5 T

   TF1 *Bfil = new TF1("Bfilament","(x*87*2*pi*0.18e-3/2)/(pi*4e-7*1e3)",0,2);   // Type A
   // TF1 *Bfil = new TF1("Bfilament","(x*87*2*pi*0.127e-3/2)/(pi*4e-7*1e3)",0,2);   // Type B
   Bfil->SetLineStyle(2);
   Bfil->SetLineColor(1);
   Bfil->SetLineWidth(2);
   Bfil->Draw("Samec");

   TLine *line = new TLine (0.5,ymin,0.5,ymax);
   line->SetLineStyle(2);
   line->SetLineColor(1);
   line->SetLineWidth(2);
   //line->Draw();

   Ic->Eval(5);

   sprintf(filename,"critical_current_c1.eps");
   c1->SaveAs(filename);
   sprintf(filename,"critical_current_c1.png");
   c1->SaveAs(filename);

}
Double_t Ic_func (Double_t *x, Double_t *par) {

   // Compute critical current in a NbTi superconductor as a function of B, at a fixed temperature T.
  // Parameters are Tc0, Bc0, alpha, beta and gamma in the parameterization.
  // 
 // plot the critical current in a NbTi SSC cable. Parameterization IEEE Trans on Applied SC, Vol 10, No 1, March 2000, p. 1054
 // Jc = (C0/B) * b**alpha * (1-b)**beta * (1 - t**nexpon)**gamma.


   Double_t T=par[0];
   Double_t Tnbp=par[1];
   Double_t C0=par[2];
   Double_t Tc0=par[3];
   Double_t Bc20=par[4];
   Double_t alpha=par[5];
   Double_t beta=par[6];
   Double_t gamma=par[7];
   Double_t nexpon=par[8];
   Double_t B=x[0];

   char string[256];

   Double_t func;

   Double_t Bc2 = Bc20 * (1 - pow(T/Tc0,nexpon));
   if (Bc2 >= Bc20) Bc2 = 0.999999*Bc20;
   Double_t br = B/Bc2;   // reduced bfield
   if (br >= 1) br = 0.9999999;
   Double_t Fp = C0 * pow(br,alpha) * pow(1-br,beta) * pow(Bc2/Bc20,gamma);

   Double_t func = Fp/B;

   /*sprintf (string,"B=%g, T=%g, br=%g, Fp=%g, func=%g\n",B,T,br,Fp,func);
   printf ("%s",string);*/

   return func;
}