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 10001;
   # define pi 3.14159; 

   // 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.1;
   Double_t xmax=1000;
   //Double_t xmin=0.5;
   //Double_t xmax=10;
   Double_t ymin=-3000;
   Double_t ymax=20000;
   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=100.;
   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
   // Double_t I1_0 = 1495;   //Value of current in solenoid loop at  t=0, units are A
   // Double_t I2_0 = 1000;      // 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");


   Double_t gamma = gamma0;

   // left side equation

   // Loop over configurations when Rl is superconducting and Rl conducting mode

   Double_t Rl0 = Rl;
   Double_t RV1_save;
   Double_t RV2_save;
   Double_t RV3_save;
   Double_t RV4_save;

   Double_t RVI1_save;
   Double_t RVI2_save;
   Double_t RVI3_save;
   Double_t RVI4_save;
   Double_t lam1_0;
   Double_t lam2_0;

   for (j=0;j<2;j++) {
   if (j == 0) {
       Rl = Rl_super;
       }
       else {
       Rl = Rl0;
       }

   Double_t A = L1;
   Double_t B = gamma*M;
   Double_t C = gamma*M;
   Double_t D = L2;

   Double_t B1 = Rd+Rl;
   Double_t B2= Rl;
   Double_t B3 = Rl;
   Double_t B4 = Rs+Rl;
   Double_t rlam1 = 0.5*(B1+B4) + 0.5*sqrt((B1-B4)*(B1-B4) + 4*B2*B3);
   Double_t rlam2 = 0.5*(B1+B4) -  0.5*sqrt((B1-B4)*(B1-B4) + 4*B2*B3);

   printf ("B1=%g, B2=%g, B3=%g, B4=%g, rlam1=%g, rlam2=%g\n",B1,B2,B3,B4,rlam1,rlam2);
  // test solution to characteristic equation

  Double_t cheq1 = (B1-rlam1)*(B4-rlam1) - B2*B3;
  Double_t cheq2 = (B1-rlam2)*(B4-rlam2) - B2*B3;

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

   Double_t Brlam1= B1 - rlam1;
   Double_t Brlam1s= Brlam1*Brlam1;
   Double_t Brlam2= B1 - rlam2;
   Double_t Brlam2s= Brlam2*Brlam2;
   Double_t Bs = B2*B2;
   printf ("Brlam1=%g, Brlam2=%g, Brlam1s=%g, Brlam2s=%g, Bs=%g\n",Brlam1,Brlam2,Brlam1s,Brlam2s,B2);

   Double_t R1 = -B2/sqrt(Brlam1s + Bs);
   Double_t R2 = -B2/sqrt(Brlam2s + Bs);
   Double_t R3 = Brlam1/sqrt(Brlam1s+Bs);
   Double_t R4 = Brlam2/sqrt(Brlam2s+Bs);
   printf (" R1=%g, R2=%g, R3=%g, R4=%g \n",R1,R2,R3,R4);

   // compute R inverse

   Double_t Rdet = R1*R4 - R3*R2;
   Double_t RI1 = R4/Rdet;
   Double_t RI2 = -R2/Rdet;
   Double_t RI3 = -R3/Rdet;
   Double_t RI4 = R1/Rdet;
   printf (" RI1=%g, RI2=%g, RI3=%g, RI4=%g Rdet=%g\n",RI1,RI2,RI3,RI4,Rdet);

   // check eigevectors

   Double_t MBR1 = B1*R1 + B2*R3;
   Double_t MBR2 = B1*R2 + B2*R4;
   Double_t MBR3 = B3*R1 + B4*R3;
   Double_t MBR4 = B3*R2 + B4*R4;

   printf (" MBR: MBR1=%g, MBR2=%g, MBR3=%g, MBR4=%g \n",MBR1,MBR2,MBR3,MBR4);
   printf (" MBR/R: MBR1=%g, MBR2=%g, MBR3=%g, MBR4=%g \n",MBR1/R1,MBR2/R2,MBR3/R3,MBR4/R4);

   Double_t RIMBR1 = RI1*MBR1 + RI2*MBR3;
   Double_t RIMBR2 = RI1*MBR2 + RI2*MBR4;
   Double_t RIMBR3 = RI3*MBR1 + RI4*MBR3;
   Double_t RIMBR4 = RI3*MBR2 + RI4*MBR4;

   printf (" Diagonal: RIMBR1=%g, RIMBR2=%g, RIMBR3=%g, RIMBR4=%g \n",RIMBR1,RIMBR2,RIMBR3,RIMBR4);

   Double_t RIM1 = RI1*A + RI2*C;
   Double_t RIM2 = RI1*B + RI2*D;
   Double_t RIM3 = RI3*A + RI4*C;
   Double_t RIM4 = RI3*B + RI4*D;

   printf (" RIM: RIM1=%g, RIM2=%g, RIM3=%g, RIM4=%g \n",RIM1,RIM2,RIM3,RIM4);

   Double_t A1 = (RIM1*R1 + RIM2*R3)/rlam1;
   Double_t A2 = (RIM1*R2 + RIM2*R4)/rlam1;
   Double_t A3 = (RIM3*R1 + RIM4*R3)/rlam2;
   Double_t A4 = (RIM3*R2 + RIM4*R4)/rlam2;



   Double_t gamma = gamma0;

   Double_t lam1 = 0.5*(A1+A4) + 0.5*sqrt((A1-A4)*(A1-A4) + 4*A2*A3);
   Double_t lam2 = 0.5*(A1+A4) - 0.5*sqrt((A1-A4)*(A1-A4) + 4*A2*A3);

   printf("\n Use right eigenvalues\n");
   printf ("A1=%g, A2=%g, A3=%g, A4=%g, lam1=%g, lam2=%g\n",A1,A2,A3,A4,lam1,lam2);

   Double_t Alam1= A1 - lam1;
   Double_t Alam1s= Alam1*Alam1;
   Double_t Alam2= A1 - lam2;
   Double_t Alam2s= Alam2*Alam2;
   Double_t As = A2*A2;
   printf ("Alam1=%g, Alam2=%g, Alam1s=%g, Alam2s=%g, As=%g\n",Alam1,Alam2,Alam1s,Alam2s,As);

   Double_t V1 = -A2/sqrt(Alam1s + As);
   Double_t V2 = -A2/sqrt(Alam2s + As);
   Double_t V3 = Alam1/sqrt(Alam1s+As);
   Double_t V4 = Alam2/sqrt(Alam2s+As);
   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*A1 + VI2*A3;
   Double_t VIM2 = VI1*A2 + VI2*A4;
   Double_t VIM3 = VI3*A1 + VI4*A3;
   Double_t VIM4 = VI3*A2 + VI4*A4;

   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 = (A1-lam1)*(A4-lam1) - A2*A3;
  Double_t cheq2 = (A1-lam2)*(A4-lam2) - A2*A3;

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

   // check eigevectors

   Double_t MV1 = A1*V1 + A2*V3;
   Double_t MV2 = A1*V2 + A2*V4;
   Double_t MV3 = A3*V1 + A4*V3;
   Double_t MV4 = A3*V2 + A4*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);

   Double_t RV1 = R1*V1+R2*V3;
   Double_t RV2 = R1*V2+R2*V4;
   Double_t RV3 = R3*V1+R4*V3;
   Double_t RV4 = R3*V2+R4*V4;

   Double_t RVI1 = VI1*RI1+VI2*RI3;
   Double_t RVI2 = VI1*RI2+VI2*RI4;
   Double_t RVI3 = VI3*RI1+VI4*RI3;
   Double_t RVI4 = VI3*RI2+VI4*RI4;

   printf ("RV1=%g, RV2=%g, RV3=%g, RV4=%g\n",RV1,RV2,RV3,RV4);
   printf ("RVI1=%g, RVI2=%g, RVI3=%g, RVI4=%g\n",RVI1,RVI2,RVI3,RVI4);

   if (j ==0) {
      Double_t RV1_save = RV1;
      Double_t RV2_save = RV2;
      Double_t RV3_save = RV3;
      Double_t RV4_save = RV4;

      Double_t RVI1_save = RVI1;
      Double_t RVI2_save = RVI2;
      Double_t RVI3_save = RVI3;
      Double_t RVI4_save = RVI4;
      Double_t lam1_0=lam1;
      Double_t lam2_0=lam2;
      }
   }

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

   TF1 *I1 = new TF1("I1_func",I1_func,xmin,xlogfactor*xmax,17);
   TF1 *I2 = new TF1("I2_func",I2_func,xmin,xlogfactor*xmax,17);
   TF1 *absI2 = new TF1("absI2_func",absI2_func,xmin,xlogfactor*xmax,17);
   TF1 *Isum = new TF1("Isum_func",Isum_func,xmin,xlogfactor*xmax,17);

   Double_t t;
   Double_t tlimit=xmax;
   Double_t I1limit=0;
   Double_t I2limit=0;

   I1->SetParameter(0,RV1);
   I1->SetParameter(1,RV2);
   I1->SetParameter(2,RV3);
   I1->SetParameter(3,RV4);
   I1->SetParameter(4,lam1);
   I1->SetParameter(5,lam2);
   I1->SetParameter(6,RVI1_save*I1_0+RVI2_save*I2_0);
   I1->SetParameter(7,RVI3_save*I1_0+RVI4_save*I2_0);
   I1->SetParameter(8,RV1_save);
   I1->SetParameter(9,RV2_save);
   I1->SetParameter(10,RV3_save);
   I1->SetParameter(11,RV4_save);
   I1->SetParameter(12,tlimit);
   I1->SetParameter(13,RVI1*I1limit+RVI2*I2limit);
   I1->SetParameter(14,RVI3*I1limit+RVI4*I2limit);
   I1->SetParameter(15,lam1_0);
   I1->SetParameter(16,lam2_0);

   I2->SetParameter(0,RV1);
   I2->SetParameter(1,RV2);
   I2->SetParameter(2,RV3);
   I2->SetParameter(3,RV4);
   I2->SetParameter(4,lam1);
   I2->SetParameter(5,lam2);
   I2->SetParameter(6,RVI1_save*I1_0+RVI2_save*I2_0);
   I2->SetParameter(7,RVI3_save*I1_0+RVI4_save*I2_0);
   I2->SetParameter(8,RV1_save);
   I2->SetParameter(9,RV2_save);
   I2->SetParameter(10,RV3_save);
   I2->SetParameter(11,RV4_save);
   I2->SetParameter(12,tlimit);
   I2->SetParameter(13,RVI1*I1limit+RVI2*I2limit);
   I2->SetParameter(14,RVI3*I1limit+RVI4*I2limit);
   I2->SetParameter(15,lam1_0);
   I2->SetParameter(16,lam2_0);

   Isum->SetParameter(0,RV1);
   Isum->SetParameter(1,RV2);
   Isum->SetParameter(2,RV3);
   Isum->SetParameter(3,RV4);
   Isum->SetParameter(4,lam1);
   Isum->SetParameter(5,lam2);
   Isum->SetParameter(6,RVI1_save*I1_0+RVI2_save*I2_0);
   Isum->SetParameter(7,RVI3_save*I1_0+RVI4_save*I2_0);
   Isum->SetParameter(8,RV1_save);
   Isum->SetParameter(9,RV2_save);
   Isum->SetParameter(10,RV3_save);
   Isum->SetParameter(11,RV4_save);
   Isum->SetParameter(12,tlimit);
   Isum->SetParameter(13,RVI1*I1limit+RVI2*I2limit);
   Isum->SetParameter(14,RVI3*I1limit+RVI4*I2limit);
   Isum->SetParameter(15,lam1_0);
   Isum->SetParameter(16,lam2_0);

   absI2->SetParameter(0,RV1);
   absI2->SetParameter(1,RV2);
   absI2->SetParameter(2,RV3);
   absI2->SetParameter(3,RV4);
   absI2->SetParameter(4,lam1);
   absI2->SetParameter(5,lam2);
   absI2->SetParameter(6,RVI1_save*I1_0+RVI2_save*I2_0);
   absI2->SetParameter(7,RVI3_save*I1_0+RVI4_save*I2_0);
   absI2->SetParameter(8,RV1_save);
   absI2->SetParameter(9,RV2_save);
   absI2->SetParameter(10,RV3_save);
   absI2->SetParameter(11,RV4_save);
   absI2->SetParameter(12,tlimit);
   absI2->SetParameter(13,RVI1*I1limit+RVI2*I2limit);
   absI2->SetParameter(14,RVI3*I1limit+RVI4*I2limit);
   absI2->SetParameter(15,lam1_0);
   absI2->SetParameter(16,lam2_0);

   // compute power dissipated

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

   // obtain limits when Rl becomes resistive

   for (j=0;j<jmax;j++) {
        t = xmin + j*(xmax*xlogfactor-xmin)/jmax;
        time[j] = t;
        Double_t current1 = I1->Eval(t);
        Double_t current2 = I2->Eval(t);
        if (abs(current1 + current2) > climit && I1limit==0) {
            tlimit = t;
            I2limit = current2;
            I1limit = current1;
           }
        }
       printf ("Values at limit: tlimit=%g, I1limit=%g, I2limit=%g\n",tlimit,I1limit,I2limit);


   // now update saved constants into function

   I1->SetParameter(12,tlimit);
   I1->SetParameter(13,RVI1*I1limit+RVI2*I2limit);
   I1->SetParameter(14,RVI3*I1limit+RVI4*I2limit);

   I2->SetParameter(12,tlimit);
   I2->SetParameter(13,RVI1*I1limit+RVI2*I2limit);
   I2->SetParameter(14,RVI3*I1limit+RVI4*I2limit);

   Isum->SetParameter(12,tlimit);
   Isum->SetParameter(13,RVI1*I1limit+RVI2*I2limit);
   Isum->SetParameter(14,RVI3*I1limit+RVI4*I2limit);

   absI2->SetParameter(12,tlimit);
   absI2->SetParameter(13,RVI1*I1limit+RVI2*I2limit);
   absI2->SetParameter(14,RVI3*I1limit+RVI4*I2limit);

   // compute power dissipated

   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 current1 = I1->Eval(t);
        Power1[j] = current1*current1*Rd;
        Energy1 = Energy1 + Power1[j]*deltat;
        Double_t current2 = I2->Eval(t);
        Power2[j] = current2*current2*Rs;
        Energy2 = Energy2 + Power2[j]*deltat;
        Double_t current3 = Isum->Eval(t);
        if (t < tlimit) {
               Power3[j] = 0;
            }
           else {
               Power3[j] = current3*current3*Rl;
           }
        Energy3 = Energy3 + Power3[j]*deltat;
        // printf("time=%g, current1=%g, current2=%g, current3=%g\n",time[j],current1,current2,current3);
        // printf("time=%g, Power1=%g, Power2=%g, Power3=%g\n",time[j],Power1[j],Power2[j],Power3[j]);
        }

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


   sprintf (string,"Solenoid current\n");
   TF1 *I1a = I1->DrawCopy("SameC");
   leg->AddEntry(I1a,string,"l");
   I1a->SetLineColor(1);
   sprintf (string,"Short current\n",gamma);
   TF1 *I2a = I2->DrawCopy("SameC");
   leg->AddEntry(I2a,string,"l");
   I2a->SetLineColor(4);
   TF1 *Isuma = Isum->DrawCopy("SameC");
   sprintf (string,"Current sum\n");
   leg->AddEntry(Isuma,string,"l");
   Isuma->SetLineColor(2);

   leg->Draw();

   sprintf (string,"L_{1}=%g H\n",L1);
   t1 = new TLatex(0.20,0.84,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.80,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.76,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf (string,"#gamma=%g\n",gamma);
   t1 = new TLatex(0.20,0.72,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.68,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.64,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf (string,"R_{loop}=%g #Omega\n",Rl);
   t1 = new TLatex(0.20,0.60,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.56,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.52,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=100000*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");

   TF1 *I1a = I1->DrawCopy("SameC");
   I1a->SetLineColor(1);
   TF1 *absI2a = absI2->DrawCopy("SameC");
   absI2a->SetLineColor(4);
   TF1 *Isumb = Isum->DrawCopy("SameC");
   Isumb->SetLineColor(2);
   leg->Draw();


   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.01;
   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(1);
   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.20,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.20,0.75,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   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);
   t1->Draw(); 
   sprintf (string,"Energy Loop=%g J\n",Energy3);
   t1 = new TLatex(0.20,0.65,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");

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

   // save canvas to file

   Int_t ig = gamma*1000;
   Int_t ir = rf*10;
   sprintf(filename,"currents_c1_%d_%d.eps",ig,ir);
   c1->SaveAs(filename);
   sprintf(filename,"currents_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 I1p_0=par[6];
   Double_t I2p_0=par[7];
   Double_t RV1_save=par[8];
   Double_t RV2_save=par[9];
   Double_t RV3_save=par[10];
   Double_t RV4_save=par[11];
   Double_t tlimit=par[12];
   Double_t I1plimit=par[13];
   Double_t I2plimit=par[14];
   Double_t lam1_0=par[15];
   Double_t lam2_0=par[16];

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

   char string[256];

   Double_t func;

   if (t1 < tlimit) {
       Double_t I1p = I1p_0*exp(-t1/lam1_0);
       Double_t I2p = I2p_0*exp(-t1/lam2_0);
       func = RV1_save*I1p + RV2_save*I2p;
       // func = RV3_save*I1p + RV4_save*I2p;
       }
   else {
       Double_t I1p = I1plimit*exp(-(t1-tlimit)/lam1);
       Double_t I2p = I2plimit*exp(-(t1-tlimit)/lam2);
       func = V1*I1p + V2*I2p;
       // func = V3*I1p + V4*I2p;
      }

   sprintf (string,"I1: t1=%f V1=%g, V2=%g, V3=%g, V4=%g, I1plimit=%g, I2plimit=%g,  lam1=%g, lam2=%g, func=%f\n",t1,V1,V2,V3,V4,I1plimit,I2plimit,lam1,lam2,func);
   //printf ("string=%s",string);
   sprintf (string,"I1: tlimit=%f RV1_save=%g, RV2_save=%g, RV3_save=%g, RV4_save=%g, I1p_0=%g, I2p_0=%g,  lam1=%g, lam2=%g, func=%f\n",tlimit,RV1_save,RV2_save,RV3_save,RV4_save,I1p_0,I2p_0,lam1_0,lam2_0,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 I1p_0=par[6];
   Double_t I2p_0=par[7];
   Double_t RV1_save=par[8];
   Double_t RV2_save=par[9];
   Double_t RV3_save=par[10];
   Double_t RV4_save=par[11];
   Double_t tlimit=par[12];
   Double_t I1plimit=par[13];
   Double_t I2plimit=par[14];
   Double_t lam1_0=par[15];
   Double_t lam2_0=par[16];

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

   char string[256];

   Double_t func;

   if (t1 < tlimit) {
       Double_t I1p = I1p_0*exp(-t1/lam1_0);
       Double_t I2p = I2p_0*exp(-t1/lam2_0);
       // func = RV1_save*I1p + RV2_save*I2p;
       func = RV3_save*I1p + RV4_save*I2p;
       }
   else {
       Double_t I1p = I1plimit*exp(-(t1-tlimit)/lam1);
       Double_t I2p = I2plimit*exp(-(t1-tlimit)/lam2);
       // func = V1*I1p + V2*I2p;
       func = V3*I1p + V4*I2p;
      }


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

   return func;
}
Double_t absI2_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 I1p_0=par[6];
   Double_t I2p_0=par[7];
   Double_t RV1_save=par[8];
   Double_t RV2_save=par[9];
   Double_t RV3_save=par[10];
   Double_t RV4_save=par[11];
   Double_t tlimit=par[12];
   Double_t I1plimit=par[13];
   Double_t I2plimit=par[14];
   Double_t lam1_0=par[15];
   Double_t lam2_0=par[16];

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

   char string[256];

   Double_t func;

   if (t1 < tlimit) {
       Double_t I1p = I1p_0*exp(-t1/lam1_0);
       Double_t I2p = I2p_0*exp(-t1/lam2_0);
       // func = RV1_save*I1p + RV2_save*I2p;
       func = RV3_save*I1p + RV4_save*I2p;
       }
   else {
       Double_t I1p = I1plimit*exp(-(t1-tlimit)/lam1);
       Double_t I2p = I2plimit*exp(-(t1-tlimit)/lam2);
       // func = V1*I1p + V2*I2p;
       func = V3*I1p + V4*I2p;
      }

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

   return fabs(func);
}
Double_t Isum_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 I1p_0=par[6];
   Double_t I2p_0=par[7];
   Double_t RV1_save=par[8];
   Double_t RV2_save=par[9];
   Double_t RV3_save=par[10];
   Double_t RV4_save=par[11];
   Double_t tlimit=par[12];
   Double_t I1plimit=par[13];
   Double_t I2plimit=par[14];
   Double_t lam1_0=par[15];
   Double_t lam2_0=par[16];

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

   char string[256];

   Double_t func;

   if (t1 < tlimit) {
       Double_t I1p = I1p_0*exp(-t1/lam1_0);
       Double_t I2p = I2p_0*exp(-t1/lam2_0);
       func = RV1_save*I1p + RV2_save*I2p + RV3_save*I1p + RV4_save*I2p;
       }
   else {
       Double_t I1p = I1plimit*exp(-(t1-tlimit)/lam1);
       Double_t I2p = I2plimit*exp(-(t1-tlimit)/lam2);
       func = V1*I1p + V2*I2p + V3*I1p + V4*I2p;
      }



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

   return func;
}