Garchupiter / Kraken-Optical-Simulator

Python - Exact ray tracing library
GNU General Public License v3.0
65 stars 14 forks source link

correct raytracing of diffraction grating? #9

Open MSLDgit opened 1 year ago

MSLDgit commented 1 year ago

Hi Garchupiter. This is a great project, thank you. I have a question regarding the diffraction gratings: I think, I noticed you use the raytrace equations from Ludwig 1973 (which are basically the equations from Spencer&Murty 1962). When I myself tried to implement these equations into the FreeCAD opticsworkbench I came across some irregularity and I couldnt find a solution yet. I see the same irregularrity with KrakenOS. Eg. I modified your grating example script to demonstrate this. A ray with wavelength 0.55µm (first orders), incident angle of +45 degree relative to the grating normal on a 1000 lpm grating should according to the grating equation ((n*lambda/d)=sin(alpha)+sin(beta), simplified for ray-plane perpendecular to grooves) be diffracted with an angle of -9.04° relative to the grating normal (opposite side of incident ray, relative to the grating normal). However using the Ludwig 1973 equations this geometry leads to an diffracted ray with an angle of 18.55° relative to the grating normal (I encounter this in my implementation and also in the KrakenOS implementation, see modified grating example script below). I am stuck here and wonder what your thoughts about this are, also because it might introduce error. Maybe I am overlooking something here.. Thank you! Miles

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""Examp Diffraction Grating Reflection"""

import numpy as np
import pkg_resources
""" Looking for if KrakenOS is installed, if not, it assumes that
an folder downloaded from github is run"""

required = {'KrakenOS'}
installed = {pkg.key for pkg in pkg_resources.working_set}
missing = required - installed

if missing:
    print("Not installed")
    import sys
    sys.path.append("../..")

import KrakenOS as Kos
### example modified for simple 45° incident ray on grating geometry
# _________________________________________#

P_Obj = Kos.surf()
P_Obj.Rc = 0.0
P_Obj.Thickness = 10
P_Obj.Glass = "AIR"
P_Obj.Diameter = 30.0

# _________________________________________#

# L1a = Kos.surf()
# L1a.Rc = 5.513435044607768E+001
# L1a.Thickness = 6.0
# L1a.Glass = "BK7"
# L1a.Diameter = 30.0

# # _________________________________________#

# L1b = Kos.surf()
# L1b.Rc = -4.408716526030626E+001
# L1b.Thickness = 3.0
# L1b.Glass = "F2"
# L1b.Diameter = 30

# # _________________________________________#

# L1c = Kos.surf()
# L1c.Rc = -2.246906271406796E+002
# L1c.Thickness = 9.737871661422000E+001 - 50.0
# L1c.Glass = "AIR"
# L1c.Diameter = 30

# _________________________________________#

Dif_Obj = Kos.surf()
Dif_Obj.Rc = 0.0
Dif_Obj.Thickness = -50
Dif_Obj.Glass = "MIRROR"
Dif_Obj.Diameter = 30.0
Dif_Obj.Grating_D = 1
Dif_Obj.Diff_Ord = -1
Dif_Obj.Grating_Angle = 90
# Dif_Obj.Surface_type = 1

# _________________________________________#

P_Ima = Kos.surf()
P_Ima.Rc = 0.0
P_Ima.Name = "Plano imagen"
P_Ima.Thickness = 0.0
P_Ima.Glass = "AIR"
P_Ima.Diameter = 300.0
P_Ima.Drawing = 0

# _________________________________________#

A = [P_Obj, Dif_Obj, P_Ima]
configuracion_1 = Kos.Setup()

# _________________________________________#

Doblete = Kos.system(A, configuracion_1)
Rayos = Kos.raykeeper(Doblete)

# _________________________________________#

tam = 1
rad = 10.0
tsis = len(A) - 1
for i in range(-tam, tam + 1):
    for j in range(-tam, tam + 1):

        x_0 = (i / tam) * rad
        y_0 = (j / tam) * rad
        r = np.sqrt((x_0 * x_0) + (y_0 * y_0))
        if r < rad:
            tet = 0.0
            pSource_0 = [x_0, y_0, 0.0]
            #dCos = [0.0, np.sin(np.deg2rad(tet)), np.cos(np.deg2rad(tet))]
            dCos = [1, 0, np.cos(np.deg2rad(tet))]
            # W = 0.4
            # Doblete.Trace(pSource_0, dCos, W)
            # Rayos.push()
            W = 0.55
            Doblete.Trace(pSource_0, dCos, W)
            Rayos.push()
            # W = 0.6
            # Doblete.Trace(pSource_0, dCos, W)
            # Rayos.push()

# ______________________________________#
print(Doblete.XYZ)
print(Doblete.LMN)
### calculate angle between grating normal and diffracted ray
print(np.rad2deg(np.arccos((np.dot((np.array([0.,0.,-1.])),Doblete.LMN[1]))/(np.linalg.norm((np.array([0.,0.,-1.])))*np.linalg.norm(Doblete.LMN[1])))))
Kos.display2d(Doblete, Rayos, 1)
### calculate angle according to typical grating equation 
### sin(beta) = (m*W/d)-sin(alpha), where beta is the angle of the diffracted ray and alpha is the angle of the incident ray
### both relative to the grating normal
print(np.rad2deg(np.arcsin((((1*0.55)/(1000/1000))-np.sin(np.deg2rad(45))))))
Garchupiter commented 1 year ago

Hi Miles, I add a file with my modification and some comments to achieve the situation you want, you will see that it is a problem in how you are considering the ray and the rotation of the diffraction grating on its own axis as 90 degrees

Line 63 changed from 90 degrees to 0 degrees

the three components (direction cosines) must must be normalized, you have two directions in these case, it can be seen in the next two images from different axis (x and y) changing line 117 (value 1 to 0).

image

image

The line 100 must be uncommented and the line 101 commented. The angle “tet” in line 98 must be 45 degrees (For the ray angle) and the diffraction order in line 62 must be 1, using these modifications now the system is the scenario that you propose. image

 

Please don't hesitate to contact me if you have any other questions. Regards, Joel H.


#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""Examp Diffraction Grating Reflection"""

