meep.hpp

麻省理工的计算光子晶体的程序

				const symmetry &S, int sn,
				void *chunkloop_data);
typedef complex<double> (*field_function)(const complex<double> *fields,
					   const vec &loc,
					   void *integrand_data_);
typedef double (*field_rfunction)(const complex<double> *fields,
				   const vec &loc,
				   void *integrand_data_);

field_rfunction derived_component_func(derived_component c, const volume &v,
				       int &nfields, component cs[6]);

class fields {
 public:
  int num_chunks;
  fields_chunk **chunks;
  src_time *sources;
  flux_vol *fluxes;
  symmetry S;

  // The following is an array that is num_chunks by num_chunks.  Actually
  // it is two arrays, one for the imaginary and one for the real part.
  double **comm_blocks[NUM_FIELD_TYPES];
  // This is the same size as each comm_blocks array, and store the sizes
  // of the comm blocks themselves for each connection-phase type
  int *comm_sizes[NUM_FIELD_TYPES][CONNECT_COPY+1];
  int comm_size_tot(int f, int pair) const {
    int sum = 0; for (int ip=0; ip<3; ++ip) sum+=comm_sizes[f][ip][pair];
    return sum;
  }

  double a, dt; // The resolution a and timestep dt=Courant/a
  volume v, user_volume;
  geometric_volume gv;
  double m;
  int t, phasein_time, is_real;
  complex<double> k[5], eikna[5];
  double coskna[5], sinkna[5];
  boundary_condition boundaries[2][5];
  bandsdata *bands;
  char *outdir;
  // fields.cpp methods:
  fields(structure *, double m=0);
  fields(const fields &);
  ~fields();
  bool equal_layout(const fields &f) const;
  void use_real_fields();
  void zero_fields();
  void remove_sources();
  void remove_polarizabilities();
  void remove_fluxes();
  void reset();
  bool disable_sources; // set to true to turn off sources (w/o deleting)
  // time.cpp
  double time_spent_on(time_sink);
  void print_times();
  // boundaries.cpp
  void set_boundary(boundary_side,direction,boundary_condition);
  void use_bloch(direction d, double k) { use_bloch(d, (complex<double>) k); }
  void use_bloch(direction, complex<double> kz);
  void use_bloch(const vec &k);
  vec lattice_vector(direction) const;

  geometric_volume total_volume(void) const;

  // h5fields.cpp:
  // low-level function:
  void output_hdf5(h5file *file, const char *dataname,
		   int num_fields, const component *components,
		   field_function fun, void *fun_data_, int reim,
		   const geometric_volume &where,
		   bool append_data = false,
		   bool single_precision = false);
  // higher-level functions
  void output_hdf5(const char *dataname,  // OUTPUT COMPLEX-VALUED FUNCTION
		   int num_fields, const component *components,
		   field_function fun, void *fun_data_,
		   const geometric_volume &where,
		   h5file *file = 0,
		   bool append_data = false,
		   bool single_precision = false,
		   const char *prefix = 0,
		   bool real_part_only = false);
  void output_hdf5(const char *dataname,  // OUTPUT REAL-VALUED FUNCTION
		   int num_fields, const component *components,
		   field_rfunction fun, void *fun_data_,
		   const geometric_volume &where,
		   h5file *file = 0,
		   bool append_data = false,
		   bool single_precision = false,
		   const char *prefix = 0);
  void output_hdf5(component c,   // OUTPUT FIELD COMPONENT (or Dielectric)
		   const geometric_volume &where,
		   h5file *file = 0,
		   bool append_data = false,
		   bool single_precision = false,
		   const char *prefix = 0);
  void output_hdf5(derived_component c,   // OUTPUT DERIVED FIELD COMPONENT
		   const geometric_volume &where,
		   h5file *file = 0,
		   bool append_data = false,
		   bool single_precision = false,
		   const char *prefix = 0);
  h5file *open_h5file(const char *name, 
		      h5file::access_mode mode = h5file::WRITE,
		      const char *prefix = NULL, bool timestamp = false);
  const char *h5file_name(const char *name,
			  const char *prefix = NULL, bool timestamp = false);

  // step.cpp methods:
  double last_step_output_wall_time;
  int last_step_output_t;
  void step();
  inline double time() const { return t*dt; };

  // cw_fields.cpp:
  bool solve_cw(double tol, int maxiters, complex<double> frequency, int L=2);
  bool solve_cw(double tol = 1e-8, int maxiters = 10000, int L=2);

  // sources.cpp:
  double last_source_time();
  void add_point_source(component c, double freq, double width, double peaktime,
                        double cutoff, const vec &, complex<double> amp = 1.0,
                        int is_continuous = 0);
  void add_point_source(component c, const src_time &src,
                        const vec &, complex<double> amp = 1.0);
  void add_volume_source(component c, const src_time &src,
			 const geometric_volume &, 
			 complex<double> A(const vec &),
			 complex<double> amp = 1.0);
  void add_volume_source(component c, const src_time &src,
			 const geometric_volume &, 
			 complex<double> amp = 1.0);
  void require_component(component c);

