fdtd-1d.txt

fdtd的一些一维到三维的程序,通常为pml边界,用txt文本编辑,无需解压密码

//************************************************************************
//    1D FDTD code modelling a pulse hitting a dielectric medium 
//    with absorbing boundary conditions at the left and right side
//***********************************************************************
//
//    Program author: Ren?Marklein
//                    Department of Electrical and Computer Engineering
//                    University of Kassel
//                    Wilhelmsh鰄er Allee 71
//                    D-34121 Kassel
//                    Tel. [phone] 
//                    [e-mail]
//
//    Date of this version:  October 31, 2001
//
//    This C program implements the finite-difference time-domain (FDTD)
//    solution of Maxwell's curl equations over a one-dimensional space
//    lattice comprised of uniform grid cells.
//
//    To illustrate the algorithm, a broadband raised cosine (RC) pulse 
//    propagating hitting an interface between vaccum (epsilon_r=1.0, mu_r=1.0) 
//    and a dielectric medium with a relative permittivity of epsilon_r=4.  
//
//    The cell size of Dx = lambda/20, G=20, is chosen to provide 20 grid 
//    cells per wavelength. The normalzied time step (Courant factor)
//    normDt= c*Dt/Dx is set to normDt=0.5.
//
//    The computational domain is truncated using a absorbing boundary 
//    condition.
//
//***********************************************************************

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<fcntl.h>
#include<io.h>
#include<sys/types.h>
#include<sys/stat.h>

#include<limits.h>

#ifndef M_PI
#define M_PI 3.1415F
#endif

#define N 200 // Number of cells in z-direction

