Astro::Sunrise - Perl extension for computing the sunrise/sunset on a given day

 use Astro::Sunrise;
#use Astro::Sunrise qw(:constants);

 ($sunrise, $sunset) = sunrise(YYYY,MM,DD,longitude,latitude,Time Zone,DST);
 ($sunrise, $sunset) = sunrise(YYYY,MM,DD,longitude,latitude,Time Zone,DST,ALT);
 ($sunrise, $sunset) = sunrise(YYYY,MM,DD,longitude,latitude,Time Zone,DST,ALT,inter);

 $sunrise = sun_rise(longitude,latitude);
 $sunset = sun_set(longitude,latitude);

 $sunrise = sun_rise(longitude,latitude,ALT);
 $sunset = sun_set(longitude,latitude,ALT);

 $sunrise = sun_rise(longitude,latitude,ALT,day_offset);
 $sunset = sun_set(longitude,latitude,ALT,day_offset);

 Eastern longitude is entered as a positive number
 Western longitude is entered as a negative number
 Northern latitude is entered as a positive number
 Southern latitude is entered as a negative number

package Astro::Sunrise;

use strict;
#use warnings;
use POSIX qw(floor);
use Math::Trig;
use Carp;
use DateTime;
use vars qw( $VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $RADEG $DEGRAD );

require Exporter;

@ISA       = qw( Exporter );
@EXPORT    = qw( sunrise sun_rise sun_set );
@EXPORT_OK = qw( DEFAULT CIVIL NAUTICAL AMATEUR ASTRONOMICAL );
%EXPORT_TAGS = ( 
	constants => [ @EXPORT_OK ],
	);
	
$VERSION =  '0.91';
$RADEG   = ( 180 / pi );
$DEGRAD  = ( pi / 180 );
my $INV360     = ( 1.0 / 360.0 );

my $upper_limb = '1';

sub sunrise  {
my ( $year, $month, $day, $lon, $lat, $TZ, $isdst, $alt, $iter ) = @_;
   my $altit      = $alt || -0.833;
   my $iteration = defined($iter) ? $iter:0 ;
   
   if ($iteration)   {
   # This is the initial start

   my $d = days_since_2000_Jan_0( $year, $month, $day ) + 0.5 - $lon / 360.0;
   my ($tmp_rise_1,$tmp_set_1) = sun_rise_set($d, $lon, $lat,$altit,15.04107);

   # Now we have the initial rise/set times next recompute d using the exact moment
   # recompute sunrise
   
   my $tmp_rise_2=9;
   my $tmp_rise_3 = 0;
   until (equal($tmp_rise_2, $tmp_rise_3, 8) )   {

         my $d_sunrise_1 = $d + $tmp_rise_1/24.0;
         ($tmp_rise_2,undef) = sun_rise_set($d_sunrise_1, $lon, $lat,$altit,15.04107);
         $tmp_rise_1 = $tmp_rise_3;
         my $d_sunrise_2 = $d + $tmp_rise_2/24.0;
         ($tmp_rise_3,undef) = sun_rise_set($d_sunrise_2, $lon, $lat,$altit,15.04107);
       
         #print "tmp_rise2 is: $tmp_rise_2 tmp_rise_3 is:$tmp_rise_3\n";
         
   }

   
#######################################################################################
# end sunrise
###################################################################################


my $tmp_set_2=9;
my $tmp_set_3=0;

   until (equal($tmp_set_2, $tmp_set_3, 8) )   {

         my $d_sunset_1 = $d + $tmp_set_1/24.0;
         (undef,$tmp_set_2) = sun_rise_set($d_sunset_1, $lon, $lat,$altit,15.04107);
         $tmp_set_1 = $tmp_set_3;
         my $d_sunset_2 = $d + $tmp_set_2/24.0;
         (undef,$tmp_set_3) = sun_rise_set($d_sunset_2, $lon, $lat,$altit,15.04107);
        
         #print "tmp_set_1 is: $tmp_set_1 tmp_set_3 is:$tmp_set_3\n";
         
   }
   
   
   return convert_hour($tmp_rise_3,$tmp_set_3,$TZ, $isdst);

   }else{
   my $d = days_since_2000_Jan_0( $year, $month, $day ) + 0.5 - $lon / 360.0;
   my ($h1,$h2) = sun_rise_set($d, $lon, $lat,$altit,15.0);
   return convert_hour($h1,$h2,$TZ, $isdst);
   }
}


