bands.cpp

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

      if (tf[j]<tf[j-1]) {
        double t1 = tf[j], t2 = td[j];
        tf[j] = tf[j-1];
        td[j] = td[j-1];
        tf[j-1] = t1;
        td[j-1] = t2;
        complex<double> temp = ta[j];
        ta[j] = ta[j-1];
        ta[j-1] = temp;
      }
    }
  }
  if (!quiet) master_printf("Looking for clusters...\n");
  int num_found = 0;
  double totwid = 0.001;
  while (freqs_so_far >= fields_considered/2 + 1) {
    int hi, lo;
    get_cluster(tf,freqs_so_far,fields_considered,deltaf,&lo,&hi);
    int mid = lo + (hi-lo)/2;
    if (tf[hi]-tf[lo] < deltaf) {
      if (!quiet) master_printf("Got a cluster from %g to %g (%d freqs)\n",
				tf[lo], tf[hi], 1+hi-lo);
      fad[num_found] = complex<double>(tf[mid],td[mid]);
      if (approx_power) {
        approx_power[num_found] = 0;
        for (int i=lo;i<=hi;i++) {
          if (abs(ta[i])*abs(ta[i]) > approx_power[num_found]) {
            approx_power[num_found] = abs(ta[i])*abs(ta[i]);
          }
        }
      }
      totwid += tf[hi]-tf[lo];
      num_found++;
      if (num_found >= maxbands) num_found--;
    } else {
      if (!quiet)
	master_printf("Rejected a cluster from %g to %g (%d/%d freqs)\n",
		      tf[lo], tf[hi], 1+hi-lo, fields_considered);
      if (verbosity > 1) master_printf("width is %g vs %g\n",
                                       tf[hi] - tf[lo], deltaf);
      lo = get_closest(tf,freqs_so_far);
      hi = lo+1;
      if (verbosity > 1) master_printf("dropping %g and %g\n", tf[hi], tf[lo]);
    }
    freqs_so_far -= 1 + hi - lo;
    for (int i=lo;i<freqs_so_far;i++) {
      tf[i] = tf[i + 1 + hi - lo];
      td[i] = td[i + 1 + hi - lo];
      ta[i] = ta[i + 1 + hi - lo];
    }
  }
  for (int i=0;i<freqs_so_far;i++) {
    if (verbosity > 1) master_printf("Have a leftover freq: %g\n", tf[i]);
  }
  return num_found;
}

complex<double> *fields::clever_cluster_bands(int maxbands, double *approx_power) {
  bands->maxbands = maxbands;
  const int max_harminvs = 120;
  const int max_freqs = max_harminvs*maxbands;
  double *tf = new double[max_freqs];
  double *td = new double[max_freqs];
  complex<double> *ta = new complex<double>[max_freqs];
  const int ntime = bands->ntime;
  if (!ta) abort("Error allocating...\n");
  int num_found = 0;
  complex<double> *fad = new complex<double>[maxbands];
  for (int i=0;i<maxbands;i++) fad[i] = 0.0;
  int freqs_so_far = 0;
  int fields_considered = 0;

  for (int p=0; p<num_bandpts && bands->index[p]!=-1; p++)
    for (int whichf = 0; whichf < 10; whichf++)
      if (v.has_field((component)whichf) && maxbands < max_freqs - freqs_so_far) {
        if (verbosity>1) master_printf("Looking at field %d\n", whichf);
        int freqs_here = bands->get_freqs(bands->f[p][whichf], ntime,
                                          ta+freqs_so_far,
                                          tf+freqs_so_far,
                                          td+freqs_so_far);
        if (freqs_here) {
          fields_considered++;
          freqs_so_far += freqs_here;
        }
        if (freqs_so_far + maxbands > max_freqs) break;
      }
  if (k == 0 && v.dim == Dcyl && m != 0) fields_considered /= 2;
  num_found = cluster_some_bands_cleverly(tf, td, ta, freqs_so_far, fields_considered,
                                          maxbands, fad, approx_power);
  delete[] ta;
  delete[] tf;
  delete[] td;
  // Get rid of bands with too little power in them...
  {
    double maxp = 0.0;
    for (int i=0;i<num_found;i++)
      maxp = max(maxp, approx_power[i]);
    double minp = maxp*bands->fpmin;
    for (int i=0;i<num_found;i++)
      if (approx_power[i] < minp) {
        for (int j=i; j<num_found-1;j++) {
          fad[j] = fad[j+1];
          approx_power[j] = approx_power[j+1];
        }
        num_found--;
        i--;
        fad[num_found] = 0.0;
      }
  }
  // Sorting by frequency again...
  for (int i=0;i<num_found;i++) {
    for (int j=i; j>0;j--) {
      if (real(fad[j])<real(fad[j-1])) {
        complex<double> t1 = fad[j];
        fad[j] = fad[j-1];
        fad[j-1] = t1;
        double temp = approx_power[j];
        approx_power[j] = approx_power[j-1];
        approx_power[j-1] = temp;
      }
    }
  }
  return fad;
}

