bessel function plot modified code

src/ipglasma.cpp

@@ -14,6 +14,7 @@
 #include <cmath>
 #include <complex>
 #include <numeric>
+#include <omp.h>
 
 typedef unsigned int uint;
 
@@ -509,7 +510,7 @@ std::complex<double> IPGlasma::ComplexAmplitude(double xpom, WilsonLine quark, i
     double result = amp.real();
     //if (result < 0) return 0;
     //if (result > 1) return 1;
-    //if (isnan(result))
+    if (isnan(result))
     {
 	cerr << "Wilson line trance NaN, antiquark coords " << q2[0] << ", " << q2[1] << "and quark wline" << quark << endl;
         exit(1);
@@ -545,6 +546,8 @@ std::complex<double> IPGlasma::SumAmplitudes()
 
     // in GeV, see LoadData
     double lattice_spacing = xcoords[1]-xcoords[0];
+    cout << "# Lattice spacing is " << lattice_spacing << " GeV^-1." << endl;
+
 
     //TODO: loop over all values
     double l = 1./(100*lattice_spacing); //[GeV] This is the magnitude of the vector ... whose direction is chosen to be in the x direction (1,0) For starters << a (lattice const)
@@ -557,19 +560,19 @@ std::complex<double> IPGlasma::SumAmplitudes()
 
 
     // quark flavors
-    const double elementary_charge = sqrt(4.*M_PI*1./137.036); 
+    const double elementary_charge = sqrt(4.*M_PI/137.036); 
 
     //TODO: check masses
     quark_flavor u_quark;
-    u_quark.mass = 2e-3; //GeV
+    u_quark.mass = 0.14; //GeV
     u_quark.e_charge = 2./3.*elementary_charge;
 
     quark_flavor d_quark;
-    d_quark.mass = 2e-3;
+    d_quark.mass = 0.14; //GeV
     d_quark.e_charge = -1./3.*elementary_charge;
 
     quark_flavor c_quark;
-    c_quark.mass = 1.;
+    c_quark.mass = 0.14; //GeV
     c_quark.e_charge = u_quark.e_charge;
 
 
@@ -577,40 +580,66 @@ std::complex<double> IPGlasma::SumAmplitudes()
     std::complex<double> sum = 0.0;
     //double sum = 0.0;
 
-
+    cout << "# r (lattice units)          WF_u" << endl << std::flush;
+   // bool is_parallel = false;
     #pragma omp declare reduction(+:std::complex<double>:omp_out+=omp_in)
     //TODO: check that this sums correctly
-    #pragma omp parallel for reduction(+:sum)
-    //error: user defined reduction not found for ‘sum’
-
+    //#pragma omp parallel for reduction(+:sum)
     for(int i =0; i < L; i += 1){
+    
+    //if (i ==0){
+      //  if (omp_get_num_threads() > 1) {
+        //    #pragma omp critical
+          //  { is_parallel = true;}
+            
+        //}
+        //if (is_parallel) { std::cout << "# OpenMP is running with " << omp_get_num_threads() << "  threads!" << std::endl;} 
+        //else { std::cout << "# OpenMP is not enabled, or running with 1 thread." << std::endl;}
+   // }
+    
+
         for(int j =0; j < L; j += 1){
     
             int q1[2] = {i,j};
-            WilsonLine quark = GetWilsonLine(WilsonLineCoordinate(i, j));
-            double r1 = DistanceToOrigin(q1);
+      //      WilsonLine quark = GetWilsonLine(WilsonLineCoordinate(i, j));
+        //    double r1 = DistanceToOrigin(q1);
     
             for(int v =0; v < L; v +=1){
-                for(int w =0; w < L; w +=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;}
+
+
                     int q2[2] = {v,w};
-                    double r2 = DistanceToOrigin(q2);
+          //          double r2 = DistanceToOrigin(q2);
 
                     // phase factor vectors
                     double diffX[2] = {X(i) - X(v), Y(j) - Y(w)};
+                    //cout << "# diffX is " << diffX << endl;
+
                     double sumZX[2] = {z0*X(i) + z1*X(v), z0*Y(j) + z1*Y(w)};
-                    
-                    // position for wave function
-                    double r = sqrt(pow(X(i)-X(v),2) + pow(Y(j)-Y(w),2));
+                    //cout << "# sumZX is " << sumZX << endl;
+
+                    //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));
+                    //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
@@ -624,27 +653,27 @@ std::complex<double> IPGlasma::SumAmplitudes()
             }
         }
     }
-    cout << "# total amplitude = " << sum << endl;
+    
+    cout << "# total amplitude = " << endl;
+    cout << sum << endl;
     return sum;
 }
 
 
-// return longitudinal photon polarization wave function
+// longitudinal photon polarization wave function
 // source arxiv.org/pdf/hep-ph/0606272
-// quark flavor f: 0=u, 1=d
 double IPGlasma::WFGammaqqbarLongitudinal(double r, double z, double Q, const quark_flavor &quark)
 { 
-    //Assuming r is the transverse size of the dipole
+    //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 lattice_spacing = xcoords[1]-xcoords[0];
     // fine structure constant set to 1/137
-    double e = sqrt(4.*M_PI*1./137.036);
+    double e = sqrt(4.*M_PI/137.036);
 
-    double epsilon = sqrt(z*(1-z)*Q*Q * pow(quark.mass, 2));
- 
+    //double epsilon = sqrt(z*(1-z)*Q*Q + pow(quark.mass, 2));
+    double epsilon = 1;
     return (quark.e_charge * e * sqrt(NC) * 2*Q*Q*z*(1-z) * std::cyl_bessel_k(0, epsilon*r)*0.5/M_PI);
 }