fdtd.cpp

应用FDTD方法计算微带传输结构的C＋＋程序

/*
* John B. Schneider
* [e-mail]
*
* Copyright (C) 2003  John B. Schneider
* 
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation (FSF) version 2
* of the License.
* 
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
* 
* The license under which this software is publish is available from
* www.fsf.org/copyleft/gpl.html or www.fsf.org/copyleft/gpl.txt.
*
* This code was provided as the solution to a homework problem I
* assigned in EE 417/517 at Washington State University in the Spring
* of 2003.  The goal was to write a program which would duplicate, at
* least to a large extend, the patch antenna work described by Sheen
* et al., IEEE Trans. MTT, 38(7):849--857, 1990.  HOWEVER, unlike
* Sheen et al., a uniform spatial step is used, thus the geometry
* modeled by this program is not identical to the geometry described
* by Sheen et al. (but it is close).  Modifying the code to realize a
* non-uniform grid should not be difficult (and is left as an
* exercise for the reader).  Furthermore, a second-order Higdon
* absorbing boundary conditions is used on four of the boundaries of
* the computational domain while a first-order one is used on the
* source-plane wall.  (Sheen et al. used a first-order ABC on all of
* the five planes.)  This code has in place the arrays to use a
* second-order ABC on the source-plane, but the actual ABC update
* equations are only first-order since they were found to cause fewer
* artifacts when switching between having the source-wall generate
* fields and absorb fields.
*
* This code can be used to model a microstrip line or a microstrip
* patch antenna (the particular problem being modeled is determined
* at compile-time via various declarations).  To compile this program
* to model a patch antenna and have each of the ABC turned on at each
* face, use a command such as this:
*
*  cc -DPATCH -DABC1 -DABC2 -DABC3 -DABC4 -DABC5 -O sheen_patch.c -o sheen_patch -lm
* 
* The switchs control things as follows:
*
*   PATCH: If present, the patch is modeled, otherwise a microstrip which
*          spans the computational domain.
*   ABC1:  ABC at the left wall, x=0.
*   ABC2:  ABC at the far wall, y=LIMY-1.
*   ABC3:  ABC at the right wall, x=LIMX-1.
*   ABC4:  ABC at the near wall, y=0.
*   ABC5:  ABC at the top wall, z=LIMZ-1.
*
* These switches are rather cumbersome but the 417 students in the
* class (i.e., the undergraduates) didn't have to implement ABC's and
* the switches allowed me to use one program for both solutions
* without having to deal with a bunch of if statements.  Note that
* the 417 students had to take a "snapshot" of the field at time step
* 300 but I have removed that portion of the code.
*
* Some of the parameters of this code are:
*    - Courant number of 1/sqrt(3.1)
*    - Computational domain size: 90 x 130 x 20 cells (in x, y, and z
*      directions, respectively)
*    - del_x = del_y = del_z = 0.265 mm
*    - Source plane at y=0
*    - Ground plane at z=0
*    - Duroid substrate with relative permittivity 2.2.  Electric
*      field nodes on interface between duroid and freespace use
*      average permittivity of media to either side.
*    - Substrate 3 cells thick
*    - Microstrip 9 cells wide
*    - Patch dimensions 47 x 60 cells
*
* As written, this program program runs for 8192 time step and writes
* a single file called "obs-point".  If you want the absolute value
* of the S11 parameter for the patch, you will have to run this
* program twice: once with the microstrip and once with patch.  
* So, on my Linux system I would issue the following commands:
*
*  % cc -DABC1 -DABC2 -DABC3 -DABC4 -DABC5 -O sheen_patch.c -o sheen_patch -lm
*  % sheen_patch
*  % mv obs-point obs-point-strip
*  % cc -DPATCH -DABC1 -DABC2 -DABC3 -DABC4 -DABC5 -O sheen_patch.c -o sheen_patch -lm
*  % sheen_patch
*  % mv obs-point obs-point-patch
*

* Now you can obtain |S11| by subtracting obs-point-strip from
* obs-point-patch to obtain the reflected field.  Take the Fourier
* transform of the incident field (i.e., obs-point-strip) and
* incident field, divide on a term by term basis, and plot the
* magnitude of the result.  In matlab you would do this with commands
* such as (this can be snipped out and saved to a separate file and
* fed to matlab):

%--------------------------------------------------------------------------
% These matlab commands can be used to plot the results for the
% Sheen et. simulation.
%
% Note that this assumes a uniform spatial step size of del_s of
% 0.265 mm, a Courant number of 1/sqrt(3.1), and 8192.
% Using those numbers we can find the temporal step size.
% 
%   c del_t/del_s = 1/sqrt(3.1)   => del_t = del_s/(c sqrt(3.1))
%                                          = 5.016996 10^-13
%
% Thus the highest frequency in the simulation is 
%
%   f_max = 1/(2 del_t) = 9.9661227 10^11
%
% Given that there are 8192 time-steps in the simulation, the spectral
% resolution is
%
%   del_f = f_max/(number_of_time_steps / 2) = 243 313 543.8 Hz
%
% Thus del_f is 0.2433 GHz and since we'll plot the spectrum in GHz, 
% this is the how we'll scale the horizontal axis for the results from 
% the FFT.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
inc=dlmread('obs-point-strip','\n');  % "incident" field
tot=dlmread('obs-point-patch','\n');  % total field
% reflected field is the difference of total and incident field
ref=tot-inc;                          
% Take FFT of incident and relected fields.  S11 transfer function
% is just reflected divided by incident field.  Note that this is
% only meaningful at frequencies where we have sufficient incident 
% energy to excite the system, but this should be fine over the range
% of frequencies considered here.
incF=fft(inc);                        
refF=fft(ref);
freq=0.24331*(0:100);
semilogy(freq(5:80),abs(refF(5:80) ./ incF(5:80)))
title('S11 for patch antenna')
xlabel('Frequency, GHz')
ylabel('|S11|')
%--------------------------------------------------------------------------

* Note that the size of the computational domain is larger than that
* used by Sheen et al., so you may have to be patient while this runs
* (or resize things to suit your needs).
*
* Enjoy!
*/

