% Adaptation for Adaptive DO solvers
% This is basically code that is common amongst the Adaptive DO solvers.

%At first timestep do preliminaries

if k==0
    if isfield(app,'addvar')
       ADDMODEVAR = app.addvar(1); 
    end
    if isfield(app,'subvar')
       SUBMODEVAR = app.subvar(1); 
    end
    if ~exist('ADDMODEVAR','var')
        %Variance of smallest mode that would require a new mode
	ADDMODEVAR = 1e-4; 
    end
    if ~exist('SUBMODEVAR','var')
        %Variance of smallest mode that would remove that mode
        %Will also require second smallest mode to be less
        %than twice this value in adaptation criterion
        SUBMODEVAR = 1e-8; 
    end

    if length(app.S)>1 %Adaptation!
        Smin=app.S(1);
        Smax=app.S(2);
        app.S=Smin;
    else %No adaptation
        Smin=app.S;
        Smax=app.S;
    end    
end

% First Step is to diagonalize the modes
CYY=cov(Yi);
clear DC
[VC,DC] = eig(CYY);  %Diagonal Covariance, and Diagonal Vectors
DC=diag(DC);
%VC = diag(sign(sum(VC)))*VC;
ui(uid,:)=fliplr(ui(uid,:)*(VC));                           % IMPORTANT CHANGE HERE
vi(vid,:)=fliplr(vi(vid,:)*(VC));
Pi(pid,:)=fliplr(Pi(pid,:)*(VC));
Yi=fliplr(Yi*(VC));   % the fliplr ensures that the first mode is the most important

%Next Step is to adapt if necessary
if ~isempty(DC)
    if isfield(app,'tellmeS')
        app.S(1)
    end
    if (DC(1)<SUBMODEVAR && app.S(1)>Smin && DC(2)<2*SUBMODEVAR)
        
         app.S(1)=app.S(1)-1;
         ui(:,app.S(1)+1)=[];
         vi(:,app.S(1)+1)=[];
         Pi(:,app.S(1)+1)=[];
         Yi(:,app.S(1)+1)=[];
         fui = zeros(length(uid),app.S(1));
         fvi = zeros(length(vid),app.S(1));
         fuij= zeros(length(uid),app.S(1)^2);
         fvij= zeros(length(vid),app.S(1)^2);
         YiYj=[];

	 %clear temporary of helper variables so that they won't be the wrong 
         %size next time through the integration loop.
         clear fui_diffusion fvi_diffusion Yi0 A1j A2j RHSY fuijmyyy fvijmyyy
         clear f1u f2u f1v f2v InProdi YiYj CYY VC
%         
%         % Don't forget here also to take into account that app.S=Smax
    elseif (min(DC)>ADDMODEVAR && app.S(1)<Smax)
        tmp=app.S(1);
        [ui,vi, S_temp] = AddMostUnstableMode(...
            u,v,ui,vi,...
            app.Nx,app.Ny,dx,dy,...
            Nodeu,Nodev,idsu,idsv,...
             -Dfs_u,-Dfs_v,uid,vid,existbasis);
        %A new mode may not have been added, so we need to check this:
        if S_temp>tmp
            app.S(1)=S_temp;
            Pi(:,app.S(1))=zeros(size(Pi,1),app.S(1)-tmp);
            Yi(:,app.S(1))=normrnd(0,sqrt(ADDMODEVAR*2/3*10^-1),size(Yi,1),1);
        end
    end
end