int bandsdata::get_freqs(complex<double> *data, int n, complex<double> *amps,
                         double *freq_re, double *freq_im) {
  
  const double total_time = n*dt;
  const double qminhere = 1.0/(1.0/qmin + 0.25/(fmin*total_time));
  return do_harminv(data, n, dt, fmin, fmax, maxbands,
		    amps,  freq_re, freq_im, NULL, 1.1, qminhere);
}

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,
	       double spectral_density, double Q_thresh, double rel_err_thresh, double err_thresh, double rel_amp_thresh, double amp_thresh) {
#ifndef HAVE_HARMINV
  abort("compiled without Harminv library, required for do_harminv");
  return 0;
#else
  int numfreqs = int(fabs(fmax-fmin)*dt*n*spectral_density); // c.f. harminv
  if (numfreqs > 150) numfreqs = 150; // prevent matrices from getting too big
  if (numfreqs < 2) numfreqs = 2;

  if (maxbands > numfreqs)
       numfreqs = maxbands;

  // check for all zeros in input
  // data is a size n array.
  {
    int i;
    for (i = 0; i < n && data[i] == 0.0; i++)
      ;
    if (i == n)
      return 0;
  }

#if 0
  // debugging: save data file and arguments for standalone harminv program
  {
    FILE *f = fopen("harminv.dat", "w");
    fprintf(f, "# -f %d -t %g %g-%g -Q %e -e %e -E %e -a %e -A %e -F\n",
	    numfreqs, dt, fmin, fmax,
	    Q_thresh, rel_err_thresh, err_thresh, rel_amp_thresh, amp_thresh);
    for (int i = 0; i < n; ++i)
      fprintf(f, "%g%+gi\n", real(data[i]), imag(data[i]));
    fclose(f);
  }
#endif

  harminv_data hd = harminv_data_create(n, data, fmin*dt, fmax*dt, numfreqs);
  harminv_solve(hd);

  int nf = harminv_get_num_freqs(hd);
  if (nf == 0) return 0;
  int *fsort = new int[nf]; // indices of frequencies, sorted as needed
  
  for (int i = 0; i < nf; ++i)
    fsort[i] = i;
  for (int i = 0; i < nf; ++i) // sort in increasing order of error
    for (int j = i + 1; j < nf; ++j) 
      if (harminv_get_freq_error(hd, fsort[i]) >
	  harminv_get_freq_error(hd, fsort[j])) {
	int k = fsort[i];
	fsort[i] = fsort[j];
	fsort[j] = k;
      }
  
  double min_err = harminv_get_freq_error(hd, fsort[0]);
  double max_amp = abs(harminv_get_amplitude(hd, 0));
  for (int i = 1; i < nf; ++i) {
    double amp = abs(harminv_get_amplitude(hd, i));
    if (max_amp < amp)
      max_amp = amp;
  }
  { // eliminate modes that fall outside the various thresholds:
    int j = 0;
    for (int i = 0; i < nf; ++i) {
      double f = abs(harminv_get_freq(hd, fsort[i]) / dt);
      double err = harminv_get_freq_error(hd, fsort[i]);
      double amp = abs(harminv_get_amplitude(hd, fsort[i]));
      if (f >= fmin && f <= fmax
	  && abs(harminv_get_Q(hd, fsort[i])) > Q_thresh
	  && err < err_thresh
	  && err < rel_err_thresh * min_err
	  && amp > amp_thresh
	  && amp > rel_amp_thresh * max_amp) {
	fsort[j++] = fsort[i];
      }
    }
    nf = j;
  }
  { // eliminate positive/negative frequency pairs
    // set indices to -1 for frequencies to be eliminated
    for (int i = 0; i < nf; ++i) 
      if (fsort[i] != -1) { // i hasn't been eliminated yet
	double f = harminv_get_freq(hd, fsort[i]);
	if (f < 0.0) {
	  double kdiff = -2 * f;
	  int kpos = i;
	  for (int k = 0; k < nf; ++k) // search for closest positive freq.
	    if (fsort[k] != -1) { // k hasn't been eliminated yet
	      double fdiff = abs(harminv_get_freq(hd, fsort[k]) + f);
	      if (fdiff < kdiff) {
		kpos = k;
		kdiff = fdiff;
	      }
	    }
	  if (kpos != i && kdiff < 2.0 / n) { // consider them the same
	    // pick the one with the smaller error
	    if (harminv_get_freq_error(hd, fsort[i]) <
		harminv_get_freq_error(hd, fsort[kpos]))
	      fsort[kpos] = -1;
	    else
	      fsort[i] = -1;
	  }
	}
      }
    int j = 0;
    for (int i = 0; i < nf; ++i) // remove the eliminated indices
      if (fsort[i] != -1)
	fsort[j++] = fsort[i];
    nf = j;
  }
  
  if (nf > maxbands)
    nf = maxbands;
  
  // sort again, this time in increasing order of freq:
  for (int i = 0; i < nf; ++i) // simple O(nf^2) sort
    for (int j = i + 1; j < nf; ++j) 
      if (abs(harminv_get_freq(hd, fsort[i])) >
	  abs(harminv_get_freq(hd, fsort[j]))) {
	int k = fsort[i];
	fsort[i] = fsort[j];
	fsort[j] = k;
      }
  
  for (int i = 0; i < nf; ++i) {
    complex<double> freq = harminv_get_omega(hd, fsort[i]) / (2*pi*dt);
    freq_re[i] = abs(real(freq));
    freq_im[i] = imag(freq);
    amps[i] = harminv_get_amplitude(hd, fsort[i]);
    if (errors)
      errors[i] = harminv_get_freq_error(hd, fsort[i]);
  }

  delete[] fsort;
  harminv_data_destroy(hd);
  return nf;
#endif
}

} // namespace meep