Untitled Document

#!/usr/bin/perl

# Written by William Patrick Arnott, 4 July 2009, from Reno NV.
# This is a PERL language script.
# The whole idea of this script is to use it to write html code for a browser, and to
# allow dynamic content where the user may change conditions and easily resubmit the calculation.
# Version control: this version was completed 25 July 2011 at 3:25 pm.
# Still left undone is to compare the Rayleigh optical depths with values from the literature across the vastness of the literature.

# Set the color scheme.
$text_color = "000000";
$page_background_color = "FFFFF0";
$page_links_color = "0000FF";
$page_active_links_color = "FF0000";
$page_visited_links_color = "0000FF";
$table_background_color= "F4E3C8";
$table_cell_color= "FFFFFF";

# The following gets the environmental variable related to the name of this script for use later.
# This is an example of how to define a variable cgu_url in the PERL language.
# CGI is defined as "Common Gateway Interface". It is a standard for interfacing external
# applications such as perl programs to web servers.
$cgi_url = "$ENV{'SCRIPT_NAME'}";

# $header is the HTML header that will be generated by the results page
# of the program. It is nothing more than a general string variable, and could have any name.
# We simply write the string to the web page later.

$header = qq[

<html><head><title>Calculate Air Mass and Rayleigh Optical Depth for SunPhotometer Data Analysis</title>
<meta name="Description" content="Calculates the solar zenith angle for Reno and the Rayleigh Scattering Coefficient for air.">
<meta name="keywords" content="Rayleigh Reno sunphotometer sun photometer handheld sunphotometer">
</head>
<body bgcolor=#$page_background_color text=#$text_color link=#$page_links_color vlink=#$page_visited_links_color alink=#$page_active_links_color>
<center>
<SCRIPT type="text/javascript" src="overlib.js"></SCRIPT>

];

# $footer is the HTML code that will go at the bottom of pages
# generated by the program.

$footer = qq[

 
<iframe src="http://free.timeanddate.com/clock/i2oq821e/n599/fn6/fs18/fcff0/tc000/ftb/bas2/bat1/bacfff/pa8/tt0/tw1/th1/ta1/tb4" frameborder="0" width="242" height="64"></iframe>
</body></html>
<center>Powered by healthy green tea smoothies ... (:-:)</center>
<center>By W. Patrick Arnott and Ben Sumlin</center>

];

# This line makes it so everything printed by the script is printed
# immediately instead of being buffered and printed whenever.

$| = 1;

# This line calls the parse_form_data subroutine, which grabs all the
# input from the form into an associative style array so we can use the
# variables. This is where data is actually read into the script from the standard input,
# which is the input form itself considering of the pressure, temperature, and wavelength.

%data = &parse_form_data();

# Use the local time PDT to initialize the time and date.
# Notice that if the pressure variable is erased at any time it will update the local time to the current time and date.
# That's a very convenient function.

$ENV{TZ} = 'America/Los_Angeles';
@timeData = localtime();

# Selectively initialize the associative data array.
if ($data{'presADA'} eq "") {
$data{'presADA'} = 850.00 ;
$data{'timezoneADA'} = -7.00 ; # Pacific daylight savings time time.
$data{'hourADA'} = $timeData[2] ;
$data{'minADA'} = $timeData[1] ;
$data{'secADA'} = $timeData[0] ;
$data{'monthADA'}= $timeData[4]+1 ;
$data{'dayADA'} = $timeData[3] ;
$data{'yearADA'} = $timeData[5]+1900 ;
$data{'latADA'} = 39.5411666 ; # Start with UNR Reno value.
$data{'lonADA'}= -119.81406 ; # Start with UNR Reno value.
$data{'wave'} = 532.00 ;
}

#################################################################
############## Here we go, make the web page. ##############
#################################################################

# This line must be in any CGI script: it tells your web browser
# that HTML stuff is coming out. If you had a CGI script that showed
# an image, you might use "Content-type: image/gif", but you have to
# have a content-type (or MIME header) of some kind. If you were using
# the Perl CGI module with "use CGI;" then you could probably do something
# like "print header;" but a lot of servers don't have the
# CGI module, so we're taking no chances.

print "Content-type: text/html\n\n";

# First, print our HTML for the top of the page.
print $header;

# This next set of things is just checking so that if they enter
# strange values. Just show a little red * next to the form
# field, and don't bother trying to calculate anything.

