void currents0(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,jshort;
   
   #define npts 3;
   #define jmax 1001;
   #define pi 3.14159; 
   #define nshorts 4;

   // 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 currents0",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.1;
   Double_t xmax=1000;
   Double_t ymin=-1000;
   Double_t ymax=30000;
   Double_t gamma0=0.34;
   Double_t turns=4600;
   Double_t climit=4000;

   Double_t L1= 26. ;  // inductance of solenoid, units are H minus one coil
   Double_t L2 =7.3e-6;  // self-inductance of loop, units are H, computed from single loop equation
   Double_t L2 =L1/(turns*turns); // scale to inductance of full solenoid divided by the number of turns.
   Double_t L2=10e-6;  // calculation by Eugene for the self-inductance.
   // 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 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 rf=0.1;
   Double_t Rs = Rs*rf;
   Double_t Rsolenoid = 1.1;    // assume 1.1 m for radius of solenoid
   Double_t conductor_area=0.008*0.008; // assume 8 mm square conductor
   Double_t Rl = rho*rfactor*2*pi*Rsolenoid/conductor_area;  //Dimension is Ohm*m , value taken from e-mail from Eugene. rho(2/75)6/0.008^2 is 4e-5 Ohm*m
   // Double_t Rl = 1e-9; // superconducting
   Double_t Rl_super = 1e-9;  // use this value for "superconducting" to avoid numerical inversion problems.
   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);
   TLegend *leg1 = new TLegend(0.16,0.65,0.55,0.85);
   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;

   // loop over values of the short resistance

  Rs = Rs /10;

   for (jshort=0; jshort<nshorts; jshort++) {

   Rs = 10*Rs;
   printf ("jshort=%d, nshorts=%d, Rs=%g\n",jshort,nshorts,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);

   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);

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


   sprintf (string,"L_{1}=%g H\n",L1);
   t1 = new TLatex(0.20,0.57,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.54,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.51,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.51,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf (string,"#gamma=%g\n",gamma);
   t1 = new TLatex(0.20,0.47,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.43,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.39,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.39,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.35,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=100;
   Double_t ymax=1000000*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");
  
    } 

   c1->cd(1);

   sprintf (string,"Short current\n",gamma);
   TF1 *I2a = I2->DrawCopy("SameC");
   if (jshort == 0) {
      I2a->SetLineColor(1);
      }
   elseif (jshort == 1) {
      I2a->SetLineColor(4);
      }
   elseif (jshort ==2) {
      I2a->SetLineColor(3);
      }
   elseif (jshort == 3) {
      I2a->SetLineColor(5);
   }
   leg->AddEntry(I2a,string,"l");

   leg1->Draw();


   c1->cd(2);

   sprintf (string,"Dump current\n");
   TF1 *I1a = I1->DrawCopy("SameC");
   I1a->SetLineColor(2);
   sprintf (string,"Short current\n",gamma);
   TF1 *I2a = I2->DrawCopy("SameC");
   if (jshort == 0) {
      I2a->SetLineColor(1);
      }
   elseif (jshort == 1) {
      I2a->SetLineColor(4);
      }
   elseif (jshort ==2) {
      I2a->SetLineColor(3);
      }
   elseif (jshort == 3) {
      I2a->SetLineColor(5);
   }
   leg1->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++) {
        Double_t 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->SetLineColor(2);
   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);
   if (jshort == 0) P1->Draw("Ap");

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

   TGraph *P2 = new TGraph (jmax,time,Power2);
   TLegend *leg = new TLegend(0.50,0.75,0.85,0.9);
   sprintf (string,"R_{s}=%g\n",Rs);
   if (jshort == 0) {
      leg1->AddEntry(P1,"Solenoid","l");
      P2->SetMarkerColor(1);
      P2->SetLineColor(1);
      leg1->AddEntry(P2,string,"l");
      }
   elseif (jshort == 1) {
      P2->SetMarkerColor(4);
      P2->SetLineColor(4);
      leg1->AddEntry(P2,string,"l");
      }
   elseif (jshort ==2) {
      P2->SetMarkerColor(3);
      P2->SetLineColor(3);
      leg1->AddEntry(P2,string,"l");
      }
   elseif (jshort == 3) {
      P2->SetMarkerColor(5);
      P2->SetLineColor(5);
      leg1->AddEntry(P2,string,"l");
   }

   P2->SetMarkerStyle(21);
   P2->SetMarkerSize(0.35);
   P2->DrawClone("Samep");

  printf ("End of Loop: jshort=%d, nshorts=%d, Rs=%g\n",jshort,nshorts,Rs);

   }

  leg1->Draw();


   // save canvas to file

   Int_t ig = gamma*1000;
   Int_t ir = rf*10;
   sprintf(filename,"currents0_c1_%d_%d.eps",ig,ir);
   c1->SaveAs(filename);
   sprintf(filename,"currents0_c1_%d_%d.png",ig,ir);
   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;
}