///////////////////////////////////////////////////////////////////////////////
//
// blackbody.C
//
// Plot curves of blackbody radiation.
//
// 2006/11/30 oxon
///////////////////////////////////////////////////////////////////////////////

#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>

using namespace std;

double T0 = 200;         // [K]
double dT = 20;         // [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 blackbody()
{
  TCanvas* can = new TCanvas("can", "can", 804, 628);

  TH2D* frame = new TH2D("frame", "", 1, 0, 35, 1, 0, 10);
  frame->SetTitle("Black Body Radiation;Wavelength [#mum];dE/dS/d#lambda [W/m^{2}/#mum]");
  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)", 1, 35);
    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");
  } // i

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

  return;
}

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

  x *= 1e-6; // [um]->[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/1e6; // [W/m^2/m]->[W/m^2/um]
}