sheen_patch.c

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

	    Ey(i,j,k) = Ey(i,j,k) + 
	      coefE0*((Hx(i,j,k)-Hx(i,j,k-1)) - (Hz(i,j,k)-Hz(i-1,j,k)));
	  else if (k == HEIGHT)
	    Ey(i,j,k) = Ey(i,j,k) + 
	      coefE01*((Hx(i,j,k)-Hx(i,j,k-1)) - (Hz(i,j,k)-Hz(i-1,j,k)));
	  else
	    Ey(i,j,k) = Ey(i,j,k) + 
	      coefE1*((Hx(i,j,k)-Hx(i,j,k-1)) - (Hz(i,j,k)-Hz(i-1,j,k)));
    
#ifdef PATCH
    /* Zero Ey 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++)
	Ey(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++)
	Ey(i,j,k) = 0.0;
#else
    /* Zero Ey on metal microstrip. */
    k=HEIGHT;
    for (i=LINE_X_START; i<=LINE_X_END; i++)
      for (j=0; j<LIMY; j++)
	Ey(i,j,k) = 0.0;
#endif

#ifdef ABC1
    /* ABC on left wall */
    /* Ey above substrate */
    i = 0;
    for (j=0; j<LIMY-1; j++) {
      for (k=HEIGHT+1; k<LIMZ; k++) {
	Ey(i,j,k) = c10*(Ey(i+2,j,k)+Eyleft(i,j,k,1))
	  + c20*(Eyleft(i,j,k,0) + Eyleft(i+2,j,k,0) - Ey(i+1,j,k)
          - Eyleft(i+1,j,k,1))  + c30*Eyleft(i+1,j,k,0) - Eyleft(i+2,j,k,1);
      
	for(ii=0; ii<3; ii++) {
	  Eyleft(ii,j,k,1) = Eyleft(ii,j,k,0);
	  Eyleft(ii,j,k,0) = Ey(ii,j,k);
	}
      }
    }
    /* Ey at substrate */
    i=0;
    for (j=0; j<LIMY-1; j++) {
      k=HEIGHT;
      Ey(i,j,k) = c101*(Ey(i+2,j,k)+Eyleft(i,j,k,1))
	+ c201*(Eyleft(i,j,k,0) + Eyleft(i+2,j,k,0) - Ey(i+1,j,k)
        - Eyleft(i+1,j,k,1))  + c301*Eyleft(i+1,j,k,0) - Eyleft(i+2,j,k,1);
      
      for(ii=0; ii<3; ii++) {
	Eyleft(ii,j,k,1) = Eyleft(ii,j,k,0);
	Eyleft(ii,j,k,0) = Ey(ii,j,k);
      }
    }
    /* Ey below substrate */
    i=0;
    for (j=0; j<LIMY-1; j++) {
      for (k=1; k<HEIGHT; k++) {
	Ey(i,j,k) = c11*(Ey(i+2,j,k)+Eyleft(i,j,k,1))
	  + c21*(Eyleft(i,j,k,0) + Eyleft(i+2,j,k,0) - Ey(i+1,j,k)
          - Eyleft(i+1,j,k,1))  + c31*Eyleft(i+1,j,k,0) - Eyleft(i+2,j,k,1);
      
	for(ii=0; ii<3; ii++) {
	  Eyleft(ii,j,k,1) = Eyleft(ii,j,k,0);
	  Eyleft(ii,j,k,0) = Ey(ii,j,k);
	}
      }
    }
#endif

