FDTD Python API

可执行文件地址：

下载后，后缀修改去掉.ra即可执行

源代码

#!/usr/bin/env python

from math import exp
from gnuplot_leon import *
imp0 = 377.0

class fdtd_leon:
# Author : Leon  Email: yangli0534@gmail.com
# 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: yangli0534@gmail.com
# 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()