<jupyter_start><jupyter_text>Derivative of $P(G_n)$ with respect to $r$: The derivative of the integrand of P(G_n) with respect to rcompared with the analytic formula for the derivative:<jupyter_code>import numpy as np
import matplotlib.pyplot as pl
%matplotlib inline
from scipy.integrate import quad
import warnings
warnings.filterwarnings("ignore")<jupyter_output><empty_output><jupyter_text>Let's define the numerical versions of the $P$ and $M$ integrals:<jupyter_code>def P(n, r, b):
    ksq = (1 - r ** 2 - b ** 2 + 2 * b * r) / (4 * b * r)
    if ksq <= 1:
        kappa = 2 * np.arcsin(np.sqrt(ksq))
    else:
        kappa = np.pi
    def integrand(xi):
        return (ksq - np.sin(xi) ** 2) ** (n / 2.0) * (r - b + 2 * b * np.sin(xi) ** 2)
    res, err = quad(integrand, -kappa / 2, kappa / 2)
    return res * 2 * r * (4 * b * r) ** (n / 2.0)
def M(n, r, b):
    ksq = (1 - r ** 2 - b ** 2 + 2 * b * r) / (4 * b * r)
    if ksq <= 1:
        kappa = 2 * np.arcsin(np.sqrt(ksq))
    else:
        kappa = np.pi
    def integrand(xi):
        return (ksq - np.sin(xi) ** 2) ** (n / 2.0)
    res, err = quad(integrand, -kappa / 2, kappa / 2)
    return res * (4 * b * r) ** (n / 2.0)<jupyter_output><empty_output><jupyter_text>Define the numerical derivatives:<jupyter_code>def dPdr(n, r, b, eps=1.e-8):
    return (P(n, r + eps, b) - P(n, r - eps, b)) / (2 * eps)
def dPdb(n, r, b, eps=1.e-8):
    return (P(n, r, b + eps) - P(n, r, b - eps)) / (2 * eps)<jupyter_output><empty_output><jupyter_text>Define the equations in the paper:<jupyter_code>def dPdr_paper(n, r, b):
    return 2 * r * ((n + 2) * M(n, r, b) - n * M(n - 2, r, b))
def dPdb_paper(n, r, b):
    return (n / b) * ((r ** 2 + b ** 2) * (M(n, r, b) - M(n - 2, r, b)) + (r ** 2 - b ** 2) ** 2 * M(n - 2, r, b))<jupyter_output><empty_output><jupyter_text>Let's plot the difference on a grid:<jupyter_code>def plot(n):
    barr = np.linspace(1e-4, 2, 30)
    rarr = np.linspace(1e-4, 2, 30)
    diff_b = np.zeros((30, 30))
    diff_r = np.zeros((30, 30))
    for i, b in enumerate(barr):
        for j, r in enumerate(rarr):
            diff_r[j, i] = abs(dPdr(n, r, b) - dPdr_paper(n, r, b))
            diff_b[j, i] = abs(dPdb(n, r, b) - dPdb_paper(n, r, b)) 
    fig, ax = pl.subplots(1, 2, figsize=(14, 5))
    ax[0].imshow(np.log10(diff_r), origin="lower", extent=(0, 2, 0, 2), vmin=-15, vmax=0)
    ax[0].set_xlabel("Impact parameter", fontsize=16)
    ax[0].set_ylabel("Radius ratio", fontsize=16)    
    ax[0].set_title("Radius derivative")
    im = ax[1].imshow(np.log10(diff_b), origin="lower", extent=(0, 2, 0, 2), vmin=-15, vmax=0)
    ax[1].set_xlabel("Impact parameter", fontsize=16)
    ax[1].set_ylabel("Radius ratio", fontsize=16)
    ax[1].set_title("Impact parameter derivative")
    cb = pl.colorbar(im)
    cb.ax.set_ylabel("Log difference", fontsize=16);
plot(2)
plot(3)
plot(4)
plot(5)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18996908278941943, 'fraction_numerical': 0.03847475094469255, 'mean_word_length': 3.022099447513812, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '196372', 'license': 'MIT', 'n_tokens_mistral': 1336, 'n_tokens_neox': 1255, 'n_words': 330, 'path': 'proofs/dpdr.ipynb', 'repo_name': 'langfzac/Limbdark.jl', 'size': '2911'}