Simulating collecting efficency of a LiDAR

  • A
  • Thread starter NielsW
  • Start date
  • #1
NielsW
2
0
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 calculate 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
• 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

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

f = 1.2         # focal length
beta = 2.4e-3   # laser beam divergence
L=0.1           # lense radius
d= 1.25         # distance of diaphrang

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

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) :        # 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,D,f,d) :       # radius of the projection on lens plane
    if b(z,f) == None or  (b(z,f) - d ) == 0:
        return None
    return (D * b(z,f) ) / (b(z,f) - d )

def AL(z,r,f,d,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,D,f,d)) == 0 )  or P(z,D,f,d) ==None ):
        return None  
    inAa = 1.1*(( M(z,r,f,d)**2 + L**2 - P(z,D,f,d)**2) / (2 * M(z,r,f,d) * L )) 
    inAb = 1.1*(( M(z,r,f,d)**2 - L**2 + P(z,D,f,d)**2) / (2 * M(z,r,f,d) * P(z,D,f,d))) 
    
    if not((inAa >-1  and inAa < 1) and (inAb >-1  and inAb < 1)):
        return None    
    teilA = (L**2 * math.acos( inAa) + P(z,D,f,d)**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,D,f,d)**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,D,f,d)**2)**2) )**(0.5)
    return teilA - teilB 

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

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

def Sensitivity (f,d,D,L,sign, beta):
    sens = []
    zValue = []
    ooz = []
    for iii in np.arange( 2, 60, 0.5):
        integrand = lambda r: IntegrationOfOmega(iii,r,f,d,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)/(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    

D = Dcalc(f, d)

#'''
Solution25p = Sensitivity (f,1.2501, D,L, sign, beta)
Solution23p = Sensitivity (f,1.2329, D,L, sign, beta)
Solution22p = Sensitivity (f,1.2245, D,L, sign, beta)

#'''

Solution25m = Sensitivity (f,1.2501, D,L,-sign, beta)
Solution23m = Sensitivity (f,1.2329, D,L,-sign, beta)
Solution22m = Sensitivity (f,1.2245, D,L,-sign, beta)

#'''

plt.plot( Solution25p[0], Solution25p[1],label='z0=30m, sign=+' )
plt.plot( Solution23p[0], Solution23p[1],label='z0=45m, sign=+' )
plt.plot( Solution22p[0], Solution22p[1],label='z0=60m, sign=+')

#'''

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")
plt.show()
 
Last edited by a moderator:
Science news on Phys.org
  • #2
Welcome to PF.

NielsW said:
Reasons, I think I am doing something wrong:
• The Graphs look different
• 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
Can you upload your graphs for comparison? Use the "Attach files" link below the Edit window. Thank you.
 
  • #3
NielsW said:
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 calculate 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
• 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

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

f = 1.2         # focal length
beta = 2.4e-3   # laser beam divergence
L=0.1           # lense radius
d= 1.25         # distance of diaphrang

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

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) :        # 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,D,f,d) :       # radius of the projection on lens plane
    if b(z,f) == None or  (b(z,f) - d ) == 0:
        return None
    return (D * b(z,f) ) / (b(z,f) - d )

def AL(z,r,f,d,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,D,f,d)) == 0 )  or P(z,D,f,d) ==None ):
        return None 
    inAa = 1.1*(( M(z,r,f,d)**2 + L**2 - P(z,D,f,d)**2) / (2 * M(z,r,f,d) * L ))
    inAb = 1.1*(( M(z,r,f,d)**2 - L**2 + P(z,D,f,d)**2) / (2 * M(z,r,f,d) * P(z,D,f,d)))
   
    if not((inAa >-1  and inAa < 1) and (inAb >-1  and inAb < 1)):
        return None   
    teilA = (L**2 * math.acos( inAa) + P(z,D,f,d)**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,D,f,d)**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,D,f,d)**2)**2) )**(0.5)
    return teilA - teilB

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

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

def Sensitivity (f,d,D,L,sign, beta):
    sens = []
    zValue = []
    ooz = []
    for iii in np.arange( 2, 60, 0.5):
        integrand = lambda r: IntegrationOfOmega(iii,r,f,d,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)/(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   

D = Dcalc(f, d)

#'''
Solution25p = Sensitivity (f,1.2501, D,L, sign, beta)
Solution23p = Sensitivity (f,1.2329, D,L, sign, beta)
Solution22p = Sensitivity (f,1.2245, D,L, sign, beta)

#'''

Solution25m = Sensitivity (f,1.2501, D,L,-sign, beta)
Solution23m = Sensitivity (f,1.2329, D,L,-sign, beta)
Solution22m = Sensitivity (f,1.2245, D,L,-sign, beta)

#'''

plt.plot( Solution25p[0], Solution25p[1],label='z0=30m, sign=+' )
plt.plot( Solution23p[0], Solution23p[1],label='z0=45m, sign=+' )
plt.plot( Solution22p[0], Solution22p[1],label='z0=60m, sign=+')

#'''

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")
plt.show()
plot.png
 

Similar threads

  • Engineering and Comp Sci Homework Help
Replies
6
Views
962
  • Engineering and Comp Sci Homework Help
Replies
3
Views
1K
Replies
3
Views
187
  • Introductory Physics Homework Help
Replies
25
Views
403
  • Advanced Physics Homework Help
Replies
8
Views
1K
  • Introductory Physics Homework Help
Replies
6
Views
280
  • Advanced Physics Homework Help
Replies
1
Views
203
Replies
0
Views
93
  • Introductory Physics Homework Help
Replies
3
Views
369
Replies
1
Views
715
Back
Top