#define mpts 101;
#define tpts 5;
Double_t position[mpts];
Double_t steady_data [mpts];
void steady_state2(void)
{
// 
// plot the number of fits to a cdc segment that result in prob > prob_cut (currently set to 0.01).
//

// #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 xpts 8;

   Double_t time[tpts];
   Double_t steady_time[tpts];
   Double_t x_t[xpts];
   Double_t steady_t100[xpts];    // steady state solution.
   Double_t steady_t0[xpts];
   Double_t steady_t1[xpts];
   Double_t steady_t2[xpts];
   Double_t steady_t3[xpts];
   Double_t steady_t4[xpts];

   // 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]={-10,0,2000};
   Double_t dummyy[npts]={-10,-10,-10};

   //
   TCanvas *c1 = new TCanvas("c1","c1 heat transfer",200,10,700,700);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);

   c1->SetGridx();
   c1->SetGridy();
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   c1->SetLogx();
   c1->SetLogy();


   Double_t xmin=0.001;
   Double_t xmax=1;
   // Double_t xmin=0.1;
   // Double_t xmax=0.5;
   // Double_t ymin=-50;
   // Double_t ymax=300;
   Double_t ymin=1;
   Double_t ymax=5000;

   Double_t x1=0.001;
   Double_t x2=0.003;
   Double_t x3=0.01;
   Double_t x4=0.02;
   Double_t x5=0.05;
   Double_t x6=0.10;
   Double_t x7=0.15;
   Double_t x8=0.30;

   // dummy to draw axes

   TGraph *dummy = new TGraph (npts,dummyx,dummyy);
   TLegend *leg = new TLegend(0.15,0.75,0.85,0.95);
   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("Position (m)");
   dummy->GetYaxis()->SetTitle("Temperature (K)");
   // dummy->Draw("Ap");
   dummy->GetXaxis()->SetNdivisions(505);

   // plot steady-state heat transfer equation with Gaussian source term with width sigma;

   TF1 *steadya = new TF1("steady_a",steady_a,xmin,xmax,3);
   TF1 *steadyb = new TF1("steady_b",steady_b,xmin,xmax,3);
   TF1 *steadyc = new TF1("steady_c",steady_c,xmin,xmax,5);
   TF1 *steadye = new TF1("steady_e",steady_e,xmin,xmax,4);
   Double_t xmean=0;

   Double_t Bath=4.2;    // units are K
   // Double_t Pv=2.1e7;   // units are W/m3  corresponds to 5 kW, maximum for power deposition in loop
   Double_t Pv=2.1e9;   // units are W/m3  corresponds to 60 W, maximum for power deposition in short
   Double_t Cooling=4.6e4;   // units are W/m3-K; Take 0.1 times estimate to account for film boiling Van Sciver Helium Cryogenics p. 187
   Double_t sigma0=0.0025;   // units are m, assume 5 mm short
   // Double_t sigma=0.0025;   // units are m, assume 5 mm short
   Double_t sigma=0.00035;   // units are m, assume 5 mm short
   Double_t k=250;  // units are W/m-K value at 4deg K

   Double_t T0=6;
   steadya->SetParameter(0,sigma);
   steadya->SetParameter(1,xmean);
   steadya->SetParameter(2,xmax);

   steadyb->SetParameter(0,sigma);
   steadyb->SetParameter(1,xmean);
   steadyb->SetParameter(2,xmax);

   steadyc->SetParameter(0,Cooling);
   steadyc->SetParameter(1,Pv);
   steadyc->SetParameter(2,sigma);
   steadyc->SetParameter(3,k);
   steadyc->SetParameter(4,Bath);

   for (j=0;j<mpts;j++) {
       Double_t x = xmin + j*(xmax-xmin)/(mpts-1);
       position[j] = x;
       steady_data[j] = steadyc->Eval(x);

      // printf ("x=%f, j=%d, data=%f\n",x,j,steady_data[j]);
       }

   
   TGraph *steady1 = new TGraph (mpts,position,steady_data);
   steady1->SetMarkerColor(2);
   steady1->SetMarkerStyle(20);
   steady1->SetMarkerSize(0.8);
   steady1->SetTitle("");
   steady1->GetXaxis()->SetRangeUser(xmin,xmax);
   steady1->GetYaxis()->SetRangeUser(ymin,ymax);
   steady1->GetXaxis()->SetTitleSize(0.04);
   steady1->GetYaxis()->SetTitleSize(0.04);
   steady1->GetYaxis()->SetTitleOffset(1.5);
   steady1->GetXaxis()->SetTitle("Position (m)");
   steady1->GetYaxis()->SetTitle("Temperature (K)");


   //
   TCanvas *c2 = new TCanvas("c2","c2 heat transfer",200,10,700,700);
   c2->SetBorderMode(0);
   c2->SetFillColor(0);

   c2->SetGridx();
   c2->SetGridy();
   c2->SetBorderMode(0);
   c2->SetFillColor(0);
   c2->SetLogx();
   c2->SetLogy();

   // compute time dependent solution

   Double_t tmin=0.0001;
   Double_t tmax=10;
   Double_t tdelta=10;
   // Double_t ymin=0.5;
   // Double_t ymax=100;
   Double_t ymin=0.1;
   Double_t ymax=1000;
   Double_t alpha=0.37;  // units are m2/s  value at 4 deg K.
   Double_t alpha=0.00011;  // units are m2/s  value at 100 deg K.

   Double_t x=x1;
   x_t[0]=x;
   Double_t t=tmin;
    steadye->SetParameter(0,x);
    steadye->SetParameter(2,Bath);

   for (j=0;j<tpts;j++) {
       Double_t t = t*tdelta;
       time[j] = t;

        steadye->SetParameter(1,alpha*t);
        Double_t A = steadye->Integral(xmin,xmax);
       steady_time[j] = fabs(A);    // take absolute to plot on log scale
       printf ("j=%d t=%g x=%g A=%g sigma=%g\n",j,t,x,A,sqrt(2*alpha*t));
       if (j == 0) {
          steady_t0[0] = A;
          }
      else if ( j == 1) {
          steady_t1[0] = A;
          }
      else if ( j == 2) {
          steady_t2[0] = A;
          }
      else if ( j == 3) {
          steady_t3[0] = A;
          }
      else if ( j == 4) {
          steady_t4[0] = A;
          }
      }

   TGraph *tdep = new TGraph (tpts,time,steady_time);
   TLegend *leg1 = new TLegend(0.25,0.8,0.85,0.95);

   // tdep->SetMarkerSize(0.8);
   tdep->SetTitle("");
   tdep->GetXaxis()->SetRangeUser(tmin,tmax);
   tdep->GetYaxis()->SetRangeUser(ymin,ymax);
   tdep->GetXaxis()->SetTitleSize(0.04);
   tdep->GetYaxis()->SetTitleSize(0.04);
   tdep->GetYaxis()->SetTitleOffset(1.5);
   tdep->GetXaxis()->SetTitle("time (s)");
   tdep->GetYaxis()->SetTitle("- T(t) + T(#infty) (K)");
   tdep->GetXaxis()->SetNdivisions(505);

   tdep->SetMarkerColor(2);
   tdep->SetMarkerStyle(20);
   tdep->Draw("Ap");
   sprintf (string,"x=%0.2f m, T_{b}=%0.1f K, #sigma=%g  m, #alpha=%gm2/s\n",x,Bath,sigma,alpha);
   leg1->AddEntry(tdep,string,"p");


   Double_t x=x2;  
   x_t[1]=x;
   Double_t t=tmin;
    steadye->SetParameter(0,x);
    steadye->SetParameter(2,Bath);

   for (j=0;j<tpts;j++) {
       Double_t t = t*tdelta;
       time[j] = t;

        steadye->SetParameter(1,alpha*t);
       Double_t A = steadye->Integral(xmin,xmax);
       steady_time[j] = fabs(A);    // take absolute to plot on log scale
       if (j == 0) {
          steady_t0[1] = A;
          }
      else if ( j == 1) {
          steady_t1[1] = A;
          }
      else if ( j == 2) {
          steady_t2[1] = A;
          }
      else if ( j == 3) {
          steady_t3[1] = A;
          }
      else if ( j == 4) {
          steady_t4[1] = A;
          }
       }

   TGraph *tdep1 = new TGraph (tpts,time,steady_time);
   tdep1->SetTitle("");

   tdep1->SetMarkerColor(4);
   tdep1->SetMarkerStyle(20);
   tdep1->Draw("Samep");
   sprintf (string,"x=%0.2f m, T_{b}=%0.1f K, #sigma=%g  m, #alpha=%gm2/s\n",x,Bath,sigma,alpha);
   leg1->AddEntry(tdep1,string,"p");

   Double_t x=x3;  
   x_t[2]=x;
   Double_t t=tmin;
    steadye->SetParameter(0,x);
    steadye->SetParameter(2,Bath);

   for (j=0;j<tpts;j++) {
       Double_t t = t*tdelta;
       time[j] = t;

       steadye->SetParameter(1,alpha*t);
       Double_t A = steadye->Integral(xmin,xmax);
       steady_time[j] = fabs(A);    // take absolute to plot on log scale
       if (j == 0) {
          steady_t0[2] = A;
          }
      else if ( j == 1) {
          steady_t1[2] = A;
          }
      else if ( j == 2) {
          steady_t2[2] = A;
          }
      else if ( j == 3) {
          steady_t3[2] = A;
          }
      else if ( j == 4) {
          steady_t4[2] = A;
          }
       }

   TGraph *tdep2 = new TGraph (tpts,time,steady_time);
   tdep2->SetTitle("");

   tdep2->SetMarkerColor(1);
   tdep2->SetMarkerStyle(20);
   tdep2->Draw("Samep");
   sprintf (string,"x=%0.2f m, T_{b}=%0.1f K, #sigma=%g  m, #alpha=%gm2/s\n",x,Bath,sigma,alpha);
   leg1->AddEntry(tdep2,string,"p");

   Double_t x=x4;  
   x_t[3]=x;
   Double_t t=tmin;
    steadye->SetParameter(0,x);
    steadye->SetParameter(2,Bath);

   for (j=0;j<tpts;j++) {
       Double_t t = t*tdelta;
       time[j] = t;

       steadye->SetParameter(1,alpha*t);
       Double_t A = steadye->Integral(xmin,xmax);
       printf ("j=%d, x=%g, t=%g, A=%g\n",j,x,t,A);
       steady_time[j] = fabs(A);    // take absolute to plot on log scale
       if (j == 0) {
          steady_t0[3] = A;
          }
      else if ( j == 1) {
          steady_t1[3] = A;
          }
      else if ( j == 2) {
          steady_t2[3] = A;
          }
      else if ( j == 3) {
          steady_t3[3] = A;
          }
      else if ( j == 4) {
          steady_t4[3] = A;
          }
       }

   TGraph *tdep3 = new TGraph (tpts,time,steady_time);
   tdep3->SetTitle("");

   tdep3->SetMarkerColor(3);
   tdep3->SetMarkerStyle(20);
   tdep3->Draw("Samep");
   sprintf (string,"x=%0.2f m, T_{b}=%0.1f K, #sigma=%g  m, #alpha=%gm2/s\n",x,Bath,sigma,alpha);
   leg1->AddEntry(tdep3,string,"p");


   Double_t x=x5;  
   x_t[4]=x;
   Double_t t=tmin;
    steadye->SetParameter(0,x);
    steadye->SetParameter(2,Bath);

   for (j=0;j<tpts;j++) {
       Double_t t = t*tdelta;
       time[j] = t;

       steadye->SetParameter(1,alpha*t);
       Double_t A = steadye->Integral(xmin,xmax);
       printf ("j=%d, x=%g, t=%g, A=%g\n",j,x,t,A);
       steady_time[j] = fabs(A);    // take absolute to plot on log scale
       if (j == 0) {
          steady_t0[4] = A;
          }
      else if ( j == 1) {
          steady_t1[4] = A;
          }
      else if ( j == 2) {
          steady_t2[4] = A;
          }
      else if ( j == 3) {
          steady_t3[4] = A;
          }
      else if ( j == 4) {
          steady_t4[4] = A;
          }
       }

   TGraph *tdep4 = new TGraph (tpts,time,steady_time);
   tdep4->SetTitle("");

   tdep4->SetMarkerColor(3);
   tdep4->SetMarkerStyle(20);
   tdep4->Draw("Samep");
   sprintf (string,"x=%0.2f m, T_{b}=%0.1f K, #sigma=%g  m, #alpha=%gm2/s\n",x,Bath,sigma,alpha);
   leg1->AddEntry(tdep4,string,"p");

   Double_t x=x6;  
   x_t[5]=x;
   Double_t t=tmin;
    steadye->SetParameter(0,x);
    steadye->SetParameter(2,Bath);

   for (j=0;j<tpts;j++) {
       Double_t t = t*tdelta;
       time[j] = t;

       steadye->SetParameter(1,alpha*t);
       Double_t A = steadye->Integral(xmin,xmax);
       printf ("j=%d, x=%g, t=%g, A=%g\n",j,x,t,A);
       steady_time[j] = fabs(A);    // take absolute to plot on log scale
       if (j == 0) {
          steady_t0[5] = A;
          }
      else if ( j == 1) {
          steady_t1[5] = A;
          }
      else if ( j == 2) {
          steady_t2[5] = A;
          }
      else if ( j == 3) {
          steady_t3[5] = A;
          }       
      else if ( j == 4) {
          steady_t4[5] = A;
          }
       }

   TGraph *tdep5 = new TGraph (tpts,time,steady_time);
   tdep5->SetTitle("");

   tdep5->SetMarkerColor(3);
   tdep5->SetMarkerStyle(20);
   tdep5->Draw("Samep");
   sprintf (string,"x=%0.2f m, T_{b}=%0.1f K, #sigma=%g  m, #alpha=%gm2/s\n",x,Bath,sigma,alpha);
   leg1->AddEntry(tdep5,string,"p");

   Double_t x=x7;  
   x_t[6]=x;
   Double_t t=tmin;
    steadye->SetParameter(0,x);
    steadye->SetParameter(2,Bath);

   for (j=0;j<tpts;j++) {
       Double_t t = t*tdelta;
       time[j] = t;

       steadye->SetParameter(1,alpha*t);
       Double_t A = steadye->Integral(xmin,xmax);
       printf ("j=%d, x=%g, t=%g, A=%g\n",j,x,t,A);
       steady_time[j] = fabs(A);    // take absolute to plot on log scale
       if (j == 0) {
          steady_t0[6] = A;
          }
      else if ( j == 1) {
          steady_t1[6] = A;
          }
      else if ( j == 2) {
          steady_t2[6] = A;
          }
      else if ( j == 3) {
          steady_t3[6] = A;
          }
      else if ( j == 4) {
          steady_t4[6] = A;
          }
       }

   TGraph *tdep6 = new TGraph (tpts,time,steady_time);
   tdep6->SetTitle("");

   tdep6->SetMarkerColor(3);
   tdep6->SetMarkerStyle(20);
   tdep6->Draw("Samep");
   sprintf (string,"x=%0.2f m, T_{b}=%0.1f K, #sigma=%g  m, #alpha=%gm2/s\n",x,Bath,sigma,alpha);
   leg1->AddEntry(tdep6,string,"p");

   Double_t x=x8;  
   x_t[7]=x;
   Double_t t=tmin;
    steadye->SetParameter(0,x);
    steadye->SetParameter(2,Bath);

   for (j=0;j<tpts;j++) {
       Double_t t = t*tdelta;
       time[j] = t;

       steadye->SetParameter(1,alpha*t);
       Double_t A = steadye->Integral(xmin,xmax);
       printf ("j=%d, x=%g, t=%g, A=%g\n",j,x,t,A);
       steady_time[j] = fabs(A);    // take absolute to plot on log scale
       if (j == 0) {
          steady_t0[7] = A;
          }
      else if ( j == 1) {
          steady_t1[7] = A;
          }
      else if ( j == 2) {
          steady_t2[7] = A;
          }
      else if ( j == 3) {
          steady_t3[7] = A;
          }
      else if ( j == 4) {
          steady_t4[7] = A;
          }
       }

   TGraph *tdep7 = new TGraph (tpts,time,steady_time);
   tdep7->SetTitle("");

   tdep7->SetMarkerColor(3);
   tdep7->SetMarkerStyle(20);
   tdep7->Draw("Samep");
   sprintf (string,"x=%0.2f m, T_{b}=%0.1f K, #sigma=%g  m, #alpha=%gm2/s\n",x,Bath,sigma,alpha);
   leg1->AddEntry(tdep7,string,"p");

   leg1->Draw();


   // sprintf(string,"#sigma =0.3f m^{2}/s\n",sigma);
   t1 = new TLatex(0.50,0.65,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   // t1->Draw(); 

   // compute asymptotic solution for future use

   for (j=0;j<xpts;j++) {
       Double_t x = x_t[j];
       steady_t100[j] = steadyc->Eval(x);
       }

   // go back to first canvas and plot x dependence

   c1->cd();

   steady1->Draw("Ap");
   steady1->GetXaxis()->SetNdivisions(505);
   sprintf (string,"Steady State T_{b}=%0.1f K, #sigma=%g  m\n",Bath,sigma);
   // leg->AddEntry(steady1,string,"p");

   for (j=0;j<xpts;j++) {
        printf ("\nj=%d x=%g, 0=%f, 1=%f, 2=%f, 3=%f, 4=%f 100=%f\n",j,x_t[j],steady_t0[j],steady_t1[j],steady_t2[j],steady_t3[j],steady_t4[j],steady_t100[j]);
        steady_t0[j] = steady_t0[j] + steady_t100[j];
        steady_t1[j] = steady_t1[j] + steady_t100[j];
        steady_t2[j] = steady_t2[j] + steady_t100[j];
        steady_t3[j] = steady_t3[j] + steady_t100[j];
        steady_t4[j] = steady_t4[j] + steady_t100[j];

        // cheat and require answer to be physical. Integration must have numerical inconsistencies

        if (steady_t0[j] < Bath) steady_t0[j]  = Bath;
        if (steady_t1[j] < Bath) steady_t1[j]  = Bath;
        if (steady_t2[j] < Bath) steady_t2[j]  = Bath;
        if (steady_t3[j] < Bath) steady_t3[j]  = Bath;
        if (steady_t4[j] < Bath) steady_t4[j]  = Bath;
        printf ("\nj=%d x=%g, 0=%f, 1=%f, 2=%f, 3=%f, 4=%f 100=%f\n",j,x_t[j],steady_t0[j],steady_t1[j],steady_t2[j],steady_t3[j],steady_t4[j],steady_t100[j]);
        }

   TGraph *steadyt0 = new TGraph (xpts,x_t,steady_t0);
   steadyt0->SetTitle("");
   steadyt0->SetMarkerColor(4);
   steadyt0->SetMarkerStyle(20);
   steadyt0->Draw("Samep");
   sprintf (string,"t=%0.4f s, T_{b}=%0.1f K, #sigma=%g  m, #alpha=%gm2/s\n",time[0],Bath,sigma,alpha);
   leg->AddEntry(steadyt0,string,"p");

   TGraph *steadyt1 = new TGraph (xpts,x_t,steady_t1);
   steadyt1->SetTitle("");
   steadyt1->SetMarkerColor(1);
   steadyt1->SetMarkerStyle(20);
   steadyt1->Draw("Samep");
   sprintf (string,"t=%0.4f s, T_{b}=%0.1f K, #sigma=%g  m, #alpha=%gm2/s\n",time[1],Bath,sigma,alpha);
   leg->AddEntry(steadyt1,string,"p");

   TGraph *steadyt2 = new TGraph (xpts,x_t,steady_t2);
   steadyt2->SetTitle("");
   steadyt2->SetMarkerColor(3);
   steadyt2->SetMarkerStyle(20);
   steadyt2->Draw("Samep");
   sprintf (string,"t=%0.4f s, T_{b}=%0.1f K, #sigma=%g  m, #alpha=%gm2/s\n",time[2],Bath,sigma,alpha);
   leg->AddEntry(steadyt2,string,"p");

   TGraph *steadyt3 = new TGraph (xpts,x_t,steady_t3);
   steadyt3->SetTitle("");
   steadyt3->SetMarkerColor(2);
   steadyt3->SetMarkerStyle(24);
   steadyt3->Draw("Samep");
   sprintf (string,"t=%0.4f s, T_{b}=%0.1f K, #sigma=%g  m, #alpha=%gm2/s\n",time[3],Bath,sigma,alpha);
   leg->AddEntry(steadyt3,string,"p");

   TGraph *steadyt4 = new TGraph (xpts,x_t,steady_t4);
   steadyt4->SetTitle("");
   steadyt4->SetMarkerColor(2);
   steadyt4->SetMarkerStyle(25);
   steadyt4->Draw("Samep");
   sprintf (string,"t=%0.4f s, T_{b}=%0.1f K, #sigma=%g  m, #alpha=%gm2/s\n",time[4],Bath,sigma,alpha);
   leg->AddEntry(steadyt4,string,"p");

   TGraph *steadyt100 = new TGraph (xpts,x_t,steady_t100);
   steadyt100->SetTitle("");
   steadyt100->SetMarkerColor(2);
   steadyt100->SetMarkerStyle(20);
   steadyt100->Draw("Samep");
   sprintf (string,"Steady State T_{b}=%0.1f K, #sigma=%g  m, #alpha=%gm2/s\n",Bath,sigma,alpha);
   leg->AddEntry(steadyt100,string,"p");

   leg->Draw();
  
   sprintf(filename,"steady_state2_c1.eps");
   c1->SaveAs(filename);
   sprintf(filename,"steady_state2_c1.png");
   c1->SaveAs(filename);
   sprintf(filename,"steady_state2_c2.eps");
   c2->SaveAs(filename);
   sprintf(filename,"steady_state2_c2.png");
   c2->SaveAs(filename);

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

   // assume that initial conditions consist of u(0) = U0, u'(0)=0
   // Pulse response function

   Double_t sigma=par[0];
   Double_t xmean=par[1];
   Double_t xtop=par[2];
   Double_t pi=3.14159;
   Double_t xmin=-10;
   Double_t xmax=10;
   Double_t x1=(x[0]-xmean)/sigma;

   char string[256];
 
   if (x1 > xmax) x1 = xmax;
   if (x1 < xmin) x1 = xmin;
   Double_t u = 1;

   sprintf (string,"xmean=%f, sigma=%f, xtop=%f, x=%f, u=%f\n",xmean,sigma,xtop,x1,u);
   printf ("string=%s",string);

   return u;
}
Double_t steady_b (Double_t *x, Double_t *par) {

   // assume that initial conditions consist of u(0) = TU, u'(0)=0
   // Pulse response function

   Double_t sigma=par[0];
   Double_t xmean=par[1];
   Double_t xtop=par[2];
   Double_t pi=3.14159;
   Double_t xmin=-10;
   Double_t xmax=10;
   Double_t x1=(x[0]-xmean)/sigma;

   char string[256];
 
   if (x1 > xmax) x1 = xmax;
   if (x1 < xmin) x1 = xmin;
   Double_t u = 1;

   /*sprintf (string,"xmean=%f, sigma=%f, x=%f, u=%f\n",xmean,sigma,x1,u);
   printf ("string=%s",string);*/

   return u;
}
Double_t steady_c (Double_t *x, Double_t *par) {

   // Compute steady state solution to heat transfer equation with constant source Pv on [-sigma,sigma]
   // in a cooling bath at temperature Bath

   Double_t Cooling=par[0];
   Double_t Pv=par[1];
   Double_t sigma=par[2];
   Double_t k=par[3];
   Double_t Bath=par[4];
   Double_t x1=x[0];
   Double_t pi=3.14159;

   char string[256];

   Double_t func;
   Double_t omega = sqrt(Cooling/k);
   Double_t PvC = Pv/Cooling;
   if (fabs(x1) < sigma) {
       func = (PvC + Bath) - PvC*cosh(omega*x1)*(cosh(omega*sigma)-sinh(omega*sigma));
       }
   else {
       Double_t xomega= omega*sigma;
       Double_t x1omega = omega*fabs(x1);
       if (xomega <20 && x1omega<20) {
           func = Bath + PvC*sinh(xomega)*(cosh(x1omega)-sinh(x1omega));
          }
       else  {
          // protect against product of large*small numbers 
          func = Bath + PvC*0.5*exp(xomega-x1omega);
       }
    }

   /*sprintf (string,"x1=%f omega=%g, PvC=%g, Bath=%g, func=%f\n",x1,omega,PvC,Bath,func);
   printf ("string=%s",string);*/

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

   // assume that initial conditions consist of u(0) = U0, u'(0)=0
   // Function called AFTER steady_data has been filled using steady_c and steady_d.

   Double_t xvar=par[0];
   Double_t alpha_t=par[1];
   Double_t Bath=par[2];
   Double_t norm=par[3];
   Double_t pi=3.14159;
   Double_t x1=x[0];
   Int_t j;

   Double_t xvarmin = position[0];
   Double_t xvarmax = position[mpts-1];
    j = (x1 - xvarmin)*(mpts-1)/(xvarmax-xvarmin);
   if (j < 0) j = 0;
   if (j >= mpts) j = mpts-1;
   Double_t deltaj = (x1 - xvarmin)*(mpts-1)/(xvarmax-xvarmin) - j;
   Double_t data = steady_data[j] + deltaj *( steady_data[j+1] - steady_data[j]);
   // printf ("xvar=%f, x1=%f, xvarmin=%f, xvarmax=%f, j=%d, deltaj=%f, data=%f\n",xvar,x1,xvarmin,xvarmax,j,deltaj,data);

   // compute logs to maintain precision

   char string[256];
   Double_t logpart1;
   Double_t logpart1min=-100;
   if (alpha_t > 0) {
      logpart1 = -(xvar-x1)*(xvar-x1)/(4*alpha_t) -0.5*log(4*pi*alpha_t);
      if (logpart1 < logpart1min) logpart1 = logpart1min;
      }
   else {
       logpart1 = logpart1min;
       }

   // normalize intergrals to compensate for finite lengths


    Double_t part2 = Bath -  data;
   
   /*if (norm == 0) { 
       Double_t part2 =1;
       }
   else if {
       Double_t part2 = Bath -  data;
        }*/

   Double_t u = exp(logpart1)*part2;

  /* if (xvar == 0) {
       sprintf (string,"x1=%g xvar=%g, data=%g, u=%f\n",x1,xvar,data,u);
       printf ("string=%s",string);
       }*/


   /*for (j=0;j<mpts;j++) {
        printf("j=%d x=%f, steady=%f\n",j,position[j],steady_data[j]);
        }*/

   return u;
}
Double_t g2_func (Double_t *x, Double_t *par) {

   //

   Double_t sigma=par[0];
   Double_t pi=3.14159;
   Double_t xmean=par[1];
   Double_t xmin=-5;
   Double_t xmax=5;
   Double_t x1=(x[0]-xmean)/sigma;

   char string[256];
 
   Double_t u = exp(-x1*x1/2)/(sqrt(2*pi)*sigma);
   
   /*sprintf (string,"x=%f, u=%f\n",x1,u);
   printf ("string=%s",string);*/

   return u;
}

Double_t steady_test (Double_t *x, Double_t *par) {

   // assume that initial conditions consist of u(0) = T0, u'(0)=0

   Double_t sigma=par[0];
   Double_t T0=par[1];
   Double_t pi=3.14159;
   Double_t xmean=1;
   Double_t xmin=-5;
   Double_t xmax=5;
   Double_t x1=(x[0]-xmean)/sigma;

   char string[256];
 
   Double_t u;
   TF1 *g2;
   // if (g2 != NULL) TF1 *g2 = new TF1("g2","TMath::Gaus(x)",xmin,xmax);
   // gint->SetParameter(0,xmean);
   // gint->SetParameter(1,sigma);
   // Double_t u = g2->Eval(x1)/(sqrt(2*pi)*sigma);
   if (x1 < xmin) {
       u = 0;
       }
   else if (x1 < xmax) {
       u = g2->Integral(xmin,x1)/(sqrt(2*pi)*sigma);
       }
   else {
       u = g2->Integral(xmin,xmax)/(sqrt(2*pi)*sigma);
   }
   
   /*sprintf (string,"x=%f, u=%f\n",x1,u);
   printf ("string=%s",string);*/

   return u;
}