Re: calculate length of day as function of space at onset of fall

[Linicks]

I found some more code that calculates sunrise/sunset. After a lot of messing about, I found it is pretty accurate (it's a day out this year though, as it doesn't understand leap years). In the code below, search for 'Nick' as I have put in a few comments of how to set it up for your lat/long/timezone etc.

#!/usr/bin/perl

use POSIX;
sub sr_main
{
    ($rs, $rise, $set) = sunriset($year, $month, $day, $lon, $lat);
    return $rs, $rise, $set;
}
sub sunriset
{
    # 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/60;
        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/3
+60.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;
#($rs,$sunrise,$sunset)=sunriset($year, $month, $day, $latitude, $logi
+tude);
#print "Sunrise: ".$sunrise."n";

$year=2016;
$month=1;
#$day=11;

# Nick - 29/10/2016
# Put your lat and long here:
$lat="+50.8198";
$lon="-1.0880";

# From the origianl code
#$lat="+43.6667";
#$lon="-79.4";
# End Nick

$counter=1;
while ($counter < 366)
{
$day=$counter;
        ($rs, $rise, $set)=sr_main;
#print $rs;
#print $rise;
$prettycounter=0;
$lines=(24-($set-$rise));
for ($lines; $lines >= 0; $lines--) {
print "-";
$prettycounter++;
}
$stars=($set-$rise)*2;
for ($stars; $stars >= 0; $stars--) {
print "*";
$prettycounter++;
}
$lines=(24-($set-$rise));
for ($lines; $lines >= 0; $lines--) {
print "-";
$prettycounter++;
}
if ($prettycounter<50) { print "-"; };

# Nick - 29/10/2016
# daylight savings (i.e. +3, -2 etc.)
my $dls = +1;

# Nick - 29/10/2016 - round down the time to minutes, not decimals
my $riseh = $rise;
my $risef = $rise;
$riseh =~ s/\..*//;
$risef =~ s/.*\././;

my $risem = (60*$risef);
$risem =~ s/\..*//;

if ($risem < 10) {
   $risem = "0".$risem;
}

my $seth = $set;
my $setf = $set;
$seth =~ s/\..*//;
$setf =~ s/.*\././;

my $setm = (60*$setf);
$setm =~ s/\..*//;

if ($setm < 10) {
   $setm = "0".$setm;
}

#print "Day: $counter , sunrise: ".($rise-4)." , sunset: ".($set-4)." 
+ ";
print "Day: $counter , sunrise: ".($riseh+$dls).":".$risem." , sunset:
+ ".($seth+$dls).":".$setm;
# End Nick
print "\n";
$counter++;
}

Nick