function [AVGuonv AVGvonu] = CoriolisAVG(Nx, Ny, Nodeu, Nodev, idsu, idsv, Nbcs)
%function [AVGuonv AVGvonu] = CoriolisAVG(Nx, Ny, Nodeu, Nodev, idsu, idsv, Nbcs)
%This function creates the averaging matrixes to compute the value of the 
% coriolis terms on the staggered u and v C-grids. 
%
% INPUTS:
%           Nx:         Number of interior nodes in x-direction (integer
%                        number)
%           Ny:         Number of interior nodes in y-direction (integer
%                        number)
%           Nodeu:  	U-velocity node-numbering matrix (see NodePad.m)
%           Nodev:  	V-velocity node-numbering matrix (see NodePad.m)
%           idsu:   	Interior U-velocity node id's (see NodePad.m)
%           idsv:   	Interior V-velocity node id's (see NodePad.m)
%           Nbcs:	Number of boundary conditions
%
% OUTPUTS:
%           AVGuonv:	Matrix that averages the u velocity on the v-grid
%           AVGvonu:	Matrix that averages the v velocity on the u-grid
%
% Authors: Themis Sapsis, Matt Ueckermann
%%
idsy = 2:Ny;
idsx = 1:Nx;
LU   = Nodeu(idsy,idsx);
RU   = Nodeu(idsy,idsx+1);
LL   = Nodeu(idsy+1,idsx);
RL   = Nodeu(idsy+1,idsx+1);

LU   = LU(idsv);
RU   = RU(idsv);
LL   = LL(idsv);
RL   = RL(idsv);

AVGuonv = spalloc(length(idsv), length(idsu) + Nbcs, 4 * length(idsv) + 5); %this is on the v grid
for i=1:length(idsv)
   AVGuonv(i,LL(i)) = 0.25;
   AVGuonv(i,LU(i)) = AVGuonv(i,LU(i))+0.25;
   AVGuonv(i,RU(i)) = AVGuonv(i,RU(i))+0.25;
   AVGuonv(i,RL(i)) = AVGuonv(i,RL(i))+0.25;
end

clear LU RU LL RL

idsy = 1:Ny;
idsx = 2:Nx;

LU   = Nodev(idsy,idsx);
RU   = Nodev(idsy,idsx+1);
LL   = Nodev(idsy+1,idsx);
RL   = Nodev(idsy+1,idsx+1);

LU   = LU(idsu);
RU   = RU(idsu);
LL   = LL(idsu);
RL   = RL(idsu);

AVGvonu = spalloc(length(idsu), length(idsv) + Nbcs, 4 * length(idsu) + 5); %this is on the u grid
for i=1:length(idsu)
   AVGvonu(i,LL(i)) = 0.25;
   AVGvonu(i,LU(i)) = AVGvonu(i,LU(i))+0.25;
   AVGvonu(i,RU(i)) = AVGvonu(i,RU(i))+0.25;
   AVGvonu(i,RL(i)) = AVGvonu(i,RL(i))+0.25;
end