int main()
{

	float *Ex, *Hy;  // Ex and Hy component
	float Jex;          // Jex component

	float *permittivity_r; // Permittivity
	float *beta;           // Inverse normalized permittivity on the dual grid

	float Ex_low_m1=0, Ex_low_m2=0;
	float Ex_high_m1=0, Ex_high_m2=0;

	int   nt,Nt;
	int   Nz;
	int   n;

	float normDt, Dt, Dt_ref;
	float normDz, Dz, Dx_ref;

	int   Nz_Source, N_Source;	

	float *pulse;
	
	float c_0, eps_0, mu_0;
	float c_ref, eps_ref, mu_ref;
	float f_0, f_ref, omega_0, omega_ref;
	float lambda_ref;
	float t;
	
	FILE *fout_log, *fout_Ex, *fout_Hy, *fout_pulse;
	FILE *fout_Ex_txt, *fout_Hy_txt, *fout_pulse_txt;
	
	Nz = N;

	c_0   = 2.99792458e8F;        // Speed of light in free space
	mu_0  = 4.0F*(float)M_PI*1.0e-7F;       // Permeability of free space
	eps_0 = 1.0F/(c_0*c_0*mu_0);  // Permittivty of free space

	c_ref = c_0;                 // Reference velocity
	eps_ref = eps_0;
	mu_ref = mu_0;

	f_0 = 10e9F;                  // Excitation frequency
	f_ref = f_0;                 // Reference frequency 

	omega_0 = 2.0F*M_PI*f_0;      // Excitation circular frequency
	omega_ref = omega_0;

	lambda_ref = c_ref/f_ref;    // Reference wavelength

	fprintf(stdout,"Reference wavelength: lambda_ref = %g\n", lambda_ref);

	Dz = lambda_ref/20;      // Reference cells width
	Dx_ref = Dz;
	normDz = Dz/Dx_ref;

	Dt_ref = Dx_ref/c_ref;       // Reference time step
	Dt = Dt_ref;

	fprintf(stdout,"Reference time step:      Dt_ref = %g\n", Dt_ref);

	normDt = 0.5;                    // Normalized time step

	Ex = calloc(sizeof(float),N);
	Hy = calloc(sizeof(float),N);
	permittivity_r = calloc(sizeof(float),N);
	beta = calloc(sizeof(float),N);

	// Initialize normalized field quantities		
	for(n=0; n<N; n++){
		Ex[n] = 0.0;
		Hy[n] = 0.0;
	}

	// Initialize normalized field quantities		
	for (n=0; n<N; n++){
		permittivity_r[n] = 1.0;		
	}
	for (n=N/2; n<N; n++){
		permittivity_r[n] = 4.0;
	}

	for (n=0; n<N; n++){
		beta[n] = normDt*2.0F/(permittivity_r[n] + permittivity_r[n+1]);
	}
	beta[N-1] = normDt/permittivity_r[N-1];		


	// Number of time steps
	printf("Number of time steps: Nt = ");
	scanf("%d",&Nt);
	printf("%d\n",Nt);

	pulse = calloc(sizeof(float),Nt);

	// Simulation time
	fprintf(stdout,"Simulation time: T = %g\n", normDt*Dt_ref*(float)Nt);
		
	// Source position
	fprintf(stdout,"Source position [%d-%d]: Nz_Source = ",0 ,Nz-1);
	scanf("%d",&Nz_Source);
	fprintf(stdout,"%d\n",Nz_Source);

	N_Source = Nz_Source;

	if ( ( fout_log = fopen("fdtd1dvac.log","w") ) == NULL ) { ; // Open log file
		fprintf(stdout,"Coun't open 'fdtd1dvac.log' output file\n");
		exit(1);
	}
	if ( ( fout_Ex    = fopen("Ex.bin","wb") ) == NULL ) {  // Ex component	
		fprintf(stdout,"Coun't open 'Ex.bin' output file\n");
		exit(1);
	}
	if ( ( fout_Hy    = fopen("Hy.bin","wb") ) == NULL ) {  // Hy component
		fprintf(stdout,"Coun't open 'Hy.bin' output file\n");
		exit(1);
	}
	if ( ( fout_pulse = fopen("Pulse.bin","wb") ) == NULL ) { // Excitation pulse
		fprintf(stdout,"Coun't open 'Pulse.bin' output file\n");
		exit(1);
	}
	if ( ( fout_Ex_txt    = fopen("Ex.txt","w") ) == NULL ) {  // Ex component	
		fprintf(stdout,"Coun't open 'Ex.txt' output file\n");
		exit(1);
	}
	if ( ( fout_Hy_txt    = fopen("Hy.txt","w") ) == NULL ) {  // Hy component
		fprintf(stdout,"Coun't open 'Hy.txt' output file\n");
		exit(1);
	}
	if ( ( fout_pulse_txt = fopen("Pulse.txt","w") ) == NULL ) { // Excitation pulse
		fprintf(stdout,"Coun't open 'Pulse.txt' output file\n");
		exit(1);
	}
	

	// Marching-on-in-time loop

	for (nt=0; nt<Nt; nt++) {

		t = Dt_ref*normDt*(float)nt; // Actual time

		// Compute Faraday's induction grid equation
		for( n=1; n<N; n++) {
			Hy[n] -= (float)normDt * ( Ex[n] - Ex[n-1] ); 
		}		

		// Initialize all Hy components outside the simulatin area to zero
		Hy[0] = 0.0;


		// Compute Ampere-Maxwell's grid equation			
		for (n=1; n<N-1; n++) {
			Ex[n] -= beta[n] * ( Hy[n+1] - Hy[n] );
		}

		// Apply open boundary condition
		Ex[0]     = Ex_low_m2;
		Ex_low_m2 = Ex_low_m1;
		Ex_low_m1 = Ex[1];

		Ex[N-1]    = Ex_high_m2;
		Ex_high_m2 = Ex_high_m1;
		Ex_high_m1 = Ex[N-2];
	
		// Apply excitation pulse
		// Raised cosine pulse with two cycles - (RC2) pulse
		if ( t < 2.0*2.0*M_PI/(omega_0) ) 
			pulse[nt] = 0.5*(float)((1.0-cos((double)(omega_0*t/2.0F)))*cos((double)(omega_0*t)));
		else 
			pulse[nt] = 0.0;
			
		// Apply excittion pulse as an impress electric current density
		Jex = pulse[nt];

        // Apply excittion pulse as an impress electric current density
		Ex[N_Source] -= (float)normDt * Jex;
		
		fprintf(fout_log,"t= %8.4f pulse= %8.4f  Ex[N_Source]= %8.4f\n", 
			t, pulse[nt], Ex[N_Source]);

		// Write output files
		fwrite(&Ex[0],sizeof(float),N,fout_Ex);		
		fwrite(&Hy[0],sizeof(float),N,fout_Hy);	
		
		for (n=0; n<N; n++) {
	      fprintf(fout_Ex_txt,"%10.5f %10.5f \n",(float)n,Ex[n]);
		}

		fwrite(&pulse[nt],sizeof(float),1,fout_pulse);

  	    fprintf(fout_pulse_txt,"%10.5f %10.5f \n",(float)nt,pulse[nt]);

	} // End of marching-on-in-time loop

	for (n=0; n<N; n++) {
	  fprintf(fout_Ex_txt,"%10.5f %10.5f \n",(float)n,Ex[n]);
	}

	for (n=0; n<N; n++) {
	  fprintf(fout_Hy_txt,"%10.5f %10.5f \n",(float)n,Hy[n]);
	}
	
	// Close all output files
	fclose(fout_Ex);
	fclose(fout_Hy);
	fclose(fout_pulse);
	
	return(0);
}