#!/usr/bin/perl

############################################################################
##	Title	et.pl
##	Usage	./et.pl siteID [-dateYYYY-MM-DD] [-output/directory]
##	Purpose	Calculates the ETo for a given day (default is yesterday),
##		and stores that information in the database.  It also
##		creates PHP pages that implement templates to display the
##		ETo information on the internet.  These files are created in
##		the '/directory' directory if the -output option is used, or
##		the current directory if it is not.
##	Date	8/1/2002
##	Author	Matt Maxwell
############################################################################

## Modules
use DBI;
use Math::Trig;
use Time::Local;

## Subroutines
sub InBounds(@_);		# Determines if the given value is within the bounds set in the database
sub Average(@_);		# Simple procedure to calculate the average of an array of values
sub NewPM(@_);			# Simple procedure to calculate the average of an array of values
sub generateHtmlMain(@_);	# Create the main html page
sub generateHtmlRad(@_);	# Create the radiation html page
sub generateHtmlSoil(@_);	# Create the soil flux html page
sub generateHtmlPsychro(@_);	# Create the psychro html page
sub generateHtmlVPD(@_);	# Create the vapor pressure deficit page
sub generateHtmlVPC(@_);	# Create the vapor pressure curve page


## Variables
$DB = "";		# Variable values removed before displaying on internet
$USERNAME = "";		# Variable values removed before displaying on internet
$PASSWORD = "";		# Variable values removed before displaying on internet
$HELP = "Syntax:\n\t./et.pl siteID [-dateYYYY-MM-DD] [-output/directory]\n";
($ROOTDIR = $0) =~ s/^(.+)\/.*$/\1/;


## Get Command-line Parameters and Variables from the Definition File

die("Error: Must supply a site ID\n$HELP") if ($ARGV[0] eq '');
$siteID = shift @ARGV;

foreach(@ARGV){
	if(/^-date/){
		# Specify the date to calculate (default is yesterday--see below)
		$dateParam=1;
		($mysqlDate = $_) =~ s/^-date//;
		if($mysqlDate !~ /^\d{4}-\d{2}-\d{2}$/){ die ("Error: Unrecognized date\n$HELP"); }
	}elsif(/^-output/){
		# Specify output directory (default is the current directory)
		($baseDir = $_) =~ s/^-output//;
		if($baseDir !~ /\/$/){ $baseDir .= '/'; }	# Add trailing slash if needed
	}
}

if(!$dateParam){
	# Default date is yesterday
	my ($year, $month, $date) = (localtime(time - 24 * 60 * 60))[5,4,3];
	$mysqlDate = sprintf("%04u-%02u-%02u",$year+1900,$month+1,$date);
}

if (-e "$ROOTDIR/$siteID.def"){
	chop(($FullName,$DBStation,$Elevation,$Latitude,$HtmlRootDir,$TemplateDir) = `$ROOTDIR/parseDef.pl $ROOTDIR/$siteID.def`);
}else{
	die "Error: $ROOTDIR/$siteID.def not found\n";
}


## Connect to database

$dbh = DBI->connect("DBI:mysql:database=$DB",$USERNAME,$PASSWORD);

## Query for the day's data

