下载后，后缀修改去掉.ra即可执行
源代码
#!/usr/bin/env python from math import exp
from gnuplot_leon import *
imp0 = 377.0 class fdtd_leon:
# Author : Leon Email: [email protected]
# fdtd simulation #initialization
def __init__(self,size=400,time=0,MaxTime=1000,delay = 30, width = 10, cdtds =1.0):
self.ez = size * [0.00]
self.hy = size * [0.00]
self.ceze = size * [0.00]
self.chye = size * [0.00]
self.cezh = size * [0.00]
self.chyh = size * [0.00]
self.size = size
self.time = 0
self.MaxTime = MaxTime
self.delay = delay
self.width = width
self.cdtds = cdtds # grid initialization
def grid_init(self,loss = 0.02, loss_layer = 180, epsr = 9.0):
for mm in range(0, self.size):
if (mm < 100):
self.ceze[mm] = 1.0
self.cezh[mm] = imp0
elif (mm < loss_layer):
self.ceze[mm] = 1.0
self.cezh[mm] = imp0/epsr
else:
self.ceze[mm] = (1.0-loss)/(1.0+loss)
self.cezh[mm] = imp0/epsr//(1.0+loss)
if(mm < loss_layer):
self.chyh[mm] = 1.0
self.chye[mm] = 1.0/imp0
else:
self.chyh[mm] = (1.0-loss)/(1.0+loss)
self.chye[mm] = 1.0/imp0/(1.0+loss) # update electric field
def update_e(self):
mm =0;
for mm in range(1,self.size-1):
self.ez[mm] = self.ez[mm]*self.ceze[mm] + (self.hy[mm]-self.hy[mm-1])*self.cezh[mm] # update magnetic field
def update_h(self):
mm =0;
for mm in range(0,self.size-1):
self.hy[mm] = self.hy[mm]*self.chyh[mm] + (self.ez[mm+1]-self.ez[mm])*self.chye[mm] def ez_inc_init(self):
#self.delay = int(raw_input('Enter delay:'))
#self.width = int(raw_input('Enter width:'))
#self.cdtds = int(raw_input('Enter cdtds:'))
#print self.delay
#print self.width
#print self.cdtds
return def ez_inc(self,time,location):
#print ''.join(['ez_inc time: ',str(time),'location: ',str(location)] )
#print exp(-((time-self.delay-location/self.cdtds)/self.width)**2)
return exp(-((time-self.delay-location/self.cdtds)/self.width)**2) def abc_init(self):
return def abc(self):
self.ez[0] = self.ez[1] def tfsf_init(self):
tfsf_boundary = raw_input('Enter location of tfsf boundary:')
self.ez_inc_init()
return int(tfsf_boundary) def tfsf_update(self,tfsf_boundary):
#print self.time
#print tfsf_boundary
tfsf_boundary = int(tfsf_boundary)
if (tfsf_boundary <= 0):
print 'tfsf boundary error \n'
return
else:
self.hy[tfsf_boundary] -= self.ez_inc(self.time,0.0)*self.chye[tfsf_boundary]
self.ez[tfsf_boundary+1] += self.ez_inc(self.time+0.5,-0.5)
例子
#!/usr/bin/env python import sys
import math
import os
from gnuplot_leon import *
from fdtd_leon import *
import threading
# Author : Leon Email: [email protected]
# fdtd simulation , plotting with gnuplot, writting in python
# python and gnuplot software packages should be installed before running this program
# 1d fdtd with absorbing boundary and TFSF boundary
# lossy dielectric material def snashot(gp, fdtd,interval):
"""Record a frame data into the gif file Parameters
----------
gp : class gnuplot_leon
fdtd : class fdtd_leon
interval : int record one every [interval] frames
"""
if(fdtd.time % interval == 0):
gp.set_frame_start('l', 1, 'green')
cnt = 0
for elem in fdtd.ez:
gp.update_point(cnt,elem)
cnt += 1
gp.set_frame_end() def report():
"""report the rate of progress Parameters
----------
none
"""
global MaxTime
global qTime
print ''.join([str(int(1000.00*int(qTime+1)/int(MaxTime))/10.0),'% has been finished!'])
if(qTime>=MaxTime-1):
return
global timer
timer = threading.Timer(2.0,report)
timer.start() gp = gnuplot_leon()
gp.set_plot_size(0.85,0.85)
gp.set_canvas_size(600,400)
#gp.set_title('fdtd simulation by leon : gnuplot class test')
title = 'fdtd simulation by leon,yangli0534\\\\@gmail.com' gp.set_title(title)
gp.set_gif()
#gp.set_png()
gp.set_file_name('demo3.gif')
gp.set_tics_color('white')
gp.set_border_color('orange')
gp.set_grid_color('orange')
gp.set_bkgr_color('gray10')
gp.set_xlabel('length','white')
gp.set_ylabel('amplitude','white')
gp.auto_scale_enable()
gp.set_key('off','sin(x)','white') size = 400#physical distance
#ez=size * [0.00]#electric field
#hy=size * [0.00]#magnetic field
#ceze=size * [0.00]#
#cezh=size * [0.00]#
#chye=size * [0.00]#
#chyh=size * [0.00]#
#sinwave=size * [0.00]#
imp0 = 377.00
LOSS = 0.01
LOSS_LAYER = 250
qTime = 0
MaxTime = 18000
delay = 30
width = 10
cdtds =1.0
epsR = 9.0
tfsf_boundary = 0
interval = 30
#cnt = 0
#elem = 0.00000
gp.set_x_range(0,size-1) fdtd = fdtd_leon(size,0,MaxTime,delay,width,cdtds)
fdtd.grid_init(LOSS, LOSS_LAYER, epsR)
fdtd.abc_init()
tfsf_boundary = fdtd.tfsf_init()
timer = threading.Timer(1,report)
timer.start()
# do time stepping
for fdtd.time in range(0, MaxTime):
qTime = fdtd.time
fdtd.update_h()
fdtd.tfsf_update(tfsf_boundary)
fdtd.abc()
fdtd.update_e()
snashot(gp,fdtd,interval) gp.set_output_valid()
gp.close()