///////////////////////////////////////////////////////////////////////////////
//
// blackbody.C
//
// Plot curves of blackbody radiation.
//
// 2006/11/30 oxon
// 2010/2/12 habig - modifed for visual
///////////////////////////////////////////////////////////////////////////////

#include "TF1.h";
#include "TH1D.h";
#include "TCanvas.h";
#include "TLatex.h";
#include "TLegend.h";
#include "TBox.h";
#include "TH2D.h";
#include "TColor.h";
#include "TROOT.h";
#include <fstream>
//#include "Math/GSLIntegrator.h"
//#include "Math/WrappedTF1.h"


using namespace std;

double T0 = 2600;         // [K]
double dT = 200;         // [K]
const int kN = 5;
double h = 6.626068e-34; // [m^2*kg/s]
double c = 2.99792458e8; // [m/s]
double k = 1.3806503e-23;// [m^2*kg/s^2/K]
double e = 1.6e-19;      // [C]

double T;

double Planck(double x);

void bb()
{
  TCanvas* can = new TCanvas("can", "can", 804, 628);

  TH2D* frame = new TH2D("frame", "", 1, 100, 5000, 1, 0, 2000);
  frame->SetTitle("Black Body Radiation;Wavelength [nm];dE/dS/d#lambda [W/m^{2}/nm]");
  frame->Draw();

  //  TText* txt = new TText(5.5, 3.8, "Sensitive Region of IR Camera");
  //  txt->SetTextSize(0.035);
  //  txt->SetTextColor(16);
  // txt->SetTextAngle(90);
  // txt->SetTextAlign(12);
  // txt->Draw();

  TBox* box = new TBox(8, 0, 12, 10);
  box->SetFillColor(16);
  box->SetLineColor(0);
  box->Draw();

  TLegend* leg = new TLegend(0.6, 0.6, 0.85, 0.85);

  TF1* planck[kN];
  for(int i=kN-1; i>=0; i--){
    T = T0 + i*dT;
    planck[i] = new TF1(Form("planck%d", i), "Planck(x)", 100, 5000);
    planck[i]->SetNpx(1200);
    planck[i]->SetLineWidth(2);
    TColor color(10000+i, (double)i/(kN-1), 0, 1-(double)i/(kN-1));
    planck[i]->SetLineColor(10000+i);
    planck[i]->GetHistogram()->DrawClone("same");
    leg->AddEntry(planck[i], Form("T = %d (K)", (int)T), "l");
    if (i == 2) {
      // Integrate the 3000k version from 400 to 700 nm
      double result;
      //     ROOT::Math::GSLIntegrator* integral = new GSLIntegrator(planck[i],&result;400.0,700.0);
      result = planck[i]->Integral(400.0,700.0);
      cout << result << " W/m^2 in visible light" << endl;
      double lambda_max = (2.898e-3/3000.0)*1e9;
      cout << "lambda_max (nm) is " << lambda_max << endl;
      cout << "ratio at 400 nm: " << planck[i]->Eval(400.0)/ planck[i]->Eval(lambda_max) << endl;
      cout << "ratio at 700 nm: " << planck[i]->Eval(700.0)/ planck[i]->Eval(lambda_max) << endl;
    }
  } // i

  leg->SetLineColor(0);
  leg->SetFillColor(0);
  leg->Draw();

  return;
}

//_____________________________________________________________________________
double Planck(double x /*[nm]*/)
{
  // Planck Distribution:
  // B(lamda) = 2hc^2/lambda^5/(exp(hc/(lamba*kT)) - 1)

  x *= 1e-9; // [nm]->[m]

  double tmp = h*c/(x*k*T); // = h\nu/kT
  double B= 2*h*pow(c, 2)/pow(x, 5)/(exp(tmp)-1);
  return B/1e9; // [W/m^2/m]->[W/m^2/nm]
}