"Fossies" - the Fresh Open Source Software Archive

Member "meep-1.10.0/src/initialize.cpp" (5 Jun 2019, 4138 Bytes) of package /linux/privat/meep-1.10.0.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "initialize.cpp" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 1.9.0_vs_1.10.0.

    1 /* Copyright (C) 2005-2019 Massachusetts Institute of Technology
    2 %
    3 %  This program is free software; you can redistribute it and/or modify
    4 %  it under the terms of the GNU General Public License as published by
    5 %  the Free Software Foundation; either version 2, or (at your option)
    6 %  any later version.
    7 %
    8 %  This program is distributed in the hope that it will be useful,
    9 %  but WITHOUT ANY WARRANTY; without even the implied warranty of
   10 %  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   11 %  GNU General Public License for more details.
   12 %
   13 %  You should have received a copy of the GNU General Public License
   14 %  along with this program; if not, write to the Free Software Foundation,
   15 %  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
   16 */
   17 
   18 #include <stdio.h>
   19 #include <stdlib.h>
   20 #include <math.h>
   21 #include <complex>
   22 
   23 #include "meep.hpp"
   24 #include "meep_internals.hpp"
   25 #include "config.h"
   26 
   27 // Cylindrical coordinates:
   28 
   29 #ifdef HAVE_LIBGSL
   30 #include <gsl/gsl_sf_bessel.h>
   31 #endif
   32 
   33 using namespace std;
   34 
   35 namespace meep {
   36 
   37 #define J BesselJ
   38 double J(int m, double kr) {
   39 #if defined(HAVE_JN)
   40   return jn(m, kr); // POSIX/BSD jn function
   41 #elif defined(HAVE_LIBGSL)
   42   return gsl_sf_bessel_Jn(m, kr);
   43 #else
   44   abort("not compiled with GSL, required for Bessel functions");
   45   return 0;
   46 #endif
   47 }
   48 static double Jprime(int m, double kr) {
   49   if (m)
   50     return 0.5 * (J(m - 1, kr) - J(m + 1, kr));
   51   else
   52     return -J(1, kr);
   53 }
   54 static double Jroot(int m, int n) {
   55 #ifdef HAVE_LIBGSL
   56   return gsl_sf_bessel_zero_Jnu(m, n + 1);
   57 #else
   58   (void)m;
   59   (void)n;
   60   abort("not compiled with GSL, required for Bessel functions");
   61   return 0;
   62 #endif
   63 }
   64 static double Jmax(int m, int n) {
   65   double rlow, rhigh = Jroot(m, n), rtry;
   66   if (n == 0)
   67     rlow = 0;
   68   else
   69     rlow = Jroot(m, n - 1);
   70   double jplow = Jprime(m, rlow), jptry;
   71   do {
   72     rtry = rlow + (rhigh - rlow) * 0.5;
   73     jptry = Jprime(m, rtry);
   74     if (jplow * jptry < 0)
   75       rhigh = rtry;
   76     else
   77       rlow = rtry;
   78   } while (rhigh - rlow > rhigh * 1e-15);
   79   return rtry;
   80 }
   81 
   82 static double ktrans, kax;
   83 static int m_for_J;
   84 static complex<double> JJ(const vec &pt) {
   85   return polar(J(m_for_J, ktrans * pt.r()), kax * pt.r());
   86 }
   87 static complex<double> JP(const vec &pt) {
   88   return polar(Jprime(m_for_J, ktrans * pt.r()), kax * pt.r());
   89 }
   90 
   91 void fields::initialize_with_nth_te(int np0) {
   92   require_component(Hz);
   93   for (int i = 0; i < num_chunks; i++)
   94     chunks[i]->initialize_with_nth_te(np0, real(k[Z]));
   95 }
   96 
   97 void fields_chunk::initialize_with_nth_te(int np0, double kz) {
   98   const int im = int(m);
   99   const int n = (im == 0) ? np0 - 0 : np0 - 1;
  100   const double rmax = Jmax(im, n);
  101   ktrans = rmax * a / gv.nr();
  102   kax = kz * 2 * pi / a;
  103   m_for_J = im;
  104   initialize_field(Hz, JJ);
  105 }
  106 
  107 void fields::initialize_with_nth_tm(int np0) {
  108   require_component(Ez);
  109   require_component(Hp);
  110   for (int i = 0; i < num_chunks; i++)
  111     chunks[i]->initialize_with_nth_tm(np0, real(k[Z]));
  112 }
  113 
  114 void fields_chunk::initialize_with_nth_tm(int np1, double kz) {
  115   const int im = int(m);
  116   const int n = np1 - 1;
  117   const double rroot = Jroot(im, n);
  118   ktrans = rroot * a / gv.nr();
  119   kax = kz * 2 * pi / a;
  120   m_for_J = im;
  121   initialize_field(Ez, JJ);
  122   initialize_field(Hp, JP);
  123 }
  124 
  125 void fields::initialize_with_n_te(int ntot) {
  126   for (int n = 0; n < ntot; n++)
  127     initialize_with_nth_te(n + 1);
  128 }
  129 
  130 void fields::initialize_with_n_tm(int ntot) {
  131   for (int n = 0; n < ntot; n++)
  132     initialize_with_nth_tm(n + 1);
  133 }
  134 
  135 void fields::initialize_field(component c, complex<double> func(const vec &)) {
  136   require_component(c);
  137   for (int i = 0; i < num_chunks; i++)
  138     chunks[i]->initialize_field(c, func);
  139   step_boundaries(type(c));
  140   if (is_D(c)) {
  141     update_eh(E_stuff);
  142     step_boundaries(E_stuff);
  143   }
  144   if (is_B(c)) {
  145     update_eh(H_stuff);
  146     step_boundaries(H_stuff);
  147   }
  148 }
  149 
  150 void fields_chunk::initialize_field(component c, complex<double> func(const vec &)) {
  151   if (f[c][0]) {
  152     LOOP_OVER_VOL(gv, c, i) {
  153       IVEC_LOOP_LOC(gv, here);
  154       complex<double> val = func(here);
  155       f[c][0][i] += real(val);
  156       if (!is_real) f[c][1][i] += imag(val);
  157     }
  158   }
  159 }
  160 
  161 } // namespace meep