Compute total solar irradiance at different latitudes on different days

solarirradiance.cpp

/*  solar.cpp: Calculates the solar energy received by places at different
 *  latitudes over the year. It is assumed that the earth revolves around the
 *  sun in a circular orbit taking exactly 365 days
 *
 *  Author: Raziman T V
 *
 *  License:
 *  You may use this document for whatever purpose you desire, as long as the
 *  following conditions are satisfied:
 *  1) You do not use it in a way that suggests that I endorse the use.
 *  2) If you redistribute this code, I would be grateful if you credit me as
 *    the author of the code. This is not a requirement.
 *  3) If this code is modified and redistributed with credit to me as the
 *    original author, the code should either contain a link to the source where
 *    the unaltered code can be found, or state explicitly that the code has
 *    been altered significantly since.
 */

#include <cmath>
#include <cstdio>

/* Constants and global variables */
const double pi = 4 * atan(1);

// Tilt of the earth's rotational axis
const double epsilon = 23.4 * pi / 180;

// Solar constant : J/day/m^2
const double solar_contant = 1361 * 86400;

// Assume that 1 year = 365 synodic days
const int days_in_year = 365;

// Number of parts to divide the day into for calculation
const int divide_day_into = 1440;

// Integration time unit
const double integration_time = 1 / (double)divide_day_into;

// Total number of times for which stuff has to be calculated
const int N_calc = days_in_year * divide_day_into;

// Assume that Vernal equinox occurs at noon on March 21
const double VE_time = 31 + 28 + 20 + 0.5;

// Hour angle and declination of the sun for each time
double h[N_calc], delta[N_calc];

/* Convert from ecliptic to equatorial coordinates */
void ecliptic_to_equatorial(double lambda, double beta, double &alpha,
                            double &delta) {
  /*               sin[delta] = sin[beta] * cos[epsilon]
   *                          + cos[beta] * sin[epsilon] * sin[lambda]
   *  cos[delta] * sin[alpha] = cos[beta] * sin[lambda] * cos[epsilon]
   *                          - sin[beta] * sin[epsilon]
   *  cos[delta] * cos[alpha] = cos[beta] * cos[lambda]
   */
  delta =
      asin(sin(beta) * cos(epsilon) + cos(beta) * sin(epsilon) * sin(lambda));
  alpha =
      atan2(cos(beta) * sin(lambda) * cos(epsilon) - sin(beta) * sin(epsilon),
            cos(beta) * cos(lambda));
}

/* Convert from equatorial to horizontal coordinates */
void equatorial_to_horizontal(double h, double delta, double phi, double &A,
                              double &a) {
  /*           sin[a] = sin[phi] * sin[delta]
   *                  + cos[phi] * cos[delta] * cos[h]
   *  cos[a] * cos[A] = cos[delta] * cos[h] * sin [phi] - sin[delta] * cos[phi]
   *  cos[a] * sin[A] = cos[delta] * sin[h]
   */
  a = asin(sin(phi) * sin(delta) + cos(phi) * cos(delta) * cos(h));
  A = atan2(cos(delta) * sin(h),
            cos(delta) * cos(h) * sin(phi) - sin(delta) * cos(phi));
}

/* Convert from time (in days) since vernal equinox to Local Sidereal Time */
double LST(double t_VE) {
  /*  days_in_year synodic days = (days_in_year+1) sidereal days
   *  Hence sidereal time progresses at the rate of 1+(1/days_in_year) times
   *  normal time
   *  Sidereal time at vernal equinox is 12h
   */
  return fmod(t_VE * (1 + 1 / (double)days_in_year) + 0.5, 1);
}

/* Calculate the Hour angle and Declination of the sun for each time value */
void init_HAdec() {
  for (int day = 0, arr = 0; day < days_in_year; day++) {
    for (int i = 0; i < divide_day_into; i++, arr++) {
      // Time, in days, since new year
      double t = day + (i + 0.5) * integration_time;

      // Time, in days, since vernal equinox
      double t_VE = fmod(t + days_in_year - VE_time, days_in_year);

      // Ecliptic longitude, latitude and Right ascension of the sun
      double lambda = 2 * pi * t_VE / days_in_year, beta = 0, alpha;

      // Find RA and declination of the sun
      ecliptic_to_equatorial(lambda, beta, alpha, delta[arr]);

      // Calculate Hour angle of sun : HA = LST -RA
      h[arr] = 2 * pi * LST(t_VE) - alpha;
    }
  }
}

int main() {
  init_HAdec();
  for (int lat = -90; lat <= 90; lat++) {
    printf("%d\t", lat);

    // Latitude in radians
    double phi = lat * pi / 180;

    for (int day = 0, arr = 0; day < days_in_year; day++) {
      // Total energy received at a place with the given latitude over the day
      double daily_energy = 0;

      for (int i = 0; i < divide_day_into; i++, arr++) {
        // Altitude and azimuth of the sun at the given place and time
        double a, A;

        // Calculate the altitude and azimuth values from HA and dec
        equatorial_to_horizontal(h[arr], delta[arr], phi, A, a);

        // If the sun is above the horizon, increment the energy received
        if (a >= 0)
          daily_energy += integration_time * sin(a) * solar_contant;
      }
      printf("%lf\t", daily_energy);
    }
    printf("\n");
  }
  return 0;
}