#!/usr/local/bin/perl

=head1

	sunrise.pm	Computes sun rise/set times, at any date
			and latitude, giving sunset times based on 0h UT.

	Translated into perl by Michael Dabrowski, from SUNRISE.C
	(c) Paul Schlyter, 1989, 1992.

	This code is Copyright (c) by www.adventist.org, All Rights
	reserved. You may freely use and distribute this library, but the
	editing of this code is forbidden. We do not give any warranty or
	guarentee as to the accuracy of the results, or the function of
	this code. This code is presented AS IS. 

	Note: This library does not account for Time Zones, nor daylight
	      savings time.

	Usage: sunriset($year, $month, $day, $latitude, $logitude);
=cut1


use POSIX;

sub sr_main
{
	# Integers used

	local($year);
	local($month);
	local($day);
	local($rs);
	local($lon); 
	local($lat);
	local($rise, $set);

	print "Longitude (+ is east) and latitude (+ is north) : ";

	$_ = <STDIN>; chomp; ($lon, $lat) = split(/ /);	

	print "Input date: ( yyyy mm dd) : ";

	$_ = <STDIN>; chomp; ($year, $month, $day) = split(/ /);

	($rs, $rise, $set) = sunriset($year, $month, $day, $lon, $lat);

	if($rs == 0)  { print "Sun rises $rise UT, sets $set UT\n"; }
	elsif($rs == +1) { print "Sun above horizon\n"; }
	elsif($rs == -1) { print "Sun below horizon\n"; }

	return $rs, $rise, $set;
}

sub sunriset
{
	# sunriset($year, $month, $day, $latitude, $logitude);
	#
	#	$year in 4 digits.
	#
	#	Return values:
	#		$rs	1 	Sun above horizon
	#		$rs	-1	Sun below horizon
	#		$rise		float
	#		$set		float

	# GLOBAL variable definitions

	local($pi) = 3.1415926535897932384;
	local($radeg) = 180.0 / $pi;		
	local($degrad) = $pi / 180.0;
	local($inv360) = 1.0 / 360.0;

	# Passed variables

	local($year) 		= $_[0]; local($month) 	= $_[1];
	local($day) 		= $_[2]; local($lon) 	= $_[3];
	local($lat) 		= $_[4]; local($altit) 	= -35.0/60.0;
	local($upper_limb) 	= 1;

	# Return Values:
		local($rc)    = 0;
		local($trise) = $_[7];
		local($tset)  = $_[8];


	local($d);	# Days since 2000 Jan 0.0
	local($sr);	# Solar distance
	local($sra);	# Sun's right ascension
	local($sdec);	# Sun's declination
	local($sradius);# Sun's apparent radius
	local($t);	# diurnam arc;
	local($tsouth); # time when Sun is at south
	local($sidtime);# local sidereal time

	$d = days_since_2000_Jan_0($year, $month, $day) + 0.5 - $lon/360.0;

	$sidtime = revolution(GMST0($d) + 180.0 + $lon);

	($sra, $sdec, $sr) = sun_RA_dec($d, $sr);

	$tsouth = 12.0 - rev180($sidtime - $sra) / 15.0;

	$sradius = 0.2666 / $sr;

	if($upper_limb)
	{ $altit -= $sradius; }

	$cost = (sind($altit) - sind($lat) * sind($sdec)) / (cosd($lat) * cosd($sdec));
	if($cost >= 1.0)
	{ $rc = -1, $t = 0.0; }
	elsif($cost <= -1.0)
	{ $rc = +1; $t = 12.0; }
	else { $t = acosd($cost) / 15.0; }

	$trise = $tsouth - $t;
	$tset = $tsouth + $t;

	return $rc, $trise, $tset;
}

sub sunpos
{
	## Input variables

	local($d) = $_[0];
	
	## Return variables

	local($lon);
	local($r);

	## Local declerations

	local($M, $w, $e, $E, $x, $y, $v);

	$M = revolution(356.0470 + 0.9856002585 * $d);
	$w = 282.9404 + 4.70935E-5 * $d;
	$e = 0.016709 - 1.151E-9 * $d;

	$E = $M + $e * $radeg * sind($M) * (1.0 + $e * cosd($M));
	$x = cosd($E) - $e;
	$y = sqrt(1.0 - $e * $e) * sind($E);
	$r = sqrt($x * $x + $y * $y); 	
	$v = atan2d($y, $x);
	$lon = $v + $w;

	if( $lon >= 360.0 ) { $lon -= 360.0; }

	return $lon, $r;
}

sub sun_RA_dec
{
	## Input variables

	local($d) = $_[0];

	## Return Variables

	local($RA, $dec);
	local($r) = $_[1];

	## Local reserved variables

	local($lon, $obl_ecl, $x, $y, $z);

	($lon, $r) = sunpos($d);

	$x = $r * cosd($lon);
	$y = $r * sind($lon);

	$obl_ecl = 23.4393 - 3.563E-7 * $d;

	$z = $y * sind($obl_ecl);
	$y = $y * cosd($obl_ecl);

	$RA = atan2d($y, $x);
	$dec = atan2d($z, sqrt($x * $x + $y * $y));

	return $RA, $dec, $r;
}

sub revolution
{ return ($_[0] - 360.0 * floor($_[0] * $inv360)); }

sub rev180
{ return ($_[0] - 360.0 * floor($_[0] * $inv360 + 0.5)); }

sub GMST0
{
 return revolution((180.0 + 356.0470 + 282.9404) + 
		    (0.9856002585 + 4.70935E-5) * $_[0]);
}

sub days_since_2000_Jan_0
{
	my($y) = $_[0]; my($m) = $_[1]; my($d) = $_[2];
	return (367 * ($y) - ((7 * (($y) + ((($m) + 9) / 12))) /4) + 
	((275 * ($m)) / 9) + ($d) - 730530);
}

sub sind { return sin($_[0] * $degrad); }
sub cosd { return cos($_[0] * $degrad); }
sub tand { return tan($_[0] * $degrad); }
sub atand { return ($radeg * atan($_[0])); }
sub asind { return ($radeg * asin($_[0])); }
sub acosd { return ($radeg * acos($_[0])); }
sub atan2d { return ($radeg * atan2($_[0], $_[1])); }

1;