fdtd.c

FDTD 可以用来创建mpeg

/*	FTDT - calculater - main program*/

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <signal.h>
#include <string.h>
#include <sys/kbd.h>
#include "engine.h"
#define  BUFSIZE 256

/*	electromagnetic constants */

#define e0 8.85e-12
#define m0 1.257e-6
#define SLight 3e8

/*	data for magnetic and electrical field 
	Hx,Hy,Hz,Ex,Ey,Ez in this order

	there is a 3D grid stored in a 1D array
	if you have a (x,y,z) point, the 1D representation
	in the array is (z*Ny+y)*Nx+x
*/

double *F[6];

/*	electrical field one timestep before
	these are needed to save the values for the
	absorbing boundary conditions
*/
	 
double *oF[6];

/*	epsilon and sigma fields (permitivity and conductivity) */

double *e,*s;

/*	coefficient fields for the fdtd formulation
	this is to calculate the coefficients before the run to save time
*/

double *cE1,*cE2,*cE3,cH;

/*	timestep in seconds */

double dt=0;

/*	final time of calculation */

double tend=1.;

/*	space discretization parameter in meters */

double dx=1;

/*	number of gridpoints in x,y and z direction */

long Nx=10,Ny=10,Nz=10;

/*	data for boxes of different values of epsilon and sigma */

struct 	{	long x1,x2,y1,y2,z1,z2;
			double sigma;
			double eps;
		} Box[100];
		
long NumBox=0;

/*	data for ramp source (voltage)
	x1,x2,y1,y2,z1,z2: the geometry of the box
	voltage : voltage-value of source
	direction : ... of the source; can be 'x','y' or 'z'
	timeon : time to grow to full value
	time : time the voltage stays until it drops back to zero
	R : inner  resistivity
*/

struct 	{	long x1,x2,y1,y2,z1,z2;
			double voltage;
			char direction;
			double timeon,time,R;
		} Volt;

/*	lumped resistors */

struct 	{	long x1,x2,y1,y2,z1,z2;
			char direction;
			double R;
		} R[100];
		
long NumR=0;

/*	lumped capacitors */

struct 	{	long x1,x2,y1,y2,z1,z2;
			char direction;
			double C;
		} C[100];
		
long NumC=0;

/*	cone of different material
	direction : of the axis
	x1,x2 : max in min on the axis
	r1,r2 : radius's of the circles at direction=x1,x2
	y,z : position af the axis
*/

struct 	{	long x1,x2,y,z,r1,r2;
			char direction;
			double eps,sigma;
		} Cone[100];
		
long NumCone=0;

/*	matlaboutput or not */

char matlab=0;

/*	data for matlab */

Engine *ep;

/*	matlab is on or not */

char matlabON=0;

/*	data for matlab outputs
	in each record there are the data for one output
	type : 0... for diferent types of output
	comp : for component output, can be 0..5 (Hx,Hy,Hz,Ex,Ey,Ez)
	plane : plane to show 'x','y','z'
	value : the plane to be shown is the plane=value plane (e.g. x=10)
	x1,..,z2 : different meanings
	when : after how many timesteps the output is updated
	F : 2D C-Array to be outputed
	MF : same array in matlab representation
	c1,c2 : values to fix the colormap during creating a movie
	num_pic : number of pictures for the mpg-movie
	Name : name of the mgp
	ind : actual number of picture
	quality : quality of movie 'h' or 'l'
*/	

struct {	char type;
			char comp;
			char plane;
			int value;
			int x1,x2,y1,y2,z1,z2;
			int when;
			double *F;
			mxArray *MF;

			int movie;
			double c1,c2;
			int num_pic;
			char Name[30];
			int ind;
			char quality;
		} Show[20];

long NumShow=0;

/*	Save whole old Field for displaying the displacement current*/

char Save_old=0;

/*	close matlab if ctrl-c or killed */

void abort_exception(int sig)
{
	if (matlabON==1)
	{
		engEvalString(ep, "close;");
		engClose(ep);
		printf("matlab closed.\n");
	}
	printf("\nfdtd aborted.\n");
	exit(0);
}

#include "parse.h"
#include "prepare.h"
#include "calc.h"




main(int argc, char **argv)
{
	char *DatafileName;
	char i;

/*	catch ctrl-c or kill to finish matlab correctly*/
	
	signal(SIGINT,abort_exception);
	signal(SIGTERM,abort_exception);
	
/*	read Inputfile*/

	if (argc==2)
	{
		DatafileName=argv[1];
		if (ReadData(DatafileName)==0)
		{
			
/* allocate the memory for all arrays */
			
			for (i=0;i<6;i++)
			{
				F [i]=(double *)malloc(sizeof(double)*Nx*Ny*Nz);
				oF[i]=(double *)malloc(sizeof(double)*Nx*Ny*Nz);
			}
			
			cE1=(double *)malloc(sizeof(double)*Nx*Ny*Nz);
			cE2=(double *)malloc(sizeof(double)*Nx*Ny*Nz);
			cE3=(double *)malloc(sizeof(double)*Nx*Ny*Nz);
			e  =(double *)malloc(sizeof(double)*Nx*Ny*Nz);
			s  =(double *)malloc(sizeof(double)*Nx*Ny*Nz);

/* call main calculation routine */

			prepare();
			calc();
		}
		else printf("there is something wrong with your datafile! \n");
	}
	else printf("\n'fdtd datafile' expected! \n");
}