bands.cpp

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

/* 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 "meep.hpp"
#include "meep_internals.hpp"

#include "config.h"
#ifdef HAVE_HARMINV
#  include <harminv.h>
#endif

namespace meep {

#define BAND(b,r,t) ((b)[(r)+(t)*nr])

bandsdata::bandsdata() {
  verbosity = 0;
  maxbands = -1;
  tstart = 0;
  for (int i=0;i<num_bandpts;i++) index[i] = -1;
  tend = -1;
  for (int i=0;i<num_bandpts;i++) for (int c=0;c<10;c++) f[i][c] = NULL;
  P = NULL;
}

bandsdata::~bandsdata() {
  for (int i=0;i<num_bandpts;i++) for (int c=0;c<10;c++) delete[] f[i][c];
  delete[] P;
}

double fields::last_source_time() {
  double last_time = 0;
    if (sources != NULL)
      last_time = max(last_time, sources->last_time_max());
  for (int i=0;i<num_chunks;i++)
    if (chunks[i]->is_mine())
      last_time = max(last_time,chunks[i]->last_source_time());
  return max_to_all(last_time);
}

double fields_chunk::last_source_time() {
  return 0;
}

void fields::prepare_for_bands(const vec &p, double endtime, double fmax,
                               double qmin, double frac_pow_min) {
  int last_source = (int)(last_source_time()/dt+0.5);
  last_source = max(last_source, t + phasein_time);
  if (!bands) bands = new bandsdata;
  bands->tstart = last_source+1;
  if (bands->tstart < t) bands->tstart = t;
  bands->tend = t + (int)(endtime/dt) - 1;

  {
    int ind[8];
    double w[8];
    int indind = 0;
    while (bands->index[indind] != -1 && indind < num_bandpts) indind++;
    for (int h=0;h<num_chunks;h++)
      if (chunks[h]->v.contains(p)) {
        chunks[h]->v.interpolate(chunks[h]->v.eps_component(), p, ind, w);
        for (int i=0;i<8&&w[i]&&indind<num_bandpts;i++) {
          bands->chunk[indind] = chunks[h];
          bands->index[indind++] = ind[i];
        }
        break;
      }
  }
  bands->fpmin = frac_pow_min;

  // Set fmin properly...
  const double epsmax = max_eps();

  double cutoff_freq = 0.0;
  if (v.dim == Dcyl) {
    cutoff_freq = 1.84*a*dt/(2*pi)/v.nr()/sqrt(epsmax);
    if (m == 0) cutoff_freq *= 0.5;
  }
  bands->fmin = sqrt(cutoff_freq*cutoff_freq + abs(k[Z])*abs(k[Z])*(a*dt)*(a*dt)/epsmax); // FIXME
  bands->fmin = cutoff_freq/dt;
  bands->qmin = qmin;
  // Set fmax and determine how many timesteps to skip over...
  bands->fmax = fmax;
  {
    // for when there are too many data points...
    double decayconst = bands->fmax*dt/qmin*8.0;
    double smalltime = 1./(decayconst + bands->fmax*dt);
    bands->scale_factor = (int)(0.06*smalltime);
    if (bands->scale_factor < 1) bands->scale_factor = 1;
    if (verbosity) master_printf("scale_factor is %d (%g,%g)\n",
                                 bands->scale_factor, bands->fmax*dt,
                                 decayconst);
  }

  if (bands->tend <= bands->tstart) {
    printf("Oi, we don't have any time to take a fourier transform!\n");
    abort("FT start is %d and end is %d\n", bands->tstart, bands->tend);
  }
  bands->ntime = (1+(bands->tend-bands->tstart)/bands->scale_factor);
  bands->dt = dt * bands->scale_factor;
  for (int c=0;c<10;c++)
    for (int i=0;i<num_bandpts;i++)
      if (v.has_field((component)c)) {
        delete[] bands->f[i][c];
        bands->f[i][c] = new complex<double>[bands->ntime];
        if (bands->f[i][c] == NULL)
          abort("Unable to allocate bandstructure array!\n");
        for (int j=0;j<bands->ntime;j++) bands->f[i][c][j] = 0.0;
      }
  bands->P = new complex<double>[bands->ntime];
  for (int i=0;i<bands->ntime;i++) bands->P[i] = 0.0;
  bands->verbosity = verbosity;
  for (int h=0;h<num_chunks;h++)
    chunks[h]->bands = bands;
}

void fields::record_bands() {
  if (!bands) return;
  if (t > bands->tend || t < bands->tstart) return;
  if (t % bands->scale_factor != 0) return;
  for (int i=0;i<num_chunks;i++)
    chunks[i]->record_bands(t);
}

void fields_chunk::record_bands(int tcount) {
  int thet = (tcount-bands->tstart)/bands->scale_factor;
  if (thet >= bands->ntime) return;
  for (int p=0; p<num_bandpts && bands->index[p]!=-1; p++)
    if (this == bands->chunk[p])
      for (int c=0;c<10;c++)
        if (v.has_field((component)c)) {
          complex<double> tmp;
          if (f[c][0] && f[c][1]) tmp = complex<double>(f[c][0][bands->index[p]],
                                                        f[c][1][bands->index[p]]);
          bands->f[p][c][thet] = broadcast(n_proc(), tmp);
        }
}

#define HARMOUT(o,n,f) ((o)[(n)+(f)*maxbands])

complex<double> fields::get_band(int nn, int maxbands) {
  //complex<double> *fad = get_the_bands(maxbands, approx_power);
  complex<double> *fad = clever_cluster_bands(maxbands);
  complex<double> thef = fad[nn-1];
  delete[] fad;
  return thef;
}

void fields::grace_bands(grace *g, int maxbands) {
  double *approx_power = new double[maxbands];
  //complex<double> *fad = get_the_bands(maxbands, approx_power);
  complex<double> *fad = clever_cluster_bands(maxbands, approx_power);

  int num_found = 0;
  for (int i=0;i<maxbands;i++) if (fad[i] != 0.0) num_found = i+1;
  for (int i = 0; i < num_found; ++i) {
    g->output_out_of_order(i, abs(k[Z]), fabs(real(fad[i])), fabs(imag(fad[i])),
                           approx_power[i]); // FIXME
  }
  delete[] fad;
  delete[] approx_power;
}

void fields::output_bands(FILE *o, const char *name, int maxbands) {
  out_bands(o, name, maxbands);
}

void fields::out_bands(FILE *o, const char *name, int maxbands) {
  double *approx_power = new double[maxbands];
  //complex<double> *fad = get_the_bands(maxbands, approx_power);
  complex<double> *fad = clever_cluster_bands(maxbands, approx_power);

  complex<double> *eigen = new complex<double>[maxbands*6];
  if (!eigen) abort("Error allocating...\n");

  for (int whichf = 0; whichf < 6; whichf++) {
    for (int n=0;n<maxbands;n++) {
      HARMOUT(eigen,n,whichf) = 0;
    }
  }
  int num_found = 0;
  for (int i=0;i<maxbands;i++) if (fad[i] != 0.0) num_found = i+1;

  if (am_master()) {
    for (int i = 0; i < num_found; ++i) {
      // k k k m index freq decay Q approx_power
      fprintf(o, "%s\t%g\t%g\t%g\t%g\t%d\t%g \t%g \t%g \t%g\n", 
	      name, 
	      real(k[0]), real(k[1]), real(k[2]),   
	      m, i, fabs(real(fad[i])), imag(fad[i]),
	      -fabs(real(fad[i])) / (2 * imag(fad[i])),
	      approx_power[i]);
    }
    fflush(o);
  }
  delete[] approx_power;
  delete[] fad;
}

static inline int get_closest(double f[], int fmax) {
  double deltamin = 1e300;
  for (int i=0;i<fmax-1;i++) {
    if (f[i+1]-f[i] < deltamin) deltamin = f[i+1]-f[i];
  }
  for (int i=0;i<fmax-1;i++) {
    if (f[i+1]-f[i] == deltamin) return i;
  }
  return 0;
}

static inline int go_higher(double f[], int fmax, int lo, int hi) {
  if (lo == 0) return 1;
  if (hi == fmax-1) return 0;
  if (f[lo]-f[lo-1] > f[hi+1]-f[hi]) return 1;
  else return 0;
}

static inline int am_done(double f[], int fmax, int lo, int hi) {
  double wid = f[hi]-f[lo] + 0.001;
  int lodone = lo == 0 || f[lo]-f[lo-1] > wid;
  int hidone = hi == fmax-1 || f[hi+1]-f[hi] > wid;
  return lodone && hidone;
}

static void get_cluster(double f[], int fmax, int maxsize, double maxwid,
                        int *out_lo, int *out_hi) {
  int lo = get_closest(f,fmax);
  int hi = lo+1;
  int minsize = maxsize/2+1;
  if (minsize < 3) minsize = 3;
  for (int i=0;i<minsize-2;i++) {
    if (go_higher(f,fmax,lo,hi)) hi++;
    else lo--;
  }
  while (hi + 1 - lo < maxsize) {
    if (am_done(f,fmax,lo,hi)) break;
    if (go_higher(f,fmax,lo,hi)) {
      if (f[hi+1]-f[lo] > maxwid) break;
      hi++;
    } else {
      if (f[hi]-f[lo-1] > maxwid) break;
      lo--;
    }
  }
  *out_lo = lo;
  *out_hi = hi;
}

int fields::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) {
  const double total_time = (bands->tend-bands->tstart)*dt;
  const double deltaf = 1.0/total_time;
  int freqs_so_far = num_freqs;
  if (!quiet) master_printf("About to sort by frequency... (%d frequencies)\n", freqs_so_far);
  // Sort by frequency...
  for (int i = 1; i < freqs_so_far; i++) {
    for (int j=i; j>0;j--) {