# Coerce values to good operating ranges.
# Pressure coersion.
$errpressure = "*error";
if ($data{'presADA'} > 0.0 and $data{'presADA'} <= 10000.0)
{ $errpressure="" }
elsif ($data{'presADA'} < 10.0) { $data{'presADA'}=10.0 }
else {$data{'presADA'}=850.0 ; }
$pressure=$data{'presADA'} ;

# Timezone coersion.
$errtimezone = "*error";
if ($data{'timezoneADA'} >= -12.0 and $data{'timezoneADA'} <= 12.0)
{ $errtimezone="" }
elsif ($data{'timezoneADA'} < -12.0) { $data{'timezoneADA'}=-12.0 }
else {$data{'timezoneADA'}=12.0 ; }
$timezonechosen=$data{'timezoneADA'} ;
$timezonechosen=sprintf("%.0f", $timezonechosen) ;

# Hour coersion.
$errhour = "*error";
if ($data{'hourADA'} >= 0.0 and $data{'hourADA'} <= 23.0)
{ $errhour="" }
elsif ($data{'hourADA'} < 0.0) { $data{'hourADA'}=12.0 }
else {$data{'hourADA'}=12.0 ; }
$hourchosen=$data{'hourADA'} ;
$hourchosen=sprintf("%.0f", $hourchosen) ;

# Minute coersion.
$errmin = "*error";
if ($data{'minADA'} >= 0.0 and $data{'minADA'} <= 59.0)
{ $errmin="" }
elsif ($data{'minADA'} < 0.0) { $data{'minADA'}=0.0 }
else {$data{'minADA'}=59.0 ; }
$min=$data{'minADA'};
$min=sprintf("%.0f",$min) ;

# Seconds coersion.
$errsec = "*error";
if ($data{'secADA'} >= 0.0 and $data{'secADA'} <= 59.0)
{ $errsec="" }
elsif ($data{'secADA'} < 0.0) { $data{'secADA'}=0.0 }
else {$data{'secADA'}=59.0 ; }
$sec=$data{'secADA'} ;

# Month coersion.
$errmon = "*error";
if ($data{'monthADA'} >= 1.0 and $data{'monthADA'} <= 12.0)
{ $errmon="" }
elsif ($data{'monthADA'} < 1.0) { $data{'monthADA'}=1.0 }
else {$data{'monthADA'}=12.0 ; }
$mon=$data{'monthADA'};
$mon=sprintf("%.0f",$mon) ;

# Day coersion. Here we have to constrain, for now, to having a month of 31 days or less.
# Sorry february. It is possible to read the month and to coerce the day accordingly.
$errday = "*error";
if ($data{'dayADA'} >= 1.0 and $data{'dayADA'} <= 31.0)
{ $errday="" }
elsif ($data{'dayADA'} < 1.0) { $data{'dayADA'}=1.0 }
else {$data{'dayADA'}=31.0 ; }
$daychosen=$data{'dayADA'};
$daychosen=sprintf("%.0f",$daychosen) ;

# Year coersion. This is a loosely needed coersion.
$erryear = "*error";
if ($data{'yearADA'} >= 1898.00 and $data{'yearADA'} <= 2300.0)
{ $erryear="" }
elsif ($data{'yearADA'} < 1898.0) { $data{'yearADA'}=1898.0 }
else {$data{'yearADA'}=2011.0 ; }
$yearchosen=$data{'yearADA'} ;
$yearchosen=sprintf("%.0f",$yearchosen) ;

# Latitude coersion.
$errlat = "*error";
if ($data{'latADA'} >= -90.0 and $data{'latADA'} <= 90.0)
{ $errlat="" }
elsif ($data{'latADA'} < -90.0) { $data{'latADA'}=-90.0 }
else {$data{'latADA'}=90.0 ; }
$lat=$data{'latADA'} ;

# Longitude coersion.
$errlon = "*error";
if ($data{'lonADA'} >= -180.0 and $data{'lonADA'} <= 180.0)
{ $errlon="" }
elsif ($data{'lonADA'} < -180.0) { $data{'lonADA'}=-180.0 }
else {$data{'lonADA'}=180.0 ; }
$lon=$data{'lonADA'} ;

