slices.cpp

采用FDTD时域有限差分法计算电磁波的传播问题等。

        ivec next = v.iloc(m,i)+iveccyl(2,0);
        vec nextrot = v[S.transform(next - iveccyl(1,0),symnum)];
        if (v.contains(next) && f[i] != f[v.index(m,next)])
          buf.printf("%g\t%g\tLH\n", nextrot.z(), nextrot.r());
        next = v.iloc(m,i)+iveccyl(0,2);
        nextrot = v[S.transform(next - iveccyl(0,1),symnum)];
        if (v.contains(next) && f[i] != f[v.index(m,next)])
          buf.printf("%g\t%g\tLV\n", nextrot.z(), nextrot.r());
        break;
      }
      case D2: {
        ivec next = v.iloc(m,i)+ivec(0,2);
        vec nextrot = v[(S.transform(next - ivec(0,1),symnum))];
        if (v.owns(next))
          if (f[i] != f[v.index(m,next)])
            buf.printf("%g\t%g\tLH\n", nextrot.x(), nextrot.y());
        next = v.iloc(m,i)-ivec(0,2);
        nextrot = v[S.transform(next + ivec(0,1),symnum)];
        if (v.owns(next))
          if (f[i] != f[v.index(m,next)])
            buf.printf("%g\t%g\tLH\n", nextrot.x(), nextrot.y());
        next = v.iloc(m,i)+ivec(2,0);
        nextrot = v[S.transform(next - ivec(1,0),symnum)];
        if (v.owns(next))
          if (f[i] != f[v.index(m,next)])
            buf.printf("%g\t%g\tLV\n", nextrot.x(), nextrot.y());
        next = v.iloc(m,i)-ivec(2,0);
        nextrot = v[S.transform(next + ivec(1,0),symnum)];
        if (v.owns(next))
          if (f[i] != f[v.index(m,next)])
            buf.printf("%g\t%g\tLV\n", nextrot.x(), nextrot.y());
        break;
      }
      case D1: {
        const ivec next = v.iloc(m,i)+ivec(2);
        const vec nextrot = v[S.transform(next - ivec(1),symnum)];
        if (v.contains(next) && f[i] != f[v.index(m,next)])
          buf.printf("%g\t%g\tLV\n", nextrot.z(), 0.0);
        break;
      }
      case D3: abort("Error in eps_outline.\n"); break;
      }
  }
}

static void output_complex_eps_body(component m, double *f[2], const volume &v,
                                    symmetry S, int symnum,
                                    const geometric_volume &what, file *out) {
  bufprint buf(out);
  if (!f[0]) return; // Field doesn't exist...
  const complex<double> ph = S.phase_shift(m, symnum);
  for (int i=0;i<v.ntot();i++) {
    const vec here = S.transform(v.loc(m,i),symnum);
    if (what.contains(here)) {
      double x = 0, y = 0;
      switch (v.dim) {
      case Dcyl: x = here.z(); y = here.r(); break;
      case D1: x = here.z(); break;
      case D2: x = here.x(); y = here.y(); break;
      case D3: x = here.x(); y = here.y(); break; // FIXME use right directions!
      }
      if (f[1]) buf.printf("%g\t%g\t%g\tP\n", x, y,
                         real(ph)*f[0][i] - imag(ph)*f[1][i]);
      else buf.printf("%g\t%g\t%g\tP\n", x, y, real(ph)*f[0][i]);
    }
  }
}

void fields_chunk::output_eps_body(component c, const symmetry &S, int sn,
                                   const geometric_volume &what, file *out,
                                   complex<double> phshift) {
  bufprint buf(out);
  if (f[c][0]) {
    if (v.a <= default_eps_resolution) {
      output_complex_eps_body(c, f[c], v, S, sn, what, out);
    } else {
      const complex<double> ph = S.phase_shift(c, sn);
      const int n = v.ntot_at_resolution(default_eps_resolution);
      for (int i=0;i<n;i++) {
        const vec here = v.loc_at_resolution(i,default_eps_resolution);
        const vec there = S.transform(here, sn);
        if (what.contains(here)) {
          double x = 0, y = 0;
          switch (v.dim) {
          case Dcyl: x = there.z(); y = there.r(); break;
          case D1: x = there.z(); break;
          case D2: x = there.x(); y = there.y(); break;
          case D3: abort("Don't support 3D o_e_b\n"); break;
          }
          ivec ilocs[8];
          double w[8];
          double fhere = 0.0;
          v.interpolate(c, here, ilocs, w);
          for (int i=0;i<8&&w[i];i++)
            if (v.contains(S.transform(ilocs[i],sn)))
              fhere += w[i]*real(phshift*get_field(c,ilocs[i]));
          if (fhere != 0.0) // save space by leaving out blanks.
            buf.printf("%g\t%g\t%g\tP\n", x, y, fhere);
        }
      }
    }
  }
}

