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

#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 coil2M_dump",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 xmax=2;
   //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=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;
   Double_t Rl_0 = 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

   Double_t gamma = gamma0;

   // left side equation

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

   // Rl = Rl_super;
   printf ("L1=%g, L2=%g, M=%g, gamma0=%g, sigma=%g, area=%g, Rd=%g, Rs=%g Rl=%g\n",L1,L2,M,gamma0,sigma,area,Rd,Rs,Rl);


   // fill matrices

      Double_t det;
      Double_t col1[2];
      Double_t col2[2];


   // Loop over time and compute currents and other quantities.

   Double_t t;
   Double_t tlimit=xmax;
   Double_t time[jmax];
   Double_t curr1[jmax];
   Double_t curr2[jmax];
   Double_t abscurr2[jmax];
   Double_t currsum[jmax];

   Double_t Power1[jmax];
   Double_t Power2[jmax];
   Double_t Power3[jmax];
   Double_t Energy1=0;
   Double_t Energy2=0;
   Double_t Energy3=0;
   Double_t Voltage1[jmax];
   Double_t Voltage2[jmax];
   Double_t Voltage3[jmax];
   Double_t Voltage4[jmax];

      // initial values
      Double_t ctemp[2];
      ctemp[0] = I1_0;
      ctemp[1] = I2_0;
      TVectorD Curr; Curr.Use(2,ctemp);
        curr1[0] = I1_0;
        curr2[0] = I2_0;
        currsum[0] = I1_0 + I2_0;
        Rl = Rl_super;
        Int_t once = 1;

   Double_t deltat = (xmax*xlogfactor-xmin)/jmax;
   for (j=0;j<jmax;j++) {
        t = xmin + j*(xmax*xlogfactor-xmin)/jmax;

      col1[0] = L1;
      col1[1] = gamma*M;
      col2[0] = gamma*M;
      col2[1] = L2;
      TVectorD x1; x1.Use(2,col1);
      TVectorD x2; x2.Use(2,col2);
      TMatrixD MatrixM(2,2);
      TMatrixDColumn(MatrixM,0) = x1;
      TMatrixDColumn(MatrixM,1) = x2;
      // MatrixM.Print();
      TMatrixD MatrixMinv = MatrixM;        // copy first to avoid overwriting M
      MatrixMinv.InvertFast(&det);
      // MatrixMinv.Print();


      col1[0] =Rd+Rl ;
      col1[1] = Rl;
      col2[0] = Rl;
      col2[1] = Rs+Rl;
      TVectorD x1; x1.Use(2,col1);
      TVectorD x2; x2.Use(2,col2);
      TMatrixD MatrixB(2,2);
      TMatrixDColumn(MatrixB,0) = x1;
      TMatrixDColumn(MatrixB,1) = x2;
      // MatrixB.Print();

       TMatrixD MatrixMinvB = MatrixMinv * MatrixB;

      TVectorD eigenvalues(2);
      TMatrixD MatrixU = MatrixMinvB.EigenVectors(eigenvalues);
      //eigenvalues.Print();
      TMatrixD MatrixUinv = MatrixU;    // copy first to avoid overwriting U
      MatrixUinv.InvertFast(&det);
      TMatrixD MatrixMinvBd = MatrixUinv * MatrixMinvB * MatrixU;


   // compute currents, power dissipated  and voltages

        time[j] = t;
        curr1[j] = Curr(0);
        curr2[j] = Curr(1);
        abscurr2[j] = fabs(curr2[j]);
        // printf (" time = %f curr1 = %f curr2 = %f\n",time[j],curr1[j],curr2[j]);
        Power1[j] = curr1[j]*curr1[j]*Rd;
        Voltage1[j] = curr1[j]*Rd;
        Energy1 = Energy1 + Power1[j]*deltat;
        Power2[j] = curr2[j]*curr2[j]*Rs;
        Voltage2[j] = curr2[j]*Rs;
        Energy2 = Energy2 + Power2[j]*deltat;
        currsum[j] = curr1[j] + curr2[j];
        if (once) {
               Power3[j] = 0;
            }
           else {
               Power3[j] = currsum[j]*currsum[j]*Rl;
           }
        Voltage3[j] = Voltage1[j] - Voltage2[j];
        Voltage4[j] = Voltage1[j] - Voltage3[j];
        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]);

       // copy currents for next step and repeat

      TVectorD CurrU  = MatrixUinv * Curr;
       // printf ("1  time = %f currU1 = %f currU2 = %f eigen = %f %f \n",time[j],CurrU(0),CurrU(1),eigenvalues(0),eigenvalues(1));
      // CurrU.Print();
      Double_t arg = -deltat*eigenvalues(0);
      if (arg > 100)  arg = 100;
      if (arg < -100) arg = -100;
      CurrU(0) = CurrU(0)* exp(arg) ;
      Double_t arg = -deltat*eigenvalues(1);
      if (arg > 100)  arg = 100;
      if (arg < -100) arg = -100;
      CurrU(1) = CurrU(1) * exp(arg);
      // eigenvalues.Print();
      // Curr.Print();
      // CurrU.Print();

      TVectorD Curr  = MatrixU * CurrU;

     // update value of resistance in loop

      if (currsum[j] < 4000 && once) {
           Rl = Rl_super;
           // printf ("1 current sum %f < 4000\n",currsum[j]);
            }
      else {
           Rl = Rl_0;
           once = 0;
           } 

        }

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


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

   TGraph *I1 = new TGraph (jmax,time,curr1);
   TLegend *leg = new TLegend(0.50,0.75,0.85,0.9);
   I1->SetMarkerColor(1);
   I1->SetMarkerStyle(21);
   I1->SetMarkerSize(0.35);
   leg->AddEntry(I1,"Solenoid Current","p");
   I1->Draw("Samep");

   TGraph *I2 = new TGraph (jmax,time,curr2);
   I2->SetMarkerColor(4);
   I2->SetMarkerStyle(21);
   I2->SetMarkerSize(0.35);
   leg->AddEntry(I2,"Short Current","p");
   I2->Draw("Samep");

   TGraph *Isum = new TGraph (jmax,time,currsum);
   Isum->SetMarkerColor(2);
   Isum->SetMarkerStyle(21);
   Isum->SetMarkerSize(0.35);
   leg->AddEntry(Isum,"Current Sum","p");
   Isum->Draw("Samep");

   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",1/eigenvalues(0));
   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",1/eigenvalues(1));
   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);
   TLegend *leg2 = new TLegend(0.50,0.75,0.85,0.9);
   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");

   I1->Draw("Samep");
   leg2->AddEntry(I1,"Solenoid Current","p");

   TGraph *aI2 = new TGraph (jmax,time,abscurr2);
   aI2->SetMarkerColor(4);
   aI2->SetMarkerStyle(21);
   aI2->SetMarkerSize(0.35);
   leg2->AddEntry(aI2,"Short Current","p");
   aI2->Draw("Samep");

   Isum->Draw("Samep");
   leg2->AddEntry(Isum,"Short Current","p");
   leg2->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");


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

   Double_t ymin=0.01;
   Double_t ymax=200000;
   TGraph *VV1 = new TGraph (jmax,time,Voltage1);
   TLegend *leg = new TLegend(0.50,0.75,0.85,0.9);
   VV1->SetMarkerColor(1);
   VV1->SetMarkerStyle(21);
   VV1->SetMarkerSize(0.35);
   VV1->SetTitle("");
   VV1->GetXaxis()->SetRangeUser(xmin,xmax*xlogfactor);
   ymin=89.5;
   ymax=90.;
   VV1->GetYaxis()->SetRangeUser(ymin,ymax);
   VV1->GetXaxis()->SetTitleSize(0.04);
   VV1->GetYaxis()->SetTitleSize(0.04);
   VV1->GetYaxis()->SetTitleOffset(1.5);
   VV1->GetXaxis()->SetTitle("Time (s)");
   VV1->GetYaxis()->SetTitle("Voltage (V)");
   VV1->GetXaxis()->SetNdivisions(505);
   VV1->Draw("Ap");
   leg->AddEntry(VV1,"Total V","p");
   leg->Draw();

   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 *VV2 = new TGraph (jmax,time,Voltage2);
   // TLegend *leg = new TLegend(0.50,0.75,0.85,0.9);
   VV2->SetMarkerColor(4);
   VV2->SetMarkerStyle(21);
   VV2->SetMarkerSize(0.35);
   ymin= -0.005;
   ymax= 0.005;
   VV2->GetYaxis()->SetRangeUser(ymin,ymax);
   VV2->GetXaxis()->SetTitleSize(0.04);
   VV2->GetYaxis()->SetTitleSize(0.04);
   VV2->GetYaxis()->SetTitleOffset(1.5);
   VV2->GetXaxis()->SetTitle("Time (s)");
   VV2->GetYaxis()->SetTitle("Voltage (V)");
   VV2->GetXaxis()->SetNdivisions(505);
   VV2->Draw("Samep");
   leg->AddEntry(VV2,"Short V","p");
   leg->Draw();

   TGraph *VV3 = new TGraph (jmax,time,Voltage3);
   // TLegend *leg = new TLegend(0.50,0.75,0.85,0.9);
   VV3->SetMarkerColor(2);
   VV3->SetMarkerStyle(21);
   VV3->SetMarkerSize(0.35);
   VV3->Draw("Samep");
   leg->AddEntry(VV3,"V coil","p");
   leg->Draw();

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


   //
   TCanvas *c2 = new TCanvas("c2","c2 coil2M_dump",200,10,700,700);

   c2->SetGridx();
   c2->SetGridy();
   c2->SetBorderMode(0);
   c2->SetFillColor(0);
   c2->Divide(2,2);


   c2->cd(1);
   c2_1->SetGridx();
   c2_1->SetGridy();
   c2_1->SetLogx();
  VV1->Draw("Ap");
   leg->Draw();

   c2->cd(2);
   c2_2->SetGridx();
   c2_2->SetGridy();
   c2_2->SetLogx();
   VV2->Draw("Ap");
   leg->Draw();

   c2->cd(3);
   c2_3->SetGridx();
   c2_3->SetGridy();
   c2_3->SetLogx();
   VV1->Draw("Ap");
   VV3->Draw("Samep");
   leg->Draw();

   c2->cd(4);
   /*c2_4->SetLogx();
   leg->Draw();
   VV4->Draw("Ap");*/

   // save canvas to file

   Int_t ig = gamma*1000;
   Int_t ir = rf*10;
   sprintf(filename,"coil2M_dump_c1_%d_%d.pdf",ig,ir);
   c1->SaveAs(filename);
   sprintf(filename,"coil2M_dump_c1_%d_%d.png",ig,ir);
   c1->SaveAs(filename);
   sprintf(filename,"coil2M_dump_c2_%d_%d.pdf",ig,ir);
   c2->SaveAs(filename);
   sprintf(filename,"coil2M_dump_c2_%d_%d.png",ig,ir);
   c2->SaveAs(filename);


}