function [ust,vst S] = AddMostUnstableMode_mex(u,v,ust,vst,Nx,Ny,dx,dy,Nodeu,Nodev,idsu,idsv,AU,AV,uid,vid,Nbcs,NbcsDu,NbcsDv,existbasis)
%function [ust,vst S] = AddMostUnstableMode_mex
%                  (u,v,ust,vst,Nx,Ny,dx,dy,Nodeu,Nodev,idsu,idsv,AU,AV,uid,vid)
%
%This function is used to add a new stochastic mode to an adaptive DO simulation. 
% it searches through the saved basis function and finds the most unstable mode
% then adds it to the stochastic modes. This occurs only if the eigenvalue for 
% this mode is larger than 0. This function uses some mexed routines, and this 
% is the only difference to AddMostUnstableMode.m
%
% INPUTS:
%	u: 	Mean u-velocity
%	v:	Mean v-velocity
%	ust:    Stochastic u-velocity modes
%	vst: 	Stochastic v-velocity modes
%       Nx:     Number of interior nodes in the x-direction
%       Ny:     Number of interior nodes in the y-direction
%       dx:     Step size in x direction
%       dy:     Step size in y direction
%       Nodeu:  Node numbering matrix for the u velocity
%       Nodev:  Node numbering matrix for the v velocity
%       idsu:   Interior U-velocity node id's (see NodePad.m)
%       idsv:   Interior V-velocity node id's (see NodePad.m)
%	AU:	Linear operator in u-velocity momentum equation
%	AV:	Linear operator in v-velocity momentum equation
%	uid:    I.d.'s of interior u-velocity nodes [u(uid)]
%       vid:    I.d.'s of interior v-velocity nodes [v(vid)]
%       Nbcs:   Number of boundary conditions
%       NbcsDu: Number of Dirichlet u-velocity boundary conditions
%       NbcsDv: Number of Dirichlet v-velocity boundary conditions
%
% OUTPUTS:
%	ust:    Stochastic u-velocity modes with (potentially) new mode
%	ust:    Stochastic v-velocity modes with (potentially) new mode
%       S:      New stochastic dimension. S = S_in or S = S_in + 1
%
% Authors: Themis Sapsis and Matt Ueckermann.

%Notes
%vst and ust are ui and vi from other frunction, and do NOT include
%boundaries! -- only interior points

%This HAS to be a function in the main code, otherwise code's ui's will be
%destroyed!

%Load the orthonormal basis set from which the new most unstable mode will be 
% constructed. After this new variables ui's and vi's will be present. 
% This basis is originally created and saved either in the setup file by the 
% user, or in the solver file. 
load(sprintf('orth_basis%0.3d.mat',existbasis))
%Find the size of this basis
Mdu=size(ui,2);

%Get the stochastic dimension of the input (S_in)
S=size(ust,2);

%% Construction of Ortonormal basis that approximates Vs 
cni=[];
for i=1:Mdu
    %Project out current stochastic basis from saved basis

    coef =dx*dy*(ust(uid,:)'*ui(uid,i)+vst(vid,:)'*vi(vid,i));
    ui(:,i) = ui(:,i)-ust*coef;
    vi(:,i) = vi(:,i)-vst*coef;
   
    nrm=sqrt(  ( ui(uid,i)'*ui(uid,i)+vi(vid,i)'*vi(vid,i) )*dy*dx  );
% No longer needed, normalize later anyway
%     ui(:,i) = ui(:,i)/nrm;
%     vi(:,i) = vi(:,i)/nrm;
    
    %Have a mode that is very close to existing element in Vs, then
    %if I subtract projection of vs, will get numerical garbage
    %therefore, only want to keep a mode if it has a large normal.
    %This is a numerical condition

    if nrm<=0.01

        cni=[cni i];        

    end
end
%Remove modes from basis that have a small magnitude
ui(:,cni)=[];
vi(:,cni)=[];

%Get new size of Vs
Mdu=size(ui,2);

QQ=zeros(Mdu,Mdu);


%% Computation of Frechet Derivative
for i=1:Mdu    
    % RHS of Navier Stokes - Calculation
    [fui fvi] = UVadvect_DO_mex(    Nx, Ny, dx, dy, ui(:,i), vi(:,i),u      ,v       ,...
        Nodeu, Nodev, idsu, idsv, Nbcs, NbcsDu, NbcsDv);

    [tmp tmp1]= UVadvect_DO_CDS_mex(Nx, Ny, dx, dy, u      , v      ,ui(:,i), vi(:,i),...
        Nodeu, Nodev, idsu, idsv, Nbcs, NbcsDu, NbcsDv);


    fui = fui + tmp + AU*ui(uid,i);
    fvi = fvi + tmp1+ AV*vi(vid,i);

    % Frechet Derivative
    QQ(i,:)=dx*dy*(ui(uid,:)'*fui+vi(vid,:)'*fvi)';   
end


%% Computation of maximum eigenvalue and corresponding direction
[V,D] = eig((QQ+QQ')/2);
%If the maximum eigenvalue is non-negative, we will, in fact, add a new 
% stochastic mode
if max(max(D))>0
    S=S+1;
    [C,I1]=max(diag(D));
   

    ust(:,S) = ui*V(:,I1);
    vst(:,S) = vi*V(:,I1);

    nrm        = sqrt(  ( ust(:,S)'*ust(:,S)+vst(:,S)'*vst(:,S) )*dy*dx  );
    ust(:,S) = ust(:,S)/nrm;
    vst(:,S) = vst(:,S)/nrm;
end

end
