function navier_stokes_visualization(frame)

n=128;
dx=1./n;

u=load(strcat(strcat('u_',int2str(frame)),'.dat'));
v=load(strcat(strcat('v_',int2str(frame)),'.dat'));
p_hat=load(strcat(strcat('p_hat_',int2str(frame)),'.dat'));
%surf(p_hat);
%figure;
particles=load(strcat(strcat('particles_',int2str(frame)),'.dat'));

for i=0:n
    for j=0:n
        x=i*dx;
        y=j*dx;
        i_periodic=i;
        if i_periodic==n
            i_periodic=0;
        end
        i_periodic_m1=i_periodic-1;
        if i_periodic_m1==-1;
            i_periodic_m1=n-1;
        end
        j_periodic=j;
        if j_periodic==n
            j_periodic=0;
        end
        j_periodic_m1=j_periodic-1;
        if j_periodic_m1==-1;
            j_periodic_m1=n-1;
        end
        
        velocity_u(j+1,i+1)=.5*(u(i_periodic+1,j_periodic+1)+u(i_periodic+1,j_periodic_m1+1));
        velocity_v(j+1,i+1)=.5*(v(i_periodic+1,j_periodic+1)+v(i_periodic_m1+1,j_periodic+1));
    end
end

[xp,yp]=meshgrid(.5*dx:dx:(1-.5*dx),.5*dx:dx:(1-.5*dx));
contour(xp,yp,p_hat),hold on

[x,y]=meshgrid(0:dx:1,0:dx:1);
quiver(x,y,velocity_u,velocity_v),hold on

s=size(particles);
for p=0:(s(1)/2)-1
    px(p+1)=particles(2*p+1);
    py(p+1)=particles(2*p+2);
end

quarter_p=s(1)/8;

plot(px(1:quarter_p),py(1:quarter_p),'b*');
plot(px((1+quarter_p):(2*quarter_p)),py((1+quarter_p):(2*quarter_p)),'r*');
plot(px((1+2*quarter_p):(3*quarter_p)),py((1+2*quarter_p):(3*quarter_p)),'g*');
plot(px((1+3*quarter_p):(4*quarter_p)),py((1+3*quarter_p):(4*quarter_p)),'m*');

%figure;
%plot(px(1:quarter_p),py(1:quarter_p),'b*'),hold on
%plot(px((1+quarter_p):(2*quarter_p)),py((1+quarter_p):(2*quarter_p)),'r*'),hold on
%plot(px((1+2*quarter_p):(3*quarter_p)),py((1+2*quarter_p):(3*quarter_p)),'g*'),hold on
%plot(px((1+3*quarter_p):(4*quarter_p)),py((1+3*quarter_p):(4*quarter_p)),'m*'),hold on