%This is version 3 with changes to decidedMaxSTDv_base_v2 are now below.
%The changes consisted of:
% 0. In getDiurnal, we now send the entire profile, including below the
% emergence height. We also set the PBLH of the profile at the emergence
% time to the emergence height
% 1. Compare each max to min below and above. If the lowest max is the
% first point in the profile, only compare to the min above. For the
% highest max, only compare to the min below unless there is a min above 
% highest max.
% 2. Set startind to 1 (instead of 2) unless the lowest max is the first
% lowest point in the profile. In that case, set it to 2.




function[PBLH,absFactor,negIndicator,not30Ind] = runAlgorithm(sigVector,heightVector,noninterpsig,hookstat,absFactorMin,absFactorMax)

PBLH = NaN;
mini = NaN.*ones(size(heightVector,2),1); 
maxi = NaN.*ones(size(heightVector,2),1);
negIndicator=0;
not30Ind=NaN;
% variablize the submitted vector as Asignal
Asignal = sigVector; 
%%disp('Asignal size');
%%disp(size(Asignal));

% differentiate the vector
Bsigdiff = diff(Asignal);
% for the length of the vector minus one, if one element is positive and
% the next element is negative we have a local maximum. The variable maxi
% contains the value of the signal at the maxima 

for i = 1: length(Bsigdiff) - 1
  if Bsigdiff(i) > 0 && Bsigdiff(i + 1) < 0
   maxi(i + 1) = Asignal(i + 1);
  end
  if Bsigdiff(i) > 0 && Bsigdiff(i + 1) == 0
%       %To deal with cases where there are several points in a row in
%       %the profile that have the same SNR value and then the profile
%       %continues to increase in SNR value. This loop prevents us from
%       %counting those plateaus as maxima.
      for j = i+2:length(Bsigdiff)-1
          if Bsigdiff(j)>0
              break
          end
          if Bsigdiff(j)<0
            maxi(i + 1) = Asignal(i + 1);
            break
          end
      end
  end
end
%%disp('maxi');
%%disp(size(maxi));

% for the length of the vector minus one, if one element is positive and
% the next element is negative we have a local maximum. The variable mini
% contains the value of the signal at the minima 

for i = 1: length(Bsigdiff) - 1
  if Bsigdiff(i) < 0 && Bsigdiff(i + 1) > 0
   mini(i + 1) = Asignal(i + 1);
  end
  if Bsigdiff(i) < 0 && Bsigdiff(i + 1) == 0
%       %To deal with cases where there are several points in a row in
%       %the profile that have the same SNR value and then the profile
%       %continues to decrease in SNR value. This loop prevents us from
%       %counting those plateaus as minima.
      for j = i+2:length(Bsigdiff)-1
          if Bsigdiff(j)<0
              break
          end
          if Bsigdiff(j)>0
            mini(i + 1) = Asignal(i + 1);
            break
          end
      end
  end
end
%%disp('mini');
%%disp(size(mini));
% MaxIndex is a vector of maxima indices only
MaxIndex = find(~isnan(maxi));
%%disp('MaxIndex');
%%disp(MaxIndex);



% maxima are the actual maximum values
maxima = Asignal(MaxIndex);


% MinIndex is a vector of minima indices only
MinIndex = find(~isnan(mini));
%%disp('MinIndex');
%%disp(MinIndex);


% minima are the actual minimum values
minima = Asignal(MinIndex);





%Find maxindex of noninterpsig
mini_ni = NaN.*ones(size(noninterpsig,2),1); 
maxi_ni = NaN.*ones(size(noninterpsig,2),1); 
% variablize the submitted vector as Asignal_ni
Asignal_ni = noninterpsig; 

%Experiment
Asignal_ni=Asignal_ni(~isnan(Asignal_ni));

% differentiate the vector
Bsigdiff_ni = diff(Asignal_ni);
% for the length of the vector minus one, if one element is positive and
% the next element is negative we have a local maximum. The variable maxi
% contains the value of the signal at the maxima 

for i = 1: length(Bsigdiff_ni) - 1
  if Bsigdiff_ni(i) > 0 && Bsigdiff_ni(i + 1) < 0
   maxi_ni(i + 1) = Asignal_ni(i + 1);
  end
  if Bsigdiff_ni(i) > 0 && Bsigdiff_ni(i + 1) == 0
