use POSIX qw(floor);
use Math::Trig;
use constant year => 365.2424; # vernal equinox year in natural days
use constant obliquity => 23.4397/180*pi; # equator-orbit angle
use constant DT => 1E-5; # equator-orbit angle
sub addToDate{
	my $date=shift(@_);
	my $days=shift(@_);
	my @monthDays=(0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31);
	my ($year, $month, $day)=$date=~/(\d+)\-(\d+)\-(\d+)/;

	while($days>0){
		$monthDays[2]=( (!($year%4)) && (($year%100)||!($year%400)) ? 29:28); # this is a leap year
		$days--; $day++;
		if($day>$monthDays[$month]){
			$day=1; $month++;
			if($month>12){
				$month=1; $year++;
			}
		}
	}
	while($days<0){
		$monthDays[2]=( (!($year%4)) && (($year%100)||!($year%400)) ? 29:28); # this is a leap year
		$days++; $day--;
		if($day<1){
			$month--;
			if($month<1){
				$month=12; $year--;
			}
			$day=$monthDays[$month];
		}
	}
	return sprintf("%04s-%02s-%02s",$year, $month, $day);
}
sub dateTime2MJD{
	# convert a string like 2005-09-14 17:30:46.12 to MJD
	my @monthDays=(31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31);
	my $dateTime=shift(@_);
	my ($year, $month, $day, $hour, $minute, $second) = $dateTime =~ /^(\d{4})\D(\d{1,2})\D(\d{1,2})\D(\d{1,2})\D(\d{1,2})\D([\d\.]+)/;
	($year, $month, $day, $hour, $minute, $second) = $dateTime =~ /^(\d{4})\D(\d{1,2})\D(\d{1,2})\D(\d{1,2})\D(\d{1,2})/ if($year==0);
	($year, $month, $day) = $dateTime =~ /^(\d{4})\D(\d{1,2})\D(\d{1,2})/ if($year==0);
	return 0 if($year==0); # couldn't parse date

	# first the integer part
	my $i=1859; # year after epoch of MJD
	my $MJD=45; # MJD of 1859-1-1T00:00:00
	while($i<$year){ $MJD+=($i%400==0?366:$i%100==0?365:$i%4==0?366:365); $i++; }
	$monthDays[1] = ($year%400==0?29:$year%100==0?28:$year%4==0?29:28);
	for($i=0; $i<$month-1; $i++) { $MJD+=$monthDays[$i]; }
	$MJD += $day-1; 
	# now the real (non integer) part
	$MJD += ($hour+($minute+$second/60)/60)/24;
	return $MJD;
}
sub MJD2dateTime{
	# convert MJD to a string like 2005-09-14 17:30:46.12
	my @monthDays=(31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31);
	my $MJD = shift(@_);
	my ($year, $month, $day, $hour, $minute, $second);

	# first the integer part
	my $year=1859; # year after epoch of MJD
	$MJD-=45; # MJD of 1859-1-1
	while($MJD>($year%400==0?366:$year%100==0?365:$year%4==0?366:365)){ $MJD-=$year%400==0?366:$year%100==0?365:$year%4==0?366:365; $year++;}
	my $month=0;
	$monthDays[1] = ($year%400==0?29:$year%100==0?28:$year%4==0?29:28);
	while($MJD>$monthDays[$month]){ $MJD-= $monthDays[$month]; $month++;} $month++; #first month = number 1
	my $day=int($MJD); $MJD-=$day; $day++; # first day is number 1

	# now the real (non integer) part
	$MJD*=24;
	my $hour=int($MJD+DT); $MJD-=$hour;
	$MJD*=60;
	my $minute=int($MJD+DT); $MJD-=$minute;
	my $second=$MJD*60;

	return sprintf("%0004d-%02d-%02d %02d:%02d:%05.2f", $year, $month, $day, $hour, $minute, abs($second));
	return $MJD;
}
sub MJD2date{
	# convert MJD to a string like 2005-09-14 17:30:46.12
	my @monthDays=(31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31);
	my $MJD = shift(@_)+.1;
	my ($year, $month, $day, $hour, $minute, $second);

	# first the integer part
	my $year=1859; # year after epoch of MJD
	$MJD-=45; # MJD of 1859-1-1
	while($MJD>($year%400==0?366:$year%100==0?365:$year%4==0?366:365)){ $MJD-=$year%400==0?366:$year%100==0?365:$year%4==0?366:365; $year++;}
	my $month=0;
	$monthDays[1] = ($year%400==0?29:$year%100==0?28:$year%4==0?29:28);
	while($MJD>$monthDays[$month]){ $MJD-= $monthDays[$month]; $month++;} $month++; #first month = number 1
	my $day=floor($MJD); $MJD-=$day; $day++; # first day is number 1

	return sprintf("%0004d-%02d-%02d", $year, $month, $day);
}
sub MJD2time{
	# convert MJD to a string like 17:30:46.12
	my $MJD = shift(@_);
	$MJD-=floor($MJD);
	my ($hour, $minute, $second);

	# now the real (non integer) part
	$MJD*=24;
	my $hour=floor($MJD); $MJD-=$hour;
	$MJD*=60;
	my $minute=floor($MJD); $MJD-=$minute;
	my $second=$MJD*60;

	return sprintf("%02d:%02d:%05.2f", $hour, $minute, $second);
}
sub nowDateTime{
	# current UTC in DATETIME
	my ($second,$minute,$hour,$day,$month,$year,$wday,$yday,$isdst) = gmtime(time);
	$month++; $year+=1900;
	return sprintf("%0004d-%02d-%02d %02d:%02d:%05.2f", $year, $month, $day, $hour, $minute, $second);
}
sub nowMJD{
	return dateTime2MJD(nowDateTime());
}
sub solarAltAzimuth{
	my $long=shift(@_);
	my $lat=shift(@_);
	my $t=shift(@_);
	#print "long:$long, lat:$lat, t:$t\n";
	my $vernal = dateTime2MJD("2005-03-20T12:33:00"); #UTC, vernal equinox
	my $tf=$t-.5; $tf-=int($tf); # part of day since mid day
	my $yf=($t-$vernal)/year; $yf-=int($yf); # part of year vernal eq.
	my $zenLong=-2*pi*$tf; # longitude of location with sun in zenith
	my $zenLat=obliquity*sin(2*pi*$yf); # longitude of location with sun in zenith
	# at location (zenLong, zenLat) sun is at (east, north, up)=(0,0,1) in A.U.
	#     rotation along negative x-axis, a=zenLat
	#   ( 1     0      0    )
	#   ( 0   cos(a) sin(a) )
	#   ( 0  -sin(a) cos(a) )
	# (zenLong,0) => (e,n,u) = (0, sin(zenLat), cos(zenLat))
	#
	#     rotation along y-axis, a=zenLong-long
	#   ( cos(a)  0   sin(a) )
	#   (   0     1     0    )
	#   (-sin(a)  0   cos(a) )
	# (long, 0) => (e,n,u) = (cos(zenLat)sin(zenLong-long),
	#                 sin(zenLat), cos(zenLat)cos(zenLong-long))
	#
	#     rotation along negative x-axis, a=-lat
	#   ( 1     0      0    )
	#   ( 0   cos(a) sin(a) )
	#   ( 0  -sin(a) cos(a) )
	# (long, lat) => (e,n,u) = ( cos(zenLat)sin(zenLong-long),
	#                  sin(zenLat)cos(lat)-cos(zenLat)cos(zenLong-long)sin(lat),
	#                  sin(zenLat)sin(lat)+cos(zenLat)cos(zenLong-long)cos(lat) )
	my $east=cos($zenLat)*sin($zenLong-$long);
	my $north=sin($zenLat)*cos($lat)-cos($zenLat)*cos($zenLong-$long)*sin($lat);
	my $up=sin($zenLat)*sin($lat)+cos($zenLat)*cos($zenLong-$long)*cos($lat);
	my $alt=atan2($up, sqrt($north**2+$east**2));
	my $azi=atan2($east, $north); $azi+=2*pi if($azi<0); # north=0, clockwise
	return ($alt, $azi, $yf);
}

1