Corrections from Dmitry Strekalov.

Fixed angular and radial overlap calculations.

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