% Last modified: 06/09/2012 by Pete

%% Stochastic Flow Past an Obstacle Set-up Script

% extract parameters from app structure
S = app.S ; MC = app.MC ;
varscale = app.varscale ;
EnKF = app.EnKF ;
subtag = app.subtag ;


%% problem framework and observation set-up

Pete_detsetup

% remove asymmetric perturbation from u- and v-velocity means
u = [ bcsDu' ; bcsNu' ; zeros(size(uid)) ] ;
v = [ bcsDv' ; bcsNv' ; zeros(size(vid)) ] ;


%% DO set-up

% initialize DO modes
ui = zeros(length(u),S) ;
vi = zeros(length(v),S) ;
Pi = zeros(length(P),S) ;
rhoi = zeros(length(rho),S) ;
var = zeros(S,1) ;

% create fall-off function around obstacle
Mask = Nodeu(2:Ny+1,2:Nx) ;
mids = (Mask == 1) ;
Mask(:) = 1 ;
Mask(mids) = 0 ;
Maskv = Nodev(2:Ny,2:Nx+1) ;
midv = (Maskv == 1) ;
Maskv(:) = 1 ;
Maskv(midv) = 0 ;

% blend mask
for i = 1:2*ceil(Ny)
   Mask(2:end-1,2:end-1) = (1/5) * ( Mask(2:end-1,2:end-1) ...
       + Mask(1:end-2, 2:end-1) + Mask(3:end, 2:end-1) ...
       + Mask(2:end-1,1:end-2) + Mask(2:end-1,3:end) ) ;
   Mask(mids) = 0 ;
   Maskv(2:end-1,2:end-1) = (1/5) * ( Maskv(2:end-1,2:end-1) ...
       + Maskv(1:end-2, 2:end-1) + Maskv(3:end, 2:end-1) ...
       + Maskv(2:end-1,1:end-2) + Maskv(2:end-1,3:end) ) ;
   Maskv(midv) = 0 ;
end

Mask = Mask(idsu) ;
Maskv = Maskv(idsv) ;

% Matt's magic DO initialization...
uinit = @(x,y,M,N) -(1/b)*sin(pi.*x/a).*sin(pi.*M.*x/a) ...
    .*( cos(pi.*y/b).*sin(pi.*N.*y/b) ...
        + N.*sin(pi.*y/b).*cos(pi.*N.*y/b) ) ;
vinit = @(x,y,M,N)  (1/a)*sin(pi.*y/b).*sin(pi.*N.*y/b) ...
    .*( cos(pi.*x/a).*sin(pi.*M.*x/a) ...
        + M.*sin(pi.*x/a).*cos(pi.*M.*x/a) ) ;

i = 1 ;
for ii = 1:S
    for iiy = 1:ceil(ii/2)
        if ( i <= S )
            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(:,i) = [ zeros(Nbcs,1) ; ...
                uinit( XU(idsu), YU(idsu), ix, iy ) / nrm .* Mask ] ;
            vi(:,i) = [ zeros(Nbcs,1) ; ...
                vinit( XV(idsv), YV(idsv), ix, iy ) / nrm .* Maskv ] ;
            var(i) = varscale * exp(-(ix+iy) + 2) ;
            i = i + 1 ;     
        end
        if ( i <= S && 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(:,i) = [ zeros(Nbcs,1) ; ...
                uinit( XU(idsu), YU(idsu), ix, iy ) / nrm .* Mask ] ;
            vi(:,i) = [ zeros(Nbcs,1) ; ...
                vinit( XV(idsv), YV(idsv), ix, iy ) / nrm .* Maskv ] ;
            var(i) = varscale * exp(-(ix+iy) + 2) ;
            i = i + 1 ;   
        end
    end
end

% orthonormalize DO modes
for i = 1:S
    for j = 1:S
        if ( i ~= j )
            coef = InProd(ui(:,i),ui(:,j),dx,dy,uid) ...
                   + InProd(vi(:,i),vi(:,j),dx,dy,vid) ;
            ui(:,i) = ui(:,i) - coef * ui(:,j) ;
            vi(:,i) = vi(:,i) - coef * vi(:,j) ;
        end
    end
    nrm = sqrt(InProd(ui(:,i),ui(:,i),dx,dy,uid) ...
               + InProd(vi(:,i),vi(:,i),dx,dy,vid)) ;
    ui(:,i) = ui(:,i)/nrm ;
    vi(:,i) = vi(:,i)/nrm ;
end

% initialize mode coefficient ensemble particles (draw half from Gaussians,
% then mirror about 0)
YYt = zeros(MC,S) ;
for i = 1:S
    tmp = normrnd( 0, sqrt(var(i)), MC/2, 1 ) ;
    tmp = [ tmp ; -tmp ] ;
    tmp = tmp * sqrt(var(i)) / std(tmp) ;
    YYt(:,i) = tmp ;
end

% enforce symmetry of mode coefficient covariances
for k = 1:4
    for i = 1:S
        for j = 1:S
            if ( i ~= j )
                coef = YYt(:,i)' * YYt(:,j) / (MC-1) ;
                YYt(:,i) = YYt(:,i) - coef * YYt(:,j) / std(YYt(:,j)).^2 ;
            end
        end
        YYt(:,i) = YYt(:,i) * sqrt(var(i)) / std(YYt(:,i)) ;
    end
end

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

