#!/usr/bin/perl
use strict;
use warnings;
use Math::Trig;
#
# SunAlign.pl (Version for Perl)
#
# Calculates position of sun in sky, as azimuth (compass bearing
# measured clockwise from True North) and angle of elevation, as
# seen from any place on earth, on any date and any time.
# Also calculates alignment of a heliostat mirror.
#
# David Williams
# P.O. Box 48512
# 3605 Lakeshore Blvd. West
# Toronto, Ontario. M8W 4Y6
# Canada
#
# Initially dated 2007 Jul 07 (Basic version)
# This version 2008 August 08
#
# Perl version ported by William K Stewart, 2008 Aug 15
# gmail username = lc([first initial][middle initial][lastname])
#
# All angles in radians except in i/o routines DegIn and DegOut
#
my $PY = 3.1415926536; # "PI" not assignable in some BASICs
my $DR = 180/$PY; # degree / radian factor
my $W = 2 * $PY/365; # earth's mean orbital angular speed in radians/day
my $WR = $PY/12; # earth's speed of rotation relative to sun (radians/hour)
my $C = -23.45/$DR; # reverse angle of earth's axial tilt in radians
my $ST = sin($C); # sine of reverse tilt
my $CT = cos($C); # cosine of reverse tilt
my $E2 = 2 * 0.0167; # twice earth's orbital eccentricity
my $SN = 10 * $W; # 10 days from December solstice to New Year (Jan 1)
my $SP = 12 * $W; # 12 days from December solstice to perihelion
my($LT, $LG, $tAZ, $tEL) = ();
while(1) {
my $S;
do{
print "1. Calculate sun's position\n";
print "2. Calculate mirror orientation\n";
print "3. Calculate both\n";
print "4. Quit program\n\n";
print "Which? (1 - 4): ";
$S = <STDIN>;
chomp $S;
} while(!($S >= 1 && $S <= 4));
print "\nYou selected: $S\n";
last if $S == 4;
print "\nUse negative numbers for directions opposite to those shown.\n\n";
DegIn ("Observer's latitude (degrees North)", \$LT);
DegIn ("Observer's longitude (degrees East)", \$LG);
my $TZN;
do {
print "Time Zone (+- hours from GMT/UT): ";
$TZN = <STDIN>;
chomp $TZN;
} while($TZN !~ /^-?\d+$/);
my $time;
do {
print "Time (HH,MM) (24-hr format): ";
$time = <STDIN>;
chomp $time;
} while($time !~ /^\d+,\d+$/);
my($HR, $MIN) = split(",", $time);
my $date;
do {
print "Date (M#,D#): ";
$date = <STDIN>;
chomp $date;
} while($date !~ /^\d+,\d+$/);
my($Mth, $Day) = split(",", $date);
print "\n";
# co-latitude
my $CL = $PY/2 - $LT;
#day of year (D = 0 on Jan 1)
my $D = (30.6 * (($Mth + 9) % 12) + 58.5 + $Day) % 365;
# orbit angle since solstice at mean speed
my $A = $W * $D + $SN;
# angle with correction for eccentricity
my $B = $A + $E2 * sin($A - $SP);
$C = ($A - atan2(tan($B),$CT))/$PY;
my $tmp = $C + .5;
my $tmp1 = $C - 1;
# solar longitude relative to mean position
#my $SL = $PY * ($C - INT($C + .5));
my $SL = $PY * ($C - int($C + .5));
$C = $ST * cos($B);
# solar declination (latitude)
my $DC = atan2($C, sqrt(1 - $C * $C));
# longitude difference
my $LD = ($HR - $TZN + $MIN / 60) * $WR + $SL + $LG;
# polar axis (perpend'r to azimuth plane)
my($sX, $sY, $sZ);
P2C($LD, $DC, \$sX, \$sY, \$sZ);
# horizontal axis
my($sAZ, $sEL);
C2P($sY, $sZ, $sX, \$sAZ, \$sEL);
# rotate by co-latitude
P2C($sAZ - $CL, $sEL, \$sY, \$sZ, \$sX);
if($sZ < 0) {
print "\n\nSun Below Horizon\n\n";
print "New Calculation\n";
next;
}
if($S != 2) {
# calculate and display sun's position
# vertical axis
C2P($sX, $sY, $sZ, \$sAZ, \$sEL);
DegOut("Sun's azimuth: ", \$sAZ);
DegOut("Sun's elevation: ", \$sEL);
print "\n";
}
if($S > 1) {
# calculate and display mirror orientation
print "For target direction of light reflected from mirror:\n";
DegIn ("Azimuth of target direction (degrees)", \$tAZ);
DegIn ("Elevation of target direction (degrees)", \$tEL);
print "\n";
my($tX, $tY, $tZ);
# target vector X,Y,Z
P2C($tAZ, $tEL, \$tX, \$tY, \$tZ);
my($mAZ, $mEL);
# angle bisection by vector addition
C2P($sX + $tX, $sY + $tY, $sZ + $tZ, \$mAZ, \$mEL);
print "Mirror aim direction (perpendicular to surface):\n";
DegOut("Azimuth: ", \$mAZ);
DegOut("Elevation: ", \$mEL);
print "\n";
}
}
sub Ang {
# calculates angle from positive X axis to vector to (X,Y)
my($X, $Y) = @_;
my $angle;
my $val = sgn($X);
if($val == 1) {
$angle = atan($Y/$X);
}
elsif($val == -1) {
$angle = atan($Y/$X) + $PY;
}
else {
$angle = sgn($Y) * $PY/2;
}
return $angle;
}
sub C2P {
# Cartesian to Polar. Convert from X,Y,Z to AZ,EL
my($X, $Y, $Z, $AZ, $EL) = @_;
$$EL = Ang(sqrt($X * $X + $Y * $Y), $Z);
my $A = Ang($Y, $X);
if($A < $PY) {
$$AZ = $A + $PY;
}
else {
$$AZ = $A - $PY;
}
return;
}
sub DegIn {
# Input angle in degrees and convert to radians
my($P, $X) = @_;
my $N;
do {
print $P,": ";
$N = <STDIN>;
chomp $N;
} while($N !~ /^-?\d+$/);
$$X = $N/$DR;
return;
}
sub DegOut {
# converts radians to degrees, rounds to nearest 0.1, and prints
my($P, $X) = @_;
#my $S = LTRIM(STR(INT(10 * abs($$X * $DR) + .5)));
my $S = 10 * abs($$X * $DR);
$S = round($S);
$S = 0 if($S eq "3600");
$S = "0" . $S if(length($S) == 1);
if($$X < 0) {
if($S) {
$S = "-" . $S;
}
}
$S =~ s/(.+)(.)$/$1\.$2/;
print "$P $S degrees\n";
return;
}
sub P2C {
# Polar to Cartesian. Convert from AZ,EL to X,Y,Z
my($AZ, $EL, $X, $Y, $Z) = @_;
$$Z = sin($EL);
$$X = -cos($EL) * sin($AZ);
$$Y = -cos($EL) * cos($AZ);
return;
}
sub sgn {
my($x) = @_;
if($x < 0) {
return -1
}
elsif($x == 0) {
return 0;
}
else {
return 1;
}
}
sub round {
my($number) = @_;
return int($number + .5 * ($number <=> 0));
}
[ Parent ]
