8478阅读
33回复

[RFEDA原创]Python的实现

UID ：63797

注册：2010-07-23
登录：2010-10-19
发帖：28
等级：仿真二级

10楼发表于: 2010-07-28 13:13:34

fd1d_2.1.c是fd1d_1.5.c的另一种实现方式，利用麦克斯韦中的散度方程引入了dx和ix：

dx[k] = dx[k]+0.5*(hy[k-1]-hy[k])
ex[k] = ga[k]*(dx[k]-ix[k])
ix[k] = ix[k]+gb[x]*ex[k]
hy[k] = hy[k]+0.5*(ex[k]-ex[k+1])

式中：

ga[k] = 1/(epsilon+sigma*dt/epsz)
gb[k] = sigma*dt/epsz

这次是把程序中的循环用向量化的方式实现的：

from pylab import *

KE = 200; NSTEPS = 500

ex = zeros(KE); hy = zeros(KE)
dx = zeros(KE); ix = zeros(KE)

epsz = 8.8e-12
ddx = 0.01; dt = ddx/6e8

kstart = 100
epsilon = 4.0; sigma = 0.04

ga = ones(KE); ga[kstart:] = 1.0/(epsilon+sigma*dt/epsz)
gb = zeros(KE); gb[kstart:] = sigma*dt/epsz

ex_low_m1 = 0.0; ex_high_m1 = 0.0
ex_low_m2 = 0.0; ex_high_m2 = 0.0

freq_in = 7e8; T = 0

for n in range(NSTEPS+1):

    T += 1

    dx[1:] = dx[1:] + 0.5*(hy[:-1]-hy[1:])
    pulse = sin(2*pi*freq_in*dt*T); dx[5] = dx[5]+pulse
    ex = ga*(dx-ix); ix = ix+gb*ex

    ex[0] = ex_low_m2; ex_low_m2 = ex_low_m1;   ex_low_m1 = ex[1]
    ex[-1] = ex_high_m2; ex_high_m2 = ex_high_m1; ex_high_m1 = ex[-2]

    hy[:-1] = hy[:-1] + 0.5*(ex[:-1]-ex[1:])

plot(ex)
show()

离线samuelp

UID ：63797

注册：2010-07-23
登录：2010-10-19
发帖：28
等级：仿真二级

11楼发表于: 2010-07-28 16:44:57

fd1d_2.2.c是计算500MHz处的频域响应

from pylab import *

KE = 200; NSTEPS = 400

ex = zeros(KE); hy = zeros(KE)
dx = zeros(KE); ix = zeros(KE)

epsz = 8.8e-12
ddx = 0.01; dt = ddx/6e8

kstart = 100
epsilon = 4.0; sigma = 0.0

ga =  ones(KE); ga[kstart:] = 1.0/(epsilon+sigma*dt/epsz)
gb = zeros(KE); gb[kstart:] = sigma*dt/epsz

ex_low_m1 = 0.0; ex_high_m1 = 0.0
ex_low_m2 = 0.0; ex_high_m2 = 0.0

freq_in = 7e8; T = 0

real_pt = zeros((KE,3)); imag_pt = zeros((KE,3))
ampn = zeros((KE,3)); phasen = zeros((KE,3))

real_in = zeros(3); imag_in = zeros(3)
amp_in = zeros(3); phase_in = zeros(3)

freq = array([1.0e8,2.0e8,5.0e8]); arg = 2*pi*freq*dt

t0 = 50.0; spread = 10.0

ex3 = zeros((KE,3))

for n in range(NSTEPS+1):

    T += 1

    dx[1:] = dx[1:] + 0.5*(hy[:-1]-hy[1:])
    pulse = exp(-0.5*((t0-T)/spread)**2); dx[5] = dx[5]+pulse
    ex = ga*(dx-ix); ix = ix+gb*ex

    ex[0]  =  ex_low_m2;  ex_low_m2 = ex_low_m1;   ex_low_m1 = ex[1]
    ex[-1] = ex_high_m2; ex_high_m2 = ex_high_m1; ex_high_m1 = ex[-2]

    for m in range(3): ex3[:,m] = ex
    real_pt = real_pt+cos(arg*T)*ex3
    imag_pt = imag_pt-sin(arg*T)*ex3

    real_in = real_in+cos(arg*T)*ex[5]
    imag_in = imag_in-sin(arg*T)*ex[5]

    hy[:-1] = hy[:-1] + 0.5*(ex[:-1]-ex[1:])

    amp_in = sqrt(real_in**2+imag_in**2)
    ampn = 1.0/amp_in*sqrt(real_pt**2+imag_pt**2)

