started writing photon wave function

src/ipglasma.cpp

@@ -11,6 +11,7 @@
 #include <fstream>
 #include <sstream>
 #include <cstdlib>
+#include <cmath>
 
 typedef unsigned int uint;
 
@@ -591,7 +592,6 @@ double IPGlasma::SumAmplitudes(int l)
     WilsonLine quark;
     for(int i =0; i < L; i += 1){
         for(int j =0; j < L; j += 1){
-            
             int q1[2] = {i,j};
             //TODO: pointer instead of copy?
             quark = GetWilsonLine(WilsonLineCoordinate(i, j));
@@ -628,17 +628,19 @@ double IPGlasma::SumAmplitudes()
     int L = xcoords.size();
     cout << "# Summing with lattice size L = " << L << endl;
     WilsonLine quark;
-    int thrownAway = 0;
-    int iWasHere = 0;
+    //int thrownAway = 0;
     #pragma omp parallel for reduction(+:sum)
     for(int i =0; i < L; i += 1){
         for(int j =0; j < L; j += 1){
+    
             int q1[2] = {i,j};
-            quark = GetWilsonLine(WilsonLineCoordinate(i, j))
+            quark = GetWilsonLine(WilsonLineCoordinate(i, j));
             double r1 = DistanceToOrigin(q1);
             //#pragma omp parallel for reduction(+:sum) collapse(2)
+    
             for(int k =0; k < L; k +=1){
                 for(int l =0; l < L; l +=1){
+    
                     int q2[2] = {k,l};
                     double r2 = DistanceToOrigin(q2);
                     if (r1 > 1. && r2 > 1.){
@@ -647,7 +649,7 @@ double IPGlasma::SumAmplitudes()
                         //cout << "#Quark coordinates " << X(q1[0]) << "," << Y(q1[1]) << " and " << X(q2[0]) << "," << Y(q2[1]) << endl;
                         //cout << "#Dipole size is " << r1-r2 << endl;
                         //cout << "#Dipole amplitude is " << Amplitude(0, q1, q2) << endl;
-                        thrownAway += 1;
+      //                  thrownAway += 1;
                         continue;
                 }
 
@@ -656,15 +658,31 @@ double IPGlasma::SumAmplitudes()
 
                    sum += Amplitude(0, q1, q2);
                    //sum += Amplitude(0, quark, q2);
-                   iWasHere += 1;  
                 }
             }
         }
     }
     cout << "# total sum = " << sum << endl;
-    cout << "# i counted " << iWasHere << " values" << endl;
-    cout << "# threw away " << thrownAway << " values (" << 100.*thrownAway*pow(L,-4) << "%)" << endl;
+ //   cout << "# i counted " << iWasHere << " values" << endl;
+   // cout << "# threw away " << thrownAway << " values (" << 100.*thrownAway*pow(L,-4) << "%)" << endl;
     return sum;
 }
 
 
+// return longitudinal photon polarization wave function
+// source arxiv.org/pdf/hep-ph/0606272
+double IPGlasma::WFGammaqqbar(double r)
+{ 
+    //Assuming r is the transverse size of the dipole
+    //constants: e_f (quark charge), e, N_C, 
+    //variables: r, Q, z, epsilon
+    //epsilon = sqrt(z(1-z)Q^2-m^2_f
+    double z = 0.5 //GeV^0
+    double QQ = 10 //GeV^2
+    double m_f = 0.16e-9 //GeV
+    double Delta = 0.5 //GeV
+    double l = a/100 // << a^-1 (a = lattice const)
+    double eps = sqrt(z*(1-z)*QQ*m_f*m_f)
+ 
+    return e_f*e*cmath.sqrt(NC)*2*Q*z*(1-z)*std::cyl_bessel_k(0, eps*r)0.5/cmath.pi 
+}