Calculating insolation

/*  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 <cstdio>
# include <cmath>

/* Constants and global variables */
const double pi=4*atan(1);
const double epsilon=23.4*pi/180;                           //Tilt of the earth's rotational axis
const double solar_contant=1361*86400;                      //Solar constant : J/day/m^2
const int days_in_year=365;                             //Assume that 1 year = 365 synodic days
const int divide_day_into=1440;                         //Number of parts to divide the day into for calculation
const double integration_time=1/(double)divide_day_into;        //Integration time unit
const int N_calc=days_in_year*divide_day_into;              //Total number of times for which stuff has to be calculated
const double VE_time=31+28+20+0.5;                          //Assume that Vernal equinox occurs at noon on March 21

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

/* 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++)
        {
            double t=day+(i+0.5)*integration_time;                      //Time, in days, since new year
            double t_VE=fmod(t+days_in_year-VE_time,days_in_year);          //Time, in days, since vernal equinox
            double lambda=2*pi*t_VE/days_in_year,beta=0,alpha;          //Ecliptic longitude, latitude and Right ascension of the sun
            ecliptic_to_equatorial(lambda,beta,alpha,delta[arr]);           //Find RA and declination of the sun
            h[arr]=2*pi*LST(t_VE)-alpha;                                //Calculate Hour angle of sun : HA = LST -RA
        }
    }
}

int main()
{
    init_HAdec();
    for(int lat=-90;lat<=90;lat++)
    {
        printf("%d\t",lat);
        double phi=lat*pi/180;                                      //Latitude in radians
        for(int day=0,arr=0;day<days_in_year;day++)
        {
            double daily_energy=0;                                  //Total energy received at a place with the given latitude over the day
            for(int i=0;i<divide_day_into;i++,arr++)
            {
                double a,A;                                         //Altitude and azimuth of the sun at the given place and time
                equatorial_to_horizontal(h[arr],delta[arr],phi,A,a);        //Calculate the altitude and azimuth values from HA and dec
                if(a>=0)daily_energy+=integration_time*sin(a)*solar_contant;//If the sun is above the horizon, increment the energy received
            }
            printf("%lf\t",daily_energy);
        }
        printf("\n");
    }
    return 0;
}