//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 NetherlandsSAstroResultsAstro(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 possibledCosLD = 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;
}