sub sun_rise_set {
    my ($d, $lon, $lat,$altit) =@_;
    #my ( $year, $month, $day, $lon, $lat, $TZ, $isdst, $alt ) = @_;
    #my $altit      = $alt || -0.833;
    #my $d = days_since_2000_Jan_0( $year, $month, $day ) + 0.5 - $lon / 360.0;

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

    my ( $sRA, $sdec ) = sun_RA_dec($d);
    my $tsouth  = 12.0 - rev180( $sidtime - $sRA ) / 15.0;
    my $sradius = 0.2666 / $sRA;

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

    # Compute the diurnal arc that the Sun traverses to reach 
    # the specified altitude altit: 

    my $cost =
      ( sind($altit) - sind($lat) * sind($sdec) ) /
      ( cosd($lat) * cosd($sdec) );

    my $t;
    if ( $cost >= 1.0 ) {
        carp "Sun never rises!!\n";
        $t = 0.0;    # Sun always below altit
    }
    elsif ( $cost <= -1.0 ) {
        carp "Sun never sets!!\n";
        $t = 12.0;    # Sun always above altit
    }
    else {
        $t = acosd($cost) / 15.0;    # The diurnal arc, hours
    }

    # Store rise and set times - in hours UT 

    my $hour_rise_ut = $tsouth - $t;
    my $hour_set_ut  = $tsouth + $t;
    return($hour_rise_ut, $hour_set_ut);
    #return convert_hour($hour_rise_ut,$hour_set_ut,$TZ, $isdst);
}

#########################################################################################################
sub GMST0 {
#
#
# FUNCTIONAL SEQUENCE for GMST0 
#
# _GIVEN
# Day number
#
# _THEN
#
# computes GMST0, the Greenwich Mean Sidereal Time  
# at 0h UT (i.e. the sidereal time at the Greenwhich meridian at  
# 0h UT).  GMST is then the sidereal time at Greenwich at any     
# time of the day..
# 
#
# _RETURN
#
# Sidtime
#
    my ($d) = @_;

    my $sidtim0 =
      revolution( ( 180.0 + 356.0470 + 282.9404 ) +
      ( 0.9856002585 + 4.70935E-5 ) * $d );
    return $sidtim0;

}

sub sunpos {

#
#
# FUNCTIONAL SEQUENCE for sunpos
#
# _GIVEN
#  day number
#
# _THEN
#
# Computes the Sun's ecliptic longitude and distance */
# at an instant given in d, number of days since     */
# 2000 Jan 0.0. 
# 
#
# _RETURN
#
# ecliptic longitude and distance
# ie. $True_solar_longitude, $Solar_distance
#
    my ($d) = @_;

    #                       Mean anomaly of the Sun 
    #                       Mean longitude of perihelion 
    #                         Note: Sun's mean longitude = M + w 
    #                       Eccentricity of Earth's orbit 
    #                       Eccentric anomaly 
    #                       x, y coordinates in orbit 
    #                       True anomaly 

    # Compute mean elements 
    my $Mean_anomaly_of_sun = revolution( 356.0470 + 0.9856002585 * $d );
    my $Mean_longitude_of_perihelion = 282.9404 + 4.70935E-5 * $d;
    my $Eccentricity_of_Earth_orbit  = 0.016709 - 1.151E-9 * $d;

    # Compute true longitude and radius vector 
    my $Eccentric_anomaly =
      $Mean_anomaly_of_sun + $Eccentricity_of_Earth_orbit * $RADEG *
      sind($Mean_anomaly_of_sun) *
      ( 1.0 + $Eccentricity_of_Earth_orbit * cosd($Mean_anomaly_of_sun) );

    my $x = cosd($Eccentric_anomaly) - $Eccentricity_of_Earth_orbit;

    my $y =
      sqrt( 1.0 - $Eccentricity_of_Earth_orbit * $Eccentricity_of_Earth_orbit )
      * sind($Eccentric_anomaly);

    my $Solar_distance = sqrt( $x * $x + $y * $y );    # Solar distance
    my $True_anomaly = atan2d( $y, $x );               # True anomaly

    my $True_solar_longitude =
      $True_anomaly + $Mean_longitude_of_perihelion;    # True solar longitude

    if ( $True_solar_longitude >= 360.0 ) {
        $True_solar_longitude -= 360.0;    # Make it 0..360 degrees
    }

    return ( $Solar_distance, $True_solar_longitude );
}

