perl/Sun.pm

   1 #!/usr/bin/perl -w
   2 #
   3 # The subroutines "julian_day" and "riseset" written by
   4 # Steve Franke, November 1999.
   5 #
   6 # The formulas used to calculate sunrise and sunset times
   7 # are described in Chapters 7, 12, 15, and 25 of
   8 # Astronomical Algorithms, Second Edition
   9 # by Jean Meeus, 1998
  10 # Published by Willmann-Bell, Inc.
  11 # P.O. Box 35025, Richmond, Virginia 23235
  12 #
  13 package Sun;
  14
  15 use POSIX;
  16
  17 require Exporter;
  18 @ISA = qw(Exporter);
  19 @EXPORT = qw($pi $d2r $r2d );
  20
  21 use strict;
  22 use vars qw($pi $d2r $r2d );
  23
  24 $pi = 3.141592653589;
  25 $d2r = ($pi/180);
  26 $r2d = (180/$pi);
  27
  28 sub julian_day
  29 {
  30         my $year = shift;
  31         my $month = shift;
  32         my $day = shift;
  33         my $julianday;
  34
  35         $year=$year-1 if( $month <= 2 );
  36         $month=$month+12 if( $month <= 2);
  37
  38         $julianday = int(365.25*($year+4716)+int(30.6001*($month+1)))+$day-13-1524.5;
  39         return $julianday;
  40 }
  41
  42 sub riseset
  43 {
  44         my $year = shift;
  45         my $month = shift;
  46         my $day = shift;
  47         my $lat = shift;
  48         my $lon = shift;
  49         my $julianday;
  50
  51         $julianday=julian_day($year,$month,$day);
  52
  53         my $tt = ($julianday-2451545)/36525.;
  54
  55         my $theta0=280.46061837+360.98564736629*($julianday-2451545.0)+
  56                 0.000387933*($tt^2)-($tt^3)/38710000;
  57         $theta0=$theta0-int($theta0/360)*360;
  58         $theta0=$theta0+360 if( $theta0 < 0 );
  59
  60         my $L0 = 280.46646+36000.76983*$tt+0.0003032*($tt^2);
  61         $L0=$L0-int($L0/360)*360;
  62         $L0=$L0+360 if( $L0 < 0 );
  63
  64         my $M = 357.52911 + 35999.05029*$tt-0.0001537*($tt^2);
  65         $M=$M-int($M/360)*360;
  66         $M=$M+360 if( $M < 0 );
  67
  68         my $C = (1.914602 - 0.004817*$tt-0.000014*($tt^2))*sin($M*$d2r) +
  69                 (0.019993 - 0.000101*$tt)*sin(2*$M*$d2r) +
  70                 0.000289*sin(3*$M*$d2r);
  71
  72         my $OMEGA = 125.04 - 1934.136*$tt;
  73
  74         my $lambda=$L0+$C-0.00569-0.00478*sin($OMEGA*$d2r);
  75
  76         my $epsilon = 23+26./60.+21.448/(60.*60.);
  77
  78         my $alpha=atan2(cos($epsilon*$d2r)*sin($lambda*$d2r),cos($lambda*$d2r))*$r2d;
  79         $alpha = $alpha-int($alpha/360)*360;
  80         $alpha=$alpha+360 if ( $alpha < 0 );
  81
  82         my $delta=asin(sin($epsilon*$d2r)*sin($lambda*$d2r))*$r2d;
  83         $delta = $delta-int($delta/360)*360;
  84         $delta = $delta+360 if ( $delta < 0 );
  85
  86         my $arg = (sin(-.8333*$d2r)-sin($lat)*sin($delta*$d2r))/(cos($lat)*cos($delta*$d2r));
  87         my $argtest = tan($lat)*tan($delta*$d2r);
  88
  89         if ( $argtest < -1. ) {
  90                 return sprintf("Sun doesn't rise.");
  91         }
  92         if ( $argtest > 1. ) {
  93                 return sprintf("Sun doesn't set.");
  94         }
  95
  96         my $H0 = acos($arg)*$r2d;
  97
  98         my $transit = ($alpha + $lon*$r2d - $theta0)/360.;
  99         $transit=$transit+1 if( $transit < 0 );
 100         $transit=$transit-1 if( $transit > 1 );
 101
 102         my $rise = $transit - $H0/360.;
 103         $rise=$rise+1 if( $rise < 0 );
 104         $rise=$rise-1 if( $rise > 1 );
 105
 106         my $set = $transit + $H0/360.;
 107         $set=$set+1 if( $set < 0 );
 108         $set=$set-1 if( $set > 1 );
 109
 110         return sprintf("Sunrise: %2.2d%2.2dZ   Sunset: %2.2d%2.2dZ",int($rise*24),
 111                                    ($rise*24-int($rise*24))*60.,
 112                                    int($set*24),($set*24-int($set*24))*60.);
 113 }
 114
 115 1;