factor de llenado

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

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

l   = 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)

Npsi  = 2000             # numero de filling fractions
psii  = 0
psif  = 1
dpsi  = (psif-psii)/Npsi
fa    = np.zeros(Npsi+1)
aux   = np.zeros((Npsi+1,Nf+1),dtype=complex)


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

data= open("python.dat","w")
for ipsi in range(0,Npsi+1):
  fa[ipsi] = psii+dpsi*ipsi
  la       = fa[ipsi]*l
  lb       = l-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[ipsi,cf] = cmath.acos(aux4)
    if abs(aux.imag[ipsi,cf]) > 0.1:
#       print(fa[ipsi],abs(aux.imag[ipsi,cf]),f[cf]/1e12)  
       data.write(str(fa[ipsi])+" "+str(abs(aux.imag[ipsi,cf]))+" "+str(f[cf]/1e12)+"\n")
data.close()