%       %To deal with cases where there are several points in a row in
%       %the profile that have the same SNR value and then the profile
%       %continues to increase in SNR value. This loop prevents us from
%       %counting those plateaus as maxima.
      for j = i+2:length(Bsigdiff_ni)-1
          if Bsigdiff_ni(j)>0
              break
          end
          if Bsigdiff_ni(j)<0
            maxi_ni(i + 1) = Asignal_ni(i + 1);
            break
          end
      end
  end
end


% for the length of the vector minus one, if one element is positive and
% the next element is negative we have a local maximum. The variable mini
% contains the value of the signal at the minima 

for i = 1: length(Bsigdiff_ni) - 1
  if Bsigdiff_ni(i) < 0 && Bsigdiff_ni(i + 1) > 0
   mini_ni(i + 1) = Asignal_ni(i + 1);
  end
  if Bsigdiff_ni(i) < 0 && Bsigdiff_ni(i + 1) == 0
%       %To deal with cases where there are several points in a row in
%       %the profile that have the same SNR value and then the profile
%       %continues to decrease in SNR value. This loop prevents us from
%       %counting those plateaus as minima.
      for j = i+2:length(Bsigdiff_ni)-1
          if Bsigdiff_ni(j)<0
              break
          end
          if Bsigdiff_ni(j)>0
            mini_ni(i + 1) = Asignal_ni(i + 1);
            break
          end
      end
  end
end

% MaxIndex is a vector of maxima indices only
MaxIndex_ni = find(~isnan(maxi_ni));
%will use for setting start index
%%disp('MaxIndex_ni');
%%disp(MaxIndex_ni);
    


% if there are no minima (MinIndex), but there is one maximum (MaxIndex), then pick
% that maximum as the PBL 

if length(MinIndex)==0 & length(MaxIndex)==1 
   PBLH = heightVector(MaxIndex(1));
   absFactor=-99;
   for i=1:MaxIndex
       if Asignal(i)<0
           negIndicator=1;
           break
       end
   end
end

if length(MinIndex)==0 & length(MaxIndex)==0
   absFactor=-99;
end

if length(MinIndex)==1 & length(MaxIndex)==0
   absFactor=-99;
end

if length(noninterpsig(~isnan(noninterpsig)))<2
    absFactor=-99;
end



% if the number of elements in each set are different then one of two
% things must happen. If there are fewer minima than maxima and the
% first element of B is negative then the 500 level signal value is
% considered a maximum




if length(MinIndex) ~= 0 && length(MaxIndex) ~= 0 && length(noninterpsig)> 2

% NOTE: we assume later that the bottom max is below the bottom min
% if this is not the case (ie., bottom level is local max), add an extra
% max at the bottom

  if MinIndex(1) < MaxIndex(1) && Bsigdiff(1) <= 0
