%% Bottom Gravity Current Test Case Setup Script

%-------------------------------------------------
% Set up deterministic portion.
%-------------------------------------------------

BGC_Setup;

app.IP = 1./((app.sigma_nd.^2).*([length(u),length(v),length(rho)]-Nbcs)*dx*dy);


%-------------------------------------------------
% Define body force for DO modes.
%-------------------------------------------------

app.Frhoi = @(t, XY, rhoi) - kappa * (rhoi .*...
    repmat(((XY(:,1) <= a2)&(XY(:,2) >= b1)),1,size(rhoi,2)));

%-------------------------------------------------
% Initialize DO modes.
%-------------------------------------------------

% extract parameters from app structure
S = app.S ; MC = app.MC ;

if isfield (app,'DO_IC_flag')
    DO_IC_flag = app.DO_IC_flag;
else
    DO_IC_flag = 'SineCosine_Plateau';
end;

if ~isfield (app,'MollifierWidth')
    app.MollifierWidth = [];
end;

ui = zeros(length(u),app.S(1)) ;
vi = zeros(length(v),app.S(1)) ;
Pi = zeros(length(P),app.S(1)) ;
rhoi = zeros(length(rho),app.S(1)) ;
Yi_Std = ones(1, app.S(1)) ;

switch DO_IC_flag
    
    case {'SineCosine_Plateau', 'SineCosine_Plateau_Sic'};
        
        % Initialize with Sines & Cosines over plateau or only over
        % plateau region with nonzero initial S.
        
        if strcmp (DO_IC_flag, 'SineCosine_Plateau_Sic')
            a_do_ic = a2;
        else
            a_do_ic = a1;
        end;
        b_do_ic = b-b1;
        
        initfn1 = @(x,y,M,N) (-(1/b_do_ic)*sin(pi.*x/a_do_ic)...
            .*sin(pi.*M.*x/a_do_ic)...
            .*(cos(pi.*(b-y)/b_do_ic).*sin(pi.*N.*(b-y)/b_do_ic)...
            + N.*sin(pi.*(b-y)/b_do_ic).*cos(pi.*N.*(b-y)/b_do_ic)))...
            .* ((x<=a_do_ic) & (y>=b1));
        initfn2 = @(x,y,M,N)  ((1/a_do_ic)*sin(pi.*(b-y)/b_do_ic)...
            .*sin(pi.*N.*(b-y)/b_do_ic) ...
            .*( cos(pi.*x/a_do_ic).*sin(pi.*M.*x/a_do_ic) ...
            + M.*sin(pi.*x/a_do_ic).*cos(pi.*M.*x/a_do_ic) ))...
            .* ((x<=a_do_ic) & (y>=b1));
        
        i = 1 ;
        i_rho = 1;
        for ii = 1:app.S(1)
            for iiy = 1:ceil(ii/2)
                if ( i <= app.S(1) )
                    ix = ii - iiy + 1 ;
                    iy = iiy ;
                    nrm = sqrt(((ix^2 + 1))/16*(3/2)^(iy==1) ...
                        + ((iy^2 + 1))/16*(3/2)^(ix==1)) ;
                    ui(Nbcs+1:length(u),i) = initfn1(XU(idsu), YU(idsu), ix, iy)...
                        / (nrm*sqrt(app.IP(1)));
                    vi(Nbcs+1:length(v),i) = initfn2(XV(idsv), YV(idsv), ix, iy)...
                        / (nrm*sqrt(app.IP(2)));
                    Yi_Std(i) = 10 * exp((1-i)/2.5) ;
                    i = i + 1;
                    if (i_rho <=app.S(1))
                        rhoi(Nbcs+1:length(rho),i_rho) =...
                            initfn1(XP(idsP),YP(idsP),ix,iy)./(nrm.*sqrt(app.IP(3)));
                        i_rho = i_rho + 1;
                    end;
                    if ( i_rho <= app.S(1) )
                        rhoi(Nbcs+1:length(rho),i_rho) =...
                            initfn2(XP(idsP),YP(idsP),ix,iy)./(nrm*sqrt(app.IP(3)));
                        i_rho = i_rho + 1;
                    end;
                end
                if ( i <= app.S(1) && ix ~= iy )
                    % ix and iy interchanged
                    iy = ii - iiy + 1 ;
                    ix = iiy ;
                    nrm = sqrt(((ix^2 + 1))/16*(3/2)^(iy==1) ...
                        + ((iy^2 + 1))/16*(3/2)^(ix==1)) ;
                    ui(Nbcs+1:length(u),i) = initfn1(XU(idsu), YU(idsu), ix, iy)...
                        / (nrm*sqrt(app.IP(1)));
                    vi(Nbcs+1:length(v),i) = initfn2(XV(idsv), YV(idsv), ix, iy)...
                        / (nrm*sqrt(app.IP(2)));
                    Yi_Std(i) = 10 * exp((1-i)/2.5) ;
                    i = i + 1;
                    if (i_rho <=app.S(1))
                        rhoi(Nbcs+1:length(rho),i_rho) =...
                            initfn1(XP(idsP),YP(idsP),ix,iy)./(nrm*sqrt(app.IP(3)));
                        i_rho = i_rho + 1;
                    end;
                    if ( i_rho <= app.S(1) )
                        rhoi(Nbcs+1:length(rho),i_rho) =...
                            initfn2(XP(idsP),YP(idsP),ix,iy)./(nrm*sqrt(app.IP(3)));
                        i_rho = i_rho + 1;
                    end;
                end
            end
        end
        
        % Gram-Schmidt Orthonormalization of DO modes
        
        uid = Nbcs+1:length(u);
        vid = Nbcs+1:length(v);
        pid = Nbcs+1:length(rho);
        
        for i = 1 : S
            for j = 1 : (i-1)
                coef = (app.IP(1)*(ui(uid,i)'*ui(uid,j))...
                    + app.IP(2)*(vi(vid,i)'*vi(vid,j))...
                    + app.IP(3)*(rhoi(pid,i)'*rhoi(pid,j)))*(dx*dy);
                ui(uid,i) = ui(uid,i) - coef * ui(uid,j) ;
                vi(vid,i) = vi(vid,i) - coef * vi(vid,j) ;
                rhoi(pid,i) = rhoi(pid,i) - coef * rhoi(pid,j);
            end;
            nrm = (app.IP(1)*(ui(uid,i)'*ui(uid,i))...
                + app.IP(2)*(vi(vid,i)'*vi(vid,i))...
                + app.IP(3)*(rhoi(pid,i)'*rhoi(pid,i)))*(dx*dy);
            ui(uid,i) = ui(uid,i)/nrm ;
            vi(vid,i) = vi(vid,i)/nrm ;
            rhoi(pid,i) = rhoi(pid,i)/nrm;
        end;
        
        % Initialize stochastic coefficients as normal random with variance
        % as calculatedin Stoch_Initial_Perturb
        Yi_Std = Yi_Std*(app.varscale/sqrt(sum(Yi_Std.^2)));
        Yi = repmat(Yi_Std, app.MC, 1).*normrnd(0,1,app.MC,app.S(1));
        
    case {'SVD_ICs'};
        
        % This option creates modes by SVD of frequent output from short run.
        
        [u,v,rho,ui,vi,rhoi,Yi] = BGC_DO_simsvdIC(app);
        
    case {'Correlation_Rho', 'Correlation_Rho_Sic' };
        
        switch (DO_IC_flag)
            case 'Correlation_Rho_Sic'
                a_limit = a2;
            otherwise
                a_limit = a1;
        end;
        
        [rhoi, Yi] = Stoch_Initial_Perturb_rho...
            (XP,YP,dx,dy,a_limit,b,b1,app,Nbcs,NodeRho,idsP);
        
    case {'Restart'};
        
        [u,v,rho,P,ui,vi,rhoi,Pi,Yi] = BGC_DO_restart...
            (app,u,v,rho,P,ui,vi,rhoi,Pi,idsu,idsv,idsP,Nodeu,Nodev,NodeP);
        
end;

%% Observation Parameters

if ~isfield (app,'truenoise')
    app.truenoise = 0.01;   % Actual standard deviation of observation noise
end;
if ~isfield (app,'modelnoise')
    app.modelnoise = 0.05;  % Modeled standard deviation of observation noise
end;

if ~isfield (app,'EnKF')
    app.EnKF = false;       % 1: Ensemble Kalman filter; 0: GMM-DO filter
end;

% Times of calling observation script (round to nearest saving time)
if isfield (app,'ObsTimes')
    ObsTimes = app.ObsTimes;
else
    ObsTimes = (30:1:120).*(60*nu_x/h^2);
end;

% round observation times down to nearest saving times
ObsTimes = round(ObsTimes/(app.PlotIntrvl*app.dt));

n_obssteps = length(ObsTimes);
obsstep = 1;

% Spatial coordinates (x,y) at which to make observations
if isfield (app,'ObsLocs')
    ObsLocs = app.ObsLocs;
else
    ObsLocs = [...
        2000, (A1-2000)*tan_theta + B1 + 25;
        2000, (A1-2000)*tan_theta + B1 + 75;
        2000, (A1-2000)*tan_theta + B1 + 125;
        3000, (A1-3000)*tan_theta + B1 + 30;
        3000, (A1-3000)*tan_theta + B1 + 90;
        3000, (A1-3000)*tan_theta + B1 + 150;
        4000, (A1-4000)*tan_theta + B1 + 30;
        4000, (A1-4000)*tan_theta + B1 + 90;
        4000, (A1-4000)*tan_theta + B1 + 150]/h;
end;

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

% One observation per location: salinity
n_obs = size(ObsLocs,1);

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

% Find indexes of observed state variables
for i = 1 : n_obs
    [~,xind] = min(abs(XP(1,:) - ObsLocs(i,1)));
    [~,yind] = min(abs(YP(:,1) - ObsLocs(i,2)));
    H(i, n_uvars + n_vvars + NodeRho(yind+1,xind+1)) = 1;
end;

% Create observation noise covariance matrix
R = (app.modelnoise^2)*eye(n_obs);

% Minimum and maximum number of components of GMM
if ~isfield(app,'M_Min');
    app.M_Min = 5;
end;
if ~isfield(app,'M_Max');
    app.M_Max = 5;
end;

if ~isfield(app, 'GMM_Fit_Cov')
    app.GMM_Fit_Cov = 'Full_NotShared';
end;

if ~isfield(app,'ObsScript');
    app.ObsScript = 'BGC_DO_GMM_Filter';
end;

if ~isfield(app,'ObsDataDir');
    app.ObsDataDir = 'Save/mat/BGC_Det_True';
end;

% initialize array of observation log-likelihoods
loglikes = double.empty(0,2) ;

