//writen by: E. van der Goot and Giorgio Cifani
//JRC, I-21020 Ispra, Italy
//date: 1996

//based on the original FORTRAN/SQL CGMS and WOFOST code
//produced by Tamme van der Wal, Kees van Diepen and Daniel van Kralingen
//Alterra, P.O. Box 6700 AA, Wageningen, The Netherlands

SAstroResultsAstro(int nDay, double dLatitude)
{
SAstroResults ar;

// This procedure does not work very well for higher latitudes, especially in winter.
// a possible fix checks for the validity of dAob (abs(dAob) <= 1.0)
double dAob;
double dAobCorrected;

double dCosLD;
double dSinLD;
float fDayL;
float fDayLP;
float fDSinB;
float fDSinBE;
float fAngot;

double dFrac;
double dSin2;
double dCos2;

// The following formula does not work very well for very high latitudes (> 66?) and
// nDay around december/january
// The procedure checks for extreme values and tries to correct as much as possible

dCosLD = CosSolDeclAngle[nDay] * cos(toRAD * dLatitude);
dSinLD = SinSolDeclAngle[nDay] * sin(toRAD * dLatitude);

dAob = dSinLD / dCosLD;

if (fabs((float)dAob) < 1.0)
        {
        dFrac = dCosLD * sqrt(1.0 - dAob * dAob) / PI;
        fDayL = (float)(12.0 * (1.0 + 2.0 * asin(dAob) / PI));
        fDSinB = (float)(3600 * (dSinLD * (fDayL) + 24.0 * dFrac));
        dSin2 = dSinLD * dSinLD;
        dCos2 = dCosLD * dCosLD;
        fDSinBE = (float)(3600.0 * (fDayL * (dSinLD + 0.4 * (dSin2 + (dCos2 / 2.0))) +
                        (12.0 * dFrac * (2.0 + 3.0 * 0.4 * dSinLD))));

        fAngot = (float)(SolarConstant[nDay] * fDSinB / 1000.0f); //(KJ/m2day)
        }
else if ((float)dAob > 1.0)
        {
        fDayL = 24.0f;
        fDSinB = (float)(3600 * (dSinLD * 24.0));
        dSin2 = dSinLD * dSinLD;
        dCos2 = dCosLD * dCosLD;
        fDSinBE = (float)(3600.0 * (24 * (dSinLD + 0.4 * (dSin2 + (dCos2 / 2.0)))));
        fAngot = (float)(SolarConstant[nDay] * fDSinB / 1000.0f); // (KJ/m2day)
        }
else if ((float)dAob < -1.0)
        {
        fDayL = 0.0f;
        fDSinBE = 0.0f;
        fAngot = 0.0f;
        }

dAobCorrected = (-sin(-4.0 * toRAD) + dSinLD) / dCosLD;

if (fabs((float)dAobCorrected) < 1.0)
        fDayLP = (float)(12.0 * (1.0 + 2.0 * asin(dAobCorrected) / PI));

        fDayLP = 0.0f;

ar.dCosLD         = dCosLD;
ar.dSinLD          = dSinLD;
ar.fDayL            = fDayL;
ar.fDayLP          = fDayLP;
ar.fDSinBE         = fDSinBE;
ar.fAngot            = fAngot;

return ar;
}