function spPlot()

%% Load data from labs to client
pmode lab2client out 1 set1
pmode lab2client out 2 set2
pmode lab2client out 3 set3 
pmode lab2client out 4 set4
pmode lab2client static 1 static

%% Move data from base workspace to function workspace
static = evalin('base', 'static');
set1 = evalin('base', 'set1');
set2 = evalin('base', 'set2');
set3 = evalin('base', 'set3');
set4 = evalin('base', 'set4');

%% Extract data from structures
xc = static.xc;
yc = static.yc;
Time = 0;
R = xc^2+yc^2;

m = static.m;
n = static.m;
M = static.M;
MM = static.MM;

plotpoints = static.plotpoints;

[pData(1:plotpoints).empty] = deal(0);

for i = 1:plotpoints
    pData(i).X = [set1(i).X set2(i).X set3(i).X set4(i).X];
    pData(i).Y = [set1(i).Y set2(i).Y set3(i).Y set4(i).Y];
    pData(i).u = [set1(i).u set2(i).u set3(i).u set4(i).u];
    pData(i).v = [set1(i).v set2(i).v set3(i).v set4(i).v];
    
    pData(i).A1 = [set1(i).A1 set2(i).A1 set3(i).A1 set4(i).A1];
    pData(i).A2 = [set1(i).A2 set2(i).A2 set3(i).A2 set4(i).A2];
    
    pData(i).t = set1(i).t;
   
end

%% Plot
for c = 1:plotpoints
  
    X = pData(c).X; % matrix of x-coordinate values
    Y = pData(c).Y; % matrix of y-coordinate values
    u = pData(c).u; % matrix of x-components of velocity
    v = pData(c).v; % matrix of y-components of velocity

    A1 = pData(c).A1;
    A2 = pData(c).A2; 

    t = pData(c).t;

    HH=0*X;
    AA=0*X;      

    AA(2:2:m,2:2:n)=(A2(2:2:m,2:2:n)+A2(2:2:m,1:2:n-1)+A2(1:2:m-1,2:2:n)+A2(1:2:m-1,1:2:n-1))/4;
    AA(2:2:m,2:2:n)=AA(2:2:m,2:2:n)+(A1(2:2:m,2:2:n)+A1(2:2:m,1:2:n-1)+A1(1:2:m-1,2:2:n)+A1(1:2:m-1,1:2:n-1))/4;
    AA(2:2:m,3:2:n-1)=(A1(2:2:m,3:2:n-1)+A1(1:2:m-1,3:2:n-1)+A1(1:2:m-1,2:2:n-2)+A1(2:2:m,2:2:n-2))/2;

    AA(3:2:m-1,2:2:n)=(A2(3:2:m-1,2:2:n)+A2(3:2:m-1,1:2:n-1)+A2(2:2:m-2,2:2:n)+A2(2:2:m-2,1:2:n-1))/2;
    AA(3:2:m-1,3:2:n-1)=(A2(3:2:m-1,3:2:n-1)+A2(2:2:m-2,3:2:n-1)+A2(3:2:m-1,2:2:n-2)+A2(2:2:m-2,2:2:n-2))/4;
    AA(3:2:m-1,3:2:n-1)=AA(3:2:m-1,3:2:n-1)+(A1(3:2:m-1,3:2:n-1)+A1(1:2:m-2,3:2:n-1)+A1(3:2:m-1,2:2:n-2)+A1(2:2:m-2,2:2:n-2))/4;
    HH(2:m,2:n)=M(2:m,2:n)./AA(2:m,2:n);


    Xc=sum(sum(M.*X))/MM;
    Yc=sum(sum(M.*Y))/MM;
    uc=sum(sum(M.*u))/MM;
    vc=sum(sum(M.*v))/MM;
    % This will do for now
    xc=[xc,Xc];
    yc=[yc,Yc];

    if max(size(R))<100;
        R=[R,Xc^2+Yc^2];    
        Time=[Time,t];
    else
        R=[R(2:100),Xc^2+Yc^2];  
        Time=[Time(2:100),t];
    end

    XX=X;
    YY=Y;

    ff=('velocity fields u-uc and v-vc -- both approximately planar fields');

    ff1=['surface plots of water elevation h+(x^2+y^2)/2, the parabolic boundary, and time  ' num2str(t)];
    clf

    % plot of water column above paraboloid
    subplot(2,2,1 )
    surf(XX,YY,HH)
    shading interp
    hold on
    contour3(XX,YY,HH,20,'k')
    hold on
    axis([-2 2 -2 2 0 2])
    axis square
    title('surface plot of water column above z=(x^2+y^2)/2')

    % plot of water comlumn absolute depth
    subplot(322)
    surf(XX-Xc,YY-Yc,u-uc)
    shading interp
    hold on
    surf(XX-Xc,YY-Yc,v-vc)
    shading interp
    contour3(XX-Xc,YY-Yc,u-uc,20,'k')
    title(ff)
    axis([-1.5 1.5 -1.5 1.5 -1.5 1.5])% -.5 .5])
    axis square

    
    
    subplot(3,2,4)
    plot(xc,yc,'r',Xc,Yc,'+k',0,0,'+k',[-2 2],[0 0],'k',[0 0],[-2 2],'k')
    axis([-1 1 -1 1])
    title('trajectory of the center of mass-initial position (1/sqrt(2) , 0) ')

    axis square

    subplot(223)
    ZZ1=((XX.^2+YY.^2)/2) ;
    ZZ2=HH+(XX.^2+YY.^2)/2;
    surf(XX,YY,ZZ1)
    shading interp
    hold on
    surf(XX,YY,ZZ2)
    shading interp
    hold on
    contour3(XX,YY,ZZ2,20,'k')
    hold on
    contour(XX,YY,ZZ2,20,'k')
    axis([-2 2 -2 2 0 2.5])
    axis square
    title(ff1)
   
    subplot(3,2,6)
    plot(Time,R)
    axis([min(Time),max(Time), min(R),max(R)])
    title('xc^2+yc^2 vs time')
    drawnow
    
    %pause(0.1)

end