#include <math.h>
#include <stdlib.h>
#include <stdio.h>

/* Size of the computational domain. */
#define LIMX  90
#define LIMY  130
#define LIMZ  20
#define HEIGHT 3  /* Number of cells for the dielectric substrate. */

/* Various start and stop point for the feed line and patch. */
#define LINE_X_START 28
#define LINE_X_END 37

#define PATCH_X_START 20
#define PATCH_X_END 67
#define PATCH_Y_START 50
#define PATCH_Y_END 110

/* Parameter to control the width of the Gaussian pulse, which is
   rather arbitrary.  We just need to ensure we have sufficient
   spectral energy at the frequencies of interest.  */
#define PPW 30

/* The time at which the source plane switches to being a regular
   (absorbing) wall. */
#define SWITCH_SRC 225

#ifndef M_PI
#define M_PI  3.14159265358979323846
#endif

double gaussian(int, double, int);
void init(void);

/* Field arrays */
double ex[LIMX][LIMY][LIMZ], ey[LIMX][LIMY][LIMZ],
       ez[LIMX][LIMY][LIMZ], hx[LIMX][LIMY][LIMZ],
       hy[LIMX][LIMY][LIMZ], hz[LIMX][LIMY][LIMZ];

/* I don't like repeated brackets so define some macros to make things
   neater. */
#define Ex(I,J,K) ex[I][J][K]
#define Ey(I,J,K) ey[I][J][K]
#define Ez(I,J,K) ez[I][J][K]
#define Hx(I,J,K) hx[I][J][K]
#define Hy(I,J,K) hy[I][J][K]
#define Hz(I,J,K) hz[I][J][K]

/* ABC arrays -- a second-oder Higdon ABC is used at most faces. */
double exfar[LIMX][3][LIMZ][2], ezfar[LIMX][3][LIMZ][2];

#define Exfar(I,J,K,N) exfar[I][J-(LIMY-3)][K][N]
#define Ezfar(I,J,K,N) ezfar[I][J-(LIMY-3)][K][N]

