Find integration limit given the value of the (2D) integral

68 Views Asked by At

(Previously asked on Stackoverflow - now deleted - it was pointed out that it makes more sense to ask my question here.)

I am trying to solve for R_j(I_0) given in Eq. 12 of this paper (see also Eq. 10; both screenshots below) given that we know the LHS; that is, given the LHS, I'd like to find the integration limit of one of the (two) integrals.

Eq. 12 Eq. 10

Currently, on Python, I am iterating through every value of R_j (i.e., j taking on a value from 0 to N) and finding the value of R_j such that the difference between the LHS and the evaluated integral is minimized. A code snippet is below. However, this is prohibitively slow (with 1600 points, this takes >3000 hours!). Is there a more efficient way of calculating the integration limit?

import numpy as np 
import mpmath as mp
from tqdm import tqdm
for i in tqdm(range(len(R_j))): #for every frequency value
    pdf = lambda r,s: 1/(2*S_j*S_j_err[i]*mp.sqrt(2*np.pi)) * s * mp.exp( -(s-S_j[i])**2/(2*S_j_err[i]**2) - r*s/S_j[i] )
    probs = np.array([ mp.quad(pdf,[R_j[j],mp.inf],[0,mp.inf]) for j in tqdm(range(len(R_j))) ])