FDTD电磁波算法的一些问题

/* FDTD-1D-2.2.c  1D FDTD simulation of a lossy dielectric medium  */
/* Simulation of a sinusoidal wave or a Gauss pulse transmiting in a lossy dielectric medium exsiting in free space */ 

/*8888888888888  Using a new formulation and using flux density  8888888888888*/

/*    The Fourier Transform has been added        */

/*    1D FDTD simulation of a frequency dependent material    */

# include <math.h>
# include <stdlib.h>
# include <stdio.h>
# define KE 200         /* KE is the number of cells to be used   */
# define NF 10          /* NF is the number of the frequence in fourier transform */
# define Pi 3.1415926
void main ()
{
   
	double ex_low_m1, ex_low_m2,ex_high_m1,ex_high_m2;
	
	double dx[KE],ix[KE],ex[KE],hy[KE],sx[KE];
	double ga[KE],gb[KE],gc[KE];

	int n, m,k, kc,NSTEPS;
	int kstart_dielectric;                                        /* 介质的左边界*/
	int kend_dielectric;                                          /* 介质的右边界*/ 
	double T;
	double t0, spread,pulse, pulse1, freq_in,  pulse2;        /*pulse1 is Gauss pulse, Pulse2 is a sinusoidal wave */
	double epsilon_r,epsilon_0;                                    /* 介质的相对介电常数 和 真空介电常数*/   

	double sigma;                                                  /* 介质的电导率 */
    double ddx, dt;                                                /* The FDTD Cell size and the time interval */
    
	double freq[NF], arg[NF], ampn[NF][KE],phasen[NF][KE];
	double real_part[NF][KE], imag_part[NF][KE]; 
	double real_in[NF], imag_in[NF], amp_in[NF], phase_in[NF];
	double mag[KE] ;
	double tau, chil, del_exp;


	FILE *fp ;


	
	/* Initialize  */
	
	ddx=0.01;                                      /* Set the cell size to 1cm */ 
    dt=ddx/(2*3e8) ;                               /* calculate the time step */
    epsilon_0=8.85419e-12;
    
	
	
     
	kc=5;           /* 激励源所引入的位置 */
    t0= 50.0;          /* Center of the incident pulse */
	spread=10.0;         /* Width of the incident pulse  */
	T=0;
	NSTEPS=1;

	ex_low_m1=0;       /* 在n+1时存储左边界处前一个格点ex[1] 结果 */
	ex_low_m2=0;       /* 在n+2时存储左边界处前一个格点ex[1] 结果 */
	ex_high_m1=0;      /* 在n+1时存储右边界处前一个格点ex[KE-2] 结果 */
	ex_high_m2=0;      /* 在n+2时存储右边界处前一个格点ex[KE-2] 结果 */
	
	
	
	
	for (k=0; k<KE; k++)
	{    ga[k]=1.0;                               /* initialize to free space */
		 gb[k]=0.;
		 gc[k]=0;                                 /* initialize to free space */
		 ex[k]=0.;
	     hy[k]=0.; 
		 dx[k]=0.;
		 ix[k]=0;
		 sx[k]=0;
		 mag[k]=0;
	

        for (m=0;m<=2;m++)
		{  real_part[m][k]=0.;           /* Real and imaginary parts of fourier transform */
	       imag_part[m][k]=0.;
	       ampn[m][k]=0.;                /* Amplituds and phase of the fourier transform */
	       phasen[m][k]=0.;
	   }
	}

	for (m=0;m<=2;m++)
	{   real_in[m]=0.;             /* Fourier Trans. of input pulse */
	    imag_in[m]=0.;
	}

/* Parameters for the Fourier Transform  */
	
	freq[0]=50.e6;
	freq[1]=200.e6;
	freq[2]=500.e6;

	
    for (m=0;m<=2;m++)
	{    arg[m]=2*Pi*freq[m]*dt;
	printf("%2d %6.2f MHz %7.5f \n ", m, freq[m]*1e-6,arg[m]);
	}

    
	/* These parameters specify the input pulse */
//	 printf("input frequence (MHz) -->");
//     scanf("%lf", &freq_in);
//     freq_in=freq_in*1e6;
//	 printf("%8.0f \n",freq_in);

	printf("Dielectric start at  kstart_dielectric -->");
     scanf("%d", &kstart_dielectric);

	 printf("Dielectric end at  kend_dielectric -->");
     scanf("%d", &kend_dielectric);

     printf("Epsilon_r -->");
     scanf("%lf", &epsilon_r);  
	 
	 printf("Conductivity -->");
     scanf("%lf", &sigma);

	 printf("chil----->");
	 scanf("%lf", &chil);

     

	 tau=1000.;                   /* Make sure tau is >0.    */
     if (chil>0.0001)
	 {printf("tau (in microseconds)--->"); scanf("%lf",&tau); del_exp=exp(-dt/tau);}
	 
	
    printf ("%d %d %f  %f %f %f \n", kstart_dielectric,kend_dielectric,epsilon_r, sigma,tau,chil);
     tau=1.e-6*tau;
	 
	  printf("del_exp= %f \n",del_exp) ;         


     for (k=kstart_dielectric; k<=kend_dielectric;k++)
	 {   ga[k]=1/(epsilon_r+sigma*dt/epsilon_0+chil*dt/tau);
		 gb[k]=sigma*dt/epsilon_0;               /* initialize to dielectric medium */
		 gc[k]=chil*dt/tau;
	 }
     
	 for (k=0;k<KE;k++)
	 { printf("%2d %4.2f %4.2f \n",k,ga[k],gb[k]); }
	
	
	 
	 /* Main part of program */
	 
	 while (NSTEPS>0) 
	{	printf("NSTEPS -->");     /* NSTEPS is the number of times that the Main loop has executed */
		scanf("%d", &NSTEPS);  
		printf ("%d \n", NSTEPS);
	
 printf("程序正在运行，请稍侯！ \n");
   
 for (n=1; n<=NSTEPS; n++)
	{
		T=T+1;                   /* T keeps track of the total number of times the main loop is executed */

	/* Main FDTD Loop  */

		/* Calculate the Dx field   */

		for (k=1; k<KE; k++)
		{ dx[k]=dx[k]+ 0.5*(hy[k-1]-hy[k]); }


		/* Put a sinusoidal wave or the Gauss pulse at the cell kc */

	    pulse1=exp(-0.5*(pow((t0-T)/spread,2.0)));
//		pulse2=sin(2*Pi*freq_in*dt*T);
		
		pulse=pulse1;                                            /* 脉冲的类型由pulse1和pulse2两个变量来选择 */
		
		dx[kc]=dx[kc]+pulse;
		printf("%5.1f %f %6.2f \n", T,pulse,dx[kc]);
        
     /* Calculate the Ex from the Dx */
	
		for (k=0;k<KE;k++)
		{     ex[k]=ga[k]*(dx[k]-ix[k]-sx[k]);
		      ix[k]=ix[k]+gb[k]*ex[k];
			  sx[k]=del_exp*sx[k]+gc[k]*ex[k];
		}
		
		/* Calculate the fourier trasform of Ex */

		for (k=0;k<KE;k++)
		{   for (m=0;m<=2;m++)
			{   real_part[m][k]=real_part[m][k]+cos(arg[m]*T)*ex[k];
			    imag_part[m][k]=imag_part[m][k]-sin(arg[m]*T)*ex[k];
			}
		}

		/* Fourier Transform of the input pulse */
		if (T<100)
		{
		      for (m=0;m<=2;m++)
			  {   real_in[m]=real_in[m]+cos(arg[m]*T)*ex[10];
			      imag_in[m]=imag_in[m]-sin(arg[m]*T)*ex[10];
			  }
		}
		      
		
		/* Absorbing boundary conditions */

		ex[0]=ex_low_m2;
		ex_low_m2=ex_low_m1;
		ex_low_m1=ex[1];


		ex[KE-1]=ex_high_m2;
		ex_high_m2=ex_high_m1;
		ex_high_m1=ex[KE-2];               
		


		/* Calculate the hy field */
		for (k=0; k<KE-1;k++)
		{ hy[k]=hy[k]+0.5*(ex[k]-ex[k+1]); }		
	}
	/* End of the Main FDTD Loop */
	    
	    /* At the end of the calculation, print out the Ex and Hy field */
	    for (k=0; k<KE; k++)
		{    printf("%3d %6.2f %6.2f \n", k, ex[k], hy[k]);
		}
        

		/* Write the E field out to a file "Ex"  */

		fp=fopen("Ex.txt","w");
		fprintf(fp, "T=%5.0f \n", T); 
		for (k=0;k<KE;k++)
		{  
			fprintf (fp," %5d %6.2f %6.2f \n", k,dx[k],ex[k]);
		}

		fclose(fp);


		/* Write the H field out to a file "Hy" */
		fp=fopen("Hy.txt","w");
		
		 fprintf(fp, "T=%5.0f \n", T); 
		for (k=0;k<KE;k++)
		{  
			fprintf (fp,"%5d  %6.2f \n",k, hy[k]);
		}

		fclose(fp);

        printf(" %5.0f \n", T);
        
		/* calculate the amplitude and phase of each frequency */

		/* Amplitude and the phase of the input pulse */

		for (m=0; m<=2;m++)
		{   amp_in[m]=sqrt(pow(imag_in[m],2.0)+pow(real_in[m],2.0));
		    phase_in[m]=atan2(imag_in[m],real_in[m]);
			printf("%d Input pulse: %8.4f %8.4f %8.4f %7.2f \n", m,real_in[m],imag_in[m],amp_in[m],(180.0/Pi)*phase_in[m]);
			for(k=0;k<KE;k++)
			{   ampn[m][k]=(1.0/amp_in[m]*sqrt(pow(real_part[m][k],2.0)+pow(imag_part[m][k],2.0)));
			    phasen[m][k]=atan2(imag_part[m][k],real_part[m][k])-phase_in[m];
			}
		}
			/* write the amplitude field out to a files "amp "  */
			fp=fopen ("amp0.txt","w");
			for (k=0;k<KE;k++)
			{ fprintf(fp,"%d %8.5f \n ", k,ampn[0][k] ); }
			fclose(fp);
			fp=fopen ("amp1.txt","w");
			for (k=0;k<KE;k++)
			{ fprintf(fp,"%d %8.5f \n ", k,ampn[1][k] ); }
			fclose(fp);
			fp=fopen ("amp2.txt","w");
			for (k=0;k<KE;k++)
			{ fprintf(fp,"%d %8.5f \n ", k,ampn[2][k] ); }
			fclose(fp);
		

		 

}
printf("程序运行完毕 \n");
}

