Return to Question

Edit: I add the code needed to calculate radiation and solar constant correction.

void getSolRadBird(INPUTS i, SOLPOS solpos, SOLRAD * o)
{
double degRad = 0.017453292519943295769236907684886;
// relative air mass
double AM=1./(cos( solpos.Z*degRad )+0.15*pow(93.885-solpos.Z,-1.25));
double AMp=AM*i.sitePressure/1013.;
// Rayleigh
double Tr=exp(-0.0903*pow(AMp,0.84)*(1+AMp-pow(AMp,1.01)));
// ozone
double Ozm=i.ozone*AM;
double Toz=1.-0.1611*Ozm*pow(1.+139.48*Ozm,-0.3035)-0.002715*Ozm/(1.+0.044*Ozm+0.0003*pow(Ozm,2.));
// mixed gases
double Tm=exp(-0.0127*pow(AMp,0.26));
// water vapor
double Wm=AM * i.water;
// total water vapor
double Tw=1.-2.4959*Wm/((1.+pow(79.034*Wm,0.6828))+6.385*Wm);
// daily turbidity
// Ta5=A*sin((Dan-B)*PI/180.)+C; aerosol optical depth at 500 nm
// Ta3=Ta5+0.1;
double Ta5=i.AOT500;
double Ta3=i.AOT380;
double Tau=0.2758*Ta3+0.35*Ta5;
double Ta=exp((-pow(Tau,0.873))*(1.+Tau-(pow(Tau,0.7088)))*pow(AM,0.9108));

double TAA=1.-0.1*(1.-AM+pow(AM,1.06))*(1.-Ta);
double TAs=Ta/TAA;
double Rs=0.0685+(1.-0.84)*(1.-TAs);
// clear irradiance
double Io=i.SolConst;
// direct
// Earth/sun distance correction, Rsq = 1/R^2
// double d=f.doy;
// double Rsq=(1.00011+0.034221*cos(6.28318*(d-1)/365)+0.00128*sin(6.28318*(d-1)/365)+0.000719*cos(2*(6.28318*(d-1)/365))+0.000077*sin(2*(6.28318*(d-1)/365)));
double Rsq=1.0/ (solpos.R*solpos.R) ;
//alert("R= "+Rsq);
double Id=Rsq*Io*.9662*Tr*Toz*Tm*Tw*Ta;
// direct on horizontal surface
double Idh=Id*cos(solpos.Z*degRad);
// diffuse (scattered)
double Ias=0.79*Io*cos(solpos.Z*degRad)*Toz*Tm*Tw*TAA;
Ias=Ias*(0.5*(1.-Tr)+0.85*(1.-TAs))/(1.-AM+pow(AM,1.02));
// total dif + dir on horizontal
double Itot=(Idh+Ias)/(1.-i.albedo*Rs);
double Idif=Itot-Idh;

// Save the values
o->air_m = AM;
o->direct = Idh;
o->diffuse = Idif;
o->total = Itot;
o->S_corrected = Rsq*Io;
}

I would like to ask you:

Earth/Sun distance correction is made in another function to calculate solar radiation. The code:

// Earth/sun distance correction, Rsq = 1/R^2
// double Rsq=(1.00011+0.034221*cos(6.28318*(d-1)/365)+0.00128*sin(6.28318*(d-1)/365)+0.000719*cos(2*(6.28318*(d-1)/365))+0.000077*sin(2*(6.28318*(d-1)/365)));
    double Rsq=1.0/ (solpos.R*solpos.R) ;

I would like to ask you:

