fdtdpre.java

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

/* 3D FDTD PACKAGE in Java by Stephen Kirkup
  ==========================================
  MARK 1. Released January 2005

  This file: FDTDpre (Main program for pre-processor)

  This is a shareware java source code that carries out the FDTD (finite difference time domain)
  algorithm. The code is available for general use. The copyright remains with the author; the
  program should not be passed on to others, even in revised form, without the permission of the
  author. The author would be pleased to hear of any interesting applications of the code, ideas
  for future development, links to the website, or bugs in the code.

  The author can be contacted on [e-mail]

  The webpages maintaining this code is www.fdtd.electromagnetics.info please refer to
  it for downloads or any news or updates.

  The java compiler/interpreter can be downloaded free from wwww.java.sun.com

  The FDTD algorithm is the most popular method for transient electromagnetic simulation.

  FDTD is also very computationally-intensive hence the program has been divided into two parts: a
  pre-processor (FDTDpre.java) that carries out the computation and stores relevant results and
  a post-processor (FDTDpost.java) that carries out a simple visualisation of the results.
  The test problem is input using the .DAT file and output to .OUT files. The .MON (monitor)
  file monitors the program (to be viewed soon after the start of the program to
  check that everything is as you expected and that the program will run in a reasonable time.

  FDTDpre requires the following class files to run: Material, PhysicalConstants, Mesh, Array,
  VectorField, Efield, Hfield, Vector, Component, readfromfile.

  */




import java.io.*;