%   if Bsigdiff(1) < 0
    MaxIndex = [1 MaxIndex']';
    maxima = [Asignal(1) maxima];
   end
%    disp('maxima');
%    disp(size(maxima));
%    disp(maxima);
%    disp('MaxIndex');
%    disp(size(MaxIndex));
%    disp(MaxIndex);
%    disp('height vector of maxindex');
%    disp(heightVector(MaxIndex));

   if max(MinIndex) < max(MaxIndex)
    addextramin = sum(~isnan(Asignal));
    MinIndex = [MinIndex' addextramin]';
    minima = [minima Asignal(addextramin)];
   end
   

%    disp('MinIndex');
%    disp(MinIndex);
%    disp('height vector of minindex');
%    disp(heightVector(MinIndex));
%    disp('minima');
%    disp(size(minima));
%    disp(minima);

%%% do not need this check - we already are inside an if-loop where it is true:    if length(MaxIndex) ~= 0

   
   %disp("Abs max index (numlow)");
   %disp(numlow);
      %Finding the abs max may have many indices because of the
   %interpolation. taking the middle index as numlow.
   numlow=find(Asignal==max(maxima));
   if length(numlow)>1
    numlow=floor(nanmean(numlow));
   end
   
   
   for i=1:numlow
       if Asignal(i)<0
           negIndicator=1;
           break
       end
   end
   
% Note - this is a tunable parameter, ie., the range of levels to use in getting the stdv
%   (it does not have to be related to the level of the column absolute max)
   
   %Checking if numlow is around 3 doesn't make sense because of the
   %interpolated profiles. This is from when profiles had 15 points. Now,
   %they often have 3000+ points. 3/15 = 0.2. Instead check if numlow is
   %less than 20% of the size of Asignal, and if it is not, set it to the
   %index that is 20% above the first entry in Asignal
   if numlow < 0.2*(length(Asignal))
    numlow = ceil(0.2*(length(Asignal)));
   end
   


   absFactor = nanstd(Asignal(1:numlow),1);
   %%disp('Initial stdv factor');
   %%disp(absFactor);


% The method employed in this study is to utilize some sort of threshold 
% value that can be multiplied by the lowest local maximum in the profile. 
% The result will be subtracted from this maximum and compared to the 
% nearest minimum signal return value located above the present level. 
% A criterion must be met. The result of the subtraction must be larger 
% than the minimum found above the local maximum in question. If this 
% criterion is met, then we can stop searching our local maximum values 
% and settle on this one as representative of the top of the PBL. But if 
% it is not met, then we will continue testing each local maximum higher 
% up in the profile until the criterion is met. If none of the maximum 
% values tested meets this criterion then that hour's PBL height is 
% undetermined. 


   totalmax = nanmax(maxima);
   stdevcap = totalmax .* absFactorMax/100;
   stdevfloor = totalmax .* absFactorMin/100 ;

% if standard deviation is larger than the cap, then
% the cap becomes the new standard deviation
   if absFactor > stdevcap
    absFactor = stdevcap;
   end

% if standard deviation is less than the floor, then
% the floor becomes the new standard deviation
   if absFactor < stdevfloor
       absFactor = stdevfloor;
       not30Ind=0;
   else
       not30Ind=1;
   end
   
   
   %%disp('final abs factor');
   %%disp(absFactor);
   

   topindex=find(Asignal==max(maxima));
   topindex=floor(nanmean(topindex));
   

   numundefsbel =  sum(isnan(Asignal(1:topindex)));
   %%disp('numundefsbel');
   %%disp(numundefsbel);
   numundefabv =  sum(isnan(Asignal(topindex+1:size(sigVector,2)))); 
   
   if numundefsbel > 0
    PBLH = NaN;
   else
     %Set startind to 1. Then set a trap in the for loop so that if the
     %fist max is the lowest point in the profile, it will only compare to
     %the min above. If i is the last element in MaxIndex
     %or last element, it will only compare to the min above
     %or below, respectively.
     startind=1;
     if MaxIndex(1)==1
         %%disp('startind:a');
         startind=2;
     end
     
     if ismember(2, MaxIndex_ni) && hookstat==1
         %%disp('startind:b');
         startind=2;
     end
     
     
     %%disp('startind');
     %%disp(startind);
     
     for i = startind: length(MaxIndex)
       if MaxIndex(i) > topindex & numundefabv > 0
        PBLH = NaN;
        break
       else
          %% min at the SAME index as the max is the min above
          %% min at the index below the max in the min below
              
          %check the min above the max to make sure it
          %is at least 1 abs factor away 
            if  maxima(i) - absFactor > minima(i)

          %%  change above line to this one to check that criterion is met for min above and below candidate max
          %%  if maxima(i) - absFactor > minima(i-1)  && maxima(i) - absFactor > minima(i)
                PBLH = heightVector(MaxIndex(i));
                %%disp('PBLH from alg');
                %%disp(PBLH);
                PBLHsnr=Asignal(MaxIndex(i));
                %%disp('PBLHsnr from alg');
                %%disp(PBLHsnr);
                break
            end
               
       end     %% end of no undefs above test
     end    % end of do loop
     if PBLH>3500
         PBLH=NaN;
     end

   end % end of no undef below test
   
   if absFactor==0
       PBLH=NaN;
   end



end   % end of check for non zero indeces

%whos


end   % end of function
