import math
import cmath
im = complex(0.0,1.0)
print(im)
ka = 0.235
ca = 3600
ra = 1500
ta = 1
CCa = math.sqrt(ka/(ra*ca*ta))
la = 0.01e-3
kb = 0.445
cb = 3300
rb = 1116
tb = 20
CCb = math.sqrt(kb/(rb*cb*tb))
lb = 0.01e-3
fi = 0
ff = 10
Nf = 100
df = (ff-fi)/Nf
data= open("gap.dat","w")
for cf in range(1,Nf):
f = fi+cf*df
w = 2.0*math.pi*f
ga = cmath.sqrt((w*w+im*(w/ta))/(CCa*CCa))
gb = cmath.sqrt((w*w+im*(w/tb))/(CCb*CCb))
na = ka/(1-im*w*ta)
nb = kb/(1-im*w*tb)
aux1 = cmath.cosh(im*ga*la)*cmath.cosh(im*gb*lb)
aux2 = ( (na*ga)/(nb*gb) + (nb*ga)/(na*ga) )
aux3 = cmath.sinh(im*ga*la)*cmath.sinh(im*gb*lb)
aux4 = aux1+0.5*aux2*aux3
aux5 = cmath.acos(aux4)
data.write(str(f)+" "+str(aux5.real)+" "+str(abs(aux5.imag)) + "\n")
data.close()
import numpy as np
import matplotlib.pyplot as plt
data = np.loadtxt('gap.dat')
x = data[:, 0]
y = data[:, 1]
z = data[:, 2]
plt.plot(y, x,'r--',z,x,'g--')
plt.show()