#ifdef ABC3
    /* ABC on right wall */
    /* Ey above substrate */
    i = LIMX-1;
    for (j=0; j<LIMY-1; j++) {
      for (k=HEIGHT+1; k<LIMZ; k++) {
	Ey(i,j,k) = c10*(Ey(i-2,j,k)+Eyright(i,j,k,1))
	  + c20*(Eyright(i,j,k,0) + Eyright(i-2,j,k,0) - Ey(i-1,j,k)
          - Eyright(i-1,j,k,1))  + c30*Eyright(i-1,j,k,0) - Eyright(i-2,j,k,1);
      
	for(ii=LIMX-3; ii<LIMX; ii++) {
	  Eyright(ii,j,k,1) = Eyright(ii,j,k,0);
	  Eyright(ii,j,k,0) = Ey(ii,j,k);
	}
      }
    }
    /* Ey at substrate */
    i=LIMX-1;
    for (j=0; j<LIMY-1; j++) {
      k=HEIGHT;
      Ey(i,j,k) = c101*(Ey(i-2,j,k)+Eyright(i,j,k,1))
	+ c201*(Eyright(i,j,k,0) + Eyright(i-2,j,k,0) - Ey(i-1,j,k)
        - Eyright(i-1,j,k,1))  + c301*Eyright(i-1,j,k,0) - Eyright(i-2,j,k,1);
      
      for(ii=LIMX-3; ii<LIMX; ii++) {
	Eyright(ii,j,k,1) = Eyright(ii,j,k,0);
	Eyright(ii,j,k,0) = Ey(ii,j,k);
      }
    }
    /* Ey below substrate */
    i=LIMX-1;
    for (j=0; j<LIMY-1; j++) {
      for (k=1; k<HEIGHT; k++) {
	Ey(i,j,k) = c11*(Ey(i-2,j,k)+Eyright(i,j,k,1))
	  + c21*(Eyright(i,j,k,0) + Eyright(i-2,j,k,0) - Ey(i-1,j,k)
          - Eyright(i-1,j,k,1))  + c31*Eyright(i-1,j,k,0) - Eyright(i-2,j,k,1);
      
	for(ii=LIMX-3; ii<LIMX; ii++) {
	  Eyright(ii,j,k,1) = Eyright(ii,j,k,0);
	  Eyright(ii,j,k,0) = Ey(ii,j,k);
	}
      }
    }
#endif

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

    /************ Ez update. *************/
    for (i=1; i<LIMX-1; i++)
      for (j=1; j<LIMY-1; j++)
	for (k=0; k<LIMZ-1; k++)
	  if (k > HEIGHT)
	    Ez(i,j,k) = Ez(i,j,k)
	      + coefE0*((Hy(i,j,k)-Hy(i-1,j,k)) - (Hx(i,j,k)-Hx(i,j-1,k)));
	  else
	    Ez(i,j,k) = Ez(i,j,k)
	      + coefE1*((Hy(i,j,k)-Hy(i-1,j,k))	- (Hx(i,j,k)-Hx(i,j-1,k)));
    
#ifdef ABC4
    if (ntime<SWITCH_SRC) {
#endif
    /* Ez update on "near" wall. */
    j = 0;
    for (i=1; i<LIMX; i++)
      for (k=0; k<LIMZ; k++)
	if (k > HEIGHT)
	  Ez(i,j,k) = Ez(i,j,k)
	    + coefE0*((Hy(i,j,k)-Hy(i-1,j,k)) - 2.0*Hx(i,j,k));
	else
	  Ez(i,j,k) = Ez(i,j,k)
	    + coefE1*((Hy(i,j,k)-Hy(i-1,j,k)) - 2.0*Hx(i,j,k));
    
    /* source at y=0 wall of computational domain */
    holdez = gaussian(ntime, cdtds,  PPW)/3.0;
    j=0;
    for (i=LINE_X_START; i<=LINE_X_END; i++)
      for (k=0; k<HEIGHT; k++)
	Ez(i,j,k) = holdez;
    /* Ez(i,j,k) = Ez(i,j,k) + holdez; */
#ifdef ABC4
    }
#endif

