Function g_eval,c,zeta ; this function returns a dimensionless G for given values of c and zeta. Zeta should be zero
pi = 3.14159265
n = 1000000
n1 = randomn(seed,n) ; these are 5 n element random normal distributions used to account for the random nature of the field
n2 = randomn(seed,n)
n3 = randomn(seed,n)
n4 = randomn(seed,n)
n5 = randomn(seed,n)
G = exp(-0.5*(zeta)^2)/(4*sqrt(2)*(pi)^(1.5))*(((sqrt(1-c^2)*n5 + c*zeta)^3) + ( (1-c^2)*((n3^2 + n4^2) -0.5*(n1^2 +n2^2)*(sqrt(1-c^2)*n5 + c*zeta)) ) + (1-c^2)*(sqrt(2)*n2*n3*n4 + (1/sqrt(2))*(n3^2 - n4^2)*n1))
G = abs(G)
return, mean(G)
end