function p = initial_rho(x,y)
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here

rho1=1.1; rho2=0.5065; rho3=1.1; rho4=0.5065;
%u1=0; u2=0.8939; u3=0.8939; u4=0;
%v1=0; v2=0; v3=0.8939; v4=0.8939;
%p1=1.1; p2=0.35; p3=1.1; p4=0.35;

p=x*0;

aa = x>=0.5 & y>=0.5  % zone 1

bb = x<0.5 & y>=0.5  % zone 2

cc = x<0.5 & y<0.5   % zone 3

dd = x>=0.5 & y<0.5   % zone 4


p = rho1*aa+rho2*bb+rho3*cc+rho4*dd;

end