Hello,
I have a program using Numeric and weave developed with Python2.3. I
just switched to Python2.4 and numpy. The Numeric/weave version is
almost a factor of 2 faster compared to numpy/weave. Is that what is to
be expected or are there options to improve the speed of numpy/weave? I
would be very appreciative for any help. Please find the source attached.
Regards
Rolf
--
------------------------------------
# Dr. Rolf Wester
# Fraunhofer Institut f. Lasertechnik
# Steinbachstrasse 15, D-52074 Aachen, Germany.
# Tel: + 49 (0) 241 8906 401, Fax: +49 (0) 241 8906 121
# EMail: [EMAIL PROTECTED]
# WWW: http://www.ilt.fraunhofer.de
import math
import sys
if sys.version[0:3] == '2.3':
from Numeric import *
import weave
else:
from numpy import *
import scipy.weave as weave
Float64 = float64
Complex64 = complex128
cc = 2.997924e14 # mum/s
mu0 = 1.256637e-12 # Vs/(A mum)
eps0 = 1.0/(mu0*cc**2) # As/(V mum)
Z0 = sqrt(mu0/eps0)
pi = math.pi
def make_spaces(n):
s = ""
for i in range(n):
s = s + ' '
return s
def split2(str, sep = ' \t\n\r', trimp = 'yes'):
vstr = []
n0 = len(str)
n = len(sep)
ipos0 = 0
ipos1 = 0
while ipos1 < n0 and ipos0 < n0:
ipos0 = ipos1
flag = 0
c2 = str[ipos0]
for i in range(n):
if(c2 == sep[i]):
flag = 1
break
while(flag == 1):
ipos0 = ipos0+1
if(ipos0 > n0-1):
break
flag = 0
c2 = str[ipos0]
for i in range(n):
if(c2 == sep[i]):
flag = 1
break
if(ipos0 > n0-1):
break
ipos1 = ipos0
flag = 0
c2 = str[ipos1]
for i in range(n):
if c2 == sep[i]:
flag = 1
break
while flag == 0:
ipos1 = ipos1+1
if ipos1 > n0-1:
break
flag = 0
c2 = str[ipos1]
for i in range(n):
if c2 == sep[i]:
flag = 1
break
if(ipos0 < ipos1):
teil = str[ipos0:ipos1]
if(trimp == 'yes'):
teil = string.strip(teil)
vstr.append(teil)
return vstr
class LSSOLV:
def __init__(self, N):
self.N = N
self.gamma = zeros(N, Complex64)
self.betai = zeros(N, Complex64)
def inline_solve(self, a, b, c, d, y):
betai = self.betai
gamma = self.gamma
N = self.N
code = """
int i;
betai[0] = 1.0/b[0];
for(i=1; i < N; i++)
betai[i] = 1.0/(b[i]-a[i]*c[i-1]*betai[i-1]);
gamma[0] = d[0]*betai[0];
for(i=1; i < N; i++)
gamma[i] = (d[i]-a[i]*gamma[i-1])*betai[i];
y[N-1] = gamma[N-1];
for(i=N-2; i >= 0; i--)
y[i] = gamma[i]-c[i]*y[i+1]*betai[i];
return_val = y[300];
"""
err = weave.inline(code,
['a', 'b', 'c', 'd', 'betai', 'gamma', 'y', 'N'],
#type_converters = weave.converters.blitz,
compiler = 'gcc',verbose=2)
return err
def blitz_solve(self, a, b, c, d, y):
betai = self.betai
gamma = self.gamma
N = self.N
code = """
int i;
betai(0) = 1.0/b(0);
for(i=1; i < N; i++)
betai(i) = 1.0/(b(i)-a(i)*c(i-1)*betai(i-1));
gamma(0) = d(0)*betai(0);
for(i=1; i < N; i++)
gamma(i) = (d(i)-a(i)*gamma(i-1))*betai(i);
y(N-1) = gamma(N-1);
for(i=N-2; i >= 0; i--)
y(i) = gamma(i)-c(i)*y(i+1)*betai(i);
return_val = y(300);
"""
err = weave.inline(code,
['a', 'b', 'c', 'd', 'betai', 'gamma', 'y', 'N'],
type_converters = weave.converters.blitz,
compiler = 'gcc')
return err
def solve(self, a, b, c, d, y):
betai = self.betai
gamma = self.gamma
N = self.N
betai[0] = 1.0/b[0]
for i in range(1,N):
betai[i] = 1.0/(b[i]-a[i]*c[i-1]*betai[i-1])
gamma[0] = d[0]*betai[0]
for i in range(1,N):
gamma[i] = (d[i]-a[i]*gamma[i-1])*betai[i]
y[self.N-1] = gamma[N-1]
for i in range(N-2,-1,-1):
y[i] = gamma[i]-c[i]*y[i+1]*betai[i]
class Parameter:
def __init__(self, filename):
self.map = {}
ein = file(filename, 'r')
inhalt = ein.readlines()
head = ''
for line in inhalt:
line = split2(line,sep = '=#', trimp = 'yes')
if len(line) > 0 and len(line[0]) > 0 and line[0][0] == '[' and line[0][len(line[0])-1] == ']':
head = line[0][1:len(line[0])-1]
if len(line) >= 2:
self.map[head + '.' + line[0]] = line[1]
def get(self, key, val):
ret = val
if self.map.has_key(key):
if type(val) == type(1.0):
ret = float(self.map[key])
elif type(val) == type(1):
ret = int(self.map[key])
else:
ret = self.map[key]
return ret
def set(self, key, val):
if type(val) == type(1.0):
self.map[key] = "%e" % val
elif type(val) == type(1):
self.map[key] = "%i" % val
else:
self.map[key] = val
def print_parameter(self, out = None):
list = []
for key in self.map:
list.append([key, self.map[key]])
list.sort()
lmax = 0
for key in self.map:
n = len(key)
if (n > lmax):
lmax = n
ausgabe = "---- parameter begin ----\n"
for l in list:
n = len(l[0])
spaces = make_spaces(lmax+2-n)
ausgabe = ausgabe + l[0] + spaces + " = " + l[1] + '\n'
ausgabe = ausgabe + "---- parameter end ----\n"
if out == None:
print ausgabe
else:
out.write("%s", ausgabe)
class Crystal:
def __init__(self):
self.d_eff = 34.1e-6 #microm/V uebliche Angabe in pm/V -> 10^-12 m/V -> 10^-9 mm/V -> 10^-6 microm/V
def get_deff(self):
return self.d_eff
def get_no(self, lamb):
#lamb = lamb * 1.0e3
l2 = lamb*lamb
n2 = self.Ae*l2/(l2-self.Be) + self.Ce*l2/(l2-self.De) + self.Ee*l2/(l2-self.Fe)
if(n2+1 >= 0.0):
return sqrt(n2+1)
return 1.0
def get_ne(self, lamb):
#lamb = lamb * 1.0e3
l2 = lamb*lamb
n2 = self.Ao*l2/(l2-self.Bo) + self.Co*l2/(l2-self.Do) + self.Eo*l2/(l2-self.Fo)
if(n2+1 >= 0.0):
return sqrt(n2+1)
return 1.0
def get_alpha(self, lamb):
if (lamb <= 1.5 ): alpha = 0.005
elif (lamb >= 1.5 and lamb < 2.0): alpha = 0.005
elif (lamb >= 2.0 and lamb < 2.5): alpha = 0.009
elif (lamb >= 2.5 and lamb < 3.0): alpha = 0.012
elif (lamb >= 3.0 and lamb < 3.5): alpha = 0.018
elif (lamb >= 3.5 and lamb < 4.0): alpha = 0.03
elif (lamb >= 4.0 and lamb < 4.5): alpha = 0.08
elif (lamb >= 4.5 and lamb < 5.0): alpha = 0.12
elif (lamb >= 5.0): alpha = 0.12
return alpha*1.0e-3
class LiNbO3(Crystal):
def __init__(self):
self.Ao = 2.6734;
self.Bo = 0.01764;
self.Co = 1.2290;
self.Do = 0.05914;
self.Eo = 12.614;
self.Fo = 474.6;
self.Ae = 2.9804;
self.Be = 0.02047;
self.Ce = 0.5981;
self.De = 0.0666;
self.Ee = 8.9543;
self.Fe = 416.08;
Crystal.__init__(self)
class LiNbO3Mg(Crystal):
def __init__(self):
self.Ae = 2.4372; self.Ao = 2.2454;
self.Be = 0.01478; self.Bo = 0.01242;
self.Ce = 1.4617; self.Co = 1.3005;
self.De = 0.05612; self.Do = 0.05313;
self.Ee = 9.6536; self.Eo = 6.8972;
self.Fe = 416.08; self.Fo = 331.33;
Crystal.__init__(self)
class NLO:
def __init__(self, par, log):
self.k1 = zeros(3, Complex64)
self.k2 = zeros(3, Complex64)
self.k3 = zeros(3, Complex64)
self.k4 = zeros(3, Complex64)
self.y0 = zeros(3, Complex64)
self.y = zeros(3, Complex64)
self.log = log
self.par = par
self.lambda_s = self.par.get("beam.lambda_s", -1.0)
self.lambda_i = self.par.get("beam.lambda_i", -1.0)
self.lambda_p = self.par.get("beam.lambda_p", -1.0)
if(self.lambda_s < 0.0): self.lambda_s = 1.0/(1.0/self.lambda_p - 1.0/self.lambda_i)
if(self.lambda_i < 0.0): self.lambda_i = 1.0/(1.0/self.lambda_p - 1.0/self.lambda_s)
if(self.lambda_p < 0.0): self.lambda_p = 1.0/(1.0/self.lambda_s + 1.0/self.lambda_i)
self.par.set("beam.lambda_s", self.lambda_s)
self.par.set("beam.lambda_i", self.lambda_i)
self.par.set("beam.lambda_p", self.lambda_p)
cr_type = self.par.get("crystal.type", "LiNbO3Mg")
cr_oe = self.par.get("crystal.oe" , "o")
if(cr_type == "LiNbO3"):
cr = LiNbO3()
else:
cr = LiNbO3Mg()
self.n_s = self.par.get("crystal.n_s", -1.5)
self.n_i = self.par.get("crystal.n_i", -1.5)
self.n_p = self.par.get("crystal.n_p", -1.5)
if(self.n_s < 0.0):
if(cr_oe == "e"):
self.n_s = cr.get_ne(self.lambda_s)
else:
self.n_s = cr.get_no(self.lambda_s)
if(self.n_i < 0.0):
if(cr_oe == "e"):
self.n_i = cr.get_ne(self.lambda_i)
else:
self.n_i = cr.get_no(self.lambda_i)
if(self.n_p < 0.0):
if(cr_oe == "e"):
self.n_p = cr.get_ne(self.lambda_p)
else:
self.n_p = cr.get_no(self.lambda_p)
self.par.set("crystal.n_s", self.n_s)
self.par.set("crystal.n_p", self.n_p)
self.par.set("crystal.n_i", self.n_i)
self.nq_s = self.par.get("waveguide.nq_s", -1.0)
self.nq_i = self.par.get("waveguide.nq_i", -1.0)
self.nq_p = self.par.get("waveguide.nq_p", -1.0)
if self.nq_s < 0.0: self.nq_s = self.n_s
if self.nq_i < 0.0: self.nq_i = self.n_i
if self.nq_p < 0.0: self.nq_p = self.n_p
self.par.set("waveguide.nq_s", self.nq_s)
self.par.set("waveguide.nq_i", self.nq_i)
self.par.set("waveguide.nq_p", self.nq_p)
self.alpha_s = self.par.get("crystal.alpha_s", -1.5) # 1/mm
self.alpha_i = self.par.get("crystal.alpha_i", -1.5) # 1/mm
self.alpha_p = self.par.get("crystal.alpha_p", -1.5) # 1/mm
if(self.alpha_s < 0.0): self.alpha_s = cr.get_alpha(self.lambda_s)
if(self.alpha_i < 0.0): self.alpha_i = cr.get_alpha(self.lambda_i)
if(self.alpha_p < 0.0): self.alpha_p = cr.get_alpha(self.lambda_p)
self.par.set("crystal.alpha_s", self.alpha_s);
self.par.set("crystal.alpha_i", self.alpha_i);
self.par.set("crystal.alpha_p", self.alpha_p);
self.alpha_s = self.alpha_s * 0.5
self.alpha_i = self.alpha_i * 0.5
self.alpha_p = self.alpha_p * 0.5
self.k_s = 2.0*pi/self.lambda_s*self.n_s;
self.k_i = 2.0*pi/self.lambda_i*self.n_i;
self.k_p = 2.0*pi/self.lambda_p*self.n_p;
self.par.set("waveguide.k_s", self.k_s)
self.par.set("waveguide.k_i", self.k_i)
self.par.set("waveguide.k_p", self.k_p)
self.kq_s = 2.0*pi/self.lambda_s*self.nq_s;
self.kq_i = 2.0*pi/self.lambda_i*self.nq_i;
self.kq_p = 2.0*pi/self.lambda_p*self.nq_p;
self.par.set("waveguide.kq_s", self.kq_s)
self.par.set("waveguide.kq_i", self.kq_i)
self.par.set("waveguide.kq_p", self.kq_p)
self.omega_s = 2.0*pi * cc/self.lambda_s;
self.omega_i = 2.0*pi * cc/self.lambda_i;
self.omega_p = 2.0*pi * cc/self.lambda_p;
self.dk0 = self.kq_p - self.kq_i - self.kq_s;
if(abs(self.dk0) < 1.0e-7):
self.dk = 0.0;
else:
self.dk = self.dk0;
self.par.set("beam.dk" , self.dk); #1/mm
self.deff = self.par.get("crystal.deff", -1.0)
if(self.deff < 0.0): self.deff = cr.get_deff()
self.par.set("crystal.deff", self.deff)
self.cs = self.deff * self.omega_s / (cc * self.nq_s);
self.ci = self.deff * self.omega_i / (cc * self.nq_i);
self.cp = self.deff * self.omega_p / (cc * self.nq_p);
self.par.set("crystal.cs", self.cs);
self.par.set("crystal.ci", self.ci);
self.par.set("crystal.cp", self.cp);
self.z_coh = par.get("crystal.z_coh" , -1.0)
if(self.z_coh < 0.0):
if(self.dk != 0.0):
self.z_coh = abs(pi/self.dk)
else:
self.z_coh = 1.0e20
self.par.set("crystal.z_coh", self.z_coh)
self.dz = self.par.get("par.dz", -1.0)
self.nsteps = self.par.get("par.nsteps", 0)
self.naus = self.par.get("par.naus", 0)
if(self.dz < 0.0):
self.dz = self.z_coh/nsteps
self.par.set("par.dz", self.dz)
self.length = par.get("crystal.length", -1.0)
h1 = -1.0e3
if(self.length > 0.0):
h1 = self.length/(self.dz*self.nsteps)
self.naus = int( h1+0.1 )
par.set("par.nsteps", self.nsteps)
par.set("par.naus", self.naus)
self.par.print_parameter()
self.z = 0.0
def get_z_coh(self):
return self.z_coh
def get_dk(self):
return dk
def get_cs(self):
return self.cs
def get_ci(self):
return self.ci
def get_cp(self):
return self.cp
def inline_step(self, Es, Ei, Ep, ix1, ix2):
k1 = self.k1
k2 = self.k2
k3 = self.k3
k4 = self.k4
y = self.y
y0 = self.y0
alpha_s = self.alpha_s
alpha_i = self.alpha_i
alpha_p = self.alpha_p
cs = self.cs
ci = self.ci
cp = self.cp
dz = self.dz
dkz = self.dk * dz
a1 = int(floor(self.z/self.z_coh))
sign_c = -1.0;
if int(floor(float(a1/2))) * 2 == a1:
sign_c = 1.0
code = """
double h = dz;
double hd3 = dz/3.0;
double hd6 = dz/6.0;
double hh = 0.5*dz;
std::complex<double> im(0.0,1.0);
std::complex<double> cs_ = im * sign_c * cs;
std::complex<double> ci_ = im * sign_c * ci;
std::complex<double> cp_ = im * sign_c * cp;
std::complex<double> exp_parg = exp( im * dkz);
std::complex<double> exp_marg = exp(- im * dkz);
double as = alpha_s*0.5;
double ai = alpha_i*0.5;
double ap = alpha_p*0.5;
int i,ix;
int n = 3;
for(ix = ix1; ix <= ix2; ix++)
{
y0[0] = Es[ix];
y0[1] = Ei[ix];
y0[2] = Ep[ix];
k1[0] = cs_ * y0[2] * conj(y0[1]) * exp_marg - as*y0[0];
k1[1] = ci_ * y0[2] * conj(y0[0]) * exp_marg - ai*y0[1];
k1[2] = cp_ * y0[0] * y0[1] * exp_parg - ap*y0[2];
for(i=0; i < n; i++) y[i] = y0[i] + hh * k1[i];
k2[0] = cs_ * y[2] * conj(y[1]) * exp_marg - as*y[0];
k2[1] = ci_ * y[2] * conj(y[0]) * exp_marg - ai*y[1];
k2[2] = cp_ * y[0] * y[1] * exp_parg - ap*y[2];
for(i=0; i < n; i++) y[i] = y0[i] + hh * k2[i];
k3[0] = cs_ * y[2] * conj(y[1]) * exp_marg - as*y[0];
k3[1] = ci_ * y[2] * conj(y[0]) * exp_marg - ai*y[1];
k3[2] = cp_ * y[0] * y[1] * exp_parg - ap*y[2];
for(i=0; i < n; i++) y[i] = y0[i] + h * k3[i];
k4[0] = cs_ * y[2] * conj(y[1]) * exp_marg - as*y[0];
k4[1] = ci_ * y[2] * conj(y[0]) * exp_marg - ai*y[1];
k4[2] = cp_ * y[0] * y[1] * exp_parg - ap*y[2];
Es[ix] = y0[0] + hd6 * (k1[0] + k4[0]) + hd3 * (k2[0] + k3[0]);
Ei[ix] = y0[1] + hd6 * (k1[1] + k4[1]) + hd3 * (k2[1] + k3[1]);
Ep[ix] = y0[2] + hd6 * (k1[2] + k4[2]) + hd3 * (k2[2] + k3[2]);
}
"""
weave.inline(code,
['k1', 'k2', 'k3', 'k4', 'y', 'y0', 'Es', 'Ei', 'Ep', 'dz', 'ix1', 'ix2', 'dkz', 'cs', 'ci', 'cp', 'alpha_s', 'alpha_i', 'alpha_p', 'sign_c'],
#type_converters = weave.converters.blitz,
compiler = 'gcc')
self.z = self.z + self.dz
def step(self, Es, Ei, Ep, ix1, ix2, dkz):
k1 = self.k1
k2 = self.k2
k3 = self.k3
k4 = self.k4
y = self.y
y0 = self.y0
alpha_s = self.alpha_s
alpha_i = self.alpha_i
alpha_p = self.alpha_p
cs = self.cs
ci = self.ci
cp = self.cp
dz = self.dz
h = dz
hd3 = dz/3.0
hd6 = dz/6.0
hh = 0.5*dz
im = 0.0+1.0j
sign_c = 1.0 + 0.0j
cs_ = im * sign_c * cs
ci_ = im * sign_c * ci
cp_ = im * sign_c * cp
exp_parg = exp( im * dkz)
exp_marg = exp(- im * dkz)
as = alpha_s*0.5
ai = alpha_i*0.5
ap = alpha_p*0.5
n = 3
#print Es[300], Ei[300], Ep[300], exp_parg, cs_, ci_, cp_
for ix in range(ix1,ix2+1):
y0[0] = Es[ix]
y0[1] = Ei[ix]
y0[2] = Ep[ix]
k1[0] = cs_ * y0[2] * conjugate(y0[1]) * exp_marg - as*y0[0]
k1[1] = ci_ * y0[2] * conjugate(y0[0]) * exp_marg - ai*y0[1]
k1[2] = cp_ * y0[0] * y0[1] * exp_parg - ap*y0[2]
for i in range(n):
y[i] = y0[i] + hh * k1[i]
k2[0] = cs_ * y[2] * conjugate(y[1]) * exp_marg - as*y[0]
k2[1] = ci_ * y[2] * conjugate(y[0]) * exp_marg - ai*y[1]
k2[2] = cp_ * y[0] * y[1] * exp_parg - ap*y[2]
for i in range(n):
y[i] = y0[i] + hh * k2[i]
k3[0] = cs_ * y[2] * conjugate(y[1]) * exp_marg - as*y[0]
k3[1] = ci_ * y[2] * conjugate(y[0]) * exp_marg - ai*y[1]
k3[2] = cp_ * y[0] * y[1] * exp_parg - ap*y[2]
for i in range(n):
y[i] = y0[i] + h * k3[i]
k4[0] = cs_ * y[2] * conjugate(y[1]) * exp_marg - as*y[0]
k4[1] = ci_ * y[2] * conjugate(y[0]) * exp_marg - ai*y[1]
k4[2] = cp_ * y[0] * y[1] * exp_parg - ap*y[2]
Es[ix] = y0[0] + hd6 * (k1[0] + k4[0]) + hd3 * (k2[0] + k3[0])
Ei[ix] = y0[1] + hd6 * (k1[1] + k4[1]) + hd3 * (k2[1] + k3[1])
Ep[ix] = y0[2] + hd6 * (k1[2] + k4[2]) + hd3 * (k2[2] + k3[2])
self.z = self.z + self.dz
class BPM1D:
def __init__(self, par, wave, log):
self.par = par
self.log = log
self.dz = self.par.get("par.dz", -1.0)
self.nsteps = self.par.get("par.nsteps", 0)
self.naus = self.par.get("par.naus", 0)
self.w_core = self.par.get("waveguide.w_core", -1.0)
self.w_cladding = self.par.get("waveguide.w_cladding", -1.0)
self.nx = self.par.get("waveguide.nx", 1)
self.alpha = self.par.get("waveguide.alpha", -1.0)
self.n_cladding = self.par.get("waveguide.n_cladding", -1.0)
self.nq = self.par.get("waveguide.nq_" + wave, -1.0)
self.n_core = self.par.get("crystal.n_" + wave, -1.5)
self.lamb = self.par.get("beam.lambda_" + wave, -1.0)
self.P = self.par.get("beam.P_" + wave, -1.0)
self.w0 = self.par.get("beam.w0_" + wave, -1.0)
self.shape = self.par.get("beam.shape_" + wave, "gauss")
self.omega = 2.0*math.pi/self.lamb*cc
self.dx = (self.w_core + 2.0*self.w_cladding) / (self.nx-1)
nkk = int(floor(self.w_cladding / self.dx));
self.E = zeros(self.nx, Complex64)
self.n2 = zeros(self.nx, Complex64)
for i in range(self.nx):
self.n2[i] = self.n_cladding*self.n_cladding
for i in range(nkk,self.nx-nkk):
self.n2[i] = self.n_core*self.n_core
self.ix1 = nkk
self.ix2 = self.nx-nkk-1
self.a = zeros(self.nx, Complex64)
self.b = zeros(self.nx, Complex64)
self.c = zeros(self.nx, Complex64)
self.d = zeros(self.nx, Complex64)
self.k0 = self.k(1.0) + 0.0j
self.k02 = self.k(1.0)*self.k(1.0) + 0.0j
self.nq2 = self.nq*self.nq + 0.0j
self.beta = self.k0*self.nq
## m_im = complex<double>(0.0, 1.0);
## m_cc0 = complex<double>(0.0, 0.0);
## m_cc1 = complex<double>(1.0, 0.0);
self.ddx2 = 1.0/self.dx**2
self.x = zeros(self.nx, Float64)
self.xmin = -(self.w_cladding+0.5*self.w_core)
for i in range(self.nx):
self.x[i] = self.xmin + i * self.dx;
# init E
for i in range(self.nx):
self.E[i] = 0.0 + 0.0j
if self.shape == "gauss":
for i in range(self.nx):
arg = self.x[i]*self.x[i]/(self.w0*self.w0)
if arg < 700.0:
self.E[i] = exp(- arg) + 0.0j
elif self.shape == "sin":
for i in range(nkk,self.nx-nkk):
a = sin((i-nkk)*self.dx/self.w_core*pi)
self.E[i] = a + 0.0j
else:
for i in range(self.nx):
if abs(self.x[i]) <= self.w0:
self.E[i] = 1.0 + 0.0j
P = self.get_P()
factor = sqrt(self.P/P)
for i in range(self.nx):
self.E[i] = self.E[i] * factor
self.z = 0.0
self.lssolv = LSSOLV(self.nx)
def k(self, n):
return 2.0*pi / self.lamb * n
def get_kq(self):
return self.nq*self.k0
def set_E(self, E):
for i in range(self.nx):
self.E[i] = E[i]
def get_E(self, E):
for i in range(self.nx):
E[i] = self.E[i]
def get_I(self):
dZq = 0.5*self.nq/Z0
I = dZq*((self.E*conjugate(self.E)).real)
return I
def get_phase(self):
phase = zeros(self.nx, Float64)
for i in range(self.nx):
x = self.E[i].real
y = self.E[i].imag
if x > 0.0 and y >= 0.0:
phase[i] = math.atan(y/x)
elif x < 0.0 and y >= 0.0:
phase[i] = pi - math.atan(-y/x)
elif x < 0.0 and y <= 0.0:
phase[i] = pi + math.atan(y/x)
elif x > 0.0 and y <= 0.0:
phase[i] = 2.0*pi - math.atan(-y/x)
elif x == 0.0 and y >= 0.0:
phase[i] = pi*0.5
elif x == 0.0 and y < 0.0:
phase[i] = pi*1.5
return phase
def get_P(self):
I = self.get_I()
P = 0.0
for i in range(self.nx):
P = P + I[i]
return P
def compute_res(self):
i = self.nx/2
ddx2 = self.ddx2
k02 = self.k02
n2 = self.n2
nq2 = self.nq2
a = ddx2*(self.E[i-1] + self.E[i+1]) - 2.0*ddx2 * self.E[i]
c = pi**2/5.0**2*self.E[i]
b = k02*(n2[i] - nq2)*self.E[i]
return a,c,(c-b)/(2.0*self.beta),(a-b)/(2.0*self.beta)
def blitz_set_matrix(self):
im = 0.0 + 1.0j
#(1,1) Pade
# ? -> siehe Randbedingungen
c1 = 1.0/(4.0 * self.k02*self.nq2) + im * self.dz* self.alpha / (2.0*self.k0*self.nq)
c2 = 1.0/(4.0 * self.k02*self.nq2) - im * self.dz*(1.0-self.alpha) / (2.0*self.k0*self.nq)
#paraxial approximation
# ? -> siehe Randbedingungen
#c1 = + im * self.dz* self.alpha / (2.0*self.k0*self.nq)
#c2 = - im * self.dz*(1.0-self.alpha) / (2.0*self.k0*self.nq)
a = self.a
b = self.b
c = self.c
d = self.d
E = self.E
ddx2 = self.ddx2
k02 = self.k02
n2 = self.n2
nq2 = self.nq2
nx = self.nx
code = """
int i;
for(i = 0; i < nx; i++)
{
a(i) = c1*ddx2;
c(i) = c1*ddx2;
b(i) = 1.0 - 2.0*c1*ddx2 + c1*k02*(n2(i) - nq2);
}
for(i = 1; i < nx-1; i++)
{
d(i) = E(i) * (1.0 - 2.0*c2*ddx2 + c2*k02*(n2(i) - nq2)) + c2*ddx2*(E(i-1) + E(i+1));
}
std::complex<double> gl(1.0,0.0);
std::complex<double> gr(1.0,0.0);
if (abs(E(1)) > 1.0e-20)
gl = E(0)/E(1);
if (abs(E(nx-2)) > 1.0e-20)
gr = E(nx-1)/E(nx-2);
std::complex<double> kr_g = log(gl);
if(kr_g.imag() > 0.0)
{
kr_g = conj(kr_g);
gl = exp(kr_g);
}
kr_g = log(gr);
if(kr_g.imag() > 0.0)
{
kr_g = conj(kr_g);
gr = exp(kr_g);
}
a(0) = std::complex<double>(0.0, 0.0);
b(0) = b(0) + c1*ddx2*gl;
d(0) = E(0) * (1.0 - 2.0*c2*ddx2 + c2*k02*(n2(0) - nq2)) + c2*ddx2*(gl*E(0) + E(1));
b(nx-1) = b(nx-1) + c1*ddx2*gr;
c(nx-1) = std::complex<double>(0.0, 0.0);
d(nx-1) = E(nx-1) * (1.0 - 2.0*c2*ddx2 + c2*k02*(n2(nx-1) - nq2)) + c2*ddx2*(gr*E(nx-1) + E(nx-2));
return_val = E(nx/2);
"""
err = weave.inline(code,
['a', 'b', 'c', 'd', 'E', 'ddx2', 'k02', 'n2', 'nq2', 'nx', 'c1', 'c2'],
type_converters = weave.converters.blitz,
compiler = 'gcc')
return err
def inline_set_matrix(self):
im = 0.0 + 1.0j
#(1,1) Pade
# ? -> siehe Randbedingungen
c1 = 1.0/(4.0 * self.k02*self.nq2) + im * self.dz* self.alpha / (2.0*self.k0*self.nq)
c2 = 1.0/(4.0 * self.k02*self.nq2) - im * self.dz*(1.0-self.alpha) / (2.0*self.k0*self.nq)
#paraxial approximation
# ? -> siehe Randbedingungen
#c1 = + im * self.dz* self.alpha / (2.0*self.k0*self.nq)
#c2 = - im * self.dz*(1.0-self.alpha) / (2.0*self.k0*self.nq)
a = self.a
b = self.b
c = self.c
d = self.d
E = self.E
ddx2 = self.ddx2
k02 = self.k02
n2 = self.n2
nq2 = self.nq2
nx = self.nx
code = """
int i;
for(i = 0; i < nx; i++)
{
a[i] = c1*ddx2;
c[i] = c1*ddx2;
b[i] = 1.0 - 2.0*c1*ddx2 + c1*k02*(n2[i] - nq2);
}
for(i = 1; i < nx-1; i++)
{
d[i] = E[i] * (1.0 - 2.0*c2*ddx2 + c2*k02*(n2[i] - nq2)) + c2*ddx2*(E[i-1] + E[i+1]);
}
std::complex<double> gl(1.0,0.0);
std::complex<double> gr(1.0,0.0);
if (abs(E[1]) > 1.0e-20)
gl = E[0]/E[1];
if (abs(E[nx-2]) > 1.0e-20)
gr = E[nx-1]/E[nx-2];
std::complex<double> kr_g = log(gl);
//wenn die Vorzeichen bei ? umgekehrt werden -> < statt >
if(kr_g.imag() > 0.0)
{
kr_g = conj(kr_g);
gl = exp(kr_g);
}
kr_g = log(gr);
//wenn die Vorzeichen bei ? umgekehrt werden -> < statt >
if(kr_g.imag() > 0.0)
{
kr_g = conj(kr_g);
gr = exp(kr_g);
}
a[0] = std::complex<double>(0.0, 0.0);
b[0] = b[0] + c1*ddx2*gl;
d[0] = E[0] * (1.0 - 2.0*c2*ddx2 + c2*k02*(n2[0] - nq2)) + c2*ddx2*(gl*E[0] + E[1]);
b[nx-1] = b[nx-1] + c1*ddx2*gr;
c[nx-1] = std::complex<double>(0.0, 0.0);
d[nx-1] = E[nx-1] * (1.0 - 2.0*c2*ddx2 + c2*k02*(n2[nx-1] - nq2)) + c2*ddx2*(gr*E[nx-1] + E[nx-2]);
return_val = E[nx/2];
"""
err = weave.inline(code,
['a', 'b', 'c', 'd', 'E', 'ddx2', 'k02', 'n2', 'nq2', 'nx', 'c1', 'c2'],
#type_converters = weave.converters.blitz,
compiler = 'gcc')
return err
def set_matrix(self):
im = 0.0 + 1.0j
#(1,1) Pade
# ? -> siehe Randbedingungen
c1 = 1.0/(4.0 * self.k02*self.nq2) + im * self.dz* self.alpha / (2.0*self.k0*self.nq)
c2 = 1.0/(4.0 * self.k02*self.nq2) - im * self.dz*(1.0-self.alpha) / (2.0*self.k0*self.nq)
#paraxial approximation
# ? -> siehe Randbedingungen
#c1 = + im * self.dz* self.alpha / (2.0*self.k0*self.nq)
#c2 = - im * self.dz*(1.0-self.alpha) / (2.0*self.k0*self.nq)
for i in range(self.nx):
self.a[i] = c1*self.ddx2
self.c[i] = c1*self.ddx2
self.b[i] = 1.0 - 2.0*c1*self.ddx2 + c1* self.k02*(self.n2[i] - self.nq2)
for i in range(1,self.nx-1):
self.d[i] = self.E[i] * (1.0 - 2.0*c2*self.ddx2 + c2* self.k02*(self.n2[i] - self.nq2) ) \
+ c2*self.ddx2*(self.E[i-1] + self.E[i+1])
if (abs(self.E[1]) > 1.0e-20):
gl = self.E[0]/self.E[1]
else:
gl = 1.0+0.0j
if (abs(self.E[self.nx-2]) > 1.0e-20):
gr = self.E[self.nx-1]/self.E[self.nx-2]
else:
gr = 1.0+0.0j
kr_g = log(gl)
if(kr_g.imag > 0.0):
kr_g = conjugate(kr_g)
gl = exp(kr_g)
kr_g = log(gr)
if(kr_g.imag > 0.0):
kr_g = conjugate(kr_g)
gr = exp(kr_g)
self.a[0] = 0.0 + 0.0j
self.b[0] = self.b[0] + c1*self.ddx2*gl
self.d[0] = self.E[0] * (1.0 - 2.0*c2*self.ddx2 + c2* self.k02*(self.n2[0] - self.nq2)) + \
c2*self.ddx2*(gl*self.E[0] + self.E[1])
self.b[self.nx-1] = self.b[self.nx-1] + c1*self.ddx2*gr
self.c[self.nx-1] = 0.0 + 0.0j
self.d[self.nx-1] = self.E[self.nx-1] * (1.0 - 2.0*c2*self.ddx2 + c2* self.k02*(self.n2[self.nx-1] - self.nq2)) + \
c2*self.ddx2*(gr*self.E[self.nx-1] + self.E[self.nx-2])
def solve(self):
self.lssolv.solve(self.a, self.b, self.c, self.d, self.E)
self.z = self.z + self.dz
def blitz_solve(self):
self.lssolv.blitz_solve(self.a, self.b, self.c, self.d, self.E)
self.z = self.z + self.dz
def inline_solve(self):
self.lssolv.inline_solve(self.a, self.b, self.c, self.d, self.E)
self.z = self.z + self.dz
class WG_OPO:
def __init__(self, parameterfile):
self.par = Parameter(parameterfile)
logfilename = self.par.get("log.logfilename", "wg_opo.log")
self.log = file(logfilename,'w')
self.nlo = NLO(self.par, self.log)
self.bpm_s = BPM1D(self.par, "s", self.log)
self.bpm_i = BPM1D(self.par, "i", self.log)
self.bpm_p = BPM1D(self.par, "p", self.log)
print self.bpm_s.get_I()[300], self.bpm_s.get_I()[350]
def run(self):
Ps = []
Pi = []
Pp = []
z = []
print self.bpm_s.get_I()[300], self.bpm_s.get_I()[350]
self.naus = self.par.get("par.naus", 100)
self.nsteps = self.par.get("par.nsteps", 1000)
for i in range(self.naus):
for j in range(self.nsteps):
self.bpm_s.inline_set_matrix()
self.bpm_s.inline_solve()
self.bpm_i.inline_set_matrix()
self.bpm_i.inline_solve()
self.bpm_p.inline_set_matrix()
self.bpm_p.inline_solve()
self.nlo.inline_step(self.bpm_s.E, self.bpm_i.E, self.bpm_p.E, self.bpm_p.ix1, self.bpm_p.ix2)
Ps.append(self.bpm_s.get_P())
Pi.append(self.bpm_i.get_P())
Pp.append(self.bpm_p.get_P())
z.append(self.bpm_s.z)
iz = len(z)-1
phase_s = self.bpm_s.get_phase()[self.bpm_s.nx/2]
phase_i = self.bpm_i.get_phase()[self.bpm_s.nx/2]
phase_p = self.bpm_p.get_phase()[self.bpm_s.nx/2]
print "%4i %10.6e %10.6e %10.6e %10.6e %10.6e %10.6e" % (i, z[iz], Ps[iz], Pi[iz], Pp[iz], phase_s, phase_p)
out = file("x_I.dat",'w')
x = self.bpm_s.x
nx = self.bpm_s.nx
Is = self.bpm_s.get_I()
Ii = self.bpm_s.get_I()
Ip = self.bpm_s.get_I()
for i in range(nx):
out.write("%10.6e %10.6e %10.6e %10.6e\n" % (x[i],Is[i],Ii[i],Ip[i]))
out.close()
out = file("z_P.dat",'w')
for i in range(len(z)):
out.write("%10.6e %10.6e %10.6e %10.6e\n" % (z[i],Ps[i],Pi[i],Pp[i]))
out.close()
wg_opo = WG_OPO("parameter.par")
wg_opo.run()
-------------------------------------------------------------------------
Take Surveys. Earn Cash. Influence the Future of IT
Join SourceForge.net's Techsay panel and you'll get the chance to share your
opinions on IT & business topics through brief surveys -- and earn cash
http://www.techsay.com/default.php?page=join.php&p=sourceforge&CID=DEVDEV
_______________________________________________
Numpy-discussion mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/numpy-discussion
