sheen_patch.c

This code can be used to model a microstrip line or a microstrip patch antenna (the particular prob

/*
 * 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++)
	  if (k > HEIGHT)