calc.h

FDTD 可以用来创建mpeg

					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 H */

				if (Show[i].type==8)
				{
					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[0][(z*Ny+y)*Nx+Show[i].value]*F[0][(z*Ny+y)*Nx+Show[i].value]+
									F[1][(z*Ny+y)*Nx+Show[i].value]*F[1][(z*Ny+y)*Nx+Show[i].value]+
									F[2][(z*Ny+y)*Nx+Show[i].value]*F[2][(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[0][(z*Ny+Show[i].value)*Nx+x]*F[0][(z*Ny+Show[i].value)*Nx+x]+
									F[1][(z*Ny+Show[i].value)*Nx+x]*F[1][(z*Ny+Show[i].value)*Nx+x]+
									F[2][(z*Ny+Show[i].value)*Nx+x]*F[2][(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[0][(Show[i].value*Ny+y)*Nx+x]*F[0][(Show[i].value*Ny+y)*Nx+x]+
									F[1][(Show[i].value*Ny+y)*Nx+x]*F[1][(Show[i].value*Ny+y)*Nx+x]+
									F[2][(Show[i].value*Ny+y)*Nx+x]*F[2][(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 current */

				if (Show[i].type==1)
				{
					if (Show[i].plane=='y')
					{
					  	for (z=0;z<Nz;z++)
							for (x=0;x<Nx;x++)
								if ((Show[i].value<Volt.y2+4) && (Show[i].value>Volt.y1-4) &&
									(x<Volt.x2+4) && (x>Volt.x1-4) &&				
									(z<Volt.z2+4) && (z>Volt.z1-4)) Show[i].F[(x*Nz+z)]=0;
								else
								Show[i].F[(x*Nz+z)]=(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])
														*s[(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=='x')
					{
					  	for (z=0;z<Nz;z++)
							for (y=0;y<Ny;y++)
								if ((Show[i].value<Volt.x2+4) && (Show[i].value>Volt.x1-4) &&
									(y<Volt.y2+4) && (y>Volt.y1-4) &&
									(z<Volt.z2+4) && (z>Volt.z1-4)) Show[i].F[(y*Nz+z)]=0;
								else
								Show[i].F[(y*Nz+z)]=(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])
														*s[(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=='z')
					{
					  	for (x=0;x<Nx;x++)
							for (y=0;y<Ny;y++)
								if ((Show[i].value<Volt.z2+4) && (Show[i].value>Volt.z1-4) &&
									(x<Volt.x2+4) && (x>Volt.x1-4) &&
									(y<Volt.y2+4) && (y>Volt.y1-4)) Show[i].F[(x*Nz+y)]=0;
								else
								Show[i].F[(x*Ny+y)]=(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])
														*s[(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);shading interp;axis([1 %d 1 %d]);colorbar;",i+1,i,Nx,Ny);
					if (Show[i].movie==1) sprintf(str,"%s caxis([%e %e]);colorbar;",str,Show[i].c1,Show[i].c2);
					engEvalString(ep, str);

				}
				
/* show voltage in a plane */

				if (Show[i].type==3)
				{
					for (t=Show[i].y1;t<=Show[i].y2;t++)
					{
						Show[i].F[t-Show[i].y1]=0;
						for (k=Show[i].x1;k<=Show[i].x2;k++)
						{
							if ((Show[i].plane=='x')&&(Show[i].comp==4)) Show[i].F[t-Show[i].y1]+=F[4][(t*Ny+k)*Nx+Show[i].value];
							if ((Show[i].plane=='x')&&(Show[i].comp==5)) Show[i].F[t-Show[i].y1]+=F[5][(k*Ny+t)*Nx+Show[i].value];
							if ((Show[i].plane=='y')&&(Show[i].comp==3)) Show[i].F[t-Show[i].y1]+=F[3][(t*Ny+Show[i].value)*Nx+k];
							if ((Show[i].plane=='y')&&(Show[i].comp==5)) Show[i].F[t-Show[i].y1]+=F[5][(k*Ny+Show[i].value)*Nx+t];
							if ((Show[i].plane=='z')&&(Show[i].comp==3)) Show[i].F[t-Show[i].y1]+=F[3][(Show[i].value*Ny+t)*Nx+k];
							if ((Show[i].plane=='z')&&(Show[i].comp==4)) Show[i].F[t-Show[i].y1]+=F[4][(Show[i].value*Ny+k)*Nx+t];
						}
						Show[i].F[t-Show[i].y1]*=dx;
						
					}
					memcpy(mxGetPr(Show[i].MF),Show[i].F, sizeof(double)*(Show[i].y2-Show[i].y1+1));
					engPutArray(ep,Show[i].MF);
					sprintf(str,"figure(%d);plot(F%d);",i+1,i);
					if (Show[i].movie==1) sprintf(str,"%s axis([%d %d %e %e]);colorbar;",str,Show[i].y1,Show[i].y2,Show[i].c1,Show[i].c2);
					engEvalString(ep, str);
				}
				
/* show voltage between two points */

				if (Show[i].type==4)
				{
					V=0;
					
					if (Show[i].x1!=Show[i].x2)
						if (Show[i].x1<Show[i].x2)
							for (k=Show[i].x1;k<=Show[i].x2;k++)
								V+=F[3][(Show[i].z1*Ny+Show[i].y1)*Nx+k];
						else
							for (k=Show[i].x2;k<=Show[i].x1;k++)
								V-=F[3][(Show[i].z1*Ny+Show[i].y1)*Nx+k];

					if (Show[i].y1!=Show[i].y2)
						if (Show[i].y1<Show[i].y2)
							for (k=Show[i].y1;k<=Show[i].y2;k++)
								V+=F[4][(Show[i].z1*Ny+k)*Nx+Show[i].x2];
						else
							for (k=Show[i].y2;k<=Show[i].y1;k++)
								V-=F[4][(Show[i].z1*Ny+k)*Nx+Show[i].x2];

					if (Show[i].z1!=Show[i].z2)
						if (Show[i].z1<Show[i].z2)
							for (k=Show[i].z1;k<=Show[i].z2;k++)
								V+=F[5][(k*Ny+Show[i].y2)*Nx+Show[i].x2];
						else
							for (k=Show[i].z2;k<=Show[i].z1;k++)
								V-=F[5][(k*Ny+Show[i].y2)*Nx+Show[i].x2];
								
					V*=dx;
					sprintf(str,"V%d=[V%d(1,2) %e];t%d=[t%d(1,2) %e];figure(%d);plot(t%d,V%d);",i,i,V,i,i,it*dt,i+1,i,i);
					engEvalString(ep, str);
				}
				
/* show displacement curents */

				if (Show[i].type==2)
				{
					if (Show[i].plane=='x')
					{
					  	for (z=0;z<Nz;z++)
							for (y=0;y<Ny;y++)
								if ((Show[i].value<Volt.x2+4) && (Show[i].value>Volt.x1-4) &&
									(y<Volt.y2+4) && (y>Volt.y1-4) &&
									(z<Volt.z2+4) && (z>Volt.z1-4)) Show[i].F[(y*Nz+z)]=0;
								else
								Show[i].F[(y*Nz+z)]=sqrt(	(F[3][(z*Ny+y)*Nx+Show[i].value]-oF[3][(z*Ny+y)*Nx+Show[i].value])*
															(F[3][(z*Ny+y)*Nx+Show[i].value]-oF[3][(z*Ny+y)*Nx+Show[i].value])+
															(F[4][(z*Ny+y)*Nx+Show[i].value]-oF[4][(z*Ny+y)*Nx+Show[i].value])*
															(F[4][(z*Ny+y)*Nx+Show[i].value]-oF[4][(z*Ny+y)*Nx+Show[i].value])+
															(F[5][(z*Ny+y)*Nx+Show[i].value]-oF[5][(z*Ny+y)*Nx+Show[i].value])*
															(F[5][(z*Ny+y)*Nx+Show[i].value]-oF[5][(z*Ny+y)*Nx+Show[i].value]))
															*e[(z*Ny+y)*Nx+Show[i].value]/dt;
						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++)
								if ((Show[i].value<Volt.y2+4) && (Show[i].value>Volt.y1-4) &&
									(x<Volt.x2+4) && (x>Volt.x1-4) &&
									(z<Volt.z2+4) && (z>Volt.z1-4)) Show[i].F[(x*Nz+z)]=0;
								else
								Show[i].F[(x*Nz+z)]=sqrt(	(F[3][(z*Ny+Show[i].value)*Nx+x]-oF[3][(z*Ny+Show[i].value)*Nx+x])*
															(F[3][(z*Ny+Show[i].value)*Nx+x]-oF[3][(z*Ny+Show[i].value)*Nx+x])+
															(F[4][(z*Ny+Show[i].value)*Nx+x]-oF[4][(z*Ny+Show[i].value)*Nx+x])*
															(F[4][(z*Ny+Show[i].value)*Nx+x]-oF[4][(z*Ny+Show[i].value)*Nx+x])+
															(F[5][(z*Ny+Show[i].value)*Nx+x]-oF[5][(z*Ny+Show[i].value)*Nx+x])*
															(F[5][(z*Ny+Show[i].value)*Nx+x]-oF[5][(z*Ny+Show[i].value)*Nx+x]))
															*e[(z*Ny+Show[i].value)*Nx+x]/dt;
						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++)
								if ((Show[i].value<Volt.z2+4) && (Show[i].value>Volt.z1-4) &&
									(x<Volt.x2+4) && (x>Volt.x1-4) &&
									(y<Volt.y2+4) && (y>Volt.y1-4)) Show[i].F[(x*Ny+y)]=0;
								else	
								Show[i].F[(x*Ny+y)]=sqrt(	(F[3][(Show[i].value*Ny+y)*Nx+x]-oF[3][(Show[i].value*Ny+y)*Nx+x])*
															(F[3][(Show[i].value*Ny+y)*Nx+x]-oF[3][(Show[i].value*Ny+y)*Nx+x])+
															(F[4][(Show[i].value*Ny+y)*Nx+x]-oF[4][(Show[i].value*Ny+y)*Nx+x])*
															(F[4][(Show[i].value*Ny+y)*Nx+x]-oF[4][(Show[i].value*Ny+y)*Nx+x])+
															(F[5][(Show[i].value*Ny+y)*Nx+x]-oF[5][(Show[i].value*Ny+y)*Nx+x])*
															(F[5][(Show[i].value*Ny+y)*Nx+x]-oF[5][(Show[i].value*Ny+y)*Nx+x]))
															*e[(Show[i].value*Ny+y)*Nx+x]/dt;
						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);shading interp;axis([1 %d 1 %d]);colorbar;",i+1,i,Nx,Ny);
					if (Show[i].movie==1) sprintf(str,"%s caxis([%e %e]);colorbar;",str,Show[i].c1,Show[i].c2);
					engEvalString(ep, str);

				}

				if ((Show[i].movie==1)&&(Show[i].ind<=Show[i].num_pic))
				{
					if (Show[i].ind==1)
					{
						sprintf(str,"M%d=moviein(%d);",i,Show[i].num_pic);
						engEvalString(ep, str);
					}
					sprintf(str,"M%d(:,%d)=getframe;",i,Show[i].ind);						
					engEvalString(ep, str);
					Show[i].ind++;
					if (Show[i].ind==Show[i].num_pic+1)
					{
					
						if (Show[i].quality=='h')
							sprintf(str,"mpgwrite(M%d,jet,'%s',[1 2 1 1])",i,Show[i].Name);
						else
							sprintf(str,"mpgwrite(M%d,jet,'%s',[])",i,Show[i].Name);
						
						engEvalString(ep, str);
					}
				}
			}
		}
	}

/* close matlab */

	if (matlab==1)
	{
		for (i=0;i<NumShow;i++)
		{
			sprintf(str,"F%d",i);
			mxDestroyArray(Show[i].MF);
		}
		
		engEvalString(ep, "close;");
		engClose(ep);
		matlabON=0;
	}
}