sub sun_RA_dec {

#
#
# FUNCTIONAL SEQUENCE for sun_RA_dec 
#
# _GIVEN
# day number, $r and $lon (from sunpos) 
#
# _THEN
#
# compute RA and dec
# 
#
# _RETURN
#
# Sun's Right Ascension (RA) and Declination (dec)
# 
#
    my ($d) = @_;

    # Compute Sun's ecliptical coordinates 
    my ( $r, $lon ) = sunpos($d);

    # Compute ecliptic rectangular coordinates (z=0) 
    my $x = $r * cosd($lon);
    my $y = $r * sind($lon);

    # Compute obliquity of ecliptic (inclination of Earth's axis) 
    my $obl_ecl = 23.4393 - 3.563E-7 * $d;

    # Convert to equatorial rectangular coordinates - x is unchanged 
    my $z = $y * sind($obl_ecl);
    $y = $y * cosd($obl_ecl);

    # Convert to spherical coordinates 
    my $RA  = atan2d( $y, $x );
    my $dec = atan2d( $z, sqrt( $x * $x + $y * $y ) );

    return ( $RA, $dec );

}    # sun_RA_dec

sub days_since_2000_Jan_0 {

#
#
# FUNCTIONAL SEQUENCE for days_since_2000_Jan_0 
#
# _GIVEN
# year, month, day
#
# _THEN
#
# process the year month and day (counted in days)
# Day 0.0 is at Jan 1 2000 0.0 UT
# Note that ALL divisions here should be INTEGER divisions
#
# _RETURN
#
# day number
#
    use integer;
    my ( $year, $month, $day ) = @_;

    my $d =
      ( 367 * ($year) -
      int( ( 7 * ( ($year) + ( ( ($month) + 9 ) / 12 ) ) ) / 4 ) +
      int( ( 275 * ($month) ) / 9 ) + ($day) - 730530 );

    return $d;

}

sub sind {
    sin( ( $_[0] ) * $DEGRAD );
}

sub cosd {
    cos( ( $_[0] ) * $DEGRAD );
}

sub tand {
    tan( ( $_[0] ) * $DEGRAD );
}

sub atand {
    ( $RADEG * atan( $_[0] ) );
}

sub asind {
    ( $RADEG * asin( $_[0] ) );
}

sub acosd {
    ( $RADEG * acos( $_[0] ) );
}

sub atan2d {
    ( $RADEG * atan2( $_[0], $_[1] ) );
}

sub revolution {
#
#
# FUNCTIONAL SEQUENCE for revolution
#
# _GIVEN
# any angle
#
# _THEN
#
# reduces any angle to within the first revolution 
# by subtracting or adding even multiples of 360.0
# 
#
# _RETURN
#
# the value of the input is >= 0.0 and < 360.0
#

    my $x = $_[0];
    return ( $x - 360.0 * floor( $x * $INV360 ) );
}

