function f = EQ(rho,phi_x)
global V; 
[Nv,~] = size(V);
[Nx,Nz] = size(rho);
rho = reshape(repmat(rho,[Nv,1]),Nx,Nv,Nz);
phi_x = reshape(repmat(phi_x,[Nv,1]),Nx,Nv,Nz);
f = (rho/sqrt(2*pi)).*exp(-(repmat(V',[Nx,1,Nz]) + phi_x).^2/2);