import numpy as np
import pkg_resources
""" Looking for if KrakenOS is installed, if not, it assumes that
an folder downloaded from github is run"""

required = {'KrakenOS'}
installed = {pkg.key for pkg in pkg_resources.working_set}
missing = required - installed

if missing:
    print("Not installed")
    import sys
    sys.path.append("../..")

import KrakenOS as Kos
### example modified for simple 45° incident ray on grating geometry
# _________________________________________#

P_Obj = Kos.surf()
P_Obj.Rc = 0.0
P_Obj.Thickness = 10
P_Obj.Glass = "AIR"
P_Obj.Diameter = 30.0

# _________________________________________#

# L1a = Kos.surf()
# L1a.Rc = 5.513435044607768E+001
# L1a.Thickness = 6.0
# L1a.Glass = "BK7"
# L1a.Diameter = 30.0

# # _________________________________________#

# L1b = Kos.surf()
# L1b.Rc = -4.408716526030626E+001
# L1b.Thickness = 3.0
# L1b.Glass = "F2"
# L1b.Diameter = 30

# # _________________________________________#

# L1c = Kos.surf()
# L1c.Rc = -2.246906271406796E+002
# L1c.Thickness = 9.737871661422000E+001 - 50.0
# L1c.Glass = "AIR"
# L1c.Diameter = 30

# _________________________________________#

Dif_Obj = Kos.surf()
Dif_Obj.Rc = 0.0
Dif_Obj.Thickness = -50
Dif_Obj.Glass = "MIRROR"
Dif_Obj.Diameter = 30.0
Dif_Obj.Grating_D = 1
Dif_Obj.Diff_Ord = 1
Dif_Obj.Grating_Angle = 0
# Dif_Obj.Surface_type = 1