void fields_chunk::output_eps_body(const polarizability_identifier &p,
                                   component c, const symmetry &S, int sn,
                                   const geometric_volume &what, file *out,
                                   complex<double> phshift) {
  bufprint buf(out);
  if (f[c][0]) {
    if (v.a <= default_eps_resolution) {
      output_complex_eps_body(c, f[c], v, S, sn, what, out);
    } else {
      const complex<double> ph = S.phase_shift(c, sn);
      const int n = v.ntot_at_resolution(default_eps_resolution);
      for (int i=0;i<n;i++) {
        const vec here = v.loc_at_resolution(i,default_eps_resolution);
        const vec there = S.transform(here, sn);
        if (what.contains(here)) {
          double x = 0, y = 0;
          switch (v.dim) {
          case Dcyl: x = there.z(); y = there.r(); break;
          case D1: x = there.z(); break;
          case D2: x = there.x(); y = there.y(); break;
          case D3: abort("Don't support 3D o_e_b\n"); break;
          }
          ivec ilocs[8];
          double w[8];
          double fhere = 0.0;
          v.interpolate(c, here, ilocs, w);
          for (int i=0;i<8&&w[i];i++)
            if (v.contains(S.transform(ilocs[i],sn)))
              fhere += w[i]*real(phshift*get_polarization_field(p,c,ilocs[i]));
          if (fhere != 0.0) // save space by leaving out blanks.
            buf.printf("%g\t%g\t%g\tP\n", x, y, fhere);
        }
      }
    }
  }
}

static void output_complex_eps_header(component m, double fmax, const volume &v,
                                      const geometric_volume &what, file *out,
                                      const char *name, component om = Hx) {
  double xmin = 1e300, ymin = 1e300, xmax = -1e300, ymax = -1e300;
  direction xdir, ydir;
  switch (v.dim) {
  case Dcyl: xdir = Z; ydir = R; break;
  case D3: xdir = X; ydir = Y; break; // FIXME: check the thin direction of what.
  case D2: xdir = X; ydir = Y; break;
  case D1: xdir = Z; ydir = Z; break; // a tad ugly...
  }
  xmin = what.in_direction_min(xdir);
  xmax = what.in_direction_max(xdir);
  ymin = what.in_direction_min(ydir);
  ymax = what.in_direction_max(ydir);
  if (v.dim == D1) { ymin = -v.inva*0.5; ymax = -ymin; }
  if (ymax == ymin) ymax = ymin + v.inva;
  if (fmax == 0.0) fmax = 0.0001;
  if (v.dim == D1) {
    // Make a 1D line plot!
    eps_1d_header(xmin, ymin, xmax, ymax, fmax, v.a, out, name);
  } else {
    // Make a 2D color plot!
    eps_header(xmin, ymin, xmax, ymax, fmax, v.a, out, name);
  }
}

static void output_eps_header(double fmax, double dx,
                              const double xmin, const double xmax,
                              const double ymin, const double ymax,
                              file *out, const char *name) {
  if (fmax == 0.0) fmax = 0.0001;
  // Make a 2D color plot!
  eps_header(xmin, ymin, xmax, ymax, fmax, 1.0/dx, out, name);
}

static void output_complex_eps_tail(file *out) {
  eps_trailer(out);
}

void structure::output_slices(const char *name) const {
  output_slices(v.surroundings(),name);
}
void structure::output_slices(const geometric_volume &what, const char *name) const {
  const int buflen = 1024;
  char nname[buflen];
  if (*name) snprintf(nname, buflen, "%s-", name);
  else *nname = 0; // No additional name!
  char *n = new char[buflen];
  if (!n) abort("Allocation failure!\n");
  snprintf(n, buflen, "%s/%sepsilon.sli", outdir, nname);
  file *out = everyone_open_write(n);
  if (!out) {
    printf("Unable to open file '%s' for slice output.\n", n);
    return;
  }
  for (int i=0;i<num_chunks;i++)
    if (chunks[i]->is_mine())
      output_slice(v.eps_component(), chunks[i]->eps, chunks[i]->v, what, out);
  everyone_close(out);

  FOR_DIRECTIONS(d)
    FOR_COMPONENTS(c) {
      snprintf(n, buflen, "%s/%ssigma-C-%s-%s.sli", outdir, nname, component_name(c), direction_name(d));
      file *out = everyone_open_write(n);
      if (!out) {
	printf("Unable to open file '%s' for slice output.\n", n);
	return;
      }
      for (int i=0;i<num_chunks;i++)
	if (chunks[i]->is_mine())
	  output_slice(c, chunks[i]->C/*decay*/[d][c]/*[component_direction(c)]*/, chunks[i]->v, what, out);
      everyone_close(out);
    }
  delete[] n;
}

void fields::output_real_imaginary_slices(const char *name) {
  output_real_imaginary_slices(v.surroundings(),name);
}

