% Last modified: 30/12/2012 by Pete

%% Deterministic Microorganism Tracer System Set-up Script

% extract parameters from app structure
truereac = app.truereac ;
reac = app.reac ;
nu = app.nu ;
a = app.a ; b = app.b ;
dx = app.dx ; dy = app.dy ;
dt = app.dt ; T = app.T ;
SaveIntrvl = app.SaveIntrvl ;
GenFig = app.GenFig ;
modelnoise = app.modelnoise ;
truenoise = app.truenoise ;
ObsTimes = app.ObsTimes ;
ObsLocs = app.ObsLocs ;
tag = app.tag ;


%% problem framework set-up

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) ;
XRho = XP ; YRho = YP ;

% define geometry mask with square obstacle
Mask = zeros(Ny,Nx) ;
xc = a/2 - 3.0 ; yc = b/2 ;
Mask( (XP-xc).^2 + (YP-yc).^2 <= (0.5)^2 ) = 1 ;
% Mask( abs(XP-xc) <= 0.5 & abs(YP-yc) <= 0.5 ) = 1 ;
Mask = double( Mask | flipud(Mask) ) ;

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

% specify tracer evolution equation parameters
rhomax = 1.0 ;      % maximum rho
rhog = 0.8 ;        % growth coefficient
rhok = 1.0 ;        % growth regularization factor
rhod = 0.1 ;        % decay coefficient
rhoL = 1.0 ;        % spatial effect length scale
rhoAL = 1.8 ;       % spatial effect amplitude
rhoT = 0.7 ;        % temporal effect time scale
rhoAT = 0.5 ;       % temporal effect amplitude

% define tracer evolution component equations
Fg  = @(rho) ...    % mean growth equation
    rhog*( rho.*(rhomax-rho)./(rhomax-rho+rhok) ) ;
Fgi = @(rho) ...    % mode growth equation
    rhog*( rho.*(rhomax-rho)./(rhomax-rho+rhok) + rhomax - 2*rho ) ;
Fd  = @(rho) ...    % mean decay equation
    rhod*rho ;
Fdi = @(rho) ...    % mode decay equation
    rhod*(rhomax-rho+rhok) ;
FL  = @(r) ...      % spatial effect
    rhoAL*exp(-((r(:,1)-xc).^2+(r(:,2)-yc).^2)/(rhoL^2)) ;
FT  = @(t) ...      % temporal effect
    rhoAT*sin(t/rhoT) ;

% specify spatial and temporal models based on selected reac:
% 
%          o-----------------o
%          |     |     |     |
%        + |  1  |  2  |  3  |
%          |     |     |     |
%          |-----------------|
%          |     |     |     |
%     T  0 |  4  |  5  |  6  |
%          |     |     |     |
%          |-----------------|
%          |     |     |     |
%        - |  7  |  8  |  9  |
%          |     |     |     |
%          o-----------------o
% 
%             -     0     +
% 
%                   L
% 
if ( 1 <= reac && reac <=9 )
    Lmodel = mod(reac-1,3) - 1 ; Tmodel = 2 - ceil(reac/3) ;
else
    Lmodel = -1/2 ; Tmodel = 1/2 ;  % implicit model
end

% define tracer evolution equations
Frho  = @(rho,t,r) (1+Tmodel*FT(t))*Fg(rho) - (1+Lmodel*FL(r)).*Fd(rho) ;
Frhoi = @(rhoi,rho,t,r) ...
    ( (1+Tmodel*FT(t))*Fgi(rho) - (1+Lmodel*FL(r)).*Fdi(rho) ) ...
    .* rhoi ./ (rhomax-rho+rhok) ;

% set boundary conditions
%#ok<*NBRAK>
bcsDP = [ 0 ] ; bcsNP = [ 0 0 0 0 ] ; bcsOP = [ 1 ] ;
bcsDu = [ 0 1 ] ; bcsNu = [ 0 0 0 ] ;
bcsDv = [ 0 0 0 0 ] ; bcsNv = [ 0 ] ;
bcsDrho = [ 0 0.1 ] ;	bcsNrho = [ 0 0 0 ] ;
Nbcs = length(bcsDP) + length(bcsNP) ;
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 ;
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 ;

% initialize P, u, v, and rho with zeros
P = [ bcsDP' ; bcsNP' ; zeros(size(Pid)) ] ;

% initialize u and v with asymmetric perturbation
xcpert = xc - 2.0 ; ycpert = yc - 1.0 ;
falloff = @(X,Y) exp( -5 * sqrt( (X-xcpert).^2 + (Y-ycpert).^2 ) ) ;
uinit = @(X,Y) -50 * (Y - ycpert) .* falloff(X,Y) ;
vinit = @(X,Y) +50 * (X - xcpert) .* falloff(X,Y) ;
u = [ bcsDu' ; bcsNu' ; uinit(XU(idsu),YU(idsu)) ] ;
v = [ bcsDv' ; bcsNv' ; vinit(XV(idsv),YV(idsv)) ] ;

% initialize rho
rhoinit = @(X,Y) 0.1 * ones(size(X)) ;
rho = [ bcsDrho' ; bcsNrho' ; rhoinit(XRho(idsrho),YRho(idsrho)) ] ;


%% 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(uid) ;       % number of u-velocity state variables
n_vvars   = length(vid) ;       % number of v-velocity state variables
n_rhovars = length(rhoid) ;     % number of rho state variables
n_vars    = n_uvars + n_vvars + n_rhovars ;

% one observation per location
n_obs = size(ObsLocs,1) ;

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

% find indexes of observed state variables
for i = 1:n_obs

    % find P indexes
    [~,xind] = min( abs(xp-ObsLocs(i,1)), [], 2 ) ;
    [~,yind] = min( abs(yp-ObsLocs(i,2)), [], 2 ) ;
    H( i , NodeP(yind+1,xind+1) - Nbcs + n_uvars + n_vvars ) = 1 ;
    
end

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

