gschwind / sg2

GNU Lesser General Public License v3.0
9 stars 2 forks source link

Fatal error in topocentric_correction_refraction_ZIM #1

Closed jrleal closed 2 years ago

jrleal commented 2 years ago

Hello, There is an error in the topocentric_correction_refraction_ZIM function that prevents it from working correctly. In that function, tan_gamma_S0 is initialized to 0.0, and then used without assigning the correct value (tan_gamma_S0 = tan(gamma_S0);). There is also another error in the calculation of gamma_S0_4. It is calculated as gamma_S0_4 = gamma_S0_4 + gamma_S0 when the correct is gamma_S0_4 = gamma_S0_3 + gamma_S0

inline double topocentric_correction_refraction_ZIM(double const gamma_S0, double const P, double const T)
{
    static const double gamma_S0_seuil = -0.010035643198967;
    static const double R = 0.029614018235657;
    /*(tan(gamma_S0_seuil + 0.0031376 / (gamma_S0_seuil+ 0.089186))) */
    double K;
    double tan_gamma_S0 = 0.0; // <= Initialized to 0.0
    double gamma_S0_2 = 0.0, gamma_S0_3 = 0.0, gamma_S0_4 = 0.0;
    unsigned long k;

    K = (P / 1013.0) * (283. / (273. + T)) * 4.848136811095360e-006;

    if (gamma_S0 <= -0.010036) {
        return gamma_S0 + (-20.774 / tan_gamma_S0) * K; // Error: Use tan_gamma_S0 without calculating its value!
    } else if (gamma_S0 <= 0.087266) {
        gamma_S0_2 = gamma_S0 * gamma_S0;
        gamma_S0_3 = gamma_S0_2 * gamma_S0;
        gamma_S0_4 = gamma_S0_4 * gamma_S0; // Error: must be gamma_S0_3 * gamma_S0
        return gamma_S0
                + (1735 - 2.969067e4 * gamma_S0 + 3.394422e5 * gamma_S0_2
                + -2.405683e6 * gamma_S0_3 + 7.66231727e6 * gamma_S0_4) * K;
    } else if (gamma_S0 <= 1.483529864195180) {
        return gamma_S0 + K * (58.1 / tan_gamma_S0 - 0.07 / pow(tan_gamma_S0, 3.0) + 0.000086 / pow(tan_gamma_S0, 5.0)); // Error: Use tan_gamma_S0 without calculating its value!
    } else {
        return NAN;
    }

}
gschwind commented 2 years ago

Thanks you for reporting the issue, I did a fix in the devel branch that I will push in the next release:

5cb9d273a9706e347be1ca75aff94aad0d6cdd63

Best regards.

gschwind commented 2 years ago

Hello,

It should be fixed in the version 2.3.0

Best regards