plot(ampn[:,2],label='T='+str(NSTEPS))
legend()
show()

图片:fd1d_2.2.JPG

离线samuelp

UID ：63797

注册：2010-07-23
登录：2010-10-19
发帖：28
等级：仿真二级

12楼发表于: 2010-07-28 17:43:45

如果是色散介质，增加参数：

gc = chil*dt/tau

fd1d_2.3.c：

from pylab import *

KE = 200; NSTEPS = 1000

ex = zeros(KE); hy = zeros(KE)
dx = zeros(KE); ix = zeros(KE)
sx = zeros(KE)

epsz = 8.8e-12
ddx = 0.01; dt = ddx/6e8

kstart = 100
epsilon = 2.0; sigma = 0.01; chil = 2.0; tau = 1.0e-9; del_exp = exp(-dt/tau)

ga = ones(KE); ga[kstart:] = 1.0/(epsilon+sigma*dt/epsz+chil*dt/tau)
gb = zeros(KE); gb[kstart:] = sigma*dt/epsz
gc = zeros(KE); gc[kstart:] = chil*dt/tau

ex_low_m1 = 0.0; ex_high_m1 = 0.0
ex_low_m2 = 0.0; ex_high_m2 = 0.0

freq_in = 7e8; T = 0

real_pt = zeros((KE,3)); imag_pt = zeros((KE,3))
ampn = zeros((KE,3)); phasen = zeros((KE,3))

real_in = zeros(3); imag_in = zeros(3)
amp_in = zeros(3); phase_in = zeros(3)

freq = array([50.0e6, 200.0e6, 500.0e6]); arg = 2*pi*freq*dt

t0 = 50.0; spread = 10.0

ex3 = zeros((KE,3))

for n in range(NSTEPS+1):

    T += 1

    dx[1:] = dx[1:] + 0.5*(hy[:-1]-hy[1:])
    pulse = exp(-0.5*((t0-T)/spread)**2); dx[5] = dx[5]+pulse
    ex = ga*(dx-ix-sx); ix = ix+gb*ex; sx = del_exp*sx+gc*ex

    ex[0] = ex_low_m2; ex_low_m2 = ex_low_m1;   ex_low_m1 = ex[1]
    ex[-1] = ex_high_m2; ex_high_m2 = ex_high_m1; ex_high_m1 = ex[-2]

    for m in range(3): ex3[:,m] = ex
    real_pt = real_pt+cos(arg*T)*ex3
    imag_pt = imag_pt-sin(arg*T)*ex3

    real_in = real_in+cos(arg*T)*ex[5]
    imag_in = imag_in-sin(arg*T)*ex[5]

    hy[:-1] = hy[:-1] + 0.5*(ex[:-1]-ex[1:])

    amp_in = sqrt(real_in**2+imag_in**2)
    ampn = 1.0/amp_in*sqrt(real_pt**2+imag_pt**2)

for m in range(3):
    plot(ampn[:,m],label='Freq: '+str(freq[m]/1e6)+'MHz')
legend()
show()

第二章完活。

图片:fd1d_2.3.JPG

离线haiyang

UID ：55297

注册：2010-03-21
登录：2023-07-03
发帖：246
等级：仿真三级

13楼发表于: 2010-07-28 23:49:35

很有用来分享

离线samuelp

UID ：63797

注册：2010-07-23
登录：2010-10-19
发帖：28
等级：仿真二级

14楼发表于: 2010-07-29 12:57:41

二维FDTD，fd2d_3.1.c：

from pylab import *
from mpl_toolkits.mplot3d import Axes3D

IE = 60;   JE = 60;      NSTEPS = 50
ic = IE/2; jc = JE/2
ddx = 0.01; dt = ddx/6.0e8; epsz = 8.8e-12

dz = zeros((IE,JE)); ez = zeros((IE,JE))
hx = zeros((IE,JE)); hy = zeros((IE,JE))
ga = ones((IE,JE))

t0 = 20.0; spread = 6.0; T = 0

