sun_earth_distance_correction returns the square of the sun-earth distance relativ to 1 AU

[depion]

Problem description

The function sun_earth_distance_correction in pyorbital/astronomy.py returns the square of the sun-earth distance relative to 1 AU. However, the function apply_earthsun_distance_correction in satpy/readers/utils.py squares the returned value again which leads to inaccurate calculations of the reflectance for the visible channels (e.g., vis_calibrate in satpy/readers/seviri_base.py).

The line

corr = 1 - 0.0334 * np.cos(2 * np.pi * (jdays2000(utc_time) - 2) / year)

in sun_earth_distance_correction calculating the squared relative distance should be replaced by

corr = 1 - 0.0167 * np.cos(2 * np.pi * (jdays2000(utc_time) - 2) / year)

as follows from the formula given in the comments of sun_earth_distance_correction, i.e.,

r = a(1-e e)/(1+e * np.cos(theta))

thus

corr = r / AU = a(1-e e)/AU 1/(1+e np.cos(theta)) =approx a(1-e e)/AU (1- e np.cos(theta))

Edit, 11 July 2022: The day of year with Earth in perihelion varies between Jan 2 and Jan 5, see http://www.astropixels.com/ephemeris/perap2001.html or https://web.archive.org/web/20080328044924/http://aa.usno.navy.mil/data/docs/EarthSeasons. Jan 2 was used up to now, but due to the variability Jan 3 would be a better choice.

[mraspaud]

@depion Thanks a lot for the heads up! It's a bit hard to read in text form, so here it is in math:

$$r = {a (1-e^2) \over (1+e \cos(\theta))}$$

thus

$$corr = {r \over AU} = {{a(1-e^2) \over AU} {1 \over (1+e \cos(\theta))}} \approx {{a(1-e^2) \over AU} (1- e \cos(\theta))}$$

Is this correct?

And do you feel like making a pull request for it?

[depion]

Hi @mraspaud! That looks correct. Yes, I will try making a pull request and correcting this.