from pdPythonLib import *
import math
import cmath
fimm = pdApp()
fimm.ConnectToApp("localhost",5101)

# functions

# Note: the Gaussian beam is defined as:
# H = exp(-((Z0-Z)/MFD)**2) [in A/m] where MFD is the mode field diameter [in um]
# E = 376.73031*H [in V/m]

def overlap(Ey,MFD,Z0,Zsteps,dZ):
    integral = 0
    for i in range (1,Zsteps,2):
        integral = integral + dZ/3*(Ey[i-1]*cmath.exp(-((Z0-(i-1)*dZ)/MFD)**2)+4*Ey[i]*cmath.exp(-((Z0-i*dZ)/MFD)**2)+Ey[i+1]*cmath.exp(-((Z0-(i+1)*dZ)/MFD)**2))
    return integral

def normaliseInput(): # from input field
    norm = 0
    fimm.Exec("{FPpath}.calcfieldprofile(0,1,1,0)")
    fimm.Exec("Set field = {FPpath}.fieldprofile.getfieldarray(2)")
    inputEy = fimm.Exec("field.fieldarray")
    fimm.Exec("Set field = {FPpath}.fieldprofile.getfieldarray(4)")
    inputHx = fimm.Exec("field.fieldarray") #given in mA/m
    for i in range (0,int(nx),1):
        inputHx[i][0]=inputHx[i][0]/1000 #now given in A/m
    for i in range (1,nx+1,2):
        norm = norm + W/float(nx)/3*(inputEy[i-1][0]*inputHx[i-1][0]+4*inputEy[i][0]*inputHx[i][0]+inputEy[i+1][0]*inputHx[i+1][0])
    norm = math.sqrt(norm.real**2+norm.imag**2)
    return norm

# main script

# Set the parameters below if you make any changes to the structure
FPpath = "app.subnodes[1].subnodes[7]"      # path of the FIMMPROP Device
Zmax = 24.61                # [initially 24.61] length of FP Device
nx = 200                    # nx for the cross-section
W = 11.22                   # [initially 11.22] width of cross-section
xpos = -0.61                # [initially -0.61] x-position at which the overlap is calculated (x = 0 being the centre of the cross-section)

# The parameter below define the resolution of the field export
Zsteps = 2001               # [MUST BE ODD NUMBER] number of points used to export the fields

# The parameters below handle the position / scanning of the Gaussian beam
Zscan = 0                   # [0 or 1] 0 for using a fixed z-position for the centre of the Gaussian beam (Z0fixed), 1 for scanning the position of the Gaussian beam
Zscansteps = 1000           # [Zscan = 1 only] number of points for the scanning of the centre of the Gaussian beam
Zscanmin = 0                # [Zscan = 1 only] min Z value for centre of Gaussian beam
Zscanmax = Zmax             # [Zscan = 1 only] min Z value for centre of Gaussian beam
Z0fixed = 11.435            # [Zscan = 0 only; initially 11.435] Z value for centre of Gaussian beam

# The parameters below handle the scanning of the wavelength
wavsteps = 3                # number of wavelength points
wavmin = 1.4                # min wavelength
wavmax = 1.7                # max wavelength

# Internal - please do not change
dZ = Zmax/(Zsteps-1)
OVres4=[]
OVres9=[]
Z0res4=[]
Z0res9=[]

