Problems from Sullivan's FDTD 3D codes

/*File: FD3D.h*/

#include "iostream.h"
#include "fstream.h"
#include "math.h"
#include "stdlib.h"
#include "malloc.h"

#define IE 40
#define JE 40
#define KE 40
#define N 40


class ThreeDFDTD
{public: 	/*float gax[IE][JE][KE],***[IE][JE][KE],gaz[IE][JE][KE];	float dx[IE][JE][KE],dy[IE][JE][KE],dz[IE][JE][KE];	float ex[IE][JE][KE],ey[IE][JE][KE],ez[IE][JE][KE];	float hx[IE][JE][KE],hy[IE][JE][KE],hz[IE][JE][KE];*/	float ***gax,******,***gaz;	float ***dx,***dy,***dz;	float ***ex,***ey,***ez;	float ***hx,***hy,***hz;	float ddx,dt,T,epsz,epsilon,sigma,eaf;	float xn,xxn,xnum,xd,curl_e;	float t0,spread,pulse,NSTEPs;	int ic,jc,kc;	void createMang(float ***&mang);	void freeMang(float ***&mang);	void initialize(void);	void fdtd(void);
};

/*File: Main.cpp*/

#include "FD3D.h"

int main()
{int mode;ThreeDFDTD c;
while(1){	cout << endl		 << "****3 Dimensional FDTD simulation*****" 		 << "\n****1> 3D FDTD. Dipole in free space"   		 << "\n****2> Exit" 	     << endl;	cin  >> mode;	switch(mode)	{	case 1:		{			c.fdtd();		break;		}	default:exit(0);	}}return 0;
}

/*File FD3D.cpp*/

#include "FD3D.h"

void ThreeDFDTD::createMang(float ***&mang)
{
      /* mang = new double**[N];
for (int i=0; i<N; i++)
    mang[i] = new double*[N];

for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
    mang[i][j] = new double[N];*/
// gan gia tri
/*for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
for (int k=0; k<N; k++)
    mang[i][j][k] = i+j+k;*/
    /*mang=(float***)calloc(IE*JE*KE,sizeof(float**));
    for(int i=0;i<N;i++)
    {
            mang[i]=(float**)calloc(IE*JE,sizeof(float*));
            for(int j=0;j<N;j++)
            {
                 mang[i][j]=(float*)calloc(IE,sizeof(float));
            }
    }*/
    mang=(float***)malloc(IE*JE*KE*sizeof(float**));
    for(int i=0;i<N;i++)
    {
            mang[i]=(float**)malloc(IE*JE*sizeof(float*));
            for(int j=0;j<N;j++)
            {
                 mang[i][j]=(float*)malloc(IE*sizeof(float));
            }
    }   
}

void ThreeDFDTD::freeMang(float ***&mang)
{
// giai phong mang
/*for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
    delete mang[i][j];
    
for (int i=0; i<N; i++)
    delete mang[i];

     delete mang;*/
     for(int i=0;i<N;i++)
     {
          for(int j=0;j<N;j++)
              free(mang[i][j]);
          free(mang[i]);
     }
     free(mang);
     
}