  // initialize.cpp:
  void initialize_field(component, complex<double> f(const vec &));
  void initialize_A(complex<double> A(component, const vec &), double freq);
  void initialize_with_nth_te(int n);
  void initialize_with_nth_tm(int n);
  void initialize_with_n_te(int n);
  void initialize_with_n_tm(int n);
  void initialize_polarizations();
  int phase_in_material(const structure *s, double time);
  int is_phasing();

  // loop_in_chunks.cpp
  void loop_in_chunks(field_chunkloop chunkloop, void *chunkloop_data,
		      const geometric_volume &where,
		      component cgrid = Dielectric,
		      bool use_symmetry = true,
		      bool snap_unit_dims = false);
  
  // integrate.cpp
  complex<double> integrate(int num_fields, const component *components,
			    field_function fun, void *fun_data_,
			    const geometric_volume &where,
			    double *maxabs = 0);
  double integrate(int num_fields, const component *components,
		   field_rfunction fun, void *fun_data_,
		   const geometric_volume &where,
		   double *maxabs = 0);
  double max_abs(int num_fields, const component *components,
		 field_function fun, void *fun_data_,
		 const geometric_volume &where);
  double max_abs(int num_fields, const component *components,
		 field_rfunction fun, void *fun_data_,
		 const geometric_volume &where);

  double max_abs(int c, const geometric_volume &where);
  double max_abs(component c, const geometric_volume &where);
  double max_abs(derived_component c, const geometric_volume &where);
  
  // dft.cpp
  dft_chunk *add_dft(component c, const geometric_volume &where,
		     double freq_min, double freq_max, int Nfreq,
		     bool include_dV_and_interp_weights = true,
		     complex<double> weight = 1.0, dft_chunk *chunk_next = 0);
  dft_chunk *add_dft_pt(component c, const vec &where,
			double freq_min, double freq_max, int Nfreq);
  dft_chunk *add_dft(const geometric_volume_list *where,
		     double freq_min, double freq_max, int Nfreq,
		     bool include_dV = true);
  void update_dfts();
  dft_flux add_dft_flux(direction d, const geometric_volume &where,
			double freq_min, double freq_max, int Nfreq);
  dft_flux add_dft_flux_box(const geometric_volume &where,
			    double freq_min, double freq_max, int Nfreq);
  dft_flux add_dft_flux_plane(const geometric_volume &where,
			      double freq_min, double freq_max, int Nfreq);
  dft_flux add_dft_flux(const geometric_volume_list *where,
			double freq_min, double freq_max, int Nfreq);
  
  // monitor.cpp
  double get_inveps(component, direction, const vec &loc) const;
  double get_eps(const vec &loc) const;
  void get_point(monitor_point *p, const vec &) const;
  monitor_point *get_new_point(const vec &, monitor_point *p=NULL) const;
  complex<double> analytic_epsilon(double freq, const vec &) const;
  
  void prepare_for_bands(const vec &, double end_time, double fmax=0,
                         double qmin=1e300, double frac_pow_min=0.0);
  void record_bands();
  complex<double> get_band(int n, int maxbands=100);
  void grace_bands(grace *, int maxbands=100);
  void output_bands(FILE *, const char *, int maxbands=100);
  complex<double> get_field(int c, const vec &loc) const;
  complex<double> get_field(component c, const vec &loc) const;
  double get_field(derived_component c, const vec &loc) const;

  // energy_and_flux.cpp
  double energy_in_box(const geometric_volume &);
  double electric_energy_in_box(const geometric_volume &);
  double magnetic_energy_in_box(const geometric_volume &);
  double thermo_energy_in_box(const geometric_volume &);
  double total_energy();
  double field_energy_in_box(const geometric_volume &);
  double field_energy_in_box(component c, const geometric_volume &);
  double field_energy();
  double flux_in_box_wrongH(direction d, const geometric_volume &);
  double flux_in_box(direction d, const geometric_volume &);
  flux_vol *add_flux_vol(direction d, const geometric_volume &where);
  flux_vol *add_flux_plane(const geometric_volume &where);
  flux_vol *add_flux_plane(const vec &p1, const vec &p2);
  double electric_energy_max_in_box(const geometric_volume &where);
  double modal_volume_in_box(const geometric_volume &where);
  double electric_sqr_weighted_integral(double (*deps)(const vec &),
				       const geometric_volume &where);
  double electric_energy_weighted_integral(double (*f)(const vec &),
					   const geometric_volume &where);

