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

# Calculate the surface intensity of a blackbody in a given band
#
# **********
# CGS units !
# **********

$h=6.63e-27;
$c=2.998e10;
$k=1.38e-16;
$hc2=$h*$c*$c;
$step=1000;		#step size as a fraction of wavelength, 1/$step = dlambda / lambda

sub count(){
	($lstart,$lend,$t)=@_;
	$hckt=$h*$c/$k/$t;
	print "Wavelength band: $lstart nm - $lend nm. ";
	$lstart*=1e-7;		#convert nm into cm
	$lend*=1e-7;		#convert nm into cm
	$dlambda=$lstart/$step;
	$lambda=$lstart;
	my $integral=0;
	while ($lambda<$lend){
		$dlambda=$lambda/$step;		#it is important that step size varies with the wavelength !
		$lambda+=$dlambda;
		$bldl=2*$hc2/$lambda**5/ (exp($hckt/$lambda) -1) *$dlambda;
		$integral+=$bldl;
	}

	printf("I_\\lambda = %E erg/cm^2/sec/sr for temperature = %d K\n", $integral, $t);
}

# main part starts here:
# call the function "count" with arguments lambda_min, lambda_max, temperature
# output is the I_lambda for a blackbody of that temperature, integrated from lambda_min to lambda_max

# First, the V band intensity of the sun
&count(505,595,5777);
# Next, the V band intensity for the O7III companion star
&count(505,595,35000);

# calculate total flux from the sun, just for fun:
&count(0.519,519000,5777);	#peak emission for 5777K is 519 nm. integrating from 519/1000 to 519*1000 should give us almost all the flux
# note that I_V (sun) is about 12% of the total energy emitted by the sun.