monitor.cpp

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

                                             const ivec &iloc) const {
  double res = 0.0;
  polarization *p = pol;
  while (is_mine() && p) {
    if (p->pb->get_identifier() == pi) res += p->local_energy(iloc);
    p = p->next;
  }
  return broadcast(n_proc(), res);
}

double fields_chunk::my_polarization_energy(const polarizability_identifier &pi,
                                            const ivec &iloc) const {
  if (!is_mine())
    abort("Can't call my_polarization_energy on someone else's chunk!\n");
  double res = 0.0;
  polarization *p = pol;
  while (p) {
    if (p->pb->get_identifier() == pi) res += p->local_energy(iloc);
    p = p->next;
  }
  return res;
}

double fields::get_eps(const ivec &origloc) const {
  ivec iloc = origloc;
  complex<double> aaack = 1.0;
  locate_point_in_user_volume(&iloc, &aaack);
  for (int sn=0;sn<S.multiplicity();sn++)
    for (int i=0;i<num_chunks;i++)
      if (chunks[i]->v.contains(S.transform(iloc,sn)))
        return chunks[i]->get_eps(S.transform(iloc,sn));
  return 0.0;
}

double fields_chunk::get_eps(const ivec &iloc) const {
  double res = 0.0;
  if (is_mine()) res = s->eps[v.index(v.eps_component(), iloc)];
  return broadcast(n_proc(), res);
}

double fields::get_inveps(component c, direction d,
			  const ivec &origloc) const {
  ivec iloc = origloc;
  complex<double> aaack = 1.0;
  locate_point_in_user_volume(&iloc, &aaack);
  for (int sn=0;sn<S.multiplicity();sn++)
    for (int i=0;i<num_chunks;i++)
      if (chunks[i]->v.contains(S.transform(iloc,sn))) {
	signed_direction ds = S.transform(d,sn);
        return chunks[i]->get_inveps(S.transform(c,sn), ds.d,
				     S.transform(iloc,sn))
	     * (ds.flipped ^ S.transform(component_direction(c),sn).flipped
		? -1 : 1);
      }
  return 0.0;
}

double fields_chunk::get_inveps(component c, direction d,
			     const ivec &iloc) const {
  double res = 0.0;
  if (is_mine()) res = s->inveps[c][d] ? s->inveps[c][d][v.index(c, iloc)] : 0;
  return broadcast(n_proc(), res);
}

double fields::get_inveps(component c, direction d,
			  const vec &loc) const {
  ivec ilocs[8];
  double w[8];
  double val[8];
  for (int i=0;i<8;i++) val[i] = 0.0;
  v.interpolate(c, loc, ilocs, w);
  for (int argh=0;argh<8&&w[argh];argh++)
    val[argh] = w[argh]*get_inveps(c,d,ilocs[argh]);
  dumbsort(val);
  double res = 0.0;
  for (int i=0;i<8;i++) res += val[i];
  return res;
}

double fields::get_eps(const vec &loc) const {
  double tr = 0;
  int nc = 0;
  FOR_ELECTRIC_COMPONENTS(c)
    if (v.has_field(c)) {
      tr += get_inveps(c, component_direction(c), loc);
      ++nc;
    }
  return nc / tr;
}

double structure::get_inveps(component c, direction d,
			     const ivec &origloc) const {
  ivec iloc = origloc;
  for (int sn=0;sn<S.multiplicity();sn++)
    for (int i=0;i<num_chunks;i++)
      if (chunks[i]->v.contains(S.transform(iloc,sn))) {
	signed_direction ds = S.transform(d,sn);
        return chunks[i]->get_inveps(S.transform(c,sn), ds.d,
				     S.transform(iloc,sn))
	     * (ds.flipped ^ S.transform(component_direction(c),sn).flipped
		? -1 : 1);
      }
  return 0.0;
}

double structure_chunk::get_inveps(component c, direction d,
				   const ivec &iloc) const {
  double res = 0.0;
  if (is_mine()) res = inveps[c][d] ? inveps[c][d][v.index(c, iloc)] : 0;
  return broadcast(n_proc(), res);
}

double structure::get_inveps(component c, direction d,
			     const vec &loc) const {
  ivec ilocs[8];
  double w[8];
  double val[8];
  for (int i=0;i<8;i++) val[i] = 0.0;
  v.interpolate(c, loc, ilocs, w);
  for (int argh=0;argh<8&&w[argh];argh++)
    val[argh] = w[argh]*get_inveps(c,d,ilocs[argh]);
  dumbsort(val);
  double res = 0.0;
  for (int i=0;i<8;i++) res += val[i];
  return res;
}

double structure::get_eps(const vec &loc) const {
  double tr = 0;
  int nc = 0;
  FOR_ELECTRIC_COMPONENTS(c)
    if (v.has_field(c)) {
      tr += get_inveps(c, component_direction(c), loc);
      ++nc;
    }
  return nc / tr;
}

monitor_point *fields::get_new_point(const vec &loc, monitor_point *the_list) const {
  monitor_point *p = new monitor_point();
  get_point(p, loc);
  p->next = the_list;
  return p;
}

