commit before editing to check ComplexAmplitude sum

src/ipglasma.cpp

@@ -601,13 +601,11 @@ std::complex<double> IPGlasma::SumAmplitudes()
         for(int j =0; j < L; j += 1){
     
             int q1[2] = {i,j};
-      //      WilsonLine quark = GetWilsonLine(WilsonLineCoordinate(i, j));
+            WilsonLine quark = GetWilsonLine(WilsonLineCoordinate(i, j));
         //    double r1 = DistanceToOrigin(q1);
     
             for(int v =0; v < L; v +=1){
                 for(int w =0; w < 100; w +=1){
-                    double r = w*lattice_spacing;
-                    if (v>0 || j>0 || i>0) {break;}
                     // check that the dipole size is not zero
                     if (i == v && j == w) { continue;}
 
@@ -625,21 +623,19 @@ std::complex<double> IPGlasma::SumAmplitudes()
                     //TODO: is this correct? 
                     //TODO: remove case where r=0
                     // dipole size for wave function
-                    //double r = sqrt(pow(X(i)-X(v),2) + pow(Y(j)-Y(w),2));
+                    double r = sqrt(pow(X(i)-X(v),2) + pow(Y(j)-Y(w),2));
                     //cout << "# r in WF is " << r << endl;
 
                     // wave functions
                     double WF_u = WFGammaqqbarLongitudinal(r, z0, Q, u_quark); 
                     //cout << "# WF_u is " << WF_u << endl;
 
-                    cout << r << " " << WF_u << endl << std::flush;
-                    continue;
                     double WF_d = WFGammaqqbarLongitudinal(r, z0, Q, d_quark); 
                     double WF_c = WFGammaqqbarLongitudinal(r, z0, Q, c_quark); 
 
-    //                sum += (WF_u + WF_d + WF_c) * ComplexAmplitude(0, quark, q2)
-      //                      * exp( -i * ( std::inner_product(std::begin(vecL), std::end(vecL), std::begin(diffX), 0)
-        //                                       + std::inner_product(std::begin(vecDelta), std::end(vecDelta), std::begin(sumZX), 0) ) );
+                    sum += (WF_u + WF_d + WF_c) * ComplexAmplitude(0, quark, q2)
+                            * exp( -i * ( std::inner_product(std::begin(vecL), std::end(vecL), std::begin(diffX), 0)
+                                               + std::inner_product(std::begin(vecDelta), std::end(vecDelta), std::begin(sumZX), 0) ) );
                     
                     //TODO: Finish eq (12) from notes
                     //          Finish amplitude
@@ -672,8 +668,7 @@ double IPGlasma::WFGammaqqbarLongitudinal(double r, double z, double Q, const qu
     // fine structure constant set to 1/137
     double e = sqrt(4.*M_PI/137.036);
 
-    //double epsilon = sqrt(z*(1-z)*Q*Q + pow(quark.mass, 2));
-    double epsilon = 1;
+    double epsilon = sqrt(z*(1-z)*Q*Q + pow(quark.mass, 2));
     return (quark.e_charge * e * sqrt(NC) * 2*Q*Q*z*(1-z) * std::cyl_bessel_k(0, epsilon*r)*0.5/M_PI);
 }