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