# Wavelength coersion.
$errwavelength = "*error";
if ($data{'wave'} > 299.999 and $data{'wave'} <= 1000.0)
{ $errwavelength="" }
elsif ($data{'wave'} < 300.0) { $data{'wave'}=300.0 }
else {$data{'wave'}=1000.0 ; }
$wavelength=$data{'wave'} ;

&RayleighCoefs ; # Calls the subroutine to compute Rayleigh Coefficients.

&AirMassCalculator ; # Calls the subroutine to compute the air mass.

$RAY = qq[
 
<table border="6" bordercolor="#336633" cellpadding="2" cellspacing="2">
<caption>Calculated Results</caption>
<tr>
<th>Air Mass at chosen date and time </th>
<td>$AirMassOut</td>
</tr>
<tr>
<th>Atmospheric Column Rayleigh Scattering Optical Depth
 at Chosen Pressure and Wavelength
<A href="javascript:void(0);" onClick="return overlib(' theory from Rayleigh-scattering Calculations for the Terrestrial Atmosphere. Applied Optics 34, 1995, 2765-2773. Fk from the table was fit to a 10th order polynomial.', STICKY, CAPTION, 'A. Bucholtz', ABOVE);"onmouseout="nd();">*</A>
</th>
<td>$Tau[2]</td>
</tr>
</table>
];

# The get_form subroutine prints out the form, saves whatever
# they entered, etc.

print &get_form ;

print "$RAY";

print $footer;

exit();

# This subroutine write the HTML code for setting up the calculation.
sub get_form
{
$form = qq[
<table border=5 bordercolor="#336699" cellpadding="2" cellspacing="2" >
<caption>
<A href="javascript:void(0);" onClick="return overlib('Authors of this page, written in PERL, and calculations contained therein, with references as noted. Enjoy!', STICKY, CAPTION, 'W. Patrick Arnott and Benjamin Sumlin', LEFT);"onmouseout="nd();">Air mass and Rayleigh Scattering Coefficient for Sunphotometer Analysis </A> 
</caption>
<tr>
<td><table><form action=$cgi_url method="POST">
<tr><td align=right>$errpressure <A href="javascript:void(0);" onClick="return overlib('10 to 2000', STICKY, CAPTION, 'Range of Values', ABOVE);"onmouseout="nd();">Pressure (mb)</A></td>
<td><input type=text name=presADA size=14 value=\"$data{'presADA'}\"></td>
</tr>
<tr> <td align=right>$errtimezone <A href="javascript:void(0);" onClick="return overlib('-12 to 12. For Reno NV, use -7 for daylight savings time in summer, and -8 otherwise.', STICKY, CAPTION, 'Range of Values', ABOVE);"onmouseout="nd();">Time zone</A> </td>
<td><input type=text name=timezoneADA size=14 value=\"$data{'timezoneADA'}\"></td>
</tr>
<tr><td align=right>$errhour <A href="javascript:void(0);" onClick="return overlib('0 to 23', STICKY, CAPTION, 'Range of Values', ABOVE);"onmouseout="nd();">Hour</A> </td>
<td><input type=text name=hourADA size=14 value=\"$data{'hourADA'}\"></td>
</tr>
<tr><td align=right>$errmin <A href="javascript:void(0);" onClick="return overlib('0 to 59', STICKY, CAPTION, 'Range of Values', ABOVE);"onmouseout="nd();">Minute</A></td>
<td><input type=text name=minADA size=14 value=\"$data{'minADA'}\"></td>
</tr>
<tr><td align=right>$errsec <A href="javascript:void(0);" onClick="return overlib('0 to 59', STICKY, CAPTION, 'Range of Values', ABOVE);"onmouseout="nd();">Seconds</A></td>
<td><input type=text name=secADA size=14 value=\"$data{'secADA'}\"></td>
</tr>
<tr><td align=right>$errmon <A href="javascript:void(0);" onClick="return overlib('1 to 12', STICKY, CAPTION, 'Range of Values', ABOVE);"onmouseout="nd();">Month</A></td>
<td><input type=text name=monthADA size=14 value=\"$data{'monthADA'}\"></td>
</tr>
<tr><td align=right>$errday <A href="javascript:void(0);" onClick="return overlib('1 to 31', STICKY, CAPTION, 'Range of Values', ABOVE);"onmouseout="nd();">Day</A></td>
<td><input type=text name=dayADA size=14 value=\"$data{'dayADA'}\"></td>
</tr>
<tr><td align=right>$erryear <A href="javascript:void(0);" onClick="return overlib('1848 to 2031', STICKY, CAPTION, 'Range of Values', ABOVE);"onmouseout="nd();">Year</A></td>
<td><input type=text name=yearADA size=14 value=\"$data{'yearADA'}\"></td>
</tr>
<tr><td align=right>$errlat <A href="javascript:void(0);" onClick="return overlib('-90 to 90', STICKY, CAPTION, 'Range of Values', ABOVE);"onmouseout="nd();">Latitude</A></td>
<td><input type=text name=latADA size=14 value=\"$data{'latADA'}\"></td>
</tr>
<tr><td align=right>$errlon <A href="javascript:void(0);" onClick="return overlib('-180 to 180', STICKY, CAPTION, 'Range of Values', ABOVE);"onmouseout="nd();">Longitude</A></td>
<td><input type=text name=lonADA size=14 value=\"$data{'lonADA'}\"></td>
</tr>
<tr><td align=right>$errwavelength <A href="javascript:void(0);" onClick="return overlib('300 nm to 1000 nm', STICKY, CAPTION, 'Range of Values', ABOVE);"onmouseout="nd();">Wavelength (nm)</A></td>
<td><input type=text name=wave size=14 value=\"$data{'wave'}\"></td>
</tr>
<tr><td align=right><input type="submit" value="Calculate"></td></tr>
</form></table><td/>
<td><img src="sunPhotometerBasics.png" width="320" height="298"></td>
</tr>
</table>
];
$form;
}

