prepare.h

任何几何形状的导体或介质的麦克斯韦方程的时域有限差分代码

/* prepares the data for the main loop */
void prepare()
{
	long x,y,z;
	char i;
	long ind,indp1,indm1,Ra;
	long it;
	char d,c;
	double div,V;
	long Nxy=Nx*Ny;
	char str[1000];

/* start matlab */

	if (matlab==1)
	{
		if (!(ep = engOpen("\0")))
		{
			fprintf(stderr, "\nCan't start MATLAB engine\n");
			exit(0);
		} else
		{
			matlabON=1;
		}


/* initialize the field neccesary for matlab output */

		for (i=0;i<NumShow;i++)
		{
			if ((Show[i].type==0) || (Show[i].type==1) || (Show[i].type==2) || (Show[i].type==7) || (Show[i].type==8))
			{
				if (Show[i].plane=='x')
				{
					Show[i].F=(double *)malloc(sizeof(double)*Nz*Ny);
					Show[i].MF = mxCreateDoubleMatrix(Nz,Ny, mxREAL);
				}
				if (Show[i].plane=='y')
				{
					Show[i].F=(double *)malloc(sizeof(double)*Nx*Nz);
					Show[i].MF = mxCreateDoubleMatrix(Nz,Nx, mxREAL);
				}
				if (Show[i].plane=='z')
				{
					Show[i].F=(double *)malloc(sizeof(double)*Nx*Ny);
					Show[i].MF = mxCreateDoubleMatrix(Ny,Nx, mxREAL);
				}
				sprintf(str,"F%d",i);
				mxSetName(Show[i].MF,str);
			}
			if (Show[i].type==3)
			{
				Show[i].F=(double *)malloc(sizeof(double)*(Show[i].y2-Show[i].y1+1));
				Show[i].MF = mxCreateDoubleMatrix(1,Show[i].y2-Show[i].y1+1, mxREAL);
				
				sprintf(str,"F%d",i);
				mxSetName(Show[i].MF,str);
			}
			if (Show[i].type==4)
			{
				sprintf(str,"V%d=[0 0];t%d=[0 0];figure(%d);hold on;",i,i,i+1);
				engEvalString(ep, str);
			}

/* Initialize the movies */
			Show[i].ind=1;
		}
		
	}

/* initialize epsilon and sigma arrays */
	
	for (x=0;x<Nx;x++)
		for (y=0;y<Ny;y++)
			for (z=0;z<Nz;z++)
			{
				ind=(z*Ny+y)*Nx+x;
				for (i=0;i<6;i++)
					F [i][ind]=0.;
				e[ind]=e0;
				s[ind]=0.;
			}
			
/* include boxes */
			
	for (i=0;i<NumBox;i++)
		for (x=Box[i].x1;x<=Box[i].x2;x++)
			for (y=Box[i].y1;y<=Box[i].y2;y++)
				for (z=Box[i].z1;z<=Box[i].z2;z++)
				{
					s[(z*Ny+y)*Nx+x]=Box[i].sigma;
					e[(z*Ny+y)*Nx+x]*=Box[i].eps;
				}
	
/* include cones */

	for (i=0;i<NumCone;i++)
	{
		if (Cone[i].r1>Cone[i].r2) Ra=Cone[i].r1;else Ra=Cone[i].r2;
		
		for (x=Cone[i].x1;x<=Cone[i].x2;x++)
		for (y=Cone[i].y-Ra;y<=Cone[i].y+Ra;y++)
		for (z=Cone[i].z-Ra;z<=Cone[i].z+Ra;z++)
		
			if (sqrt((y-Cone[i].y)*(y-Cone[i].y)+(z-Cone[i].z)*(z-Cone[i].z))
				<=Cone[i].r1+(int)((double)(x)*(double)(Cone[i].r2-Cone[i].r1)/(double)(Cone[i].x2-Cone[i].x1)))
			
			{
				if (Cone[i].direction=='z')
				{
					s[(x*Ny+z)*Nx+y]=Cone[i].sigma;
					e[(x*Ny+z)*Nx+y]*=Cone[i].eps;
				}
				if (Cone[i].direction=='x')
				{
					s[(y*Ny+z)*Nx+x]=Cone[i].sigma;
					e[(y*Ny+z)*Nx+x]*=Cone[i].eps;
				}
				if (Cone[i].direction=='y')
				{
					s[(y*Ny+x)*Nx+z]=Cone[i].sigma;
					e[(y*Ny+x)*Nx+z]*=Cone[i].eps;
				}
			}
	}
	
/* choose timestep (50% of magical timestep)
	otherwice you get some instabilities maybe caused by the source or ABC's
*/

	if (dt==0.) dt=0.5*dx/SLight/sqrt(3.);

	printf("Timestep : %e\n",dt);
	
/* create coefficient fields, cE1 for E, cE2 for dH, cE3 for Voltage,cH for dE */
	
	for (x=0;x<Nx;x++)
		for (y=0;y<Ny;y++)
			for (z=0;z<Nz;z++)
			{
				ind=(z*Ny+y)*Nx+x;
				cE1[ind]=(1.-s[ind]*dt/2./e[ind])/(1.+s[ind]*dt/2./e[ind]);
				cE2[ind]=            dt/e[ind]/dx/(1.+s[ind]*dt/2./e[ind]);
				cE3[ind]=0.;
			}
			
	cH=dt/m0/dx;
	
/* include lumped resistors */

	for (i=0;i<NumR;i++)
	{
		div=R[i].z2-R[i].z1+1.;
		if (R[i].direction=='x') div=R[i].x2-R[i].x1+1;
		if (R[i].direction=='y') div=R[i].y2-R[i].y1+1;
		
		for (x=R[i].x1;x<=R[i].x2;x++)
		for (y=R[i].y1;y<=R[i].y2;y++)
		for (z=R[i].z1;z<=R[i].z2;z++)
		{
			ind=(z*Ny+y)*Nx+x;
			cE1[ind]=(1.-dt/(2.*R[i].R/div*e0*dx))	/(1.+dt/(2.*R[i].R/div*e0*dx));
			cE2[ind]=dt/e0/dx						/(1.+dt/(2.*R[i].R/div*e0*dx));
			cE3[ind]=0;
		}
	}

/* include lumped capacitors */

	for (i=0;i<NumC;i++)
	{
		div=C[i].z2-C[i].z1+1.;
		if (C[i].direction=='x') div=C[i].x2-C[i].x1+1;
		if (C[i].direction=='y') div=C[i].y2-C[i].y1+1;
		
		for (x=C[i].x1;x<=C[i].x2;x++)
		for (y=C[i].y1;y<=C[i].y2;y++)
		for (z=C[i].z1;z<=C[i].z2;z++)
		{
			ind=(z*Ny+y)*Nx+x;
			cE1[ind]=1.;
			cE2[ind]=dt/e0/dx						/(1.+C[i].C*div/(e0*dx));
			cE3[ind]=0;
		}
	}

/* apply source (means change coefficient fields cE for Voltage) */
		
	d=5; div=Volt.z2-Volt.z1+1.;
	if (Volt.direction=='x') {d=3;div=Volt.x2-Volt.x1+1;}
	if (Volt.direction=='y') {d=4;div=Volt.y2-Volt.y1+1;}

		for (x=Volt.x1;x<=Volt.x2;x++)
		for (y=Volt.y1;y<=Volt.y2;y++)
		for (z=Volt.z1;z<=Volt.z2;z++)
		{
			ind=(z*Ny+y)*Nx+x;
			cE1[ind]=(1.-dt/(2.*Volt.R/div*e0*dx))	/(1.+dt/(2.*Volt.R/div*e0*dx));
			cE2[ind]=dt/e0/dx						/(1.+dt/(2.*Volt.R/div*e0*dx));
			cE3[ind]=dt/(Volt.R*e0*dx*dx)			/(1.+dt/(2.*Volt.R/div*e0*dx));
		}

}