pdtype.py 5.04 KB
Newer Older
Andreas Freise's avatar
Andreas Freise committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals

import pykat.exceptions as pkex
import numpy as np
import math
import copy
import warnings
import cmath
from math import factorial
from scipy.integrate import quad
from scipy.special import hermite
from pykat.SIfloat import SIfloat
16
17
18
19
#import scipy.special
from pykat.optics.gaussian_beams import HG2LG
from pykat.math.jacobi import jacobi
from pykat.math.laguerre import laguerre
Andreas Freise's avatar
Andreas Freise committed
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76


def finesse_split_photodiode(maxtem, xy='x'):
    """Prints beat coefficients for HG modes on a split photo detector
    in the format for insertion into the kat.ini file of Finesse.
    Example use:
    finesse_split_photodiode(5, xy='x')
    The default kat.ini numbers have been generated with:
    finesse_split_photodiode(40, xy='x')
    finesse_split_photodiode(40, xy='y')
    """
    assert (xy=="x" or xy=="y"), "xy function parameter must be 'x' or 'y'" 
    assert (isinstance(maxtem, int) and maxtem >=0), "maxtem must be integer >=0" 

    if (xy=="x"):
        print('PDTYPE x-split')
    else:
        print('PDTYPE y-split')

    # running through even modes
    for n1 in np.arange(0,maxtem+1,2):
        # running through odd modes
        for n2 in np.arange(1,maxtem+1,2) :
            c_n1n2= HG_split_diode_coefficient_numerical(n1,n2)
            if (xy=="x"):       
                print("{:2d} x {:2d} x {:.15g}".format(n1,n2,c_n1n2))
            else:
                print("x {:2d} x {:2d} {:.15g}".format(n1,n2,c_n1n2))
    print('END')

                
def HG_split_diode_coefficient_analytical(n1,n2):
    """ Using an analytical equation to compute beat coefficients
    betwween mode n1 and n2 on a split photo detector. Uses arbitrary
    precision (from gmpy) because otherwise errors get very large
    for n1,n2>30. This is for comparison with the numerical method
    HG_split_diode_coefficient_numerical only. """
    import gmpy
    
    temp=gmpy.mpq(0.0)
    for l in np.arange(0,n1/2+1):
        for k in np.arange(0,(n2-1)/2+1):
            temp +=  gmpy.mpq(pow((-1.0/4.0),l+k)*gmpy.mpz(factorial((n2+n1-1)/2-l-k))/gmpy.mpz((factorial(l)*factorial(k)*factorial(n1-2*l)*factorial(n2-2*k))))

    c_n1n2=gmpy.mpq(temp*gmpy.mpq(math.sqrt(2.0**(n1+n2)*gmpy.mpz(factorial(n1)*factorial(n2))/np.pi)))
    return float(gmpy.mpf(c_n1n2))
    
    
def HG_split_diode_coefficient_numerical(n1,n2):
    """ Compute beam coefficient between mode n1 and n2 on a split
    photo detector using numerical integration, tested up to n1,n2 = 40.
    This is primarily used by finesse_split_photodiode()"""
    A = 2.0 * np.sqrt(2.0/np.pi) * np.sqrt(1.0 / (2.0**(n1+n2) * factorial(n1) * factorial(n2)))
    f = lambda x: hermite(n1)(np.sqrt(2.0)*x) * math.exp(-2.0*x*x) * hermite(n2)(np.sqrt(2.0)*x)    
    c_n1n2= quad(f, 0.0, np.Inf, epsabs=1e-10, epsrel=1e-10, limit=200)
    return A*c_n1n2[0]
    
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133

def LG_bullseye_coefficients_numerical(l1,p1, l2, p2, w, r):
    """
    w = beam radisu [m]
    r = radius of inner element
    """
    if (l1 != l2):
        return 0.0
        
    l = np.abs(l1)
    # computer pre-factor of integral
    A = math.sqrt(factorial(p1)*factorial(p2)/(factorial(l+p1)*factorial(l+p2)))

    # definng integrand
    f = lambda x: x**l * laguerre(p1,l,x) * laguerre(p2,l,x) * math.exp(-x)
    # defining integration limit
    int_limit = 2.0 * r**2 / w**2
    # perform numerical integrations
    val1, res= quad(f, 0.0, int_limit, epsabs=1e-10, epsrel=1e-10, limit=500)
    val2, res= quad(f, int_limit, np.Inf, epsabs=1e-10, epsrel=1e-10, limit=500)
    return A*(val2-val1)

def HG_bullseye_coefficients_numerical(n1, m1, n2, m2, w, r):
    mparity=np.mod(m1+m2,2)
    nparity=np.mod(n1+n2,2)
                
    k = 0.0 + 0.0*1j
    if (mparity==0 and nparity==0):            
        a,p1,l1 = HG2LG(n1,m1)
        b,p2,l2 = HG2LG(n2,m2)
        for idx1, aa in enumerate(l1):
            for idx2, bb in enumerate(l2):
                if l1[idx1]==l2[idx2]:
                    c = LG_bullseye_coefficients_numerical(l1[idx1],p1[idx1], l2[idx2], p2[idx2], w, r)
                    k = k + a[idx1] * np.conj(b[idx2]) * c
    return float(np.real(k))
    
def finesse_bullseye_photodiode(maxtem, w=1.0, r=0.5887050, name='bullseye'):
    """Prints beat coefficients for HG modes on a split photo detector
    in the format for insertion into the kat.ini file of Finesse.
    Example use:
    finesse_split_photodiode(5, 0.01, 0.03, name="bull1:3")
    The default kat.ini numbers have been generated with:
    ...
    """
    assert (isinstance(maxtem, int) and maxtem >=0), "maxtem must be integer >=0" 

    print('PDTYPE {0}'.format(name))

    for n1 in np.arange(0,maxtem+1):
        for m1 in np.arange(0,maxtem+1):
            for n2 in np.arange(n1,maxtem+1):
                for m2 in np.arange(m1,maxtem+1):
                    c = HG_bullseye_coefficients_numerical(n1,m1, n2, m2, w, r)
                    if (np.abs(c)>10e-13 and not np.isnan(c)):
                        print("{:2d} {:2d} {:2d} {:2d} {:.15g}".format(n1,m1,n2,m2, c))
    print('END')