#ifdef ABC1
    /* ABC on left wall */
    /* Ez above substrate */
    i = 0;
    for (j=0; j<LIMY-1; j++) {
      for (k=HEIGHT; k<LIMZ; k++) {
	Ez(i,j,k) = c10*(Ez(i+2,j,k)+Ezleft(i,j,k,1))
	  + c20*(Ezleft(i,j,k,0) + Ezleft(i+2,j,k,0) - Ez(i+1,j,k)
          - Ezleft(i+1,j,k,1))  + c30*Ezleft(i+1,j,k,0) - Ezleft(i+2,j,k,1);
      
      for(ii=0; ii<3; ii++) {
	  Ezleft(ii,j,k,1) = Ezleft(ii,j,k,0);
	  Ezleft(ii,j,k,0) = Ez(ii,j,k);
	}
      }
    }
    /* Ez below substrate */
    i=0;
    for (j=0; j<LIMY-1; j++) {
      for (k=0; k<HEIGHT; k++) {
	Ez(i,j,k) = c11*(Ez(i+2,j,k)+Ezleft(i,j,k,1))
	  + c21*(Ezleft(i,j,k,0) + Ezleft(i+2,j,k,0) - Ez(i+1,j,k)
          - Ezleft(i+1,j,k,1))  + c31*Ezleft(i+1,j,k,0) - Ezleft(i+2,j,k,1);
      
	for(ii=0; ii<3; ii++) {
	  Ezleft(ii,j,k,1) = Ezleft(ii,j,k,0);
	  Ezleft(ii,j,k,0) = Ez(ii,j,k);
	}
      }
    }
#endif

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

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

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

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

#ifdef ABC3
    /* ABC on right wall */
    /* Ez above substrate */
    i = LIMX-1;
    for (j=0; j<LIMY-1; j++) {
      for (k=HEIGHT; k<LIMZ; k++) {
	Ez(i,j,k) = c10*(Ez(i-2,j,k)+Ezright(i,j,k,1))
	  + c20*(Ezright(i,j,k,0) + Ezright(i-2,j,k,0) - Ez(i-1,j,k)
          - Ezright(i-1,j,k,1))  + c30*Ezright(i-1,j,k,0) - Ezright(i-2,j,k,1);
      
	for(ii=LIMX-3; ii<LIMX; ii++) {
	  Ezright(ii,j,k,1) = Ezright(ii,j,k,0);
	  Ezright(ii,j,k,0) = Ez(ii,j,k);
	}
      }
    }
    /* Ez below substrate */
    i=LIMX-1;
    for (j=0; j<LIMY-1; j++) {
      for (k=0; k<HEIGHT; k++) {
	Ez(i,j,k) = c11*(Ez(i-2,j,k)+Ezright(i,j,k,1))
	  + c21*(Ezright(i,j,k,0) + Ezright(i-2,j,k,0) - Ez(i-1,j,k)
          - Ezright(i-1,j,k,1))  + c31*Ezright(i-1,j,k,0) - Ezright(i-2,j,k,1);
      
	for(ii=LIMX-3; ii<LIMX; ii++) {
	  Ezright(ii,j,k,1) = Ezright(ii,j,k,0);
	  Ezright(ii,j,k,0) = Ez(ii,j,k);
	}
      }
    }
#endif


    /************ Hx update. ************/
    for (i=0; i<LIMX; i++)
      for (j=0; j<LIMY-1; j++)
	for (k=0; k<LIMZ-1; k++)
	  Hx(i,j,k) = Hx(i,j,k) 
	    + coefH*((Ey(i,j,k+1)-Ey(i,j,k)) - (Ez(i,j+1,k)-Ez(i,j,k)));
    
    /************ Hy update. ************/
    for (i=0; i<LIMX-1; i++)
      for (j=0; j<LIMY; j++)
	for (k=0; k<LIMZ-1; k++)
	  Hy(i,j,k) = Hy(i,j,k) 
	    + coefH*((Ez(i+1,j,k)-Ez(i,j,k)) - (Ex(i,j,k+1)-Ex(i,j,k)));
    
    /************ Hz update. ************/
    for (i=0; i<LIMX-1; i++)
      for (j=0; j<LIMY-1; j++)
	for (k=0; k<LIMZ; k++)
	  Hz(i,j,k) = Hz(i,j,k) 
	    + coefH*((Ex(i,j+1,k)-Ex(i,j,k)) - (Ey(i+1,j,k)-Ey(i,j,k)));
    
    fprintf(obs,"%f\n",Ez((LINE_X_START+LINE_X_END)/2,PATCH_Y_START-10,2));

  }

  fclose(obs);
  return 0;

}
/* ------------------------- end of main ------------------------- */

