bryancole / raypier_optics

A raytracing toolkit for optical design
Other
57 stars 8 forks source link

Modelling an aspheric lens: AFL12-15 #10

Open Jonas231 opened 2 years ago

Jonas231 commented 2 years ago

Hello Bryan,

I have this code (which basically stems from here: https://raypier-optics.readthedocs.io/en/latest/introduction.html#basic-usage). The question is: How can I define some aspherical surface?

I attached a test file. Do you see what is wrong?

Best regards, Jonas

import numpy as np

coeffs = [0.683740494, 0.420323613, 0.58502748, 0.00460352869, 0.0133968856, 64.4932732] Lmin = 0.18 Lmax = 2.3 # micrometer from raypier.api import OpticalMaterialFromFormula m1 = OpticalMaterialFromFormula(formula_id = 2, coefs = coeffs, wavelen_min = Lmin, wavelen_max = Lmax)

from raypier.api import GeneralLens, AxiconFace, PlanarFace, OpticalMaterial, CircleShape,\ RayTraceModel, CollimatedGaussletSource, EFieldPlane, GaussletCapturePlane, IntensitySurface, ParallelRaySource, SphericalFace, AsphericFace

Shape = CircleShape(radius=6.25) # APE-Y [mm]
cuy_L = [0.0, 0.13537295248409, 0.0, 0.0] #curvatures [1/mm] thi_L = [1e+20, 4.0, 12.32032337, 0.0] # thickness [mm]

aspheric coefficients

asp_L = [np.array([], dtype='<U118'), np.array([' ASP -0.6700000000 0.3440000000E-04 0.1430000000E-06 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000'], dtype='<U181'), np.array([], dtype='<U118'), np.array([], dtype='<U118')]

sut_L = ['S', 'A', 'S', 'S'] # surface type

semi-apertures:

ape_L = [" 1 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 1 0 0 1 ''", " 1 5.625000000 5.625000000 0.000000000 0.000000000 0.000000000 1 0 0 1 ''", " 1 6.250000000 6.250000000 0.000000000 0.000000000 0.000000000 1 0 0 1 ''", " 1 0.1676174929E-04 0.1676174929E-04 0.000000000 0.000000000 0.000000000 1 0 0 1 ''"]

Nsurfaces = 4

if 1:

Build a couple of lenses

shape_L = []
radius_L = []
f_L = []
for i in range(1,Nsurfaces):
    radius_L.append(float(ape_L[i][5:5+10]))
    shape_L.append(CircleShape(radius=radius_L[i-1]))
    print("dia: ", float(ape_L[i][5:5+10]))
    if sut_L[i] == "S":
        print("i: ",i , " sphere")

        f_L.append(SphericalFace(curvature=cuy_L[i], z_height=thi_L[i]))
    if sut_L[i] == "A":
        print("i: ",i , " asphere")
        asp_coeffs = asp_L[i][0].split()
        asphere = AsphericFace(z_height=thi_L[i],
                              curvature=cuy_L[i],            
                              conic_const= float(asp_coeffs[1]),
                              A4 = float(asp_coeffs[2]),  # A
                              A6 = float(asp_coeffs[3]), #I'm not 100% sure if the signs of the polynomial terms are right.
                              A8 = float(asp_coeffs[4]), # C
                              A10 = float(asp_coeffs[5]), # D
                              A12 = float(asp_coeffs[6]), # E
                              A14 = float(asp_coeffs[7]), # F
                              A16 = float(asp_coeffs[8])#, # G
                              #A18 = asp_coeffs[9]  # H
                              )
        f_L.append(asphere)

lens = GeneralLens(centre=(0,0,0), direction=(0,0,1), shape=Shape, surfaces=[f_L[0], f_L[1]], materials=[m1])

Add a source

src = ParallelRaySource(origin=(0,0,-5.0), direction=(0,0,1), rings=5, show_normals=True, display="wires", opacity=0.1)

model = RayTraceModel(optics=[lens], sources=[src])

model.ipython_view(800,600)

model.configure_traits()

Jonas231 commented 2 years ago

raypier_simple.txt

bryancole commented 2 years ago

Hi Jonas, I think the problem is the definition of "curvature". For the General Lens surface types, I have defined "curvature" as a shorthand for "radius of curvature", rather than the arguably more logical 1./radius_of_curvature.This also means that it is not possible to specify a spherical surface with infinite radius of curvature (i.e. curvature of zero, meaning a plane). Hence, in your list of curvatures, you need to give real radius's or if youactually want a plane, use a PlanarFace class. I would agree the ambiguity between curvature vs. radius-of-curvature is annoying and confusing. It would have been nice to specify curvature as a quantity than can vary smoothly through zero. A second problem is that one of your surfaces has a thickness of 1e20. This is unlikely to give good results, although at present, you are not including it in the for-loop (because you skip index 0). Bryan On Tue, 2022-06-14 at 07:11 -0700, Jonas231 wrote:

Hello Bryan, I have this code (which basically stems from here: https://raypier-optics.readthedocs.io/en/latest/introduction.html#basic-usage ).

The question is: How can I define some aspherical surface? I attached a test file. Do you see what is wrong? Best regards,

Jonas import numpy as np coeffs = [0.683740494, 0.420323613, 0.58502748, 0.00460352869, 0.0133968856, 64.4932732]

Lmin = 0.18

Lmax = 2.3 # micrometer

from raypier.api import OpticalMaterialFromFormula

m1 = OpticalMaterialFromFormula(formula_id = 2, coefs = coeffs, wavelen_min = Lmin, wavelen_max = Lmax) from raypier.api import GeneralLens, AxiconFace, PlanarFace, OpticalMaterial, CircleShape,

RayTraceModel, CollimatedGaussletSource, EFieldPlane, GaussletCapturePlane, IntensitySurface, ParallelRaySource, SphericalFace, AsphericFace Shape = CircleShape(radius=6.25) # APE-Y [mm]

cuy_L = [0.0, 0.13537295248409, 0.0, 0.0] #curvatures [1/mm]

thi_L = [1e+20, 4.0, 12.32032337, 0.0] # thickness [mm] aspheric coefficients asp_L = [np.array([], dtype='<U118'),

np.array([' ASP -0.6700000000 0.3440000000E-04 0.1430000000E- 06 0.000000000 0.000000000 0.000000000 0.00000000 0 0.000000000 0.000000000 0.000000000'],

dtype='<U181'),

np.array([], dtype='<U118'),

np.array([], dtype='<U118')] sut_L = ['S', 'A', 'S', 'S'] # surface type semi-apertures: ape_L = [" 1 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 1 0 0 1 ''",

" 1 5.625000000 5.625000000 0.000000000 0.000000000 0.000000000 1 0 0 1 ''",

" 1 6.250000000 6.250000000 0.000000000 0.000000000 0.000000000 1 0 0 1 ''",

" 1 0.1676174929E-04 0.1676174929E- 04 0.000000000 0.000000000 0.000000000 1 0 0 1 ''"] Nsurfaces = 4 if 1:

