#!/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 = [ 60,140] # C
R0 = 0.35 # cm, disk radius at 20 C
r = R0/4 # cm, rim radius
wls = 532.0 # nm
wlRange1 = [1063.5, 1064.5] # signal min wavelength in nm, max wl in nm
MgO = 6.2 # %  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!"