下载后，后缀修改去掉.ra即可执行
源代码
1 #!/usr/bin/env python 2 3 from math import exp 4 from gnuplot_leon import * 5 imp0 = 377.0 6 7 class fdtd_leon: 8 # Author : Leon Email: yangli0534@gmail.com 9 # fdtd simulation 10 11 #initialization 12 def __init__(self,size=400,time=0,MaxTime=1000,delay = 30, width = 10, cdtds =1.0): 13 self.ez = size * [0.00] 14 self.hy = size * [0.00] 15 self.ceze = size * [0.00] 16 self.chye = size * [0.00] 17 self.cezh = size * [0.00] 18 self.chyh = size * [0.00] 19 self.size = size 20 self.time = 0 21 self.MaxTime = MaxTime 22 self.delay = delay 23 self.width = width 24 self.cdtds = cdtds 25 26 # grid initialization 27 def grid_init(self,loss = 0.02, loss_layer = 180, epsr = 9.0): 28 for mm in range(0, self.size): 29 if (mm < 100): 30 self.ceze[mm] = 1.0 31 self.cezh[mm] = imp0 32 elif (mm < loss_layer): 33 self.ceze[mm] = 1.0 34 self.cezh[mm] = imp0/epsr 35 else: 36 self.ceze[mm] = (1.0-loss)/(1.0+loss) 37 self.cezh[mm] = imp0/epsr//(1.0+loss) 38 if(mm < loss_layer): 39 self.chyh[mm] = 1.0 40 self.chye[mm] = 1.0/imp0 41 else: 42 self.chyh[mm] = (1.0-loss)/(1.0+loss) 43 self.chye[mm] = 1.0/imp0/(1.0+loss) 44 45 # update electric field 46 def update_e(self): 47 mm =0; 48 for mm in range(1,self.size-1): 49 self.ez[mm] = self.ez[mm]*self.ceze[mm] + (self.hy[mm]-self.hy[mm-1])*self.cezh[mm] 50 51 # update magnetic field 52 def update_h(self): 53 mm =0; 54 for mm in range(0,self.size-1): 55 self.hy[mm] = self.hy[mm]*self.chyh[mm] + (self.ez[mm+1]-self.ez[mm])*self.chye[mm] 56 57 def ez_inc_init(self): 58 #self.delay = int(raw_input('Enter delay:')) 59 #self.width = int(raw_input('Enter ')) 60 #self.cdtds = int(raw_input('Enter cdtds:')) 61 #print self.delay 62 #print self.width 63 #print self.cdtds 64 return 65 66 67 def ez_inc(self,time,location): 68 #print ''.join(['ez_inc time: ',str(time),'location: ',str(location)] ) 69 #print exp(-((time-self.delay-location/self.cdtds)/self.width)**2) 70 return exp(-((time-self.delay-location/self.cdtds)/self.width)**2) 71 72 def abc_init(self): 73 return 74 75 def abc(self): 76 self.ez[0] = self.ez[1] 77 78 def tfsf_init(self): 79 tfsf_boundary = raw_input('Enter location of tfsf boundary:') 80 self.ez_inc_init() 81 return int(tfsf_boundary) 82 83 def tfsf_update(self,tfsf_boundary): 84 #print self.time 85 #print tfsf_boundary 86 tfsf_boundary = int(tfsf_boundary) 87 if (tfsf_boundary <= 0): 88 print 'tfsf boundary error ' 89 return 90 else: 91 self.hy[tfsf_boundary] -= self.ez_inc(self.time,0.0)*self.chye[tfsf_boundary] 92 self.ez[tfsf_boundary+1] += self.ez_inc(self.time+0.5,-0.5)
例子
1 #!/usr/bin/env python 2 3 import sys 4 import math 5 import os 6 from gnuplot_leon import * 7 from fdtd_leon import * 8 import threading 9 # Author : Leon Email: yangli0534@gmail.com 10 # fdtd simulation , plotting with gnuplot, writting in python 11 # python and gnuplot software packages should be installed before running this program 12 # 1d fdtd with absorbing boundary and TFSF boundary 13 # lossy dielectric material 14 15 def snashot(gp, fdtd,interval): 16 """Record a frame data into the gif file 17 18 Parameters 19 ---------- 20 gp : class gnuplot_leon 21 fdtd : class fdtd_leon 22 interval : int record one every [interval] frames 23 """ 24 if(fdtd.time % interval == 0): 25 gp.set_frame_start('l', 1, 'green') 26 cnt = 0 27 for elem in fdtd.ez: 28 gp.update_point(cnt,elem) 29 cnt += 1 30 gp.set_frame_end() 31 32 def report(): 33 """report the rate of progress 34 35 Parameters 36 ---------- 37 none 38 """ 39 global MaxTime 40 global qTime 41 print ''.join([str(int(1000.00*int(qTime+1)/int(MaxTime))/10.0),'% has been finished!']) 42 if(qTime>=MaxTime-1): 43 return 44 global timer 45 timer = threading.Timer(2.0,report) 46 timer.start() 47 48 gp = gnuplot_leon() 49 gp.set_plot_size(0.85,0.85) 50 gp.set_canvas_size(600,400) 51 #gp.set_title('fdtd simulation by leon : gnuplot class test') 52 title = 'fdtd simulation by leon,yangli0534\\@gmail.com' 53 54 gp.set_title(title) 55 gp.set_gif() 56 #gp.set_png() 57 gp.set_file_name('demo3.gif') 58 gp.set_tics_color('white') 59 gp.set_border_color('orange') 60 gp.set_grid_color('orange') 61 gp.set_bkgr_color('gray10') 62 gp.set_xlabel('length','white') 63 gp.set_ylabel('amplitude','white') 64 gp.auto_scale_enable() 65 gp.set_key('off','sin(x)','white') 66 67 size = 400#physical distance 68 #ez=size * [0.00]#electric field 69 #hy=size * [0.00]#magnetic field 70 #ceze=size * [0.00]# 71 #cezh=size * [0.00]# 72 #chye=size * [0.00]# 73 #chyh=size * [0.00]# 74 #sinwave=size * [0.00]# 75 imp0 = 377.00 76 LOSS = 0.01 77 LOSS_LAYER = 250 78 qTime = 0 79 MaxTime = 18000 80 delay = 30 81 width = 10 82 cdtds =1.0 83 epsR = 9.0 84 tfsf_boundary = 0 85 interval = 30 86 #cnt = 0 87 #elem = 0.00000 88 gp.set_x_range(0,size-1) 89 90 fdtd = fdtd_leon(size,0,MaxTime,delay,width,cdtds) 91 fdtd.grid_init(LOSS, LOSS_LAYER, epsR) 92 fdtd.abc_init() 93 tfsf_boundary = fdtd.tfsf_init() 94 timer = threading.Timer(1,report) 95 timer.start() 96 # do time stepping 97 for fdtd.time in range(0, MaxTime): 98 qTime = fdtd.time 99 fdtd.update_h() 100 fdtd.tfsf_update(tfsf_boundary) 101 fdtd.abc() 102 fdtd.update_e() 103 snashot(gp,fdtd,interval) 104 105 106 107 gp.set_output_valid() 108 gp.close()