Build a couple of lenses

shape_L = []

radius_L = []

f_L = []

for i in range(1,Nsurfaces):

radius_L.append(float(ape_L[i][5:5+10]))

shape_L.append(CircleShape(radius=radius_L[i-1]))

print("dia: ", float(ape_L[i][5:5+10]))

if sut_L[i] == "S":

print("i: ",i , " sphere") f_L.append(SphericalFace(curvature=cuy_L[i], z_height=thi_L[i])) if sut_L[i] == "A": print("i: ",i , " asphere") asp_coeffs = asp_L[i][0].split() asphere = AsphericFace(z_height=thi_L[i], curvatur e=cuy_L[i], conic_const= float(asp_coeffs[1]), A4 = float(asp_coeffs[2]), # A A6 = float(asp_coeffs[3]), #I'm not 100% sure if the signs of the polynomial terms are right. A8 = float(asp_coeffs[4]), # C A10 = float(asp_coeffs[5]), # D A12 = float(asp_coeffs[6]), # E A14 = float(asp_coeffs[7]), # F A16 = float(asp_coeffs[8])#, # G #A18 = asp_coeffs[9] # H ) f_L.append(asphere) lens = GeneralLens(centre=(0,0,0),

direction=(0,0,1),

shape=Shape,

surfaces=[f_L[0], f_L[1]],

materials=[m1]) Add a source src = ParallelRaySource(origin=(0,0,-5.0),

direction=(0,0,1),

rings=5,

show_normals=True,

display="wires",

opacity=0.1) model = RayTraceModel(optics=[lens],

sources=[src])

model.ipython_view(800,600)

model.configure_traits()

— Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you are subscribed to this thread.Message ID: @.***>

Jonas231 commented 2 years ago

