#!/usr/bin/perl
#Varun Bhalerao
#http://www.astro.caltech.edu/~varun/ay20/ps2/blackbodysun.pl

#cgs units !

$h=6.63e-27;
$c=2.998e10;
$k=1.38e-16;
$hc2=$h*$c*$c*2; #per steradian is wanted, so no pi !

#print "Enter the temperature: ";
$t=5777;
#$t=shift;

$hckt=$h*$c/$k/$t;
$lmax=0.3/$t;
$factor=100;		#factor above and below lambda_max in which to integrate
$step=100;		#step size as a fraction of wavelength, 1/$step = dlambda / lambda
#$lstart=$lmax/$factor;
#$lend=$lmax*$factor;
$lstart=505e-7;
$lend=595e-7;
print "Calculating from $lstart to $lend \n";
$dlambda=$lstart/$step;

$lambda=$lstart;
while ($lambda<$lend){
	$bldl=$hc2/$lambda**5/ (exp($hckt/$lambda) -1) *$dlambda;
	$integral+=$bldl;
	$dlambda=$lambda/$step;
	$lambda+=$dlambda;
	#print "$lambda $bldl\n";
}

$sigma=5.67e-5;
$exp=$sigma*$t**4/3.1416; #per steradian is wanted, so divide by pi !
print "$t $integral\n";
print "expected $exp\n";
printf("%E\n",$integral);
printf("%E\n",$exp);