#!/usr/bin/perl
use strict;
use Math::Trig;
# References:
# (a) Michalsky, J. J., 1988, The Astronomical Almanac's algorithm for
# approximate solar position (1950-2050), Solar Energy, 227---235, 1988
#
# (b) Spencer, J. W., 1989, Comments on The Astronomical
# Almanac's algorithm for approximate solar position (1950-2050)
# Solar Energy, 42, 353

# Input:
# year - the year number (e.g. 1977)
# day  - the day number of the year starting with 1 for
#        January 1
# time - decimal time. E.g. 22.89 (8.30am eastern daylight time is
#        equal to 8.5+5(hours west of Greenwich) -1 (for daylight savings
#        time correction
# lat -  local latitude in degrees (north is positive)
# lon -  local longitude (east of Greenwich is positive.
#                         i.e. Honolulu is 15.3, -157.8)
#
# Output:
# a - azimuth angle of the sun (measured east from north 0 to 360)
# e - elevation of the sun
#
# Spencer correction introduced and 3 lines of Michalsky code
# commented out (after calculation of az)
#
my $mxtim=100;
my ($cozena, $dlat, $dlon, $dday, $dtime, $val, $hour);
my ($musun,$a,$e,$az,$el,$refrac,$time); 

my $year=2007;
my $day=65;
my $hour=10.5;
my $lon=22.0;
my $lat=52.0;
 
my $pi=4.*atan(1.0);
my $twopi=2.*$pi;
my $rad=$pi/180.0;
# get the current Julian date
my $delta=$year-1949.;
my $leap=int($delta/4.);
my $jd=32916.5+$delta*365.+$leap+$day+$hour/24.;
# calculate ecliptic coordinates
$time=$jd-51545.0;
print $time,"TIME\n";
# force mean longitude between 0 and 360 degs
my $mnlong=280.460+0.9856474*$time;
$mnlong=mod($mnlong,360.0);  
if($mnlong<0) {$mnlong=$mnlong+360};
# mean anomaly in radians between 0, 2*pi
my $mnanom=357.528+0.9856003*$time;
$mnanom=mod($mnanom,360);
if ($mnanom<0) {$mnanom=$mnanom+360};
$mnanom=$mnanom*$rad;
print $mnlong,"MNLONG\n";
print $mnanom,"MNANOM\n";
# compute ecliptic longitude and obliquity of ecliptic
#       eclong=mnlong+1.915*(mnanom)+0.20*sin(2.*mnanom)
my $eclong=$mnlong+1.915*sin($mnanom)+0.020*sin(2.0*$mnanom);
$eclong=mod($eclong,360);
if ($eclong<0) {$eclong=$eclong+360}
my $oblqec=23.429-0.0000004*$time;
$eclong=$eclong*$rad;
$oblqec=$oblqec*$rad;
print $eclong,"ECLONG\n";
print $oblqec,"OBLEQEC\n";
# calculate right ascention and declination
my $num=cos($oblqec)*sin($eclong);
my $den=cos($eclong);
my $ra=atan($num/$den);
print $num,"NUM\n";
print $den,"DEN\n";
print $ra,"RA\n";
# force ra between 0 and 2*pi
if($den<0) {$ra=$ra+$pi}
 else {
  if ($num<0) {$ra=$ra+$twopi};
      }
print $num,"\n";
print $den,"\n";
print $ra,"\n";      
# dec in radians
my $dec=asin(sin($oblqec)*sin($eclong));
# calculate Greenwich mean sidereal time in hours
my $gmst=6.697375+0.0657098242*$time+$hour;
print $time,$hour,"AAA\n";
# hour not changed to sidereal sine "time" includes the fractional day
$gmst=mod($gmst,24);
if($gmst<0) {$gmst=$gmst+24};
# calculate local mean sidereal time in radians
print $gmst,"\n";
my $lmst=$gmst+$lon/15;
print $lmst,"LMST\n";
$lmst=mod($lmst,24);
print $lmst,"LMST\n";
if($lmst<0) {$lmst=$lmst+24};
$lmst=$lmst*15*$rad;
# calculate hour angle in radians between -pi, pi
my $ha=$lmst-$ra;

print $ha,"HA\n";
print $lmst,"LMST\n";
print $ra,"RA\n";
if($ha<-$pi) {$ha=$ha+$twopi};
if($ha>$pi) {$ha=$ha-$twopi};
print $ha,"HA\n";
$lat=$lat*$rad;
# calculate azimuth and elevation
$el=asin(sin($dec)*sin($lat)+cos($dec)*cos($lat)*cos($ha));
$az=asin(-cos($dec)*sin($ha)/cos($el));
print $el,"EL\n";
print $az,"AZ\n";
# add J. W. Spencer code (next 5 lines)
if(sin($dec)-sin($el)*sin($lat)>=0) {
 if(sin($az)<0) {$az=$az+$twopi}
}
if(sin($dec)-sin($el)*sin($lat)<0)  
{$az=$pi-$az}
 
  
# end Spencer's corrections

# this puts azimuth between 0 and 2*pi radians
# calculate refraction correction for US stand. atm.
$el=$el/$rad;
if($el>-0.56)
 {$refrac=3.51561*(0.1594+0.0196*$el+0.00002*$el**2)/(1.+0.505*$el+0.0845*$el**2)}
 else {$refrac=0.56};
$el=$el+$refrac;
$az=$az/$rad;
$lat=$lat/$rad;

print $refrac,"\n";
print $el,"\n";
print $az,"\n";

sub mod
{
my $n;
my ($x,$y)=@_;
$n=int($x/$y);
my $res=$x-$n*$y;
return $res;
}