有问题的程序.doc

#pragma hdrstop
#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
  
//源程序
double g(double t)
{   
  const  double pi=3.1415926;
  const  double w=2*pi*1e9;
  if((t>0)&&(t<1e-9))
   return ((1-cos(w*t))/2);
  else
   return 0;
}
int main(int argc, char* argv[])
{   
  double dz[101][101];
  double ez[101][101],hx[101][101],hy[101][101],hxrec[101][4][2],hyrec[4][101][2];
  double murl_coef1,murl_coef2;
  double pi,c,u0,epsilon0;
     double t,dx,dt,dy;

  int maxgridsx,maxgridsy;
  int m,n,i,m1,m2,m3,m4,n1,n2,n3,n4;
  int msrc,nsrc;
  int nxwatch,nywatch,ntwin,ndtwin;
  double twin1,twin2,tunit,tcalc;
  double coef1,coef2,coef3,coef4;
  int ntotal,ncnt;

  FILE  *fp;
  //物理常数
  pi=cos(-1.0);
  c=3e8;
  u0=4*pi*1e-7;
  epsilon0=1e-9/(36*pi);
  //剖分间隔
  dx=0.01;
  dt=dx/c/sqrt(2.0);
  dy=dx;
  //输入网格数及线源位置
  puts("Please input max grids of x: ");
  scanf("%i",&maxgridsx);
  puts("Please input max grids of y:");
  scanf("%i",&maxgridsy);
  puts("Please input x position of line source");
  scanf("%i",&msrc);
  puts("Please input y position of line source");
  scanf("%i",&nsrc);
  
   
  
  /*金属导体边界(m3,n3)-(m4,n4)*/
  m3=80;
  m4=m3+10;
  n3=10;
  n4=n3+10;
  /*观察点参数*/
  nxwatch=msrc;
  nywatch=nsrc+5;
  twin1=0;
  twin2=20e-9;
  ntwin=200;
  tunit=1e-9;
  tcalc=twin2;
  /*初始化数组*/
  for(m=0;m<maxgridsx;m++)
   for(n=0;n<maxgridsy;n++)
   {   dz[m][n]=0;
       ez[m][n]=0;
   }
   
   for(m=0;m<maxgridsx;m++)
    for(n=0;n<maxgridsy+1;n++)
     hx[m][n]=0;
  for(m=0;m<maxgridsx+1;m++)
  for(n=0;n<maxgridsy;n++)
   hy[m][n]=0;
      for(m=0;m<maxgridsx;m++)
   for(n=0;n<4;n++)
    for(i=0;i<2;i++)
     hxrec[m][n]=0;
      
   for(m=0;m<4;m++)
    for(n=0;n<maxgridsy;n++)
     for(i=0;i<2;i++)
      hyrec[m][n]=0;
/*有关系数*/
    coef1=dt/dx;
    coef2=dt/dy;
    coef3=dt/dx/u0;
    coef4=dt/dy/u0;
    murl_coef1=(c*dt-dx)/(c*dt+dx);
    murl_coef2=2*dx/(c*dt+dx);

    ndtwin=(int)((twin2-twin1)/ntwin/dt);
    if(0==ndtwin)
     ndtwin=1;
    /*打开数据文件*/
       fp=fopen("FDTD.dat","w");
    fputs("time data\n",fp);
    ntotal=(int)(tcalc/dt);
    ncnt=-1;
   /*开始FDTD迭代*/
    for(t=0;t<=tcalc;t+=dt)
    {   
          ncnt=ncnt+1;
    /*dz的更新方程*/
    for(m=0;m<maxgridsx;m++)
     for(n=0;n<maxgridsy;n++)
      dz[m][n]+=coef1*(hy[m+1][n]-hy[m][n])-coef2*(hx[m][n+1]-hx[m][n]);
      /*更新ez*/
    for(m=0;m<maxgridsx;m++)
     for(n=0;n<maxgridsy;n++)
      ez[m][n]=dz[m][n]/epsilon0;
     /*加入线源*/
      
     
   /*设置金属导体*/
    for(m=m3;m<m4;m++)
     for(n=n3;n<n4;n++)
      ez[m][n]=0;
   /*更新hx*/
   for(m=0;m<maxgridsx;m++)
    for(n=0;n<maxgridsy;n++)
     hx[m][n]-=coef4*(ez[m][n]-ez[m][n-1]);
   /*hx1阶Mur吸收边界条件*/
             for(m=0;m<maxgridsx;m++)
    {
      hx[m][0]=-hxrec[m][1][1]+murl_coef1*(hx[m][1]+hxrec[m][0][1])+murl_coef2*(hxrec[m][0][0]+hxrec[m][1][0]);
      hx[m][maxgridsy]=-hxrec[m][2][1]+murl_coef1*(hx[m][maxgridsy-1]+hxrec[m][3][1])+murl_coef2*(hxrec[m][3][0]+hxrec[m][2][0]);
      for(i=0;i<2;i++)
      {
           hxrec[m][1]=hxrec[m][0];
              hxrec[m][0]=hx[m];
           hxrec[m][2+i][1]=hxrec[m][2+i][0];
           hxrec[m][2+i][0]=hx[m][maxgridsy-1+i];
      }
    }
    /*更新hy*/
    for(m=1;m<maxgridsx;m++)
     for(n=0;n<maxgridsy;n++)
      hy[m][n]+=coef3*(ez[m][n]-ez[m-1][n]);
    /*hy吸收边界条件*/
     for(n=0;n<maxgridsy;n++)
     {
      hy[0][n]=-hyrec[1][n][1]+murl_coef1*(hy[1][n]+hyrec[0][n][1])+murl_coef2*(hyrec[0][n][0]+hyrec[1][n][0]);
         hy[maxgridsx][n]=-hyrec[2][n][1]+murl_coef1*(hy[maxgridsx-1][n]+hyrec[3][n][1])+murl_coef2*(hyrec[3][n][0]+hyrec[2][n][0]);
           for(i=0;i<2;i++)
     {
      hyrec[n][1]=hyrec[n][0];
         hyrec[n][0]=hy[n];
      hyrec[2+i][n][1]=hyrec[2+i][n][0];
      hyrec[2+i][n][0]=hy[maxgridsx-1+i][n];
     }   
     
     
     }
     /*输出电场分量twin1~twin1+twatch*/
     if((t>=twin1)&&(t<=twin2)&&(0==(ncnt%ndtwin)))
     {   
        fprintf(fp,"%e%e\n",t/tunit,ez[nxwatch][nywatch]);
        printf("%i/%i\n",ncnt,ntotal);
                  ez[msrc][nsrc]=g(t);
      
         
     }
  }

        return 0;
}