public class FDTDpre
{ public static void main (String[] args) throws IOException
  {

    int MaxNtout;
    String fileName="EXPERIMENT.DAT";
    readfromfile read = new readfromfile(fileName);

	// Total time of simulation
	double totalt;
	// The speed of electromagnetic radiation in the outer layer (Material 1)
	double c1;

	// Variables in FDTD package
	int Nx, Ny, Nz;

    Array sigma = new Array();
    Array mu = new Array();
    Array epsilon = new Array();

	// Tolerance
	final double eps=0.000001;
	final int ndlsplit=3;

	String expname;
	final int MAXMAT=20;

	// nmat, the actual number of distinct materials
	int nmat=0;

    // material object
    Material[] material = new Material[MAXMAT+1];

	// The counter through the distinct materials
	int imat;

	// Dimensions of the domain
	double xlength,ylength,zlength;

	// Cell size (Yee cell)
	double dl;

	// Set MAXCOMP, the maximum number of distinct materials
	final int MAXCOMP=50;

	// Components
	Component[] component = new Component[MAXCOMP+1];

	// ncomp, icomp the number, counter and label of the components
	int ncomp, icomp, jcomp;


	// Excitation
	Vector[] E_excite = new Vector[MAXCOMP+1];

	boolean changemat;

    // Time
    double time;
	// time step
	double dt;
	// number of time steps in simulation
	int Nt;
	// the number of time steps that are output
	int Ntout;
	// the number of time steps jumped over
	int Ntjump;
	int n;
	double time_appear, time_disappear;

	// Output Control
	//  Number of times that the electromagnetic data is saved
	int nsave;
	double sumcond,sumrel_perm,sumsusc;
	double avcond, avrel_perm, avsusc;
	double xlow,ylow,zlow,xx,yy,zz;
	double div;

	// General counters
	int i,j,k,ii,jj,kk;

    // Setting Excitation
    double sumEx, sumEy, sumEz;

	String s,monitor,stemfile,savefile;
	int    itemp;
    double temp1,temp2;
    double coeff;

  // coordinates corresponding to E and H
    double[][] Ecube = new double[4][4],
			   Hcube = new double[4][4];

    double[]    Hold = new double[4],
    	        Hnew = new double[4],
		        Eold = new double[4],
		        Enew = new double[4];

  // material properties averaged in x,y,z directions
    double[]  sigval = new double[4],
		      rhoval = new double[4],
			   muval = new double[4],
		      epsval = new double[4];

  // counters
    int
	         Jx, // Jx = 1..Nx
        	 Jy, // Jy = 1..Ny
	         Jz; // Jz = 1..Nz

  // the voltage across contacts on the capacitor and the current
    double Vab, Iab;

  // position on Graphics screen
    int IX, IZ, IXP, IZP, JXP, JZP;

    Vector point =new Vector();
    double xl,xu,yl,yu,zl,zu;


    double Exval, Eyval, Ezval, Val;

    double current;
    double	dtomu, temp, con, con1, con2;

    String name;
    double conductivity, susceptibility, rel_permittivity;
    boolean appeared;

	Mesh M;

    PhysicalConstants Const = new PhysicalConstants();


    double epsilon0=Const.epsilon0;
	double mu0=Const.mu0;
	double c0=Const.c0;

	// Voltage Probes
	final int maxpd=20;
	Vector startpoint[] = new Vector[maxpd+1];
	Vector endpoint[] = new Vector[maxpd+1];
	String pdfile[] = new String[maxpd+1];

    time=0.0;



 // Set up file names for output
    itemp=fileName.indexOf('.');
    stemfile=fileName.substring(0,itemp);
 // Monitor file: Check this file if there are errors
    monitor=stemfile+".mon";
 // Output file
    savefile=stemfile+".out";

 //setup the streams
    File outFileX = new File("x.out");
    FileOutputStream outFileStreamX = new FileOutputStream(outFileX);
    DataOutputStream outDataStreamX  = new DataOutputStream(outFileStreamX);

    File outFileY = new File("y.out");
    FileOutputStream outFileStreamY = new FileOutputStream(outFileY);
    DataOutputStream outDataStreamY  = new DataOutputStream(outFileStreamY);

    File outFileZ = new File("z.out");
    FileOutputStream outFileStreamZ = new FileOutputStream(outFileZ);
    DataOutputStream outDataStreamZ  = new DataOutputStream(outFileStreamZ);

   System.out.println("Reading Input File");
  // Get experiment name
    expname=read.get_string(5);

  // Get material properties
  //  Get number of materials and validate
    nmat=read.get_int(10,1);
    if (nmat>MAXMAT)
    { System.out.println("Number of Materials is too Large (>MAXMAT) ");
    }
    if (nmat<1)
    { System.out.println("There must be at least one material");
    }
    // Get material properties
    for (imat=1; imat<=nmat; imat++)
    {
	  name=read.get_string(12+3*imat);
      conductivity=read.get_double(13+3*imat,1);
      rel_permittivity=read.get_double(13+3*imat,2);
      susceptibility=read.get_double(13+3*imat,3);
      material[imat] = new Material(name, conductivity, rel_permittivity, susceptibility);
    }

  // Get domain and mesh
  // Get size of Yee cell
  dl=read.get_double(22+3*nmat,1);
  // Validate dl
  if (dl<=0.0)
  { System.out.println("The cell size must be positive");
    System.out.println("Input terminated");
  }
  // Get domain dimensions
  // Get xlength and validate
  xlength=read.get_double(17+3*nmat,1);
  if (xlength<dl)
  { System.out.println("x-dimension is less than the cell size");
    System.out.println("Input terminated");
  }
  ylength=read.get_double(17+3*nmat,2);
  if (ylength<dl)
  { System.out.println("y-dimension is less than the cell size");
    System.out.println("Input terminated");
  }
  zlength=read.get_double(17+3*nmat,3);
  if (zlength<dl)
  { System.out.println("z-dimension is less than the cell size");
    System.out.println("Input terminated");
  }


  ncomp=read.get_int(29+3*nmat,1);
  if (ncomp>MAXCOMP)
  { System.out.println("Number of material components is too Large (>MAXMAT) ");
  }
  if (ncomp<0)
  { System.out.println("There must be at least zero material components");
  }

 // Set time
  // Set totalt
  totalt=read.get_double(35+3*nmat+2*ncomp,1);
  if (totalt<=0.0)
  { System.out.println("Total simulation time must be positive");
  }


  for (icomp=1; icomp<=ncomp; icomp++)
  { imat=read.get_int(32+3*nmat+2*icomp,1);
    xl = read.get_double(32+3*nmat+2*icomp,2);
    xu = read.get_double(32+3*nmat+2*icomp,3);
    yl = read.get_double(32+3*nmat+2*icomp,4);
    yu = read.get_double(32+3*nmat+2*icomp,5);
    zl = read.get_double(32+3*nmat+2*icomp,6);
    zu = read.get_double(32+3*nmat+2*icomp,7);
    time_appear=read.get_double(45+3*nmat+3*ncomp+icomp,2);
    time_disappear=read.get_double(45+3*nmat+3*ncomp+icomp,3);
	component[icomp]=new Component(imat,xl,xu,yl,yu,zl,zu,time_appear,time_disappear);
  }

  // Set Excitation
  for (icomp=1; icomp<=ncomp; icomp++)
	{
		double x,y,z;
		x=read.get_double(41+3*nmat+2*ncomp+icomp,2);
		y=read.get_double(41+3*nmat+2*ncomp+icomp,3);
		z=read.get_double(41+3*nmat+2*ncomp+icomp,4);

		E_excite[icomp] = new Vector(x,y,z);
	}

    M = new Mesh(xlength,ylength,zlength,dl,totalt);

    HField H = new HField();

    EField E = new EField();

  // Set Nx, Ny and Nz and check that mesh is OK
	Nx=M.Nx;
	Ny=M.Ny;
	Nz=M.Nz;
	dt=M.dt;
	Nt=M.Nt;

  //Set MaxNtout
    MaxNtout=read.get_int(54+3*nmat+4*ncomp,1);
    System.out.println("MaxNtout = "+MaxNtout);
    Ntjump=(Nt-1)/MaxNtout+1;
    Ntout=Nt/Ntjump;

    System.out.println("Nt = "+Nt);
    System.out.println("Ntout = "+ Ntout);


// Set the speed of electromagnetic radiation in the outer field
    c1=c0/Math.sqrt((1.0+material[1].susceptibility)*material[1].rel_perm);

// Send information to the monitoring file *.mon
    System.out.println("Sending information to the monitoring file "+monitor);
    PrintWriter out =
    new PrintWriter(
		new BufferedWriter(
			new FileWriter(monitor,false)));
			   out.print("Monitoring file for ");
			   out.println(read.getFileName());
               out.println("Nt=");
               out.println(Nt);
               out.println("dt=  value           E  (exponent)");
               temp2=Math.log(dt)/Math.log(10.0);
               itemp=(int)temp2-1;
               temp1=dt/Math.pow(10.0,itemp);
               out.print(temp1);
               out.print("  ");
               out.println(itemp);
               out.println("Nx=");
               out.println(Nx);
               out.println("Ny");
               out.println(Ny);
               out.println("Nz=");
               out.println(Nz);
               out.println("Experiment name: "+ expname);
               out.print("Number of Materials ");
               out.println(nmat);
               for (imat=1; imat<=nmat; imat++)
               { out.println(material[imat].name);
                 out.print(": cond=");
                 out.print(material[imat].conductivity);
                 out.print("   rel_perm=");
                 out.print(material[imat].rel_perm);
                 out.print("   susc=");
                 out.println(material[imat].susceptibility);
               }
               out.print("Dimensions: ");
               out.print(xlength);
               out.print("   ");
               out.print(ylength);
               out.print("   ");
               out.println(zlength);
               out.print("Yee Cell Size:");
               out.println(dl);

               out.print("Number of material components:");
               out.println(ncomp);

               for (icomp=1; icomp<=ncomp; icomp++)
               { out.print("Material no ");
                 out.print(component[icomp].imat);
                 out.print(" (");
                 out.println(material[component[icomp].imat].name+") ");
                 out.print("   ");
                 out.print(component[icomp].xl);
                 out.print("   ");
                 out.print(component[icomp].xu);
                 out.print("   ");
                 out.print(component[icomp].yl);
                 out.print("   ");
                 out.print(component[icomp].yu);
                 out.print("   ");
                 out.print(component[icomp].zl);
                 out.print("   ");
                 out.print(component[icomp].zu);

                 out.println();
               }

               out.println();
               out.println("Excitation");
               for (icomp=1; icomp<=ncomp; icomp++)