adifdtd1.cpp

2D ADI FDTD code.采用不同的三对角矩阵解法

  for (it=0;it<=itmax;it++)
  {  //-------------Time interative---------->
    if(it%snapStep == 0)
      printf("it=%d\n",it);
    
    for(i=0;i<=imax;i++)
	{
      for(j=0;j<=jmax;j++)
	  {
		if(i==i0 && j==j0 && it<=nstop){/*excitation*/
		  HZINC1[i][j]=0.5*CFLN*exp(-((it)*tt-t0)*((it)*tt-t0)/(t*t));
		  HZINC2[i][j]=0.5*CFLN*exp(-((it+0.5)*tt-t0)*((it+0.5)*tt-t0)/(t*t));
		}
		else 
		{
		  HZINC1[i][j]=0.0;
		  HZINC2[i][j]=0.0;
		}
      }
	}
	
////////////////////////////////////////////////////////////////////    
    /*The start of the 1st-updating*/
    
    /*EX field component -- Y_CUT*/
    for(i=0;i<=imax;i++)
	{
      for(j=0;j<=jmax;j++)
	  {
		if(j ==0)
			aa[i][j] = -(t2e*t2h)/(yy*yy);
	    else
			aa[i][j] = -(tt*tt)/((4*epsArray[i][j] + thetaArray[i][j]*tt)*muArray[i][j-1]*yy*yy);
		
		cc[i][j] = -(tt*tt)/((4*epsArray[i][j] + thetaArray[i][j]*tt)*muArray[i][j]*yy*yy); 
		bb[i][j] = 1 -aa[i][j] - cc[i][j];

//		aa[i][j]=-(t2e*t2h)/(yy*yy);
//		bb[i][j]=1.0+2*(t2e*t2h)/(yy*yy);
//		cc[i][j]=-(t2e*t2h)/(yy*yy);
      }
	}
    
    for(i=0;i<imax;i++)
	{
      for(j=1;j<jmax;j++)
	  {
		r[i][j] = (4*epsArray[i][j] - thetaArray[i][j]*tt)/(4*epsArray[i][j] + thetaArray[i][j]*tt)*EX[i][j] +
			      + (2*tt)/((4*epsArray[i][j] + thetaArray[i][j]*tt)*yy)*(HZ[i][j] - HZ[i][j-1])
				  - (tt*tt)/(xx*yy*muArray[i][j]*(4*epsArray[i][j] + thetaArray[i][j]*tt))*(EY[i+1][j]-EY[i][j]) 
				  + (tt*tt)/(xx*yy*muArray[i][j-1]*(4*epsArray[i][j] + thetaArray[i][j]*tt))*(EY[i+1][j-1]-EY[i][j-1]);


//		r[i][j]=EX[i][j]+(t2e/yy)*HZ[i][j]-(t2e/yy)*HZ[i][j-1]
//		  -((t2e*t2h)/(xx*yy))*(EY[i+1][j]-EY[i][j]) +((t2e*t2h)/(xx*yy))*(EY[i+1][j-1]-EY[i][j-1]);
      }
	}
    
    for(i=0;i<imax;i++)
    {                   /*------------------------------------------>*/
      for(j=1;j<jmax;j++)
      {
		if(j==1)
		{
		  bet1 = bb[i][j];
		  IEX[i][j] = r[i][j]/bet1;

		}
		else// if (j>=2 && j<=jmax-1)
		{
		  g1[i][j] = cc[i][j-1]/bet1;       /*追赶法解3对角*/
		  bet1 = bb[i][j] - aa[i][j]*g1[i][j];
		  IEX[i][j] = (r[i][j] - aa[i][j]*IEX[i][j-1])/bet1;
		}
      }

      for(j=jmax-2;j>=1;j--)
      {
		IEX[i][j]=IEX[i][j]-g1[i][j+1]*IEX[i][j+1];
      }
    }               /*-------------------------------------------->*/
    
    /*EY field component -- Z_CUT*/
    
    for (i=1;i<imax;i++)
	{
      for (j=0;j<jmax;j++)
	  {
		IEY[i][j] = (4*epsArray[i][j] - thetaArray[i][j]*tt)/(4*epsArray[i][j] + thetaArray[i][j]*tt)*EY[i][j] 
			       - (2*tt)/(4*epsArray[i][j] + thetaArray[i][j]*tt)*(HZ[i][j]-HZ[i-1][j])/xx; 

//		IEY[i][j]=EY[i][j]-t2e*(HZ[i][j]-HZ[i-1][j])/xx;
	  }
	}
    
    
    /*HZ field components*/
    
    for (i=0;i<imax;i++)
	{
      for (j=0;j<jmax;j++)
	  {
		IHZ[i][j] = HZ[i][j] + HZINC1[i][j]
			+ tt/(2*muArray[i][j])*((IEX[i][j+1]-IEX[i][j])/yy - (EY[i+1][j]-EY[i][j])/xx);

//		IHZ[i][j]=HZ[i][j]+HZINC1[i][j] + t2h*((IEX[i][j+1]-IEX[i][j])/yy - (EY[i+1][j]-EY[i][j])/xx);
      }
	}
    
    /*end of the 1st-updating*/

///////////////////////////////////////////////////////////////////////    
    /*start of the 2nd-updating*/
 
    /*EX field component -- Z_CUT*/   
    for (i=0;i<imax;i++)
	{
		for (j=1;j<jmax;j++)
		{
			EX[i][j] = (4*epsArray[i][j] - thetaArray[i][j]*tt)/(4*epsArray[i][j] + thetaArray[i][j]*tt)*IEX[i][j] 
				      + (2*tt)/(4*epsArray[i][j] + thetaArray[i][j]*tt)*(IHZ[i][j]-IHZ[i][j-1])/yy;

//			EX[i][j]=IEX[i][j]+t2e*(IHZ[i][j]-IHZ[i][j-1])/yy;
		}
	}
    
    
    /*EY field component -- X_CUT*/
    
    for(i=0;i<=imax;i++)
	{
      for(j=0;j<=jmax;j++)
	  {
		  if(i==0)
			  aa[i][j]=-(t2e*t2h)/(xx*xx);
		  else
			  aa[i][j] = -(tt*tt)/((4*epsArray[i][j] + thetaArray[i][j]*tt)*muArray[i-1][j]*xx*xx);

		cc[i][j] = -(tt*tt)/((4*epsArray[i][j] + thetaArray[i][j]*tt)*muArray[i][j]*xx*xx); 
		bb[i][j] = 1 -aa[i][j] - cc[i][j];
		
//		aa[i][j]=-(t2e*t2h)/(xx*xx);
//		bb[i][j]=1.0+2*(t2e*t2h)/(xx*xx);
//		cc[i][j]=-(t2e*t2h)/(xx*xx);
      }
	}
    
    for (j=0;j<jmax;j++)
	{
      for(i=1;i<imax;i++)
	  {

		  r[i][j] = (4*epsArray[i][j] - thetaArray[i][j]*tt)/(4*epsArray[i][j] + thetaArray[i][j]*tt)*IEY[i][j] +
			      + (2*tt)/((4*epsArray[i][j] + thetaArray[i][j]*tt)*xx)*(IHZ[i-1][j] - IHZ[i][j])
				  - (tt*tt)/(xx*yy*muArray[i][j]*(4*epsArray[i][j] + thetaArray[i][j]*tt))*(IEX[i][j+1]-IEX[i][j]) 
				  + (tt*tt)/(xx*yy*muArray[i-1][j]*(4*epsArray[i][j] + thetaArray[i][j]*tt))*(IEX[i-1][j+1]-IEX[i-1][j]);


//		  r[i][j]=IEY[i][j]-(t2e/xx)*(IHZ[i][j] - IHZ[i-1][j])
//		  -((t2e*t2h)/(xx*yy))*(IEX[i][j+1]-IEX[i][j]) +((t2e*t2h)/(xx*yy))*(IEX[i-1][j+1]-IEX[i-1][j]);
      }
	}
    
    
    for (j=0;j<jmax;j++)//--------------------------------------->
	{  
		for(i=1;i<imax;i++)
		{
			if(i==1)
			{
			  bet2 = bb[i][j];
			  EY[i][j] = r[i][j]/bet2;

			}
			else// if(i>=2&&i<=imax-1)
			{
			  g2[i][j] = cc[i-1][j]/bet2;       /*追赶法解3对角*/
			  bet2 = bb[i][j] - aa[i][j]*g2[i][j];
			  EY[i][j] = (r[i][j] - aa[i][j]*EY[i-1][j])/bet2;
			}
      }

	  for(i=imax-2;i>=1;i--)
      {
		EY[i][j]=EY[i][j]-g2[i+1][j]*EY[i+1][j];
      }
    }               /*-------------------------------------------->*/

    
    /*HZ field components*/

    for (i=0;i<imax;i++)
	{
      for (j=0;j<jmax;j++)
	  {
		HZ[i][j] = IHZ[i][j] + HZINC2[i][j]
			+ tt/(2*muArray[i][j])*((IEX[i][j+1]-IEX[i][j])/yy - (EY[i+1][j]-EY[i][j])/xx);

//		HZ[i][j]=IHZ[i][j]+HZINC2[i][j]
//		  +t2h*((IEX[i][j+1]-IEX[i][j])/yy
//		  -(EY[i+1][j]-EY[i][j])/xx);
      }
	}
    
    /*end of the 2nd-updating*/
/*	
	char *strEXtemp = "EX_";
	char *strEYtemp = "EY_";
	char *strHZtemp = "HZ_";
	
	char stimulusFilename[30]; 
	char addString[40];
	if(fmod(it, snapStep)< 1e-8)
	{
		strcpy(stimulusFilename, strEXtemp);
		sprintf(addString, "%d.dat", it);
		strcat(stimulusFilename, addString);
		fp3 = fopen(stimulusFilename, "w");
		for(int ik = 0; ik < imax; ik++)
		{
			for(int jk =0; jk < jmax; jk++)
			{
				fprintf(fp3,"%15.7f   ", EX[ik][jk]);   
			}
			fprintf(fp3,"\n"); 
		}
		fclose(fp3);

		strcpy(stimulusFilename, strEYtemp);
		sprintf(addString, "%d.dat", it);
		strcat(stimulusFilename, addString);
		fp3 = fopen(stimulusFilename, "w");
		for( ik = 0; ik < imax; ik++)
		{
			for(int jk =0; jk < jmax; jk++)
			{
				fprintf(fp3,"%15.7f   ", EY[ik][jk]);   
			}
			fprintf(fp3,"\n"); 
		}
		fclose(fp3);

		strcpy(stimulusFilename, strHZtemp);
		sprintf(addString, "%d.dat", it);
		strcat(stimulusFilename, addString);
		fp3 = fopen(stimulusFilename, "w");
		for( ik = 0; ik < imax; ik++)
		{
			for(int jk =0; jk < jmax; jk++)
			{
				fprintf(fp3,"%15.7f   ", HZ[ik][jk]);   
			}
			fprintf(fp3,"\n"); 
		}
		fclose(fp3);
	}
*/

   fprintf(fp1,"%15.6f, %15.7f, %15.7f, %15.7f, %15.7f, %15.7f\n",it*tt*(10e9),EX[nobsX][nobsY],EY[nobsX][nobsY],
	   HZ[nobsX][nobsY], HZINC1[i0][j0], HZINC2[i0][j0]);   
	/*recording the field component*/   
	H1[it]=HZ[int(imax*0.9)][int(jmax*0.5)];

  }  //-------------Time interative---------->

  printf("Time used before FT  = %5d (seconds)\n",(time(NULL)-myt));
  for(k=1;k<1499;k++){
  xre1=xim1=0.0;
 ff[k]=k/(tt*itmax);
  w=2.0*pi*ff[k];

  for(i=0;i<=itmax;i++)
    {
      xre1=xre1+H1[i]*cos(w*tt*i);
      xim1=xim1-H1[i]*sin(w*tt*i);
    }
  rr1[k]=(xre1*xre1+xim1*xim1);
  fprintf(fp2,"%15.7f, %15.7f\n",ff[k]/10e+9, rr1[k]);
  }

  fclose (fp1);
  fclose (fp2);
 
  /*stuff related to C++*/
  // delete EX, EY, HZ
  for(i=0;i<=imax;i++){
    delete [] EX[i];
  }
  delete [] EX;

  for(i=0;i<=imax;i++){
    delete [] EY[i];
  }
  delete [] EY;

  for(i=0;i<=imax;i++){
    delete [] HZ[i];
  }
  delete [] HZ;

 for(i=0;i<=imax;i++){
    delete [] IEX[i];
  }
  delete [] IEX;

  for(i=0;i<=imax;i++){
    delete [] IEY[i];
  }
  delete [] IEY;

  for(i=0;i<=imax;i++){
    delete [] IHZ[i];
  }
  delete [] IHZ;

  // delete HZINC1, HZINC2
  for(i=0;i<=imax;i++)
    delete [] HZINC1[i];
  delete [] HZINC1;
 
  for(i=0;i<=imax;i++)
    delete [] HZINC2[i];
  delete [] HZINC2;
  
/////////////////////
  for(i=0;i<=imax;i++){
    delete [] thetaArray[i];
  }
  delete [] thetaArray;

  for(i=0;i<=imax;i++){
    delete [] muArray[i];
  }
  delete [] muArray;
  
  for(i=0;i<=imax;i++){
    delete [] epsArray[i];
  }
  delete [] epsArray;
/////////////////////

 for(i=0;i<=imax;i++){
    delete [] aa[i];
  }
  delete [] aa;

  for(i=0;i<=imax;i++){
    delete [] bb[i];
  }
  delete [] bb;

  for(i=0;i<=imax;i++){
    delete [] cc[i];
  }
  delete [] cc;
  
 for(i=0;i<=imax;i++){
    delete [] r[i];
  }
  delete [] r;

  for(i=0;i<=imax;i++){
    delete [] g1[i];
  }
  delete [] g1;

  for(i=0;i<=imax;i++){
    delete [] g2[i];
  }
  delete [] g2;
  
  delete [] H1;


  printf("Total CPU Time used = %5d (seconds)\n",(time(NULL)-myt));
  return 0;
}