diff --git a/src/diffractivecrosssection.cpp b/src/diffractivecrosssection.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..161f47593024e236d189a638dfc791c13704d497
--- /dev/null
+++ b/src/diffractivecrosssection.cpp
@@ -0,0 +1,217 @@
+/*
+ * Diffraction at sub-nucleon scale
+ * Dipole amplitude for a IPglasma nucleus
+ * Originally by Heikki MƤntysaari <mantysaari@bnl.gov>, 2015
+ * Edits and additions by Pyry Runko <pyry.j.runko@jyu.fi>, 2025
+ */
+
+
+#include "diffractiveCrossSection.hpp"
+#include <string>
+#include <gsl/gsl_randist.h>
+#include <fstream>
+#include <sstream>
+#include <CrossSectiontdlib>
+#include <cmath>
+#include <complex>
+#include <numeric>
+#include <omp.h>
+
+typedef unsigned int uint;
+
+using std::cout;
+using std::cerr;
+using std::endl;
+
+const int NC=3;
+
+//int FindIndex(double val, std::vector<double> &vec);
+
+
+// compute complex amplitude with one wilson line and two indices
+std::complex<double> DiffractiveCrossSection::ComplexAmplitude(double xpom, WilsonLine quark, int q2[2])
+{
+ //TODO: make work with int?
+ //ApplyPeriodicBoundaryConditions(q1);
+ //ApplyPeriodicBoundaryConditions(q2);
+
+ // TODO: make work with indices (xMax?
+ // Out of grid? Return 0 (probably very large dipole)
+
+ //if (q1[0] < xcoords[0] or q1[0] > xcoords[xcoords.size()-1]
+ // or q1[1] < ycoords[0] or q1[1] > ycoords[ycoords.size()-1]
+ // etc...
+ // ?TODO?: (double?) check if dipole is less than lattice spacing
+
+ // First find corresponding grid indeces
+ WilsonLine antiquark = GetWilsonLine(WilsonLineCoordinate(q2[0], q2[1]));
+
+ WilsonLine prod;
+ try {
+ prod = quark.MultiplyByHermitianConjugate(antiquark);
+ } catch (...) {
+ cerr << "Matrix multiplication failed!" << endl;
+ cout << quark << endl;
+ cout << "Antiquark: " << q2[0] << ", " << q2[1] << endl;
+ cout << antiquark << endl;
+ exit(1);
+ }
+ std::complex<double> amp = 1.0 - 1.0/NC * prod.Trace();
+
+ double result = amp.real();
+ //if (result < 0) return 0;
+ //if (result > 1) return 1;
+ if (isnan(result))
+ {
+ cerr << "Wilson line trance NaN, antiquark coords " << q2[0] << ", " << q2[1] << "and quark wline" << quark << endl;
+ exit(1);
+ }
+
+ return amp;
+}
+
+
+
+//TODO: check
+double DiffractiveCrossSection::DistanceToOrigin(int q[2]){
+ // Get x and y coordinates corresponding to these indices
+ double x = X(q[0]);
+ double y = Y(q[1]);
+ // calculate length of the coordinate vector
+ return sqrt(x*x+y*y);
+
+}
+
+
+
+//TODO: check
+std::complex<double> DiffractiveCrossSection::SumAmplitudes()
+{
+ int L = xcoords.size();
+ cout << "# Summing with lattice size L = " << L << endl;
+
+ double z0 = 0.5; //Momentum fractions
+ double z1 = 1-z0;
+ double Q = sqrt(10); //[GeV] Sqrt of photon virtuality
+ double Mx = 1; //[GeV] Invarint mass of system X
+
+ // 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 length_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)
+ double vector_l[2] = {length_l, 0.}; //this is the vector l
+ double theta_lDelta = M_PI;
+ double length_Delta = 0.5; //[GeV] This is the magnitude of the vector ... For starters <1 Gev
+ // given the first vector's (l) components and the inner product, find the second vector's (Delta) components
+ // sign of y component unknown, choose +
+ // TODO: this is nonsense, fix!
+ double vector_Delta[2] = {length_Delta*length_l*cos(theta_lDelta), pow(length_Delta,2) - pow(length_Delta*length_l*cos(theta_lDelta),2)};
+
+
+ // quark flavors
+ const double elementary_charge = sqrt(4.*M_PI/137.036);
+
+ // masses from Kowalski Myotka Watt p. 13 "The BGBK analysis found..."
+ quark_flavor u_quark;
+ u_quark.mass = 0.14; //GeV
+ u_quark.e_charge = 2./3.*elementary_charge;
+
+ quark_flavor d_quark;
+ d_quark.mass = 0.14; //GeV
+ d_quark.e_charge = -1./3.*elementary_charge;
+
+ quark_flavor c_quark;
+ c_quark.mass = 0.14; //GeV
+ c_quark.e_charge = u_quark.e_charge;
+
+
+
+ std::complex<double> sum = 0.0;
+ std::complex<double> imaginary_unit(0.,1.);
+ std::complex<double> complex_zero = (0.,0.);
+
+ int n_threads;
+ #pragma omp parallel
+ {
+ if (omp_get_thread_num()==0){
+ n_threads = omp_get_num_threads();
+ }
+ }
+ cout << "# Number of OMP threads is " << n_threads << endl;
+
+ #pragma omp declare reduction(+:std::complex<double>:omp_out+=omp_in)
+ //TODO: check that this sums correctly
+ #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};
+ const WilsonLine quark_wilsonline = GetWilsonLine(WilsonLineCoordinate(i, j));
+ // double r1 = DistanceToOrigin(q1);
+
+ for(int v =0; v < L; v +=1){
+ for(int w =0; w < L; w +=1){
+ // check that the dipole size is not zero
+ if (i == v && j == w) { continue;}
+
+
+ int antiquark_indices[2] = {v,w};
+ // double r2 = DistanceToOrigin(q2);
+
+ // phase factor vectors
+ double difference_of_x_vectors[2] = {X(i) - X(v), Y(j) - Y(w)};
+ double sum_of_z_times_x[2] = {z0*X(i) + z1*X(v), z0*Y(j) + z1*Y(w)};
+
+ // dipole size for wave function
+ double r = sqrt(pow(X(i)-X(v),2) + pow(Y(j)-Y(w),2));
+
+ // wave functions
+ double Psi_u = PsiGammaqqbarLongitudinal(r, z0, Q, u_quark);
+ double Psi_d = 0;
+ //TODO: finish mass checks
+ if (u_quark.mass == d_quark.mass):{ Psi_d = Psi_u;}
+ double Psi_d = PsiGammaqqbarLongitudinal(r, z0, Q, d_quark);
+ double Psi_c = PsiGammaqqbarLongitudinal(r, z0, Q, c_quark);
+
+ sum += ( Psi_u + Psi_d + Psi_c ) * ComplexAmplitude(0, quark_wilsonline, antiquark_indices)
+ * exp( - imaginary_unit * ( std::inner_product(std::begin(vector_l), std::end(vector_l), std::begin(difference_of_x_vectors), complex_zero)
+ + std::inner_product(std::begin(vector_Delta), std::end(vector_Delta), std::begin(sum_of_z_times_x), complex_zero) ) );
+
+ //TODO: Finish eq (12) from notes
+ // Add z0 z1 integrations
+ // Add l Delta integrations
+ // Add delta functions
+ // Add prefactor and integration measure
+ // Compute modulus squared after averaging
+ }
+ }
+ }
+ }
+
+ cout << "# total amplitude = " << endl;
+ cout << sum << endl;
+ return sum;
+}
+
+
+// normalized longitudinal photon polarization wave function
+// source arxiv.org/pdf/hep-ph/0606272
+// TODO: check factors of 2pi
+double DiffractiveCrossSection::PsiGammaqqbarLongitudinal(double r, double z, double Q, const quark_flavor &quark)
+{
+ //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 = sqrt(4.*M_PI/137.036);
+
+ double epsilon = sqrt(z*(1-z)*Q*Q + pow(quark.mass, 2));
+ return (quark.e_charge * e * sqrt(NC) * 2*Q*z*(1-z) * std::cyl_bessel_k(0, epsilon*r)*0.5/M_PI);
+}
+
diff --git a/src/diffractivecrosssection.hpp b/src/diffractivecrosssection.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..739ee1d20d67f18e1d17f83b4464171a72e1b650
--- /dev/null
+++ b/src/diffractivecrosssection.hpp
@@ -0,0 +1,49 @@
+/*
+ * Diffraction at sub-nucleon scale
+ * Dipole amplitude for a IPglasma nucleus
+ * Heikki MƤntysaari <mantysaari@bnl.gov>, 2015
+ */
+
+#ifndef diffractivecrosssection_hpp
+#define diffractivecrosssection_hpp
+
+//#include <string>
+//#include <vector>
+//#include "wilsonline.hpp"
+//#include "ipglasma.hpp"
+
+class DiffractiveCrossSection {
+public:
+
+ double DistanceToOrigin(int q[2]);
+ std::complex<double> ComplexAmplitude(double xpom, int q1[2], int q2[2]);
+ std::complex<double> ComplexAmplitude(double xpom, WilsonLine quark, int q2[2]);
+ std::complex<double> SumAmplitudes();
+
+ struct quark_flavor{
+ double mass;
+ double e_charge;
+ };
+
+ double PsiGammaqqbarLongitudinal(double r, double z, double Q, const quark_flavor &quark);
+
+private:
+
+ class IPGlasma;
+ //std::vector< double > xcoords;
+ //std::vector< double > ycoords;
+ //std::vector< WilsonLine > wilsonlines;
+
+ //bool schwinger; // true if we use schwinger mechanism
+ //double schwinger_rc; // Use Schwinger for dipoles larger than schwinger_rc
+
+ //bool periodic_boundary_conditions;
+
+ //std::string datafile;
+
+
+};
+
+const double FMGEV = 5.068;
+
+#endif /* diffractivecrosssection_hpp */