void currents2(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 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]={0,10,500};
   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->SetLogy();


   Double_t xmin=0;
   Double_t xmax=200;
   Double_t ymin=-6000;
   Double_t ymax=6000;
   Double_t gamma0=0.1;
   Double_t L1= 10.3;  // inductance of solenoid, units are H
   Double_t L2 =7.3e-6;  // self-inductance of loop, units are H
   Double_t M=sqrt(L1*L2); // mutual inductance between solenoid and one coil/gamma . 0 < gamma < 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 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->Draw("Ap");
   dummy->GetXaxis()->SetNdivisions(505);

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

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

   Double_t gamma = gamma0;
   Double_t A = L1/Rd;
   Double_t B = gamma*M/Rd;
   Double_t C = L1/Rd + gamma*M/Rs;
   // Double_t B = C/2;
   Double_t D = L2/Rs + gamma*M/Rd;
   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,A);
   I1->SetParameter(4,I1_0);
   I1->SetParameter(5,I2_0);
   I1->SetParameter(1,lam1);
   I1->SetParameter(2,lam2);
   I1->SetParameter(3,B);

   sprintf (string,"I1a: #gamma=%g \n",gamma);
   TF1 *I1a = I1->DrawCopy("SameC");
   leg->AddEntry(I1a,string,"l");
   I1a->SetLineColor(2);

   I2->SetParameter(0,A);
   I2->SetParameter(4,I1_0);
   I2->SetParameter(5,I2_0);
   I2->SetParameter(1,lam1);
   I2->SetParameter(2,lam2);
   I2->SetParameter(3,B);

   sprintf (string,"I2a: #gamma=%g \n",gamma);
   TF1 *I2a = I2->DrawCopy("SameC");
   leg->AddEntry(I2a,string,"l");
   I2a->SetLineColor(4);

   leg->Draw();

   sprintf (string,"#lambda_{1}=%g s\n",lam1);
   t1 = new TLatex(0.20,0.85,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.81,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
  
   sprintf(filename,"currents_c1.pdf");
   c1->SaveAs(filename);
   sprintf(filename,"currents_c1.png");
   c1->SaveAs(filename);

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

   // Compute current in solenoid as a function of time

   Double_t A=par[0];
   Double_t lam1=par[1];
   Double_t lam2=par[2];
   Double_t B=par[3];
   Double_t I1_0=par[4];
   Double_t I2_0=par[5];

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

   char string[256];

   Double_t func;
   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;

   Double_t I1p = I1_0*exp(-t1/lam1);
   Double_t I2p = I2_0*exp(-t1/lam2);
   Double_t factor = sqrt(1+Alam1s/B2)/(1-Alam1/Alam2);

   func = factor*(I1p - B*I2p/Alam2);

   sprintf (string,"I1: t1=%f A=%g, Alam1=%g, Alam2=%g, B=%g, I1_0=%g, I2_0=%g,  factor=%g, func=%f\n",t1,A,Alam1,Alam2,B,I1_0,I2_0,factor,func);
   //printf ("string=%s",string);

   // return func;
   // printf("I1: I1p=%g, I2p=%g\n\n",I1p,I2p); 
   return I1p;
}
Double_t I2_func (Double_t *x, Double_t *par) {

   // Compute current in solenoid as a function of time

   Double_t A=par[0];
   Double_t lam1=par[1];
   Double_t lam2=par[2];
   Double_t B=par[3];
   Double_t I1_0=par[4];
   Double_t I2_0=par[5];

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

   char string[256];

   Double_t func1,func2;
   Double_t Alam1= fabs(A - lam1);
   Double_t Alam1s= Alam1*Alam1;
   Double_t Alam2= fabs(A - lam2);
   Double_t Alam2s= Alam2*Alam2;
   Double_t B2 = B*B;

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

   Double_t factor1 = sqrt(1+Alam1s/B2)/(1-Alam1/Alam2);
   func1 = factor1*(I1p - B*I2p/Alam2);
   Double_t factor2 = sqrt(1+Alam2s/B2)/(1-Alam2/Alam1);
   func2 = factor2*(I1p - B*I2p/Alam1);

   sprintf (string,"I2: t1=%f A=%g, Alam1=%g, Alam2=%g, B=%g, I1_0=%g, I2_0=%g,  factor2=%g, func=%f\n",t1,A,Alam1,Alam2,B,I1_0,I2_0,factor2,func2-func1);
   //printf ("string=%s",string);

  // return I2 - I1, which is the current flowing through Rs

   // compute determinant 

   Double_t U1 = B/sqrt(Alam1s + B2);
   Double_t U2 = B/sqrt(Alam2s + B2);
   Double_t U3 = Alam1/sqrt(Alam1s+B2);
   Double_t U4 = Alam2/sqrt(Alam2s+B2);
   Double_t logU4 = log(Alam2) - 0.5*log(Alam2s+B2);
   Double_t U4b = exp(logU4);
   Double_t det = U1*U4 - U2*U3;
   //printf (" U1=%g, U2=%g, U3=%g, U4=%g, U4b=%g, det=%g\n",U1,U2,U3,U4,U4b,det);

   // return func2 - func1;
   // printf("I1: I1p=%g, I2p=%g\n\n",I1p,I2p); 
   return I2p;
}