1
$\begingroup$

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 David R. Brooks. The code is derived from this link 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. 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;
}
Improve this question
$\endgroup$
8
  • $\begingroup$ 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. $\endgroup$
    – David Hammen
    Commented Jul 27, 2017 at 15:11
  • $\begingroup$ 1.000001018*(1.0-e*e)/(1.0+e*cos(f)) taken from here: [url]books.google.cz/… 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. $\endgroup$
    – user1141649
    Commented Jul 27, 2017 at 15:42
  • $\begingroup$ 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. $\endgroup$
    – David Hammen
    Commented Jul 27, 2017 at 16:05
  • $\begingroup$ @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 $\endgroup$
    – user1141649
    Commented Jul 27, 2017 at 16:35
  • $\begingroup$ 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". $\endgroup$
    – Santiago
    Commented Jul 27, 2017 at 17:44

2 Answers 2

Reset to default
3
$\begingroup$

This doesn't answer your question entirely, but if by "sun-earth radius" you mean the distance between the Earth and the Sun, you can visit https://ssd.jpl.nasa.gov/?horizons with these settings:

enter image description here

to see that the length of the Earth's semimajor axis (the "A" value in each row) doesn't change much over the years (it remains very close to 1 AU, as expected).

However, the semimajor axis is only the "average" distance from the Sun for a given year, which may not be what you want.

If you use these settings:

enter image description here

you can see the actual coordinates of the Earth as viewed from the Sun. From these, you can either compute the distance yourself, or you can look at the "LT" value which is light travel time, and thus gives the distance in light days. Note that I intentionally chose the odd interval of 500 days to demonstrate how the Sun's distance changes during the year.

Additionally, these are "computed" distances, not actually historical records. However, given the inaccuracy of measuring distance to the Sun historically, the computed numbers are probably more accurate.

This question is fine here, but you may want to visit astronomy.stackexchange.com if you have other questions that are more directly based on astronomy.

Improve this answer
$\endgroup$
3
  • $\begingroup$ I tried to look for the coordinates of day 1849-Jul-28 and I got 5,87E-03 light seconds which should be 1,76E+03km. Compared with my program where R=1,01677040986377 AU ... My resultis 1,52107E+11 . 1.76E+03km is 1760 km - too small number... My calculation from coords. d=sqrt(xx + yy + zz) = 6,08E-01*6,08E-01 + -8,15E-01-8,15E-01 + -2,17E-04*-2,17E-04 = 1,02E+00ld. Roughly 155000 km. $\endgroup$
    – user1141649
    Commented Jul 28, 2017 at 10:18
  • $\begingroup$ OK, what I am gonna do is that I export all days from 1600 to 2400 AD. I will calculate the diagonal of the cube box coordinates (d=sqrt(xx*+yy+z*z)) to get distance and diagonal of the cube box velocity vectors to get speed. Then I will save the distances and speed to file and I will add it to my program so it will be able to get the data needed to calculate exact radius. $\endgroup$
    – user1141649
    Commented Jul 28, 2017 at 10:47
  • $\begingroup$ My mistake. At the end of the results, it says "LT = One-way down-leg Newtonian light-time (day)", so it's actually in light days, not light days, not light seconds. I've corrected the post. $\endgroup$
    – user967
    Commented Jul 28, 2017 at 15:18
2
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ The 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. $\endgroup$
    – Fred
    Commented Jan 21, 2018 at 1:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.