void fields::output_real_imaginary_slices(const geometric_volume &what,
                                          const char *name) {
  const int buflen = 1024;
  char nname[buflen];
  if (*name) snprintf(nname, buflen, "%s-", name);
  else *nname = 0; // No additional name!
  char *n = (char *)malloc(buflen);
  if (!n) abort("Allocation failure!\n");
  char *r_or_i = "-re";
  DOCMP {
    FOR_COMPONENTS(c)
      if (v.has_field(c)) {
        snprintf(n, buflen, "%s/%s%s%s-%09.2f.sli",
                 outdir, nname, component_name(c),
                 r_or_i, time());
        file *out = everyone_open_write(n);
        if (!out) {
          printf("Unable to open file '%s' for slice output.\n", n);
          return;
        }
        for (int i=0;i<num_chunks;i++)
          if (chunks[i]->is_mine())
            output_slice(c, chunks[i]->f[c][cmp],
                         chunks[i]->v, what, out);
        everyone_close(out);
      }
    r_or_i = "-im";
  }

  free(n);
}

void fields::output_slices(const char *name) {
  output_slices(v.surroundings(), name);
}
void fields::output_slices(const geometric_volume &what, const char *name) {
  am_now_working_on(FieldOutput);
  const int buflen = 1024;
  char nname[buflen];
  if (*name) snprintf(nname, buflen, "%s-", name);
  else *nname = 0; // No additional name!
  char *n = (char *)malloc(buflen);
  if (!n) abort("Allocation failure!\n");
  char time_step_string[buflen];
  if (a == 1)
    snprintf(time_step_string, buflen, "%08.0f", time());
  else
    snprintf(time_step_string, buflen, "%09.2f", time());
  FOR_COMPONENTS(c)
    if (v.has_field(c)) {
      snprintf(n, buflen, "%s/%s%s-%s.sli", outdir, nname,
               component_name(c), time_step_string);
      file *out = everyone_open_write(n);
      if (!out) {
        printf("Unable to open file '%s' for slice output.\n", n);
        return;
      }
      complex<double> phshift = optimal_phase_shift(c);
      for (int i=0;i<num_chunks;i++)
        if (chunks[i]->is_mine())
          output_complex_slice(c, chunks[i]->f[c], phshift,
                               chunks[i]->v, what, out);
      everyone_close(out);
    }
  //if (new_s) {
  //  snprintf(n, buflen, "%s/%sepsilon-%s.sli", outdir, nname, time_step_string);
  //  output_slice(v.eps_component(), s->eps, v, what, n);
  //}
  free(n);
  finished_working();
}

void fields::eps_energy_slice(const char *name) {
  if (v.dim == D3) abort("FIXME need to support 3D energy slices...\n");
  geometric_volume what = user_volume.surroundings();
  if (v.dim == Dcyl) what.set_direction_min(R,-what.in_direction_max(R));
  eps_energy_slice(what,name);
}

void fields::eps_energy_slice(const geometric_volume &what, const char *name) {
  am_now_working_on(FieldOutput);
  const int buflen = 1024;
  char nname[buflen];
  if (*name) snprintf(nname, buflen, "%s-", name);
  else *nname = 0; // No additional name!
  char *n = new char[buflen];
  if (!n) abort("Allocation failure!\n");
  char time_step_string[buflen];
  snprintf(time_step_string, buflen, "%09.2f", time());
  
  snprintf(n, buflen, "%s/%senergy-%s.eps", outdir, nname, time_step_string);
  file *out = everyone_open_write(n);
  if (!out) abort("Unable to open file '%s' for slice output.\n", n);
  const double fmax = max(maxpolenergy_to_master(), -minpolenergy_to_master());
  if (am_master())
    output_complex_eps_header(v.eps_component(), fmax, user_volume,
                              what, out, n, v.eps_component());
  if (v.dim != D1) abort("Still only works in 1D.  :( \n");
  {
    bufprint buf(out);
    for (int i=0;i<num_chunks;i++)
      if (chunks[i]->is_mine())
        for (int sn=0;sn<S.multiplicity();sn++)
          for (int n=0;n<chunks[i]->v.ntot();n++) {
            const ivec here = chunks[i]->v.iloc(v.eps_component(), n);
            buf.printf("%g\t0\t%g\tP\n", chunks[i]->v[S.transform(here,sn)].z(),
                       chunks[i]->get_polarization_energy(here));
          }
  }
  all_wait();
  for (int i=0;i<num_chunks;i++)
    if (chunks[i]->is_mine())
      for (int sn=0;sn<S.multiplicity();sn++)
        for (int otherc=0;otherc<10;otherc++)
          if (S.transform((component)otherc,sn) == 0.5) // fixme -- what is this?
            eps_outline(v.eps_component(), chunks[i]->s->eps,
                        chunks[i]->v, what, S, sn, out);
  all_wait();
  if (am_master()) output_complex_eps_tail(out);
  everyone_close(out);
  delete[] n;
  finished_working();
}

void fields::eps_energy_slice(const polarizability_identifier &p, const char *name) {
  if (v.dim == D3) abort("FIXME need to support 3D energy slices...\n");
  geometric_volume what = user_volume.surroundings();
  if (v.dim == Dcyl) what.set_direction_min(R,-what.in_direction_max(R));
  eps_energy_slice(p, what,name);
}

void fields::eps_energy_slice(const polarizability_identifier &p,
                              const geometric_volume &what, const char *name) {
  am_now_working_on(FieldOutput);