void steady_state3(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]={-1,0,1};
   Double_t dummyy[npts]={-1,-1,-1};

   //
   TCanvas *c1 = new TCanvas("c1","c1 steady_state3",200,10,700,700);

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


   Double_t xmin=0.001;
   Double_t xmax=10;
   // Double_t xmin=0.1;
   // Double_t xmax=0.5;
   Double_t ymin=1;
   Double_t ymax=5000;
   Double_t Bath=4.2;    // units are K
   // Double_t Pv0=2.1e7;   // units are W/m3  corresponds to 5 kW, maximum for power deposition in loop
   Double_t Pv0=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

   // 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("Position (m)");
   dummy->GetYaxis()->SetTitle("Temperature (K)");
   dummy->Draw("Ap");
   dummy->GetXaxis()->SetNdivisions(505);

   // plot temperature function. Consider x as a parameter and plot time dependence.

   TF1 *steady = new TF1("steady_func",steady_func,xmin,xmax,5);
   steady->SetParameter(0,Cooling);
   Double_t Pv = Pv0;
   steady->SetParameter(1,Pv);
   steady->SetParameter(2,sigma);
   steady->SetParameter(3,k);
   steady->SetParameter(4,Bath);

   sprintf (string,"k=%g W/m-K\n",k);
   TF1 *steady1 = steady->DrawCopy("SameC");
   leg->AddEntry(steady1,string,"l");
   steady1->SetLineColor(2);

   // Pv = 2*Pv0;
   // steady->SetParameter(1,Pv);
   k = 500;
   steady->SetParameter(3,k);
   sprintf (string,"k=%g W/m-K\n",k);
   TF1 *steady2 = steady->DrawCopy("SameC");
   leg->AddEntry(steady2,string,"l");
   steady2->SetLineColor(4);

   // Pv = 0.5*Pv0;
   // steady->SetParameter(1,Pv);
   k = 1200;
   steady->SetParameter(3,k);
   sprintf (string,"k=%g W/m-K\n",k);
   TF1 *steady3 = steady->DrawCopy("SameC");
   leg->AddEntry(steady3,string,"l");
   steady3->SetLineColor(3);

   leg->Draw();

   sprintf (string,"Cooling=%g W/m^{3}-K\n",Cooling);
   t1 = new TLatex(0.20,0.85,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf (string,"#sigma=%g m\n",sigma);
   t1 = new TLatex(0.20,0.81,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf (string,"Pv=%g W/m3\n",Pv0);
   t1 = new TLatex(0.20,0.77,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf (string,"Bath=%g K\n",Bath);
   t1 = new TLatex(0.20,0.73,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
  
   sprintf(filename,"steady_state3_c1.eps");
   c1->SaveAs(filename);
   sprintf(filename,"steady_state3_c1.png");
   c1->SaveAs(filename);

}
Double_t steady_func (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;
}