/* ########################## gaussian ########################### */
/* Gaussian pulse.
 */
double gaussian(int ntime, double cdtds, int ippw) {
  double arg, time;
  time = (double)(ntime);
  arg = pow(M_PI*(cdtds*time/(double)(ippw)-1.),2);
  /* If the argument is greater than 70.0, the result is so small we
     may as well skip calculating the exponential and return zero. */
  if (arg > 70.0) {
    return 0.0;
  } else {
    return exp(-arg);
  }
}
/* ------------------------ end of gaussian ------------------------ */

/* ########################### init ############################# */
void init() {

  int i, j, k;

  /* Initialize all fields. */
  /* Ex. */
  for (k=0;k<LIMZ;k++)
    for (j=0;j<LIMY;j++)
      for (i=0;i<LIMX;i++)
	Ex(i,j,k) = 0.;

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

  /* Ez. */
  for (k=0;k<LIMZ;k++)
    for (j=0;j<LIMY;j++)
      for (i=0;i<LIMX;i++)
	Ez(i,j,k) = 0.;

  /* Hx. */
  for (k=0;k<LIMZ;k++)
    for (j=0;j<LIMY;j++)
      for (i=0;i<LIMX;i++)
	Hx(i,j,k) = 0.;

  /* Hy. */
  for (i=0;i<LIMX;i++)
    for (j=0;j<LIMY; j++)
      for (k=0; k<LIMZ; k++)
	Hy(i,j,k) = 0.;

  /* Hz. */
  for (i=0; i<LIMX; i++)
    for (j=0; j<LIMY; j++)
      for (k=0; k<LIMZ; k++)
	Hz(i,j,k) = 0.;

  /* ABC arrays */
  for (i=0; i<LIMX; i++)
    for(j=LIMY-3; j<LIMY; j++)
      for(k=0; k<LIMZ; k++) {
	Exfar(i,j,k,0) = 0.0;
	Exfar(i,j,k,1) = 0.0;
	Ezfar(i,j,k,0) = 0.0;
	Ezfar(i,j,k,1) = 0.0;
	Exnear(i,j,k,0) = 0.0;
	Exnear(i,j,k,1) = 0.0;
	Eznear(i,j,k,0) = 0.0;
	Eznear(i,j,k,1) = 0.0;
      }

  for (i=0; i<3; i++)
    for(j=0; j<LIMY; j++)
      for(k=0; k<LIMZ; k++) {
	Eyleft(i,j,k,0) = 0.0;
	Eyleft(i,j,k,1) = 0.0;
	Ezleft(i,j,k,0) = 0.0;
	Ezleft(i,j,k,1) = 0.0;
      }

  for (i=0; i<LIMX; i++)
    for(j=0; j<LIMY; j++)
      for(k=0; k<3; k++) {
	Extop(i,j,k,0) = 0.0;
	Extop(i,j,k,1) = 0.0;
	Eytop(i,j,k,0) = 0.0;
	Eytop(i,j,k,1) = 0.0;
      }

  for (i=LIMX-3; i<LIMX; i++)
    for(j=0; j<LIMY; j++)
      for(k=0; k<LIMZ; k++) {
	Eyright(i,j,k,0) = 0.0;
	Eyright(i,j,k,1) = 0.0;
	Ezright(i,j,k,0) = 0.0;
	Ezright(i,j,k,1) = 0.0;
      }

}
/* ------------------------- end of init ------------------------ */