erweifd2d3[1].3.1.cpp

计算二维光子晶体带隙的源程序

#include <iostream>
#include "math.h"

using namespace std ;
#define PRECISION double  
//the precision of the data used in this program.It can be float or double

/************************************************************************
Fundamental constants
************************************************************************/
PRECISION pi=3.1415926535 ;
PRECISION cc=2.99792458e8 ;
//speed of light in free space
PRECISION muz=4.0*pi*1.0e-7 ;
//permeability of free space
PRECISION epsz=1.0/(cc*cc*muz);
//permittivity of free space

PRECISION freq=5.0e+14 ;
//center frequency of source excitation
PRECISION lambda=cc/freq ;
//center wavelength of source excitation
PRECISION omega=2.0*pi*freq ;

/*************************************************************************
*       Function declaration
**************************************************************************/

PRECISION **MyAllocPrecision2D(int dimensionX,int dimensionY,PRECISION value);
PRECISION *MyAllocPrecision1D(int N,PRECISION value);

int main()
{
    
/***********************************************************************
*     declaration parameters
    ************************************************************************/
    int i,j,n,a;
    
    //***********************************************************************
    //     Grid parameters
    //***********************************************************************/
    int ie=900 ;
    //number of grid cells in x-direction
    int je=900 ;
    //number of grid cells in y-direction
    
    int ib=ie+1 ;
    int jb=je+1 ;
    
    int is=10 ;
    //location of z-directed hard source
    int js=420 ;
    //location of z-directed hard source
    
	int ip=880; //location of probe
	int jp=420; //location of probe
	
	/************************************************************************/
	/*                                                                      */
	/************************************************************************/
	
    PRECISION dx=1.0e-8 ;
    //space increment of square lattice
    PRECISION dt=dx/(2.0*cc);
    //time step
    
    int nmax;
	cout<<"Input total Time steps:";
	cin>>nmax;
    //total number of time steps
    
    int iebc=8 ;
    //thickness of left and right PML region
    int jebc=8 ;
    //thickness of front and back PML region
    PRECISION rmax=0.00001 ;
    int orderbc=2 ;
    int ibbc=iebc+1 ;
    int jbbc=jebc+1 ;
    int iefbc=ie+2*iebc ;
    int jefbc=je+2*jebc ;
    int ibfbc=iefbc+1 ;
    int jbfbc=jefbc+1 ;
    
    /***********************************************************************
	Material parameters
    ***********************************************************************/
    
    int media=2 ;
    //第一列留空，第二列为背景介质参数，第三列为介质柱参数
    PRECISION eps[3]=
    {
        0.0,2.1,1.0 
    }
    ;
    PRECISION sig[3]=
    {
        0.0,0.0,0.0 
    }
    ;
    PRECISION mur[3]=
    {
        0.0,1.0,1.0 
    }
    ;
    PRECISION sim[3]=
    {
        0.0,0.0,0.0 
    }
    ;
    
    /***********************************************************************
	Wave excitation
    ***********************************************************************/
	PRECISION rtau=160.0e-17 ;
    PRECISION tau=rtau/dt ;
    PRECISION delay=3*tau ;

    PRECISION*source=MyAllocPrecision1D(nmax,0.0);
    int sm;
	sm=7*tau;
	if (nmax<sm) sm=nmax;
    
	for(n=1;n<=sm;n++)
    {
        source[n]=sin(omega*(n-delay)*dt)*exp(-(pow((n-delay),2)/pow(tau,2)));
		//source[n]=exp(-(4*pi*(n*dt-T0)*(n*dt-T0))/(TAO*TAO));
    }
    
	FILE *FpSource;
	FpSource=fopen("source.txt","w");
	for(i=1;i<=nmax;i++) fprintf(FpSource,"%d %e\n",i,source[i]);
	fclose(FpSource);
	
	

    /***********************************************************************
	Field arrays
    /***********************************************************************/
    
    PRECISION**ex=MyAllocPrecision2D(ie,jb,0.0);
    // fields in main grid 
    PRECISION**ey=MyAllocPrecision2D(ib,je,0.0);
    PRECISION**hz=MyAllocPrecision2D(ie,je,0.0);
    
    PRECISION**exbcf=MyAllocPrecision2D(iefbc,jebc,0.0);
    //fields in front PML region
    PRECISION**eybcf=MyAllocPrecision2D(ibfbc,jebc,0.0);
    PRECISION**hzxbcf=MyAllocPrecision2D(iefbc,jebc,0.0);
    PRECISION**hzybcf=MyAllocPrecision2D(iefbc,jebc,0.0);
    
    PRECISION**exbcb=MyAllocPrecision2D(iefbc,jbbc,0.0);
    //fields in back PML region
    PRECISION**eybcb=MyAllocPrecision2D(ibfbc,jebc,0.0);
    PRECISION**hzxbcb=MyAllocPrecision2D(iefbc,jebc,0.0);
    PRECISION**hzybcb=MyAllocPrecision2D(iefbc,jebc,0.0);
    
    PRECISION**exbcl=MyAllocPrecision2D(iebc,jb,0.0);
    //fields in left PML region
    PRECISION**eybcl=MyAllocPrecision2D(iebc,je,0.0);
    PRECISION**hzxbcl=MyAllocPrecision2D(iebc,je,0.0);
    PRECISION**hzybcl=MyAllocPrecision2D(iebc,je,0.0);
    
    PRECISION**exbcr=MyAllocPrecision2D(iebc,jb,0.0);
    //fields in right PML region
    PRECISION**eybcr=MyAllocPrecision2D(ibbc,je,0.0);
    PRECISION**hzxbcr=MyAllocPrecision2D(iebc,je,0.0);
    PRECISION**hzybcr=MyAllocPrecision2D(iebc,je,0.0);
    
    /***********************************************************************
	Updating coefficients
    ***********************************************************************/
    
    PRECISION eaf ;
    PRECISION*ca=MyAllocPrecision1D(media,0.0);
    PRECISION*cb=MyAllocPrecision1D(media,0.0);
    PRECISION haf ;
    PRECISION*da=MyAllocPrecision1D(media,0.0);
    PRECISION*db=MyAllocPrecision1D(media,0.0);
    
    for(i=1;i<=media;i++)
    {
        eaf=dt*sig[i]/(2.0*epsz*eps[i]);
        ca[i]=(1.0-eaf)/(1.0+eaf);
        cb[i]=dt/epsz/eps[i]/dx/(1.0+eaf);
        haf=dt*sim[i]/(2.0*muz*mur[i]);
        da[i]=(1.0-haf)/(1.0+haf);
        db[i]=dt/muz/mur[i]/dx/(1.0+haf);
    }
    
    /***********************************************************************
	Geometry specification (main grid)
    ***********************************************************************/
    
    //     Initialize entire main grid to free space
    
    PRECISION**caex=MyAllocPrecision2D(ie,jb,ca[1]);
    PRECISION**cbex=MyAllocPrecision2D(ie,jb,cb[1]);
    
    PRECISION**caey=MyAllocPrecision2D(ib,je,ca[1]);
    PRECISION**cbey=MyAllocPrecision2D(ib,je,cb[1]);
    
    PRECISION**dahz=MyAllocPrecision2D(ie,je,da[1]);
    PRECISION**dbhz=MyAllocPrecision2D(ie,je,db[1]);
    
    //     Add metal cylinder
/*
    int diam=40 ;
    // diameter of cylinder: 6 cm
    int rad=diam/2.0 ;
    // radius of cylinder: 3 cm
    int icenter ;
    // i-coordinate of cylinder's center
    int jcenter;
    // j-coordinate of cylinder's center
	
	int ldist=20;
	int lattice=50;
	int p,q;    
    PRECISION dist2 ;

	for(p=0;p<16;p++)
		for(q=0;q<16;q++){
			icenter=(q+0.5)*lattice+ldist;
			jcenter=(p+0.5)*lattice+ldist;
			
			for(i=1;i<=ie;i++)
			{
				for(j=1;j<=je;j++)
				{
					dist2=pow((i+0.5-icenter),2.0)+pow((j-jcenter),2.0);
					if(dist2<=rad*rad)
					{
						caex[i][j]=ca[2];
						cbex[i][j]=cb[2];
					}
					dist2=pow((i-icenter),2.0)+pow((j+0.5-jcenter),2.0);
					if(dist2<=rad*rad)
					{
						caey[i][j]=ca[2];
						cbey[i][j]=cb[2];
					}
				}
			}
		}
						
	/********************************************************************/
	/*					  Fill the PML regions							*/
    /********************************************************************/
    
    PRECISION delbc=iebc*dx ;
    PRECISION sigmam=-log(rmax/100.0)*epsz*cc*(orderbc+1)/(2*delbc);
    PRECISION bcfactor=eps[1]*sigmam/(dx*(pow(delbc,orderbc))*(orderbc+1));
    
    //    FRONT region 
    PRECISION**caexbcf=MyAllocPrecision2D(iefbc,jebc,1.0);
    PRECISION**cbexbcf=MyAllocPrecision2D(iefbc,jebc,0.0); //////////////modified
    
    PRECISION y1,y2,sigmay,sigmays,ca1,cb1 ;
    for(j=2;j<=jebc;j++)
    {
        y1=(jebc-j+1.5)*dx ;
        y2=(jebc-j+0.5)*dx ;
        sigmay=bcfactor*(pow(y1,(orderbc+1))-pow(y2,(orderbc+1)));
        ca1=exp(-sigmay*dt/(epsz*eps[1]));
        cb1=(1.0-ca1)/(sigmay*dx);
        for(i=1;i<=iefbc;i++)
        {
            caexbcf[i][j]=ca1 ;
            cbexbcf[i][j]=cb1 ;
        }
    }
    
    sigmay=bcfactor*pow((0.5*dx),(orderbc+1));
    ca1=exp(-sigmay*dt/(epsz*eps[1]));
    cb1=(1-ca1)/(sigmay*dx);
    
    for(i=1;i<=ie;i++)
    {
        caex[i][1]=ca1 ;
        cbex[i][1]=cb1 ;
    }
    
    PRECISION**caexbcl=MyAllocPrecision2D(iebc,jb,ca1);
    PRECISION**cbexbcl=MyAllocPrecision2D(iebc,jb,cb1);
    PRECISION**caexbcr=MyAllocPrecision2D(iebc,jb,ca1);
    PRECISION**cbexbcr=MyAllocPrecision2D(iebc,jb,cb1);
    
    PRECISION da1,db1 ;
    PRECISION**dahzybcf=MyAllocPrecision2D(iefbc,jebc,0.0);
    PRECISION**dbhzybcf=MyAllocPrecision2D(iefbc,jebc,0.0);
    PRECISION**caeybcf=MyAllocPrecision2D(ibfbc,jebc,0.0);
    PRECISION**cbeybcf=MyAllocPrecision2D(ibfbc,jebc,0.0);
    PRECISION**dahzxbcf=MyAllocPrecision2D(iefbc,jebc,0.0);
    PRECISION**dbhzxbcf=MyAllocPrecision2D(iefbc,jebc,0.0);
    
    for(j=1;j<=jebc;j++)
    {
        y1=(jebc-j+1)*dx ;
        y2=(jebc-j)*dx ;
        sigmay=bcfactor*(pow(y1,(orderbc+1))-pow(y2,(orderbc+1)));
        sigmays=sigmay*(muz/(epsz*eps[1]));
        da1=exp(-sigmays*dt/muz);
        db1=(1-da1)/(sigmays*dx);
        for(i=1;i<=iefbc;i++)
        {
            dahzybcf[i][j]=da1 ;
            dbhzybcf[i][j]=db1 ;
            
            dahzxbcf[i][j]=da[1];
            dbhzxbcf[i][j]=db[1];
        }
        
        for(i=1;i<=ibfbc;i++)
        {
            caeybcf[i][j]=ca[1];
            cbeybcf[i][j]=cb[1];
        }
        
        
    }
    
    
    
    //    BACK region 
    
    PRECISION**caexbcb=MyAllocPrecision2D(iefbc,jbbc,1.0);
    PRECISION**cbexbcb=MyAllocPrecision2D(iefbc,jbbc,0.0);
    for(j=2;j<=jebc;j++)
    {
        y1=(j-0.5)*dx ;
        y2=(j-1.5)*dx ;
        sigmay=bcfactor*(pow(y1,(orderbc+1))-pow(y2,(orderbc+1)));
        ca1=exp(-sigmay*dt/(epsz*eps[1]));
        cb1=(1-ca1)/(sigmay*dx);
        for(i=1;i<=iefbc;i++)
        {
            caexbcb[i][j]=ca1 ;
            cbexbcb[i][j]=cb1 ;
        }
        
    }
    
    sigmay=bcfactor*pow((0.5*dx),(orderbc+1));
    ca1=exp(-sigmay*dt/(epsz*eps[1]));
    cb1=(1-ca1)/(sigmay*dx);
    for(i=1;i<=ie;i++)
    {
        caex[i][jb]=ca1 ;
        cbex[i][jb]=cb1 ;
    }
    
    for(i=1;i<=iebc;i++)
    {
        caexbcl[i][jb]=ca1 ;
        cbexbcl[i][jb]=cb1 ;
        caexbcr[i][jb]=ca1 ;
        cbexbcr[i][jb]=cb1 ;
    }
    
    
    PRECISION**dahzybcb=MyAllocPrecision2D(iefbc,jebc,0.0);
    PRECISION**dbhzybcb=MyAllocPrecision2D(iefbc,jebc,0.0);
    PRECISION**caeybcb=MyAllocPrecision2D(ibfbc,jebc,0.0);
    PRECISION**cbeybcb=MyAllocPrecision2D(ibfbc,jebc,0.0);
    PRECISION**dahzxbcb=MyAllocPrecision2D(iefbc,jebc,0.0);
    PRECISION**dbhzxbcb=MyAllocPrecision2D(iefbc,jebc,0.0);
    
    for(j=1;j<=jebc;j++)
    {
        y1=j*dx ;
        y2=(j-1)*dx ;
        sigmay=bcfactor*(pow(y1,(orderbc+1))-pow(y2,(orderbc+1)));
        sigmays=sigmay*(muz/(epsz*eps[1]));
        da1=exp(-sigmays*dt/muz);
        db1=(1-da1)/(sigmays*dx);
        
        for(i=1;i<=iefbc;i++)
        {
            dahzybcb[i][j]=da1 ;
            dbhzybcb[i][j]=db1 ;
            dahzxbcb[i][j]=da[1];
            dbhzxbcb[i][j]=db[1];
        }
        
        for(i=1;i<=ibfbc;i++)
        {
            caeybcb[i][j]=ca[1];
            cbeybcb[i][j]=cb[1];
        }
        
    }
    
    //     LEFT region 
    
    PRECISION**caeybcl=MyAllocPrecision2D(iebc,je,1.0);
    PRECISION**cbeybcl=MyAllocPrecision2D(iebc,je,0.0);
    PRECISION x1,x2,sigmax,sigmaxs ;
    
    for(i=2;i<=iebc;i++)
    {
        x1=(iebc-i+1.5)*dx ;
        x2=(iebc-i+0.5)*dx ;
        sigmax=bcfactor*(pow(x1,(orderbc+1))-pow(x2,(orderbc+1)));
        ca1=exp(-sigmax*dt/(epsz*eps[1]));
        cb1=(1-ca1)/(sigmax*dx);
        for(j=1;j<=je;j++)
        {
            caeybcl[i][j]=ca1 ;
            cbeybcl[i][j]=cb1 ;
        }
        
        for(j=1;j<=jebc;j++)
        {
            caeybcf[i][j]=ca1 ;
            cbeybcf[i][j]=cb1 ;
            caeybcb[i][j]=ca1 ;
            cbeybcb[i][j]=cb1 ;
        }
    }
    sigmax=bcfactor*pow((0.5*dx),(orderbc+1));
    ca1=exp(-sigmax*dt/(epsz*eps[1]));
    cb1=(1-ca1)/(sigmax*dx);
    for(j=1;j<=je;j++)
    {
        caey[1][j]=ca1 ;
        cbey[1][j]=cb1 ;
    }
    
    for(j=1;j<=jebc;j++)
    {
        caeybcf[iebc+1][j]=ca1 ;
        cbeybcf[iebc+1][j]=cb1 ;
        caeybcb[iebc+1][j]=ca1 ;
        cbeybcb[iebc+1][j]=cb1 ;
    }
    
    PRECISION**dahzxbcl=MyAllocPrecision2D(iebc,je,0.0);
    PRECISION**dbhzxbcl=MyAllocPrecision2D(iebc,je,0.0);
    PRECISION**dahzybcl=MyAllocPrecision2D(iebc,je,0.0);
    PRECISION**dbhzybcl=MyAllocPrecision2D(iebc,je,0.0);
    
    for(i=1;i<=iebc;i++)
    {
        x1=(iebc-i+1)*dx ;
        x2=(iebc-i)*dx ;
        sigmax=bcfactor*(pow(x1,(orderbc+1))-pow(x2,(orderbc+1)));
        sigmaxs=sigmax*(muz/(epsz*eps[1]));
        da1=exp(-sigmaxs*dt/muz);
        db1=(1-da1)/(sigmaxs*dx);
        
        for(j=1;j<=je;j++)
        {
            dahzxbcl[i][j]=da1 ;
            dbhzxbcl[i][j]=db1 ;
            dahzybcl[i][j]=da[1];
            dbhzybcl[i][j]=db[1];
        }
        
        for(j=2;j<=je;j++)
        {
            caexbcl[i][j]=ca[1];
            cbexbcl[i][j]=cb[1];
        }
        
        for(j=1;j<=jebc;j++)
        {
            dahzxbcf[i][j]=da1 ;
            dbhzxbcf[i][j]=db1 ;
            dahzxbcb[i][j]=da1 ;
            dbhzxbcb[i][j]=db1 ;
        }
        
    }
    
    
    
    //     RIGHT region 
    
    PRECISION**caeybcr=MyAllocPrecision2D(ibbc,je,1.0);
    PRECISION**cbeybcr=MyAllocPrecision2D(ibbc,je,0.0);
    
    for(i=2;i<=iebc;i++)
    {
        x1=(i-0.5)*dx ;
        x2=(i-1.5)*dx ;
        sigmax=bcfactor*(pow(x1,(orderbc+1))-pow(x2,(orderbc+1)));
        ca1=exp(-sigmax*dt/(epsz*eps[1]));