update_e_from_d_update.hpp

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

/* This file was automatically generated --- DO NOT EDIT */

if (have_d_minus_p) {
  if (f[Er][0]) {
    // Am in cylindrical coordinates.
    FOR_E_AND_D(ec,dc) if (f[ec][0]) {
      const int yee_idx = v.yee_index(ec);
      const int d_ec = component_direction(ec);
      const int s_ec = stride_any_direction[d_ec];
      const direction d_1 = (direction)(((d_ec-2)+1)%3+2);
      const component ec_1 = direction_component(ec,d_1);
      const component dc_1 = direction_component(dc,d_1);
      const int s_1 = stride_any_direction[d_1];
      const direction d_2 = (direction)(((d_ec-2)+2)%3+2);
      const component ec_2 = direction_component(ec,d_2);
      const component dc_2 = direction_component(dc,d_2);
      const int s_2 = stride_any_direction[d_2];
      if (s->inveps[ec][d_1]) {
        if (s->inveps[ec][d_2]) {
          if (s->kerr[ec]) {
            if (is_real) {
              if (num_any_direction[Z]==1) {
                const int sR = stride_any_direction[R];
                const int sZ = 1;
                for (int iR=0; iR<num_any_direction[R]; iR++) {
                  const int i = yee_idx + iR*sR;
                  const double dsqr_here = d_minus_p[ec][0][i]*d_minus_p[ec][0][i]
                     + (0.25*(d_minus_p[ec_1][0][i]*d_minus_p[ec_1][0][i] + ((d_minus_p[ec_1][0][i-s_1])*(d_minus_p[ec_1][0][i-s_1]))
                        + ((d_minus_p[ec_1][0][i+s_ec])*(d_minus_p[ec_1][0][i+s_ec]))
                           + (d_minus_p[ec_1][0][i-s_1+s_ec])*(d_minus_p[ec_1][0][i-s_1+s_ec])))
                              + 0.25*(d_minus_p[ec_2][0][i]*d_minus_p[ec_2][0][i] + ((d_minus_p[ec_2][0][i-s_2])*(d_minus_p[ec_2][0][i-s_2]))
                                 + ((d_minus_p[ec_2][0][i+s_ec])*(d_minus_p[ec_2][0][i+s_ec]))
                                    + (d_minus_p[ec_2][0][i-s_2+s_ec])*(d_minus_p[ec_2][0][i-s_2+s_ec]));
                  const double frac_change = calc_nonlinear_inveps(dsqr_here,
                     s->inveps[ec][d_ec][i], s->kerr[ec][i]);
                  f[ec][0][i] = frac_change*(((s->inveps[ec][d_ec][i])*d_minus_p[ec][0][i])
                     + (0.25*((s->inveps[ec][d_1][i])*(d_minus_p[ec_1][0][i] + (d_minus_p[ec_1][0][i-s_1]) + (d_minus_p[ec_1][0][i+s_ec]) + d_minus_p[ec_1][0][i-s_1+s_ec]))) + 0.25*((s->inveps[ec][d_2][i])*(d_minus_p[ec_2][0][i] + (d_minus_p[ec_2][0][i-s_2]) + (d_minus_p[ec_2][0][i+s_ec]) + d_minus_p[ec_2][0][i-s_2+s_ec])));
                }
              } else { // not num_any_direction[Z]==1
                if (num_any_direction[R]==1) {
                  const int sR = stride_any_direction[R];
                  const int sZ = 1;
                  for (int iZ=0; iZ<num_any_direction[Z]; iZ++) {
                    const int i = yee_idx + iZ;
                    const double dsqr_here = d_minus_p[ec][0][i]*d_minus_p[ec][0][i]
                       + (0.25*(d_minus_p[ec_1][0][i]*d_minus_p[ec_1][0][i] + ((d_minus_p[ec_1][0][i-s_1])*(d_minus_p[ec_1][0][i-s_1]))
                          + ((d_minus_p[ec_1][0][i+s_ec])*(d_minus_p[ec_1][0][i+s_ec]))
                             + (d_minus_p[ec_1][0][i-s_1+s_ec])*(d_minus_p[ec_1][0][i-s_1+s_ec])))
                                + 0.25*(d_minus_p[ec_2][0][i]*d_minus_p[ec_2][0][i] + ((d_minus_p[ec_2][0][i-s_2])*(d_minus_p[ec_2][0][i-s_2]))
                                   + ((d_minus_p[ec_2][0][i+s_ec])*(d_minus_p[ec_2][0][i+s_ec]))
                                      + (d_minus_p[ec_2][0][i-s_2+s_ec])*(d_minus_p[ec_2][0][i-s_2+s_ec]));
                    const double frac_change = calc_nonlinear_inveps(dsqr_here,
                       s->inveps[ec][d_ec][i], s->kerr[ec][i]);
                    f[ec][0][i] = frac_change*(((s->inveps[ec][d_ec][i])*d_minus_p[ec][0][i])
                       + (0.25*((s->inveps[ec][d_1][i])*(d_minus_p[ec_1][0][i] + (d_minus_p[ec_1][0][i-s_1]) + (d_minus_p[ec_1][0][i+s_ec]) + d_minus_p[ec_1][0][i-s_1+s_ec]))) + 0.25*((s->inveps[ec][d_2][i])*(d_minus_p[ec_2][0][i] + (d_minus_p[ec_2][0][i-s_2]) + (d_minus_p[ec_2][0][i+s_ec]) + d_minus_p[ec_2][0][i-s_2+s_ec])));
                  }
                } else { // not num_any_direction[R]==1
                  const int sR = stride_any_direction[R];
                  const int sZ = 1;
                  for (int iR=0; iR<num_any_direction[R]; iR++) {
                    for (int iZ=0; iZ<num_any_direction[Z]; iZ++) {
                      const int i = yee_idx + iZ + iR*sR;
                      const double dsqr_here = d_minus_p[ec][0][i]*d_minus_p[ec][0][i]
                         + (0.25*(d_minus_p[ec_1][0][i]*d_minus_p[ec_1][0][i] + ((d_minus_p[ec_1][0][i-s_1])*(d_minus_p[ec_1][0][i-s_1]))
                            + ((d_minus_p[ec_1][0][i+s_ec])*(d_minus_p[ec_1][0][i+s_ec]))
                               + (d_minus_p[ec_1][0][i-s_1+s_ec])*(d_minus_p[ec_1][0][i-s_1+s_ec])))
                                  + 0.25*(d_minus_p[ec_2][0][i]*d_minus_p[ec_2][0][i] + ((d_minus_p[ec_2][0][i-s_2])*(d_minus_p[ec_2][0][i-s_2]))
                                     + ((d_minus_p[ec_2][0][i+s_ec])*(d_minus_p[ec_2][0][i+s_ec]))
                                        + (d_minus_p[ec_2][0][i-s_2+s_ec])*(d_minus_p[ec_2][0][i-s_2+s_ec]));
                      const double frac_change = calc_nonlinear_inveps(dsqr_here,
                         s->inveps[ec][d_ec][i], s->kerr[ec][i]);
                      f[ec][0][i] = frac_change*(((s->inveps[ec][d_ec][i])*d_minus_p[ec][0][i])
                         + (0.25*((s->inveps[ec][d_1][i])*(d_minus_p[ec_1][0][i] + (d_minus_p[ec_1][0][i-s_1]) + (d_minus_p[ec_1][0][i+s_ec]) + d_minus_p[ec_1][0][i-s_1+s_ec]))) + 0.25*((s->inveps[ec][d_2][i])*(d_minus_p[ec_2][0][i] + (d_minus_p[ec_2][0][i-s_2]) + (d_minus_p[ec_2][0][i+s_ec]) + d_minus_p[ec_2][0][i-s_2+s_ec])));
                    }
                  }
                }
              }
            } else { // not is_real
              if (num_any_direction[Z]==1) {
                const int sR = stride_any_direction[R];
                const int sZ = 1;
                for (int iR=0; iR<num_any_direction[R]; iR++) {
                  const int i = yee_idx + iR*sR;
                  const double dsqr_here = (d_minus_p[ec][1][i]*d_minus_p[ec][1][i]
                     + d_minus_p[ec][0][i]*d_minus_p[ec][0][i]) + (0.25*((d_minus_p[ec_1][1][i]*d_minus_p[ec_1][1][i]
                        + d_minus_p[ec_1][0][i]*d_minus_p[ec_1][0][i]) + (((d_minus_p[ec_1][1][i-s_1])*(d_minus_p[ec_1][1][i-s_1]))
                           + (d_minus_p[ec_1][0][i-s_1])*(d_minus_p[ec_1][0][i-s_1]))
                              + (((d_minus_p[ec_1][1][i+s_ec])*(d_minus_p[ec_1][1][i+s_ec]))
                                 + (d_minus_p[ec_1][0][i+s_ec])*(d_minus_p[ec_1][0][i+s_ec]))
                                    + ((d_minus_p[ec_1][1][i-s_1+s_ec])*(d_minus_p[ec_1][1][i-s_1+s_ec]))
                                       + (d_minus_p[ec_1][0][i-s_1+s_ec])*(d_minus_p[ec_1][0][i-s_1+s_ec])))
                                          + 0.25*((d_minus_p[ec_2][1][i]*d_minus_p[ec_2][1][i] + d_minus_p[ec_2][0][i]*d_minus_p[ec_2][0][i])
                                             + (((d_minus_p[ec_2][1][i-s_2])*(d_minus_p[ec_2][1][i-s_2]))
                                                + (d_minus_p[ec_2][0][i-s_2])*(d_minus_p[ec_2][0][i-s_2]))
                                                   + (((d_minus_p[ec_2][1][i+s_ec])*(d_minus_p[ec_2][1][i+s_ec]))
                                                      + (d_minus_p[ec_2][0][i+s_ec])*(d_minus_p[ec_2][0][i+s_ec]))
                                                         + ((d_minus_p[ec_2][1][i-s_2+s_ec])*(d_minus_p[ec_2][1][i-s_2+s_ec]))
                                                            + (d_minus_p[ec_2][0][i-s_2+s_ec])*(d_minus_p[ec_2][0][i-s_2+s_ec]));
                  const double frac_change = calc_nonlinear_inveps(dsqr_here,
                     s->inveps[ec][d_ec][i], s->kerr[ec][i]);
                  f[ec][0][i] = frac_change*(((s->inveps[ec][d_ec][i])*d_minus_p[ec][0][i])
                     + (0.25*((s->inveps[ec][d_1][i])*(d_minus_p[ec_1][0][i] + (d_minus_p[ec_1][0][i-s_1]) + (d_minus_p[ec_1][0][i+s_ec]) + d_minus_p[ec_1][0][i-s_1+s_ec]))) + 0.25*((s->inveps[ec][d_2][i])*(d_minus_p[ec_2][0][i] + (d_minus_p[ec_2][0][i-s_2]) + (d_minus_p[ec_2][0][i+s_ec]) + d_minus_p[ec_2][0][i-s_2+s_ec])));
                  f[ec][1][i] = frac_change*(((s->inveps[ec][d_ec][i])*d_minus_p[ec][1][i])
                     + (0.25*((s->inveps[ec][d_1][i])*(d_minus_p[ec_1][1][i] + (d_minus_p[ec_1][1][i-s_1]) + (d_minus_p[ec_1][1][i+s_ec]) + d_minus_p[ec_1][1][i-s_1+s_ec]))) + 0.25*((s->inveps[ec][d_2][i])*(d_minus_p[ec_2][1][i] + (d_minus_p[ec_2][1][i-s_2]) + (d_minus_p[ec_2][1][i+s_ec]) + d_minus_p[ec_2][1][i-s_2+s_ec])));
                }
              } else { // not num_any_direction[Z]==1
                if (num_any_direction[R]==1) {
                  const int sR = stride_any_direction[R];
                  const int sZ = 1;
                  for (int iZ=0; iZ<num_any_direction[Z]; iZ++) {
                    const int i = yee_idx + iZ;
                    const double dsqr_here = (d_minus_p[ec][1][i]*d_minus_p[ec][1][i]
                       + d_minus_p[ec][0][i]*d_minus_p[ec][0][i]) + (0.25*((d_minus_p[ec_1][1][i]*d_minus_p[ec_1][1][i]
                          + d_minus_p[ec_1][0][i]*d_minus_p[ec_1][0][i]) + (((d_minus_p[ec_1][1][i-s_1])*(d_minus_p[ec_1][1][i-s_1]))
                             + (d_minus_p[ec_1][0][i-s_1])*(d_minus_p[ec_1][0][i-s_1]))
                                + (((d_minus_p[ec_1][1][i+s_ec])*(d_minus_p[ec_1][1][i+s_ec]))
                                   + (d_minus_p[ec_1][0][i+s_ec])*(d_minus_p[ec_1][0][i+s_ec]))
                                      + ((d_minus_p[ec_1][1][i-s_1+s_ec])*(d_minus_p[ec_1][1][i-s_1+s_ec]))
                                         + (d_minus_p[ec_1][0][i-s_1+s_ec])*(d_minus_p[ec_1][0][i-s_1+s_ec])))
                                            + 0.25*((d_minus_p[ec_2][1][i]*d_minus_p[ec_2][1][i] + d_minus_p[ec_2][0][i]*d_minus_p[ec_2][0][i])
                                               + (((d_minus_p[ec_2][1][i-s_2])*(d_minus_p[ec_2][1][i-s_2]))
                                                  + (d_minus_p[ec_2][0][i-s_2])*(d_minus_p[ec_2][0][i-s_2]))
                                                     + (((d_minus_p[ec_2][1][i+s_ec])*(d_minus_p[ec_2][1][i+s_ec]))
                                                        + (d_minus_p[ec_2][0][i+s_ec])*(d_minus_p[ec_2][0][i+s_ec]))
                                                           + ((d_minus_p[ec_2][1][i-s_2+s_ec])*(d_minus_p[ec_2][1][i-s_2+s_ec]))
                                                              + (d_minus_p[ec_2][0][i-s_2+s_ec])*(d_minus_p[ec_2][0][i-s_2+s_ec]));
                    const double frac_change = calc_nonlinear_inveps(dsqr_here,
                       s->inveps[ec][d_ec][i], s->kerr[ec][i]);
                    f[ec][0][i] = frac_change*(((s->inveps[ec][d_ec][i])*d_minus_p[ec][0][i])
                       + (0.25*((s->inveps[ec][d_1][i])*(d_minus_p[ec_1][0][i] + (d_minus_p[ec_1][0][i-s_1]) + (d_minus_p[ec_1][0][i+s_ec]) + d_minus_p[ec_1][0][i-s_1+s_ec]))) + 0.25*((s->inveps[ec][d_2][i])*(d_minus_p[ec_2][0][i] + (d_minus_p[ec_2][0][i-s_2]) + (d_minus_p[ec_2][0][i+s_ec]) + d_minus_p[ec_2][0][i-s_2+s_ec])));
                    f[ec][1][i] = frac_change*(((s->inveps[ec][d_ec][i])*d_minus_p[ec][1][i])
                       + (0.25*((s->inveps[ec][d_1][i])*(d_minus_p[ec_1][1][i] + (d_minus_p[ec_1][1][i-s_1]) + (d_minus_p[ec_1][1][i+s_ec]) + d_minus_p[ec_1][1][i-s_1+s_ec]))) + 0.25*((s->inveps[ec][d_2][i])*(d_minus_p[ec_2][1][i] + (d_minus_p[ec_2][1][i-s_2]) + (d_minus_p[ec_2][1][i+s_ec]) + d_minus_p[ec_2][1][i-s_2+s_ec])));
                  }
                } else { // not num_any_direction[R]==1
                  const int sR = stride_any_direction[R];
                  const int sZ = 1;
                  for (int iR=0; iR<num_any_direction[R]; iR++) {
                    for (int iZ=0; iZ<num_any_direction[Z]; iZ++) {
                      const int i = yee_idx + iZ + iR*sR;
                      const double dsqr_here = (d_minus_p[ec][1][i]*d_minus_p[ec][1][i]
                         + d_minus_p[ec][0][i]*d_minus_p[ec][0][i]) + (0.25*((d_minus_p[ec_1][1][i]*d_minus_p[ec_1][1][i]
                            + d_minus_p[ec_1][0][i]*d_minus_p[ec_1][0][i]) + (((d_minus_p[ec_1][1][i-s_1])*(d_minus_p[ec_1][1][i-s_1]))
                               + (d_minus_p[ec_1][0][i-s_1])*(d_minus_p[ec_1][0][i-s_1]))
                                  + (((d_minus_p[ec_1][1][i+s_ec])*(d_minus_p[ec_1][1][i+s_ec]))
                                     + (d_minus_p[ec_1][0][i+s_ec])*(d_minus_p[ec_1][0][i+s_ec]))
                                        + ((d_minus_p[ec_1][1][i-s_1+s_ec])*(d_minus_p[ec_1][1][i-s_1+s_ec]))
                                           + (d_minus_p[ec_1][0][i-s_1+s_ec])*(d_minus_p[ec_1][0][i-s_1+s_ec])))
                                              + 0.25*((d_minus_p[ec_2][1][i]*d_minus_p[ec_2][1][i] + d_minus_p[ec_2][0][i]*d_minus_p[ec_2][0][i])
                                                 + (((d_minus_p[ec_2][1][i-s_2])*(d_minus_p[ec_2][1][i-s_2]))
                                                    + (d_minus_p[ec_2][0][i-s_2])*(d_minus_p[ec_2][0][i-s_2]))
                                                       + (((d_minus_p[ec_2][1][i+s_ec])*(d_minus_p[ec_2][1][i+s_ec]))
                                                          + (d_minus_p[ec_2][0][i+s_ec])*(d_minus_p[ec_2][0][i+s_ec]))
                                                             + ((d_minus_p[ec_2][1][i-s_2+s_ec])*(d_minus_p[ec_2][1][i-s_2+s_ec]))
                                                                + (d_minus_p[ec_2][0][i-s_2+s_ec])*(d_minus_p[ec_2][0][i-s_2+s_ec]));
                      const double frac_change = calc_nonlinear_inveps(dsqr_here,
                         s->inveps[ec][d_ec][i], s->kerr[ec][i]);
                      f[ec][0][i] = frac_change*(((s->inveps[ec][d_ec][i])*d_minus_p[ec][0][i])
                         + (0.25*((s->inveps[ec][d_1][i])*(d_minus_p[ec_1][0][i] + (d_minus_p[ec_1][0][i-s_1]) + (d_minus_p[ec_1][0][i+s_ec]) + d_minus_p[ec_1][0][i-s_1+s_ec]))) + 0.25*((s->inveps[ec][d_2][i])*(d_minus_p[ec_2][0][i] + (d_minus_p[ec_2][0][i-s_2]) + (d_minus_p[ec_2][0][i+s_ec]) + d_minus_p[ec_2][0][i-s_2+s_ec])));
                      f[ec][1][i] = frac_change*(((s->inveps[ec][d_ec][i])*d_minus_p[ec][1][i])
                         + (0.25*((s->inveps[ec][d_1][i])*(d_minus_p[ec_1][1][i] + (d_minus_p[ec_1][1][i-s_1]) + (d_minus_p[ec_1][1][i+s_ec]) + d_minus_p[ec_1][1][i-s_1+s_ec]))) + 0.25*((s->inveps[ec][d_2][i])*(d_minus_p[ec_2][1][i] + (d_minus_p[ec_2][1][i-s_2]) + (d_minus_p[ec_2][1][i+s_ec]) + d_minus_p[ec_2][1][i-s_2+s_ec])));
                    }
                  }
                }
              }
            }
          } else { // not s->kerr[ec]
            DOCMP {
              if (num_any_direction[Z]==1) {
                const int sR = stride_any_direction[R];
                const int sZ = 1;