diff --git a/src/ipglasma.cpp b/src/ipglasma.cpp
index 7f9a084be209e6c34c32e8e536a2b173c8d8cd87..97df1099db81a442f17547710ade0292b2f74b98 100644
--- a/src/ipglasma.cpp
+++ b/src/ipglasma.cpp
@@ -671,18 +671,26 @@ double IPGlasma::SumAmplitudes()
 
 // return longitudinal photon polarization wave function
 // source arxiv.org/pdf/hep-ph/0606272
-double IPGlasma::WFGammaqqbar(double r)
+// quark flavor f: 0=u, 1=d
+double IPGlasma::WFGammaqqbar(double r, int f)
 { 
     //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
+
+    // fine structure constant set to 1/137
+    double e = cmath.sqrt(4*cmath.pi*1/137)
+    double m = 0.16e-9 //GeV
+    double e_f;
+    if (f){ e_f = -1/3};
+    else { e_f = 2/3};
+
     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)
+    double eps = cmath.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 
 }