运行之后出现如下错误

root@gallup-virtual-machine:/home/source# gcc -Wall  li.c -o li 

li.c:7:6: 警告： ‘main’的返回类型不是‘int’ [-Wmain]
li.c: 在函数‘main’中:
li.c:28:9: 警告： 变量‘mag’被设定但未被使用 [-Wunused-but-set-variable]
li.c:25:41: 警告： 变量‘phasen’被设定但未被使用 [-Wunused-but-set-variable]
li.c:19:45: 警告： 未使用的变量‘pulse2’ [-Wunused-variable]
li.c:19:35: 警告： 未使用的变量‘freq_in’ [-Wunused-variable]
/tmp/ccwrACyT.o: In function `main':
li.c:(.text+0x3fb): undefined reference to `exp'
li.c:(.text+0x6a3): undefined reference to `exp'
li.c:(.text+0x860): undefined reference to `cos'
li.c:(.text+0x8cb): undefined reference to `sin'
li.c:(.text+0x986): undefined reference to `cos'
li.c:(.text+0x9d0): undefined reference to `sin'
li.c:(.text+0xcfe): undefined reference to `sqrt'
li.c:(.text+0xd3a): undefined reference to `atan2'
li.c:(.text+0xe1e): undefined reference to `sqrt'
li.c:(.text+0xe7e): undefined reference to `atan2'
collect2: ld 返回 1

root@gallup-virtual-machine:/home/source# 

解决方法：

1、 root@gallup-virtual-machine:/home/source# gcc   li.c  -lm -o li