opticspy / lightpipes

LightPipes for Python, "Pure Python version"
https://opticspy.github.io/lightpipes/
BSD 3-Clause "New" or "Revised" License
227 stars 52 forks source link

LG doughnut mode #71

Open linmanshuo opened 2 years ago

linmanshuo commented 2 years ago

Could you provide the mathematical expression for the LG doughnut mode in LightPipes? Thanks in advance!

ldoyle commented 2 years ago

Dear linmanshuo,

I have never worked with the higher order Gauss beams in LightPipes (or outside), maybe Fred has some more insight. For the moment I can suggest to you to study the source code of the functions to be sure.

The GaussBeam(...,doughnut=True) method will mix (i.e. add the complex fields) two Laguerre Gauss modes of order (n=0,m=1) where one is rotated and phase shifted:

Fout = GaussLaguerre(Fin, w0, p=n, l=m, A = 1.0 ) #first copy
Fout = Interpol( Fout, Fout.siz, Fout.N, 0, 0, 360 / (4 * m), 1) #rotate by 90°/m, e.g. 90° for m=1
Fout = MultPhase(Fout, PI/2 ) #shift phase by pi/2=90°
Fout = BeamMix(GaussLaguerre(Fin, w0, p=n, l=m, A=1 ), Fout) #add complex fields to create spiral phase and doughnut intensity

The mathematical basis used for the GaussLaguerre function is listed in the command reference in addition to the code. https://opticspy.github.io/lightpipes/command-reference.html#LightPipes.GaussLaguerre

I must admit it looks like there is a discrepancy in the prefactors ($\rho$ instead of $\rho/2$) between the command reference and the actual code. This will lead to a wrong amplitude/intensity for a given $A$, but correct shape.

GubioGL commented 1 year ago

I worked with laguer modes in my master's research. I tried using your laguerre command but it has some problems with it.

First, which only has the real part of the" iltheta" phase term, it should have the others as well. According to the phase plots it is only correct for a certain region of the simulation. I looked at all the source code and this is partly due to the scipy functions you call to generate the laguerre-gauss functions

See https://en.wikipedia.org/wiki/Gaussian_beam in Laguerre-Gaussian modes.

If you compare the phase of laguerre-gauss doughnut from the page https://opticspy.github.io/lightpipes/DoughnutModes.html with fig 1 of paper Physical meaning of the radial index of Laguerre-Gauss beams and fig 1 page 3 of Coupling of atmospheric perturbed optical beams with guided modes of propagation you will notice some difference .

I have some code that I used in my master's degree that I can try to implement in lights pipes. I used light pipes to do all my research simulations, but that part had problems. I have other suggestions for implementation for the generation of partially coherent beams (subject of my research).

thanks

FredvanGoor commented 12 months ago

I propose the following modification of the GaussLaguerre command (in core.py):

