{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":101849,"databundleVersionId":13093295,"sourceType":"competition"}],"dockerImageVersionId":31089,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Limb darkening\n\nThis notebook explains how ExoSim2 computes limb darkening.\nStar is dimmer near the edge (or limb) than the center.\n\n\n[ExoSim2](https://github.com/arielmission-space/ExoSim2-public)\nis using [batman](https://github.com/lkreidberg/batman) for transit curve calculation:\n\n```python\nm = batman.TransitModel(params, timeline, transittype=transittype)\n```\n\nsrc/exosim/tasks/astrosignal/estimatePlanetarySignal.py\n\n\nSee also the original code: `calc_limb_darkening()`, `area()`, `intensity()` in\n`batman/c_src/_nonlinear_ld.c`, and `common.h`.\n\n### Remark\n\n* I am not computing the Kepler orbit, just using straight lines in front of the star.\n* This is how the data are generated, but not necessarily how you should fit the data. The computation could be too complicated and slow, and you can use simpler and faster formula for inference.\n\n**Version 2 and later:**\n\nFix minor texts.","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2025-08-01T21:40:54.408929Z","iopub.execute_input":"2025-08-01T21:40:54.409336Z","iopub.status.idle":"2025-08-01T21:40:54.418205Z","shell.execute_reply.started":"2025-08-01T21:40:54.409307Z","shell.execute_reply":"2025-08-01T21:40:54.416675Z"}}},{"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\nimport math\n\nRs = 1     # radius of the star, i.e., lengths are in units of the stellar radius\nRp = 0.3   # radius of the planet","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-10T22:31:25.899375Z","iopub.execute_input":"2025-08-10T22:31:25.899791Z","iopub.status.idle":"2025-08-10T22:31:25.907506Z","shell.execute_reply.started":"2025-08-10T22:31:25.899755Z","shell.execute_reply":"2025-08-10T22:31:25.90585Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Limb darkening coefficients\n\n\nThere are several empirical formulae for the limb darkening,\nand I do not know which one is used for this competition,\nbut [Claret (2000)](https://ui.adsabs.harvard.edu/abs/2000A%26A...363.1081C)\nis one of the most flexible one:\n\n$$\nI(\\mu) / I_1 = 1 - a_1(1 - \\mu^{1/2}) - a_2(1 - \\mu) - a_3(1 - \\mu^{3/2}) - a_4(1 - \\mu^2),\n$$\n\nwhere $\\mu = \\cos \\theta$, and $\\theta$ is the angle between the stellar surface normal and the line of sight.\n\n","metadata":{}},{"cell_type":"code","source":"# Explain mu\n# written by ChatGPT\n\nR = 1                   # radius of the star\ntheta = np.radians(45)\nL1 = 0.5                # length normal\nL2 = 0.5                # length tangent\n\n# Figure\nfig, ax = plt.subplots()\norigin = np.array([0, 0])\n\n# Half circle\ntheta_vals = np.linspace(np.pi, 2 * np.pi, 300)\nx_circle = R * np.cos(theta_vals)\ny_circle = R * np.sin(theta_vals)\nax.plot(x_circle, y_circle, 'k')\n\n# Angle theta\nxP = R * np.cos(theta)\nyP = -R * np.sin(theta)\nax.plot([0, xP], [0, yP], 'b')\n\n# Line of sight\nax.plot([xP, xP], [yP, yP - L1], 'g')\n\n# Tangent space\ndx = -np.sin(theta)\ndy = -np.cos(theta)\nxL = xP + L2 * dx\nyL = yP + L2 * dy\nxR = xP - L2 * dx\nyR = yP - L2 * dy\nax.plot([xL, xR], [yL, yR], 'r')\n\n# Normal vector\nn_len = 0.5\nxn = xP + n_len * np.cos(theta)\nyn = yP - n_len * np.sin(theta)\nax.annotate('', xy=(xn, yn), xytext=(xP, yP),\n            arrowprops=dict(arrowstyle='->', color='black'))\nax.text(xn + 0.05, yn, 'n', fontsize=12)\n\n# Annotate theta\n# You can also use patches.Arc. I would not do this, ChatGPT is doing.\narc_theta = np.linspace(-np.pi / 2, -np.pi / 2 + theta, 100)\nr_arc = 0.3\nx_arc = xP + r_arc * np.cos(arc_theta)\ny_arc = yP + r_arc * np.sin(arc_theta)\nax.plot(x_arc, y_arc, 'k')\n\ndx = 0.1\ntheta_label_angle = -np.pi / 2 + theta / 2\nx_theta_text = xP + (r_arc + dx) * np.cos(theta_label_angle)\ny_theta_text = yP + (r_arc + dx) * np.sin(theta_label_angle)\nax.text(x_theta_text, y_theta_text, r'$\\theta$', fontsize=12, ha='center', va='center')\n\n# Annotate line of sight\nlos_arrow_len = 0.2\nlos_start = yP - L1 - 0.3\nlos_end = los_start + los_arrow_len\nax.annotate('', xy=(xP, los_end), xytext=(xP, los_start),\n            arrowprops=dict(arrowstyle='->', color='darkgreen'))\n\nax.text(xP + 0.05, los_start + los_arrow_len / 2, 'line of sight', \n        fontsize=10, verticalalignment='center', color='darkgreen')\n\n# Axis lines\nx_arrow_len = 1.3\nz_arrow_len = 1.3\narrow_style = dict(arrowstyle='->', color='lightgray', linewidth=1.5)\n\n# x axis\nax.annotate('', xy=(x_arrow_len, 0), xytext=(0, 0), arrowprops=arrow_style)\nax.text(x_arrow_len + 0.05, 0, 'x', fontsize=12, color='gray', va='center')\n\n# z axis (line of sight)\nz_arrow_start = -R - 0.2\nz_arrow_end = 0.2\nax.annotate('', xy=(0, z_arrow_end), xytext=(0, z_arrow_start),\n            arrowprops=arrow_style)\n\nax.text(0, z_arrow_end + 0.05, 'z', fontsize=12, color='gray', ha='center')\n\n# Plot range\nax.set_aspect('equal')\nax.set_xlim(-1.5, 1.5)\nax.set_ylim(-1.5, 0.5)\nax.axis('off')\n\nplt.show()","metadata":{"trusted":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2025-08-10T22:31:25.9091Z","iopub.execute_input":"2025-08-10T22:31:25.90947Z","iopub.status.idle":"2025-08-10T22:31:26.041901Z","shell.execute_reply.started":"2025-08-10T22:31:25.909434Z","shell.execute_reply":"2025-08-10T22:31:26.040966Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"In the sky, i.e., the projection to xy plane,\nthe distance from the center of star to a point in star, $d$,\nis related to $\\mu$ by\n\n$$ d / R_s = \\sin \\theta, $$\n\nwhere $R_s$ is the radius of the star.\n\n$$ \\mu^2 = 1 - (d / R_s)^2 $$\n\n","metadata":{}},{"cell_type":"markdown","source":"# Step 1: Planet size effect without limb darkening\n\nBefore computing limb darking, I compute the edge effect when the planet eclipses the star partially near the\nedge of the star.","metadata":{"execution":{"iopub.status.busy":"2025-08-01T21:39:55.164537Z","iopub.execute_input":"2025-08-01T21:39:55.164968Z","iopub.status.idle":"2025-08-01T21:39:55.173833Z","shell.execute_reply.started":"2025-08-01T21:39:55.164942Z","shell.execute_reply":"2025-08-01T21:39:55.171976Z"}}},{"cell_type":"markdown","source":"\n","metadata":{}},{"cell_type":"code","source":"R = Rs    # radius of the star\nr = Rp    # radius of the planet\nx = -0.8   # position of the star\n\nfig, ax = plt.subplots()\n\n# Parameter for drawing circle\n# Again, ChatGPT is doing this, instead of patches.Circle\ntheta_vals = np.linspace(0, 2 * np.pi, 300)\n\nx_circle = R * np.cos(theta_vals)\ny_circle = R * np.sin(theta_vals)\nplt.plot(x_circle, y_circle, color='orange')\n\n# Circle r (planet)\nx_circle = x + r * np.cos(theta_vals)\ny_circle = r * np.sin(theta_vals)\nax.plot(x_circle, y_circle, color='black')\n\n# Length x\narrow_style = dict(arrowstyle='-', color='blue', linewidth=1.5, shrinkA=0, shrinkB=0)\nplt.annotate('', xy=(x, 0), xytext=(0, 0), arrowprops=arrow_style)\n\nplt.text(x / 2, -0.15, 'd', fontsize=12, color='blue', ha='center')\n\n# Length R\ncost = (x**2 + R**2 - r**2) / (2 * x * R)\nsint = math.sqrt(1 - cost ** 2)\n\narrow_style = dict(arrowstyle='-', color='green', linewidth=1.5)\nplt.annotate('', xy=(cost, sint), xytext=(0, 0), arrowprops=arrow_style)\nplt.text(x / 2, 0.3, '$R_s$', fontsize=12, color='green', ha='center', va='center')\n\n# Length r\narrow_style = dict(arrowstyle='-', color='red', linewidth=1.5, shrinkA=0, shrinkB=0)\nplt.annotate('', xy=(cost, sint), xytext=(x, 0), arrowprops=arrow_style)\nplt.text(x -0.12 , 0.0, '$R_p$', fontsize=12, color='red', ha='center', va='center')\n\n# Angle alpha\ntheta = math.atan2(sint, cost)\narc_theta = np.linspace(theta, np.pi, 100)\nr_arc = 0.3\nx_arc = r_arc * np.cos(arc_theta)\ny_arc = r_arc * np.sin(arc_theta)\nax.plot(x_arc, y_arc, 'cyan')\n\n# Angle beta\ncos_beta = (x**2 + r**2 - R**2) / abs(2 * x * r)\nsin_beta = math.sqrt(1 - cos_beta ** 2)\nbeta = math.atan2(sin_beta, cos_beta)\narc_theta = np.linspace(0, beta, 100)\nr_arc = 0.1\nx_arc = r_arc * np.cos(arc_theta)\ny_arc = r_arc * np.sin(arc_theta)\nax.plot(x + x_arc, y_arc, 'magenta')\n\nplt.text(0, 0.7, 'star', fontsize=12, color='orange', ha='center', va='center')\nplt.text(x, -0.4, 'planet', fontsize=12, color='black', ha='left', va='center')\nplt.text(0, 0, '$\\\\alpha$', fontsize=12, color='cyan', ha='left', va='center')\nplt.text(x, -0.12, '$\\\\beta$', fontsize=12, color='magenta', ha='center', va='center')\n\n\n# Figure range\nax.set_aspect('equal')\nax.set_xlim(-1.5, 1.5)\nax.set_ylim(-1.5, 1.5)\nax.axis('off')\n\nplt.show()","metadata":{"trusted":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2025-08-10T22:31:26.042602Z","iopub.execute_input":"2025-08-10T22:31:26.042858Z","iopub.status.idle":"2025-08-10T22:31:26.159837Z","shell.execute_reply.started":"2025-08-10T22:31:26.042834Z","shell.execute_reply":"2025-08-10T22:31:26.158751Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"The area of two intesecting disks is:\n\n1. Sector of angle $\\beta$ and radius $R_p$, which is entirely inside the stellar cicle,\n2. add the remaining part intersecting with the stellar edge,\n    - sector of angle $\\alpha$ and radius $R_s$,\n    - subtract the triangle of lengths $R_s, R_p, d$,\n3. times 2 for the other half.\n\n$$\nA = \\beta R_p^2 + \\alpha R_s^2 - d R_s \\sin \\alpha\n$$","metadata":{}},{"cell_type":"code","source":"def area(R, r, d):\n    \"\"\"\n\tCompute the area of overlapping disks with radii R and r; separated by a distance d.\n    Pure mathematics.\n\n    Args:\n      R, r (array): Radii of two disks\n      x (array): Distance betwen two centers\n\n    Three arrays must have the same shape\n    \"\"\"\n    assert R.shape == r.shape == d.shape\n\n    d = np.abs(d)  # distance (d >= 0)\n\n    # Case 1. Default: no overlap\n    out = np.zeros_like(d)\n\n    # Case 2. Planet completely inside the sterllar circle\n    idx_inside = d <= np.abs(R - r)\n    out[idx_inside] = np.pi * np.minimum(r, R)[idx_inside] ** 2\n    \n    # Case 3. Star and planet are overlapping\n    idx = np.logical_and(np.abs(R - r) < d, d < R + r)\n    d = d[idx]\n    r = r[idx]\n    R = R[idx]\n\n    # Law of cosines\n    cos_alpha = (d ** 2 + R ** 2 - r ** 2) / (2 * R * d)\n    cos_beta = (r ** 2 + d ** 2 - R ** 2) / (2 * r * d)\n\n    assert np.all(np.abs(cos_alpha) < 1.00001)\n    assert np.all(np.abs(cos_beta) < 1.00001)\n    cos_alpha = np.clip(cos_alpha, -1, 1)\n    cos_beta = np.clip(cos_beta, -1, 1)\n\n    # overlap = r x sin(beta) / 2\n    overlap_sq = 2 * (R ** 2 * r ** 2 + r ** 2 * d ** 2 + d ** 2 * R ** 2) - (R ** 4 + r ** 4 + d ** 4)\n    overlap_sq = np.clip(overlap_sq, 0, None)\n\n    out[idx] = R ** 2 * np.arccos(cos_alpha) + r ** 2 * np.arccos(cos_beta) - 0.5 * np.sqrt(overlap_sq)\n    return out\n\n\nd = np.linspace(-2, 2, 201)\nR = np.ones_like(d)\nr = Rp * np.ones_like(d)\n\na = area(R, r, d)\ndelta = 1 - a / math.pi  # Fraction compared to star area\n\nplt.figure(figsize=(6, 3))\nplt.ylabel('Relative flux')\nplt.title('No limb darkening; r/R=0.3')\nplt.plot(d, delta, label='b=0')\nplt.text(0, 0.92, '↓ Flat', ha='center')\n\n# delta\ndelta1 = (Rp / Rs) ** 2\nplt.axhline(1 - delta1, ls=':', color='gray', alpha=0.5)\narrow_style = dict(arrowstyle='<->', color='red', linewidth=1.5, shrinkA=0, shrinkB=0)\n\nplt.annotate('', xy=(1.5, 1), xytext=(1.5, 1 - delta1), arrowprops=arrow_style)\nplt.text(1.55, 1 - delta1 / 2, '$\\\\delta$', va='center', fontsize=14)\n\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-10T22:31:26.161023Z","iopub.execute_input":"2025-08-10T22:31:26.161354Z","iopub.status.idle":"2025-08-10T22:31:26.369773Z","shell.execute_reply.started":"2025-08-10T22:31:26.161326Z","shell.execute_reply":"2025-08-10T22:31:26.368905Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Without limb darkening, the light curve is flat when the planet is entirely inside the stellar circle.\n\nBelow, I introduce the impact parameter $b$ for a linear trajectory that does not pass the stellar center.","metadata":{"execution":{"iopub.status.busy":"2025-08-01T22:21:50.322491Z","iopub.execute_input":"2025-08-01T22:21:50.322812Z","iopub.status.idle":"2025-08-01T22:21:50.329883Z","shell.execute_reply.started":"2025-08-01T22:21:50.322791Z","shell.execute_reply":"2025-08-01T22:21:50.328317Z"}}},{"cell_type":"code","source":"# Visualize impact parameter b\nR = Rs\nr = Rp\ns = -0.5\nb = -0.6\n\nfig, ax = plt.subplots()\n\ntheta_vals = np.linspace(0, 2 * np.pi, 300)\n\n# Star\nx_circle = R * np.cos(theta_vals)\ny_circle = R * np.sin(theta_vals)\nplt.plot(x_circle, y_circle, color='orange')\n\n# Planet\nx_circle = s + r * np.cos(theta_vals)\ny_circle = b + r * np.sin(theta_vals)\nax.plot(x_circle, y_circle, color='black')\n\n# Orbit (linear approximation)\narrow_style = dict(arrowstyle='<-', color='gray', linewidth=1.5, shrinkA=0, shrinkB=0)\nplt.annotate('s', xy=(-1.5, b), xytext=(1.5, b), arrowprops=arrow_style)\nplt.text(s, b - 0.1, '(s, b)', fontsize=12, color='gray', ha='center', va='center')\n\narrow_style = dict(arrowstyle='-', color='red', linewidth=1.5, shrinkB=0)\nplt.annotate('', xy=(0, 0), xytext=(0, b), arrowprops=arrow_style)\nplt.text(0.04, b / 2, 'b', fontsize=12, color='red', ha='left', va='center')\n\narrow_style = dict(arrowstyle='-', color='blue', linewidth=1.5, shrinkB=0)\nplt.annotate('', xy=(0, 0), xytext=(s, b), arrowprops=arrow_style)\nplt.text(s / 2 - 0.1, b / 2 + 0.06, 'd', fontsize=12, color='blue', ha='center', va='center')\n\n# Figure\nplt.title('Impact parameter b: $d^2 = s^2 + b^2$')\nax.set_aspect('equal')\nax.set_xlim(-1.5, 1.5)\nax.set_ylim(-1.5, 1.5)\nax.axis('off')\n\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-10T22:31:26.370556Z","iopub.execute_input":"2025-08-10T22:31:26.370809Z","iopub.status.idle":"2025-08-10T22:31:26.527004Z","shell.execute_reply.started":"2025-08-10T22:31:26.37079Z","shell.execute_reply":"2025-08-10T22:31:26.525803Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.figure(figsize=(6, 3))\nplt.ylabel('Relative flux')\nplt.title('Impact paramter b for off-center trajectory')\n\nd = np.linspace(-2, 2, 201)\nR = np.ones_like(d)\nr = Rp * np.ones_like(d)\n\na = area(R, r, d)\nplt.plot(d, 1 - a / math.pi, label='b=0')\n\n# Add impact paramber b\nb = 0.6\ns = np.linspace(-2, 2, 201)  # coordinate along the planetary orbit (linear approximation)\nxx_b = np.sqrt(b ** 2 + d ** 2) * np.sign(d)\na_b = area(R, r, xx_b)\n\nplt.plot(s, 1 - a_b / math.pi, '--', label=('b=%.1f' % b))\nplt.xlabel('s')\nplt.legend()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-10T22:31:26.527753Z","iopub.execute_input":"2025-08-10T22:31:26.528Z","iopub.status.idle":"2025-08-10T22:31:26.739822Z","shell.execute_reply.started":"2025-08-10T22:31:26.527981Z","shell.execute_reply":"2025-08-10T22:31:26.738821Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# 2. Limb darkening function\n\nThe intensity function is normalized by the total.\n\n$$\n\\int_0^{R_s} I \\, 2\\pi r dr = 2 \\pi R_s^2 \\int_0^1 I(\\mu) \\mu d\\mu = \\mathcal{N}\n$$\n\nusing,\n\n$$ (r / R_s)^2 = 1 - \\mu^2, $$\n\nas mentioned near the begenning. We can integrate analytically,\n\n$$\n\\int_0^1 \\mu^{k/2} \\mu d\\mu = \\frac{1}{2 + k /2}.\n$$","metadata":{}},{"cell_type":"code","source":"def intensity(d, c):\n    \"\"\"\n    Compute the stelar intensity with Claret 4-parameter non-linear\n    \n    Args:\n      d (array): Distance between stellar center and planet center in units of stellar radius Rs=1\n      c (4, ): Limb darkening coefficent, 4-parameter nonlinear + normalization\n    \"\"\"\n    assert len(c) == 4\n\n    # Overall normalization, such that, the total intensity is 1\n    norm = (0.5 - c[0] / 10 - c[1] / 6 - 3 * c[2] / 14 - c[3] / 4) * 2 * np.pi\n    \n    # From batmann batman/c_src/_nonlinear_ld.c; to avoid negative in sqrt?\n    d = np.minimum(d, 0.99995)\n\n    sqrtmu = (1 - d ** 2) ** 0.25  # d = Rs sinθ, Rs = 1, mu = cos θ -> mu**2 = 1 - d**2\n    \n    return (1 - c[0] * (1 - sqrtmu)\n              - c[1] * (1 - sqrtmu ** 2)\n              - c[2] * (1 - sqrtmu ** 3)\n              - c[3] * (1 - sqrtmu ** 4)) / norm\n\n\nc = (0.4, 0.3, 0.2, 0.1)  # I don't know what typcial coefficients are\nd = np.cos(np.linspace(0, 1, 101))\nll = intensity(d, c)\n\nplt.figure(figsize=(6, 3))\nplt.title('Limb darkening function')\nplt.xlabel('$d / R_s$')\nplt.ylabel('$I(\\mu) / \\\\mathcal{N}$')\nplt.plot(d, ll)\nplt.show()\n\nprint('Limb darkening coefficients:', c)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-10T22:31:26.740847Z","iopub.execute_input":"2025-08-10T22:31:26.741195Z","iopub.status.idle":"2025-08-10T22:31:26.929916Z","shell.execute_reply.started":"2025-08-10T22:31:26.741172Z","shell.execute_reply":"2025-08-10T22:31:26.929109Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Check total\nc = (0.4, 0.3, 0.2, 0.1)\nr = np.linspace(0, 1, 1001)\ndr = np.diff(r)\nr_mid = 0.5 * (r[:-1] + r[1:])\nI = intensity(r_mid, c)\n\nA = np.pi * r ** 2\ndA = np.diff(A)\n\ntotal = np.sum(I * dA)\nprint('Check total:', total)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-10T22:31:26.931621Z","iopub.execute_input":"2025-08-10T22:31:26.93188Z","iopub.status.idle":"2025-08-10T22:31:26.939466Z","shell.execute_reply.started":"2025-08-10T22:31:26.931859Z","shell.execute_reply":"2025-08-10T22:31:26.938554Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Integrate the blocked intensity\n\nThis is the final step: Compute the intensity blocked by the planet.\n\nThe integration variable $r$ is the radius from the stellar center,\nbecause the intensity $I(\\mu)$ is a function of $r$.\n\n$$\n\\delta = \\int_\\mathrm{overlap} I\\, 2\\pi r dr \\approx \\sum_n I_{n + 1/2} \\Delta A_n\n$$\n\n**Because the planet is usually much smaller than the star, I am not sure such integration is necessary.**","metadata":{"execution":{"iopub.status.busy":"2025-08-01T23:01:10.69831Z","iopub.execute_input":"2025-08-01T23:01:10.698715Z","iopub.status.idle":"2025-08-01T23:01:10.706615Z","shell.execute_reply.started":"2025-08-01T23:01:10.698688Z","shell.execute_reply":"2025-08-01T23:01:10.70514Z"}}},{"cell_type":"code","source":"R = Rs    # radius of the star\nr = Rp    # radius of the planet\nx = -0.8   # position of the star\n\nfig, ax = plt.subplots()\n\n# The for circle drawing\ntheta_vals = np.linspace(0, 2 * np.pi, 300)\n\nx_circle = R * np.cos(theta_vals)\ny_circle = R * np.sin(theta_vals)\nplt.plot(x_circle, y_circle, color='orange')\n\n# Circle r (planet)\nx_circle = x + r * np.cos(theta_vals)\ny_circle = r * np.sin(theta_vals)\nax.plot(x_circle, y_circle, color='black')\n\nplt.text(0, 1.15, 'star', fontsize=12, color='orange', ha='center', va='center')\nplt.text(x, -0.5, 'planet', fontsize=12, color='black', ha='center', va='center')\n\n# Integration circles\nx_circle = np.cos(theta_vals)\ny_circle = np.sin(theta_vals)\nplt.plot(0.65 * x_circle, 0.65 * y_circle, color='lightgray')\nplt.plot(0.75 * x_circle, 0.75 * y_circle, color='lightgray')\n\nplt.text(-0.45, 0.3, '$r_{n}$', fontsize=12, color='gray', ha='center', va='center')\nplt.text(-0.70, 0.5, '$r_{n+1}$', fontsize=12, color='gray', ha='center', va='center')\nplt.text(-0.70, 0.0, '$\\\\Delta A_{n}$', fontsize=12, color='black', ha='center', va='center')\n\n# Figure range\nax.set_aspect('equal')\nax.set_xlim(-1.5, 1.5)\nax.set_ylim(-1.5, 1.5)\nax.axis('off')\n\nplt.show()\n\nprint('Integrate the radius-dependent intensity on star-planet overlap.')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-10T22:31:26.940555Z","iopub.execute_input":"2025-08-10T22:31:26.941133Z","iopub.status.idle":"2025-08-10T22:31:27.03975Z","shell.execute_reply.started":"2025-08-10T22:31:26.941107Z","shell.execute_reply":"2025-08-10T22:31:27.038784Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def limb_darkening(d, rp, c, *, nstep=101):\n    \"\"\"\n    Compute the intensity from the star\n    \n    d  (float): Distance between stellar center and planet center in units of stellar radius\n    rp (float): Planet radius\n    c  (array): Limb darkening coefficients (4, )\n    \"\"\"\n    assert isinstance(d, float)\n    assert isinstance(rp, float)\n    \n    # Range of integration\n    r_min = max(d - rp, 0)     # lower bound for integration\n    r_max = min(d + rp, 1.0)   # upper bound for integration\n\n    if r_min >= 1:\n        return 1  # No overlap between star and planet\n    elif r_max - r_min < 1e-7:\n        return 1  # Overlap is too small\n\n    # This is different from the batman code.\n    # See the original code for accurate error control.\n    r = np.linspace(r_min, r_max, nstep)\n    dr = np.diff(r)\n    r_mid = r[:-1] + dr\n\n    rp_array = np.full_like(r, rp)\n    d_array = np.full_like(r, d)\n    A = area(r, rp_array, d_array)\n\n    dA = np.diff(A)\n    I = intensity(r_mid, c)\n\n    assert len(dA) == len(I)\n\n    # \\int I(x) dA\n    integ = np.sum(I * dA)  # intensity blocked by the planet\n\n    return 1 - integ  # intensity not blocked\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-10T22:31:27.040795Z","iopub.execute_input":"2025-08-10T22:31:27.041113Z","iopub.status.idle":"2025-08-10T22:31:27.048644Z","shell.execute_reply.started":"2025-08-10T22:31:27.041079Z","shell.execute_reply":"2025-08-10T22:31:27.047674Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def plot_limb_darkening(rp, b, c, *, title=''):\n    \"\"\"\n    rp (float): Planet radius in units of stellar radius (Rs=1)\n    \"\"\"\n    # = 0.1  # imbact parameter\n    s = np.linspace(-2, 2, 201)  # parameter along the planetary orbit (linear)\n    d = np.sqrt(b ** 2 + s ** 2) #* np.sign(s)\n    R = np.ones_like(d)\n    r = np.full_like(d, rp)\n\n    # Area of star-planet overlap\n    area_b = area(R, r, d)\n    \n    ll = []\n    for d1 in d:\n        l = limb_darkening(d1, rp, c)\n        ll.append(l)\n\n    plt.figure(figsize=(6, 2))\n    plt.title(title)\n    plt.xlabel('$s / R_s$')\n    plt.ylabel('Relative flux')\n    plt.plot(s, 1 - area_b / math.pi, label='uniform')\n    plt.plot(s, ll)\n    plt.axhline(1 - rp ** 2, ls=':', color='gray', alpha=0.5)\n    plt.show()\n\nb = 0.1\nc = (0.4, 0.3, 0.2, 0.1)\nplot_limb_darkening(Rp, b, c, title='$R_p/R_s=0.3, b/R_s=0.1$')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-10T22:31:27.050052Z","iopub.execute_input":"2025-08-10T22:31:27.050414Z","iopub.status.idle":"2025-08-10T22:31:27.37333Z","shell.execute_reply.started":"2025-08-10T22:31:27.050383Z","shell.execute_reply":"2025-08-10T22:31:27.37242Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Limb darkening makes the light curve round, because the planet blocks more light when the planet passes the brighter part of the star nearer to the center.","metadata":{}},{"cell_type":"code","source":"# Very large Rp/Rs = 0.3 is for visualization purpose.\n# Rp is much smaller in the competition\nb = 0.1\nc = (0.1, 0.1, 0.1, 0.1)\nplot_limb_darkening(0.03, b, c, title='$R_p/R_s=0.03, b/R_s=0.1$')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-10T22:31:27.374229Z","iopub.execute_input":"2025-08-10T22:31:27.374507Z","iopub.status.idle":"2025-08-10T22:31:27.607589Z","shell.execute_reply.started":"2025-08-10T22:31:27.374487Z","shell.execute_reply":"2025-08-10T22:31:27.606644Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Note that the intensity is normalized by the total.\nThe minimum delta is not $1 - (R_p/R_s)^2$ with limb darkening.","metadata":{}},{"cell_type":"code","source":"# No limb darkening for c=0\nb = 0.1\nc = (0, 0, 0, 0)\nplot_limb_darkening(Rp, b, c, title='Check c=0 case')\n\nprint('Area calculation agree when no limb darkening c=0.')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-10T22:31:27.61132Z","iopub.execute_input":"2025-08-10T22:31:27.611585Z","iopub.status.idle":"2025-08-10T22:31:27.821397Z","shell.execute_reply.started":"2025-08-10T22:31:27.611564Z","shell.execute_reply":"2025-08-10T22:31:27.820552Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Reference\n\nExoSim2: https://github.com/arielmission-space/ExoSim2-public\n\nbatman: https://github.com/lkreidberg/batman\n\nClaret (2000): https://ui.adsabs.harvard.edu/abs/2000A%26A...363.1081C\n","metadata":{}}]}