# _________________________________________#

P_Ima = Kos.surf()
P_Ima.Rc = 0.0
P_Ima.Name = "Plano imagen"
P_Ima.Thickness = 0.0
P_Ima.Glass = "AIR"
P_Ima.Diameter = 300.0
P_Ima.Drawing = 0

# _________________________________________#

A = [P_Obj, Dif_Obj, P_Ima]
configuracion_1 = Kos.Setup()

# _________________________________________#

Doblete = Kos.system(A, configuracion_1)
Rayos = Kos.raykeeper(Doblete)

# _________________________________________#

tam = 1
rad = 10.0
tsis = len(A) - 1
for i in range(-tam, tam + 1):
    for j in range(-tam, tam + 1):

        x_0 = (i / tam) * rad
        y_0 = (j / tam) * rad
        r = np.sqrt((x_0 * x_0) + (y_0 * y_0))
        if r < rad:
            tet = 45
            pSource_0 = [x_0, y_0, 0.0]
            dCos = [0.0, np.sin(np.deg2rad(tet)), np.cos(np.deg2rad(tet))]
            # dCos = [1, 0, np.cos(np.deg2rad(tet))]
            # W = 0.4
            # Doblete.Trace(pSource_0, dCos, W)
            # Rayos.push()
            W = 0.55
            Doblete.Trace(pSource_0, dCos, W)
            Rayos.push()
            # W = 0.6
            # Doblete.Trace(pSource_0, dCos, W)
            # Rayos.push()

# ______________________________________#
print(Doblete.XYZ)
print(Doblete.LMN)
### calculate angle between grating normal and diffracted ray
print(np.rad2deg(np.arccos((np.dot((np.array([0.,0.,-1.])),Doblete.LMN[1]))/(np.linalg.norm((np.array([0.,0.,-1.])))*np.linalg.norm(Doblete.LMN[1])))))
Kos.display2d(Doblete, Rayos, 0)
### calculate angle according to typical grating equation
### sin(beta) = (m*W/d)-sin(alpha), where beta is the angle of the diffracted ray and alpha is the angle of the incident ray
### both relative to the grating normal
print(np.rad2deg(np.arcsin((((1*0.55)/(1000/1000))-np.sin(np.deg2rad(45))))))

El 20 ene 2023, a las 8:14, MSLDgit @.***> escribió:

Hi Garchupiter. This is a great project, thank you. I have a question regarding the diffraction gratings: I think, I noticed you use the raytrace equations from Ludwig 1973 (which are basically the equations from Spencer&Murty 1962). When I myself tried to implement these equations into the FreeCAD opticsworkbench I came across some irregularity and I couldnt find a solution yet. I see the same irregularrity with KrakenOS. Eg. I modified your grating example script to demonstrate this. A ray with wavelength 0.55µm (first orders), incident angle of +45 degree relative to the grating normal on a 1000 lpm grating should according to the grating equation ((n*lambda/d)=sin(alpha)+sin(beta), simplified for ray-plane perpendecular to grooves) be diffracted with an angle of -9.04° relative to the grating normal (opposite side of incident ray, relative to the grating normal). However using the Ludwig 1973 equations this geometry leads to an diffracted ray with an angle of 18.55° relative to the grating normal (I encounter this in my implementation and also in the KrakenOS implementation, see modified grating example script below). I am stuck here and wonder what your thoughts about this are, also because it might introduce error. Maybe I am overlooking something here.. Thank you! Miles

!/usr/bin/env python3

-- coding: utf-8 --

"""Examp Diffraction Grating Reflection"""

import numpy as np import pkg_resources """ Looking for if KrakenOS is installed, if not, it assumes that an folder downloaded from github is run"""

