计算二维TE波的ADI-FDTD方法的程序.rar

/* this is the program for 2D ADI-FDTD method (TE wave) - a point source in a cavity (dimensions of the cavity are defined with a and b) that is closed with PECs at all the four outer boundaries. The 1st and 2nd resonance frequencies of the cavity are 150GHz and 300 GHz*/ /*the final output of this program is the resonance frequency of the 2-D cavity*/ #include <stdio.h> #include <stdlib.h> #include <math.h> #include <time.h> #include <iostream> typedef double *PTR; typedef double **DPTR; typedef double **TPTR; int main () { unsigned int myt = time(NULL); // #define CFLN 5 /*you can change this number to 1,2,3,4,5,6,...*/ // #define imax 100 /*number of cells in x direction*/ // #define jmax 50 /*number of cells in y direction*/ // #define itmax 200000 /*total iterations*/ // #define snapStep 1000 double eps0=8.854e-12; double xmu0,xx,yy; // double a=1.0e-3,b=0.5e-3; /*a is the the length of the cavity and b is the width*/ double e,m,ttt,tt,t,t0,c0,pi,t2e,t2h; int i,j,i0,j0,it,k, nobsX, nobsY; int nstop; float a, b, fmax, sourceX, sourceY, obsX, obsY; float leftX, leftY, rightX, rightY, thetaA; int nleftX, nleftY, nrightX, nrightY; int nBool; double xre1,xim1; /*used for F-transform*/ double w,rr1[1500],ff[1500];/*used for F-transform*/ float CFLN; int imax, jmax, itmax, snapStep; printf("==================Object Grid Configuration==============\n"); printf("The X size of object(mm): "); scanf("%f", &a); printf("The Y size of object(mm): "); scanf("%f", &b); printf("The number of cells in X: "); scanf("%d", &imax); printf("The number of cells in Y:"); scanf("%d", &jmax); xx = a*(1e-3)/imax; yy = b*(1e-3)/jmax; ///////////////////////////////// printf("If the domain contain the inner conductor (1 or 0): "); scanf("%d", &nBool); if(nBool == 1) { /// conductor object printf("The conductivity of the inner conductor: "); scanf("%f", &thetaA); printf("The inner conductor coordiation, bottom left X (mm): "); scanf("%f", &leftX); printf("The inner conductor coordiation, bottom left Y (mm): : "); scanf("%f", &leftY); printf("The inner conductor coordiation, top right X (mm): "); scanf("%f", &rightX); printf("The inner conductor coordiation, top right Y (mm): : "); scanf("%f", &rightY); nleftX = int(leftX*(1e-3)/xx); nleftY = int(leftY*(1e-3)/yy); nrightX = int(rightX*(1e-3)/xx); nrightY = int(rightY*(1e-3)/yy); } printf("==============Excitation and Time Configuration==========\n"); printf("The source location (mm): X = "); scanf("%f", &sourceX); printf("The source location (mm): Y = "); scanf("%f", &sourceY); i0 = int(sourceX*(1e-3)/xx); j0 = int(sourceY*(1e-3)/yy); printf("The CFLN: "); scanf("%f", &CFLN); printf("The total iterations in time: "); scanf("%d", &itmax); printf("The max frequency(GHz): "); scanf("%f", &fmax); fmax = fmax*1e9; printf("==============Observation Point==========\n"); printf("The number of snapstep: "); scanf("%d", &snapStep); printf("The observation point (mm): X = "); scanf("%f", &obsX); printf("The observation point (mm): Y = "); scanf("%f", &obsY); nobsX = int(obsX*(1e-3)/xx); nobsY = int(obsY*(1e-3)/yy); FILE *fp1,*fp2, *fp3; char stimulusFilename1[40], stimulusFilename2[40], stimulusFilename3[40]; char addString1[40], addString2[40], addString3[40]; if(nBool == 1) sprintf(addString1, "DoolittleField_Cond_CFLN=%f.csv", CFLN); else sprintf(addString1, "DoolittleField_CFLN=%f.csv", CFLN); strcpy(stimulusFilename1, addString1); fp1=fopen(stimulusFilename1,"w"); if(nBool == 1) sprintf(addString2, "DoolittleFreq_Cond_CFLN=%f.csv", CFLN); else sprintf(addString2, "DoolittleFreq_CFLN=%f.csv", CFLN); strcpy(stimulusFilename2, addString2); fp2=fopen(stimulusFilename2,"w"); if(nBool == 1) sprintf(addString3, "DoolittleConfig_Cond_CFLN=%f.csv", CFLN); else sprintf(addString3, "DoolittleConfig_CFLN=%f.csv", CFLN); strcpy(stimulusFilename3, addString3); fp3=fopen(stimulusFilename3,"w"); fprintf(fp3, "The size of object: X = %f mm, Y = %f mm\n",a, b); fprintf(fp3, "imax = %d, jmax = %d\n", imax,jmax); fprintf(fp3, "xx = %f(mm),yy = %(mm), ITMAX = %d\n",xx/1.0e-3,yy/1.0e-3,itmax); fprintf(fp3, "The max frequency(GHz): %10.7f GHz\n", fmax/1e9); fprintf(fp3, "The CFLN = %f\n", CFLN); fprintf(fp3, "The source location (mm): X = %f mm , Y = %f mm\n", sourceX, sourceY); fprintf(fp3, "The observation point (mm): X = %f mm , Y = %f mm\n", obsX, obsY); fclose(fp3); printf("imax=%5d,jmax=%5d\n",imax,jmax); printf("xx=%12.6f(mm),yy=%12.6f(mm), ITMAX=%d\n",xx/1.0e-3,yy/1.0e-3,itmax); c0=3.0e8; xmu0=1.0/(eps0*c0*c0); e=eps0; m=xmu0; /*position of the excitation, you can use almost arbitrary number for i0 and j0*/ // i0 = int(imax/10); // j0 = int(jmax/2); ttt=1.0/(c0*sqrt(1.0/(xx*xx)+1.0/(yy*yy))); tt=CFLN*ttt;/*CFLN number applied here*/ printf("ttt=%e,tt=%e,tt/ttt=%10.5f\n",ttt,tt,tt/ttt); //fmax=4e9;/*maximum frequcncy used for the gaussian pulse*/ t=1.0/(2.0*fmax);/*used for gaussian pulse*/ t0=4.0*t;/*used for gaussian pulse*/ nstop=(int)(2*t0/tt); pi=4.0*atan(1.0); printf("c0=%e,nstop=%d\n",c0,nstop); /*some coefficients*/ t2e=tt/(2.0*e); t2h=tt/(2.0*m); /*stuff related to C++*/ // all required variables DPTR EX,EY,HZ,HZINC1,HZINC2,IEX,IEY,IHZ; /*note- EX, EY and HZ are used to store the field components at time step n and n+1; IEX,IEY,and IHZ are used to store the field components at time step n+1/2; HZINC1 and HZINC2 are used for the excitation*/ double bet1, bet2; DPTR muArray, thetaArray, epsArray; /*used for materials array*/ DPTR aa, bb, cc, r, g1, g2; /*used for solving the tri-diagonal matrix*/ PTR H1;/*used to store the recording field components at a certain point to calculate the resonance frequency of the cavity*/ // allocate memories for each variable H1=new double [itmax+1]; for (i=0;i<itmax+1;i++){ H1[i]=0.0; } //// EX, EY, HZ //////// EX = new PTR[(imax+1)]; for(i=0;i<=imax; i++ ) { EX[i] = new double[jmax+1]; for(j=0;j<=jmax;j++) { EX[i][j] = 0.0; } } EY = new PTR[(imax+1)]; for(i=0;i<=imax;i++) { EY[i] = new double[jmax+1]; for(j=0;j<=jmax;j++) { EY[i][j] = 0.0; } } HZ = new PTR[(imax+1)]; for(i=0;i<=imax;i++) { HZ[i] = new double[jmax+1]; for(j=0;j<=jmax;j++) { HZ[i][j] = 0.0; } } IEX = new PTR[(imax+1)]; for(i=0;i<=imax; i++ ) { IEX[i] = new double[jmax+1]; for(j=0;j<=jmax;j++) { IEX[i][j] = 0.0; } } IEY = new PTR[(imax+1)]; for(i=0;i<=imax;i++) { IEY[i] = new double[jmax+1]; for(j=0;j<=jmax;j++) { IEY[i][j] = 0.0; } } IHZ = new PTR[(imax+1)]; for(i=0;i<=imax;i++) { IHZ[i] = new double[jmax+1]; for(j=0;j<=jmax;j++) { IHZ[i][j] = 0.0; } } HZINC1 = new PTR[(imax+1)]; for(i=0;i<=imax;i++) { HZINC1[i] = new double[jmax+1]; for(j=0;j<=jmax;j++) { HZINC1[i][j] = 0.0; } } HZINC2 = new PTR[(imax+1)]; for(i=0;i<=imax;i++) { HZINC2[i] = new double[jmax+1]; for(j=0;j<=jmax;j++) { HZINC2[i][j] = 0.0; } } //////materials ////////// muArray = new PTR[(imax+1)]; for(i=0;i<=imax; i++ ) { muArray[i] = new double[jmax+1]; for(j=0;j<=jmax;j++) { muArray[i][j] = xmu0; } } thetaArray = new PTR[(imax+1)]; for(i=0;i<=imax; i++ ) { thetaArray[i] = new double[jmax+1]; for(j=0;j<=jmax;j++) { thetaArray[i][j] = 0.0; } } epsArray = new PTR[(imax+1)]