update_from_e.hpp

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

                      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 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][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 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 num_any_direction[X]==1
                      if (num_any_direction[Y]==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 iZ=0; iZ<num_any_direction[Z]; iZ++) {
                            const int i = yee_idx + iX*sX + 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][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 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 num_any_direction[Y]==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++) {
                            for (int iZ=0; iZ<num_any_direction[Z]; iZ++) {
                              const int i = yee_idx + iX*sX + iY*sY + 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][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 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]);
                            }
                          }
                        }
                      }
                    }