complex<double> monitor_point::get_component(component w) {
  return f[w];
}
  
double monitor_point::poynting_in_direction(direction d) {
  int direction_shift = 0;
  if (loc.dim == Dcyl)
    direction_shift = 2;

  direction d1 = (direction) (((d+1) - direction_shift)%3 + direction_shift);
  direction d2 = (direction) (((d+2) - direction_shift)%3 + direction_shift);

  // below Ex and Hx are used just to say that we want electric or magnetic component
  complex<double> E1 = get_component(direction_component(Ex, d1));
  complex<double> E2 = get_component(direction_component(Ex, d2));
  complex<double> H1 = get_component(direction_component(Hx, d1));
  complex<double> H2 = get_component(direction_component(Hx, d2));

  return (real(E1)*real(H2) - real(E2)*real(H1)) + (imag(E1)*imag(H2) - imag(E2)*imag(H1));
}

double monitor_point::poynting_in_direction(vec v) {
  if (v.dim != loc.dim)
    abort("poynting_in_direction: v.dim != loc.dim\n");
  v = v / abs(v);
  double result = 0.0;
  LOOP_OVER_DIRECTIONS(v.dim, d)
    result += v.in_direction(d) * poynting_in_direction(d);
  return result;
}

void monitor_point::fourier_transform(component w,
                                      complex<double> **a, complex<double> **f,
                                      int *numout, double fmin, double fmax,
                                      int maxbands) {
  int n = 1;
  monitor_point *p = next;
  double tmax = t, tmin = t;
  while (p) {
    n++;
    if (p->t > tmax) tmax = p->t;
    if (p->t < tmin) tmin = p->t;
    p = p->next;
  }
  p = this;
  complex<double> *d = new complex<double>[n];
  for (int i=0;i<n;i++,p=p->next) {
    d[i] = p->get_component(w);
  }
  if (fmin > 0.0) { // Get rid of any static fields_chunk!
    complex<double> mean = 0.0;
    for (int i=0;i<n;i++) mean += d[i];
    mean /= n;
    for (int i=0;i<n;i++) d[i] -= mean;
  }
#if HAVE_SOME_FFTW
  if ((fmin > 0.0 || fmax > 0.0) && maxbands > 0) {
#else
  if ((fmin <= 0.0 && fmax <= 0.0) || maxbands <= 0) {
    maxbands = n;
    fmin = 0; fmax = (n-1)*(1.0/(tmax-tmin));
  }
#endif
    *a = new complex<double>[maxbands];
    *f = new complex<double>[maxbands];
    *numout = maxbands;
    delete[] d;
    for (int i = 0;i<maxbands;i++) {
      double df = (maxbands == 1) ? 0.0 : (fmax-fmin)/(maxbands-1);
      (*f)[i] = fmin + i*df;
      (*a)[i] = 0.0;
      p = this;
      while (p) {
        double inside = 2*pi*real((*f)[i])*p->t;
        (*a)[i] += p->get_component(w)*complex<double>(cos(inside),sin(inside));
        p = p->next;
      }
      (*a)[i] /= (tmax-tmin);
    }
#if HAVE_SOME_FFTW
  } else {
    *numout = n;
    *a = new complex<double>[n];
    *f = d;
    fftw_complex *in = (fftw_complex *) d, *out = (fftw_complex *) *a;
    fftw_plan p;
#ifdef HAVE_LIBFFTW3
    p = fftw_plan_dft_1d(n, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(p);
    fftw_destroy_plan(p);
#else
    p = fftw_create_plan(n, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_one(p, in, out);
    fftw_destroy_plan(p);
#endif
    for (int i=0;i<n;i++) {
      (*f)[i] = i*(1.0/(tmax-tmin));
      if (real((*f)[i]) > 0.5*n/(tmax-tmin)) (*f)[i] -= n/(tmax-tmin);
      (*a)[i] *= (tmax-tmin)/n;
    }
  }
#endif
}

void monitor_point::harminv(component w,
                            complex<double> **a, complex<double> **f,
                            int *numout, double fmin, double fmax,
                            int maxbands) {
  int n = 1;
  monitor_point *p = next;
  double tmax = t, tmin = t;
  while (p) {
    n++;
    if (p->t > tmax) tmax = p->t;
    if (p->t < tmin) tmin = p->t;
    p = p->next;
  }
  p = this;
  complex<double> *d = new complex<double>[n];
  for (int i=0;i<n;i++,p=p->next) {
    d[i] = p->get_component(w);
  }
  *a = new complex<double>[n];
  double *f_re = new double[n];
  double *f_im = new double[n];
  *numout = do_harminv(d, n, (tmax-tmin)/(n-1), fmin, fmax, maxbands,
                       *a, f_re, f_im, NULL);
  *f = new complex<double>[*numout];
  for (int i=0;i<*numout;i++)
    (*f)[i] = complex<double>(f_re[i],f_im[i]);
  delete[] f_re;
  delete[] f_im;
  delete[] d;
}

} // namespace meep