periodo

import numpy as np
import math
import cmath
import matplotlib.pyplot as plt
im  = complex(0.0,1.0)

######################  -> parametros materiales

l0   = 2.0e-9             # periodo CU

# Si
ka  = 150.0             
ca  = 710.0
ra  = 2330.0
ta  = 150.0e-12
CCa = math.sqrt(ka/(ra*ca*ta))

#Ge
kb  = 60.0
cb  = 310.0
rb  = 5300.0
tb  = 200.0e-12
CCb = math.sqrt(kb/(rb*cb*tb))

######################  -> parametros para el ciclo

Nf    = 1000              # Numero de frecuencias
fi    = 0
ff    = 800.0e9
df    = (ff-fi)/Nf
f     = np.zeros(Nf+1)

Nl    = 100             # numero de filling fractions
li    = 1*l0
lf    = 2*l0
dl    = (lf-li)/Nl
l     = np.zeros(Nl+1)
aux   = np.zeros((Nl+1,Nf+1),dtype=complex)

fa    = 0.5
  
for cf in range(1,Nf+1):
  f[cf]    = fi+cf*df

data= open("python.dat","w")
for il in range(0,Nl+1):
  l[il]    = li + dl*il
  la       = fa*l[il]
  lb       = l[il]-la
  for cf in range(1,Nf+1):
    w    = 2.0*math.pi*f[cf]
    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*gb)/(na*ga) )
    aux3 = cmath.sinh(im*ga*la)*cmath.sinh(im*gb*lb)
    aux4 = aux1+0.5*aux2*aux3
    aux[il,cf] = cmath.acos(aux4)
    if abs(aux.imag[il,cf]) > 0.1:
       print(l[il]/1e-9,abs(aux.imag[il,cf]),f[cf]/1e12)  
       data.write(str(l[il]/1e-9)+" "+str(abs(aux.imag[il,cf]))+" "+str(f[cf]/1e12)+"\n")
data.close()