calc.h

FDTD 可以用来创建mpeg

/*	main calculation and matlab output */

void calc()
{
	long x,y,z;
	int i;
	long ind,indp1,indm1;
	long it,k,t;
	char d,c;
	double div,Vx,Vy,Vz;
	long Nxy=Nx*Ny;
	char str[1000];
	double Jc,Jd;
	double V,l;



/* start iteration - main loop*/

	it=0;
	
	while (it*dt<tend)
	{
		
/* save electrical field boundaries for ABC */
		
		ind=0;
		for (y=0;y<Ny;y++)
			for (x=0;x<Nx;x++)
			{
				oF[3][ind+Nxy]       =F[3][ind+Nxy];
				oF[3][ind+(Nz-2)*Nxy]=F[3][ind+(Nz-2)*Nxy];
				oF[4][ind+Nxy]       =F[4][ind+Nxy];
				oF[4][ind+(Nz-2)*Nxy]=F[4][ind+(Nz-2)*Nxy];
				oF[3][ind]           =F[3][ind];
				oF[3][ind+(Nz-1)*Nxy]=F[3][ind+(Nz-1)*Nxy];
				oF[4][ind]           =F[4][ind];
				oF[4][ind+(Nz-1)*Nxy]=F[4][ind+(Nz-1)*Nxy];
				ind++;
			}
		
		ind=0;
		for (z=0;z<Nz;z++)
			for (y=0;y<Ny;y++)
			{
				oF[4][ind+1]   =F[4][ind+1];
				oF[4][ind+Nx-2]=F[4][ind+Nx-2];
				oF[5][ind+1]   =F[5][ind+1];
				oF[5][ind+Nx-2]=F[5][ind+Nx-2];
				oF[4][ind]     =F[4][ind];
				oF[4][ind+Nx-1]=F[4][ind+Nx-1];
				oF[5][ind]     =F[5][ind];
				oF[5][ind+Nx-1]=F[5][ind+Nx-1];
				ind+=Nx;
			}
		
		ind=0;
		for (z=0;z<Nz;z++)
		{
			for (x=0;x<Nx;x++)
			{
				oF[3][ind+Nx]       =F[3][ind+Nx];
				oF[3][ind+(Ny-2)*Nx]=F[3][ind+(Ny-2)*Nx];
				oF[5][ind+Nx]       =F[5][ind+Nx];
				oF[5][ind+(Ny-2)*Nx]=F[5][ind+(Ny-2)*Nx];
				oF[3][ind]          =F[3][ind];
				oF[3][ind+(Ny-1)*Nx]=F[3][ind+(Ny-1)*Nx];
				oF[5][ind]          =F[5][ind];
				oF[5][ind+(Ny-1)*Nx]=F[5][ind+(Ny-1)*Nx];
				ind++;
			}
			ind+=Nxy-Nx;
		}
		
/* save whole E-fields for dD/dt (displacement current) */

		if (Save_old == 1)
			for (i=3;i<=5;i++)
				for (ind=0;ind<Nx*Ny*Nz;ind++)
					oF[i][ind]=F[i][ind];

/* set voltage value */

		Vx=Vy=Vz=0.;
		if (it*dt<Volt.timeon)
			if (Volt.direction=='x') Vx=it*dt/Volt.timeon*Volt.voltage;
			else
			if (Volt.direction=='y') Vy=it*dt/Volt.timeon*Volt.voltage;
			else
			if (Volt.direction=='z') Vz=it*dt/Volt.timeon*Volt.voltage;
			else;
		else
		if (it*dt<Volt.time+Volt.timeon)
			if (Volt.direction=='x') Vx=Volt.voltage;
			else
			if (Volt.direction=='y') Vy=Volt.voltage;
			else
			if (Volt.direction=='z') Vz=Volt.voltage;
			else;

/* fdtd formulation a la Lee */
		
		ind=0;
		indp1=(1*Ny+1)*Nx+1;
		
		for (z=0;z<Nz-1;z++)
		{
			for (y=0;y<Ny-1;y++)
			{
				for (x=0;x<Nx-1;x++)
				{
					F[0][ind]+=cH*(F[4][indp1]-F[4][indp1-Nxy]-F[5][indp1]+F[5][indp1-Nx]);
					F[1][ind]+=cH*(F[5][indp1]-F[5][indp1-1]  -F[3][indp1]+F[3][indp1-Nxy]);
					F[2][ind]+=cH*(F[3][indp1]-F[3][indp1-Nx] -F[4][indp1]+F[4][indp1-1]);

					ind++;
					indp1++;
				}
				ind++;
				indp1++;
			}
			ind+=Nx;
			indp1+=Nx;
		}
		
		indm1=0;
		ind=(1*Ny+1)*Nx+1;
		
		for (z=0;z<Nz-2;z++)
		{
			for (y=0;y<Ny-2;y++)
			{
				for (x=0;x<Nx-2;x++)
				{
					F[3][ind]=cE1[ind]*F[3][ind]+cE2[ind]*(F[2][indm1+Nx] -F[2][indm1]-F[1][indm1+Nxy]+F[1][indm1])+cE3[ind]*Vx;
					F[4][ind]=cE1[ind]*F[4][ind]+cE2[ind]*(F[0][indm1+Nxy]-F[0][indm1]-F[2][indm1+1]  +F[2][indm1])+cE3[ind]*Vy;
					F[5][ind]=cE1[ind]*F[5][ind]+cE2[ind]*(F[1][indm1+1]  -F[1][indm1]-F[0][indm1+Nx] +F[0][indm1])+cE3[ind]*Vz;

					ind++;
					indm1++;
				}
				ind+=2;
				indm1+=2;
			}
			ind+=Nx+Nx;
			indm1+=Nx+Nx;
		}
	
/* absorbing boundary conditions */

		ind=0;
		for (y=0;y<Ny;y++)
			for (x=0;x<Nx;x++)
			{
				F[3][ind]           =oF[3][ind+Nxy]       +(SLight*dt-dx)/(SLight*dt+dx)*(F[3][ind+Nxy]       -oF[3][ind]             );
				F[3][ind+(Nz-1)*Nxy]=oF[3][ind+(Nz-2)*Nxy]+(SLight*dt-dx)/(SLight*dt+dx)*(F[3][ind+(Nz-2)*Nxy]-oF[3][ind+(Nz-1)*Nxy]);
				F[4][ind]           =oF[4][ind+Nxy]       +(SLight*dt-dx)/(SLight*dt+dx)*(F[4][ind+Nxy]       -oF[4][ind]             );
				F[4][ind+(Nz-1)*Nxy]=oF[4][ind+(Nz-2)*Nxy]+(SLight*dt-dx)/(SLight*dt+dx)*(F[4][ind+(Nz-2)*Nxy]-oF[4][ind+(Nz-1)*Nxy]);
				ind++;
			}
		
		ind=0;
		for (z=0;z<Nz;z++)
			for (y=0;y<Ny;y++)
			{
				F[4][ind]     =oF[4][ind+1]   +(SLight*dt-dx)/(SLight*dt+dx)*(F[4][ind+1]   -oF[4][ind]     );
				F[4][ind+Nx-1]=oF[4][ind+Nx-2]+(SLight*dt-dx)/(SLight*dt+dx)*(F[4][ind+Nx-2]-oF[4][ind+Nx-1]);
				F[5][ind]     =oF[5][ind+1]   +(SLight*dt-dx)/(SLight*dt+dx)*(F[5][ind+1]   -oF[5][ind]     );
				F[5][ind+Nx-1]=oF[5][ind+Nx-2]+(SLight*dt-dx)/(SLight*dt+dx)*(F[5][ind+Nx-2]-oF[5][ind+Nx-1]);
				ind+=Nx;
			}
		
		ind=0;
		for (z=0;z<Nz;z++)
		{
			for (x=0;x<Nx;x++)
			{
				F[3][ind]          =oF[3][ind+Nx]       +(SLight*dt-dx)/(SLight*dt+dx)*(F[3][ind+Nx]       -oF[3][ind]          );
				F[3][ind+(Ny-1)*Nx]=oF[3][ind+(Ny-2)*Nx]+(SLight*dt-dx)/(SLight*dt+dx)*(F[3][ind+(Ny-2)*Nx]-oF[3][ind+(Ny-1)*Nx]);
				F[5][ind]          =oF[5][ind+Nx]       +(SLight*dt-dx)/(SLight*dt+dx)*(F[5][ind+Nx]       -oF[5][ind]          );
				F[5][ind+(Ny-1)*Nx]=oF[5][ind+(Ny-2)*Nx]+(SLight*dt-dx)/(SLight*dt+dx)*(F[5][ind+(Ny-2)*Nx]-oF[5][ind+(Ny-1)*Nx]);
				ind++;
			}
			ind+=Nxy-Nx;
		}
			

		printf("%d %e\n",it,it*dt);
		it++;
		
/* matlab output */

		if (matlab==1)
		{
			for (i=0;i<NumShow;i++)
			if (it % Show[i].when==0)
			{
			
/* show (any component) */

				if (Show[i].type==0)
				{
					if (Show[i].plane=='x')
					{
					  	for (z=0;z<Nz;z++)
							for (y=0;y<Ny;y++)
								Show[i].F[(y*Nz+z)]=F[Show[i].comp][(z*Ny+y)*Nx+Show[i].value];
						memcpy((void *)mxGetPr(Show[i].MF),(void *)Show[i].F, sizeof(double)*Ny*Nz);
					}

					if (Show[i].plane=='y')
					{
					  	for (z=0;z<Nz;z++)
							for (x=0;x<Nx;x++)
								Show[i].F[(x*Nz+z)]=F[Show[i].comp][(z*Ny+Show[i].value)*Nx+x];
						memcpy((void *)mxGetPr(Show[i].MF),(void *)Show[i].F, sizeof(double)*Nx*Nz);
					}
					if (Show[i].plane=='z')
					{
					  	for (x=0;x<Nx;x++)
							for (y=0;y<Ny;y++)
								Show[i].F[(x*Ny+y)]=F[Show[i].comp][(Show[i].value*Ny+y)*Nx+x];
						memcpy((void *)mxGetPr(Show[i].MF),(void *)Show[i].F, sizeof(double)*Nx*Ny);
					}
					engPutArray(ep,Show[i].MF);
					sprintf(str,"figure(%d);surf(F%d);axis([1 %d 1 %d]);shading interp;colorbar;",i+1,i,Nx,Ny);/*axis([1 %d 1 %d]);*/
					if (Show[i].movie==1) sprintf(str,"%s caxis([%e %e]);colorbar;",str,Show[i].c1,Show[i].c2);
					engEvalString(ep, str);
				}
			
/* show magnitude of E */

				if (Show[i].type==7)
				{
					if (Show[i].plane=='x')
					{
					  	for (z=0;z<Nz;z++)
							for (y=0;y<Ny;y++)
								Show[i].F[(y*Nz+z)]=log10(sqrt(
									F[3][(z*Ny+y)*Nx+Show[i].value]*F[3][(z*Ny+y)*Nx+Show[i].value]+
									F[4][(z*Ny+y)*Nx+Show[i].value]*F[4][(z*Ny+y)*Nx+Show[i].value]+
									F[5][(z*Ny+y)*Nx+Show[i].value]*F[5][(z*Ny+y)*Nx+Show[i].value]));
						memcpy((void *)mxGetPr(Show[i].MF),(void *)Show[i].F, sizeof(double)*Ny*Nz);
					}

					if (Show[i].plane=='y')
					{
					  	for (z=0;z<Nz;z++)
							for (x=0;x<Nx;x++)
								Show[i].F[(x*Nz+z)]=log10(sqrt(
									F[3][(z*Ny+Show[i].value)*Nx+x]*F[3][(z*Ny+Show[i].value)*Nx+x]+
									F[4][(z*Ny+Show[i].value)*Nx+x]*F[4][(z*Ny+Show[i].value)*Nx+x]+
									F[5][(z*Ny+Show[i].value)*Nx+x]*F[5][(z*Ny+Show[i].value)*Nx+x]));
						memcpy((void *)mxGetPr(Show[i].MF),(void *)Show[i].F, sizeof(double)*Nx*Nz);
					}
					if (Show[i].plane=='z')
					{
					  	for (x=0;x<Nx;x++)
							for (y=0;y<Ny;y++)
								Show[i].F[(x*Ny+y)]=log10(sqrt(
									F[3][(Show[i].value*Ny+y)*Nx+x]*F[3][(Show[i].value*Ny+y)*Nx+x]+
									F[4][(Show[i].value*Ny+y)*Nx+x]*F[4][(Show[i].value*Ny+y)*Nx+x]+
									F[5][(Show[i].value*Ny+y)*Nx+x]*F[5][(Show[i].value*Ny+y)*Nx+x]));
						memcpy((void *)mxGetPr(Show[i].MF),(void *)Show[i].F, sizeof(double)*Nx*Ny);
					}