  void set_output_directory(const char *name);
  void verbose(int v=1);
  double count_volume(component);
  // polarization.cpp
  void initialize_polarization_energy(const polarizability_identifier &,
                                      double energy(const vec &));

  // fields.cpp
  bool have_component(component);
  void set_rshift(double rshift);
  // material.cpp
  double max_eps() const;
  // step.cpp
  void force_consistency(field_type ft);

  bool nosize_direction(direction d) const;
  direction normal_direction(const geometric_volume &where) const;

 private: 
  int verbosity; // Turn on verbosity for debugging purposes...
  double last_wall_time;
  time_sink working_on, was_working_on;
  double times_spent[Other+1];
  // time.cpp
  void am_now_working_on(time_sink);
  void finished_working();
  // boundaries.cpp
  bool chunk_connections_valid;
  void find_metals();
  void disconnect_chunks();
  void connect_chunks();
  void connect_the_chunks(); // Intended to be ultra-private...
  bool on_metal_boundary(const ivec &);
  ivec ilattice_vector(direction) const;
  bool locate_component_point(component *, ivec *, complex<double> *) const;
  bool locate_point_in_user_volume(ivec *, complex<double> *phase) const;
  void locate_volume_source_in_user_volume(const vec p1, const vec p2, vec newp1[8], vec newp2[8],
                                           complex<double> kphase[8], int &ncopies) const;
  // step.cpp
  void phase_material();
  void step_h();
  void step_h_source();
  void step_d();
  void step_d_source();
  void update_e_from_d();
  void update_from_e();
  void step_boundaries(field_type);
  void calc_sources(double tim);
  int cluster_some_bands_cleverly(double *tf, double *td, complex<double> *ta,
                                  int num_freqs, int fields_considered, int maxbands,
                                  complex<double> *fad, double *approx_power);
  void out_bands(FILE *, const char *, int maxbands);
  complex<double> *clever_cluster_bands(int maxbands, double *approx_power = NULL);
  // monitor.cpp
  complex<double> get_field(component c, const ivec &iloc) const;
  double get_polarization_energy(const ivec &) const;
  double get_polarization_energy(const vec &) const;
  double get_polarization_energy(const polarizability_identifier &, const ivec &) const;
  double get_polarization_energy(const polarizability_identifier &, const vec &) const;
  double get_inveps(component, direction, const ivec &iloc) const;
  double get_eps(const ivec &iloc) const;
};

class flux_vol {
 public:
  flux_vol(fields *f_, direction d_, const geometric_volume &where_) : where(where_) {
    f = f_; d = d_; cur_flux = cur_flux_half = 0; 
    next = f->fluxes; f->fluxes = this;
  }
  ~flux_vol() { delete next; }

  void update_half() { cur_flux_half = flux_wrongE(); 
                       if (next) next->update_half(); }
  void update() { cur_flux = (flux_wrongE() + cur_flux_half) * 0.5;
                  if (next) next->update(); }

  double flux() { return cur_flux; }

  flux_vol *next;
 private:
  double flux_wrongE() { return f->flux_in_box_wrongH(d, where); }
  fields *f;
  direction d;
  geometric_volume where;
  double cur_flux, cur_flux_half;
};

class grace_point;
enum grace_type { XY, ERROR_BARS };

class grace {
 public:
  grace(const char *fname, const char *dirname = ".");
  ~grace();
  
  void new_set(grace_type t = XY);
  void new_curve();
  void set_legend(const char *);
  void set_range(double xmin, double xmax, double ymin, double ymax);
  void output_point(double x, double y,
                    double dy = -1.0, double extra = -1.0);
  void output_out_of_order(int n, double x, double y,
                           double dy = -1.0, double extra= -1.0);
 private:
  void flush_pts();
  FILE *f;
  char *fn, *dn;
  grace_point *pts;
  int set_num,sn;
};

// The following is a utility function to parse the executable name use it
// to come up with a directory name, avoiding overwriting any existing
// directory, unless the source file hasn't changed.

const char *make_output_directory(const char *exename, const char *jobname = NULL);
void trash_output_directory(const char *dirname);
FILE *create_output_file(const char *dirname, const char *fname);

// The following allows you to hit ctrl-C to tell your calculation to stop
// and clean up.
void deal_with_ctrl_c(int stop_now = 2);
// When a ctrl_c is called, the following variable (which starts with a
// zero value) is incremented.
extern int interrupt;

int do_harminv(complex<double> *data, int n, double dt,
	       double fmin, double fmax, int maxbands,
	       complex<double> *amps, double *freq_re, double *freq_im,
	       double *errors = NULL,
	       double spectral_density = 1.1, double Q_thresh = 50,
	       double rel_err_thresh = 1e20, double err_thresh = 0.01, 
	       double rel_amp_thresh = -1, double amp_thresh = -1);

} /* namespace meep */

#endif /* MEEP_H */