one2.c

这几个程序是用于计算一维光子晶体反射系数和透射系数的

/* this is a one-dimensional FDTD program using mur's first-order ABC

   where kmax is the total number of cells used, a dilectric layer

   (thickness d=0.5e-3 and er = 4.0, nl=20 is the number of cells in the
 
 dilectric layer,so, cell size zz=0.5e-3/nl),itmax=800 is number of
 
 iterations. */

/* Author: Anping Zhao, Nokia Research Center, Helsinki, Finland, 
           Email: [e-mail]  */

#include <stdio.h>

#include <stdlib.h>

#include <math.h>


int main ()

{
#define kmax 50


double  eps0=8.854e-12, xmu0=1.2566e-6, er=1.0;

int itmax=800,nl=20;
int k,i,n1,n2,nstop,mm;

double c0,zz,ttt,tt,t,t0,pi,c1,c2,c3,fmax;

double e[kmax+1][2],h[kmax+1][2],einc[kmax+1][2],eps[kmax+1];

FILE *fp;

fp=fopen("cin","w");


/*  mm=5 means that excitation is placed five cells inside the

    computational region. The value of mm can be "arbitrary", but

    the source excitation must be palced between the open boundary

    and the discontinuity.*/

	mm=5;

	c0=1.0/(sqrt(eps0*xmu0));

	zz=0.5e-3/nl;


/*  ttt is the maximum value of time step */

	ttt=zz/c0;


/*  tt is the time step used in the simulation */

	tt=0.95*ttt;


/*  c1 is the coefficient for applying mur's abc */

	c1=(c0*tt-zz)/(c0*tt+zz);


/*  t and t0 are parameters used in Gaussian pulse.*/

	t=30.0*tt;

	t0=3.0*t;

        fmax=1.0/(2.0*t);

	pi=4.0*atan(1.0);


/*  nstop means that after nstop time steps the excitation should be
 
    stop (or the Gaussian pulse is totally lanched).*/

	nstop=2*t0/tt;


/*  er is the relative permittivity of the layer. if you chose er=4 then
 
    you get total (incident + reflected) fields, but if you chose er=1

    than you only get incident fields */


printf("kmax = %d, er = %e, tt = %e\n",kmax,er, tt);

printf("nstop = %d, zz = %e, fmax=%e\n",nstop, zz,fmax/1.0e9);


/*  c2 and c3 are the coefficients used in maxwell's equations.*/

	c2=tt/zz;

	c3=c2/xmu0;

	for (k=0;k<kmax+1;k++){

	  for (i=0;i<2;i++){

	e[k][i]=0.0;

	einc[k][i]=0.0;

	h[k][i]=0.0;

	  }

	}

	for (i=0;i<kmax+1;i++){

	  if(i>=(kmax-nl)/2 && i<=(kmax-nl)/2+nl){

	    eps[i]=er*eps0;}

	  else{

	eps[i]=eps0;

	  }

	}


/*  the following two line is to get average values for er at
 
    interfaces. */

	eps[(kmax-nl)/2]=((er+1)/2)*eps0;

	eps[(kmax-nl)/2+nl]=((er+1)/2)*eps0;


/*  updating process control */

	n1=1;

	n2=0;

	for (i=1;i<itmax+1;i++){

	  if(i%100==0){

	    printf("i = %d\n",i);}

	  if(n1==0){

	n1=1;

	n2=0;

	  }

	  else{

	n1=0;

	n2=1;

	  }

	  for (k=0;k<kmax+1;k++){
/*  this is for giving excitation field */

	    if (i <= nstop && k == mm) {

	      einc[k][n1]=exp(-((i-1)*tt-t0)*((i-1)*tt-t0)/(t*t));

	    }


/*  for maxwell's equation of H field */
 
	    if(k!=kmax){

	h[k][n2]=h[k][n1]-c3*(e[k+1][n1]-e[k][n1]);

	    }


/*  for maxwell's equation of E field */

	    if(k!=0 && k!=kmax){

	e[k][n2]=e[k][n1]-(c2/eps[k])*(h[k][n2]-h[k-1][n2])+einc[k][n1];

	    }


/*  mur's first-order abc for e field at k=0.
 
    but in 'if'  k must be 1 (not 0)! */

	    if(k == 1) {

	e[0][n2]=e[1][n1]+c1*(e[1][n2]-e[0][n1]);

	    }


/*  mur's abc for e field at k=kmax */

	    if(k == kmax) {

	e[kmax][n2]=e[kmax-1][n1]+c1*(e[kmax-1][n2]-e[kmax][n1]);

	    }

	  }

fprintf(fp,"%d%15.10f%15.10f\n",i,e[(kmax-nl)/2][n2],e[(kmax-nl)/2+nl][n2]);

	    }


    fclose (fp);

    return 0;

}