pipcalc3.c

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

/* program pipcalc.c - last modified 4/2/94 */
/* 
   This subroutine is used to generate an output file "pipdim.par" which  
   contains a 3-D array of numbers that describe the geometry of a buried
   pipe. the file "pipdim.par" is read into fdtd3d during program execution.
*/
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "findif.h"
#define PIPDIMH 16
#define GRIDXY 11
FILE *fsph,*fspar,*fsphdim;
/*main() */
int Epipx[GRIDXY][GRIDXY][PIPDIMH],Epipy[GRIDXY][GRIDXY][PIPDIMH];
int Epipz[GRIDXY][GRIDXY][PIPDIMH];
void pipcalc3() 

{
short int i,j,k,radint,length,xrd,yrd,zrd,len,xc,yc,zc,iflag,xaint,yaint,zaint;
float y,x,a1,a2,alpha,delx,flarang,angrad,costheta;
float dist,radflt,pi;
float xa,ya,za,theta,phi,diprad;
int horizdim,vertdim;
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("pipdim.par","wb")) == NULL) {
       printf("CAN'T OPEN pipdim.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<GRIDXY;i++) {
      for(j=0;j<GRIDXY;j++) {
         for(k=0;k<PIPDIMH;k++) {
            Epipx[i][j][k]=1;
            Epipy[i][j][k]=1;
            Epipz[i][j][k]=1;
         }
      }
   }
   xc=GRIDXY/2;
   yc=GRIDXY/2;
   zc=PIPDIMH/2;
   printf("xc,yc,zc = %d,%d,%d\n",xc,yc,zc);
   printf("dx = %f\n",fd1.dx);
   printf("radius= %f\n",fd1.rad);
   printf("phi (degrees from x axis)= %f\n",fd1.cylphi);
   phi=(fd1.cylphi)*pi/180.;
   length=fd1.objlen/fd1.dy;
   printf("length in grid cells = %d\n",length);
   radflt=fd1.rad/fd1.dx;
   if((radflt<1.0)&&(radflt>=0.5)) printf("DIAMETER LESS THAN 2 POINTS - TRUNCATING TO 1 CUBE\n");
   
   if(radflt<0.5) {
     printf("DIAMETER LESS THAN 1 POINT - TRUNCATING TO CUBE EDGE\n");
     printf("phi (relative to x-axis) = %d\n",(int) fd1.cylphi);
     if(((int)fd1.cylphi!=0)&&((int)fd1.cylphi!=90)) {
       printf("PHI MUST BE SPECIFIED AS 0 OR 90 DEGREES - EXITING\n");
       exit();
     }
   }
   printf("radius in pts= %f\n",radflt);
   radint=radflt+3;
   theta=pi/2.;
   costheta=0.0;
   diprad=0.0;
   if(radflt>=1.0) {
      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*costheta;
         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>=GRIDXY-1)) {
                     printf("xrd = %d\n",xrd);
                     printf("xrd EXCEEDING ALLOWABLE ARRAY SIZE, EXITING\n");
                     exit();
                  }
                  if((yrd<0)||(yrd>=GRIDXY-1)) {
                     printf("yrd = %d\n",yrd);
                     printf("yrd EXCEEDING ALLOWABLE ARRAY SIZE, EXITING\n");
                     exit();
                  }
                  if((zrd<0)||(zrd>=PIPDIMH-1)) {
                    printf("zrd = %d\n",zrd);
                    printf("zrd EXCEEDING ALLOWABLE ARRAY SIZE, EXITING\n");
                    exit();
                  }
                  printf(" xa,ya,za,xrd,yrd,zrd = %f,%f,%f %d %d %d\n",xa,ya,za,xrd,yrd,zrd);
                  Epipx[xrd+1][yrd][zrd]=0;
                  Epipx[xrd+1][yrd+1][zrd]=0;
                  Epipx[xrd+1][yrd+1][zrd+1]=0;
                  Epipx[xrd+1][yrd][zrd+1]=0;
                  Epipy[xrd][yrd+1][zrd]=0;
                  Epipy[xrd][yrd+1][zrd+1]=0;
                  Epipy[xrd+1][yrd+1][zrd]=0;
                  Epipy[xrd+1][yrd+1][zrd+1]=0;
                  Epipz[xrd][yrd][zrd+1]=0;
                  Epipz[xrd][yrd+1][zrd+1]=0;
                  Epipz[xrd+1][yrd+1][zrd+1]=0;
                  Epipz[xrd+1][yrd][zrd+1]=0;
           
               }
            }
         }
      }
   }
   else if((radflt>=0.5)&&(radflt<1.0)) {
      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);
         xaint=xa+0.5;
         yaint=ya+0.5;
         zaint=za+0.5;
                  if((xaint<0)||(xaint>=GRIDXY-1)) {
                     printf("xaint = %d\n",xaint);
                     printf("xaint EXCEEDING ALLOWABLE ARRAY SIZE, EXITING\n");
                     exit();
                  }
                  if((yaint<0)||(yaint>=GRIDXY-1)) {
                     printf("yaint = %d\n",yaint);
                     printf("yaint EXCEEDING ALLOWABLE ARRAY SIZE, EXITING\n");
                     exit();
                  }
                  if((zaint<0)||(zaint>=PIPDIMH-1)) {
                    printf("zaint = %d\n",zaint);
                    printf("zaint EXCEEDING ALLOWABLE ARRAY SIZE, EXITING\n");
                    exit();
                  }
         printf(" xaint,yaint,zaint = %d,%d,%d\n",xaint,yaint,zaint);
         Epipx[xaint+1][yaint][zaint]=0;
         Epipx[xaint+1][yaint+1][zaint]=0;
         Epipx[xaint+1][yaint+1][zaint+1]=0;
         Epipx[xaint+1][yaint][zaint+1]=0;
         Epipy[xaint][yaint+1][zaint]=0;
         Epipy[xaint][yaint+1][zaint+1]=0;
         Epipy[xaint+1][yaint+1][zaint]=0;
         Epipy[xaint+1][yaint+1][zaint+1]=0;
         Epipz[xaint][yaint][zaint+1]=0;
         Epipz[xaint][yaint+1][zaint+1]=0;
         Epipz[xaint+1][yaint+1][zaint+1]=0;
         Epipz[xaint+1][yaint][zaint+1]=0;
       }
    }
   else if(radflt<0.5) {
      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);
         xaint=xa+0.5;
         yaint=ya+0.5;
         zaint=za+0.5;
                  if((xaint<0)||(xaint>=GRIDXY-1)) {
                     printf("xaint = %d\n",xaint);
                     printf("xaint EXCEEDING ALLOWABLE ARRAY SIZE, EXITING\n");
                     exit();
                  }
                  if((yaint<0)||(yaint>=GRIDXY-1)) {
                     printf("yaint = %d\n",yaint);
                     printf("yaint EXCEEDING ALLOWABLE ARRAY SIZE, EXITING\n");
                     exit();
                  }
                  if((zaint<0)||(zaint>=PIPDIMH-1)) {
                    printf("zaint = %d\n",zaint);
                    printf("zaint EXCEEDING ALLOWABLE ARRAY SIZE, EXITING\n");
                    exit();
                  }
         if((int)fd1.cylphi==0) {
              Epipx[xaint+1][yaint][zaint]=0;
              printf(" xaint,yaint,zaint = %d,%d,%d\n",xaint,yaint,zaint);
         }
         else if((int)fd1.cylphi==90){
              Epipy[xaint][yaint+1][zaint]=0;
              printf(" xaint,yaint,zaint = %d,%d,%d\n",xaint,yaint,zaint);
         }
       }
    }
   
   horizdim=GRIDXY;
   vertdim =PIPDIMH;
   fwrite(&horizdim,sizeof(int),1,fsphdim);
   fwrite(&vertdim,sizeof(int),1,fsphdim);
   fwrite(Epipx,sizeof(int),GRIDXY*GRIDXY*PIPDIMH,fsphdim);
   fwrite(Epipy,sizeof(int),GRIDXY*GRIDXY*PIPDIMH,fsphdim);
   fwrite(Epipz,sizeof(int),GRIDXY*GRIDXY*PIPDIMH,fsphdim);
   fclose(fsphdim);
}