<p>Something that I noticed in your formula is that you are using 365 for the amount of days in a year.  Perhaps trying it using 365.25 would be more correct as it accounts for the leap year.  Please let me know if this helps/corrects the issue.</p>

CC BY-SA 3.0
2018-01-21T01:46:14.213Z
The [Gregorian calendar](https://en.wikipedia.org/wiki/Gregorian_calendar) is more accurate than the Julian. Unlike the Julian calendar, which has 365.25 days each years, the Georgian calendar has 365.2425 (365 days 5 hours 49 minutes & 12 seconds) each year.
Fred
CC BY-SA 3.0
2018-01-20T20:58:55.933Z
2018-01-20T20:58:55.933Z
Peter Denton
<p><strong>Purpose of my question:</strong>
I create program to calculate solar radiation and I need to calculate radius between sun and earth.</p>

<p>The based on book "Guide to HTML, JavaScript and PHP"
For Scientists and Engineers, By David R. Brooks. The code is derived from this <a href="http://www.instesre.org/Solar/BirdModelNew.htm" rel="nofollow noreferrer">link</a> which is a calculator. I edited the code to C.</p>

<p><strong>Radius calculation</strong></p>

<p>In my code the radius calculation is defined as
<code>o-&gt;R =1.000001018*(1.0-e*e)/(1.0+e*cos(f));</code>
where <em>e</em> - eccentricity of the earth's orbit:
      <em>f</em>- true anomaly of the sun:</p>

<pre><code>void getSolarPosition(INPUTS * i, SOLPOS *o) {
if ( !i-&gt;JulianDate )
   i-&gt;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-&gt;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 &lt; 0.)
   L_true = (L0+C) - ceil(L_save)*360.;
else
   L_true = (L0+C) - floor(L_save)*360.;
if (L_true &lt; 0.) L_true+=360.;
// f - true anomaly of the sun:
double  f = M_rad + C * degRad;
// Earth-sun distance:
o-&gt;R =1.000001018*(1.0-e*e)/(1.0+e*cos(f));
// Sidereal time (Theta0)
double Sidereal_time=280.46061837+ 
  360.98564736629*( i-&gt;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 &lt; 0.) Sidereal_time=Sidereal_time-ceil(S_save)*360.;
   else Sidereal_time=Sidereal_time-floor(S_save)*360.;

if (Sidereal_time &lt; 0.) Sidereal_time+=360.;

// Obliquity - (Axial tilt)
o-&gt;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-&gt;right_ascension = atan2(sin(L_true*degRad)*cos( o-&gt;obliquity*degRad ),
                                                   cos(L_true*degRad));
// declination: sin(delta)
o-&gt;declination = asin(sin( o-&gt;obliquity*degRad )*sin(L_true*degRad));
// hour angle H of the sun with respect to the observer's longitude Lobs
o-&gt;hour_angle=Sidereal_time + i-&gt;lon - o-&gt;right_ascension / degRad;
o-&gt;elevation = (asin(sin( i-&gt;lat*degRad )*sin( o-&gt;declination)+cos(i-&gt;lat*degRad)*cos( o-&gt;declination )*cos( o-&gt;hour_angle*degRad )))/degRad;
// Solar Zenit Angle
o-&gt;Z = 90.-o-&gt;elevation;
}
</code></pre>

<p>This is function to calculate solar position. It works exactly the same as Bird and Hulstrom's Solar Irradiance Model refered the calculator (see link above). Here I use <em>i</em> input object where input data are are saved and <em>o</em> ouput object where the calculated data regarding solar position are saved after they are calculated. atan2 - Returns the principal value of the arc tangent of y/x, expressed in radians (whatever it means, this is taken from C/C++ manual - </p>

<ul>
<li>I am not mathematician</li>
</ul>

<p>). Ceil rounds up; floor rounds down.</p>