$sth = $dbh->prepare("SELECT 	data_date,
				h_mx_we_airt_1, h_mn_we_airt_1,
				h_mx_we_relh_1, h_mn_we_relh_1,
				h_av_we_wnds_1, h_av_we_srad_1,
				h_to_we_pcpt_1
			FROM ".$DBStation."_hour WHERE data_date LIKE '$mysqlDate%' ORDER BY data_date");
$result = $sth->execute;

# Ensure that the values are changed in the while loop
$tempMax=-9999;
$tempMin=9999;
$relHumMax=-9999;
$relHumMin=9999;
while(my ($hDateTime,$hTempMax,$hTempMin,$hRelHumMax,$hRelHumMin,$hWindSpeed,$hSolRad,$hPrecip)=$sth->fetchrow_array){
	$tempMax = $hTempMax if($hTempMax>$tempMax);	# Get the highest temperature
	$tempMin = $hTempMin if($hTempMin<$tempMin);	# Get the lowest temperature
	$relHumMax = $hRelHumMax if($hRelHumMax>$relHumMax);	# Get the largest relative humidity
	$relHumMin = $hRelHumMin if($hRelHumMin<$relHumMin);	# Get the smallest relative humidity

	if($hSolRad<0){ $hSolRad=0; }

	push(@windSpeeds,$hWindSpeed);	# These two arrays will be used to find the averages
	push(@solRads,$hSolRad);

	$totPrecip += $hPrecip;

	# Store all the values in the @details array
	push(@details,($hDateTime,$hTempMin,$hTempMax,$hRelHumMax,$hRelHumMin,$hWindSpeed,$hSolRad,$hPrecip));
}
$aveWindSpeed = Average(@windSpeeds);
$aveSolRad = Average(@solRads);
$sth->finish;


## Check to see if all of the values are within the appropriate ranges
$badData=0;
# If any of these fields are out of bounds the $badData variable get set to 1
InBounds("h_av_we_airt_1",$tempMax);
InBounds("h_av_we_airt_1",$tempMin);
InBounds("h_av_we_relh_1",$relHumMax);
InBounds("h_av_we_relh_1",$relHumMin);
InBounds("h_av_we_wnds_1",$aveWindSpeed);
InBounds("h_av_we_srad_1",$aveSolRad);


## Calculate ETo
# Get ETo in millimeters and the intermediate values needed to display the html scripts, those values are explained more at the end of the NewPM subroutine
($ET{mm},@IntermediateValues) = NewPM($tempMax,$tempMin,$relHumMax,$relHumMin,$aveWindSpeed,$aveSolRad,$Elevation,$Latitude);
$ET{in} = $ET{mm} * .039370078;		# Convert millimeters to inches


## Store ETo in the database
if($badData){
	$sth = $dbh->prepare("REPLACE INTO ".$DBStation."_day (data_date, insert_date, d_av_we_eto_1, d_av_we_eto_2) VALUES ('$mysqlDate',NOW(),-6666,-6666);");
}else{
	$sth = $dbh->prepare("REPLACE INTO ".$DBStation."_day (data_date, insert_date, d_av_we_eto_1, d_av_we_eto_2) VALUES ('$mysqlDate',NOW(),".$ET{mm}.",".$ET{in}.");");
}
$sth->execute;
$sth->finish;
$dbh->disconnect;


## Create Web Pages
generateHtmlMain(@IntermediateValues);
generateHtmlRad(@IntermediateValues);
generateHtmlSoil(@IntermediateValues);
generateHtmlPsychro(@IntermediateValues);
generateHtmlVPD(@IntermediateValues);
generateHtmlVPC(@IntermediateValues);


#########################################################################
## Done with main program
#########################################################################


## Subrotine Implementations
sub Average(@_){
	my $count = 0;
    	my $total = 0;
    	foreach(@_){
		$count++;
		$total += $_;
    	}
	if($count==0){
		return 0;
	}else{
		return ($total/$count);
	}
}

sub InBounds(@_){
	## Variables
    	#  $DBStation  = a global from the Main section of code
    	#  $dbh        = The database handle 
    	
    	my $tempvar = $_[1];
	
    	my $ibh = $dbh->prepare("SELECT hi_val, low_val  FROM " . $DBStation ."_param WHERE param_id = '$_[0]'");
    	$ibh->execute;
    	my ($upperbound, $lowerbound) = $ibh->fetchrow_array ;
	
    	if($tempvar > $upperbound || $tempvar < $lowerbound){
		$baddata = 1;
		return 0;
    	}
    	return 1;
}

sub NewPM(@_){
    	# 0 - max temp in Celcius 
    	# 1 - min temp in Celcius 
    	# 2 - max relative humity in percent
    	# 3 - min relative humity in percent 
    	# 4 - Wind Speed in miles per hour
    	# 5 - Radiation in kilowatts per sqaure meters
    	# 6 - Elevation in feet
    	# 7 - Latitude
	
	my $maxtemp    = $_[0]; #celcius
	my $mintemp    = $_[1]; #celcius
	my $maxhum     = $_[2]; #percent
	my $minhum     = $_[3]; #percent
	my $windspdtmp = $_[4]; #mph
	my $solrad     = $_[5]; # MJ / m^2
	my $elevation  = $_[6]; #ft
	my $lat        = $_[7];
	
	my $pi = atan2(1,1) *4 ;
	my $albedo = .23;

	# Average Wind Speed In Meter Per Second
	# 1 mph = .44704 m/ss
	my $avg_wind_mps = $windspdtmp*.44704; # m/s
	
	# Mean Temperature
	my $meantemp = ($maxtemp + $mintemp)/2 ;  # celcius
	
	# Elevation
	$elevation*=0.3048;	# Convert to meters
	
	# Julian Day    (yesterday)
	my $jday = `/usr/local/bin/yearday` -1 ; 
	if($jday == 0)  {
		$jday=364;
	}
	
	# dr (inverse relative distance Earth-Sun)
	my $dr = 1+0.033*cos(2*$pi*$jday/365);
	
	# Solar Declination  
	my $soldec =0.409*sin(2*$pi*$jday/365-1.39); 
	
	# Latitude in Radians
	my $radlat = $lat*$pi/180;
	
	# Sunset Hour Angle
	my $sunangle = acos(-tan($radlat)*tan($soldec)) ;
	
	# E0 tmax
	my $eomax = 0.6108 * exp(17.27*$maxtemp/($maxtemp+237.3));
	
	# E0 tmin
	my $eomin =0.6108 * exp(17.27*$mintemp/($mintemp+237.3));
	
	# Mean saturation vapour pressure at the air temp
	my $es = ($eomax +$eomin) /2 ;
	
	# Actual vapour pressure derived from relative humditiy data
	my $ea = ($eomax*$minhum/100+$eomin*$maxhum/100)/2 ;
	
	# Atmospheric Pressure
	my $p = 101.3*((293-0.0065*$elevation)/293) **  5.26 ;
	
	# Lambda (psychrometric constant)  =.00163*$p/lhv   where lhv=latent heat of vaporization=2.45
	my $lambda =0.000665*$p;
	
	# Delta (slope vapour pressure)
	my $delta = 4098*(0.6108*exp(17.27*$meantemp/($meantemp+237.3)))/($meantemp+237.3)**2 ;
	
	# Ra - extraterrestrial radiation
	my $ra = 24*60*0.082/$pi*$dr*($sunangle*sin($radlat)*sin($soldec)+cos($radlat)*cos($soldec)*sin($sunangle));
	
	# Rso (clear sky solar radiation)
	my $rso = (0.75 + 2*10**-5*$elevation) * $ra ;
	
	# Rns (Net solar or net shortware radiation)
	my $rns = (1-$albedo)*$solrad ;
	
	# Rnl (net outgoing longwave radiation)
	my $rnl  = 4.903*10** -9*((($maxtemp+273.16)**4)+($mintemp+273.16)**4) / 2 * (0.34-0.14*$ea**0.5) * (1.35*$solrad/$rso-0.35) ;
	
	#Rn (net solar radiation)
	my $rn = $rns-$rnl ;
	
	#et (reference evapotranspiration - mm)
	
	my $et = (0.408*$delta*($rn-0)+$lambda*(900/($meantemp+273))*$avg_wind_mps*($es-$ea))/ ($delta+$lambda*(1+0.34*$avg_wind_mps)) ;
	#in case we try to take the root of a negitive number
	if($et eq "nan"){$et = -6666;}

	## Return values:
	#	$et		Evapotranspiration in millimeters
	#	$rn		Net Radiation
	#	$lambda		Psychrometic Constant
	#	$es-$ea		Vapor Pressure Deficit
	#	$delta		Slope Vapor Pressure Curve
	#	$rnl		Longwave Radiation
	#	$ea		Dew Point Vapor Pressure
	#	$rso		Clear Sky Shortwave Radiation
	#	$ra		Extra-Terrestrial Radiation
	#	$dr		Earth-Sun Relative Distance
	#	$jday		Julian Day of the Year
	#	$sunangle	Angle at Sunset
	#	$radlat		Latitude in radians	
	#	$soldec		Solar Declination
	#	$p		Atmospheric Pressure
	#	$es		Saturation Vapor Pressure
	#	$ea		Actual Vapor Pressure	
	#	$eomax		eMax
	#	$eomin		eMin

	return ($et,$rn,$lambda,$es-$ea,$delta,$rnl,$ea,$rso,$ra,$dr,$jday,$sunangle,$radlat,$soldec,$p,$es,$ea,$eomax,$eomin);
}

sub generateHtmlMain(@_){
	my($NetRadiation,$PsychrometicConstant,$VaporPressureDeficit,$SlopeVaporPressureCurve) = @_[0,1,2,3];

	open HTML, ">".$baseDir."index.php" or die "Error. Can't create ".$baseDir."index.php: $!\n";
	print HTML "<?\n";
	print HTML "\t\$FullName = \"$FullName\";\n";
	print HTML "\t\$BadData = \"$badData\";\n";
	print HTML "\t\$Date = \"$mysqlDate\";\n";
	my $max = scalar(@details) /8;
	for(my $i = 0; $i < $max; $i++){
		my $j = $i *8;
		print HTML "\t\$DataArray[$i]{datetime} = \"".$details[$j]."\";\n";
		print HTML "\t\$DataArray[$i]{lowTemp} = \"".sprintf("%0.1f",$details[$j+1])."\";\n";
		print HTML "\t\$DataArray[$i]{highTemp} = \"".sprintf("%0.1f",$details[$j+2])."\";\n";
		print HTML "\t\$DataArray[$i]{lowHumidity} = \"".sprintf("%0.1f",$details[$j+3])."\";\n";
		print HTML "\t\$DataArray[$i]{highHumidity}= \"".sprintf("%0.1f",$details[$j+4])."\";\n";
		print HTML "\t\$DataArray[$i]{wind} = \"".sprintf("%0.1f",$details[$j+5])."\";\n";
		print HTML "\t\$DataArray[$i]{solRad} = \"".sprintf("%0.3f",$details[$j+6])."\";\n";
		print HTML "\t\$DataArray[$i]{precip} = \"".sprintf("%0.2f",$details[$j+7])."\";\n";
	}
	print HTML "\t\$MinTemp = \"".sprintf("%0.1f",$tempMin)."\";\n";
	print HTML "\t\$MaxTemp = \"".sprintf("%0.1f",$tempMax)."\";\n";
	print HTML "\t\$MinRelHumidity = \"".sprintf("%0.1f",$relHumMin)."\";\n";
	print HTML "\t\$MaxRelHumidity = \"".sprintf("%0.1f",$relHumMax)."\";\n";
	print HTML "\t\$MeanTemp = \"".sprintf("%0.1f",($tempMax+$tempMin)/2)."\";\n";
	print HTML "\t\$TotalRad = \"".sprintf("%0.3f",$aveSolRad)."\";\n";
	print HTML "\t\$ET{in} = \"".sprintf("%0.2f",$ET{in})."\";\n";
	print HTML "\t\$ET{mm} = \"".sprintf("%0.2f",$ET{mm})."\";\n";
	print HTML "\t\$NetRadiation = \"".sprintf("%0.4f",$NetRadiation)."\";\n";
	print HTML "\t\$AverageTemp = \"".sprintf("%0.1f",($tempMax+$tempMin)/2)."\";\n";
	print HTML "\t\$AverageWind= \"".sprintf("%0.4f",$aveWindSpeed*.44704)."\";\n";		# Convert to m/s
	print HTML "\t\$PsychrometicConstant= \"".sprintf("%0.4f",$PsychrometicConstant)."\";\n";
	print HTML "\t\$VaporPressureDeficit= \"".sprintf("%0.4f",$VaporPressureDeficit)."\";\n";
	print HTML "\t\$SlopeVaporPressureCurve= \"".sprintf("%0.4f",$SlopeVaporPressureCurve)."\";\n";
	print HTML "\t\$RootDir= \"".$HtmlRootDir."\";\n";
	print HTML "\n";
	print HTML "\tinclude \"$TemplateDir/main.inc\";\n";
	print HTML "?>\n";
	close HTML;
}

sub generateHtmlRad(@_){
	my($NetRadiation,$LongwaveRad,$DewPointVP,$ClearSkyShortwaveRad,$ExtraTerrRad,$RelativeDistance,$DayOfYear,$SunsetAngle,$radLat,$SolarDeclination)=@_[0,4,5,6,7,8,9,10,11,12];

	open HTML, ">".$baseDir."radiation.php" or die "Error. Can't create ".$baseDir."radiation.php: $!\n";
	print HTML "<?\n";
	print HTML "\t\$FullName = \"$FullName\";\n";
	print HTML "\t\$NetRadiation= \"".sprintf("%0.4f",$NetRadiation)."\";\n";
	print HTML "\t\$ShortwaveRad= \"".sprintf("%0.3f",$aveSolRad)."\";\n";
	print HTML "\t\$LongwaveRad= \"".sprintf("%0.4f",$LongwaveRad)."\";\n";
	print HTML "\t\$MaxTempK= \"".sprintf("%0.2f",$tempMax+273.16)."\";\n";
	print HTML "\t\$MinTempK= \"".sprintf("%0.2f",$tempMin+273.16)."\";\n";
	print HTML "\t\$DewPointVP= \"".sprintf("%0.4f",$DewPointVP)."\";\n";
	print HTML "\t\$ClearSkyShortwaveRad= \"".sprintf("%0.4f",$ClearSkyShortwaveRad)."\";\n";
	print HTML "\t\$Elevation= \"".sprintf("%0.0f",$Elevation * .3048)."\";\n";
	print HTML "\t\$ExtraTerrRad= \"".sprintf("%0.3f",$ExtraTerrRad)."\";\n";
	print HTML "\t\$RelativeDistance= \"".sprintf("%0.4f",$RelativeDistance)."\";\n";
	print HTML "\t\$DayOfYear= \"".$DayOfYear."\";\n";
	print HTML "\t\$SunsetAngle= \"".sprintf("%0.4f",$SunsetAngle)."\";\n";
	print HTML "\t\$Latitude= \"".sprintf("%0.4f",$radLat)."\";\n";
	print HTML "\t\$SolarDeclination= \"".sprintf("%0.4f",$SolarDeclination)."\";\n";
	print HTML "\t\$RootDir= \"".$HtmlRootDir."\";\n";
	print HTML "\n";
	print HTML "\tinclude \"$TemplateDir/radiation.inc\";\n";
	print HTML "?>\n";
	close HTML;
}	

sub generateHtmlSoil(@_){
	open HTML, ">".$baseDir."soil.php" or die "Error. Can't create ".$baseDir."soil.php: $!\n";
	print HTML "<?\n";
	print HTML "\t\$FullName = \"$FullName\";\n";
	print HTML "\n";
	print HTML "\tinclude \"$TemplateDir/soil.inc\";\n";
	print HTML "?>\n";
	close HTML;
}

sub generateHtmlPsychro(@_){
	my($PsychrometricConstant,$AtmosphericPressure) = @_[1,13];
	
	open HTML, ">".$baseDir."psychrometic.php" or die "Error. Can't create ".$baseDir."psychrometic.php: $!\n";
	print HTML "<?\n";
	print HTML "\t\$FullName = \"$FullName\";\n";
	print HTML "\t\$PsychrometricConstant= \"".sprintf("%0.4f",$PsychrometricConstant)."\";\n";
	print HTML "\t\$AtmosphericPressure= \"".sprintf("%0.3f",$AtmosphericPressure)."\";\n";
	print HTML "\t\$Elevation= \"".sprintf("%0.0f",$Elevation*.3048)."\";\n";
	print HTML "\t\$RootDir= \"".$HtmlRootDir."\";\n";
	print HTML "\n";
	print HTML "\tinclude \"$TemplateDir/psychrometic.inc\";\n";
	print HTML "?>\n";
	close HTML;
}	

sub generateHtmlVPD(@_){
	my($VPDeficit,$SaturationVP,$ActualVP,$eMax,$eMin) = @_[2,14,15,16,17];

	open HTML, ">".$baseDir."vpd.php" or die "Error. Can't create ".$baseDir."vpd.php: $!\n";
	print HTML "<?\n";
	print HTML "\t\$FullName = \"$FullName\";\n";
	print HTML "\t\$VPDeficit= \"".sprintf("%0.4f",$VPDeficit)."\";\n";
	print HTML "\t\$SaturationVP= \"".sprintf("%0.4f",$SaturationVP)."\";\n";
	print HTML "\t\$ActualVP= \"".sprintf("%0.4f",$ActualVP)."\";\n";
	print HTML "\t\$eMax= \"".sprintf("%0.4f",$eMax)."\";\n";
	print HTML "\t\$eMin= \"".sprintf("%0.4f",$eMin)."\";\n";
	print HTML "\t\$TMax= \"".sprintf("%0.1f",$tempMax)."\";\n";
	print HTML "\t\$TMin= \"".sprintf("%0.1f",$tempMin)."\";\n";
	print HTML "\t\$RHMax= \"".sprintf("%0.3f",$relHumMax)."\";\n";
	print HTML "\t\$RHMin= \"".sprintf("%0.3f",$relHumMin)."\";\n";
	print HTML "\t\$RootDir= \"".$HtmlRootDir."\";\n";
	print HTML "\n";
	print HTML "\tinclude \"$TemplateDir/vpd.inc\";\n";
	print HTML "?>\n";
	close HTML;
}	
sub generateHtmlVPC(@_){
	my($SlopeVPCurve) = @_[3];

	open HTML, ">".$baseDir."vpc.php" or die "Error. Can't create ".$baseDir."vpc.php: $!\n";
	print HTML "<?\n";
	print HTML "\t\$FullName = \"$FullName\";\n";
	print HTML "\t\$SlopeVPCurve= \"".sprintf("%0.4f",$SlopeVPCurve)."\";\n";
	print HTML "\t\$MeanTemp= \"".sprintf("%0.2f",($tempMax+$tempMin)/2)."\";\n";
	print HTML "\t\$RootDir= \"".$HtmlRootDir."\";\n";
	print HTML "\n";
	print HTML "\tinclude \"$TemplateDir/vpc.inc\";\n";
	print HTML "?>\n";
	close HTML;
}