Correcting for sun diameter and refraction in sunrise/sunset equations

[irotas]

I am using the following C library for a variety of solar position calculations:
http://rredc.nrel.gov/solar/codesandalgorithms/solpos/

In the library, they use the function srss() to compute sunrise and sunset times.

Unfortunately, it seems that their equation is for the time that the center of the sun reaches 0 elevation, and completely ignores the diameter of the sun and atmospheric refraction.

I would like to adjust the srss() function to account for these effects. I found a good explanation online for the required adjustment:
"Sunrise or Sunset is defined for a geocentric altitude of -50 arc minutes; 34 arc minutes is to correct for atmospheric refraction and 16 arc minutes is to correct for the Sun's semidiameter" (although this seems a little suspicious since the refraction depends on atmospheric conditions and the sun is not always the same distance from the earth, but it's probably a decent approximation)

I also found online the angular speed of the earth, which is 7.2921159 × 10-5 radians/second.

My first thought was to simply calculate how long it takes the Earth to rotate 50 arc minutes at the given angular speed, but this isn't correct since the sun doesn't rise/set orthogonally to the horizon.

If there was a way to calculate the angle of the tangent of the trajectory of the sun just as it hits 0 elevation, it may be possible to use a linear approximation to determine how much of an correction to make.

It may also be possible that I'm making this entirely too complicated and there's an easier solution. Or maybe it's much harder than I think. In any case, it's been fun to ponder about and I'm hoping someone here can offer some advice.

Thanks!

-Adam

References:
http://www.adeptscience.co.uk/products/mathsim/mathcad/add-ons/free_ebooks/astro_form_samp.htm
http://hypertextbook.com/facts/2002/JasonAtkins.shtml

 

[irotas]

I found a solution:
http://en.wikipedia.org/wiki/Sunrise_equation#Hour_Angle

The correction can by made by modifying how ‘solpos’ computes the cosine of the sunset hour angle. In function ‘ssha()’, simply modify the computation of ‘cssha’ to account for the sine of -50/60 degrees
Original: cssha = -tdat->sl * tdat->sd / cdcl;
Corrected: cssha = (-0.014543898 - tdat->sl * tdat->sd) / cdcl;

The resulting sunrise/sunset calculations have been compared with Weather.com for a variety of locations at all seasons. Computed times are never more than 1 minute off, and are usually exactly the same.