Merge branch 'limb_darkening'

- Add limb darkening effect. Currently only suppurt the linear limb darkening law and not in an adaptive scheme.
- Add tests for contour integral and limb darkening effect.
- Update docstring format.
- Rename some functions and refactor the file structures.

Author: CoastEgo

docs/api/model.md

@@ -1,10 +1,13 @@
 # Magnification Model
 
 Light curve for the binary lens with adaptive contour integration and point sorce approximation.
-::: microlux.extended_light_curve
+::: microlux.binary_mag
 ---
 Light curve with point source approximation.
 ::: microlux.point_light_curve
 ---
+Adaptive contour integration for the light curve.
+::: microlux.extended_light_curve
+---
 Adaptive contour integration for the single epoch.
-::: microlux.contour_integral
+::: microlux.contour_integral

example/KB0371.ipynb

@@ -209,7 +209,7 @@
     "import numpyro.distributions as dist\n",
     "from numpyro.infer import MCMC, NUTS\n",
     "from numpyro.diagnostics import print_summary\n",
-    "from BinaryJax import model\n",
+    "from microlux import binary_mag\n",
     "\n",
     "import VBBinaryLensing\n",
     "from MulensModel import caustics\n",
@@ -699,7 +699,7 @@
     "    alpha_deg = parms[4]\n",
     "    s = 10.**parms[5]\n",
     "    q = 10.**parms[6]\n",
-    "    mag_Jax = model(t0, u0, tE, rho, q, s, alpha_deg, times)\n",
+    "    mag_Jax = binary_mag(t0, u0, tE, rho, q, s, alpha_deg, times)\n",
     "    return mag_Jax\n",
     "initial_guess = [8.59238794e+03, 1.42915228e-01, 6.61567944e+00, -2.23131913e+00, 2.71714918e+02, -7.73128397e-02, -1.14229367e+00]\n",
     "times,flux,ferr = HJD,flux,ferr\n",
@@ -788,7 +788,7 @@
     "    alpha_deg = parms[4]\n",
     "    s = 10.**parms[5]\n",
     "    q = 10.**parms[6]\n",
-    "    mag_Jax = model(t0, u0, tE, rho, q, s, alpha_deg, times)\n",
+    "    mag_Jax = binary_mag(t0, u0, tE, rho, q, s, alpha_deg, times)\n",
     "    return mag_Jax\n",
     "def light_curve_Jax_pmap(times,parms,i):\n",
     "    times = jnp.reshape(times,(-1,N_pmap),order='C')\n",

mkdocs.yml

@@ -74,7 +74,7 @@ plugins:
                     show_root_heading: true
                     show_root_full_path: true
                     show_if_no_docstring: false
-                    show_signature_annotations: false
+                    show_signature_annotations: true
                     members_order: source
                     heading_level: 2
                 # rendering:

src/microlux/__init__.py

@@ -1,8 +1,9 @@
 # # -*- coding: utf-8 -*-
 all = [
-    "model",
     "point_light_curve",
+    "extended_light_curve",
     "contour_integral",
+    "binary_mag",
     "Iterative_State",
     "Error_State",
     "to_lowmass",
@@ -13,8 +14,9 @@
     to_centroid as to_centroid,
     to_lowmass as to_lowmass,
 )
+from .countour import contour_integral as contour_integral
 from .model import (
-    contour_integral as contour_integral,
+    binary_mag as binary_mag,
     extended_light_curve as extended_light_curve,
     point_light_curve as point_light_curve,
 )

src/microlux/basic_function.py

@@ -38,6 +38,11 @@ def to_lowmass(s, q, x):
 
 
 def Quadrupole_test(rho, s, q, zeta, z, cond, tol=1e-2):
+    """
+    The quadrupole test, ghost image test, and planetary caustic test proposed by Bozza 2010 to check the validity of the point source approximation.
+    The coefficients are fine-tuned in our implementation.
+
+    """
     m1 = 1 / (1 + q)
     m2 = q / (1 + q)
     cQ = 2
@@ -93,6 +98,9 @@ def Quadrupole_test(rho, s, q, zeta, z, cond, tol=1e-2):
 
 
 def get_poly_coff(zeta_l, s, m2):
+    """
+    get the polynomial cofficients of the polynomial equation of the lens equation. The low mass object is at the origin and the primary is at s.
+    """
     zeta_conj = jnp.conj(zeta_l)
     c0 = s**2 * zeta_l * m2**2
     c1 = -s * m2 * (2 * zeta_l + s * (-1 + s * zeta_l - 2 * zeta_l * zeta_conj + m2))
@@ -140,6 +148,11 @@ def dot_product(a, b):
 
 
 def basic_partial(z, theta, rho, q, s, caustic_crossing):
+    """
+
+    basic partial derivatives of the lens equation with respect to zeta, z, and theta used in the error estimation.
+
+    """
     z_c = jnp.conj(z)
     parZetaConZ = 1 / (1 + q) * (1 / (z_c - s) ** 2 + q / z_c**2)
     par2ConZetaZ = -2 / (1 + q) * (1 / (z - s) ** 3 + q / (z) ** 3)
@@ -198,8 +211,8 @@ def refine_gradient(zeta_l, q, s, z):
 @refine_gradient.defjvp
 def refine_gradient_jvp(primals, tangents):
     """
-    use the custom jvp to refine the gradient of roots respect to zeta_l, based on the equation on V.Bozza 2010 eq 20.
-    The necessity of this function is still under investigation.
+    use the custom jvp to refine the gradient of roots respect to zeta_l, based on the equation on V.Bozza 2010 eq 20 and also see our paper for the details
+    This will simplify the computational graph and accelerate the gradient calculation
     """
     zeta, q, s, z = primals
     tangent_zeta, tangent_q, tangent_s, tangent_z = tangents