/*	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);					}