void currents(void)
{
// 
// plot the steady state solution to the heat transer problem with a source from [-sigma,sigma] in a uniform cooling bath
//

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

gROOT->Reset();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
gStyle->SetPalette(1,0);
gStyle->SetOptStat(kFALSE);
gStyle->SetOptFit(kFALSE);
// 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 jmax 1001;

   // 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 xlogfactor=1;
   Double_t dummyx[npts]={0,10,500*xlogfactor};
   Double_t dummyy[npts]={-10000,-10000,-10000};

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

   c1->SetGridx();
   c1->SetGridy();
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   c1->Divide(2,2);
   // c1->SetLogy();

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


   Double_t xmin=0;
   Double_t xmax=200;
   Double_t ymin=-3000;
   Double_t ymax=3000;
   Double_t gamma0=0.1;
   Double_t L1= 10.3*(4599./4600);  // inductance of solenoid, units are H minus one coil
   Double_t L2 =7.3e-6;  // self-inductance of loop, units are H
   // Double_t L3 = L1/4600; // self-inductance of one coil due to solenoid
   Double_t L3 = 0;
   Double_t M=sqrt(L1*L2); // mutual inductance between solenoid and one coil/gamma . 0 < gamma < 1.
   // Double_t M=1;
   Double_t Rd = 0.06; // solenoid dump resistor, units are Ohms
   Double_t sigma = 0.0025; // half length of shorted region, units are m
   Double_t area = 4*sigma*sigma;  // area of shorted region, units are m2
   Double_t rho = 1.7e-8;  // cu resistivity at room temperature, units are Ohm-m
   Double_t rfactor = 2./75; // magneto-resistive effect x /1/RRR (cf Brindza e-mail of 4/30/2010
   Double_t Rs = rho*rfactor*2*sigma/area; // resistance of short, depends on resistivity at 4.4 K, length and area of shorted region.
   Double_t Rs = Rs/10;
   Double_t I1_0 = 1500;   //Value of current in solenoid loop at  t=0, units are A
   Double_t I2_0 = 0;      // Value of induced current in shorted loop at t=0, units are A
  
   printf ("L1=%g, L2=%g, M=%g, gamma0=%g, sigma=%g, area=%g, Rd=%g, Rs=%g\n",L1,L2,M,gamma0,sigma,area,Rd,Rs);

   // dummy to draw axes

   TGraph *dummy = new TGraph (npts,dummyx,dummyy);
   TLegend *leg = new TLegend(0.50,0.75,0.85,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("Time (s)");
   dummy->GetYaxis()->SetTitle("Current (A)");
   dummy->GetXaxis()->SetNdivisions(505);
   dummy->Draw("Ap");

   // plot current in each loop as a function of time

   TF1 *I1 = new TF1("I1_func",I1_func,xmin,xlogfactor*xmax,8);
   TF1 *I2 = new TF1("I2_func",I2_func,xmin,xlogfactor*xmax,8);

   Double_t gamma = gamma0;
   /*Double_t A = (L1+L2)/Rd;
   Double_t B = (L2 - gamma*M)/Rd;
   Double_t C = (L2 - gamma*M)/Rs;
   Double_t D = (L2+L3)/Rs;*/

   Double_t A = L1/Rd;
   Double_t B = gamma*M/Rd;
   Double_t C = gamma*M/Rs;
   Double_t D = L2/Rs;
   Double_t lam1 = 0.5*(A+D) + 0.5*sqrt((A-D)*(A-D) + 4*B*C);
   Double_t lam2 = 0.5*(A+D) - 0.5*sqrt((A-D)*(A-D) + 4*B*C);

   printf ("A=%g, B=%g, C=%g, D=%g, lam1=%g, lam2=%g\n",A,B,C,D,lam1,lam2);

   // check transformation matrix, use left eigenvalues

   Double_t Alam1= A - lam1;
   Double_t Alam1s= Alam1*Alam1;
   Double_t Alam2= A - lam2;
   Double_t Alam2s= Alam2*Alam2;
   Double_t C2 = C*C;

   Double_t U1 = -C/sqrt(Alam1s + C2);
   Double_t U2 = Alam1/sqrt(Alam1s+C2);
   Double_t U3 = -C/sqrt(Alam2s + C2);
   Double_t U4 = Alam2/sqrt(Alam2s+C2);
   printf (" U1=%g, U2=%g, U3=%g, U4=%g \n",U1,U2,U3,U4);

   // compute U  inverse

   Double_t Udet = U1*U4 - U3*U2;
   Double_t UI1 = U4/Udet;
   Double_t UI2 = -U2/Udet;
   Double_t UI3 =  -U3/Udet;
   Double_t UI4 = U1/Udet;
   /*Double_t UI1 = U1;
   Double_t UI2 = U3;
   Double_t UI3 =  U2;
   Double_t UI4 = U4;*/
   printf (" UI1=%g, UI2=%g, UI3=%g, UI4=%g Udet=%g\n",UI1,UI2,UI3,UI4,Udet);

   // unity

   Double_t Id1 = U1*UI1 + U2*UI3;
   Double_t Id2 = U1*UI2 + U2*UI4;
   Double_t Id3 = U3*UI1 + U4*UI3;
   Double_t Id4 = U3*UI2 + U4*UI4;

   printf (" Unit: Id1=%g, Id2=%g, Id3=%g, Id4=%g \n",Id1,Id2,Id3,Id4);

   Double_t UM1 = U1*A + U2*C;
   Double_t UM2 = U1*B + U2*D;
   Double_t UM3 = U3*A + U4*C;
   Double_t UM4 = U3*B + U4*D;

   printf (" UM: UM1=%g, UM2=%g, UM3=%g, UM4=%g \n",UM1,UM2,UM3,UM4);

   Double_t D1 = UM1*UI1 + UM2*UI3;
   Double_t D2 = UM1*UI2 + UM2*UI4;
   Double_t D3 = UM3*UI1 + UM4*UI3;
   Double_t D4 = UM3*UI2 + UM4*UI4;

   printf(" Diagonal: D1=%g, D2=%g, D3=%g, D4=%g \n",D1,D2,D3,D4);


   // check transformation matrix use right eigenvalues

   printf("\n Use right eigenvalues\n");

   Double_t Alam1= A - lam1;
   Double_t Alam1s= Alam1*Alam1;
   Double_t Alam2= A - lam2;
   Double_t Alam2s= Alam2*Alam2;
   Double_t B2 = B*B;
   printf ("Alam1=%g, Alam2=%g, Alam1s=%g, Alam2s=%g, B2=%g\n",Alam1,Alam2,Alam1s,Alam2s,B2);

   Double_t V1 = -B/sqrt(Alam1s + B2);
   Double_t V2 = -B/sqrt(Alam2s + B2);
   Double_t V3 = Alam1/sqrt(Alam1s+B2);
   Double_t V4 = Alam2/sqrt(Alam2s+B2);
   /*Double_t V1 = 1;
   Double_t V2 = 1;
   Double_t V3 = -Alam1/B;
   Double_t V4 = -Alam2/B;*/
   printf (" V1=%g, V2=%g, V3=%g, V4=%g \n",V1,V2,V3,V4);

   // compute V inverse

   Double_t Vdet = V1*V4 - V3*V2;
   Double_t VI1 = V4/Vdet;
   Double_t VI2 = -V2/Vdet;
   Double_t VI3 = -V3/Vdet;
   Double_t VI4 = V1/Vdet;
   printf (" VI1=%g, VI2=%g, VI3=%g, VI4=%g Vdet=%g\n",VI1,VI2,VI3,VI4,Vdet);

   // unity

   Double_t Id1 = V1*VI1 + V2*VI3;
   Double_t Id2 = V1*VI2 + V2*VI4;
   Double_t Id3 = V3*VI1 + V4*VI3;
   Double_t Id4 = V3*VI2 + V4*VI4;

   printf (" Unit: Id1=%g, Id2=%g, Id3=%g, Id4=%g \n",Id1,Id2,Id3,Id4);

   Double_t VIM1 = VI1*A + VI2*C;
   Double_t VIM2 = VI1*B + VI2*D;
   Double_t VIM3 = VI3*A + VI4*C;
   Double_t VIM4 = VI3*B + VI4*D;

   printf (" VIM: VIM1=%g, VIM2=%g, VIM3=%g, VIM4=%g \n",VIM1,VIM2,VIM3,VIM4);

   Double_t D1 = VIM1*V1 + VIM2*V3;
   Double_t D2 = VIM1*V2 + VIM2*V4;
   Double_t D3 = VIM3*V1 + VIM4*V3;
   Double_t D4 = VIM3*V2 + VIM4*V4;

   printf(" Diagonal: D1=%g, D2=%g, D3=%g, D4=%g \n",D1,D2,D3,D4);

  // test solution to characteristic equation

  Double_t cheq1 = (A-lam1)*(D-lam1) - B*C;
  Double_t cheq2 = (A-lam2)*(D-lam2) - B*C;

  printf ("Characteristic equation check: cheq1=%g, cheq2=%g\n\n",cheq1,cheq2);

   // test orthogonality of left and right eigenvectors

  Double_t test12 = U1*V2 + U2*V4;
  Double_t test21 = U3*V1 + U4*V3;

   // check eigevectors

   Double_t MV1 = A*V1 + B*V3;
   Double_t MV2 = A*V2 + B*V4;
   Double_t MV3 = C*V1 + D*V3;
   Double_t MV4 = C*V2 + D*V4;

   printf (" MV: MV1=%g, MV2=%g, MV3=%g, MV4=%g \n",MV1,MV2,MV3,MV4);
   printf (" MV/V: MV1=%g, MV2=%g, MV3=%g, MV4=%g \n",MV1/V1,MV2/V2,MV3/V3,MV4/V4);

  printf ("Eigenvector test: test12=%g, test21=%g\n",test12,test21);

   I1->SetParameter(0,V1);
   I1->SetParameter(1,V2);
   I1->SetParameter(2,V3);
   I1->SetParameter(3,V4);
   I1->SetParameter(4,lam1);
   I1->SetParameter(5,lam2);
   I1->SetParameter(6,VI1*I1_0+VI2*I2_0);
   I1->SetParameter(7,VI3*I1_0+VI4*I2_0);


   sprintf (string,"Solenoid current\n");
   TF1 *I1a = I1->DrawCopy("SameC");
   leg->AddEntry(I1a,string,"l");
   I1a->SetLineColor(2);

   I2->SetParameter(0,V1);
   I2->SetParameter(1,V2);
   I2->SetParameter(2,V3);
   I2->SetParameter(3,V4);
   I2->SetParameter(4,lam1);
   I2->SetParameter(5,lam2);
   I2->SetParameter(6,VI1*I1_0+VI2*I2_0);
   I2->SetParameter(7,VI3*I1_0+VI4*I2_0);

   sprintf (string,"Single loop current\n",gamma);
   TF1 *I2a = I2->DrawCopy("SameC");
   leg->AddEntry(I2a,string,"l");
   I2a->SetLineColor(4);

   leg->Draw();

   sprintf (string,"L_{1}=%g H\n",L1);
   t1 = new TLatex(0.20,0.49,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf (string,"L_{2}=%g H\n",L2);
   t1 = new TLatex(0.20,0.45,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf (string,"L_{3}=%g H\n",L3);
   t1 = new TLatex(0.20,0.41,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf (string,"M=%g H\n",gamma*M);
   t1 = new TLatex(0.20,0.37,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf (string,"#gamma=%g\n",gamma);
   t1 = new TLatex(0.20,0.33,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf (string,"R_{dump}=%g #Omega\n",Rd);
   t1 = new TLatex(0.20,0.29,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf (string,"R_{short}=%g #Omega\n",Rs);
   t1 = new TLatex(0.20,0.25,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf (string,"#lambda_{1}=%g s\n",lam1);
   t1 = new TLatex(0.20,0.21,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf (string,"#lambda_{2}=%g s\n",lam2);
   t1 = new TLatex(0.20,0.17,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
  

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

   Double_t ymin=10;
   Double_t ymax=500000*xlogfactor;
   TGraph *dummy2 = new TGraph (npts,dummyx,dummyy);
   dummy2->SetMarkerColor(1);
   dummy2->SetMarkerStyle(21);
   dummy2->SetTitle("");
   dummy2->GetXaxis()->SetRangeUser(xmin,xlogfactor*xmax);
   dummy2->GetYaxis()->SetRangeUser(ymin,ymax);
   dummy2->GetXaxis()->SetTitleSize(0.04);
   dummy2->GetYaxis()->SetTitleSize(0.04);
   dummy2->GetYaxis()->SetTitleOffset(1.5);
   dummy2->GetXaxis()->SetTitle("Time (s)");
   dummy2->GetYaxis()->SetTitle("Current (A)");
   dummy2->GetXaxis()->SetNdivisions(505);
   dummy2->Draw("Ap");

   sprintf (string,"Solenoid current\n");
   TF1 *I1a = I1->DrawCopy("SameC");
   I1a->SetLineColor(2);
   sprintf (string,"Single loop current\n",gamma);
   TF1 *I2a = I2->DrawCopy("SameC");
   I2a->SetLineColor(4);
   leg->Draw();

   // compute power dissipated

   Double_t Power1[jmax];
   Double_t Power2[jmax];
   Double_t time[jmax];
   Double_t Energy1=0;
   Double_t Energy2=0;

   Double_t deltat = (xmax*xlogfactor-xmin)/jmax;
   for (j=0;j<jmax;j++) {
        t = xmin + j*(xmax*xlogfactor-xmin)/jmax;
        time[j] = t;
        Double_t current = I1->Eval(t);
        Power1[j] = current*current*Rd;
        Energy1 = Energy1 + Power1[j]*deltat;
        Double_t current = I2->Eval(t);
        Power2[j] = current*current*Rs;
        Energy2 = Energy2 + Power2[j]*deltat;
        // printf("time=%g, Power1=%g, Power2=%g\n",time[j],Power1[j],Power2[j]);
        }

   Double_t Estored=0.5*L1*I1_0*I1_0;
   printf ("Estored=%g, Energy1=%g, Energy2=%g\n",Estored,Energy1,Energy2);

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

   Double_t ymin=0.1;
   Double_t ymax=200000;
   TGraph *P1 = new TGraph (jmax,time,Power1);
   TLegend *leg = new TLegend(0.50,0.75,0.85,0.9);
   P1->SetMarkerColor(2);
   P1->SetMarkerStyle(21);
   P1->SetMarkerSize(0.35);
   P1->SetTitle("");
   P1->GetXaxis()->SetRangeUser(xmin,xmax*xlogfactor);
   P1->GetYaxis()->SetRangeUser(ymin,ymax);
   P1->GetXaxis()->SetTitleSize(0.04);
   P1->GetYaxis()->SetTitleSize(0.04);
   P1->GetYaxis()->SetTitleOffset(1.5);
   P1->GetXaxis()->SetTitle("Time (s)");
   P1->GetYaxis()->SetTitle("Power (W)");
   P1->GetXaxis()->SetNdivisions(505);
   P1->Draw("Ap");

   sprintf (string,"Stored Energy=%g J\n",Estored);
   t1 = new TLatex(0.30,0.80,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw(); 
   sprintf (string,"Energy Solenoid=%g J\n",Energy1);
   t1 = new TLatex(0.30,0.75,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw(); 
   sprintf (string,"Energy Loop=%g J\n",Energy2);
   t1 = new TLatex(0.30,0.70,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw(); 

   TGraph *P2 = new TGraph (jmax,time,Power2);
   TLegend *leg = new TLegend(0.50,0.75,0.85,0.9);
   P2->SetMarkerColor(4);
   P2->SetMarkerStyle(21);
   P2->SetMarkerSize(0.35);
   P2->Draw("Samep");


   // save canvas to file

   Int_t ig = gamma*1000;
   sprintf(filename,"currents_c1_%d.pdf",ig);
   c1->SaveAs(filename);
   sprintf(filename,"currents_c1_%d.png",ig);
   c1->SaveAs(filename);


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

   // Compute current in solenoid as a function of time

   Double_t V1=par[0];
   Double_t V2=par[1];
   Double_t V3=par[2];
   Double_t V4=par[3];
   Double_t lam1=par[4];
   Double_t lam2=par[5];
   Double_t I1_0=par[6];
   Double_t I2_0=par[7];

   Double_t t1=x[0];
   Double_t pi=3.14159;

   char string[256];

   Double_t func;

   Double_t I1p = I1_0*exp(-t1/lam1);
   Double_t I2p = I2_0*exp(-t1/lam2);

   func = V1*I1p + V2*I2p;
   // func = V3*I1p + V4*I2p;

   sprintf (string,"I1: t1=%f V1=%g, V2=%g, V3=%g, V4=%g, I1_0=%g, I2_0=%g,  lam1=%g, lam2=%g, func=%f\n",t1,V1,V2,V3,V4,I1_0,I2_0,lam1,lam2,func);
   // printf ("string=%s",string);

   return func;
}
Double_t I2_func (Double_t *x, Double_t *par) {

   // Compute current in solenoid as a function of time

   Double_t V1=par[0];
   Double_t V2=par[1];
   Double_t V3=par[2];
   Double_t V4=par[3];
   Double_t lam1=par[4];
   Double_t lam2=par[5];
   Double_t I1_0=par[6];
   Double_t I2_0=par[7];

   Double_t t1=x[0];
   Double_t pi=3.14159;

   char string[256];

   Double_t func;

   Double_t I1p = I1_0*exp(-t1/lam1);
   Double_t I2p = I2_0*exp(-t1/lam2);

   // func = V1*I1p + V2*I2p;
   func = V3*I1p + V4*I2p;

   sprintf (string,"I2: t1=%f V1=%g, V2=%g, V3=%g, V4=%g, I1_0=%g, I2_0=%g,  lam1=%g, lam2=%g, func=%f\n",t1,V1,V2,V3,V4,I1_0,I2_0,lam1,lam2,func);
   // printf ("string=%s",string);

   return func;
}