fdtd-2d.txt

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

//************************************************************************
//    Electromagnetic 2D TM FDTD code with PEC boundary conditions
//***********************************************************************
//    C Program author: 
//    Dr.-Ing. Ren?Marklein
//    University of Kassel, Department of Electrical and 
//    Computer Engineering, Electromagnetic Theory
//    Wilhelmsh鰄er Allee 71, D-34121 Kassel, Germany
//    Tel. [phone], [e-mail]
//
//    Date of this version:  December 19, 2001
//
//    This C program implements the finite-difference time-domain 
//    (FDTD) solution of Maxwell's eqautions in differential form over a 
//    two-dimensional space in TM_y case using a uniform grid complex.
//    To illustrate the algorithm, a broadband raised cosine (RC) pulse 
//    with a center frequency of f_0 = 10 GHz propagating in vaccum 
//    (epsilon_r=1.0, mu_r=1.0) is modeled.  
//    The cell size of Dx = lambda_0/20, G = 20, is chosen to provide 20 
//    grid cells per wavelength at center frequency. The normalzied time 
//    step (Courant factor) normDt = c_0*Dt/Dx is set to normDt = 0.5.
//    The computational domain is truncated using the perfectly 
//    electrically conducting (PEC) boundary condition.
//
//***********************************************************************

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<malloc.h>

#ifndef M_PI
#define M_PI 3.1415F
#endif

int main()
{
    float *Hx, *Hz, *Ey;  // Hx, Hz, and Ey component
    float Jey;            // Jey component
    int   nt,Nt;
    int   Nx, Nz;
	int   Mx, Mz;
    int   n, N;
	int   nx, nz;

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

    int   Nx_Source, 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;

	int   delta_snapshot, generate_snapshot;
    
    FILE *fout_log, *fout_Hx, *fout_Hz, *fout_Ey, *fout_pulse;
        
    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

    // Number of cells in z-direction
	Nx = Nz = 0;
	while ( Nx < 1 ) {
		printf("Number of cells in x-direction: Nx = ");
		scanf("%d",&Nx);
	}
	while ( Nz < 1 ) {
		printf("Number of cells in z-direction: Nz = ");
		scanf("%d",&Nz);
	}

	Mx = 1;
	Mz = Nx;
    N  = Nx*Nz;

	delta_snapshot    = 10;
	generate_snapshot = 0;
    
    Hx = calloc(sizeof(float),N);
	Hz = calloc(sizeof(float),N);
    Ey = calloc(sizeof(float),N);

    // Initialize field quantities      
    for(n=0; n<N; n++){
        Hx[n] = 0.0;
		Hz[n] = 0.0;
        Ey[n] = 0.0;
    }

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

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

    // Time intervall to simulate
    fprintf(stdout,"Simulation time: T = %g\n", normDt*Dt_ref*(float)Nt);
        
    // Source position
    fprintf(stdout,"Source position [%d-%d]: Nx_Source = ",2 ,Nx-3);
    scanf("%d",&Nx_Source);

	fprintf(stdout,"Source position [%d-%d]: Nz_Source = ",2 ,Nz-3);
    scanf("%d",&Nz_Source);

    N_Source = Nx_Source + Nz_Source * Nx;

    // Open log file
    if ( ( fout_log = fopen("fdtd2dtmvac.log","w") ) == NULL ) { ; 
        fprintf(stdout,"Coun't open 'fdtd1dvac.log' output file\n");
        exit(1);
    }
    // Open binary output files
    if ( ( fout_Hx    = fopen("Hx.bin","wb") ) == NULL ) {  // Ex component 
        fprintf(stdout,"Coun't open 'Hx.bin' output file\n");
        exit(1);
    }
	if ( ( fout_Hz    = fopen("Hz.bin","wb") ) == NULL ) {  // Ex component 
        fprintf(stdout,"Coun't open 'Hz.bin' output file\n");
        exit(1);
    }
    if ( ( fout_Ey    = fopen("Ey.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);
    }
   
    // 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+Mz; n<N; n++) {
            Hx[n] -= (float)normDt * ( Ey[n-Mz] - Ey[n   ] );     
            Hz[n] -= (float)normDt * ( Ey[n]    - Ey[n-Mx] ); 
		}
              
        // Initialize all Hx and Hz components outside the simulatin area to zero
		// Upper boundary
		for (nx=0; nx<Nx; nx++) {
			n = nx;
		    Hx[n] = 0.0F;
		}
		// Lower boundary
		for (nx=0; nx<Nx; nx++) {
			n = nx + (Nz-1)*Nx;
		    Hx[n] = 0.0F;
			Hz[n] = 0.0F;
		}
		// Left boundary
		for (nz=0; nz<Nz; nz++) {
			n = nz*Nx;
		    Hz[n] = 0.0F;
		}
		// Right boundary
		for (nz=0; nz<Nz; nz++) {
			n = (Nx-1) + nz*Nx;
		    Hx[n] = 0.0F;
			Hz[n] = 0.0F;
		}
        
        // Compute Ampere-Maxwell's grid equation           
        for (n=1; n<N-Mz; n++) {
		    Ey[n] += (float)normDt * ( Hx[n+Mz] - Hx[n] - Hz[n+Mx] + Hz[n] );            
		}
        
        // Apply PEC boundary condition
		// Upper boundary
		for (nx=0; nx<Nx; nx++) {
			n = nx;
		    Ey[n] = 0.0F;
		}
		// Lower boundary
		for (nx=0; nx<Nx; nx++) {
			n = nx + (Nz-1)*Nx;
		    Ey[n] = 0.0F;
			n = nx + (Nz-2)*Nx;
		    Ey[n] = 0.0F;
		}
		// Left boundary
		for (nz=0; nz<Nz; nz++) {
			n = nz*Nx;
		    Ey[n] = 0.0F;
		}
		// Right boundary
		for (nz=0; nz<Nz; nz++) {
			n = (Nx-1) + nz*Nx;
		    Ey[n] = 0.0F;
			n = (Nx-2) + nz*Nx;
		    Ey[n] = 0.0F;
		}
      
    
        // 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
        Jey = pulse[nt];

        // Apply excittion pulse as an impress electric current density
		Ey[N_Source] -= (float)normDt * Jey;
	       
        printf("t= %8.4g pulse= %8.4g  Ey[N_Source]= %8.4g\n", t, pulse[nt], Ey[N_Source]);

		++generate_snapshot;

		if ( generate_snapshot == delta_snapshot ) {
			
			// Write binary output files
			fwrite(&Hx[0],sizeof(float),N,fout_Hx);     
			fwrite(&Hz[0],sizeof(float),N,fout_Hz);     
			fwrite(&Ey[0],sizeof(float),N,fout_Ey);
			fwrite(&pulse[nt],sizeof(float),1,fout_pulse);      

			generate_snapshot = 0;
		}
		

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

    // Close all output files
    fclose(fout_Hx);
	fclose(fout_Hz);
    fclose(fout_Ey);
    fclose(fout_pulse);
   
    return(0);
}