%skyhawk#flyinghawk
%計算最后一個相點的最近相點的位置及最短距離
function [idx,min_d,idx1,min_d1]=nearest_point(tau,m,whlsj,whlsl,P)
%參數說明：
%輸入：m - 嵌入維數， whlsj － 待分析數據， whlsl - 待分析的數據個數， P - 平均循環周期
%      idx - 最后一個相點的最近相點的位置， min_d - 最后一個相點與其最近相點間的距離 (考慮相角)
%      idx1 - 最后一個相點的最近相點的位置， min_d1 - 最后一個相點與其最近相點間的距離 (不考慮相角)
% *可調參數
% P = 5  ; %  &&選擇演化相點距當前點的位置差，即若當前相點為I，則演化相點只能在|I－J|<P的相點中搜尋
min_point=5  ; %&&要求最少搜索到的點數
MAX_CISHU=5 ;  %&&最大增加搜索范圍次數
% global lmd;

% 相空間重構
LAST_POINT = whlsl-(m-1)*tau;  %相點個數
for j=1:LAST_POINT            
    for k=1:m
        Y(k,j)=whlsj((k-1)*tau+j);
    end
end

% 求最大、最小和平均相點距離
max_d = 0.;
min_d = 1.0e+100;
avg_d = 0.;
for i = 1 : LAST_POINT-1
    avg_dd = 0.;
    for j = i+1 : LAST_POINT
        d = 0.;
        for k = 1 : m
            d = d + (Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));
        end
        d = sqrt(d);
        if max_d < d
           max_d = d;
        end
        if min_d > d
           min_d = d;
        end
        avg_dd = avg_dd + d;
    end
    avg_dd = avg_dd/(LAST_POINT-i-1+1);
    avg_d = avg_d + avg_dd;
end
avg_d = avg_d/(LAST_POINT-1);
  
dlt_eps = (avg_d - min_d) * 0.08 ;         % 若在min_eps～max_eps中找不到演化相點時，對max_eps的寬大增幅
min_eps = min_d + dlt_eps / 8 ;            % 演化相點與當前相點距離的最小限
max_eps = min_d + 2 * dlt_eps  ;           % 演化相點與當前相點距離的最大限
    
% 從P～N-m個相點中找與第一個相點最近的相點位置(Loc_DK)及其最短距離DK
DK = 1.0e+100;
Loc_DK = LAST_POINT-P;
for i = 1 : LAST_POINT-P
    d = 0.;
    for k = 1 : m
        d = d + (Y(k,i)-Y(k,LAST_POINT-1))*(Y(k,i)-Y(k,LAST_POINT-1));
    end
    d = sqrt(d);
        
    if (d < DK) & (d > min_eps) 
        DK = d;
        Loc_DK = i;
    end
end

DK1 = 0.;
for k = 1 : m
    DK1 = DK1 + (Y(k,LAST_POINT)-Y(k,Loc_DK+1))*(Y(k,LAST_POINT)-Y(k,Loc_DK+1));
end
DK1 = sqrt(DK1);

old_Loc_DK=Loc_DK;
  
% 以下程序計算最后一個相點的最近距離點  （考慮相角）：要求距離在指定距離范圍內盡量短，與DK1的角度最小

max_eps = min_d + 2 * dlt_eps ;            % 演化相點與當前相點距離的最大限

point_num = 0  ; % 在指定距離范圍內找到的候選相點的個數
cos_sita = 0.  ; % 夾角余弦的比較初值 ——要求一定是銳角
zjfwcs=0       ; % 增加范圍次數
        
while (point_num == 0)
    % 搜索相點
    for j = 1 : LAST_POINT-1
        if abs(j-LAST_POINT) <=( P-1)      % 候選點距當前點太近，跳過！
           continue;      
        end
                
        % 計算候選點與當前點的距離
        dnew = 0.;
        for k = 1 : m
            dnew = dnew + (Y(k,LAST_POINT)-Y(k,j))*(Y(k,LAST_POINT)-Y(k,j));
        end
        dnew = sqrt(dnew);
                
        if (dnew < min_eps)|( dnew > max_eps )   % 不在距離范圍，跳過！
           continue;              
        end                
                
        % 計算夾角余弦及比較
        DOT = 0.;
        for k = 1 : m
            DOT = DOT+(Y(k,LAST_POINT)-Y(k,j))*(Y(k,LAST_POINT)-Y(k,old_Loc_DK+1));
        end
        CTH = DOT/(dnew*DK1);
                
        if acos(CTH) > (3.14151926/4)      % 不是小于45度的角，跳過！
           continue;
        end
                
        if CTH > cos_sita   % 新夾角小于過去已找到的相點的夾角，保留
            cos_sita = CTH;
            Loc_DK = j;
            DK = dnew;
        end
        point_num = point_num +1;
    end  % end of for j = 1 : LAST_POINT-1     
        
    if point_num < min_point
        point_num = 0          ;     %&&擴大距離范圍后創新搜索
        cos_sita = 0.;
        max_eps = max_eps + dlt_eps;
        zjfwcs =zjfwcs +1;
        if zjfwcs > MAX_CISHU    %&&超過最大放寬次數，改找最近的點
           DK = 1.0e+100;
           for ii = 1 : LAST_POINT-1
               if abs(LAST_POINT-ii) <= P-1      %&&候選點距當前點太近，跳過！
                   continue;     
               end
               d = 0.;
               for k = 1 : m
                   d = d + (Y(k,LAST_POINT)-Y(k,ii))*(Y(k,LAST_POINT)-Y(k,ii));
               end
               d = sqrt(d);
        
               if (d < DK) & (d > min_eps) 
                   DK = d;
                   Loc_DK = ii;
               end
           end  % end of for ii = 1 : LAST_POINT-1
           break  
        end %end of if zjfwcs > MAX_CISHU
    end %end of if point_num <= min_point
end  % end of while (point_num == 0)
idx=Loc_DK;  %返回中心點最近相點的位置
min_d=DK;    %返回中心點到其最近相點的距離
point_num;
% 以下程序計算最后一個相點的最近距離點  （不考慮相角）
 
% 求最小距離點
min_d1 = 1e+100;
idx1 = LAST_POINT-1;

for jj = 1:LAST_POINT-1         

    if abs(jj-LAST_POINT) <=( P-1)      % 候選點距當前點太近，跳過！
       continue;      
    end

    sum_d=0.;
    for k=1:m
       sum_d = sum_d+(Y(k,jj)-Y(k,LAST_POINT))^2;
    end
    sum_d = sqrt(sum_d);
    
    if (min_d1 > sum_d)&(sum_d > 0)
        min_d1 = sum_d;
        idx1 = jj;
    end
end