slices.cpp

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

/* Copyright (C) 2006 Massachusetts Institute of Technology
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2, or (at your option)
%  any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software Foundation,
%  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdarg.h>

#include "meep.hpp"
#include "meep_internals.hpp"

namespace meep {

complex<double> fields::optimal_phase_shift(component c) const {
  complex<double> mean = field_mean(c, false, false);
  complex<double> meanr= field_mean(c, true , false);
  complex<double> meani= field_mean(c, false, true );
  if (abs(mean) > abs(meanr) && abs(mean) > abs(meani)) return abs(mean)/mean;
  if (abs(meanr) > abs(meani)) return abs(meanr)/meanr;
  if (abs(meani) > 0.0) return abs(meani)/meani;
  return 1.0;
}

complex<double> fields::field_mean(component c, bool abs_real,
                                   bool abs_imag) const {
  complex<double> themean = 0.0;
  for (int i=0;i<num_chunks;i++)
    themean += chunks[i]->field_mean(c, abs_real, abs_imag);
  return sum_to_all(themean);
}

complex<double> fields_chunk::field_mean(component c, bool abs_real,
                                         bool abs_imag) const {
  complex<double> themean = 0.0;
  if (f[c][0]) for (int i=0;i<v.ntot();i++)
    themean += (v.dV(c,i) & gv).full_volume()*
      (abs_real)?fabs(f[c][0][i]):f[c][0][i];
  if (f[c][1]) for (int i=0;i<v.ntot();i++)
    themean += (v.dV(c,i) & gv).full_volume()*
      complex<double>(0,abs_imag?fabs(f[c][1][i]):f[c][1][i]);
  return themean;
}

const int bufsize = 1<<16;

class bufprint {
public:
  bufprint(file *out);
  ~bufprint();
  void flushme();
  void printf(char *fmt, ...);
private:
  file *f;
  char *buf;
  int inbuf;
};

bufprint::bufprint(file *out) {
  f = out;
  buf = new char[bufsize];
  if (!buf) abort("not enough memory for a buffer!");
  inbuf = 0;
}

bufprint::~bufprint() {
  flushme();
  delete[] buf;
}

void bufprint::flushme() {
  if (inbuf) {
    buf[inbuf] = 0;
    i_fprintf(f, "%s", buf);
    inbuf = 0;
  }
}

void bufprint::printf(char *fmt, ...) {
  va_list ap;
  va_start(ap, fmt);
  int written = vsnprintf(buf + inbuf, bufsize - inbuf - 1, fmt, ap);
  if (written <= 0 || written > bufsize - inbuf - 1) {
    flushme();
    written = vsnprintf(buf + inbuf, bufsize - inbuf - 1, fmt, ap);
  }
  inbuf += written;
  if (inbuf > bufsize/2) flushme();
}

static void output_complex_slice(component m, double *f[2],
                                 complex<double> phshift, const volume &v,
                                 const geometric_volume &what, file *out) {
  if (!f[0] || ! f[1]) return; // Field doesn't exist...
  double c = real(phshift), s = imag(phshift);
  bufprint buf(out);
  for (int i=0;i<v.ntot();i++) {
    vec here = v.loc(m,i);
    if (what.contains( here)) {
      if (v.dim == Dcyl) {
        buf.printf("%lg\t%lg\t0\t%lg\n", here.r(), here.z(), c*f[0][i]+s*f[1][i]);
      } else if (v.dim == D3  ) {
        buf.printf("%lg\t%lg\t%lg\t%lg\n", here.x(), here.y(), here.z(), c*f[0][i]+s*f[1][i]);
      } else if (v.dim == D2  ) {
        buf.printf("%lg\t%lg\t0\t%lg\n", here.x(), here.y(), c*f[0][i]+s*f[1][i]);
      } else if (v.dim == D1  ) {
        buf.printf("0\t0\t%lg\t%lg\n", here.z(), c*f[0][i]+s*f[1][i]);
      }
    }
  }  
}

static void output_slice(component m, const double *f, const volume &v,
                         const geometric_volume &what, file *out) {
  if (!f) return; // Field doesn't exist...
  bufprint buf(out);
  for (int i=0;i<v.ntot();i++) {
    vec here = v.loc(m,i);
    if (what.contains(here)) {
      if (v.dim == Dcyl) {
        buf.printf("%lg\t%lg\t0\t%lg\n", here.r(), here.z(), f[i]);
      } else if (v.dim == D3  ) {
        buf.printf("%lg\t%lg\t%lg\t%lg\n", here.x(), here.y(), here.z(), f[i]);
      } else if (v.dim == D2  ) {
        buf.printf("%lg\t%lg\t0\t%lg\n", here.x(), here.y(), f[i]);
      } else if (v.dim == D1  ) {
        buf.printf("0\t0\t%lg\t%lg\n", here.z(), f[i]);
      }
    }
  }
}

static const double default_eps_resolution = 20.0;
static const double default_eps_size = 1000.0;

static void eps_header(double xmin, double ymin, double xmax, double ymax,
                       double fmax, double a, file *out, const char *name) {
  bufprint buf(out);
  buf.printf("%%!PS-Adobe-3.0 EPSF\n");
  double size = xmax - xmin + ymax - ymin;
  buf.printf("%%%%BoundingBox: 0 0 %g %g\n",
           (xmax-xmin)*default_eps_size/size, (ymax-ymin)*default_eps_size/size);
  buf.printf("gsave\n");
  buf.printf("gsave\n");
  buf.printf("/title (%s) def\n", name);
  buf.printf("%g %g scale\n", default_eps_size/size, default_eps_size/size);
  buf.printf("%g %g translate\n", -xmin, -ymin);
  buf.printf("/max %g def\n", fmax);
  const double dx = 1.0/min(a, default_eps_resolution);
  buf.printf("/dx %g def\n", dx);
  buf.printf("/hdx %g def\n", dx*0.5);
  buf.printf("dx 10 div setlinewidth\n");
  buf.printf("/dotrad %g def\n", dx*0.25);
  buf.printf("1 setlinecap\n");
  buf.printf("/P {\n\
    max div\n\
    dup 0 lt {\n\
        1 add\n\
        dup 1\n\
    }{\n\
        neg 1 add\n\
        dup 1 3 1 roll\n\
    } ifelse\n\
    setrgbcolor\n\
    newpath\n\
    moveto\n");
  buf.printf("    %g %g rmoveto\n", dx*0.5, dx*0.5);
  buf.printf("    0 %g rlineto\n", -dx);
  buf.printf("    %g 0 rlineto\n", -dx);
  buf.printf("    0 dx rlineto\n", dx);
  buf.printf("    dx 0 rlineto\n", dx);
  buf.printf("    gsave\n\
    fill\n\
    grestore\n");
  buf.printf("    %g setlinewidth\n", dx*0.1);
  buf.printf("    stroke\n\
} def\n\
/LV {\n\
    0 0 0 setrgbcolor\n\
    moveto\n");
  buf.printf("    %g setlinewidth\n", dx*0.1);
  buf.printf("    0 %g rmoveto\n", dx*0.5);
  buf.printf("    0 %g rlineto\n", -dx);
  buf.printf("    stroke\n\
} def\n\
/LH {\n");
  buf.printf("    %g setlinewidth\n", dx*0.1);
  buf.printf("    0 0 0 setrgbcolor\n\
    moveto\n");
  buf.printf("    %g 0 rmoveto\n", dx*0.5);
  buf.printf("    %g 0 rlineto\n", -dx);
  buf.printf("    stroke\n\
} def\n");
  buf.printf("    /DV { [0 %g] 0 setdash LV } def\n", dx/4);
  buf.printf("    /DH { [0 %g] 0 setdash LH } def\n", dx/4);
  buf.printf("    /D { moveto\n\
    %g setlinewidth\n\
    0 0.8 0 setrgbcolor [0 %g] 0 setdash\n\
    lineto stroke } def\n", 0.6*dx, 3*dx);
  // B for Boundary...
  buf.printf("/B {\n\
    0 0 0 setrgbcolor \n\
    newpath dotrad 0 360 arc fill \n\
} def\n");
}

static void eps_1d_header(double xmin, double ymin, double xmax, double ymax,
                          double fmax, double a, file *out, const char *name) {
  bufprint buf(out);
  buf.printf("%%!PS-Adobe-3.0 EPSF\n");
  const double size = xmax - xmin;
  const double fsize = (5.0 < 0.2*size)?0.2*size:5.0;
  const double dx = 1.0/a;
  buf.printf("%%%%BoundingBox: 0 0 %g %g\n",
           (xmax-xmin)*default_eps_size/size, default_eps_size*fsize/size);
  buf.printf("gsave\n");
  buf.printf("gsave\n");
  buf.printf("/title (%s) def\n", name);
  buf.printf("%g %g scale\n", default_eps_size/size, default_eps_size/size);
  buf.printf("%g %g translate\n", -xmin, 0.5*fsize);
  buf.printf("%g 0 moveto %g 0 lineto %g setlinewidth stroke\n",
            xmin, xmax, dx*0.1); 
  buf.printf("/height %g def\n", (default_eps_size*0.5*fsize)/size);
  buf.printf("/max %g def\n", fmax);
  buf.printf("/fscale %g def\n", 2.2*fmax/fsize);
  buf.printf("/dotrad %g def\n", dx);
  buf.printf("/dx %g def\n", dx);
  buf.printf("/hdx %g def\n", dx*0.5);
  buf.printf("dx 10 div setlinewidth\n");
  buf.printf("1 setlinecap\n");
  buf.printf("/P {\n\
    fscale div exch pop \n\
    newpath dotrad 0 360 arc fill \n\
} def\n");
  buf.printf("/LV {\n\
    0.8 0.8 0 setrgbcolor\n\
    pop dup height moveto height neg lineto\n");
  buf.printf("    %g setlinewidth\n", dx*0.5);
  buf.printf("    stroke\n\
} def\n\
/LH {\n");
  buf.printf("    %g setlinewidth\n", dx*0.1);
  buf.printf("    0 0 0 setrgbcolor\n\
    moveto\n");
  buf.printf("    %g 0 rmoveto\n", 10*dx*0.5);
  buf.printf("    %g 0 rlineto\n", -10*dx);
  buf.printf("    stroke\n\
} def\n");
  buf.printf("    /DV { [0 %g] 0 setdash LV } def\n", dx/4);
  buf.printf("    /D { moveto\n\
    %g setlinewidth\n\
    0 0.8 0 setrgbcolor [0 %g] 0 setdash\n\
    lineto stroke } def\n", 0.6*dx, 3*dx);
}

static void eps_trailer(file *out) {
  bufprint buf(out);
  buf.printf("grestore\n");
  buf.printf(" 0.25 0 0.5 setrgbcolor\n");
  buf.printf("/Times-Roman findfont 16 scalefont setfont\n");
  buf.printf("newpath 5 5 moveto title show\n");
  buf.printf("grestore\n");
  buf.printf("showpage\n");
  buf.printf("%%%%Trailer\n");
  buf.printf("%%%%EOF\n");
}

// static void eps_dotted(file *out, component m, const double *f, const volume &v,
//                        const geometric_volume &what) {
//   bufprint buf(out);
//   if (!f) return; // Field doesn't exist...
//   for (int i=0;i<v.ntot();i++)
//     if (what.contains(v.loc(m,i)))
//       switch (v.dim) {
//       case Dcyl:
//         {
//           ivec next = v.iloc(m,i)+iveccyl(2,0);
//           if (v.contains(next) && f[i] + f[v.index(m,next)] != 0.0 &&
//               f[i]*f[v.index(m,next)] == 0.0)
//             buf.printf("%g\t%g\tDH\n", v[next].z(), v[next].r() - 0.5/v.a);
//           next = v.iloc(m,i)+iveccyl(0,2);
//           if (v.contains(next) && f[i] + f[v.index(m,next)] != 0.0 &&
//               f[i]*f[v.index(m,next)] == 0.0)
//             buf.printf("%g\t%g\tDV\n", v[next].z() - 0.5/v.a, v[next].r());
//           break;
//         }
//       case D2:
//         {
//           ivec next = v.iloc(m,i)+ivec(2,0);
//           if (v.contains(next) && f[i] + f[v.index(m,next)] != 0.0 &&
//               f[i]*f[v.index(m,next)] == 0.0)
//             buf.printf("%g\t%g\tDH\n", v[next].x() - 0.5/v.a, v[next].y());
//           next = v.iloc(m,i)+ivec(0,2);
//           if (v.contains(next) && f[i] + f[v.index(m,next)] != 0.0 &&
//               f[i]*f[v.index(m,next)] == 0.0)
//             buf.printf("%g\t%g\tDV\n", v[next].x(), v[next].y() - 0.5/v.a);
//           break;
//         }
//       case D1: case D3: break;
//       }
// }

void fields::outline_chunks(file *out) {
  bufprint buf(out);
  if (v.dim == D1) return;
  if (my_rank()) return;
  for (int i=0;i<num_chunks;i++) {
    double xlo, xhi, ylo, yhi;
    switch (v.dim) {
    case Dcyl:
      xlo = chunks[i]->v.boundary_location(Low,Z);
      xhi = chunks[i]->v.boundary_location(High,Z);
      ylo = chunks[i]->v.boundary_location(Low,R);
      yhi = chunks[i]->v.boundary_location(High,R);
    break;
    case D2:
      xlo = chunks[i]->v.boundary_location(Low,X);
      ylo = chunks[i]->v.boundary_location(Low,Y);
      xhi = chunks[i]->v.boundary_location(High,X);
      yhi = chunks[i]->v.boundary_location(High,Y);
    break;
    case D3: // FIXME make this smart about plane direction.
      xlo = chunks[i]->v.boundary_location(Low,X);
      ylo = chunks[i]->v.boundary_location(Low,Y);
      xhi = chunks[i]->v.boundary_location(High,X);
      yhi = chunks[i]->v.boundary_location(High,Y);
    case D1: abort("Error in outline chunks 1D.\n"); break;
    }
    buf.printf("%g\t%g\t%g\t%g\tD\n", xlo, yhi, xhi, yhi);
    buf.printf("%g\t%g\t%g\t%g\tD\n", xhi, yhi, xhi, ylo);
  }
}

static void eps_outline(component m, const double *f,
                        const volume &v, const geometric_volume &what,
                        symmetry S, int symnum, file *out) {
  bufprint buf(out);
  if (!f) return; // Field doesn't exist...
  for (int i=0;i<v.ntot();i++) {
    const vec here = S.transform(v.loc(m,i),symnum);
    if (what.contains(here))
      switch (v.dim) {
      case Dcyl: {