void ThreeDFDTD::initialize(void)
{int i,j,k;ic=IE/2;jc=JE/2;kc=KE/2;ddx=.01;dt=ddx/(6e8);epsz=8.8e-12;/*createMang(ex);createMang(ey);createMang(ez);createMang(hx);createMang(hy);createMang(hz);createMang(dx);createMang(dy);createMang(dz);createMang(gax);createMang(***);createMang(gaz);for(k=0;k<KE;k++)	for(j=0;j<JE;j++)		for(i=0;i<IE;i++)		{			ex[i][j][k]=0.0;			ey[i][j][k]=0.0;			ez[i][j][k]=0.0;			dx[i][j][k]=0.0;			dy[i][j][k]=0.0;			dz[i][j][k]=0.0;			hx[i][j][k]=0.0;			hy[i][j][k]=0.0;			hz[i][j][k]=0.0;			gax[i][j][k]=1.0;			***[i][j][k]=1.0;			gaz[i][j][k]=1.0;		}
for(k=11;k<30;k++){	gaz[ic][jc][k]=0.;}gaz[ic][jc][kc]=0.;t0=20.0;spread=6.0;T=0;NSTEPs=1;*/
}
void ThreeDFDTD::fdtd(void)
{int n;int i,j,k;initialize();createMang(ex);createMang(ey);createMang(ez);createMang(hx);createMang(hy);createMang(hz);createMang(dx);createMang(dy);createMang(dz);createMang(gax);createMang(***);createMang(gaz);for(k=0;k<KE;k++)	for(j=0;j<JE;j++)		for(i=0;i<IE;i++)		{			ex[i][j][k]=0.0;			ey[i][j][k]=0.0;			ez[i][j][k]=0.0;			dx[i][j][k]=0.0;			dy[i][j][k]=0.0;			dz[i][j][k]=0.0;			hx[i][j][k]=0.0;			hy[i][j][k]=0.0;			hz[i][j][k]=0.0;			gax[i][j][k]=1.0;			***[i][j][k]=1.0;			gaz[i][j][k]=1.0;		}
for(k=11;k<30;k++){	gaz[ic][jc][k]=0.;}//gaz[ic][jc][kc]=0.;t0=20.0;spread=6.0;T=0;NSTEPs=1;cout << endl		 << "Please input NSTEPs= ";	cin  >> NSTEPs;		for(n=1;n<=NSTEPs;n++)	{		T=T+1;		/****Start of the main FDTD loop****/		/****Calculate the Dx field****/		for(k=1;k<KE;k++)			for(j=1;j<JE;j++)				for(i=1;i<IE;i++)				{				dx[i][j][k]=dx[i][j][k]+.5*(hz[i][j][k]-hz[i][j-1][k]-hy[i][j][k]+hy[i][j][k-1]);				}		/****Calculate the Dy field****/		for(k=1;k<KE;k++)			for(j=1;j<JE;j++)				for(i=1;i<IE;i++)				{				dy[i][j][k]=dy[i][j][k]+.5*(hx[i][j][k]-hx[i][j][k-1]-hz[i][j][k]+hz[i-1][j][k]);				}
		/****Calculate the Dz field****/		for(k=1;k<KE;k++)			for(j=1;j<JE;j++)				for(i=1;i<IE;i++)				{				dz[i][j][k]=dz[i][j][k]+.5*(hy[i][j][k]-hy[i-1][j][k]-hx[i][j][k]+hx[i][j-1][k]);				}		/****Source****/		pulse=exp(-.5*pow((t0-T)/spread,2));		dz[ic][jc][kc]=pulse;		/****Calculate the E field from D****/		for(k=1;k<KE-1;k++)			for(j=1;j<JE-1;j++)				for(i=1;i<IE-1;i++)				{
ex[i][j][k]=gax[i][j][k]*dx[i][j][k];
ey[i][j][k]=***[i][j][k]*dy[i][j][k];
ez[i][j][k]=gaz[i][j][k]*dz[i][j][k];				}		/****Calculate the Hx field****/		for(k=1;k<KE-1;k++)			for(j=1;j<JE-1;j++)				for(i=1;i<IE;i++)				{				hx[i][j][k]=hx[i][j][k]+.5*(ey[i][j][k+1]-ey[i][j][k]-ez[i][j+1][k]-ez[i][j][k]);				}		/****Calculate the Hy field****/		for(k=1;k<KE-1;k++)			for(j=1;j<JE;j++)				for(i=1;i<IE-1;i++)				{				hy[i][j][k]=hy[i][j][k]+.5*(ez[i+1][j][k]-ez[i][j][k]-ex[i][j][k+1]+ex[i][j][k]);				}		/****Calculate the Hx field****/		for(k=1;k<KE;k++)			for(j=1;j<JE-1;j++)				for(i=1;i<IE-1;i++)				{				hz[i][j][k]=hz[i][j][k]+.5*(ex[i][j+1][k]-ex[i][j][k]-ey[i+1][j][k]+ey[i][j][k]);				}		/****End of the main FDTD loop****/		/****Write the Ez field into Ez.txt****/				ofstream Ez("Ez.txt");		for(j=0;j<JE;j++)		{			for(i=0;i<IE;i++)			{				 Ez <<   ez[i][j][kc] <<" ";			}			Ez << "\n";		}		Ez.close();			for(k=1;k<KE;k++)			for(i=1;i<IE-1;i++)			{				cout << ez[i][jc][k] << " ";			}			cout << endl;	}	
    printf("\n No error there");
    freeMang(ex);freeMang(ey);freeMang(ez);freeMang(hx);freeMang(hy);freeMang(hz);freeMang(dx);freeMang(dy);freeMang(dz);freeMang(gax);freeMang(***);freeMang(gaz);
}

//Allocate dynamic memory
#define ALLOC_3D(name,nx,ny,nz,type){ \
        name=(type*)calloc(nx*ny*nz,sizeof(type));\
        if(!name){                        \
                 perror("ALLOC_3D");      \
                 exit(-1);                \
        };                                \
}                                         

#define Ez(M,N,P) ez[(M*IE+N)*JE+P]

int main()
{
   
    double *ez;
    ALLOC_3D(ez,IE,JE,KE,double);
  ......
}