Below I report the corresponding functions with the usual approach, however a change is necessary to account for the fact that log is not defined in 0. My proposal is:
np.log(mu24) -> np.max([ np.log(mu24), min_log ])
np.log(1.0 - rsq) -> np.max([ np.log(1.0 - rsq), 2*min_log ])
Where min_log could depend on precision.
Below I report the corresponding functions with the usual approach, however a change is necessary to account for the fact that log is not defined in 0. My proposal is: np.log(mu24) -> np.max([ np.log(mu24), min_log ]) np.log(1.0 - rsq) -> np.max([ np.log(1.0 - rsq), 2*min_log ]) Where min_log could depend on precision.
integral definitions for classic_log method
def integral_r_classic_log(limb_darkening_coefficients, r): a1, a2 = limb_darkening_coefficients rsq = r r mu44 = 1.0 - rsq mu24 = np.sqrt(mu44) return 0.5 (a1 - 1.0) mu44 - (a1/3.0 + a2/9.0) mu44 mu24 + (a2/3.0) mu44 mu24 np.log(mu24)
def num_classic_log(r, limb_darkening_coefficients, rprs, z): a1, a2 = limb_darkening_coefficients rsq = r * r
mu = np.sqrt(1.0 - rsq)
def integral_r_f_classic_log(limb_darkening_coefficients, rprs, z, r1, r2, precision=3): x1 = (r2 - r1) / 2.0 x2 = (r2 + r1) / 2.0
return x1 np.sum(gausstab[precision][0][:, None] num_classic_log(x1[None, :] * gausstab[precision][1][:, None] + x2[None, :], limb_darkening_coefficients, rprs, z), 0)