# This is one amazing subroutine. It unpacks data from the standard input, and makes it useful
# to proceed with the calculation. Note that variable derive their definition in the form
# subroutine.
sub parse_form_data {
my($string,%data,@data);

# get data
if ($ENV{'REQUEST_METHOD'} eq 'GET') {
$_ = $string = $ENV{'QUERY_STRING'};
tr/\"~;/_/;
$string = $_;
}

else { read(STDIN, $string, $ENV{'CONTENT_LENGTH'});
$_ = $string;
$OK_CHARS='a-zA-Z0-9=&%\n\/_\-\.@';
tr/\"~;/_/;
$string = $_;
}

# split data into name=value pairs
# Numerical values are in variable string. They are separated by &.
@data = split(/&/, $string);

# split into name=value pairs in associative array
foreach (@data) {
split(/=/, $_);
$_[0] =~ s/\+/ /g; # plus to space
$_[0] =~ s/%(..)/pack("c", hex($1))/ge; # hex to alphanumeric
$data{"$_[0]"} = $_[1];
}

# translate special characters
foreach (keys %data) {
$data{"$_"} =~ s/\+/ /g; # plus to space
$data{"$_"} =~ s/%(..)/pack("c", hex($1))/ge; # hex to alphanumeric
}

%data; # return associative array of name=value
}