<p>The problem is that if I set old date like
<strong>1849</strong>/06/31 11:15 The Solar constant corrected to Radius does not fit the <a href="http://lasp.colorado.edu/lisird/tsi/historical_tsi.html" rel="nofollow noreferrer">historical records</a>. In the case the result would 1322.3 be for SolConst 1367. Which is crazy. According historical data it should be 1361.035. So I expect the radius is wrong calculated.</p>

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

<pre><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) ;
</code></pre>

<p>I would like to ask you:</p>

<p>1) where can I get historical records of the sun-earth radius</p>

<p>2) what is wrong with this formula or this calculation? Can you suggest better formula?</p>

<p>Edit:
I add the code needed to calculate radiation and solar constant correction.</p>

<pre><code>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-&gt;air_m = AM;
o-&gt;direct = Idh;
o-&gt;diffuse = Idif;
o-&gt;total = Itot;
o-&gt;S_corrected = Rsq*Io;
}
</code></pre>

CC BY-SA 3.0
2017-07-27T15:11:51.010Z
Your code only shows how you are computing `o->R` and other parameters related to the Earth's orbit. You are not showing how you are correcting the solar constant.
David Hammen
CC BY-SA 3.0
2017-07-27T15:42:28.763Z
`1.000001018*(1.0-e*e)/(1.0+e*cos(f))` taken from here: [url]https://books.google.cz/books?id=NTh1O5X47FEC&pg=PA190&lpg=PA190&dq=1.000001018&source=bl&ots=yFm62MILCT&sig=qUu90Ws4qS7MKHQUZkvP-69OCBk&hl=en&sa=X&ved=0ahUKEwjr167p46nVAhXQh7QKHQOQDhwQ6AEIKzAB#v=onepage&q=1.000001018&f=false[/url] You can also to open the calculator and view the source code in JavaScript. I think the calculation of the Solar Constant correction or radiation is not needed for calculating R.
user1141649
CC BY-SA 3.0
2017-07-27T16:05:05.257Z
That is not the question I asked. You are showing how you calculate `0->R` rather than the solar constant. Your calculation of the radius looks correct. Your bug most likely is somewhere else. That said, since your code is in C, learn to use your debugger.
David Hammen
CC BY-SA 3.0
2017-07-27T16:35:02.247Z
@David Hammen: Q updated. See o->S_corrected = Rsq*Io;. The program in C works exactly same as the HTML/Javascript calculator. The output value for SolConst corrected in the calculator is wrong. Wait a moment I will add the function with the solar constant and radiation calculation
user1141649
CC BY-SA 3.0
2017-07-27T17:44:36.170Z
First: You should verify your formula using more than a single point. Second: Is irradiance affected by sunspots? At least your reference includes a link to "The Impact of the Revised Sunspot Record on Solar Irradiance Reconstructions".
Santiago
CC BY-SA 3.0
2017-07-27T18:54:18.983Z
@Barry Carter: The main problem is not that the historical curve is non-smooth, the problem is that the calculations result in Solar Constant corrected which is non-sense because it goes down under the 1360 W*m2 (I got 1322.28). This is wrong. I thought it is cause by "Radius" R, but as David Hammen notes, the problem probably lays in the function getSolRadBird - in the Bird model calculation. The good news is that the Radisus R should be fine. Also, thank you for your answer I will try it.
user1141649
CC BY-SA 3.0
2018-01-22T20:25:43.383Z
I think there is a more fundamental problem: before the advent of space-based measurements (ca 1970), it simply wasn't possible to measure the solar constant with any accuracy.
jamesqf
<earth-history><astronomy><orbit>
Earth's orbit and historical sun-earth distance
user1141649
# Earth's orbit and historical sun-earth distance

(Asked by: username_0 on 2017-07-27)

**Purpose of my question:**
I create program to calculate solar radiation and I need to calculate radius between sun and earth.

