Dmitri

#include <stdio.h>
#include <math.h>

const double PI = 3.1415926535897932384626433832795028841971;
const double DR = 180 / PI; //degree / radian factor;
const double W = 2 * PI / 365; // earth's mean orbital angular speed in radians/day
const double WR = PI / 12; // earth's speed of rotation relative to sun (radians/hour)
double C = -23.45 / DR; // reverse angle of earth's axial tilt in radians
const double ST = sin(C); // sine of reverse tilt
const double CT = cos(C); // cosine of reverse tilt
const double E2 = 2 * 0.0167; // twice earth's orbital eccentricity
const double SN = 10 * W; // 10 days from December solstice to New Year (Jan 1)
const double SP = 12 * W; // 12 days from December solstice to perihelion

double rad2deg(double radians)
{
    return (180.0/PI)*radians;
}

double deg2rad(double degrees)
{
    return (PI/180.0)*degrees;
}

double ang(double x, double y)
{
        double ang;
        if (x > 0)
          ang = atan(y/x);
        else if (x < 0)
          ang = atan(y/x) + PI;
        else if (y > 0)
          ang = 1 * PI / 2;
        else if (y < 0)
          ang = -1 * PI / 2;
        else
          ang = 0;

        return ang;
}

void P2C(double AZ, double EL, double *x, double *y, double *z)
{
        double c;
        *z = sin(EL);
        c = cos(EL);
        *x = c * sin(AZ);
        *y = c * cos(AZ);
}

void C2P(double x, double y, double z, double *AZ, double *EL)
{
        *EL = ang(sqrt(x * x + y * y), z);
        *AZ = ang(y, x);
        if (*AZ < 0)
          *AZ = *AZ + 2 * PI;
}

int main (int argc, char *argv[])
{
        double A, B, SL, D, CL, DC, LD, sAZ, sEL, mAZ, mEL, tAZ, tEL, tX, tY, tZ, sX, sY, sZ, LT, LG;
        int TZN, Mth, Day, Hr, Mn;

        //User input. Note: For brevity, no error checks on user inputs
        printf("Use negative numbers for directions opposite to those shown.\n\n");
        printf("Observer's latitude (degrees North): ");
        scanf("%lf", &LT);
        LT = deg2rad(LT);

        printf("Observer's longitude (degrees East): ");
        scanf("%lf", &LG);
        LG = deg2rad(LG);

        printf("\nFor target direction of light reflected from mirror:\n");
        printf("Azimuth of target direction (degrees): ");
        scanf("%lf", &tAZ);
        tAZ = deg2rad(tAZ);

        printf("Elevation of target direction (degrees): ");
        scanf("%lf", &tEL);
        tEL = deg2rad(tEL);

        printf("\nTime Zone (+/- hours from GMT/UT): ");
        scanf("%d", &TZN);

        printf("Month: ");
        scanf("%d", &Mth);

        printf("Day: ");
        scanf("%d", &Day);

        printf("Hour (24-hr format): ");
        scanf("%d", &Hr);

        printf("Minutes: ");
        scanf("%d", &Mn);


        CL = PI / 2 - LT; // co-latitude
        D = int(30.6 * ((Mth + 9) % 12) + 58.5 + Day) % 365; // day of year (D = 0 on Jan 1)
        A = W * D + SN; // orbit angle since solstice at mean speed
        B = A + E2 * sin(A - SP); // angle with correction for eccentricity
        C = (A - atan(tan(B) / CT)) / PI;
        SL = PI * (C - floor(C + .5)); // solar longitude relative to mean position
        C = ST * cos(B);
        DC = atan(C / sqrt(1 - C * C)); // solar declination (latitude)
        LD = (Hr - TZN + Mn / 60.0) * WR + SL + LG; // longitude difference

        P2C(LD, DC, &sX, &sY, &sZ); // polar axis (perpend'r to azimuth plane)
        C2P(sY, sZ, sX, &sAZ, &sEL); // horizontal axis
        P2C(sAZ + CL, sEL, &sY, &sZ, &sX); // rotate by co-latitude

        // Sun's position
        C2P(sX, sY, sZ, &sAZ, &sEL); // vertical axis

        printf("\nSun Azimuth:          %.1f degrees\n",rad2deg(sAZ));
        printf("Sun Elevation:        %.1f degrees\n\n",rad2deg(sEL));

        // Mirror aim direction
        P2C(tAZ, tEL, &tX, &tY, &tZ);
        C2P(sX + tX, sY + tY, sZ + tZ, &mAZ, &mEL);
        printf("Mirror aim direction (perpendicular to surface): \n");
        printf("Mirror Azimuth:       %.1f degrees\n",rad2deg(mAZ));
        printf("Mirror Elevation:     %.1f degrees\n\n",rad2deg(mEL));

        return 0;
}