update_e_from_d_update.hpp

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

                  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])));
                      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])));
                    }
                  }
                }
              }
            }
          } 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])
                     + 0.25*((s->inveps[ec][d_1][i])*(d_minus_p[ec_1][cmp][i]
                        + (d_minus_p[ec_1][cmp][i-s_1]) + (d_minus_p[ec_1][cmp][i+s_ec])
                           + d_minus_p[ec_1][cmp][i-s_1+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;
                    f[ec][cmp][i] = ((s->inveps[ec][d_ec][i])*d_minus_p[ec][cmp][i])
                       + 0.25*((s->inveps[ec][d_1][i])*(d_minus_p[ec_1][cmp][i]
                          + (d_minus_p[ec_1][cmp][i-s_1]) + (d_minus_p[ec_1][cmp][i+s_ec])
                             + d_minus_p[ec_1][cmp][i-s_1+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;
                      f[ec][cmp][i] = ((s->inveps[ec][d_ec][i])*d_minus_p[ec][cmp][i])
                         + 0.25*((s->inveps[ec][d_1][i])*(d_minus_p[ec_1][cmp][i]
                            + (d_minus_p[ec_1][cmp][i-s_1]) + (d_minus_p[ec_1][cmp][i+s_ec])
                               + d_minus_p[ec_1][cmp][i-s_1+s_ec]));
                    }
                  }
                }
              }
            }
          }
        }
      } else { // not 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_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_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_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_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_2][i])*(d_minus_p[ec_2][1][i] +