The based on book "Guide to HTML, JavaScript and PHP"
For Scientists and Engineers, By username_1 R. Brooks. The code is derived from this [link](http://www.instesre.org/Solar/BirdModelNew.htm) which is a calculator. I edited the code to C.

**Radius calculation**

In my code the radius calculation is defined as
`o->R =1.000001018*(1.0-e*e)/(1.0+e*cos(f));`
where *e* - eccentricity of the earth's orbit:
*f*- true anomaly of the sun:

```
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;
}

```

This is function to calculate solar position. It works exactly the same as Bird and Hulstrom's Solar Irradiance Model refered the calculator (see link above). Here I use *i* input object where input data are are saved and *o* ouput object where the calculated data regarding solar position are saved after they are calculated. atan2 - Returns the principal value of the arc tangent of y/x, expressed in radians (whatever it means, this is taken from C/C++ manual -

* I am not mathematician

). Ceil rounds up; floor rounds down.

The problem is that if I set old date like
**1849**/06/31 11:15 The Solar constant corrected to Radius does not fit the [historical records](http://lasp.colorado.edu/lisird/tsi/historical_tsi.html). In the case the result would 1322.3 be for SolConst 1367. Which is crazy. According historical data it should be 1361.035. So I expect the radius is wrong calculated.

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:

1) where can I get historical records of the sun-earth radius

2) what is wrong with this formula or this calculation? Can you suggest better formula?

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;
}

```

**Tags:** earth-history, astronomy, orbit

> **Comments:**
> **username_1**: Your code only shows how you are computing `o->R` and other parameters related to the Earth's orbit. You are not showing how you are correcting the solar constant.  
> **username_0**: `1.000001018*(1.0-e*e)/(1.0+e*cos(f))` taken from here: [url]https://books.google.cz/books?id=NTh1O5X47FEC&pg=PA190&lpg=PA190&dq=1.000001018&source=bl&ots=yFm62MILCT&sig=qUu90Ws4qS7MKHQUZkvP-69OCBk&hl=en&sa=X&ved=0ahUKEwjr167p46nVAhXQh7QKHQOQDhwQ6AEIKzAB#v=onepage&q=1.000001018&f=false[/url] You can also to open the calculator and view the source code in JavaScript. I think the calculation of the Solar Constant correction or radiation is not needed for calculating R.  
> **username_1**: That is not the question I asked. You are showing how you calculate `0->R` rather than the solar constant. Your calculation of the radius looks correct. Your bug most likely is somewhere else. That said, since your code is in C, learn to use your debugger.  
> **username_0**: @username_1: Q updated. See o->S_corrected = Rsq*Io;. The program in C works exactly same as the HTML/Javascript calculator. The output value for SolConst corrected in the calculator is wrong. Wait a moment I will add the function with the solar constant and radiation calculation  
> **username_2**: First: You should verify your formula using more than a single point. Second: Is irradiance affected by sunspots? At least your reference includes a link to "The Impact of the Revised Sunspot Record on Solar Irradiance Reconstructions".  
> **username_0**: @Barry Carter: The main problem is not that the historical curve is non-smooth, the problem is that the calculations result in Solar Constant corrected which is non-sense because it goes down under the 1360 W*m2 (I got 1322.28). This is wrong. I thought it is cause by "Radius" R, but as username_1 username_1 notes, the problem probably lays in the function getSolRadBird - in the Bird model calculation. The good news is that the Radisus R should be fine. Also, thank you for your answer I will try it.  
> **username_3**: I think there is a more fundamental problem: before the advent of space-based measurements (ca 1970), it simply wasn't possible to measure the solar constant with any accuracy.  

---

## 1 answers

### Answer by username_4 on 2018-01-20. Score: 2

Something that I noticed in your formula is that you are using 365 for the amount of days in a year. Perhaps trying it using 365.25 would be more correct as it accounts for the leap year. Please let me know if this helps/corrects the issue.

> **Comments:**
> **username_5**: The [Gregorian calendar](https://en.wikipedia.org/wiki/Gregorian_calendar) is more accurate than the Julian. Unlike the Julian calendar, which has 365.25 days each years, the Georgian calendar has 365.2425 (365 days 5 hours 49 minutes & 12 seconds) each year.  

---

