function [F,norm_F] = Global_Equil(phi)
global V; global eps;
[Nx,Nz] = size(phi);
[Nv,~] = size(V);
M = exp(-V'.^2/2)/sqrt(2*pi);
M = repmat(M,[Nx,1,Nz]);
F = exp(-phi/eps);
F = repmat(reshape(F,Nx*Nz,1),[1,Nv]);
F = permute(reshape(F, [Nx,Nz,Nv]),[1,3,2]);
F = F.*M;
global delx;global delv;
norm_F = permute(sum(sum(F))*delx*delv,[3,1,2]);
norm_F = permute(repmat(norm_F,[1,Nx,Nv]),[2,3,1]);
F = F./norm_F;