required = {'KrakenOS'} installed = {pkg.key for pkg in pkg_resources.working_set} missing = required - installed

if missing: print("Not installed") import sys sys.path.append("../..")

import KrakenOS as Kos

example modified for simple 45° incident ray on grating geometry

_____

P_Obj = Kos.surf() P_Obj.Rc = 0.0 P_Obj.Thickness = 10 P_Obj.Glass = "AIR" P_Obj.Diameter = 30.0

_____

L1a = Kos.surf()

L1a.Rc = 5.513435044607768E+001

L1a.Thickness = 6.0

L1a.Glass = "BK7"

L1a.Diameter = 30.0

_____

L1b = Kos.surf()

L1b.Rc = -4.408716526030626E+001

L1b.Thickness = 3.0

L1b.Glass = "F2"

L1b.Diameter = 30

_____

L1c = Kos.surf()

L1c.Rc = -2.246906271406796E+002

L1c.Thickness = 9.737871661422000E+001 - 50.0

L1c.Glass = "AIR"

L1c.Diameter = 30

_____

Dif_Obj = Kos.surf() Dif_Obj.Rc = 0.0 Dif_Obj.Thickness = -50 Dif_Obj.Glass = "MIRROR" Dif_Obj.Diameter = 30.0 Dif_Obj.Grating_D = 1 Dif_Obj.Diff_Ord = -1 Dif_Obj.Grating_Angle = 90

Dif_Obj.Surface_type = 1

_____

P_Ima = Kos.surf() P_Ima.Rc = 0.0 P_Ima.Name = "Plano imagen" P_Ima.Thickness = 0.0 P_Ima.Glass = "AIR" P_Ima.Diameter = 300.0 P_Ima.Drawing = 0

_____

A = [P_Obj, Dif_Obj, P_Ima] configuracion_1 = Kos.Setup()

_____

Doblete = Kos.system(A, configuracion_1) Rayos = Kos.raykeeper(Doblete)

_____

tam = 1 rad = 10.0 tsis = len(A) - 1 for i in range(-tam, tam + 1): for j in range(-tam, tam + 1):

    x_0 = (i / tam) * rad
    y_0 = (j / tam) * rad
    r = np.sqrt((x_0 * x_0) + (y_0 * y_0))
    if r < rad:
        tet = 0.0
        pSource_0 = [x_0, y_0, 0.0]
        #dCos = [0.0, np.sin(np.deg2rad(tet)), np.cos(np.deg2rad(tet))]
        dCos = [1, 0, np.cos(np.deg2rad(tet))]
        # W = 0.4
        # Doblete.Trace(pSource_0, dCos, W)
        # Rayos.push()
        W = 0.55
        Doblete.Trace(pSource_0, dCos, W)
        Rayos.push()
        # W = 0.6
        # Doblete.Trace(pSource_0, dCos, W)
        # Rayos.push()

__

print(Doblete.XYZ) print(Doblete.LMN)

calculate angle between grating normal and diffracted ray

print(np.rad2deg(np.arccos((np.dot((np.array([0.,0.,-1.])),Doblete.LMN[1]))/(np.linalg.norm((np.array([0.,0.,-1.])))*np.linalg.norm(Doblete.LMN[1]))))) Kos.display2d(Doblete, Rayos, 1)

calculate angle according to typical grating equation

sin(beta) = (m*W/d)-sin(alpha), where beta is the angle of the diffracted ray and alpha is the angle of the incident ray

both relative to the grating normal

print(np.rad2deg(np.arcsin((((1*0.55)/(1000/1000))-np.sin(np.deg2rad(45)))))) — Reply to this email directly, view it on GitHub https://github.com/Garchupiter/Kraken-Optical-Simulator/issues/9, or unsubscribe https://github.com/notifications/unsubscribe-auth/AED73RRIW7KKLBUUOCYESQTWTK2YFANCNFSM6AAAAAAUBWJYSI. You are receiving this because you are subscribed to this thread.