def GaussLaguerre(Fin, w0, p = 0, l = 0, A = 1.0, ecs = 1 ):
    """
    *Substitutes a Laguerre-Gauss mode (beam waist) in the field.*

        :math:`F_{p,l}(x,y,z=0) = A \\left(\\rho\\right)^{\\frac{|l|}{2} }L^p_l(\\rho)e^{-\\frac{\\rho}{2}}\\cos(l\\theta)`,

        with: :math:`\\rho=\\frac{2(x^2+y^2)}{w_0^2}`

        :math:`\\theta=atan(y/x)`

        if :math:`sincos = 0` replace :math:`cos(l\\theta)` by :math:`exp(-il\\theta)`

        if :math:`sincos = 2` replace :math:`cos(l\\theta)` by :math:`sin(l\\theta)`

    :param Fin: input field
    :type Fin: Field
    :param w0: Gaussian spot size parameter in the beam waist (1/e amplitude point)
    :type w0: int, float
    :param p: mode index (default = 0.0)
    :param l: mode index (default = 0.0)
    :type p: int, float
    :type l: int, float
    :param A: amplitude (default = 1.0)
    :type A: int, float
    :param ecs: 0 = exp, 1 = cos, 2 = sin (default = 1)
    :type ecs: int, float
    :return: output field (N x N square array of complex numbers).
    :rtype: `LightPipes.field.Field`
    :Example:

    >>> F = GaussLaguerre(F, 3*mm) # Fundamental Gauss mode, LG0,0 with a beam radius of 3 mm
    >>> F = GaussLaguerre(F, 3*mm, p=3) # Idem, LG3,0
    >>> F = GaussLaguerre(F, 3*mm, p=3, l=1, A=2.0) # Idem, LG3,1, amplitude 2.0
    >>> F = GaussLaguerre(F, 3*mm, 3, 1, 2.0) # Idem
    >>> F = GaussLaguerre(F, 3*mm, p=3, l=2, ecs=0) LG3,1 with exponential phase factor
    >>> F = GaussLaguerre(F, 3*mm, p=3, l=2, ecs=2) LG3,1 with sine phase factor

    .. seealso::

        * :ref:`Examples: Laguerre Gauss modes.<Laguerre Gauss modes.>`

    Reference::

        A. Siegman, "Lasers", p. 642
    """
    # ************* Backward compatibility section ****************
    #The general backward_compatible decorator does not work for this command,
    #because of the positional argument w0.
    #Old style: GaussLaguerre(p, l, A, w0,Fin)
    #New style: GaussLaguerre(Fin, w0, p=0, l=0, A=1.0)
    _using_oldstyle = False
    if not isinstance(Fin, Field):
        #first arg is not a field, either backward compat syntax or
        # complete usage error -> find out if Field is last, else error
        if isinstance(A, Field):
            #found field in last arg
            _using_oldstyle = True #just in case code wants to know this later
            # in function
            Fin, w0, p, l, A = A, l, Fin, w0, p
            #caution: python can swap the values only if written on single
            # line, if split up a temporary assignment is necessary
            # (since a=b, b=a would not work, only temp=a, a=b, b=temp)
            #-> now all the variables contain what is expected in new style
        else:
            raise ValueError('GaussLaguerre: Field is neither first nor '
                             + 'last parameter (backward compatibility check)'
                             + ', please check syntax/usage.')
    # ************* end of Backward compatibility section *********
    if  ecs not in [0, 1, 2]:
        raise ValueError('GaussLaguerre: invalid value for ecs. ecs must be: 0, 1 or 2')
    Fout = Field.copy(Fin)
    R, Phi = Fout.mgrid_polar
    w02=w0*w0
    la=abs(l)
    rho = 2*R*R/w02
    if ecs == 1:
        Fout.field = A * rho**(la/2) * genlaguerre(p,la)(rho) * _np.exp(-rho/2) * _np.cos(l*Phi)
    if ecs == 2 and l == 0:
        Fout.field = A * rho**(la/2) * genlaguerre(p,la)(rho) * _np.exp(-rho/2)
    if ecs == 2 and l != 0:
        Fout.field = A * rho**(la/2) * genlaguerre(p,la)(rho) * _np.exp(-rho/2) * _np.sin(l*Phi)  
        #################################
    # """
    # Corrections by: GubioGL 23-8-2023:
    # I changed only that part, not calculating the R and PHI with the function,
    # because in the part of the phase of the .mgrid_polar function there is a "+ _np.pi",
    # which shouldn't be here.
    # In order not to change that function, it is simpler to just change it here.
    # optional exponential(..), cos( ..)  or sin( .. ) (esc = 0,1,2)
    # Thank you GubioGL!
    # Fred van Goor, 26-8-2023
    # """
    if ecs == 0:
        Y, X    =   Fout.mgrid_cartesian
        R       =   _np.sqrt(X**2+Y**2)
        rho     =   2*R**2/w0**2
        Fout.field = A * rho**(la/2) * genlaguerre(p,la)(rho) * _np.exp(-rho/2) * _np.exp(-1j*l*_np.arctan2(Y, X))
    Fout._IsGauss=False
    return Fout

Now you can choose your phase factor: exp, cosine or sine with ecs (default = 1: cosine). All these are eigenmodes. Please comment. Fred

GubioGL commented 12 months ago

I think it's a good way to solve it.