sub rev180 {
#
#
# FUNCTIONAL SEQUENCE for rev180
#
# _GIVEN
# 
# any angle
#
# _THEN
#
# Reduce input to within +180..+180 degrees
# 
#
# _RETURN
#
# angle that was reduced
#
    my ($x) = @_;
    
    return ( $x - 360.0 * floor( $x * $INV360 + 0.5 ) );
}

sub equal {
    my ($A, $B, $dp) = @_;

    return sprintf("%.${dp}g", $A) eq sprintf("%.${dp}g", $B);
  }


sub convert_hour   {

#
#
# FUNCTIONAL SEQUENCE for convert_hour 
#
# _GIVEN
# Hour_rise, Hour_set, Time zone offset, DST setting
# hours are in UT
#
# _THEN
#
# convert to local time
# 
#
# _RETURN
#
# hour:min rise and set 
#

  my ($hour_rise_ut, $hour_set_ut, $TZ, $isdst) = @_;

  my $rise_local = $hour_rise_ut + $TZ;
  my $set_local = $hour_set_ut + $TZ;
  if ($isdst) {
    $rise_local +=1;
    $set_local +=1;
  }

  # Rise and set should be between 0 and 24;
  if ($rise_local<0) {
    $rise_local+=24;
  } elsif ($rise_local>24) {
    $rise_local -=24;
  }
  if ($set_local<0) {
    $set_local+=24;
  } elsif ($set_local>24) {
    $set_local -=24;
  }

  my $hour_rise =  int ($rise_local);
  my $hour_set  =  int($set_local);

  my $min_rise  = floor(($rise_local-$hour_rise)*60+0.5);
  my $min_set   = floor(($set_local-$hour_set)*60+0.5);

  if ($min_rise>=60) {
    $min_rise -=60;
    $hour_rise+=1;
    $hour_rise-=24 if ($hour_rise>=24);
  }
  if ($min_set>=60) {
    $min_set -=60;
    $hour_set+=1;
    $hour_set-=24 if ($hour_set>=24);
  }

  if ( $min_rise < 10 ) {
    $min_rise = sprintf( "%02d", $min_rise );
  }
  if ( $min_set < 10 ) {
    $min_set = sprintf( "%02d", $min_set );
  }
  $hour_rise = sprintf( "%02d", $hour_rise );
  $hour_set  = sprintf( "%02d", $hour_set );
  return ( "$hour_rise:$min_rise", "$hour_set:$min_set" );

}

sub sun_rise
   {
   my $longitude = shift;
   my $latitude = shift;
   my $alt = shift || -0.833;
   my $offset = int( shift || 0 );

   my $today = DateTime->today->set_time_zone( 'local' );
   $today->add( days => $offset );

   my( $sun_rise, undef ) = sunrise( $today->year, $today->mon, $today->mday,
                                     $longitude, $latitude,
                                     ( $today->offset / 3600 ),
                                     #
                                     # DST is always 0 because DateTime
                                     # currently (v 0.16) adds one to the
                                     # offset during DST hours
                                     0,
                                     $alt );
   return $sun_rise;
   }

sub sun_set
   {
   my $longitude = shift;
   my $latitude = shift;
   my $alt = shift || -0.833;
   my $offset = int( shift || 0 );

   my $today = DateTime->today->set_time_zone( 'local' );
   $today->add( days => $offset );

   my( undef, $sun_set ) = sunrise( $today->year, $today->mon, $today->mday,
                                    $longitude, $latitude,
                                    ( $today->offset / 3600 ),
                                    #
                                    # DST is always 0 because DateTime
                                    # currently (v 0.16) adds one to the
                                    # offset during DST hours
                                    0,
                                    $alt );
   return $sun_set;
   }

sub DEFAULT      () { -0.833 }
sub CIVIL        () { - 6 }
sub NAUTICAL     () { -12 }
sub AMATEUR      () { -15 }
sub ASTRONOMICAL () { -18 }

1;