update_from_e.hpp

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

                  }
                } 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 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 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[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 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[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][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]