Hi Bryan, OpTaliX_raypier_simple.txt thank you. I did not see that. Okay I converted to radius of curvature and replaced inf by 1e10. I tried the Spherical Face class with radius of curvature of 1e10. This should also work instead of PlanarFace class, right? 1e20 did not work, though. One could make spherical class fall back to Planar if radius of curvature exceeds a certain value.

cuy_L = [0.0, -0.13537295248409, 0.0, 0.0] # curvatures [1/mm] rad_L = 1/np.array(cuy_L) # radius of curvature [mm] rad_L[rad_L == np.inf] = 1e10 thi_L = [1e+20, 0, 4, 12.32032337]

Now the layout plot in raypier looks like this: grafik

Compare to ZEMAX:

grafik

A sphere is produced on the left side. I am not sure why. How can I limit the number of rays shown only to the transmitted ones (like in ZEMAX)? Jonas

bryancole commented 2 years ago

Hi Jonas, The problem with using a SphericalFace with very very large radius is that the intersection calculation subtracts the incident rays start/end points with the sphere centre point. Substracting two very large values creates a large error in the floating-point maths. Really, I should re-write the SphericalFace class to use a different intersection formula. I think putting in a fall-back to substitute in a PlanarFace if the radius gets too big is the best plan (until I can re-write the underlying SphericalFace). To avoid having so many weak reflections, you can set the "reflection_threshold" of the material attribute of the len faces. E.g. for face in lens.faces.faces: face.material.reflection_threshold=1.1 #a value >1.0 ensures no reflections are produced. Unfortunately, the reflection_threshold and transmission_threshold attributes only exist on the underlying cmaterials.InterfaceMaterial subclasses, and are not exposed through the OpticalMaterial object, so this is a little inconvenient using the GeneralOptic framework. I've just added these attributes to the GeneralLens class and pushed to github. As second problem with your model is that the refractive index values that your custom glass is evaluating seem rather high (I found a value of >3.5 at 0.78um wavelength). Is that what you expect?The very high refractive index results in a lot of total-internal reflections. If you swap in as standard glass e.g. m1 = OpticalMaterial(glass_name="N-BK7") then things look better. Maybe the glass dispersion coefficients ordering is wrong for your glass-type. Bryan On Wed, 2022-06-15 at 02:30 -0700, Jonas231 wrote:

Hi Bryan,

OpTaliX_raypier_simple.txt

thank you. I did not see that. Okay I converted to radius of curvature and replaced inf by 1e10.

I tried the Spherical Face class with radius of curvature of 1e10. This should also work instead of PlanarFace class, right?

1e20 did not work, though. One could make spherical class fall back to Planar if radius of curvature exceeds a certain value. cuy_L = [0.0, -0.13537295248409, 0.0, 0.0] # curvatures [1/mm]

rad_L = 1/np.array(cuy_L) # radius of curvature [mm]

rad_L[rad_L == np.inf] = 1e10

thi_L = [1e+20, 0, 4, 12.32032337] Now the layout plot in raypier looks like this:

Compare to ZEMAX:

A sphere is produced on the left side. I am not sure why.

How can I limit the number of rays shown only to the transmitted ones (like in ZEMAX)?

Jonas

— Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you commented.Message ID: < @.***>

bryancole commented 2 years ago

Yes, just checked. Your dispersion coefficients are in the wrong order (and missing the first zero). Check the formula definitions for your glass data. I think you want: coeffs = [0.0, 0.683740494, 0.00460352869, 0.420323613, 0.0133968856, 0.58502748, 64.4932732] It's basically fused silica. Bryan On Wed, 2022-06-15 at 02:30 -0700, Jonas231 wrote:

Hi Bryan,

OpTaliX_raypier_simple.txt

thank you. I did not see that. Okay I converted to radius of curvature and replaced inf by 1e10.

I tried the Spherical Face class with radius of curvature of 1e10. This should also work instead of PlanarFace class, right?

1e20 did not work, though. One could make spherical class fall back to Planar if radius of curvature exceeds a certain value. cuy_L = [0.0, -0.13537295248409, 0.0, 0.0] # curvatures [1/mm]

rad_L = 1/np.array(cuy_L) # radius of curvature [mm]

rad_L[rad_L == np.inf] = 1e10

thi_L = [1e+20, 0, 4, 12.32032337] Now the layout plot in raypier looks like this:

Compare to ZEMAX:

A sphere is produced on the left side. I am not sure why.

How can I limit the number of rays shown only to the transmitted ones (like in ZEMAX)?

Jonas

— Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you commented.Message ID: < @.***>