for n in range(NSTEPS):
    T += 1
    for i in range(1,IE):
        for j in range(1,JE):
            dz[j] = dz[j]+0.5*(hy[j]-hy[i-1][j]-hx[j]+hx[j-1]) .A|@?p[
    pulse = exp(-0.5*((t0-T)/spread)**2); dz[ic][ic] = pulse HE\K@3-
    ez = ga*dz wKY_Bo/d
    for i in range(IE-1): _','9|
        for j in range(JE-1): 4 H&#q>
            hx[j] = hx[j]+0.5*(ez[j]-ez[j+1]) 2>59q$|
            hy[j] = hy[j]+0.5*(ez[i+1][j]-ez[j]) O33`+UV"W
fig = plt.figure() -ze J#B)C
ax = Axes3D(fig) f,Ghb~y
x = arange(IE); y = arange(JE) O&hTNIfi
X,Y = meshgrid(x,y) Gp\ kU:}&
ax.plot_surface(X,Y,ez,rstride=1, cstride=1, cmap=cm.jet) Kf-JcBsrT
show() IvNT6]6 P

图片:fd2d_3.1.JPG

离线caocheng82

UID ：10116

注册：2008-03-28
登录：2025-05-26
发帖：697
等级：积极交流六级

15楼发表于: 2010-07-29 15:00:31

这个图是python画的么？

感觉很不错，与MATLAB相比还不错

不知道能否写动画。

离线caocheng82

UID ：10116

注册：2008-03-28
登录：2025-05-26
发帖：697
等级：积极交流六级

16楼发表于: 2010-07-29 15:02:05

还有就是，不知道楼主对scipy，numpy等python的数值包用得怎样？这些需要配置lapack等fortran包，当时配置了好长时间以失败而告终。

离线samuelp

UID ：63797

注册：2010-07-23
登录：2010-10-19
发帖：28
等级：仿真二级

17楼发表于: 2010-07-29 16:00:18

是用绘图库matplotlib画的，也可用mayavi，都是针对科学绘图的。

安装pythonxy或者enthought就什么都不用再配置了呀，都是针对科学计算做的软件包，numpy等全都在里面，pythonxy是免费的，并且包括了enthought。enthought虽然是商业软件，但是可以免费使用。出了pythonxy后，应该很少有人自己配置numpy，scipy，gcc等软件包了。还有一款叫sage的，适合做计算平台，不适合做软件开发。

动画：
01    from pylab import  *
02    import os
03    # set the interactive mode of pylab to ON
04    ion()
05    # opens a new figure to plot into
06    fig_hndl = figure()
07    # make an empty list into which we'll append
08    # the filenames of the PNGs that compose each
09    # frame.
10    files=[]
11    # filename for the name of the resulting movie
12    filename = 'animation'
13    number_of_frames = 100
14    for t in range(number_of_frames):
15        # draw the frame
16        imshow(rand(100,100))
17        # form a filename
18        fname = '_tmp%03d.png'%t
19        # save the frame
20        savefig(fname)
21        # append the filename to the list
22        files.append(fname)
23    # call mencoder
24    os.system("mencoder 'mf://_tmp*.png' -mf type=png:fps=10
25        -ovc lavc -lavcopts vcodec=wmv2 -oac copy
26        -o " + filename + ".mpg")
27    # cleanup
28    for fname in files: os.remove(fname)

图片:fd2d_3.1.JPG

离线caocheng82

UID ：10116

注册：2008-03-28
登录：2025-05-26
发帖：697
等级：积极交流六级

18楼发表于: 2010-07-29 16:26:01

用过很多开源的东西，比如scilab，gnuplot，maxima等，由于很多都是从LINUX泊来的，因此在WINDOWS下的配置令人绝望。
希望这个东西不错吧。

还有能否介绍一下python和C++的联合编程么？想用python做前后处理，fdtd计算用C++。各发挥各自的优点。两者怎样接口，能否给大家培训一下。谢谢！

离线samuelp

UID ：63797

注册：2010-07-23
登录：2010-10-19
发帖：28
等级：仿真二级

19楼发表于: 2010-07-29 17:23:46

我在windows下装的，spyder相当于matlab，eclipse pydev相当于vc ide，完全傻瓜式的，不需要计算机背景知识。混合编程可参考《Python Scripting for Computational Science》的相应章节。
ithov.com上有下载

发帖回复