update_from_e.hpp

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

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

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 (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];
            }
            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][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[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][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[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])
                             + (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])