function data = spOneStep(data)

%% Data in from last step

% This data is constant
m = data.m; % number of grid lines in x-direction
n = data.n; % number of grid lines in y-direction
b = data.b; % rotation coefficient
a = data.a; % bottom friction coefficient
mu = data.mu; % eddy viscocity
dt = data.dt; % time step
%dx = data.dx;

m1 = data.m1; % matrices of masses
m2 = data.m2; 
M = data.M;

% This data is changed
X = data.X; % matrix of x-coordinate values
Y = data.Y; % matrix of y-coordinate values
u = data.u; % matrix of x-components of velocity
v = data.v; % matrix of y-components of velocity


%% Calculate properties along edges of triangles
% positions differences

DxX1(1:m,1:2:n-1)=X(2:m+1,1:2:n-1)-X(1:m,1:2:n-1);
DxY1(1:m,1:2:n-1)=Y(2:m+1,1:2:n-1)-Y(1:m,1:2:n-1);
DxX1(1:m,2:2:n)=X(2:m+1,3:2:n+1)-X(1:m,3:2:n+1);
DxY1(1:m,2:2:n)=Y(2:m+1,3:2:n+1)-Y(1:m,3:2:n+1);

DxX2(1:m,1:2:n-1)=X(2:m+1,2:2:n)-X(1:m,2:2:n);
DxY2(1:m,1:2:n-1)=Y(2:m+1,2:2:n)-Y(1:m,2:2:n);
DxX2(1:m,2:2:n)=X(2:m+1,2:2:n)-X(1:m,2:2:n);
DxY2(1:m,2:2:n)=Y(2:m+1,2:2:n)-Y(1:m,2:2:n);

DyX1(2:2:m,1:n)=X(2:2:m,2:n+1)-X(2:2:m,1:n);
DyY1(2:2:m,1:n)=Y(2:2:m,2:n+1)-Y(2:2:m,1:n);
DyX1(1:2:m-1,1:n)=X(2:2:m,2:n+1)-X(2:2:m,1:n);
DyY1(1:2:m-1,1:n)=Y(2:2:m,2:n+1)-Y(2:2:m,1:n);

DyX2(1:2:m-1,1:n)=X(1:2:m-1,2:n+1)-X(1:2:m-1,1:n);
DyY2(1:2:m-1,1:n)=Y(1:2:m-1,2:n+1)-Y(1:2:m-1,1:n);

DyX2(2:2:m,1:n)=X(3:2:m+1,2:n+1)-X(3:2:m+1,1:n);
DyY2(2:2:m,1:n)=Y(3:2:m+1,2:n+1)-Y(3:2:m+1,1:n);

% calculate this data to avoid storing it

A1=(DxX1.*DyY1-DyX1.*DxY1)/2;
A2=(DxX2.*DyY2-DyX2.*DxY2)/2;

h1=m1./A1;
h2=m2./A2;

% velocity differences
Dxu1(1:m,1:2:n-1)=u(2:m+1,1:2:n-1)-u(1:m,1:2:n-1);
Dxv1(1:m,1:2:n-1)=v(2:m+1,1:2:n-1)-v(1:m,1:2:n-1);
Dxu1(1:m,2:2:n)=u(2:m+1,3:2:n+1)-u(1:m,3:2:n+1);
Dxv1(1:m,2:2:n)=v(2:m+1,3:2:n+1)-v(1:m,3:2:n+1);

Dxu2(1:m,1:2:n-1)=u(2:m+1,2:2:n)-u(1:m,2:2:n);
Dxv2(1:m,1:2:n-1)=v(2:m+1,2:2:n)-v(1:m,2:2:n);
Dxu2(1:m,2:2:n)=u(2:m+1,2:2:n)-u(1:m,2:2:n);
Dxv2(1:m,2:2:n)=v(2:m+1,2:2:n)-v(1:m,2:2:n);

Dyu1(2:2:m,1:n)=u(2:2:m,2:n+1)-u(2:2:m,1:n);
Dyv1(2:2:m,1:n)=v(2:2:m,2:n+1)-v(2:2:m,1:n);
Dyu1(1:2:m-1,1:n)=u(2:2:m,2:n+1)-u(2:2:m,1:n);
Dyv1(1:2:m-1,1:n)=v(2:2:m,2:n+1)-v(2:2:m,1:n);

