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
In reply to Re: calculate length of day as function of space at onset of fall
by Linicks
in thread calculate length of day as function of space at onset of fall
by Aldebaran
| For: | Use: | ||
| & | & | ||
| < | < | ||
| > | > | ||
| [ | [ | ||
| ] | ] |