%% Bottom Gravity Current Test Case Setup Script

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

Bottom_Gravity_Current_Setup

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

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

%-------------------------------------------------
% Set scalings for DO inner products.
%-------------------------------------------------

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

app.IP = 1./(dx*dy*(app.sigma_nd.*app.sigma_nd));

 % want norm to be 1, i.e. (dx*dy*var_nd_rho*rhoi*rhoi)=1
 % rescale = sqrt(app.IP(3))

%-------------------------------------------------
% 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)) ;
var = ones(app.S(1),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. (Not tested with revises app.IP)

      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)));
                  var(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)));
                  var(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

      % orthonormalize DO modes

      uid = Nbcs+1:length(u);
      vid = Nbcs+1:length(v);
      pid = Nbcs+1:length(rho);

      for i = 1:S
          for j = 1:S
              if ( i ~= j )
                  coef = app.IP(1) * InProd(ui(:,i),ui(:,j),dx,dy,uid) ...
                         + app.IP(2) * InProd(vi(:,i),vi(:,j),dx,dy,vid) ...
                         + app.IP(3) * InProd(rhoi(:,i), rhoi(:,i), dx, dy, pid);
                  ui(:,i) = ui(:,i) - coef * ui(:,j) ;
                  vi(:,i) = vi(:,i) - coef * vi(:,j) ;
                  rhoi(:,i) = rhoi(:,i) - coef * rhoi(:,j);
              end
          end
          nrm = sqrt(app.IP(1)*InProd(ui(:,i),ui(:,i),dx,dy,uid) ...
                     + app.IP(2)*InProd(vi(:,i),vi(:,i),dx,dy,vid) ...
                     + app.IP(3)*InProd(rhoi(:,i), rhoi(:,i), dx, dy, pid)) ;
          ui(:,i) = ui(:,i)/nrm ;
          vi(:,i) = vi(:,i)/nrm ;
          rhoi(:,i) = rhoi(:,i) / nrm;
      end

      % Redimensionalize modes

      ui = ui .*app.varscale.*app.sigma_nd(1);
      vi = vi .*app.varscale.*app.sigma_nd(2);
      rhoi = rhoi .*app.varscale.*app.sigma_nd(3);

      %Initialize stochastic coefficients as normal random with variance as calculated
      %in Stoch_Initial_Perturb
      Yi=zeros(app.MC,app.S(1));
      for i=1:app.S(1)
          Yi(:,i)=var(i).*normrnd(0,1,app.MC,1);
      end

   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;

%keyboard
      [rhoi, Yi] = Stoch_Initial_Perturb_rho (XP,YP,dx,dy,a_limit,b,b1,app,Nbcs,NodeP,idsP);

   case {'Restart'};

   % This option creates ICs from output of previous DO simulation.
   % (even from previous simulation in truncated domain).

   [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;