Dyu2(1:2:m-1,1:n)=u(1:2:m-1,2:n+1)-u(1:2:m-1,1:n);
Dyv2(1:2:m-1,1:n)=v(1:2:m-1,2:n+1)-v(1:2:m-1,1:n);

Dyu2(2:2:m,1:n)=u(3:2:m+1,2:n+1)-u(3:2:m+1,1:n);
Dyv2(2:2:m,1:n)=v(3:2:m+1,2:n+1)-v(3:2:m+1,1:n);

%% Pressure terms in x momentum eqn.

p1RLx=(h1.^2).*(DyY1)/4;
p2RLx=(h2.^2).*(DyY2)/4;
p1UDx=(h1.^2).*(DxY1)/4;
p2UDx=(h2.^2).*(DxY2)/4;

P1rxx=0*X;
P1rxx(1:m,1:n)=p1RLx;

P1Lxx=0*X;
P1Lxx(2:m+1,1:n)=p1RLx;

P2rxx=0*X;
P2rxx(1:m,1:n)=p2RLx;

P2Lxx=0*X;
P2Lxx(2:m+1,1:n)=p2RLx;

P2Uxy=0*Y;
P2Uxy(1:m,1:n)=p2UDx;

P2Dxy=0*Y;
P2Dxy(1:m,2:n+1)=p2UDx;

P1Uxy=0*Y;
P1Uxy(1:m,1:n)=p1UDx;

P1Dxy=0*Y;
P1Dxy(1:m,2:n+1) = p1UDx;

% Pressure Forces in x eqn 

PPx=0*X;
PPx(:,1)=PPx(:,1)+ dt*(P1Lxx(:,1)-P1rxx(:,1));
PPx(:,3:2:n+1)=PPx(:,3:2:n+1)+dt*( P1Lxx(:,3:2:n+1)-P1rxx(:,3:2:n+1));
PPx(:,3:2:n+1)=PPx(:,3:2:n+1)+dt*(P1Lxx(:,2:2:n)-P1rxx(:,2:2:n));
PPx(:,2:2:n)=PPx(:,2:2:n)+dt*(P2Lxx(:,1:2:n-1)-P2rxx(:,1:2:n-1));
PPx(:,2:2:n)=PPx(:,2:2:n)+dt*(P2Lxx(:,2:2:n)-P2rxx(:,2:2:n));
PPx(1,:)=PPx(1,:)-dt*(P2Dxy(1,:)-P2Uxy(1,:));
PPx(3:2:m+1,:)=PPx(3:2:m+1,:)-dt*(P2Dxy(3:2:m+1,:)-P2Uxy(3:2:m+1,:));
PPx(3:2:m+1,:)=PPx(3:2:m+1,:)-dt*(P2Dxy(2:2:m,:)-P2Uxy(2:2:m,:));
PPx(2:2:m,:)=PPx(2:2:m,:)-dt*(P1Dxy(1:2:m-1,:)-P1Uxy(1:2:m-1,:));
PPx(2:2:m,:)=PPx(2:2:m,:)-dt*(P1Dxy(2:2:m,:)-P1Uxy(2:2:m,:));


%% Pressure terms in y momentum eqn.

p2udy=(h2.^2).*(DxX2)/4;
p1udy=(h1.^2).*(DxX1)/4;
p2RLy=(h2.^2).*(DyX2)/4;
p1RLy=(h1.^2).*(DyX1)/4;

P2uyy=0*Y;
P2uyy(1:m,1:n)= p2udy;

P2dyy=0*Y;
P2dyy(1:m,2:n+1)=p2udy;

P1uyy=0*Y;
P1uyy(1:m,1:n)= p1udy;

P1dyy=0*Y;
P1dyy(1:m,2:n+1)= p1udy;


P1ryx=0*X;
P1ryx(1:m,1:n)= p1RLy ;

P1Lyx=0*X;
P1Lyx(2:m+1,1:n)=p1RLy; 

P2ryx=0*X;
P2ryx(1:m,1:n)=p2RLy;

P2Lyx=0*X;
P2Lyx(2:m+1,1:n)= p2RLy;


% Pressure Forces in y eqn

PPy=0*Y;
PPy(1,:)=PPy(1,:)+dt*(P2dyy(1,:)-P2uyy(1,:));

