fdtdpre.java

该内容是关于fdtd的源码,我用过,ok!不错!

               { out.print(icomp);
                 out.print("  ");
            	 out.print(E_excite[icomp].x+" "+E_excite[icomp].y+" "+E_excite[icomp].z);
                 out.println();
                }
               out.println();
               out.println("Times of appearance and disappearance");
               for (icomp=1; icomp<=ncomp; icomp++)
               { out.print(icomp);
                 out.print("  ");
            	 out.print(component[icomp].time_appear);
                 out.print("  ");
            	 out.println(component[icomp].time_disappear);
		       }
		       out.println("MaxNtout = "+MaxNtout);


		           int npd;

			       npd=read.get_int(58+3*nmat+4*ncomp,1);

               out.close();

			   // Open output files
			   	PrintWriter pds =  new PrintWriter(new BufferedWriter(new FileWriter("pds.out",false)));


  System.out.println("Initialise E,H data structures to zero");



 outDataStreamX.writeInt(Ntout);
 outDataStreamX.writeInt(Ny);
 outDataStreamX.writeInt(Nz);

 outDataStreamY.writeInt(Ntout);
 outDataStreamY.writeInt(Nx);
 outDataStreamY.writeInt(Nz);

 outDataStreamZ.writeInt(Ntout);
 outDataStreamZ.writeInt(Nx);
 outDataStreamZ.writeInt(Ny);




  // Loop (n) through the time steps
  for ( n=1; n<=Ntout*Ntjump+1; n++ )
  { // compute H from E

    System.out.println(n);

    time=((double)n-0.5)*dt;

    changemat=false;
    for (icomp=1; icomp<=ncomp; icomp++)
      {
      if (!component[icomp].exists(time-dt)==component[icomp].exists(time))
        {
			changemat=true;
		}
	  }
    if (changemat)
    {
	  System.out.println("Setting Up Material Property Data Structures");
	div=dl/(double)(ndlsplit+1);
    System.out.println("Material changes");
    for (i=1; i<=Nx-1; i++)
    { xlow=i*dl+div/2;
	  for (j=1; j<=Ny-1; j++)
 	  { ylow=j*dl+div/2;
	    for (k=1; k<=Nz-1; k++)
	    { zlow=k*dl+div/2;
		sumcond=0.0;
		sumrel_perm=0.0;
		for (ii=0; ii<=ndlsplit-1; ii++)
		{ point.x=xlow+ii*div;
		  for (jj=0; jj<=ndlsplit-1; jj++)
		  { point.y=ylow+jj*div;
			for (kk=0; kk<=ndlsplit-1; kk++)
			{ point.z=zlow+kk*div;
			  imat=1;
			  for (icomp=1; icomp<=ncomp; icomp++)
			    if (component[icomp].intime(point,time)) imat=component[icomp].imat;
			  sumcond=sumcond+material[imat].conductivity;
			  sumrel_perm=sumrel_perm+material[imat].rel_perm;
			}
		  }
		}
		avcond=sumcond/(double)(ndlsplit*ndlsplit*ndlsplit);
		sigma.grid[i][j][k]=avcond;
		avrel_perm=sumrel_perm/(double)(ndlsplit*ndlsplit*ndlsplit);
		epsilon.grid[i][j][k]=avrel_perm*epsilon0;
	    }
	  }
    }

    // Set mu,rho
    for (i=1; i<=Nx-1; i++)
    { xlow=i*dl+div/2-dl/2;
	  for (j=1; j<=Ny-1; j++)
	  { ylow=j*dl+div/2-dl/2;
	    for (k=1; k<=Nz-1; k++)
	    { zlow=k*dl+div/2-dl/2;
		sumsusc=0.0;
		for (ii=0; ii<=ndlsplit-1; ii++)
		{ point.x=xlow+ii*div;
		  for (jj=0; jj<=ndlsplit-1; jj++)
		  { point.y=ylow+jj*div;
			for (kk=0; kk<=ndlsplit-1; kk++)
			{ point.z=zlow+kk*div;
			  imat=1;
			  for (icomp=1; icomp<=ncomp; icomp++)
			    if (component[icomp].intime(point,time)) imat=component[icomp].imat;
			  sumsusc=sumsusc+material[imat].susceptibility;
			}
		  }
		}
		avsusc=sumsusc/(double)(ndlsplit*ndlsplit*ndlsplit);
		mu.grid[i][j][k]=mu0*(1.0+avsusc);
	    }
	  }
    }
    }
    System.out.println("Computing H");

    for ( Jx=1; Jx<=Nx-2; Jx++ )
    { for ( Jy=1; Jy<=Ny-2; Jy++ )
      { for ( Jz=1; Jz<=Nz-2; Jz++ )
        { Hold[1] = H.x.grid[Jx][Jy][Jz];
          Hold[2] = H.y.grid[Jx][Jy][Jz];
          Hold[3] = H.z.grid[Jx][Jy][Jz];
          muval[1] = (mu.grid[Jx][Jy][Jz]+mu.grid[Jx+1][Jy][Jz])/2.0;
          muval[2] = (mu.grid[Jx][Jy][Jz]+mu.grid[Jx][Jy+1][Jz])/2.0;
          muval[3] = (mu.grid[Jx][Jy][Jz]+mu.grid[Jx][Jy][Jz+1])/2.0;

          rhoval[1]=1.0e-20;
          rhoval[2]=1.0e-20;
          rhoval[3]=1.0e-20;

          Ecube[1][2] = (E.x.grid[Jx-1][Jy][Jz]-E.x.grid[Jx-1][Jy-1][Jz])/dl;
          Ecube[1][3] = (E.x.grid[Jx-1][Jy][Jz]-E.x.grid[Jx-1][Jy][Jz-1])/dl;
          Ecube[2][1] = (E.y.grid[Jx][Jy-1][Jz]-E.y.grid[Jx-1][Jy-1][Jz])/dl;
          Ecube[2][3] = (E.y.grid[Jx][Jy-1][Jz]-E.y.grid[Jx][Jy-1][Jz-1])/dl;
          Ecube[3][1] = (E.z.grid[Jx][Jy][Jz-1]-E.z.grid[Jx-1][Jy][Jz-1])/dl;
          Ecube[3][2] = (E.z.grid[Jx][Jy][Jz-1]-E.z.grid[Jx][Jy-1][Jz-1])/dl;
          Ecube[1][1] = 0.0;
          Ecube[2][2] = 0.0;
          Ecube[3][3] = 0.0;


		  // equation 3.26a
		    dtomu = dt/muval[1];
		     temp = rhoval[1]*dtomu/2.0;
			  con = 1.0+temp;
			 con1 = (1.0-temp)/con;
			 con2 = dtomu/con;
		  Hnew[1] = con1*Hold[1]+con2*(Ecube[2][3]-Ecube[3][2]);

		  // equation 3.26b
			dtomu = dt/muval[2];
			 temp = rhoval[2]*dtomu/2.0;
			  con = 1.0+temp;
			 con1 = (1.0-temp)/con;
			 con2 = dtomu/con;
	      Hnew[2] = con1*Hold[2]+con2*(Ecube[3][1]-Ecube[1][3]);

	      // equation 3.26c
			dtomu = dt/muval[3];
			 temp = rhoval[3]*dtomu/2.0;
		      con = 1.0+temp;
			 con1 = (1.0-temp)/con;
			 con2 = dtomu/con;
          Hnew[3] = con1*Hold[3]+con2*(Ecube[1][2]-Ecube[2][1]);

          // renewH( Hold, Ecube, dt, rhoval, muval, Hnew );
          H.x.grid[Jx][Jy][Jz] = Hnew[1];
          H.y.grid[Jx][Jy][Jz] = Hnew[2];
          H.z.grid[Jx][Jy][Jz] = Hnew[3];
		}
	  }
	}

// First Order Mur Correction
	System.out.println("Mur Correction in H");
	coeff=(c1*dt-dl)/(c1*dt+dl);
	H.x.Mur1(Nx,Ny,Nz,coeff);
	H.y.Mur1(Nx,Ny,Nz,coeff);
	H.z.Mur1(Nx,Ny,Nz,coeff);


	// excitation
	System.out.println("Apply Excitation");
    for (i=1; i<=Nx-1; i++)
    { for (j=1; j<=Ny-1; j++)
	  { for (k=1; k<=Nz-1; k++)
		{ jcomp=0;
		  point.x=(i+1)*dl;
		  point.y=j*dl+dl/2;
		  point.z=k*dl+dl/2;
		  for (icomp=1; icomp<=ncomp; icomp++)
		  { if (Math.abs(E_excite[icomp].x)>=eps)
			{if (component[icomp].intime(point,time)) jcomp=icomp;
			}
		  }
          if (jcomp!=0) E.x.grid[i][j][k]=E_excite[jcomp].x;
          jcomp=0;
		  point.x=i*dl+dl/2;
		  point.y=j*dl+dl;
		  point.z=k*dl+dl/2;
		  for (icomp=1; icomp<=ncomp; icomp++)
		  { if (Math.abs(E_excite[icomp].y)>=eps)
			{if (component[icomp].intime(point,time)) jcomp=icomp;
			}
		  }
		  if (jcomp!=0) E.y.grid[i][j][k]=E_excite[jcomp].y;
          jcomp=0;
		  point.x=i*dl+dl/2;
		  point.y=j*dl+dl/2;
		  point.z=k*dl+dl;
		  for (icomp=1; icomp<=ncomp; icomp++)
		  { if (Math.abs(E_excite[icomp].z)>=eps)
			{if (component[icomp].intime(point,time)) jcomp=icomp;
			}
		  }
		  if (jcomp!=0) E.z.grid[i][j][k]=E_excite[jcomp].z;
		}
	  }
    }



	pds.print(time+" ");
    for (int ipd=1; ipd<=npd; ipd++)
    {
	Vector start = new Vector();
	Vector end = new Vector();
    float pd;
    int line=60+3*nmat+4*ncomp+ipd*2;
    start.set(read.get_double(line,1),read.get_double(line,2),read.get_double(line,3));
    end.set(read.get_double(line,4),read.get_double(line,5),read.get_double(line,6));
    pd=(float)E.getPD(start,end);
    pds.print(pd+" ");
    if (ipd==1)
    {
    Vector middle = new Vector();
    start.mid(end);
    System.out.println(E.getvalue(start).y);
	}
	}
	pds.println();




	// compute E from H
    System.out.println("Computing E");
	for ( Jx=1; Jx<=(Nx-1); Jx++ )
	{ for ( Jy=1; Jy<=(Ny-1); Jy++ )
	  { for ( Jz=1; Jz<=(Nz-1); Jz++ )
	    {
		  Eold[1] = E.x.grid[Jx][Jy][Jz];
	      Eold[2] = E.y.grid[Jx][Jy][Jz];
	      Eold[3] = E.z.grid[Jx][Jy][Jz];

	      sigval[1] = (sigma.grid[Jx][Jy][Jz]+sigma.grid[Jx+1][Jy][Jz])/2.0;
	      sigval[2] = (sigma.grid[Jx][Jy][Jz]+sigma.grid[Jx][Jy+1][Jz])/2.0;
	      sigval[3] = (sigma.grid[Jx][Jy][Jz]+sigma.grid[Jx][Jy][Jz+1])/2.0;
	      epsval[1] = (epsilon.grid[Jx][Jy][Jz]+epsilon.grid[Jx+1][Jy][Jz])/2.0;
	      epsval[2] = (epsilon.grid[Jx][Jy][Jz]+epsilon.grid[Jx][Jy+1][Jz])/2.0;
	      epsval[3] = (epsilon.grid[Jx][Jy][Jz]+epsilon.grid[Jx][Jy][Jz+1])/2.0;
		  Hcube[1][2] = (H.x.grid[Jx][Jy+1][Jz+1]-H.x.grid[Jx][Jy][Jz+1])/dl;
		  Hcube[1][3] = (H.x.grid[Jx][Jy+1][Jz+1]-H.x.grid[Jx][Jy+1][Jz])/dl;
		  Hcube[2][1] = (H.y.grid[Jx+1][Jy][Jz+1]-H.y.grid[Jx][Jy][Jz+1])/dl;
		  Hcube[2][3] = (H.y.grid[Jx+1][Jy][Jz+1]-H.y.grid[Jx+1][Jy][Jz])/dl;
		  Hcube[3][1] = (H.z.grid[Jx+1][Jy+1][Jz]-H.z.grid[Jx][Jy+1][Jz])/dl;
		  Hcube[3][2] = (H.z.grid[Jx+1][Jy+1][Jz]-H.z.grid[Jx+1][Jy][Jz])/dl;
		  Hcube[1][1] = 0.0;
		  Hcube[2][2] = 0.0;
		  Hcube[3][3] = 0.0;

          double	dtoeps;

          // equation 3.27a
           dtoeps = dt/epsval[1];
             temp = sigval[1]*dtoeps/2.0;
              con = 1.0+temp;
             con1 = (1.0-temp)/con;
             con2 = dtoeps/con;
          Enew[1] = con1*Eold[1]-con2*(Hcube[2][3]-Hcube[3][2]);

          // equation 3.27b
           dtoeps = dt/epsval[2];
	         temp = sigval[2]*dtoeps/2.0;
	          con = 1.0+temp;
	         con1 = (1.0-temp)/con;
	         con2 = dtoeps/con;
          Enew[2] = con1*Eold[2]-con2*(Hcube[3][1]-Hcube[1][3]);

          // equation 3.27c
	       dtoeps = dt/epsval[3];
	         temp = sigval[3]*dtoeps/2.0;
		      con = 1.0+temp;
	         con1 = (1.0-temp)/con;
	         con2 = dtoeps/con;
          Enew[3] = con1*Eold[3]-con2*(Hcube[1][2]-Hcube[2][1]);

          // renewE( Eold, Hcube, dt, sigval, epsval, Enew );
		  E.x.grid[Jx][Jy][Jz] = Enew[1];
		  E.y.grid[Jx][Jy][Jz] = Enew[2];
		  E.z.grid[Jx][Jy][Jz] = Enew[3];

		}
	  }
	}


    for (icomp=1; icomp<=ncomp; icomp++)
   {
	   appeared=component[icomp].exists(time);
       System.out.print(appeared);
   }
    System.out.println();


	coeff=(c1*dt-dl)/(c1*dt+dl);
	E.x.Mur1(Nx,Ny,Nz,coeff);
	E.y.Mur1(Nx,Ny,Nz,coeff);
	E.z.Mur1(Nx,Ny,Nz,coeff);


// number of planar records
    int Nrec;


	// plot E-field strength on x-plane
	// E-Field Strength = sqrt ( (E.x*E.x) + (E.y*E.y) + (E.z*E.z) )
	    if (n%Ntjump==0)
		    {
				Jx=Nx/2;
				for ( Jy = 1; Jy <= Ny-2; Jy++ )
				{ for ( Jz = 1; Jz <= Nz-2; Jz++ )
			 	 {
					Exval=(E.x.grid[Jx][Jy][Jz]+E.x.grid[Jx-1][Jy][Jz])/2;
					Eyval=(E.y.grid[Jx][Jy][Jz]+E.y.grid[Jx][Jy-1][Jz])/2;
					Ezval=(E.z.grid[Jx][Jy][Jz]+E.z.grid[Jx][Jy][Jz-1])/2;
		          	Val=Math.sqrt(Exval*Exval+Eyval*Eyval+Ezval*Ezval);
				  	outDataStreamX.writeFloat((float)Val);
				  }
				}

				Jy=Ny/2;
			// plot E-field strength on y-plane
				for ( Jx = 1; Jx <= Nx-2; Jx++ )
				{ for ( Jz = 1; Jz <= Nz-2; Jz++ )
				  {
					Exval=(E.x.grid[Jx][Jy][Jz]+E.x.grid[Jx-1][Jy][Jz])/2;
					Eyval=(E.y.grid[Jx][Jy][Jz]+E.y.grid[Jx][Jy-1][Jz])/2;
					Ezval=(E.z.grid[Jx][Jy][Jz]+E.z.grid[Jx][Jy][Jz-1])/2;
		          	Val=Math.sqrt(Exval*Exval+Eyval*Eyval+Ezval*Ezval);
				  outDataStreamY.writeFloat((float)Val);
			  	}
				}

				Jz=Nz/2;
			// plot E-field strength on z-plane
				for ( Jx = 1; Jx <= Nx-2; Jx++ )
				{ for ( Jy = 1; Jy <= Ny-2; Jy++ )
				  {
					Exval=(E.x.grid[Jx][Jy][Jz]+E.x.grid[Jx-1][Jy][Jz])/2;
					Eyval=(E.y.grid[Jx][Jy][Jz]+E.y.grid[Jx][Jy-1][Jz])/2;
					Ezval=(E.z.grid[Jx][Jy][Jz]+E.z.grid[Jx][Jy][Jz-1])/2;
		          	Val=Math.sqrt(Exval*Exval+Eyval*Eyval+Ezval*Ezval);
				  outDataStreamZ.writeFloat((float)Val);
				  }
				}

	}
  }

  pds.close();

  //output done, so close the streams
  outDataStreamX.close();
  outDataStreamY.close();
  outDataStreamZ.close();

}



}