view release on metacpan or  search on metacpan

lib/Drought/PET/Thornthwaite.pm  view on Meta::CPAN

sub pet_thornthwaite {
    # Define a NaN value for later use
    my $inf           = exp(~0 >> 1);
    my $nan           = $inf / $inf;

    # Make sure the caller supplied sufficient arguments
    unless(@_ and @_ >= 5) { croak "Insufficient arguments supplied"; }

    # Get the arguments
    my $yday        = shift;
    my $ndays       = shift;
    my $lat         = shift;
    my $temperature = shift;
    my $tei         = shift;
    my $missing_val = $nan;

    # Replace the default missing value if the caller supplied one
    if(@_) { $missing_val = shift; }

    # Validate the non-data arguments
    if(not looks_like_number($yday))       { croak "YDAY arg must be numeric";           }
    $yday  = int($yday);
    if(not $yday >= 1 or not $yday <= 366) { croak "YDAY arg must be between 1 and 366"; }
    if(not looks_like_number($ndays))      { croak "NDAYS arg must be numeric";          }
    $ndays = int($ndays);
    if(not $ndays >= 1)                    { croak "NDAYS arg must be 1 or greater";     }
    if(not looks_like_number($lat))        { croak "LAT arg must be numeric";            }
    if(not $lat >= 0 or not $lat <= 90)    { croak "LAT arg must be between 0 and 90";   }

    # Return the missing value if the data arguments are missing or invalid
    if(
        not defined($temperature)           or
        not defined($tei)                   or
        not looks_like_number($temperature) or
        not looks_like_number($tei)         or
        $temperature == $missing_val        or
        $tei         == $missing_val        or
        not defined($temperature <=> 0)     or
        not defined($tei <=> 0)
    ) { return $missing_val; }

    # Calculate the PET!

    if($temperature <= 0) { return 0; }

    elsif($temperature > 0 and $temperature < 26.5) {
        # Calculate the hours of sunlight
        my $lat_radians       = (pi/180)*$lat;
        my $solar_declination = 0.409*sin((2*pi/365)*$yday - 1.39);
        my $sunset_angle      = acos(-tan($lat_radians)*tan($solar_declination));
        my $daylength         = (24/pi)*$sunset_angle;
        # Calculate alpha term
        my $alpha             = (6.75e-7)*($tei**3) - (7.71e-5)*($tei**2) + 0.01792*$tei + 0.49239;
        # Apply the original Thornthwaite Equation
        my $pet_gross         = 16*((10*$temperature)/($tei))**$alpha;
        return $pet_gross*($daylength/12)*($ndays/30);
    }

    else {
        # Apply the hot temperature formula described in Huang et al. 1996
        return (-415.85 + 32.25*$temperature - 0.43*($temperature**2))*($ndays/30);
    }

}

=head2 tei_thornthwaite

 my $tei = tei_thornthwaite($jan,$feb,$mar,$apr,$may,$jun,$jul,$aug,$sep,$oct,$nov,$dec);

Calculates and returns the temperature efficiency index (TEI, sometimes called the 
Thornthwaite heat index) based on a monthly average temperature climatology for a given 
location. Since the equation involves a summation of adjusted monthly temperatures, 12 
arguments are required, consisting of the average temperature for each calendar month of 
the year. An optional 13th argument can also be supplied to define a numeric value 
interpreted as missing data. If not supplied, the default missing value is NaN.

The missing value will be returned if any of the temperature values are missing, 
undef, or invalid (e.g., non-numeric).

=cut

sub tei_thornthwaite {
    # Define a NaN value for later use
    my $inf           = exp(~0 >> 1);
    my $nan           = $inf / $inf;

    # Make sure the caller supplied sufficient arguments
    unless(@_ and @_ >= 12) { croak "Insufficient arguments supplied"; }

    # Get the arguments
    my @monthly_temperatures;
    my $invalid = 0;
    for(my $i=0; $i<12; $i++) { push(@monthly_temperatures,shift); } 
    my $missing_val = $nan;

    # Replace the default missing value if the caller supplied one
    if(@_) { $missing_val = shift; }

    # Calculate the TEI!
    my $tei = 0;

    for(my $i=0; $i<12; $i++) {
        # Return a missing value if any of the monthly temps are missing or invalid
        if(
            not defined($monthly_temperatures[$i])           or
            not looks_like_number($monthly_temperatures[$i]) or
            $monthly_temperatures[$i] == $missing_val        or
            not defined($monthly_temperatures[$i] <=> 0)
        ) { return $missing_val; }
        # Set any monthly mean temp below 0 to 0
        $monthly_temperatures[$i] = $monthly_temperatures[$i] > 0 ? $monthly_temperatures[$i] : 0;
        $tei += ($monthly_temperatures[$i]/5)**1.514;
    }

    return $tei;
}

=head1 INSTALLATION

The best way to install this module is with a CPAN client, which will resolve and