Calculation of LiDAR Sensitivity

Question

I have the following paper "Collecting Performance of a LiDAR Telescope at Short Distances": http://earsel.org/wp-content/uploads/2016/11/3-3_03_Ohm.pdf

I am supposed to plot the efficiency of the LiDAR, as shown in Fig. 4 in the paper. However, my graphs do not look at all as they do in the paper. I calculate b, B, M and P and with this AL as shown in the paper. Then I calculate Omega.

Then I integrate, as shown in Eq. (7). I integrate of the implicit variable r from 0 to R= beta *z / 2. The python code I wrote is shown below.

Reasons, I think I am doing something wrong:

• The Graphs look different

• At the distance z0, the graphs are zero.They should not be zero there, instead they should have a “kink”

• There is a removable discontinuity when z = z0

• A lot of the graph is undefined, even when the graph is shown, parts of the makeup is undefined.

The graph I get

import math 
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import numpy as np
from scipy.integrate import quad



def R(z,beta):   
        return z * beta /2


def z0 (d,f):            # particular depth
        if (d-f)==0:
            return None
        return (d*f)/(d-f)


def Dcalc(f, d, beta) :        # radius of diaphrang
        return R(z0(d,f), beta) * d / z0(d,f)

def b(z, f):             #  image distance
        if z-f==0:
            return None
        return (z * f ) / (z - f ) 

def B (z, r, f):        # image size
        if z-f==0:
            return None
        return (r * f ) / (z - f ) 



def M(z, r,f,d) :       # center of projection of diaphrgm to the plan of the lens
        if b(z,f) == None or  (b(z,f) - d ) == 0:
            return None
        return (B(z,r,f) * d ) / (b(z,f) - d ) 


def P(z,f,d,beta) :       # radius of the projection on lens plane
        if b(z,f) == None or  (b(z,f) - d ) == 0:
            return None
        return (Dcalc(f, d, beta) * b(z,f) ) / (b(z,f) - d )


def AL(z,r,f,d,L,sign):    # Area
        if ( M(z,r,f,d)==None or  ((M(z,r,f,d) * L) == 0) or ((M(z,r,f,d) *P(z ,f, d, beta)) == 0 )  or P(z ,f, d, beta) ==None ):
            return None  
        inAa = (( M(z,r,f,d)**2 + L**2 - P(z ,f, d, beta)**2) / (2 * M(z,r,f,d) * L )) 
        inAb = (( M(z,r,f,d)**2 - L**2 + P(z ,f, d, beta)**2) / (2 * M(z,r,f,d) * P(z ,f, d, beta))) 


        if not((inAa >-1  and inAa < 1) and (inAb >-1  and inAb < 1)):
            return None    
        teilA = (L**2 * math.acos( inAa) + P(z ,f, d, beta)**2 *        math.acos( inAb))
        if sign >0:
            teilB = 0.5 * (( 4 * L**2 * M(z,r,f,d)**2 + ( M(z,r,f,d)**2 + L**2 -P(z ,f, d, beta)**2)**2) )**(0.5) 
        else:
            teilB = 0.5 * (( 4 * L**2 * M(z,r,f,d)**2 - ( M(z,r,f,d)**2 + L**2 -P(z ,f, d, beta)**2)**2) )**(0.5)
        return teilA - teilB 




def Omega(z,r,f,d,L,sign):
        if AL(z,r,f,d,L,sign) == None or z==0:
            return None
        return AL(z,r,f,d,L,sign)/ z**2 



def IntegrationOfOmega(z,r,f,d,L,sign):
        if Omega(z,r,f,d,L,sign) ==None:
            return 0
        return Omega(z,r,f,d,L,sign)*r



def Sensitivity (f,d,L,sign, beta):
        sens = []
        zValue = []
        ooz = []
        for iii in np.arange( 2, 60, 0.1):
            integrand = lambda r: IntegrationOfOmega(iii,r,f,d,L,sign)
            zValue.append(iii)
            if iii==0:
                    ooz.append(None)
            else:
                    ooz.append(1.82*L**2/(iii**2))
        

            if (R(iii,beta)) == 0 or R(iii,beta) == None:
                    sens.append (None)
            elif  d-(iii*f)==0   or   (iii-f)==0:      
                    sens.append (None)
            else:
                    result, error = quad(integrand, 0, R(iii,beta))
                    sens.append(2 / (R(iii,beta) ** 2) * result)
        return  zValue , sens , ooz    


def Beta(D,d):
        return 2* D/d


f = 1.2         # 1.2 ... focal length
beta = 2.4e-3   # 2.4e-3 ... laser beam divergence
L=0.1           # 0.1 ... lense radius

sign = 1        # 1 ... in the plus minus part, the sign defines wether plus or minus

d1 = 1.2501     # 1.2501  corresponds to z0=30m
d2 = 1.2329     # 1.2329  corresponds to z0=45m
d3 = 1.2245     # 1.2245  corresponds to z0=60m


Solution25m = Sensitivity (f, d1, L, -sign, beta)
Solution23m = Sensitivity (f, d2, L, -sign, beta)
Solution22m = Sensitivity (f, d3, L, -sign, beta)

plt.plot( Solution25m[0], Solution25m[1], label='z0=30m, sign=-' )
plt.plot( Solution23m[0], Solution23m[1], label='z0=30m, sign=-' )
plt.plot( Solution22m[0], Solution22m[1], label='z0=30m, sign=-' ) 



plt.legend()
plt.xlabel("depth z/m")
plt.ylabel("Sensitivity S")
plt.savefig('plot.png')
plt.show()

I am not entirely sure, if this is a coding problem, or a problem with my understanding of the physics.