#!/usr/bin/perl

#Run this with "perl pdldemo.perl"

use PDL;
use PDL::Graphics::PLplot;
use strict;                          # I recommend always using this

# Global constants

$main::h=6.626e-27;              # erg s
$main::hbar=1.055e-27;           # erg s
$main::k=1.381e-16;              # erg K^-1
$main::ergpereV=1.602e-12;       # erg eV^-1
$main::c=3e10;                   # cm s^-1
$main::pi=3.14159265359;         # I bet there's a built-in for this


# main routine

my ($nu,$Lnu,$starR,$T);

$T=40000;                  # K
$starR=7e10;               # cm -- just a guess, don't use this

($nu,$Lnu)=blackbody($T,$starR);

# Just for fun, plot this using a PDL::Graphics::PLplot command

my $pl = PDL::Graphics::PLplot->new(DEV=>"xwin" , FILE=>$ENV{'DISPLAY'});
$pl->xyplot($nu,$Lnu);
print STDERR "Hit RETURN in the graphics window...\n";
$pl->close();

exit;

# **********************************************************************
# **********************************************************************
# This function creates two PDL arrays, $nu and $Lnu, and fills them with
#   the appropriate values for a blackbody 

sub blackbody {
    my $T=shift;          # temperature in Kelvin
    my $starR=shift;      # star radius in cm

    my ($Lnu,$nu,$i);       # Declare your local variables

    # These next two statements create PDL double arrays of length 1000

    $nu = zeroes(double,1000);
    $Lnu = zeroes(double,1000);

    # Fill the arrays with the blackbody values for 0..100 eV in 0.1 eV

    foreach $i (0..999) {

        # Note that standard Perl arrays are subscripted with [],
        #   but PDL arrays are subscripted with ->index() ; what's more,
        #   for various esoteric reasons, to *set* PDL varaibles using
        #   ->index(), you assign with .= rather than just = ; 
        #   see examples below.
        # I *think* that ->() is a shortcut for ->index()

        $nu->index($i).=$i/10. * $main::ergpereV / $main::h;

        # (Units are $Hz)

        # Below is the slow way to calculate Lnu-- index by index.
        #  It may be clearer, however.  It's comment out here in 
        #  favor of the fast way below.

        # $Lnu->index($i) .= 4 * $main::pi * $starR*$starR *
        #     2*$main::h*($nu->index($i)**3)/($main::c*$main::c) /
        #     ( exp( $main::h * $nu->index($i) / ($main::k * $T) ) -1 );

        # (Units are ergs s^-1 Hz^-1)
    }

    # The beauty of PDL is that it can do things to arrays all at once.
    # This will be a lot faster, too, since the internal code inside PDL
    # is compiledd (while most of Perl is interpreted)

    $Lnu = 4 * $main::pi * $starR*$starR *
        2*$main::h*($nu**3)/($main::c*$main::c) /
        ( exp( $main::h * $nu / ($main::k * $T) ) -1 );

    $Lnu->index(0).=0;         # eliminate a nan
    
    # (Units are ergs s^-1 Hz^-1)


    # The function returns a 2-element array of PDL variables, which
    #   are extracted in the function call in the main program above

    return ($nu,$Lnu);
}