sphcalc.c

matlab编写的三维FDTD程序 pml边界条件 实现电磁波传播数值模拟和仿真

/* program sphcalc.c - calculates array of 1.0's and 0.0's for sphere  */
/* 7/18/93 */
/* 10/15/93 - modified to change array to interger array */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "findif.h"
#define SPHDIM 34
FILE *fsph,*fspar,*fsphdim;
/*main() */
int Espx[SPHDIM][SPHDIM][SPHDIM],Espy[SPHDIM][SPHDIM][SPHDIM],Espz[SPHDIM][SPHDIM][SPHDIM];
void sphcalc() 

{
short int i,j,k,radint;
int sphnum;
float y,x,a1,a2,alpha,delx,flarang,angrad;
float dist,radflt;
   if ((fspar=fopen("fdmenu.par","rb")) == NULL) {
     printf("CAN'T OPEN fdmenu.PAR FILE\n");
   }
   else {
      fread(&fd1,sizeof(struct fddat),1,fspar);
   }
   fclose(fspar);
   if ((fsphdim=fopen("sphdim.par","wb")) == NULL) {
       printf("CAN'T OPEN sphdim.par OUTPUT FILE: \n");
   }
/* MAKE 2-D ARRAY OF VALUES (0'S AND 1'S TO MULTIPLY WITH E-FIELD IN FD3DSA) */
   for (i=0;i<=SPHDIM;i++) {
      for(j=0;j<SPHDIM;j++) {
         for(k=0;k<SPHDIM;k++) {
            Espx[i][j][k]=1;
            Espy[i][j][k]=1;
            Espz[i][j][k]=1;
         }
      }
   }
   printf("dx = %f\n",fd1.dx);
   printf("radius= %f\n",fd1.rad);
   radflt=fd1.rad/fd1.dx;
   radint=radflt+3;
   printf("radius in pts= %f\n",radflt);
   for(i=0;i<2*radint+6;i++) {
      for(j=0;j<2*radint+6;j++) {
         for(k=0;k<2*radint+6;k++) {
            dist=pow((i-radint+0.5)*(i-radint+0.5)+(j-radint+0.5)*(j-radint+0.5)
                     +(k-radint+0.5)*(k-radint+0.5),0.5);
            if(radflt-dist>=0.0) {
               Espx[i+1][j][k]=0;
               Espx[i+1][j+1][k]=0;
               Espx[i+1][j+1][k+1]=0;
               Espx[i+1][j][k+1]=0;
               Espy[i][j+1][k]=0;
               Espy[i][j+1][k+1]=0;
               Espy[i+1][j+1][k]=0;
               Espy[i+1][j+1][k+1]=0;
               Espz[i][j][k+1]=0;
               Espz[i][j+1][k+1]=0;
               Espz[i+1][j+1][k+1]=0;
               Espz[i+1][j][k+1]=0;
            }
         }
      }
   }
   sphnum=SPHDIM;
   fwrite(&sphnum,sizeof(int),1,fsphdim);
   fwrite(Espx,sizeof(int),SPHDIM*SPHDIM*SPHDIM,fsphdim);
   fwrite(Espy,sizeof(int),SPHDIM*SPHDIM*SPHDIM,fsphdim);
   fwrite(Espz,sizeof(int),SPHDIM*SPHDIM*SPHDIM,fsphdim);
   fclose(fsphdim);
}