#!/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.