double extop[LIMX][LIMY][3][2], eytop[LIMX][LIMY][3][2];

#define Extop(I,J,K,N) extop[I][J][K-(LIMZ-3)][N]
#define Eytop(I,J,K,N) eytop[I][J][K-(LIMZ-3)][N]

double exnear[LIMX][3][LIMZ][2], eznear[LIMX][3][LIMZ][2];

#define Exnear(I,J,K,N) exnear[I][J][K][N]
#define Eznear(I,J,K,N) eznear[I][J][K][N]

double eyleft[3][LIMY][LIMZ][2], ezleft[3][LIMY][LIMZ][2];

#define Eyleft(I,J,K,N) eyleft[I][J][K][N]
#define Ezleft(I,J,K,N) ezleft[I][J][K][N]

double eyright[3][LIMY][LIMZ][2], ezright[3][LIMY][LIMZ][2];

#define Eyright(I,J,K,N) eyright[I-(LIMX-3)][J][K][N]
#define Ezright(I,J,K,N) ezright[I-(LIMX-3)][J][K][N]

int main() {
  double mu0, clight, eps0;
  double coefH, coefE0, coefE1, coefE01, holdez, cdtds;
  int i, j, k, ii, jj, kk, ntime, ntmax=8192;

  /* ABC parameters. */
  double c10, c20, c30, c40, c101, c201, c301, c401, c11, c21, c31, c41, temp;

  /* output file */
  FILE *obs;
  obs=fopen("obs-point","w");

  mu0=M_PI*4.e-7;
  clight=2.99792458e8;
  eps0=1.0/(mu0*clight*clight);

  /* Run simulation close to Courant limit. */
  cdtds   = 1./sqrt(3.1);
  coefH   = cdtds/clight/mu0;
  coefE0  = cdtds/clight/eps0;
  coefE1  = coefE0/2.2;
  coefE01 = coefE0/((1.0+2.2)/2.0);

  /* ABC coefficients. */
  temp = cdtds;
  c10 =    -(1.0/temp - 2.0 + temp)/(1.0/temp + 2.0 + temp);
  c20 = 2.0*(1.0/temp - temp)/(1.0/temp + 2.0 + temp);
  c30 = 4.0*(1.0/temp + temp)/(1.0/temp + 2.0 + temp);
  c40 = (temp-1.0)/(temp+1.0);
  temp = cdtds/sqrt((1.0 + 2.2)/2.0);
  c101 =    -(1.0/temp - 2.0 + temp)/(1.0/temp + 2.0 + temp);
  c201 = 2.0*(1.0/temp - temp)/(1.0/temp + 2.0 + temp);
  c301 = 4.0*(1.0/temp + temp)/(1.0/temp + 2.0 + temp);
  c401 = (temp-1.0)/(temp+1.0);
  temp = cdtds/sqrt(2.2);
  c11 =    -(1.0/temp - 2.0 + temp)/(1.0/temp + 2.0 + temp);
  c21 = 2.0*(1.0/temp - temp)/(1.0/temp + 2.0 + temp);
  c31 = 4.0*(1.0/temp + temp)/(1.0/temp + 2.0 + temp);
  c41 = (temp-1.0)/(temp+1.0);

  /* Initialize all the fields to zero. */
  init();

  for (ntime=0; ntime<ntmax; ntime++) {
    printf("Working on time step %d ...\n",ntime);

    /********** Ex update. *************/
    for (i=0; i<LIMX-1; i++)
      for (j=1; j<LIMY-1; j++)
        for (k=1; k<LIMZ-1; k++)
          if (k > HEIGHT)
            Ex(i,j,k) = Ex(i,j,k) + 
              coefE0*((Hz(i,j,k)-Hz(i,j-1,k)) - (Hy(i,j,k)-Hy(i,j,k-1)));
          else if (k == HEIGHT)
            Ex(i,j,k) = Ex(i,j,k) + 
              coefE01*((Hz(i,j,k)-Hz(i,j-1,k)) - (Hy(i,j,k)-Hy(i,j,k-1)));
          else
            Ex(i,j,k) = Ex(i,j,k) + 
              coefE1*((Hz(i,j,k)-Hz(i,j-1,k)) - (Hy(i,j,k)-Hy(i,j,k-1)));

#ifdef ABC4
    if (ntime<SWITCH_SRC) {
#endif
    /* Ex nodes on y=0 PMC wall have special updates. */
    j=0;
    for (i=0;i<LIMX-1;i++)
      for (k=1;k<LIMZ-1;k++)
        if (k > HEIGHT)
          Ex(i,j,k) = Ex(i,j,k) + 
            coefE0*(2.0*Hz(i,j,k) - (Hy(i,j,k)-Hy(i,j,k-1)));
        else if (k == HEIGHT)
          Ex(i,j,k) = Ex(i,j,k) + 
            coefE01*(2.0*Hz(i,j,k) - (Hy(i,j,k)-Hy(i,j,k-1)));
        else
          Ex(i,j,k) = Ex(i,j,k) + 
            coefE1*(2.0*Hz(i,j,k) - (Hy(i,j,k)-Hy(i,j,k-1)));
#ifdef ABC4
    }
#endif

#ifdef PATCH
    /* Zero Ex on metal. */
    k=HEIGHT;
    /* First the feed strip. */
    for (i=LINE_X_START; i<LINE_X_END; i++)
      for (j=0; j<LIMY/2-10; j++)
        Ex(i,j,k) = 0.0;
    /* Next the patch. */
    for (i=PATCH_X_START; i<PATCH_X_END; i++)
      for (j=PATCH_Y_START; j<=PATCH_Y_END; j++)
        Ex(i,j,k) = 0.0;
#else
    /* Zero Ex on metal microstrip. */
    k=HEIGHT;
    for (i=LINE_X_START; i<LINE_X_END; i++)
      for (j=0; j<LIMY; j++)
        Ex(i,j,k) = 0.0;
#endif

#ifdef ABC2
    /* ABC at far wall */
    /* Ex above substrate */
    for (i=0; i<LIMX-1; i++) {
      j=LIMY-1;
      for (k=HEIGHT+1; k<LIMZ-1; k++) {
        Ex(i,j,k) = c10*(Ex(i,j-2,k)+Exfar(i,j,k,1))
          + c20*(Exfar(i,j,k,0) + Exfar(i,j-2,k,0) - Ex(i,j-1,k)
          - Exfar(i,j-1,k,1))  + c30*Exfar(i,j-1,k,0) - Exfar(i,j-2,k,1);
      
        for(jj=LIMY-3; jj<LIMY; jj++) {
          Exfar(i,jj,k,1) = Exfar(i,jj,k,0);
          Exfar(i,jj,k,0) = Ex(i,jj,k);
        }
      }
    }
    /* Ex at substrate */
    for (i=0; i<LIMX-1; i++) {
      j=LIMY-1;
      k=HEIGHT;
      Ex(i,j,k) = c101*(Ex(i,j-2,k)+Exfar(i,j,k,1))
        + c201*(Exfar(i,j,k,0) + Exfar(i,j-2,k,0) - Ex(i,j-1,k)
        - Exfar(i,j-1,k,1))  + c301*Exfar(i,j-1,k,0) - Exfar(i,j-2,k,1);
      
      for(jj=LIMY-3; jj<LIMY; jj++) {
        Exfar(i,jj,k,1) = Exfar(i,jj,k,0);
        Exfar(i,jj,k,0) = Ex(i,jj,k);
      }
    }
    /* Ex below substrate */
    for (i=0; i<LIMX-1; i++) {
      j=LIMY-1;
      for (k=1; k<HEIGHT; k++) {
        Ex(i,j,k) = c11*(Ex(i,j-2,k)+Exfar(i,j,k,1))
          + c21*(Exfar(i,j,k,0) + Exfar(i,j-2,k,0) - Ex(i,j-1,k)
          - Exfar(i,j-1,k,1))  + c31*Exfar(i,j-1,k,0) - Exfar(i,j-2,k,1);
      
        for(jj=LIMY-3; jj<LIMY; jj++) {
          Exfar(i,jj,k,1) = Exfar(i,jj,k,0);
          Exfar(i,jj,k,0) = Ex(i,jj,k);
        }
      }
    }
#endif

#ifdef ABC4
    /* ABC at near wall -- only apply ABC after source introduced */
    if (ntime >= SWITCH_SRC) {
      /* Ex above substrate */
      for (i=0; i<LIMX-1; i++) {
        j=0;
        for (k=HEIGHT+1; k<LIMZ-1; k++) {
          Ex(i,j,k) = Exnear(i,j+1,k,0) + c40*(Ex(i,j+1,k)-Ex(i,j,k));
          /* Ex(i,j,k) = c10*(Ex(i,j+2,k)+Exnear(i,j,k,1))
            + c20*(Exnear(i,j,k,0) + Exnear(i,j+2,k,0) - Ex(i,j+1,k)
              - Exnear(i,j+1,k,1)) + c30*Exnear(i,j+1,k,0) - Exnear(i,j+2,k,1);
          */

          for(jj=2; jj>=0; jj--) {
            Exnear(i,jj,k,1) = Exnear(i,jj,k,0);
            Exnear(i,jj,k,0) = Ex(i,jj,k);
          }
        }
      }
      /* Ex at substrate */
      for (i=0; i<LIMX-1; i++) {
        j=0;
        k=HEIGHT;
        Ex(i,j,k) = Exnear(i,j+1,k,0) + c401*(Ex(i,j+1,k)-Ex(i,j,k));
        /* Ex(i,j,k) = c101*(Ex(i,j+2,k)+Exnear(i,j,k,1))
          + c201*(Exnear(i,j,k,0) + Exnear(i,j+2,k,0) - Ex(i,j+1,k)
          - Exnear(i,j+1,k,1)) + c301*Exnear(i,j+1,k,0) - Exnear(i,j+2,k,1);
        */

        for(jj=2; jj>=0; jj--) {
          Exnear(i,jj,k,1) = Exnear(i,jj,k,0);
          Exnear(i,jj,k,0) = Ex(i,jj,k);
        }
      }
      /* Ex below substrate */
      for (i=0; i<LIMX-1; i++) {
        j=0;
        for (k=1; k<HEIGHT; k++) {
          Ex(i,j,k) = Exnear(i,j+1,k,0) + c41*(Ex(i,j+1,k)-Ex(i,j,k));
          /* Ex(i,j,k) = c11*(Ex(i,j+2,k)+Exnear(i,j,k,1))
            + c21*(Exnear(i,j,k,0) + Exnear(i,j+2,k,0) - Ex(i,j+1,k)
            - Exnear(i,j+1,k,1)) + c31*Exnear(i,j+1,k,0) - Exnear(i,j+2,k,1);
          */

          for(jj=2; jj>=0; jj--) {
            Exnear(i,jj,k,1) = Exnear(i,jj,k,0);
            Exnear(i,jj,k,0) = Ex(i,jj,k);
          }
        }
      }
    }
#endif

#ifdef ABC5
    /* ABC at top */
    for (i=0; i<LIMX-1; i++)
      for (j=0; j<LIMY; j++) {
        k=LIMZ-1;
        Ex(i,j,k) = c10*(Ex(i,j,k-2)+Extop(i,j,k,1))
          + c20*(Extop(i,j,k,0) + Extop(i,j,k-2,0) - Ex(i,j,k-1)
          - Extop(i,j,k-1,1))  + c30*Extop(i,j,k-1,0) - Extop(i,j,k-2,1);
      
        for(kk=LIMZ-3; kk<LIMZ; kk++) {
          Extop(i,j,kk,1) = Extop(i,j,kk,0);
          Extop(i,j,kk,0) = Ex(i,j,kk);
        }
    }
#endif

    /*********** Ey update. ***********/
    for (i=1; i<LIMX-1; i++)
      for (j=0; j<LIMY-1; j++)
        for (k=1; k<LIMZ-1; k++)