update_e_from_d_update.hpp

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

                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]);
                  f[ec][1][i] = frac_change*((s->inveps[ec][d_ec][i])*d_minus_p[ec][1][i]);
                }
              } 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]);
                    f[ec][1][i] = frac_change*((s->inveps[ec][d_ec][i])*d_minus_p[ec][1][i]);
                  }
                } 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]);
                      f[ec][1][i] = frac_change*((s->inveps[ec][d_ec][i])*d_minus_p[ec][1][i]);
                    }
                  }
                }
              }
            }
          } else { // not s->kerr[ec]
            DOCMP {
              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;
                  f[ec][cmp][i] = (s->inveps[ec][d_ec][i])*d_minus_p[ec][cmp][i];
                }
              } 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;
                    f[ec][cmp][i] = (s->inveps[ec][d_ec][i])*d_minus_p[ec][cmp][i];
                  }
                } 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;
                      f[ec][cmp][i] = (s->inveps[ec][d_ec][i])*d_minus_p[ec][cmp][i];
                    }
                  }
                }
              }
            }
          }
        }
      }
    }
  } else { // not f[Er][0]
    if (f[Ey][0]) {
      if (f[Ez][0]) {
        if (stride_any_direction[Z]) {
          // Am in 3D
          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+1)%3);
            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)%3);
            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 sX = stride_any_direction[X];
                      const int sY = stride_any_direction[Y];
                      const int sZ = 1;
                      for (int iX=0; iX<num_any_direction[X]; iX++) {
                        for (int iY=0; iY<num_any_direction[Y]; iY++) {
                          const int i = yee_idx + iX*sX + iY*sY;
                          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[X]==1) {
                        const int sX = stride_any_direction[X];
                        const int sY = stride_any_direction[Y];
                        const int sZ = 1;
                        for (int iY=0; iY<num_any_direction[Y]; iY++) {
                          for (int iZ=0; iZ<num_any_direction[Z]; iZ++) {
                            const int i = yee_idx + iY*sY + 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][