PPy(3:2:m+1,:)=PPy(3:2:m+1,:)+dt*(P2dyy(3:2:m+1,:)-P2uyy(3:2:m+1,:));
PPy(3:2:m+1,:)=PPy(3:2:m+1,:)+dt*(P2dyy(2:2:m,:)-P2uyy(2:2:m,:));

PPy(2:2:m,:)=PPy(2:2:m,:)+dt*(P1dyy(1:2:m-1,:)-P1uyy(1:2:m-1,:));
PPy(2:2:m,:)=PPy(2:2:m,:)+dt*(P1dyy(2:2:m,:)-P1uyy(2:2:m,:));

PPy(:,1)=PPy(:,1)-dt*(P1Lyx(:,1)-P1ryx(:,1));

PPy(:,3:2:n+1)=PPy(:,3:2:n+1)-dt*(P1Lyx(:,3:2:n+1)-P1ryx(:,3:2:n+1));
PPy(:,3:2:n+1)=PPy(:,3:2:n+1)-dt*(P1Lyx(:,2:2:n)-P1ryx(:,2:2:n));

PPy(:,2:2:n)=PPy(:,2:2:n)-dt*(P2Lyx(:,2:2:n)-P2ryx(:,2:2:n));
PPy(:,2:2:n)=PPy(:,2:2:n)-dt*(P2Lyx(:,1:2:n-1)-P2ryx(:,1:2:n-1));



% this ends the pressure computation

%% viscous stresses computation

s1=mu*A1.*(h1.^2).*(-(DxX1).*(Dyv1)+(Dxu1).*(DyY1)+(DyX1).*(Dxv1)-(DxY1).*(Dyu1))/dx/dx;
s2=mu*A2.*(h2.^2).*(-(DxX2).*(Dyv2)+(Dxu2).*(DyY2)+(DyX2).*(Dxv2)-(DxY2).*(Dyu2))/dx/dx;

T1=mu*A1.*(h1.^2).*((DxX1).*(Dyu1)+(DyY1).*(Dxv1)-(DyX1).*(Dxu1)-(DxY1).*(Dyv1))/dx/dx;
T2=mu*A2.*(h2.^2).*((DxX2).*(Dyu2)+(DyY2).*(Dxv2)-(DyX2).*(Dxu2)-(DxY2).*(Dyv2))/dx/dx;

V1rxx=0*X;
V1rxx(1:m,1:n)=(s1.*DyY1-T1.*DyX1)/2;

V1Lxx=0*X;
V1Lxx(2:m+1,1:n)=(s1.*DyY1-T1.*DyX1)/2;


V2rxx=0*X;
V2rxx(1:m,1:n)=(s2.*DyY2-T2.*DyX2)/2;

V2Lxx=0*X;
V2Lxx(2:m+1,1:n)=(s2.*DyY2-T2.*DyX2)/2 ;


V1uxy=0*X;
V1uxy(1:m,1:n)=(T1.*DxX1-s1.*DxY1)/2;

V1dxy=0*X;
V1dxy(1:m,2:n+1)=(T1.*DxX1-s1.*DxY1)/2;

V2uxy=0*X;
V2uxy(1:m,1:n)=(T2.*DxX2-s2.*DxY2)/2;

V2dxy=0*X;
V2dxy(1:m,2:n+1)=(T2.*DxX2-s2.*DxY2)/2;



% Viscous Forces in x eqn
VVx=0*X;

VVx(:,1)=VVx(:,1)+(V1rxx(:,1)-V1Lxx(:,1));

VVx(:,3:2:n+1)=VVx(:,3:2:n+1)+(V1rxx(:,3:2:n+1)-V1Lxx(:,3:2:n+1));
VVx(:,3:2:n+1)=VVx(:,3:2:n+1)+(V1rxx(:,2:2:n)-V1Lxx(:,2:2:n));

VVx(:,2:2:n)=VVx(:,2:2:n)+ (V2rxx(:,2:2:n)-V2Lxx(:,2:2:n));
VVx(:,2:2:n)=VVx(:,2:2:n)+ (V2rxx(:,1:2:n-1)-V2Lxx(:,1:2:n-1));

VVx(1,:)=VVx(1,:)+(V2uxy(1,:)-V2dxy(1,:));

