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;
    }
}