void getSolarPosition(INPUTS * i, SOLPOS *o) {
if ( !i->JulianDate )
   i->JulianDate = getJulianDate(i); // date must be set to calculate Julian Date
double degRad = 0.017453292519943295769236907684886; // PI/180.0
// 2451545.0 - January 1, 2000, at 12:00:00 UT
double T=( i->JulianDate-2451545.0)/36525.0;
    
// L0 - geometric mean longitude of the sun:
double L0=280.46645+36000.76983*T+0.0003032*T*T;
// M - Mean anomaly of the sun:
double M = 357.52910+35999.05030*T-0.0001559*T*T-0.00000048*T*T*T;

double M_rad = M * degRad;
// e - eccentricity of the earth's orbit:
double e=0.016708617-0.000042037*T-0.0000001236*T*T;
// C - sun's center
double C=(1.914600-0.004817*T-0.000014*T*T) * sin(M_rad)
        +(0.019993-0.000101*T)*sin(2.*M_rad)+0.000290*sin(3.*M_rad);
double L_save=(L0+C)/360.;
// L_true - True longitude of the sun
double L_true;
if (L_save < 0.)
   L_true = (L0+C) - ceil(L_save)*360.;
else
   L_true = (L0+C) - floor(L_save)*360.;
if (L_true < 0.) L_true+=360.;
// f - true anomaly of the sun:
double  f = M_rad + C * degRad;
// Earth-sun distance:
o->R =1.000001018*(1.0-e*e)/(1.0+e*cos(f));
// Sidereal time (Theta0)
double Sidereal_time=280.46061837+ 
  360.98564736629*( i->JulianDate-2451545.) + 
  0.000387933*T*T - T*T*T/38710000.;
// Replacement code for Sidereal=fmod(Sidereal,360.)
double S_save=Sidereal_time/360.;
if (S_save < 0.) Sidereal_time=Sidereal_time-ceil(S_save)*360.;
   else Sidereal_time=Sidereal_time-floor(S_save)*360.;
    
if (Sidereal_time < 0.) Sidereal_time+=360.;
    
// Obliquity - Axial tilt - is an astronomical term describing the angle of tilt of the Earth's axis of rotation.
o->obliquity=23.0+26./60.+21.448/3600.-46.8150/3600.*T-0.00059/3600.*T*T + 0.001813/3600.*T*T*T;
    
// right_ascension: tan(alpha)
o->right_ascension = atan2(sin(L_true*degRad)*cos( o->obliquity*degRad ),
                                                   cos(L_true*degRad));
// declination: sin(delta)
o->declination = asin(sin( o->obliquity*degRad )*sin(L_true*degRad));
// hour angle H of the sun with respect to the observer's longitude Lobs
o->hour_angle=Sidereal_time + i->lon - o->right_ascension / degRad;
o->elevation = (asin(sin( i->lat*degRad )*sin( o->declination)+cos(i->lat*degRad)*cos( o->declination )*cos( o->hour_angle*degRad )))/degRad;
// Solar Zenit Angle
o->Z = 90.-o->elevation;
}

).

void getSolarPosition(INPUTS * i, SOLPOS *o) {
if ( !i->JulianDate )
   i->JulianDate = getJulianDate(i); // date must be set to calculate Julian Date
double degRad = 0.017453292519943295769236907684886; // PI/180.0
// 2451545.0 - January 1, 2000, at 12:00:00 UT
double T=( i->JulianDate-2451545.0)/36525.0;
    
// L0 - geometric mean longitude of the sun:
double L0=280.46645+36000.76983*T+0.0003032*T*T;
// M - Mean anomaly of the sun:
double M = 357.52910+35999.05030*T-0.0001559*T*T-0.00000048*T*T*T;

double M_rad = M * degRad;
// e - eccentricity of the earth's orbit:
double e=0.016708617-0.000042037*T-0.0000001236*T*T;
// C - sun's center
double C=(1.914600-0.004817*T-0.000014*T*T) * sin(M_rad)
        +(0.019993-0.000101*T)*sin(2.*M_rad)+0.000290*sin(3.*M_rad);
double L_save=(L0+C)/360.;
// L_true - True longitude of the sun
double L_true;
if (L_save < 0.)
   L_true = (L0+C) - ceil(L_save)*360.;
else
   L_true = (L0+C) - floor(L_save)*360.;
if (L_true < 0.) L_true+=360.;
// f - true anomaly of the sun:
double  f = M_rad + C * degRad;
// Earth-sun distance:
o->R =1.000001018*(1.0-e*e)/(1.0+e*cos(f));
// Sidereal time (Theta0)
double Sidereal_time=280.46061837+ 
  360.98564736629*( i->JulianDate-2451545.) + 
  0.000387933*T*T - T*T*T/38710000.;
// Replacement code for Sidereal=fmod(Sidereal,360.)
double S_save=Sidereal_time/360.;
if (S_save < 0.) Sidereal_time=Sidereal_time-ceil(S_save)*360.;
   else Sidereal_time=Sidereal_time-floor(S_save)*360.;
    
if (Sidereal_time < 0.) Sidereal_time+=360.;
    
// Obliquity - (Axial tilt)
o->obliquity=23.0+26./60.+21.448/3600.-46.8150/3600.*T-0.00059/3600.*T*T + 0.001813/3600.*T*T*T;
    
// right_ascension: tan(alpha)
o->right_ascension = atan2(sin(L_true*degRad)*cos( o->obliquity*degRad ),
                                                   cos(L_true*degRad));
// declination: sin(delta)
o->declination = asin(sin( o->obliquity*degRad )*sin(L_true*degRad));
// hour angle H of the sun with respect to the observer's longitude Lobs
o->hour_angle=Sidereal_time + i->lon - o->right_ascension / degRad;
o->elevation = (asin(sin( i->lat*degRad )*sin( o->declination)+cos(i->lat*degRad)*cos( o->declination )*cos( o->hour_angle*degRad )))/degRad;
// Solar Zenit Angle
o->Z = 90.-o->elevation;
}

). Ceil rounds up; floor rounds down.