#!/usr/bin/perl
# Varun Bhalerao
# Calculate the energy given as kinetic energy by ionizong photons
# cgs units !
# Ay 126 problem set 1.
# refer to the problem set at http://www.astro.caltech.edu/~varun/ay126/hw1.pdf
# Varun Bhalerao

# usage: ./q-energy.pl temperature

$h=6.63e-27;
$c=2.998e10;
$k=1.38e-16;
$hc2=$h*$c*$c*2*3.14159; 

($t)=@ARGV;

if (not $t){
	print "Usage: ./q-energy.pl temperature\n";
	exit;
}

$hckt=$h*$c/$k/$t;
$lstart = $hckt/10000.0;     # calculate from ~1000 times lower than peak lambda
$lend   = 9.14e-6;              # upper limit is the ionization potential
$e0     = $h*$c/$lend;          # ionization potential
$step=100;		# step size as a fraction of wavelength, 1/$step = dlambda / lambda

$lambda=$lstart;
while ($lambda<$lend){
	$dlambda=$lambda/$step;
	$bldl=2*3.1416*$c/$lambda**4/ (exp($hckt/$lambda) -1) *$dlambda;
	$energy+=$bldl *($h*$c/$lambda - $e0);
	$count+=$bldl;
	$lambda+=$dlambda;
	#print "$lambda $bldl\n";
}

$sigma=5.67e-5;
$exp=$sigma*$t**4; #per steradian is wanted, so divide by pi !
$ratio=$energy/$exp;

print  "Temperature      : \t$t (Kelvin)\n";
printf("Flux emitted     : \t%1.3E (ergs /cm^2/ sec) [not per sterad]\n", $energy);
printf("No of photons    : \t%1.3E (cm^2/ sec) [not per sterad]\n", $count);
printf("Fraction         : \t%1.3f (=%1.3E) [Fraction of total blackbody flux]\n",$ratio,$ratio);