Double Concave Lens Help

[jiftesam]

Me along with a team created a double concave lens that simulate the refraction of light passing through it. However, I would like to find the volume in between the space of the concave lens via plotly. How should I go about this?

[mjhoptics]

Hello @jiftesam,

A general way to calculate the volume of an (arbitrary) solid is by ray casting. A parallel grid of rays is traced through the 2 surfaces with the lens material set to air. The length of the ray segment between the surfaces is summed over all of the rays in the grid to get the volume.

Here is an example jupyter notebook demonstrating the approach:

Investigate issue #137, volume of lens

%matplotlib inline
# use widget to enable interactive figures.
#%matplotlib widget

isdark = False

from rayoptics.environment import *

def vol_cyl(h, a):
    ''' compute volume of a right cylinder of height h and radius a. '''
    return np.pi*h*a**2

Create a model with a single bi-concave lens element

opm = OpticalModel()


opm.radius_mode = True
sm  = opm['seq_model']
osp = opm['optical_spec']
pm = opm['parax_model']
em = opm['ele_model']
pt = opm['part_tree']
ar = opm['analysis_results']

osp['pupil'] = PupilSpec(osp, key=['object', 'pupil'], value=4.0)
osp['fov'] = FieldSpec(osp, key=['object', 'angle'], value=2.5, flds=[0.,], is_relative=True)
osp['wvls'] = WvlSpec([('F', 0.5), ('d', 1.0), ('C', 0.5)], ref_wl=1)

sm.gaps[0].thi = 1e+11

opm.add_lens(power=-0.1, th=.3, sd=2.,med='N-BK7, Schott', t=3.)
sm.stop_surface = 1

opm.update_model(scr_model=sm)

sm.list_model()

              r            t        medium     mode   zdr      sd
  Obj:     0.000000  1.00000e+11       air             1      2.0000
 Stop:   -10.386858     0.300000     N-BK7             1      2.0000
    2:    10.386858      3.00000       air             1      2.0465
  Img:     0.000000      0.00000                       1      2.6422

pm.first_order_data()

efl                 -10
f                   -10
f'                  -10
ffl                10.1
pp1             0.09841
bfl               -10.1
ppk            -0.09841
pp sep           0.1032
f/#                -2.5
m                 1e-10
red               1e+10
obj_dist          1e+11
obj_ang               1
enp_dist             -0
enp_radius            2
na obj            2e-11
n obj                 1
img_dist          -10.1
img_ht          -0.1746
exp_dist         -13.29
exp_radius        1.981
na img           0.1961
n img                 1
optical invariant      0.03491

layout_plt = plt.figure(FigureClass=InteractiveLayout, opt_model=opm, 
                        do_draw_ray_fans=True, do_draw_edge_rays=False, is_dark=False).plot()

Change the lens material to air, so that no refraction takes place

sm.gaps[1].medium = Air()
opm.update_model(scr_model=sm)

ray_casting_plt = plt.figure(FigureClass=InteractiveLayout, opt_model=opm, 
                        do_draw_ray_fans=True, do_draw_edge_rays=False, is_dark=False).plot()

f0 = osp['fov'].fields[0]
wvl = osp['wvls'].central_wvl

num_rays=21
r2g_pkg = (sampler.create_generator, (sampler.R_2_quasi_random_generator, num_rays**2),
            dict(mapper=sampler.concentric_sample_disk))
ray_list = analyses.RayList(opm, pupil_gen=r2g_pkg, f=f0, wl=wvl)

Plot the sampling points in the pupil

csd_r2_grid_x = []
csd_r2_grid_y = []
pupil_coords = r2g_pkg[0](*r2g_pkg[1], **r2g_pkg[2])
for csd in pupil_coords:
    csd_r2_grid_x.append(csd[0])
    csd_r2_grid_y.append(csd[1])

fig, ax = plt.subplots()
ax.scatter(csd_r2_grid_x, csd_r2_grid_y, s=4)
ax.set_aspect('equal')
plt.show()

Calculate the lens volume by ray casting

num_samples = len(ray_list.ray_list)

r = abs(sm.ifcs[1].profile.r)
a = osp['pupil'].value/2

sum_dst = 0
for ray_pkg in ray_list.ray_list:
    p_x, p_y, ray_pkg = ray_pkg
    ray, op, wv = ray_pkg
    dst = RaySeg(*ray[1]).dst
    sum_dst += dst
A = np.pi*a**2
dA = A/num_samples
V_cast = sum_dst * dA

calculate exact volume for a lens

assumes equi-radii surfaces

h = r - np.sqrt(r**2 - a**2)
th = sm.gaps[1].thi + 2*h

v_cap = np.pi*h*(3*a**2 + h**2)/6
v_th = np.pi*th*a**2
V_exact = v_th - 2*v_cap

print(f"{V_exact=:8.4f}   {V_cast=:8.4f}   diff={V_exact-V_cast:8.4g}")

V_exact=  6.2047   V_cast=  6.2052   diff=-0.0004978

Hope this helps.
Mike Hayford