1 files changed, 71 insertions, 0 deletions
diff --git a/phasematching/PDC_phasematchT2.py b/phasematching/PDC_phasematchT2.py
new file mode 100644
index 0000000..4bbb60d
--- /dev/null
+++ b/phasematching/PDC_phasematchT2.py
@@ -0,0 +1,71 @@
+#!/usr/bin/env mpython_q
+"""
+This program finds various {qs, q1, q2} PDC channels near a given pump wavelength and for the a range of signal wavelemgths.
+Equatorial modes only. Version 2 also allows for nonzero p1 p2 p3 but then the angular overlap is not accurate.
+"""
+
+from WGM_lib import Get_frac_L, nLNO1, WGM_freq, T_phasematch,RadOvlp
+import math
+import os
+import pylab
+import bisect
+
+#os.chdir("F:\Users docs\Dima\Projects\WGM switch\Numeric simulations")
+os.chdir("./")
+
+Trange = [10,180] # C
+R0 = 0.15 # cm, disk radius at 20 C
+r = R0/4 # cm, rim radius
+wls = 532.0 # nm
+wlRange1 = [1000, 1100] # signal min wavelength in nm, max wl in nm
+MgO = 5.7 # %  5.83% new wafer in Erlangen; 5.54% JPL wafer by SH, 5.02% JPL wafer by 1560+780
+MgO_th = 5.0 # %
+dL = 5 # signal FSR step
+qlist = [[1,1,1],[1,1,2],[2,1,1]]
+p1max, p2max, dp = 2,2,1 # dp = max(ps-p1-p2); angular overlap rapidly decreases with dp 
+
+p1range = range(0, p1max+1) # must start from an even number!
+p2range = range(0, p2max+1) # must start from an even number!
+dprange = range(0,dp/2+1)
+def R(T): return R0*(1+1.7e-5*(T-20)) # expantion coefficient of Lithium Niobate
+T0 = 0.5*(Trange[0]+Trange[1])
+           
+f = open("PhaseMatchEqv"+str(10*R0)+"mm_doped"+str(MgO)+"_"+str(wls)+"nm2.dat", "w")
+f.write("T wl1 wl2 L1 L2 Ls p1 p2 ps q1 q2 qs rad_ovlp ang_ovlp ovlp_sq \n" )
+for qq in qlist:
+    q1,q2,qs = qq[0],qq[1],qq[2]
+    T = T0
+    for p1 in p1range:
+        for p2 in p2range:
+            psrange = [] # ps must not exceed p1+p2 and have the same parity - a selection rule!
+            for x in dprange:
+                ps = p1+p2-2*x
+                if ps < 0: break
+                psrange.append(ps)
+            for ps in psrange:
+                L1 = int(math.floor(Get_frac_L(R(Trange[0]),r,wlRange1[1],Trange[0],0,q1,p1,MgO,MgO_th))) # starting (the minimal) L1
+                wl1 = wlRange1[1]
+                while wl1 > wlRange1[0]:
+                    wl1, freq1 = WGM_freq(R(T),r,L1,T,0,q1,p1,MgO,MgO_th)
+                    Ls = int(Get_frac_L(R(T),r,wls,T,1,qs,ps,MgO,MgO_th))
+                    L2 = Ls-L1+p1+p2-ps # orbital phase matching for m's !
+                    Tcorr = T_phasematch(L1,p1,q1,L2,p2,q2,Ls,ps,qs,R(T),r,T,MgO,MgO_th)
+                    if Tcorr >= Trange[0] and Tcorr <= Trange[1]:
+                        T = Tcorr
+                        wl1, freq1 = WGM_freq(R(T),r,L1,T,0,q1,p1,MgO,MgO_th) # find exact wavelengths
+                        Ls = int(Get_frac_L(R(T),r,wls,T,1,qs,ps,MgO,MgO_th))
+                        L2 = Ls-L1+p1+p2-ps
+                        T = T_phasematch(L1,p1,q1,L2,p2,q2,Ls,ps,qs,R(T),r,T,MgO,MgO_th) # second iteration
+                        wl2 = 1.0/(1.0/wls-1.0/wl1)
+                        # wl20, freq2 = WGM_freq(R(T),r,L2,T,0,q2,0,MgO,MgO_th) # test
+                        ovlpr = (RadOvlp(R(T),L1,q1,L2,q2,Ls,qs))**2 #only radial overlap-square in cm^-3
+                        if p1+p2+ps: # the formula below is only correct for al equatorial modes; -1 for others.
+                            ovlpa = -1
+                        else:
+                            ovlpa = math.pi**(-1.5)*math.sqrt((L1*L2/float(L1+L2))) # angular overlap-square
+                        f.write("%.12g %.12g %.12g %.12g %.12g %.12g %.12g %.12g %.12g %.12g %.12g %.12g %.12g %.12g %.12g \n" \
+                        % (T, wl1, wl2, L1, L2, Ls, p1, p2, ps, q1, q2, qs, ovlpr, ovlpa, ovlpa*ovlpr))
+                        # print 2.99792E8/wls-freq1-freq2 # test: this should not exceed half fsr
+                    L1 += dL                
+f.close()
+print "\n Done!"