# First we will scan the wavelength and for each value of the MFD calculate the overlap with the fibre
# If Zscan is enabled then the position of the Gaussian beam will be scanned and the maximum overlap selected
for i in range (0,wavsteps,1):
    wav = wavmin+i*(wavmax-wavmin)/(wavsteps-1)
    fimm.Exec("{FPpath}.lambda={wav}")
    fimm.Exec("{FPpath}.update")
    fimm.Exec("{FPpath}.calczfield({Zsteps},0,{Zmax},{xpos},0,1)")
    Ey = fimm.Exec("{FPpath}.zfieldcomp.EyTotal") # first term Ey[0] is None
    Ey = Ey[1:] # Ey[0] -> z = 0, Ey[i] -> z = i*Zmax/(Zsteps-1)
    normaliseInp = normaliseInput()
    for MFD in [4,9]: # mode field diameter in microns
        OVnormmax = 0
        if (Zscan==1): # scanning the position of the Gaussian beam
            Z0opt = 0
            for j in range (0,Zscansteps,1):
                Z0 = Zscanmin+j*(Zscanmax-Zscanmin)/(Zscansteps-1)
                OV = overlap(Ey,MFD,Z0,Zsteps,dZ)
                OVmod = math.sqrt(OV.real**2+OV.imag**2)
                if (MFD==4):
                    normaliseGaussian4 = 0
                    for i in range (1,Zsteps,2):
                        normaliseGaussian4 = normaliseGaussian4 + 376.73031*dZ/3*(math.exp(-2*((Z0-(i-1)*dZ)/4)**2)+4*math.exp(-2*((Z0-i*dZ)/4)**2)+math.exp(-2*((Z0-(i+1)*dZ)/4)**2))
                    OVnorm = OVmod/cmath.sqrt(normaliseInp*normaliseGaussian4)
                    if (OVnorm.real**2+OVnorm.imag**2>OVnormmax.real**2+OVnormmax.imag**2):
                        OVnormmax = OVnorm
                        Z0opt = Z0
                elif (MFD==9):
                    normaliseGaussian9 = 0
                    for i in range (1,Zsteps,2):
                        normaliseGaussian9 = normaliseGaussian9 + 376.73031*dZ/3*(math.exp(-2*((Z0-(i-1)*dZ)/9)**2)+4*math.exp(-2*((Z0-i*dZ)/9)**2)+math.exp(-2*((Z0-(i+1)*dZ)/9)**2))
                    OVnorm = OVmod/cmath.sqrt(normaliseInp*normaliseGaussian9)
                    if (OVnorm.real**2+OVnorm.imag**2>OVnormmax.real**2+OVnormmax.imag**2):
                        OVnormmax = OVnorm
                        Z0opt = Z0
            if (MFD==4):
                OVres4.append(OVnormmax.real)
                Z0res4.append(Z0opt)
            elif (MFD==9):
                OVres9.append(OVnormmax.real)
                Z0res9.append(Z0opt)
        elif (Zscan==0): # using a fixed position for the Gaussian beam
            OV = overlap(Ey,MFD,Z0fixed,Zsteps,dZ)
            OVmod = math.sqrt(OV.real**2+OV.imag**2)
            if (MFD==4):
                normaliseGaussian4 = 0
                for i in range (1,Zsteps,2):
                    normaliseGaussian4 = normaliseGaussian4 + 376.73031*dZ/3*(math.exp(-2*((Z0fixed-(i-1)*dZ)/4)**2)+4*math.exp(-2*((Z0fixed-i*dZ)/4)**2)+math.exp(-2*((Z0fixed-(i+1)*dZ)/4)**2))
                OVnorm = OVmod/cmath.sqrt(normaliseInp*normaliseGaussian4)
                OVres4.append(OVnorm.real)
            elif (MFD==9):
                normaliseGaussian9 = 0
                for i in range (1,Zsteps,2):
                    normaliseGaussian9 = normaliseGaussian9 + 376.73031*dZ/3*(math.exp(-2*((Z0fixed-(i-1)*dZ)/9)**2)+4*math.exp(-2*((Z0fixed-i*dZ)/9)**2)+math.exp(-2*((Z0fixed-(i+1)*dZ)/9)**2))
                OVnorm = OVmod/cmath.sqrt(normaliseInp*normaliseGaussian9)
                OVres9.append(OVnorm.real)

# Displaying results
if (Zscan==1):
    print "wavelength(um) coupling_norm(MFD=4) optimumZ0(MFD=4)(um) coupling_norm(MFD=9) optimumZ0(MFD=9)(um)"
    for i in range (0,wavsteps,1):
        wav = wavmin+i*(wavmax-wavmin)/(wavsteps-1)
        print str(wav) + " " + str(OVres4[i]) + " " +str(Z0res4[i]) + " " + str(OVres9[i]) + " " +str(Z0res9[i])
elif (Zscan==0):
    print "wavelength(um) coupling_norm(MFD=4) coupling_norm(MFD=9)"
    for i in range (0,wavsteps,1):
        wav = wavmin+i*(wavmax-wavmin)/(wavsteps-1)
        print str(wav) + " " + str(OVres4[i]) + " " + str(OVres9[i]) 

del fimm
