Sun's position, equation of time: WRONG BY 27 DAYS

Question

I'm trying to improve the U.S. Naval Observatory's algorithm for calculating the position of the sun, by writing it in a simpler way where all the numbers used have been properly identified - so normal people can understand it as well.

But somehow, the equation of time is off by 27 days. Maybe someone here can spot what's wrong?

Test run:

1) Failure:
test_equation_of_time(TestSolarCalculations):
Expected 2011-10-01 10:23:00 UTC, not 2011-10-28 10:59:31 UTC.

2) Failure:
test_suns_declination(TestSolarCalculations):
Expected -3.18, not -3.2087920449753007.

New algorithm:

# Constants for J2000.0 / 1 Jan 2000 12:00Z:
#
EPOCHS_JULIAN_DATE = 2451545
ANOMALISTIC_YEAR_IN_DAYS = 365.259636
TROPICAL_YEAR_IN_DAYS = 365.2421897
SUNS_MEAN_ANOMALY_AT_EPOCH = degrees_to_radians 357.5291
SUNS_MEAN_LONGITUDE_AT_EPOCH = degrees_to_radians 280.459
SUNS_GEODETIC_PRECESSION = degrees_to_radians 1.915
EARTHS_ORBITAL_ECCENTRICITY = 0.020
EARTHS_ADJUSTED_AVERAGE_RADIUS_IN_AU = 1.00014
EARTHS_APPROXIMATE_ATMOSPHERIC_REFRACTION = degrees_to_radians 0.01671
EARTHS_ECLIPTIC_MEAN_OBLIQUITY = degrees_to_radians 23.439
EARTHS_ECLIPTIC_OBLIQUITY_CHANGE_RATE = degrees_to_radians 0.00000036

def days_from_epoch
  @todays_julian_date - EPOCHS_JULIAN_DATE
end

def suns_daily_mean_anomaly_change_rate
  degrees_to_radians(360 / ANOMALISTIC_YEAR_IN_DAYS)
end

def suns_mean_anomaly
  SUNS_MEAN_ANOMALY_AT_EPOCH + suns_daily_mean_anomaly_change_rate * days_from_epoch
end

def suns_daily_mean_longitude_change_rate
  degrees_to_radians(360 / TROPICAL_YEAR_IN_DAYS)
end

def suns_mean_longitude
  SUNS_MEAN_LONGITUDE_AT_EPOCH + suns_daily_mean_longitude_change_rate * days_from_epoch
end

def suns_apparent_ecliptic_longitude
  suns_mean_longitude + SUNS_GEODETIC_PRECESSION * sin(suns_mean_anomaly) + EARTHS_ORBITAL_ECCENTRICITY * sin(2 * suns_mean_anomaly)
end

def suns_distance_from_earth_in_au
  EARTHS_ADJUSTED_AVERAGE_RADIUS_IN_AU - EARTHS_APPROXIMATE_ATMOSPHERIC_REFRACTION * cos(suns_mean_anomaly) - (EARTHS_APPROXIMATE_ATMOSPHERIC_REFRACTION ^ 2 / 2) * cos(2 * suns_mean_anomaly)
end

def earths_ecliptic_mean_obliquity
  EARTHS_ECLIPTIC_MEAN_OBLIQUITY - EARTHS_ECLIPTIC_OBLIQUITY_CHANGE_RATE * days_from_epoch
end

def suns_right_ascension
  atan2(cos(suns_apparent_ecliptic_longitude) * sin(suns_apparent_ecliptic_longitude), cos(suns_apparent_ecliptic_longitude)) / 15
end

# Time.utc(Time.now.year) => 2011-01-01 00:00:00 UTC
#
def equation_of_time
  Time.utc(Time.now.year) + (radians_to_degrees(suns_mean_longitude) / 15 - suns_right_ascension) * 60 * 60 * 24
end

def suns_declination
  radians_to_degrees(asin(sin(earths_ecliptic_mean_obliquity) * sin(suns_apparent_ecliptic_longitude)))
end

Thanks!

Mats

Answer 1

May be it is an issue with the precision of Ruby floating point numbers, not an issue with your program's logic? I haven't used ruby, but it seems there are some issues in ruby floating point operations exists according to thisposts. Issue with precision of Ruby math operations

According to the post, you may need to use BigDecimal class for your floating point operations. I actually didn't go through your program, this is just a guess.

Answer 2

Some of my code.

# From angles.rb:<br>
# eccentricity of elliptical Earth orbit around Sun
# Horner calculation method  
def eccentricity_Earth( ta = A2000 )
  ta = check_jct_zero( ta )      
  # 0.016708617 - ta[ 0 ] * ( 0.000042037 + ta[ 0 ] * 0.0000001235 )
  [-0.0000001235, -0.000042037, 0.016708617].inject(0.0) {|p, a| p * ta[0] + a}
end
alias_method :eccentricity_earth_orbit, :eccentricity_Earth

From: Equation of Time Ruby gem

Where it's used?

# From angles.rb:<br>
# equation of centre
# added to mean anomaly to get true anomaly. 
def center( ta = A2000)
  ta = check_jct_zero( ta )    
  sine_1M = sin( 1.0 * deg_to_rad( @ma ) )
  sine_2M = sin( 2.0 * deg_to_rad( @ma ) )
  sine_3M = sin( 3.0 * deg_to_rad( @ma ) )
  sine_4M = sin( 4.0 * deg_to_rad( @ma ) )
  sine_5M = sin( 5.0 * deg_to_rad( @ma ) )
  e = eccentricity_Earth( ta )
  rad_to_deg( sine_1M * (     2.0  * e    - e**3/4.0 + 5/96.0 * e**5 ) +  
              sine_2M * (   5/4.0  * e**2 - 11/24.0 * e**4 )           + 
              sine_3M * ( 13/12.0  * e**3 - 43/64.0 * e**5 )           +
              sine_4M *  103/96.0  * e**4                              +
              sine_5M * 1097/960.0 * e**5                              )

  # sine_1M *( 1.914602 - ta[ 0 ] * ( 0.004817 + ta[ 0 ] * 0.000014 )) +
  # sine_2M *( 0.019993 - ta[ 0 ] * 0.000101 )                         +
  # sine_3M *  0.000289
end
alias_method :equation_of_center, :center

Answer 3

Here:

def suns_apparent_ecliptic_longitude
  suns_mean_longitude 
  + SUNS_GEODETIC_PRECESSION * sin(suns_mean_anomaly) 
  + EARTHS_ORBITAL_ECCENTRICITY * sin(2 * suns_mean_anomaly)
end

you are mixing units. suns_mean_longitude is radians, as is SUNS_GEODETIC_PRECESSION (and therefore so is the second term), but EARTHS_ORBITAL_ECCENTRICITY is 0.020 degrees.

Additionally:

EARTHS_APPROXIMATE_ATMOSPHERIC_REFRACTION = degrees_to_radians 0.01671

should just be

EARTHS_APPROXIMATE_ATMOSPHERIC_REFRACTION = 0.01671

I think, since this is used in a formula that determines a distance.

Additionally additionally,

def suns_right_ascension
  atan2(cos(suns_apparent_ecliptic_longitude) * sin(suns_apparent_ecliptic_longitude), cos(suns_apparent_ecliptic_longitude)) / 15
end

I haven't checked what units this atan2 function returns - are you sure it's degrees, which it would have to be in order for dividing by 15 to make sense?

There may be other problems; these are only those that I can see.