summaryrefslogtreecommitdiff
path: root/phasematching
diff options
context:
space:
mode:
Diffstat (limited to 'phasematching')
-rw-r--r--phasematching/WGM_lib.py824
1 files changed, 444 insertions, 380 deletions
diff --git a/phasematching/WGM_lib.py b/phasematching/WGM_lib.py
index 3e19fee..ee03b3d 100644
--- a/phasematching/WGM_lib.py
+++ b/phasematching/WGM_lib.py
@@ -1,380 +1,444 @@
-#!/usr/bin/env mpython_q
-def ntest(lda,T,pol=0,MgO=0,MgO_th=5.0,E = 0): return 1.5
-
-def nLNO1(lda,T,pol,MgO=0,MgO_th=5.0,E = 0):
- """
- This is Sellmeyer eq. for Lithium Niobate from "Influence of the defect structure on the refractive indices
- of undoped and Mg-doped lithium niobate", U.Schlarb and K.Betzler, Phys. Rev. B 50, 751 (1994).
- Temperature in Centigrades, wavelength in microns. pol=0 "ordinary", pol=1 "extraordinary" MgO
- =0 MgO concentration in % MgO_th the threshold concentration. External electric field E in V/cm
- """
- import math
- if pol: # extraordinary
- w0 = 218.203
- mu = 6.4047E-6
- A0 = 3.9466E-5
- A_NbLi = 23.727E-8
- A_Mg = 7.6243E-8
- A_IR = 3.0998E-8
- A_UV = 2.6613
- rr = 33E-10 # r_33 cm/V
- else: # ordinary
- w0 = 223.219
- mu = 1.1082E-6
- A0 = 4.5312E-5
- A_NbLi = -1.4464E-8
- A_Mg = -7.3548E-8
- A_IR = 3.6340E-8
- A_UV = 2.6613
- rr = 11E-10 # r_13 cm/V
- B = A0+A_Mg*MgO
- if MgO < MgO_th: B+=(MgO_th-MgO)*A_NbLi
- def f(T):
- return (T + 273)**2 + 4.0238E5*(1/math.tanh(261.6/(T + 273)) - 1)
- def n(lda,T):
- return math.sqrt(B/((w0+(f(T)-f(24.5))*mu)**(-2)-lda**(-2))-A_IR*lda**2+A_UV)
- return n(lda,T)+0.5*rr*E*n(lda,T)**3
-
-def nBBO(lda,T,pol=0,MgO=0,MgO_th=5.0,E=0):
- """
- This is Sellmeyer eq. for BBO from http://refractiveindex.info/
- (Handbook of Optics, 3rd edition, Vol. 4. McGraw-Hill 2009);
- T is the angle (in degrees) between the optical axis and polarization:
- T = 0 extraordinary, T = 90 ordinary.
- """
- import math
- a = 2.7405,2.3730 # ordinary, extraordinary
- b = 0.0184,0.0128
- c = 0.0179,0.0156
- d = 0.0155,0.0044
- lda = lda/1000.0 # nm -> microns
- no = math.sqrt(a[0]-d[0]*lda*lda+b[0]/(lda*lda-c[0]))
- ne = math.sqrt(a[1]-d[1]*lda*lda+b[1]/(lda*lda-c[1]))
- o = math.sin(math.pi*T/180)/no
- e = math.cos(math.pi*T/180)/ne
- return (o*o+e*e)**(-0.5)
-
-def nBBO1(lda,T,pol=0,MgO=0,MgO_th=5.0,E=0):
- """
- This is Sellmeyer eq. for BBO from Optics Communications 184 (2000) 485-491;
- T is the angle (in degrees) between the optical axis and polarization:
- T = 0 extraordinary, T = 90 ordinary.
- """
- import math
- a = 2.7359,2.3753 # ordinary, extraordinary
- b = 0.01878,0.01224
- c = 0.01822,0.01667
- d = 0.01471,0.01627
- x = 0.0006081,0.0005716
- y = 0.00006740,0.00006305
- lda = lda/1000.0 # nm -> microns
- no = math.sqrt(a[0]-d[0]*lda*lda+b[0]/(lda*lda-c[0])+x[0]*lda**4-y[0]*lda**6)
- ne = math.sqrt(a[1]-d[1]*lda*lda+b[1]/(lda*lda-c[1])+x[1]*lda**4-y[1]*lda**6)
- o = math.sin(math.pi*T/180)/no
- e = math.cos(math.pi*T/180)/ne
- return (o*o+e*e)**(-0.5)
-
-def FSR_simple(R,lda,T,pol,L,q = 1,MgO=0,MgO_th=5.0,E = 0, n = nLNO1):
- import math
- #from WGM_lib import nLNO1
- from scipy.special import ai_zeros
- q -= 1
- AiRoots = -ai_zeros(40)[0]
- a = 29.9792/(2*math.pi*R*n(lda,T,pol,MgO,MgO_th,E))
- b = 1+0.5*AiRoots[q]*((L/2.0+0.5)**(1.0/3.0)-(L/2.0-0.5)**(1.0/3.0))
- return a*b
-
-def Get_frac_L(R,r,lda,T,pol,q = 1,p = 0,MgO=0,MgO_th=5.0,E = 0, n = nLNO1):
- """
- This iteratively solves the WGM dispersion equation for a spheroid at the target wavelength lda
- and returs the fractional orbital number L
- pol=0 "ordinary", pol=1 "extraordinary" MgO =0 MgO concentration in %
- MgO_th the threshold concentration. External electric field E in V/cm
- R, r are big and small radia in mm
- """
- import math
- #from WGM_lib import nLNO1
- from scipy.special import ai_zeros
- q -= 1
- r = math.sqrt(r*R) # redefine the rim radius to spheroid semi-axis
- AiRoots = -ai_zeros(40)[0]
- freq_term = 2*math.pi*R*n(lda,T,pol,MgO,MgO_th,E)*1E7/lda
- geom_term = (2*p*(R-r)+R)/2.0/r
- pol_term = n(lda,T,pol,MgO,MgO_th,E)**(1-2*pol)/math.sqrt(n(lda,T,pol,MgO,MgO_th,E)**2-1)
- L = freq_term
- dL = 1
- while math.fabs(dL) >= 1E-9:
- L1 = freq_term - AiRoots[q]*(L/2.0)**(1.0/3.0) - geom_term + pol_term\
- - 3*AiRoots[q]**2*(L/2.0)**(-1.0/3.0)/20
- dL = L1-L
- L = L1
- return L
-
-def WGM_freq(R,r,L,T,pol,q = 1,p = 0,MgO=0,MgO_th=5.0,E = 0, n = nLNO1):
- """
- This iteratively solves the WGM dispersion equation for a spheroid at the target
- orbital number L and returs the wavelength (nm) and frequency (GHz)
- pol=0 "ordinary", pol=1 "extraordinary" MgO =0 MgO concentration in %
- MgO_th the threshold concentration. External electric field E in V/cm
- R, r are big and small radia in mm
- """
- import math
- #from WGM_lib import nLNO1
- from scipy.special import ai_zeros
- q -= 1
- r = math.sqrt(r*R) # redefine the rim radius to spheroid semi-axis
- AiRoots = -ai_zeros(40)[0] # >0
- geom_term = (2*p*(R-r)+R)/2.0/r
- def nm2GHz(lda): return 2.99792E8/lda
- def wl1(lda):
- freq_term = 2*math.pi*R*n(lda,T,pol,MgO,MgO_th,E)*1E7
- pol_term = n(lda,T,pol,MgO,MgO_th,E)**(1-2*pol)/math.sqrt(n(lda,T,pol,MgO,MgO_th,E)**2-1)
- return freq_term/(L+AiRoots[q]*(L/2.0)**(1.0/3.0)\
- +geom_term-pol_term+3*AiRoots[q]**2*(L/2.0)**(-1.0/3.0)/20)
- lda0 = 2*math.pi*R*n(1000,T,pol,MgO,MgO_th,E)*1E7/L # initial wavelength guess based on n(1000 nm)
- lda = wl1(lda0)
- #print lda0, lda
- df = 1
- while math.fabs(df) > 1E-6: #1 kHz accuracy
- lda1 = wl1(lda)
- df = nm2GHz(lda1)-nm2GHz(lda)
- lda = lda1
- #print df,lda1,n(lda,T,pol,MgO,MgO_th,E)
- return lda, nm2GHz(lda)
-
-def Gaunt(L1,m1,L2,m2,L3,m3):
- """
- Gaunt formula for Legendre polynomials overlap. Large factorials are treated in Stirling's approximation.
- """
- import math
- def fa(n):
- if n<=10: fac = float(math.factorial(n))
- else: fac = math.sqrt(2*math.pi*n)*(n/math.e)**n
- return fac
- params = [(L1,m1),(L2,m2),(L3,m3)]
- params = sorted(params,key=lambda list: list[1])
- l,u = params[-1]
- prms = params[:-1]
- prms = sorted(prms,key=lambda list: list[0])
- n,w = prms[0]
- m,v = prms[1]
- if not v+w == u: return 0
- else:
- s = (l+m+n)/2.0
- if not s-round(s) == 0 : return 0
- elif m+n <l : return 0
- elif m-n >l : return 0
- else:
- p = int(max(0, n-m-u))
- q = int(min(n+m-u,l-u,n-w))
- a = (-1)**(s-m-w)*math.sqrt((2*l+1)*(2*m+1)*(2*n+1)*fa(l-u)*fa(n-w)/2.0/fa(m-v))/4/math.pi/math.pi
- b = a*math.sqrt(math.sqrt(2*math.pi*(m+v)*(n+w)/(l+u))/(s-m)/(s-n))/2.0/math.pi
- c = b*math.exp(0.5*(m+n-l+u-v-w))*(s-l)**(l-s)
- d = c*((n+w)/float(l+n-m))**(s-m)*((n+w)/2.0)**((w+m-l)/2.0)
- t0 = d*2**((w-n)/2.0)*((m+v)/float(l+m-n))**((l+m-n)/2.0)*(m+v)**((v-l+n)/2.0)*(l+m-n)**(l-u)
- sum = 0
- for t in range (p,q+1):
- f1 = (-1)**t*fa(m+n-u-t)/fa(t)/fa(l-u-t)/fa(n-w-t)
- f2 = f1*((l+u+t)/float(m-n+u+t))**(0.5+t)
- f3 = f2*((l+u+t)/float(l+u))**(0.5*l+0.5*u)
- f4 = f3*((l+u+t)/float(l+m+n))**((l+m+n)/2.0)*(l+u+t)**((u-m-n)/2.0)
- sum+= f4*((l+m-n)/float(u+m-n+t))**(u+m-n)
- return t0*sum
-
-def Gaunt_low(L1,m1,L2,m2,L3,m3):
- """
- Gaunt formula for Legendre polynomials overlap computed exactly. Good only for low orders.
- """
- import math
- def fa(n):
- if n<=10: fac = float(math.factorial(n))
- else: fac = math.sqrt(2*math.pi*n)*(n/math.e)**n
- return fac
- def fa_ratio(n1, n2): #calculates n1!/n2!
- p = 1.0
- for j in range (min(n1, n2)+1, max(n1, n2)+1):
- p = p*j
- if n1 >= n2: return p
- else: return 1/p
- params = [(L1,m1),(L2,m2),(L3,m3)]
- params = sorted(params,key=lambda list: list[1])
- l,u = params[-1]
- prms = params[:-1]
- prms = sorted(prms,key=lambda list: list[0])
- n,w = prms[0]
- m,v = prms[1]
- if not v+w == u: return 0
- else:
- s = (l+m+n)/2.0
- if not s-round(s) == 0 : return 0
- elif m+n <l : return 0
- elif m-n >l : return 0
- else:
- s = int(s)
- p = int(max(0, n-m-u))
- q = int(min(n+m-u,l-u,n-w))
- #print l,m,n,s
- #print fa(l-u),fa(n-w), fa(m-v),fa(s-l)
- a1 = (-1)**(s-m-w)*math.sqrt((2*l+1)*(2*m+1)*(2*n+1)*fa(l-u)*fa(n-w)/fa(m-v)/math.pi)/fa(s-l)/float(2*s+1)/4.0
- a2 = math.sqrt(fa(m+v)*fa(n+w)/fa(l+u))*fa(s)/fa(s-m)/fa(s-n)
- sum = 0
- for t in range (p,q+1):
- b1 = (-1)**t*fa(m+n-u-t)/fa(t)/fa(l-u-t)/fa(n-w-t)
- b2 = fa_ratio(l+u+t,2*s)*fa_ratio(2*s-2*n,m-n+u+t)
- sum+= b1*b2
- #print a1, a2, sum
- return a1*a2*sum
-
-def RadOvlp(R,L1,q1,L2,q2,L3,q3):
- """
- Calculates radial overlap of three Airy functions
- """
- import math
- from scipy.special import airy, ai_zeros
- from scipy.integrate import quad
- AiRoots = -ai_zeros(40)[0]
- q1, q2, q3 = q1-1, q2-1, q3-1
- def nk(L,q): return L*(1+AiRoots[q]*(2*L**2)**(-1.0/3.0)+0.5/L)
- def AiArg(x,L,q): return -AiRoots[q]*(L+0.5-nk(L,q)*x)/(L+0.5-nk(L,q))
- def RadFunc(x,L,q): return airy(AiArg(x,L,q))[0]/math.sqrt(x)
- L0 = min(L1,L2,L3)
- q0 = max(q1,q2,q3)
- Rmin = (L0-3*(L0/2.0)**(1.0/3.0))/nk(L0,q0)
- N1 = R**3*quad(lambda x: RadFunc(x,L1,q1)**2, Rmin, 1)[0]
- N2 = R**3*quad(lambda x: RadFunc(x,L2,q2)**2, Rmin, 1)[0]
- N3 = R**3*quad(lambda x: RadFunc(x,L3,q3)**2, Rmin, 1)[0]
- def Core(x): return RadFunc(x,L1,q1)*RadFunc(x,L2,q2)*RadFunc(x,L3,q3)*x*x
- return R**3*quad(lambda x: Core(x), Rmin, 1)[0]/math.sqrt(N1*N2*N3)
-
-def AngOvlp(L1,p1,L2,p2,L3,p3,limit=0.05):
- """
- Calculates angular overlap of three spherical functions. Default limit is sufficient for up to p = 10
- """
- import math
- from scipy.special import hermite
- from scipy.integrate import quad
- def HG(p,L,x): return (2**L*math.factorial(p))**-0.5*hermite(p)(x*math.sqrt(L))*(L/math.pi)**0.25*math.exp(-L*x**2/2.0)
- def HGnorm(p,L,x): return HG(p,L,x)*(4*math.pi*quad(lambda x: HG(p,L,x)**2, 0, limit*(14000.0/L)**0.5)[0])**(-0.5)
- return 4*math.pi*quad(lambda x: HGnorm(p1,L1,x)*HGnorm(p2,L2,x)*HGnorm(p3,L3,x), 0, limit*(14000.0/min(L1,L2,L3))**0.5)[0]
-
-def T_phasematch(L1,p1,q1,L2,p2,q2,L3,p3,q3,R,r,T0=70,MgO=0,MgO_th=5.0,E = 0,Tstep=0.1, n = nLNO1):
- '''
- Finds 1/wl2 + 1/wl2 = 1/wl3, L1+L2 -> L3 phase matching temperature down to the accuracy deltaf (GHz)
- based on the initial guess T0 with initial search step Tstep
- Requires importing Sellmeyer equations. This program is unaware of the orbital selection rules!
- '''
- import math
- from WGM_lib import WGM_freq
- def Df(T): # freqency detuning in GHz
- lda_1, freq_1 = WGM_freq(R,r,L1,T,0,q1,p1,MgO,MgO_th,E,n)
- lda_2, freq_2 = WGM_freq(R,r,L2,T,0,q2,p2,MgO,MgO_th,E,n)
- lda_3, freq_3 = WGM_freq(R,r,L3,T,1,q3,p3,MgO,MgO_th,E,n)
- return freq_3-freq_1-freq_2
- T = T0
- df = Df(T)
- while math.fabs(df)>1E-6: #1 kHz accuracy
- #print "T = ",T," df = ", df,
- df1 = Df(T+Tstep)
- ddfdt = (df1-df)/Tstep
- dT = df/ddfdt
- T -= dT
- Tstep = math.fabs(dT/2)
- df = Df(T)
- #print " dT = ", dT, " T -> ",T, " df -> ", df
- if T<-200 or T> 500: break
- return T
-
-def T_phasematch_blk(wl1,wl2, T0=70,MgO=0,MgO_th=5.0,E = 0,Tstep=0.1, n = nLNO1):
- '''
- Finds 1/wl2 + 1/wl2 = 1/wl3 (wl in nanometers), collinear bulk phase matching temperature
- based on the initial guess T0 with initial search step Tstep
- Requires importing Sellmeyer equations.
- '''
- import math
- wl3 = 1/(1.0/wl1+1.0/wl2)
- def Dk(T):
- k1 = 2*math.pi*n(wl1,T,0,MgO,MgO_th,E)*1e7/wl1 # cm^-1
- k2 = 2*math.pi*n(wl2,T,0,MgO,MgO_th,E)*1e7/wl2
- k3 = 2*math.pi*n(wl3,T,1,MgO,MgO_th,E)*1e7/wl3
- return k3-k2-k1
- T = T0
- dk = Dk(T)
- while math.fabs(dk)>0.1: #10 cm beat length
- dk1 = Dk(T+Tstep)
- ddkdt = (dk1-dk)/Tstep
- dT = dk/ddkdt
- T -= dT
- Tstep = math.fabs(dT/2)
- dk = Dk(T)
- #print " dT = ", dT, " T -> ",T, " dk -> ", dk
- if T<-200 or T> 500:
- print "No phase matching!"
- break
- return T
-
-def SH_phasematch_blk(angle, wl0=1000,step=1.01, n = nBBO):
- '''
- Finds the fundamental wavelength (in nm) for frequency doubling o+o->e at a given angle between the optical axis and the e.
- Requires importing Sellmeyer equations.
- '''
- import math
- def Dk(wl):
- k1 = 2*math.pi*n(wl,90)*1e7/wl # ordinary cm^-1
- k2 = 4*math.pi*n(wl/2.0,angle)*1e7/wl # extraordinary with angle
- return k2-2*k1
- wl = wl0
- dk = Dk(wl)
- while math.fabs(dk)>0.1: #10 cm beat length
- dk1 = Dk(wl*step)
- ddkdt = (dk1-dk)/step
- dwl = dk/ddkdt
- wl = wl/dwl
- #step = min(step,math.fabs(dwl/2) )
- dk = Dk(wl)
- print " dwl = ", dwl, " wl -> ",wl, " dk -> ", dk
- return wl
-
-
-def Phasematch(wlP,p1,q1,p2,q2,p3,q3,R,r,T,MgO=0,MgO_th=5.0,E = 0, n = nLNO1):
- '''
-
- NOT FULLY FUNCTIONAL!!
-
- Finds wl1 and wl2 such that 1/wl2 + 1/wl2 = 1/wlP and L1-p1+L2-p2 = L3-p3 at a given temperature T.
- Requires importing Sellmeyer equations. R is evaluated at T. This program is unaware of the other orbital
- selection rules! L1, L2, L3 are allowed to be non-integer.
- '''
- import math
- from WGM_lib import Get_frac_L
- def Dm(wl): # m3-m1-m2, orbital detuning
- L3 = Get_frac_L(R,r,wlP,T,1,q3,p3,MgO,MgO_th)
- L1 = Get_frac_L(R,r,wl,T,0,q1,p1,MgO,MgO_th)
- wl2 = 1.0/(1.0/wlP-1.0/wl)
- L2 = Get_frac_L(R,r,wl2,T,0,q2,p2,MgO,MgO_th)
- return L2-L3+p3+L1-p1-p2
- wl = 2*wlP #try at degeneracy
- dm = Dm(wl)
- print "dm = ", dm
- ddm = 0
- wlstep = 1.0 # chosing an appropriate wavelength step so the m variation is not too large or too small
- while math.fabs(ddm)<1:
- wlstep = 2*wlstep
- ddm = Dm(wl+wlstep)-dm
- print "wlstep = ", wlstep, " initial slope = ", ddm/wlstep
- while math.fabs(dm)>1E-6: # equiv *FSR accuracy
- dm1 = Dm(wl+wlstep)
- print "dm1 = ", dm1
- slope = (dm1-dm)/wlstep
- print "slope = ", slope, "nm^-1"
- dwl = dm/slope
- wl -= dwl
- print "new wl = ", wl
- wlstep = min(math.fabs(dwl/2), wlstep)
- dm = Dm(wl)
- print " d wl = ", dwl, " wl -> ",wl, " dm -> ", dm
- if math.fabs(dm)> Get_frac_L(R,r,2*wlP,T,0,q1,p1,MgO,MgO_th): break
- L1 = Get_frac_L(R,r,wl,T,0,q1,p1,MgO,MgO_th)
- wl2 = 1.0/(1.0/wlP-1.0/wl)
- L2 = Get_frac_L(R,r,wl2,T,0,q2,p2,MgO,MgO_th)
- return L1, L2, wl, wl2
-
+#!/usr/bin/env mpython_q
+def ntest(lda,T,pol=0,MgO=0,MgO_th=5.0,E = 0): return 1.5
+
+def nLNO1(lda,T,pol,MgO=0,MgO_th=5.0,E = 0):
+ """
+ This is Sellmeyer eq. for Lithium Niobate from "Influence of the defect structure on the refractive indices
+ of undoped and Mg-doped lithium niobate", U.Schlarb and K.Betzler, Phys. Rev. B 50, 751 (1994).
+ Temperature in Centigrades, wavelength in microns. pol=0 "ordinary", pol=1 "extraordinary" MgO
+ =0 MgO concentration in % MgO_th the threshold concentration. External electric field E in V/cm
+ """
+ import math
+ if pol: # extraordinary
+ w0 = 218.203
+ mu = 6.4047E-6
+ A0 = 3.9466E-5
+ A_NbLi = 23.727E-8
+ A_Mg = 7.6243E-8
+ A_IR = 3.0998E-8
+ A_UV = 2.6613
+ rr = 33E-10 # r_33 cm/V
+ else: # ordinary
+ w0 = 223.219
+ mu = 1.1082E-6
+ A0 = 4.5312E-5
+ A_NbLi = -1.4464E-8
+ A_Mg = -7.3548E-8
+ A_IR = 3.6340E-8
+ A_UV = 2.6613
+ rr = 11E-10 # r_13 cm/V
+ B = A0+A_Mg*MgO
+ if MgO < MgO_th: B+=(MgO_th-MgO)*A_NbLi
+ def f(T):
+ return (T + 273)**2 + 4.0238E5*(1/math.tanh(261.6/(T + 273)) - 1)
+ def n(lda,T):
+ return math.sqrt(B/((w0+(f(T)-f(24.5))*mu)**(-2)-lda**(-2))-A_IR*lda**2+A_UV)
+ return n(lda,T)+0.5*rr*E*n(lda,T)**3
+
+def nLNO1_disp(lda,T,pol,MgO=0,MgO_th=5.0):
+ """
+ dispersion dn/dlda (in nm) for Lithium Niobate dispersion nLNO1
+ """
+ import math
+ if pol: # extraordinary
+ w0 = 218.203
+ mu = 6.4047E-6
+ A0 = 3.9466E-5
+ A_NbLi = 23.727E-8
+ A_Mg = 7.6243E-8
+ A_IR = 3.0998E-8
+ A_UV = 2.6613
+ else: # ordinary
+ w0 = 223.219
+ mu = 1.1082E-6
+ A0 = 4.5312E-5
+ A_NbLi = -1.4464E-8
+ A_Mg = -7.3548E-8
+ A_IR = 3.6340E-8
+ A_UV = 2.6613
+ B = A0+A_Mg*MgO
+ if MgO < MgO_th: B+=(MgO_th-MgO)*A_NbLi
+ def f(T):
+ return (T + 273)**2 + 4.0238E5*(1/math.tanh(261.6/(T + 273)) - 1)
+ def n(lda,T):
+ return math.sqrt(B/((w0+(f(T)-f(24.5))*mu)**(-2)-lda**(-2))-A_IR*lda**2+A_UV)
+ return -(B*lda**(-3)*((w0+(f(T)-f(24.5))*mu)**(-2)-lda**(-2))**(-2)+A_IR*lda)/n(lda,T)
+
+def nBBO(lda,T,pol=0,MgO=0,MgO_th=5.0,E=0):
+ """
+ This is Sellmeyer eq. for BBO from http://refractiveindex.info/
+ (Handbook of Optics, 3rd edition, Vol. 4. McGraw-Hill 2009);
+ T is the angle (in degrees) between the optical axis and polarization:
+ T = 0 extraordinary, T = 90 ordinary.
+ """
+ import math
+ a = 2.7405,2.3730 # ordinary, extraordinary
+ b = 0.0184,0.0128
+ c = 0.0179,0.0156
+ d = 0.0155,0.0044
+ lda = lda/1000.0 # nm -> microns
+ no = math.sqrt(a[0]-d[0]*lda*lda+b[0]/(lda*lda-c[0]))
+ ne = math.sqrt(a[1]-d[1]*lda*lda+b[1]/(lda*lda-c[1]))
+ o = math.sin(math.pi*T/180)/no
+ e = math.cos(math.pi*T/180)/ne
+ return (o*o+e*e)**(-0.5)
+
+def nBBO1(lda,T,pol=0,MgO=0,MgO_th=5.0,E=0):
+ """
+ This is Sellmeyer eq. for BBO from Optics Communications 184 (2000) 485-491;
+ T is the angle (in degrees) between the optical axis and polarization:
+ T = 0 extraordinary, T = 90 ordinary.
+ """
+ import math
+ a = 2.7359,2.3753 # ordinary, extraordinary
+ b = 0.01878,0.01224
+ c = 0.01822,0.01667
+ d = 0.01471,0.01627
+ x = 0.0006081,0.0005716
+ y = 0.00006740,0.00006305
+ lda = lda/1000.0 # nm -> microns
+ no = math.sqrt(a[0]-d[0]*lda*lda+b[0]/(lda*lda-c[0])+x[0]*lda**4-y[0]*lda**6)
+ ne = math.sqrt(a[1]-d[1]*lda*lda+b[1]/(lda*lda-c[1])+x[1]*lda**4-y[1]*lda**6)
+ o = math.sin(math.pi*T/180)/no
+ e = math.cos(math.pi*T/180)/ne
+ return (o*o+e*e)**(-0.5)
+
+def FSR_simple(R,lda,T,pol,L,q = 1,MgO=0,MgO_th=5.0,E = 0, n = nLNO1):
+ import math
+ #from WGM_lib import nLNO1
+ from scipy.special import ai_zeros
+ q -= 1
+ AiRoots = -ai_zeros(40)[0]
+ a = 29.9792/(2*math.pi*R*n(lda,T,pol,MgO,MgO_th,E))
+ b = 1+0.5*AiRoots[q]*((L/2.0+0.5)**(1.0/3.0)-(L/2.0-0.5)**(1.0/3.0))
+ return a*b
+
+def Get_frac_L(R,r,lda,T,pol,q = 1,p = 0,MgO=0,MgO_th=5.0,E = 0, n = nLNO1):
+ """
+ This iteratively solves the WGM dispersion equation for a spheroid at the target wavelength lda
+ and returs the fractional orbital number L
+ pol=0 "ordinary", pol=1 "extraordinary" MgO =0 MgO concentration in %
+ MgO_th the threshold concentration. External electric field E in V/cm
+ R, r are big and small radia in mm
+ """
+ import math
+ #from WGM_lib import nLNO1
+ from scipy.special import ai_zeros
+ q -= 1
+ r = math.sqrt(r*R) # redefine the rim radius to spheroid semi-axis
+ AiRoots = -ai_zeros(40)[0]
+ freq_term = 2*math.pi*R*n(lda,T,pol,MgO,MgO_th,E)*1E7/lda
+ geom_term = (2*p*(R-r)+R)/2.0/r
+ pol_term = n(lda,T,pol,MgO,MgO_th,E)**(1-2*pol)/math.sqrt(n(lda,T,pol,MgO,MgO_th,E)**2-1)
+ L = freq_term
+ dL = 1
+ while math.fabs(dL) >= 1E-9:
+ L1 = freq_term - AiRoots[q]*(L/2.0)**(1.0/3.0) - geom_term + pol_term\
+ - 3*AiRoots[q]**2*(L/2.0)**(-1.0/3.0)/20
+ dL = L1-L
+ L = L1
+ return L
+
+def WGM_freq(R,r,L,T,pol,q = 1,p = 0,MgO=0,MgO_th=5.0,E = 0, n = nLNO1):
+ """
+ This iteratively solves the WGM dispersion equation for a spheroid at the target
+ orbital number L and returs the wavelength (nm) and frequency (GHz)
+ pol=0 "ordinary", pol=1 "extraordinary" MgO =0 MgO concentration in %
+ MgO_th the threshold concentration. External electric field E in V/cm
+ R, r are big and small radia in mm
+ """
+ import math
+ #from WGM_lib import nLNO1
+ from scipy.special import ai_zeros
+ q -= 1
+ r = math.sqrt(r*R) # redefine the rim radius to spheroid semi-axis
+ AiRoots = -ai_zeros(40)[0] # >0
+ geom_term = (2*p*(R-r)+R)/2.0/r
+ def nm2GHz(lda): return 2.99792E8/lda
+ def wl1(lda):
+ freq_term = 2*math.pi*R*n(lda,T,pol,MgO,MgO_th,E)*1E7
+ pol_term = n(lda,T,pol,MgO,MgO_th,E)**(1-2*pol)/math.sqrt(n(lda,T,pol,MgO,MgO_th,E)**2-1)
+ return freq_term/(L+AiRoots[q]*(L/2.0)**(1.0/3.0)\
+ +geom_term-pol_term+3*AiRoots[q]**2*(L/2.0)**(-1.0/3.0)/20)
+ lda0 = 2*math.pi*R*n(1000,T,pol,MgO,MgO_th,E)*1E7/L # initial wavelength guess based on n(1000 nm)
+ lda = wl1(lda0)
+ #print lda0, lda
+ df = 1
+ while math.fabs(df) > 1E-6: #1 kHz accuracy
+ lda1 = wl1(lda)
+ df = nm2GHz(lda1)-nm2GHz(lda)
+ lda = lda1
+ #print df,lda1,n(lda,T,pol,MgO,MgO_th,E)
+ return lda, nm2GHz(lda)
+
+def Gaunt(L1,m1,L2,m2,L3,m3):
+ """
+ Gaunt formula for Legendre polynomials overlap. Large factorials are treated in Stirling's approximation.
+ """
+ import math
+ def fa(n):
+ if n<=10: fac = float(math.factorial(n))
+ else: fac = math.sqrt(2*math.pi*n)*(n/math.e)**n
+ return fac
+ params = [(L1,m1),(L2,m2),(L3,m3)]
+ params = sorted(params,key=lambda list: list[1])
+ l,u = params[-1]
+ prms = params[:-1]
+ prms = sorted(prms,key=lambda list: list[0])
+ n,w = prms[0]
+ m,v = prms[1]
+ if not v+w == u: return 0
+ else:
+ s = (l+m+n)/2.0
+ if not s-round(s) == 0 : return 0
+ elif m+n <l : return 0
+ elif m-n >l : return 0
+ else:
+ p = int(max(0, n-m-u))
+ q = int(min(n+m-u,l-u,n-w))
+ a = (-1)**(s-m-w)*math.sqrt((2*l+1)*(2*m+1)*(2*n+1)*fa(l-u)*fa(n-w)/2.0/fa(m-v))/4/math.pi/math.pi
+ b = a*math.sqrt(math.sqrt(2*math.pi*(m+v)*(n+w)/(l+u))/(s-m)/(s-n))/2.0/math.pi
+ c = b*math.exp(0.5*(m+n-l+u-v-w))*(s-l)**(l-s)
+ d = c*((n+w)/float(l+n-m))**(s-m)*((n+w)/2.0)**((w+m-l)/2.0)
+ t0 = d*2**((w-n)/2.0)*((m+v)/float(l+m-n))**((l+m-n)/2.0)*(m+v)**((v-l+n)/2.0)*(l+m-n)**(l-u)
+ sum = 0
+ for t in range (p,q+1):
+ f1 = (-1)**t*fa(m+n-u-t)/fa(t)/fa(l-u-t)/fa(n-w-t)
+ f2 = f1*((l+u+t)/float(m-n+u+t))**(0.5+t)
+ f3 = f2*((l+u+t)/float(l+u))**(0.5*l+0.5*u)
+ f4 = f3*((l+u+t)/float(l+m+n))**((l+m+n)/2.0)*(l+u+t)**((u-m-n)/2.0)
+ sum+= f4*((l+m-n)/float(u+m-n+t))**(u+m-n)
+ return t0*sum
+
+def Gaunt_low(L1,m1,L2,m2,L3,m3):
+ """
+ Gaunt formula for Legendre polynomials overlap computed exactly. Good only for low orders.
+ """
+ import math
+ def fa(n):
+ if n<=10: fac = float(math.factorial(n))
+ else: fac = math.sqrt(2*math.pi*n)*(n/math.e)**n
+ return fac
+ def fa_ratio(n1, n2): #calculates n1!/n2!
+ p = 1.0
+ for j in range (min(n1, n2)+1, max(n1, n2)+1):
+ p = p*j
+ if n1 >= n2: return p
+ else: return 1/p
+ params = [(L1,m1),(L2,m2),(L3,m3)]
+ params = sorted(params,key=lambda list: list[1])
+ l,u = params[-1]
+ prms = params[:-1]
+ prms = sorted(prms,key=lambda list: list[0])
+ n,w = prms[0]
+ m,v = prms[1]
+ if not v+w == u: return 0
+ else:
+ s = (l+m+n)/2.0
+ if not s-round(s) == 0 : return 0
+ elif m+n <l : return 0
+ elif m-n >l : return 0
+ else:
+ s = int(s)
+ p = int(max(0, n-m-u))
+ q = int(min(n+m-u,l-u,n-w))
+ #print l,m,n,s
+ #print fa(l-u),fa(n-w), fa(m-v),fa(s-l)
+ a1 = (-1)**(s-m-w)*math.sqrt((2*l+1)*(2*m+1)*(2*n+1)*fa(l-u)*fa(n-w)/fa(m-v)/math.pi)/fa(s-l)/float(2*s+1)/4.0
+ a2 = math.sqrt(fa(m+v)*fa(n+w)/fa(l+u))*fa(s)/fa(s-m)/fa(s-n)
+ sum = 0
+ for t in range (p,q+1):
+ b1 = (-1)**t*fa(m+n-u-t)/fa(t)/fa(l-u-t)/fa(n-w-t)
+ b2 = fa_ratio(l+u+t,2*s)*fa_ratio(2*s-2*n,m-n+u+t)
+ sum+= b1*b2
+ #print a1, a2, sum
+ return a1*a2*sum
+
+def RadOvlp(R,L1,q1,L2,q2,L3,q3):
+ """
+ Calculates radial overlap of three Airy functions
+ """
+ import math
+ from scipy.special import airy, ai_zeros
+ from scipy.integrate import quad
+ AiRoots = -ai_zeros(40)[0]
+ q1, q2, q3 = q1-1, q2-1, q3-1
+ def nk(L,q): return L*(1+AiRoots[q]*(2*L**2)**(-1.0/3.0))
+ def AiArg(x,L,q): return -AiRoots[q]*(L-nk(L+0.5,q)*x)/(L-nk(L+0.5,q))
+ def RadFunc(x,L,q): return airy(AiArg(x,L,q))[0]/math.sqrt(x)
+ L0 = min(L1,L2,L3)
+ q0 = max(q1,q2,q3)+1
+ Rmin = 1-5*math.sqrt(q0)*L0**(-2.0/3.0)
+ N1 = R**3*quad(lambda x: x*x*RadFunc(x,L1,q1)**2, Rmin, 1)[0]
+ N2 = R**3*quad(lambda x: x*x*RadFunc(x,L2,q2)**2, Rmin, 1)[0]
+ N3 = R**3*quad(lambda x: x*x*RadFunc(x,L3,q3)**2, Rmin, 1)[0]
+ #print N1, N2, N3
+ def Core(x): return RadFunc(x,L1,q1)*RadFunc(x,L2,q2)*RadFunc(x,L3,q3)*x*x
+ return R**3*quad(lambda x: Core(x), Rmin, 1)[0]/math.sqrt(N1*N2*N3)
+
+def AngOvlp(L1,p1,L2,p2,L3,p3,limit=3):
+ """
+ Calculates angular overlap of three spherical functions. Default limit is sufficient for up to p = 10
+ """
+ import math
+ from scipy.special import hermite
+ from scipy.integrate import quad
+ #def HG(p,L,x): return (2**L*math.factorial(p))**-0.5*hermite(p)(x*math.sqrt(L))*(L/math.pi)**0.25*math.exp(-L*x**2/2.0)
+ def HG(p,L,x): return hermite(p)(x*math.sqrt(L))*math.exp(-L*x**2/2.0) # removed a large factor that would drop out in normalization anyway
+ def HGnorm(p,L,x): return HG(p,L,x)*(4*math.pi*quad(lambda x: HG(p,L,x)**2, 0, limit*((p+1)*math.pi/L/2)**0.5)[0])**(-0.5)
+ return 2*quad(lambda x: HGnorm(p1,L1,x)*HGnorm(p2,L2,x)*HGnorm(p3,L3,x), 0, limit*((max(p1,p2,p3)+1)*math.pi/min(L1,L2,L3)/2)**0.5)[0]
+
+def T_phasematch(L1,p1,q1,L2,p2,q2,L3,p3,q3,R,r,T0=70,MgO=0,MgO_th=5.0,E = 0,Tstep=0.1, n = nLNO1):
+ '''
+ Finds 1/wl2 + 1/wl2 = 1/wl3, L1+L2 -> L3 phase matching temperature down to the accuracy deltaf (GHz)
+ based on the initial guess T0 with initial search step Tstep
+ Requires importing Sellmeyer equations. This program is unaware of the orbital selection rules!
+ '''
+ import math
+ from WGM_lib import WGM_freq
+ def Df(T): # freqency detuning in GHz
+ lda_1, freq_1 = WGM_freq(R,r,L1,T,0,q1,p1,MgO,MgO_th,E,n)
+ lda_2, freq_2 = WGM_freq(R,r,L2,T,0,q2,p2,MgO,MgO_th,E,n)
+ lda_3, freq_3 = WGM_freq(R,r,L3,T,1,q3,p3,MgO,MgO_th,E,n)
+ return freq_3-freq_1-freq_2
+ T = T0
+ df = Df(T)
+ while math.fabs(df)>1E-6: #1 kHz accuracy
+ #print "T = ",T," df = ", df,
+ df1 = Df(T+Tstep)
+ ddfdt = (df1-df)/Tstep
+ dT = df/ddfdt
+ T -= dT
+ Tstep = math.fabs(dT/2.0)
+ df = Df(T)
+ #print " dT = ", dT, " T -> ",T, " df -> ", df
+ if T<-200 or T> 500: break
+ return T
+
+def T_phasematch_blk(wl1,wl2, T0=70,MgO=0,MgO_th=5.0,E = 0,Tstep=0.1, n = nLNO1):
+ '''
+ Finds 1/wl2 + 1/wl2 = 1/wl3 (wl in nanometers), collinear bulk phase matching temperature
+ based on the initial guess T0 with initial search step Tstep
+ Requires importing Sellmeyer equations.
+ '''
+ import math
+ wl3 = 1/(1.0/wl1+1.0/wl2)
+ def Dk(T):
+ k1 = 2*math.pi*n(wl1,T,0,MgO,MgO_th,E)*1e7/wl1 # cm^-1
+ k2 = 2*math.pi*n(wl2,T,0,MgO,MgO_th,E)*1e7/wl2
+ k3 = 2*math.pi*n(wl3,T,1,MgO,MgO_th,E)*1e7/wl3
+ return k3-k2-k1
+ T = T0
+ dk = Dk(T)
+ while math.fabs(dk)>0.1: #10 cm beat length
+ dk1 = Dk(T+Tstep)
+ ddkdt = (dk1-dk)/Tstep
+ dT = dk/ddkdt
+ T -= dT
+ Tstep = math.fabs(dT/2)
+ dk = Dk(T)
+ #print " dT = ", dT, " T -> ",T, " dk -> ", dk
+ if T<-200 or T> 500:
+ print "No phase matching!"
+ break
+ return T
+
+def SH_phasematch_blk(angle, wl0=1000,step=1.01, n = nBBO):
+ '''
+ Finds the fundamental wavelength (in nm) for frequency doubling o+o->e at a given angle between the optical axis and the e.
+ Requires importing Sellmeyer equations.
+ '''
+ import math
+ def Dk(wl):
+ k1 = 2*math.pi*n(wl,90)*1e7/wl # ordinary cm^-1
+ k2 = 4*math.pi*n(wl/2.0,angle)*1e7/wl # extraordinary with angle
+ return k2-2*k1
+ wl = wl0
+ dk = Dk(wl)
+ while math.fabs(dk)>0.1: #10 cm beat length
+ dk1 = Dk(wl*step)
+ ddkdt = (dk1-dk)/step
+ dwl = dk/ddkdt
+ wl = wl/dwl
+ #step = min(step,math.fabs(dwl/2) )
+ dk = Dk(wl)
+ print " dwl = ", dwl, " wl -> ",wl, " dk -> ", dk
+ return wl
+
+
+def Phasematch(wlP,p1,q1,p2,q2,p3,q3,R,r,T,MgO=0,MgO_th=5.0,E = 0, n = nLNO1):
+ '''
+
+ NOT FULLY FUNCTIONAL!!
+
+ Finds wl1 and wl2 such that 1/wl2 + 1/wl2 = 1/wlP and L1-p1+L2-p2 = L3-p3 at a given temperature T.
+ Requires importing Sellmeyer equations. R is evaluated at T. This program is unaware of the other orbital
+ selection rules! L1, L2, L3 are allowed to be non-integer.
+ '''
+ import math
+ from WGM_lib import Get_frac_L
+ def Dm(wl): # m3-m1-m2, orbital detuning
+ L3 = Get_frac_L(R,r,wlP,T,1,q3,p3,MgO,MgO_th)
+ L1 = Get_frac_L(R,r,wl,T,0,q1,p1,MgO,MgO_th)
+ wl2 = 1.0/(1.0/wlP-1.0/wl)
+ L2 = Get_frac_L(R,r,wl2,T,0,q2,p2,MgO,MgO_th)
+ return L2-L3+p3+L1-p1-p2
+ wl = 2*wlP #try at degeneracy
+ dm = Dm(wl)
+ print "dm = ", dm
+ ddm = 0
+ wlstep = 1.0 # chosing an appropriate wavelength step so the m variation is not too large or too small
+ while math.fabs(ddm)<1:
+ wlstep = 2*wlstep
+ ddm = Dm(wl+wlstep)-dm
+ print "wlstep = ", wlstep, " initial slope = ", ddm/wlstep
+ while math.fabs(dm)>1E-6: # equiv *FSR accuracy
+ dm1 = Dm(wl+wlstep)
+ print "dm1 = ", dm1
+ slope = (dm1-dm)/wlstep
+ print "slope = ", slope, "nm^-1"
+ dwl = dm/slope
+ wl -= dwl
+ print "new wl = ", wl
+ wlstep = min(math.fabs(dwl/2), wlstep)
+ dm = Dm(wl)
+ print " d wl = ", dwl, " wl -> ",wl, " dm -> ", dm
+ if math.fabs(dm)> Get_frac_L(R,r,2*wlP,T,0,q1,p1,MgO,MgO_th): break
+ L1 = Get_frac_L(R,r,wl,T,0,q1,p1,MgO,MgO_th)
+ wl2 = 1.0/(1.0/wlP-1.0/wl)
+ L2 = Get_frac_L(R,r,wl2,T,0,q2,p2,MgO,MgO_th)
+ return L1, L2, wl, wl2
+
+def n_eff_FSR(wl,R,T,FSR,pol,MgO=0,MgO_th=5.0,E = 0,n = nLNO1,n_disp=nLNO1_disp):
+ '''
+ Finds effective refraction index (material+geometrical) based on the FSR measuremert (in GHz).
+ R in cm, wl in nm. Also tries to guess the appropriate q and L. The first order approximation.
+ '''
+ import math
+ from scipy.special import ai_zeros
+ import bisect
+ nr = n(wl,T,pol,MgO,MgO_th,E)
+ wl_deriv = n_disp(wl,T,pol,MgO,MgO_th)
+ denom = 2*math.pi*3*R*FSR*(nr-wl_deriv*wl)/29.9792-2
+ neff = nr/denom
+ AiRoots = -ai_zeros(40)[0]
+ L = 2*math.pi*R*neff*1e7/wl
+ a_q = (nr/neff-1)*(2*L*L)**(1.0/3.0)
+ j = bisect.bisect(AiRoots, a_q)
+ #q = 0
+ #if a_q>AiRoots[0]/2.0: q = 1
+ if j: q = j+(a_q-AiRoots[j-1])/(AiRoots[j]-AiRoots[j-1])
+ else: q = 1+(a_q-AiRoots[0])/(AiRoots[1]-AiRoots[0])
+ return neff,q,L
+
+def n_eff(wl,R,r,T,pol,q,p,MgO=0,MgO_th=5.0,E = 0,n = nLNO1):
+ '''
+ Finds effective refraction index (material+geometrical) for a given WGM.
+ R,r in cm, wl in nm.
+ '''
+ import math
+ from WGM_lib import Get_frac_L
+ L = Get_frac_L(R,r,wl,T,pol,q,p,MgO,MgO_th,E,n)
+ neff = L*wl*1e-7/(2*math.pi*R)
+ return neff,L
+ \ No newline at end of file
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-16 08:50:28 +0000