sub RayleighCoefs {

# Arrays are filled as follows: 0=Argon, 1=Carbon Dioxide, 2=Air .
@Rgas=(8.3143*1000.0/39.95,8.3143*1000.0/44.01,8.3143*1000.0/28.98) ; # Gas constants.
@conc=(0.0093,0.000380,1.0) ; # Partial Pressure of Argon, CO2, and air for optical depth calculation.
$g=9.80665 ; # acceleration due to gravity at the surface.
$l=$wavelength ;
$T=$temperature ;
$P=$pressure ;

$nu=1.e-2/($l*1.e-9) ; # Wavenumber in 1/cm.
$nu2=$nu**2 ;
$wl2=$l/1000.0 ; $wl2=$wl2**2 ; # Wavelength in microns, and it squared.

$T=$T+273.15 ; $T0=273.15+15.0 ; $P0=1013.25 ;

$fac=$P*$T0/($P0*$T); # Used to get the gas concentration at T and P.

$N0=25.4743e18 ; # Gas concentration at T0 and P0, 1/cm^3.

$pi=4.0*atan2(1.0,1.0) ;

@Fkk=(1.0,1.0,1.0) ; # Depolarization factor initilization.
@sig=(0.0, 0.0, 0.0) ; # Rayleigh cross section in cm^2 / molecule initilization.
@beta_par=(0.0, 0.0, 0.0) ; # Forward and Backward Scatter parallel in Mm-1.
@beta_per=(0.0, 0.0, 0.0) ; # perpendicular .
@beta=(0.0, 0.0, 0.0) ; # Total scattering cross section.
@nm1_10to4th=(0.0, 0.0, 0.0) ; # Refractive index at 15 C and 1013.25 mb.
@tau=(0.0,0.0,0.0) ; # optical depth of the atmosphere for each gas according to its concentration.

for ($gas=1; $gas <= 3; $gas++) {

if ($gas==1) { # Case of Argon!!!
$f1=6432.135 + 286.06021e12/(14.4e9-$nu2) ;
$ref=1.0 + 1.e-8 * $f1 ;
$Fk=1.0 ;

} elsif ($gas==2) { # case of Carbon Dioxide!!!
$f1= 5799.25/(128908.9**2 - $nu2)
+ 120.05/(89223.8**2 - $nu2)
+ 5.3334/(75037.5**2-$nu2)
+ 4.3244/(67837.7**2-$nu2)
+ 0.1218145e-4/(2418.136**2-$nu2) ;
$ref=1.0 + 1116.63027063992 * $f1 ;
# print("co2 refractive Index From adjusted model = ") ; print($ref) ; print("\n") ;
$Fk=1.1364 + 25.3e-12 * $nu2 ;
$sigCO2= 27.2e-45 * $nu**4.1343 ; # Emperical Rayleigh for CO2 from Eq. 15 of Sneep's paper.
$a=sqrt($sigCO2*$N0**2/($Fk*24.0*$pi**3*$nu**4)) ;
$ref=sqrt((2.0*$a+1.0)/(1.0-$a)) ;
$factor=($ref-1.0)/$f1 ;
# print("emperical factor for ref index = ") ; print($factor) ; print("\n") ;
# print("emperical co2 Rayleigh Cross Section = ") ; print($sigCO2) ; print(" cm-1 \n") ;

$betas=$N0*$sigCO2 ;
$beta=$betas*$fac*1.0e8 ;
# print("emperical co2 Rayleigh Coefficient = ") ; print($beta) ; print(" Mm-1 \n") ;

} else { # assume ($gas==3) # case of Air!!!
if ($l => 230.0) {
$f1=5791817.0/(238.0185-1.0/$wl2) ;
$f1=$f1+167909.0 / (57.362 - 1.0/$wl2 ) ;
} else {
$f1=8060.51+2480990.0/(132.274-1/$wl2)+17455.7/(39.32957-1.0/$wl2);
}
$ref=1.0 + 1.0e-8*$f1 ;
$a=$l;
$Fk=5.427964E-23*$a**8 - 2.806622E-19*$a**7 + 6.207856E-16*$a**6
-7.662177E-13*$a**5 + 5.767936E-10*$a**4 - 2.712549E-7*$a**3 +
7.799646E-5*$a**2 - 1.261146E-2*$a + 1.937987E+0 ;
# $Fk = 1.049 ; # From Table 2 of Bucholtz et. al. for 532 nm.
}

$rhon=6.0*($Fk-1.0)/(3.0+7.0*$Fk) ;
$gamma=$rhon/(2.0-$rhon) ;
$r=(1.0 + 5.0*$gamma) / (3.0*(1.0 + $gamma)) ; # Fpar/Fper = Bpar/Bper .

# print("Refractive Index= ") ; print($ref) ; print(" \n") ;
# print("rho_n= ") ; print($rhon) ; print(" \n") ;
# print("gamma= ") ; print($gamma) ; print(" \n") ;
# print("r=Fparallel/Fperpendicular= ") ; print($r) ; print(" \n") ;

$nm1_10to4th[$gas-1]=($ref-1)*10000.0 ;
$sig[$gas-1]=24.0* $pi**3 * $nu**4 * (($ref**2 - 1.0)/($ref**2 + 2.0))**2 * $Fk / $N0**2 ;
# print("Rayleigh Cross Section STP q= ") ; print($sig) ; print("\n") ;

$Fpar=(1.0+5.0*$gamma)/(4.0*(1.0+2.0*$gamma)) ;
$Bpar=$Fpar ;
$Fper=3.0*(1.0+$gamma)/(4.0*(1.0+2.0*$gamma)) ;
$Bper=$Fper ;

$betas=$N0*$sig[$gas-1] ;

$beta[$gas-1]=$betas*$fac*1.0e8 ;
$beta_par[$gas-1]=$beta[$gas-1]*$Fpar ;
$beta_per[$gas-1]=$beta[$gas-1]*$Fper ;

$Tau[$gas-1]=$conc[$gas-1]*$beta[$gas-1]*$T*$Rgas[$gas-1]/($g*1.0e6) ;

$Fkk[$gas-1] = $Fk ;
} # End the for loop

# $Tau[0]=sprintf("%1.4e",$Tau[0]) ; $Tau[1]=sprintf("%1.4e",$Tau[1]) ;

$Tau[2]=&round($Tau[2]) ;
foreach $sig (@sig) { $sig=$sig*1.0e27 ; $sig=&round($sig) ; }
foreach $Fkk (@Fkk) { $Fkk=&round($Fkk) ; }
foreach $beta (@beta) { $beta=&round($beta) ; }
foreach $beta_par (@beta_par) { $beta_par=&round($beta_par) ; }
foreach $beta_per (@beta_per) { $beta_per=&round($beta_per) ; }
foreach $nm1_10to4th (@nm1_10to4th) { $nm1_10to4th=&round($nm1_10to4th) ; }

@Tau ;
@Fkk ;
@sig ;
@beta ;
@beta_par ;
@beta_per ;
@nm1_10to4th ;
}

