cylcalc.c

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

/* program cylcalc.c - last modified 10/15/93 */
/*
 This subroutine is used to generate an output file, called "sphdim.par"
 which contains a 3-D array of numbers that describe the geometry of a
 burried cylinder.  The file "sphdim.par' is read into FDTD3d durning program
 execution. 
*/
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "findif.h"
#define SPHDIM 10
FILE *fsph,*fspar,*fsphdim;
/*main() */
int Ecylx[SPHDIM][SPHDIM][SPHDIM],Ecyly[SPHDIM][SPHDIM][SPHDIM],Ecylz[SPHDIM][SPHDIM][SPHDIM];
void cylcalc() 

{
short int i,j,k,radint,length,xrd,yrd,zrd,len,xc,yc,zc,iflag;
float y,x,a1,a2,alpha,delx,flarang,angrad;
float dist,radflt,pi;
float xa,ya,za,theta,phi,diprad;
int cylnum;
pi=3.1415927;
   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++) {
            Ecylx[i][j][k]=1;
            Ecyly[i][j][k]=1;
            Ecylz[i][j][k]=1;
         }
      }
   }
   xc=SPHDIM/2;
   yc=xc;
   zc=xc;
   printf("xc,yc,zc = %d\n",xc);
   printf("dx = %f\n",fd1.dx);
   printf("radius= %f\n",fd1.rad);
   printf("phi (degrees from x axis)= %f\n",fd1.cylphi);
   printf("dip (degrees from horizontal)= %f\n",fd1.cyldip);
   theta=(fd1.cyldip+90.)*pi/180.;
   diprad=fd1.cyldip*pi/180.;
   phi=(fd1.cylphi)*pi/180.;
   length=fd1.objlen/fd1.dy;
   printf("length in grid cells = %d\n",length);
   radflt=fd1.rad/fd1.dx;
   printf("radius in pts= %f\n",radflt);
   radint=radflt+3;
   for(len=(-length/2);len<length/2;len++) {
      xa = xc+len*sin(theta)*cos(phi);
      ya = yc+len*sin(theta)*sin(phi);
      za = zc+len*cos(theta);
      printf(" xa,ya,za = %f,%f,%f\n",xa,ya,za);
       for(i=(-radint-3);i<radint+3;i++) {
           iflag=1;
           if(i<0) iflag=(-1);
        for(k=(-radint-3);k<radint+3;k++) {
            dist=pow((i+0.5)*(i+0.5)+(k+0.5)*(k+0.5),0.5);
            if(radflt-dist>=0.0) {
               xrd=xa+(i*sin(phi))+k*sin(diprad)*cos(phi)+0.5;
               yrd=ya+((-1.)*i*cos(phi))+k*sin(diprad)*sin(phi)+0.5;
               zrd=za+k*cos(diprad);
               if((xrd<0)||(xrd>=SPHDIM-1)) {
                 printf("xrd = %d\n",xrd);
                 printf("xrd EXCEEDING ALLOWABLE ARRAY SIZE, EXITING\n");
                 exit();
               }
               if((yrd<0)||(yrd>=SPHDIM-1)) {
                 printf("yrd = %d\n",yrd);
                 printf("yrd EXCEEDING ALLOWABLE ARRAY SIZE, EXITING\n");
                 exit();
               }
               if((zrd<0)||(zrd>=SPHDIM-1)) {
                 printf("zrd = %d\n",zrd);
                 printf("zrd EXCEEDING ALLOWABLE ARRAY SIZE, EXITING\n");
                 exit();
               }
               Ecylx[xrd+1][yrd][zrd]=0;
               Ecylx[xrd+1][yrd+1][zrd]=0;
               Ecylx[xrd+1][yrd+1][zrd+1]=0;
               Ecylx[xrd+1][yrd][zrd+1]=0;
               Ecyly[xrd][yrd+1][zrd]=0;
               Ecyly[xrd][yrd+1][zrd+1]=0;
               Ecyly[xrd+1][yrd+1][zrd]=0;
               Ecyly[xrd+1][yrd+1][zrd+1]=0;
               Ecylz[xrd][yrd][zrd+1]=0;
               Ecylz[xrd][yrd+1][zrd+1]=0;
               Ecylz[xrd+1][yrd+1][zrd+1]=0;
               Ecylz[xrd+1][yrd][zrd+1]=0;
           
            }
         }
      }
   }
   cylnum=SPHDIM;
   fwrite(&cylnum,sizeof(int),1,fsphdim);
   fwrite(Ecylx,sizeof(int),SPHDIM*SPHDIM*SPHDIM,fsphdim);
   fwrite(Ecyly,sizeof(int),SPHDIM*SPHDIM*SPHDIM,fsphdim);
   fwrite(Ecylz,sizeof(int),SPHDIM*SPHDIM*SPHDIM,fsphdim);
   fclose(fsphdim);
}