update_from_e.hpp

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

                             + op->P[ec_1][0][i-s_1])) + (((np->P[ec_1][0][i+s_ec]) +
                                op->P[ec_1][0][i+s_ec])*((np->P[ec_1][0][i+s_ec]) + op->P[ec_1][0][i+s_ec]))
                                   + ((np->P[ec_1][0][i-s_1+s_ec]) + op->P[ec_1][0][i-s_1+s_ec])*((np->P[ec_1][0][i-s_1+s_ec])
                                      + op->P[ec_1][0][i-s_1+s_ec]))) + 0.25*((((np->P[ec_2][0][i])
                                         + op->P[ec_2][0][i])*((np->P[ec_2][0][i]) + op->P[ec_2][0][i]))
                                            + (((np->P[ec_2][0][i-s_2]) + op->P[ec_2][0][i-s_2])*((np->P[ec_2][0][i-s_2])
                                               + op->P[ec_2][0][i-s_2])) + (((np->P[ec_2][0][i+s_ec]) +
                                                  op->P[ec_2][0][i+s_ec])*((np->P[ec_2][0][i+s_ec]) + op->P[ec_2][0][i+s_ec]))
                                                     + ((np->P[ec_2][0][i-s_2+s_ec]) + op->P[ec_2][0][i-s_2+s_ec])*((np->P[ec_2][0][i-s_2+s_ec])
                                                        + op->P[ec_2][0][i-s_2+s_ec]));
                    const double Pdot_sqr = (((np->P[ec][0][i]) - (op->P[ec][0][i]))*((np->P[ec][0][i])
                       - (op->P[ec][0][i]))) + (0.25*((((np->P[ec_1][0][i]) - (op->P[ec_1][0][i]))*((np->P[ec_1][0][i])
                          - (op->P[ec_1][0][i]))) + (((np->P[ec_1][0][i-s_1]) - (op->P[ec_1][0][i-s_1]))*((np->P[ec_1][0][i-s_1])
                             - (op->P[ec_1][0][i-s_1]))) + (((np->P[ec_1][0][i+s_ec])
                                - (op->P[ec_1][0][i+s_ec]))*((np->P[ec_1][0][i+s_ec]) - (op->P[ec_1][0][i+s_ec]))) + ((np->P[ec_1][0][i-s_1+s_ec]) - (op->P[ec_1][0][i-s_1+s_ec]))*((np->P[ec_1][0][i-s_1+s_ec]) - (op->P[ec_1][0][i-s_1+s_ec])))) + 0.25*((((np->P[ec_2][0][i]) - (op->P[ec_2][0][i]))*((np->P[ec_2][0][i]) - (op->P[ec_2][0][i]))) + (((np->P[ec_2][0][i-s_2]) - (op->P[ec_2][0][i-s_2]))*((np->P[ec_2][0][i-s_2]) - (op->P[ec_2][0][i-s_2]))) + (((np->P[ec_2][0][i+s_ec]) - (op->P[ec_2][0][i+s_ec]))*((np->P[ec_2][0][i+s_ec]) - (op->P[ec_2][0][i+s_ec]))) + ((np->P[ec_2][0][i-s_2+s_ec]) - (op->P[ec_2][0][i-s_2+s_ec]))*((np->P[ec_2][0][i-s_2+s_ec]) - (op->P[ec_2][0][i-s_2+s_ec])));
                    const double energy_P_here = (om_sqr*P_sqr + (1.0/(Courant*Courant))*Pdot_sqr)*oo_ep_om_psqr;
                    np->s[ec][i] = -(energy_here - energy_P_here)*(om_psqr*oo_e_sat);
                    op->P[ec][0][i] = funinv*((2-om_sqr)*np->P[ec][0][i] + (0.5*g-1)*op->P[ec][0][i]
                       + np->s[ec][i]*f[ec][0][i]);
                  }
                }
              }
            }
          } else { // not np->saturation_factor
            for (int i=0;i<ntot;i++) {
              np->energy[ec][i] += (0.5)*(np->P[ec][0][i] - op->P[ec][0][i])*f[ec][0][i];
              op->P[ec][0][i] = funinv*((2-om_sqr)*np->P[ec][0][i] + (0.5*g-1)*op->P[ec][0][i]
                 + np->s[ec][i]*f[ec][0][i]);
            }
          }
        }
      }
    } else { // not is_real
      if (pol) {
        for (polarization *np=pol,*op=olpol; np; np=np->next,op=op->next) {
          const double om_psqr = fabs(op->pb->saturated_sigma);
          const double oo_ep_om_psqr = 0.5/(om_psqr);
          const double oo_e_sat = 1.0/op->pb->energy_saturation;
          const double om_sqr = op->pb->omeganot*op->pb->omeganot;
          const double g = op->pb->gamma;
          const double funinv = 1.0/(1+0.5*g);
          if (np->saturation_factor) {
            for (int i=0;i<ntot;i++) {
              np->energy[ec][i] += (0.5)*(np->P[ec][0][i] - op->P[ec][0][i])*f[ec][0][i];
              np->energy[ec][i] += (0.5)*(np->P[ec][1][i] - op->P[ec][1][i])*f[ec][1][i];
            }
            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 energy_here = (np->energy[ec][i]) + (0.25*((np->energy[ec_1][i])
                   + (np->energy[ec_1][i-s_1]) + (np->energy[ec_1][i+s_ec])
                      + np->energy[ec_1][i-s_1+s_ec])) + 0.25*((np->energy[ec_2][i])
                         + (np->energy[ec_2][i-s_2]) + (np->energy[ec_2][i+s_ec])
                            + np->energy[ec_2][i-s_2+s_ec]);
                const double P_sqr = ((((np->P[ec][1][i]) + op->P[ec][1][i])*((np->P[ec][1][i])
                   + op->P[ec][1][i])) + ((np->P[ec][0][i]) + op->P[ec][0][i])*((np->P[ec][0][i])
                      + op->P[ec][0][i])) + (0.25*(((((np->P[ec_1][1][i]) + op->P[ec_1][1][i])*((np->P[ec_1][1][i])
                         + op->P[ec_1][1][i])) + ((np->P[ec_1][0][i]) + op->P[ec_1][0][i])*((np->P[ec_1][0][i])
                            + op->P[ec_1][0][i])) + ((((np->P[ec_1][1][i-s_1]) + op->P[ec_1][1][i-s_1])*((np->P[ec_1][1][i-s_1])
                               + op->P[ec_1][1][i-s_1])) + ((np->P[ec_1][0][i-s_1]) + op->P[ec_1][0][i-s_1])*((np->P[ec_1][0][i-s_1])
                                  + op->P[ec_1][0][i-s_1])) + ((((np->P[ec_1][1][i+s_ec]) + op->P[ec_1][1][i+s_ec])*((np->P[ec_1][1][i+s_ec]) + op->P[ec_1][1][i+s_ec])) + ((np->P[ec_1][0][i+s_ec]) + op->P[ec_1][0][i+s_ec])*((np->P[ec_1][0][i+s_ec]) + op->P[ec_1][0][i+s_ec])) + (((np->P[ec_1][1][i-s_1+s_ec]) + op->P[ec_1][1][i-s_1+s_ec])*((np->P[ec_1][1][i-s_1+s_ec]) + op->P[ec_1][1][i-s_1+s_ec])) + ((np->P[ec_1][0][i-s_1+s_ec]) + op->P[ec_1][0][i-s_1+s_ec])*((np->P[ec_1][0][i-s_1+s_ec]) + op->P[ec_1][0][i-s_1+s_ec]))) + 0.25*(((((np->P[ec_2][1][i]) + op->P[ec_2][1][i])*((np->P[ec_2][1][i]) + op->P[ec_2][1][i])) + ((np->P[ec_2][0][i]) + op->P[ec_2][0][i])*((np->P[ec_2][0][i]) + op->P[ec_2][0][i])) + ((((np->P[ec_2][1][i-s_2]) + op->P[ec_2][1][i-s_2])*((np->P[ec_2][1][i-s_2]) + op->P[ec_2][1][i-s_2])) + ((np->P[ec_2][0][i-s_2]) + op->P[ec_2][0][i-s_2])*((np->P[ec_2][0][i-s_2]) + op->P[ec_2][0][i-s_2])) + ((((np->P[ec_2][1][i+s_ec]) + op->P[ec_2][1][i+s_ec])*((np->P[ec_2][1][i+s_ec]) + op->P[ec_2][1][i+s_ec])) + ((np->P[ec_2][0][i+s_ec]) + op->P[ec_2][0][i+s_ec])*((np->P[ec_2][0][i+s_ec]) + op->P[ec_2][0][i+s_ec])) + (((np->P[ec_2][1][i-s_2+s_ec]) + op->P[ec_2][1][i-s_2+s_ec])*((np->P[ec_2][1][i-s_2+s_ec]) + op->P[ec_2][1][i-s_2+s_ec])) + ((np->P[ec_2][0][i-s_2+s_ec]) + op->P[ec_2][0][i-s_2+s_ec])*((np->P[ec_2][0][i-s_2+s_ec]) + op->P[ec_2][0][i-s_2+s_ec]));
                const double Pdot_sqr = ((((np->P[ec][1][i]) - (op->P[ec][1][i]))*((np->P[ec][1][i])
                   - (op->P[ec][1][i]))) + ((np->P[ec][0][i]) - (op->P[ec][0][i]))*((np->P[ec][0][i])
                      - (op->P[ec][0][i]))) + (0.25*(((((np->P[ec_1][1][i]) - (op->P[ec_1][1][i]))*((np->P[ec_1][1][i]) - (op->P[ec_1][1][i]))) + ((np->P[ec_1][0][i]) - (op->P[ec_1][0][i]))*((np->P[ec_1][0][i]) - (op->P[ec_1][0][i]))) + ((((np->P[ec_1][1][i-s_1]) - (op->P[ec_1][1][i-s_1]))*((np->P[ec_1][1][i-s_1]) - (op->P[ec_1][1][i-s_1]))) + ((np->P[ec_1][0][i-s_1]) - (op->P[ec_1][0][i-s_1]))*((np->P[ec_1][0][i-s_1]) - (op->P[ec_1][0][i-s_1]))) + ((((np->P[ec_1][1][i+s_ec]) - (op->P[ec_1][1][i+s_ec]))*((np->P[ec_1][1][i+s_ec]) - (op->P[ec_1][1][i+s_ec]))) + ((np->P[ec_1][0][i+s_ec]) - (op->P[ec_1][0][i+s_ec]))*((np->P[ec_1][0][i+s_ec]) - (op->P[ec_1][0][i+s_ec]))) + (((np->P[ec_1][1][i-s_1+s_ec]) - (op->P[ec_1][1][i-s_1+s_ec]))*((np->P[ec_1][1][i-s_1+s_ec]) - (op->P[ec_1][1][i-s_1+s_ec]))) + ((np->P[ec_1][0][i-s_1+s_ec]) - (op->P[ec_1][0][i-s_1+s_ec]))*((np->P[ec_1][0][i-s_1+s_ec]) - (op->P[ec_1][0][i-s_1+s_ec])))) + 0.25*(((((np->P[ec_2][1][i]) - (op->P[ec_2][1][i]))*((np->P[ec_2][1][i]) - (op->P[ec_2][1][i]))) + ((np->P[ec_2][0][i]) - (op->P[ec_2][0][i]))*((np->P[ec_2][0][i]) - (op->P[ec_2][0][i]))) + ((((np->P[ec_2][1][i-s_2]) - (op->P[ec_2][1][i-s_2]))*((np->P[ec_2][1][i-s_2]) - (op->P[ec_2][1][i-s_2]))) + ((np->P[ec_2][0][i-s_2]) - (op->P[ec_2][0][i-s_2]))*((np->P[ec_2][0][i-s_2]) - (op->P[ec_2][0][i-s_2]))) + ((((np->P[ec_2][1][i+s_ec]) - (op->P[ec_2][1][i+s_ec]))*((np->P[ec_2][1][i+s_ec]) - (op->P[ec_2][1][i+s_ec]))) + ((np->P[ec_2][0][i+s_ec]) - (op->P[ec_2][0][i+s_ec]))*((np->P[ec_2][0][i+s_ec]) - (op->P[ec_2][0][i+s_ec]))) + (((np->P[ec_2][1][i-s_2+s_ec]) - (op->P[ec_2][1][i-s_2+s_ec]))*((np->P[ec_2][1][i-s_2+s_ec]) - (op->P[ec_2][1][i-s_2+s_ec]))) + ((np->P[ec_2][0][i-s_2+s_ec]) - (op->P[ec_2][0][i-s_2+s_ec]))*((np->P[ec_2][0][i-s_2+s_ec]) - (op->P[ec_2][0][i-s_2+s_ec])));
                const double energy_P_here = (om_sqr*P_sqr + (1.0/(Courant*Courant))*Pdot_sqr)*oo_ep_om_psqr;
                np->s[ec][i] = -(energy_here - energy_P_here)*(om_psqr*oo_e_sat);
                op->P[ec][0][i] = funinv*((2-om_sqr)*np->P[ec][0][i] + (0.5*g-1)*op->P[ec][0][i]
                   + np->s[ec][i]*f[ec][0][i]);
                op->P[ec][1][i] = funinv*((2-om_sqr)*np->P[ec][1][i] + (0.5*g-1)*op->P[ec][1][i]
                   + np->s[ec][i]*f[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 energy_here = (np->energy[ec][i]) + (0.25*((np->energy[ec_1][i])
                     + (np->energy[ec_1][i-s_1]) + (np->energy[ec_1][i+s_ec])
                        + np->energy[ec_1][i-s_1+s_ec])) + 0.25*((np->energy[ec_2][i])
                           + (np->energy[ec_2][i-s_2]) + (np->energy[ec_2][i+s_ec])
                              + np->energy[ec_2][i-s_2+s_ec]);
                  const double P_sqr = ((((np->P[ec][1][i]) + op->P[ec][1][i])*((np->P[ec][1][i])
                     + op->P[ec][1][i])) + ((np->P[ec][0][i]) + op->P[ec][0][i])*((np->P[ec][0][i])
                        + op->P[ec][0][i])) + (0.25*(((((np->P[ec_1][1][i]) + op->P[ec_1][1][i])*((np->P[ec_1][1][i])
                           + op->P[ec_1][1][i])) + ((np->P[ec_1][0][i]) + op->P[ec_1][0][i])*((np->P[ec_1][0][i])
                              + op->P[ec_1][0][i])) + ((((np->P[ec_1][1][i-s_1]) + op->P[ec_1][1][i-s_1])*((np->P[ec_1][1][i-s_1])
                                 + op->P[ec_1][1][i-s_1])) + ((np->P[ec_1][0][i-s_1]) + op->P[ec_1][0][i-s_1])*((np->P[ec_1][0][i-s_1])
                                    + op->P[ec_1][0][i-s_1])) + ((((np->P[ec_1][1][i+s_ec]) + op->P[ec_1][1][i+s_ec])*((np->P[ec_1][1][i+s_ec]) + op->P[ec_1][1][i+s_ec])) + ((np->P[ec_1][0][i+s_ec]) + op->P[ec_1][0][i+s_ec])*((np->P[ec_1][0][i+s_ec]) + op->P[ec_1][0][i+s_ec])) + (((np->P[ec_1][1][i-s_1+s_ec]) + op->P[ec_1][1][i-s_1+s_ec])*((np->P[ec_1][1][i-s_1+s_ec]) + op->P[ec_1][1][i-s_1+s_ec])) + ((np->P[ec_1][0][i-s_1+s_ec]) + op->P[ec_1][0][i-s_1+s_ec])*((np->P[ec_1][0][i-s_1+s_ec]) + op->P[ec_1][0][i-s_1+s_ec]))) + 0.25*(((((np->P[ec_2][1][i]) + op->P[ec_2][1][i])*((np->P[ec_2][1][i]) + op->P[ec_2][1][i])) + ((np->P[ec_2][0][i]) + op->P[ec_2][0][i])*((np->P[ec_2][0][i]) + op->P[ec_2][0][i])) + ((((np->P[ec_2][1][i-s_2]) + op->P[ec_2][1][i-s_2])*((np->P[ec_2][1][i-s_2]) + op->P[ec_2][1][i-s_2])) + ((np->P[ec_2][0][i-s_2]) + op->P[ec_2][0][i-s_2])*((np->P[ec_2][0][i-s_2]) + op->P[ec_2][0][i-s_2])) + ((((np->P[ec_2][1][i+s_ec]) + op->P[ec_2][1][i+s_ec])*((np->P[ec_2][1][i+s_ec]) + op->P[ec_2][1][i+s_ec])) + ((np->P[ec_2][0][i+s_ec]) + op->P[ec_2][0][i+s_ec])*((np->P[ec_2][0][i+s_ec]) + op->P[ec_2][0][i+s_ec])) + (((np->P[ec_2][1][i-s_2+s_ec]) + op->P[ec_2][1][i-s_2+s_ec])*((np->P[ec_2][1][i-s_2+s_ec]) + op->P[ec_2][1][i-s_2+s_ec])) + ((np->P[ec_2][0][i-s_2+s_ec]) + op->P[ec_2][0][i-s_2+s_ec])*((np->P[ec_2][0][i-s_2+s_ec]) + op->P[ec_2][0][i-s_2+s_ec]));
                  const double Pdot_sqr = ((((np->P[ec][1][i]) - (op->P[ec][1][i]))*((np->P[ec][1][i])
                     - (op->P[ec][1][i]))) + ((np->P[ec][0][i]) - (op->P[ec][0][i]))*((np->P[ec][0][i])
                        - (op->P[ec][0][i]))) + (0.25*(((((np->P[ec_1][1][i]) - (op->P[ec_1][1][i]))*((np->P[ec_1][1][i]) - (op->P[ec_1][1][i]))) + ((np->P[ec_1][0][i]) - (op->P[ec_1][0][i]))*((np->P[ec_1][0][i]) - (op->P[ec_1][0][i]))) + ((((np->P[ec_1][1][i-s_1]) - (op->P[ec_1][1][i-s_1]))*((np->P[ec_1][1][i-s_1]) - (op->P[ec_1][1][i-s_1]))) + ((np->P[ec_1][0][i-s_1]) - (op->P[ec_1][0][i-s_1]))*((np->P[ec_1][0][i-s_1]) - (op->P[ec_1][0][i-s_1]))) + ((((np->P[ec_1][1][i+s_ec]) - (op->P[ec_1][1][i+s_ec]))*((np->P[ec_1][1][i+s_ec]) - (op->P[ec_1][1][i+s_ec]))) + ((np->P[ec_1][0][i+s_ec]) - (op->P[ec_1][0][i+s_ec]))*((np->P[ec_1][0][i+s_ec]) - (op->P[ec_1][0][i+s_ec]))) + (((np->P[ec_1][1][i-s_1+s_ec]) - (op->P[ec_1][1][i-s_1+s_ec]))*((np->P[ec_1][1][i-s_1+s_ec]) - (op->P[ec_1][1][i-s_1+s_ec]))) + ((np->P[ec_1][0][i-s_1+s_ec]) - (op->P[ec_1][0][i-s_1+s_ec]))*((np->P[ec_1][0][i-s_1+s_ec]) - (op->P[ec_1][0][i-s_1+s_ec])))) + 0.25*(((((np->P[ec_2][1][i]) - (op->P[ec_2][1][i]))*((np->P[ec_2][1][i]) - (op->P[ec_2][1][i]))) + ((np->P[ec_2][0][i]) - (op->P[ec_2][0][i]))*((np->P[ec_2][0][i]) - (op->P[ec_2][0][i]))) + ((((np->P[ec_2][1][i-s_2]) - (op->P[ec_2][1][i-s_2]))*((np->P[ec_2][1][i-s_2]) - (op->P[ec_2][1][i-s_2]))) + ((np->P[ec_2][0][i-s_2]) - (op->P[ec_2][0][i-s_2]))*((np->P[ec_2][0][i-s_2]) - (op->P[ec_2][0][i-s_2]))) + ((((np->P[ec_2][1][i+s_ec]) - (op->P[ec_2][1][i+s_ec]))*((np->P[ec_2][1][i+s_ec]) - (op->P[ec_2][1][i+s_ec]))) + ((np->P[ec_2][0][i+s_ec]) - (op->P[ec_2][0][i+s_ec]))*((np->P[ec_2][0][i+s_ec]) - (op->P[ec_2][0][i+s_ec]))) + (((np->P[ec_2][1][i-s_2+s_ec]) - (op->P[ec_2][1][i-s_2+s_ec]))*((np->P[ec_2][1][i-s_2+s_ec]) - (op->P[ec_2][1][i-s_2+s_ec]))) + ((np->P[ec_2][0][i-s_2+s_ec]) - (op->P[ec_2][0][i-s_2+s_ec]))*((np->P[ec_2][0][i-s_2+s_ec]) - (op->P[ec_2][0][i-s_2+s_ec])));
                  const double energy_P_here = (om_sqr*P_sqr + (1.0/(Courant*Courant))*Pdot_sqr)*oo_ep_om_psqr;
                  np->s[ec][i] = -(energy_here - energy_P_here)*(om_psqr*oo_e_sat);
                  op->P[ec][0][i] = funinv*((2-om_sqr)*np->P[ec][0][i] + (0.5*g-1)*op->P[ec][0][i]
                     + np->s[ec][i]*f[ec][0][i]);
                  op->P[ec][1][i] = funinv*((2-om_sqr)*np->P[ec][1][i] + (0.5*g-1)*op->P[ec][1][i]
                     + np->s[ec][i]*f[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 energy_here = (np->energy[ec][i]) + (0.25*((np->energy[ec_1][i])
                       + (np->energy[ec_1][i-s_1]) + (np->energy[ec_1][i+s_ec])
                          + np->energy[ec_1][i-s_1+s_ec])) + 0.25*((np->energy[ec_2][i])