sub round {
my($number) = @_ ; # The incoming number(s) is(are) stored in the @_ variable.
my($c) = abs($number) ;
if ($c ge 1.0)
{
$number=$number*1000.000 ;
return( int($number + .5 * ($number <=> 0))/1000.0 );
} elsif ($c ge 0.1)
{
$number=$number*10000.000 ;
return( int($number + .5 * ($number <=> 0))/10000.0 );
} elsif ($c ge 0.01)
{
$number=$number*100000.000 ;
return( int($number + .5 * ($number <=> 0))/100000.0 );
} elsif ($c ge 0.001)
{
$number=$number*1000000.000 ;
return( int($number + .5 * ($number <=> 0))/1000000.0 );
}
else
{
$number=$number*10000000.000 ;
return( int($number + .5 * ($number <=> 0))/10000000.0 );
}
}

sub AirMassCalculator {

use Math::Trig;
use POSIX "fmod";

# Air Mass Calculator v1.0 - Input Date and Time and get Air Mass for Reno, NV.
# Written by Ben Sumlin.
# Adapted for the subroutine by Pat Arnott.

$month = $mon ; # This value is taken from the calling program.
$day = $daychosen ; # This value is taken from the calling program.
$year = $yearchosen ; # This value is taken from the calling program.
$minute = $min ; # This value is taken from the calling program.
$second = $sec ; # This value is taken from the calling program.
$hour = $hourchosen ; # This value is taken from the calling program.
$timezone = $timezonechosen ; # This value is taken from the calling program.
$Latitude = $lat ; # This value is taken from the calling program.
$Longitude = $lon ; # This value is taken from the calling program.

# $JDN = Delta_Days(1899,12,30,$year,$month,$day); #Should match excel sheet

# Code fragment to compute JDN.
# Input variable definitions: $y=year, $m=month, $day.
# Output variable = JDN.
my $y=$year ;
my $m=$month ;
$rem=$y%4.0 ; # The symbol % implies modulus taken here.
if ( $m==4 || $m==6 || $m==9 || $m==11 ) { $daymax=30.0 ; }
elsif ( $m==2 ) {
$daymax=28.0 ; if ($rem<0.1) { $daymax=29.0 ; }
}
else { $daymax=31.0 ; }
# work on determining if this is a leap year.
@d=(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0) ;
for($i = 1; $i <= 12; $i++) { $d[$i]=31.0 ; }
if ($rem>0.1) { $d[2]=28 ; }
else { $d[2]=29 ; }
$d[4]=30 ; $d[6]=30 ; $d[9]=30 ; $d[11]=30 ;
if ($m<1.1) { $daysum=0. }
else
{
$daysum=0.0 ;
for($i = 1; $i <= $m-1; $i++) { $daysum=$daysum+$d[$i] }
}
$JDN=($y-1900.0)*365.0 + $daysum + $day ;
$cor2=0.0 ;
if ( $rem>0.5 && $rem<1.5) { $cor2 = 1 ; }
if ( $rem>1.5 && $rem<2.5 && ($y-2)%8<0.1) { $cor2 = 1 ; }

if ( $rem<0.1 )
{ # process the leap year case.
if ($y >= 1904) { $cor = sprintf("%.0f",($y-1904.0)/4.00) + 1 ; }
else { $cor=0 ; }
}
else # Not a leap year.
{
if ($y >= 1904) { $cor = sprintf("%.0f",($y-1904.0)/4.00) + 1 + $cor2 ; }
else { $cor=0 ; }
}

$JDN=$JDN+$cor ;
# End code fragment used to compute $JDN.

$time = ($hour + ($minute + ($second / 60)) / 60) / 24;
$JD = $JDN + 2415018.5 + $time - ($timezone / 24); #Julian Day
$JC = ($JD - 2451545) / 36525; #Julian Century

$GeomMeanLongSun = fmod((280.46646 + $JC*(36000.76983 + $JC*0.0003032)),360);
$GeomMeanAnomSun = (357.52911+$JC*(35999.05029-(0.0001537 * $JC)));
$EccentEarthOrbit = (0.016708637-$JC*(0.000042037+0.0001537*$JC));
$SunEqOfCtr = (sin((deg2rad($GeomMeanAnomSun)))*(1.914602-$JC*(0.004817+0.000014*$JC)))+(sin(deg2rad(2*$GeomMeanAnomSun))*(0.019993-0.000101*$JC))+(sin(deg2rad(3*$GeomMeanAnomSun))*0.000289);
$SunTrueLong = $GeomMeanLongSun + $SunEqOfCtr;
$SunTrueAnom = $GeomMeanAnomSun + $SunEqOfCtr;
$SunRadVector = (1.0000010108*(1-$EccentEarthOrbit*$EccentEarthOrbit))/(1+$EccentEarthOrbit*cos(deg2rad($SunTrueAnom)));
$SunAppLong = $SunTrueLong-0.00569-0.00478*sin(deg2rad(125.04-1934.136*$JC));
$MeanObliqEcliptic = 23+(26+((21.448-$JC*(46.815+$JC*(0.00059-$JC*0.001813))))/60)/60;
$ObliqCorr = $MeanObliqEcliptic+0.00256*cos(deg2rad(125.04-1934.136*$JC));
$SunRtAscen = rad2deg(atan2(cos(deg2rad($ObliqCorr))*sin(deg2rad($SunAppLong)),cos(deg2rad($SunAppLong))));
$SunDeclin = rad2deg(asin(sin(deg2rad($ObliqCorr))*sin(deg2rad($SunAppLong))));
$vary = tan(deg2rad($ObliqCorr/2))*tan(deg2rad($ObliqCorr/2));
$EqOfTime = 4*rad2deg($vary*sin(2*deg2rad($GeomMeanLongSun))-2*$EccentEarthOrbit*sin(deg2rad($GeomMeanAnomSun))+4*$EccentEarthOrbit*$vary*sin(deg2rad($GeomMeanAnomSun))*cos(2*deg2rad($GeomMeanLongSun))-0.5*$vary*$vary*sin(4*deg2rad($GeomMeanLongSun))-1.25*$EccentEarthOrbit*$EccentEarthOrbit*sin(2*deg2rad($GeomMeanAnomSun)));
$HASunrise = rad2deg(acos(cos(deg2rad(90.833))/(cos(deg2rad($Latitude))*cos(deg2rad($SunDeclin)))-tan(deg2rad($Latitude))*tan(deg2rad($SunDeclin))));
$TrueSolarTime = fmod((($time*1440) + $EqOfTime + ($Longitude*4) - ($timezone*60)),1440);
$HourAngle = 0;
if (($TrueSolarTime/4) < 0)
{
$HourAngle = ($TrueSolarTime/4)+180;
}
else
{
$HourAngle = ($TrueSolarTime/4)-180;
}
$SolarZenithAngle = rad2deg(acos(sin(deg2rad($Latitude))*sin(deg2rad($SunDeclin))+cos(deg2rad($Latitude))*cos(deg2rad($SunDeclin))*cos(deg2rad($HourAngle))));
$CosineSolarZenithAngle = cos(deg2rad($SolarZenithAngle));
$AirMass = 1/$CosineSolarZenithAngle;
$AirMass=&round($AirMass) ;
$AirMassOut = $AirMass ;
$AirMassOut ;

}