VVx(2:2:m,:)=VVx(2:2:m,:)+(V1uxy(2:2:m,:)-V1dxy(2:2:m,:));
VVx(2:2:m,:)=VVx(2:2:m,:)+(V1uxy(1:2:m-1,:)-V1dxy(1:2:m-1,:));

VVx(3:2:m+1,:)=VVx(3:2:m+1,:)+(V2uxy(3:2:m+1,:)-V2dxy(3:2:m+1,:));
VVx(3:2:m+1,:)=VVx(3:2:m+1,:)+(V2uxy(2:2:m,:)-V2dxy(2:2:m,:));



V1ryx=0*X;
V1ryx(1:m,1:n)=(T1.*DyY1+s1.*DyX1)/2;

V1Lyx=0*X;
V1Lyx(2:m+1,1:n)=(T1.*DyY1+s1.*DyX1)/2;

V2ryx=0*X;
V2ryx(1:m,1:n)=(T2.*DyY2+s2.*DyX2)/2;

V2Lyx=0*X;
V2Lyx(2:m+1,1:n)=(T2.*DyY2+s2.*DyX2)/2; 


V1uyy=0*X;
V1uyy(1:m,1:n)=(T1.*DxY1+s1.*DxX1)/2;

V1dyy=0*X;
V1dyy(1:m,2:n+1)=(T1.*DxY1+s1.*DxX1)/2; 

V2uyy=0*Y;
V2uyy(1:m,1:n)=(T2.*DxY2+s2.*DxX2)/2;

V2dyy=0*X;
V2dyy(1:m,2:n+1)=(T2.*DxY2+s2.*DxX2)/2;


% Viscous Forces in y eqn
VVy=0*Y;

VVy(:,1)=VVy(:,1)+(V1ryx(:,1)-V1Lyx(:,1));

VVy(:,3:2:n+1)=VVy(:,3:2:n+1)+(V1ryx(:,3:2:n+1) -V1Lyx(:,3:2:n+1));
VVy(:,3:2:n+1)=VVy(:,3:2:n+1)+(V1ryx(:,2:2:n)-V1Lyx(:,2:2:n));

VVy(:,2:2:n)=VVy(:,2:2:n)+(V2ryx(:,2:2:n)-V2Lyx(:,2:2:n));
VVy(:,2:2:n)=VVy(:,2:2:n)+(V2ryx(:,1:2:n-1)-V2Lyx(:,1:2:n-1));

VVy(1,:)=VVy(1,:)-(V2uyy(1,:)-V2dyy(1,:));

VVy(3:2:m+1,:)=VVy(3:2:m+1,:)-(V2uyy(3:2:m+1,:)-V2dyy(3:2:m+1,:));
VVy(3:2:m+1,:)=VVy(3:2:m+1,:)-(V2uyy(2:2:m,:)-V2dyy(2:2:m,:));

VVy(2:2:m,:)=VVy(2:2:m,:)-(V1uyy(2:2:m,:)-V1dyy(2:2:m,:));
VVy(2:2:m,:)=VVy(2:2:m,:)-(V1uyy(1:2:m-1,:)-V1dyy(1:2:m-1,:));

% this concludes viscous stress computation

%%  Computation of x and y Forces

Fx=PPx+VVx;
Fy=PPy+VVy;


% this concludes the force computation

%%  Velocity update

u=(1-a*dt)*u+(Fx./M-dt*(X-b*v));
v=(1-a*dt)*v+(Fy./M-dt*(Y+b*u));


%% Leap Frog

X=X+dt*u;
Y=Y+dt*v;


%% Store data to pass to next step
%update some data
A1=(DxX1.*DyY1-DyX1.*DxY1)/2;
A2=(DxX2.*DyY2-DyX2.*DxY2)/2;


%  only store data that has been changed
data.X = X; % matrix of x-coordinate values
data.Y = Y; % matrix of y-coordinate values
data.u = u; % matrix of x-components of velocity
data.v = v; % matrix of y-components of velocity

data.A1 = A1;
data.A2 = A2;

% This data is stored for plotting - DEPRECATED FOR NOW
%SS = s1.*s1+T1.*T1+s2.^2 +T2.^2;
%data.SS1_1 = max(max(SS));
%data.SS1_2 = max(max(max(A1)),max(max(A2)))/dx/dx;
%data.SS1_3 = min(min(A2))/dx/dx;

end