% Last modified: 30/12/2012 by Pete

%% Deterministic Flow Past an Obstacle Set-up Script

% extract parameters from app structure
trueshape = app.trueshape ;
shape = app.shape ;
nu = app.nu ;
a = app.a ; b = app.b ;
dx = app.dx ; dy = app.dy ;
dt = app.dt ; T = app.T ;
SaveIntrvl = app.PlotIntrvl ;
GenFig = app.GenFig ;
modelnoise = app.modelnoise ;
truenoise = app.truenoise ;
ObsTimes = app.ObsTimes ;
ObsLocs = app.ObsLocs ;
tag = app.tag ;
trunk = app.trunk;


%% problem framework set-up

Re = 1/nu ;             % Reynolds number
Nx = a/dx ;             % number of discrete elements in x-direction
Ny = b/dy ;             % number of discrete elements in y-direction
Nt = ceil(T/dt) ;       % number of time steps
app.Nx = Nx ; app.Ny = Ny ;

% define grid coordinates
xp = linspace( dx/2, a-dx/2, Nx ) ;
yp = fliplr( linspace( dy/2, b-dy/2, Ny ) ) ;
xu = linspace( dx, a-dx, Nx-1 ) ;
yu = fliplr( linspace( dy/2, b-dy/2, Ny ) ) ;
xv = linspace( dx/2, a-dx/2, Nx ) ;
yv = fliplr( linspace( dy, b-dy, Ny-1 ) ) ;
[XP,YP] = meshgrid(xp,yp) ;
[XU,YU] = meshgrid(xu,yu) ;
[XV,YV] = meshgrid(xv,yv) ;

% define geometry mask
Mask = zeros(Ny,Nx) ;
xc = a/2 - 3.0 ; yc = b/2 ;
switch shape
    
    % circle
    case 1
        Mask( (XP-xc).^2 + (YP-yc).^2 <= (0.5)^2 ) = 1 ;
        
    % square
    case 2
    		%Mask( abs(XP-xc) <= 0.5 & abs(YP-yc) <= 0.5 ) = 1 ;
    		%CM Hack for Idealized Elizabeth Islands
    	  % Cuttyhunk
        Mask( abs(XP-31/6) <= 2/3 & abs(YP-yc) <= 0.5 ) = 1 ;
        % Nashawena
        Mask( abs(XP-25/3) <= 1.0 & abs(YP-yc) <= 0.5 ) = 1 ;
        % Pasque
        Mask( abs(XP-11) <= 2/3 & abs(YP-yc) <= 0.5 ) = 1 ;
        
    % upstream-pointing triangle
    case 3
        Mask( XP-xc-0.5 <= 0 ...
              & (YP-yc) - 0.5*(XP-xc+0.5) <= 0 ...
              & (YP-yc) + 0.5*(XP-xc+0.5) >= 0 ) = 1 ;
        
    % downstream-pointing triangle
    case 4
        Mask( XP-xc+0.5 >= 0 ...
              & (YP-yc) + 0.5*(XP-xc-0.5) <= 0 ...
              & (YP-yc) - 0.5*(XP-xc-0.5) >= 0 ) = 1 ;
        
    % diamond
    case 5
        Mask( (YP-yc) - (XP-xc+0.5) <= 0 ...
              & (YP-yc) + (XP-xc+0.5) >= 0 ...
              & (YP-yc) + (XP-xc-0.5) <= 0 ...
              & (YP-yc) - (XP-xc-0.5) >= 0 ) = 1 ;
        
end
Mask = double( Mask | flipud(Mask) ) ;

% define forcing functions
Fu = @(t,r) 0 ; Fv = @(t,r) 0 ;

% set boundary conditions
%#ok<*NBRAK>
bcsDq = [ 0 ] ; bcsNq = [ 0 0 0 0 ]' ; bcsOq = [ 1 ] ;
bcsDu = [ 0 1 ]' ; bcsNu = [ 0 0 0 ]' ;
bcsDv = [ 0 0 0 0 ]' ; bcsNv = [ 0 ] ;
bcsDrho = [ 0 0 0 0 ]' ;	bcsNrho = [ 0 ] ;
Nbcs = length(bcsDq) + length(bcsNq) ;
bnd = [ 1 1 1 1 ] ;

% create node matrices
[ NodeP, Nodeu, Nodev, idsP, idsu, idsv ] = NodePad(Mask,bnd,1) ;
NodeRho = NodeP ; idsRho = idsP ;
NodeP(NodeP == 1) = -1 ; NodeP(NodeP == 3) = 1 ; NodeP(NodeP == -1) = 3 ;
Nodev(Nodev == 3) = -1 ; Nodev(Nodev == 5) = 3 ; Nodev(Nodev == -1) = 5 ;
NodeRho(NodeRho == 3) = -1 ; NodeRho(NodeRho == 5) = 3 ; ...
    NodeRho(NodeRho == -1) = 5 ;
pid = (1:length(idsP))' + Nbcs*ones(size(idsP)) ;
uid = (1:length(idsu))' + Nbcs*ones(size(idsu)) ;
vid = (1:length(idsv))' + Nbcs*ones(size(idsv)) ;
rhoid = pid ;

app.RmSing = 1;

% initialize P, u, and v with zeros
P = [ bcsDq ; bcsNq ; zeros(size(pid)) ] ;
u = [ bcsDu ; bcsNu ; zeros(size(uid)) ] ;
v = [ bcsDv ; bcsNv ; zeros(size(vid)) ] ;

% initialize rho with asymmetric perturbation
perturb = 1 ;
rhoinit = @(X,Y) double( abs( X-(0.25+abs(perturb)) ) <= 0.25 ...
    & abs( Y-(yc+0.5*sign(perturb)) ) <= 0.25 ) * sign(perturb) ;
rho = [ bcsDrho ; bcsNrho ; rhoinit(XP(idsP),YP(idsP)) ] ;


%% observation set-up

% round observation times down to nearest saving times
ObsTimes = ObsTimes - mod( ObsTimes , SaveIntrvl*dt ) ;
n_obssteps = size(ObsTimes,2) ;
obsstep = 1 ;

n_uvars = length(u) ;           % number of u-velocity state variables
n_vvars = length(v) ;           % number of v-velocity state variables
n_vars  = n_uvars + n_vvars ;	% total number of state variables

% two observations per location, one of u-velocity and one of v-velocity
n_obs = 2*size(ObsLocs,1) ;

% create observation matrix
H = zeros(n_obs,n_vars) ;

% find indexes of observed state variables (the u-velocity grid points
% immediately to the North and the v-velocity grid points immediately
% to the East of the sampling locations are observed)
for i = 1:n_obs/2
    
    % find u-velocity indexes
    [~,xind] = min( abs(xu-ObsLocs(i,1)), [], 2 ) ;
    [~,yind] = min( abs(xu-(ObsLocs(i,2)+dy/2)), [], 2 ) ;
    H( 2*i - 1 , Nodeu(yind+1,xind+1) ) = 1 ;
    
    % find v-velocity indexes
    [~,xind] = min( abs(xv-(ObsLocs(i,1)+dx/2)), [], 2 ) ;
    [~,yind] = min( abs(yv-ObsLocs(i,2)), [], 2 ) ;
    H( 2*i , Nodev(yind+1,xind+1) + n_uvars ) = 1 ;
    
end

% create observation noise matrix
R = (modelnoise^2)*eye(n_obs) ;

