%%   資料邊界
clear all; close all;clc;
% %%   資料
data = load ('ftp_75_gear.txt');            %匯入工況資料
data=data(1:end,:);
time=data(:,1);
speed=data(:,2);


%% 求出曲線的垂直距離為2的包絡線
 x=time;
 y=speed;
 r=2;
ratio=atan(gradient(y)./gradient(x));
for i=1:length(x)
    k=pi/2+ratio(i);
    a(i)=x(i)+r*cos(k);
    b(i)=y(i)+r*sin(k);
end
%% 求每點豎直線與包絡線的交點
a_m=zeros(length(x),1);
for i=4:length(x)-4
    record=-ones(20,1)*2;
    for j=i-3:i+3
        if(    (  x(i)>min(a(j),a(j+1))|| x(i)==min(a(j),a(j+1))  )  ...
                &&  (x(i)<max(a(j),a(j+1))|| x(i)==max(a(j),a(j+1))   )  )
                m_x=a(j:j+1);
                m_y=b(j:j+1);
                fit_par=polyfit(m_x,m_y,1);% 計算一次函式的斜率 與 b
                record(j-i+9+1)= x(i)*fit_par(1)+fit_par(2);
        end  
    end
    a_m(i)=max(record);
end




%% 補全前面後面的五個點
a_m(1:5)=2;
for i=length(a_m)-7:length(a_m)
     record=-ones(10,1)*2;
    for j=length(a_m)-15:length(a_m)
        if(    (  x(i)>min(a(j-1),a(j))|| x(i)==min(a(j-1),a(j))  )  ...
                &&  (x(i)<max(a(j-1),a(j))|| x(i)==max(a(j-1),a(j))   )  )
                m_x=a(j-1:j);
                m_y=b(j-1:j);
                fit_par=polyfit(m_x,m_y,1);% 計算一次函式的斜率 與 b
                record(j-length(a_m)+15+1)= x(i)*fit_par(1)+fit_par(2);
        end  
    end
    a_m(i)=max(record);
end


%% 求出曲線的垂直距離為2的包絡線
 x=time;
 y=speed;
 r=2;
ratio=atan(gradient(y)./gradient(x));
for i=1:length(x)
    k=pi/2+ratio(i);
    a(i)=x(i)-r*cos(k);
    b(i)=y(i)-r*sin(k);
end
%% 求每點豎直線與包絡線的交點
a_mm=zeros(length(x),1);
for i=4:length(x)-4
    record=ones(20,1)*2000;
    for j=i-3:i+3
        if(    (  x(i)>min(a(j),a(j+1))|| x(i)==min(a(j),a(j+1))  )  ...
                &&  (x(i)<max(a(j),a(j+1))|| x(i)==max(a(j),a(j+1))   )  )
                m_x=a(j:j+1);
                m_y=b(j:j+1);
                fit_par=polyfit(m_x,m_y,1);% 計算一次函式的斜率 與 b
                record(j-i+9+1)= x(i)*fit_par(1)+fit_par(2);
        end  
    end
    a_mm(i)=min(record);
end




%% 補全前面後面的五個點
a_mm(1:5)=-2;
for i=length(a_mm)-7:length(a_mm)
     record=ones(10,1)*2000;
    for j=length(a_mm)-15:length(a_mm)
        if(    (  x(i)>min(a(j-1),a(j))|| x(i)==min(a(j-1),a(j))  )  ...
                &&  (x(i)<max(a(j-1),a(j))|| x(i)==max(a(j-1),a(j))   )  )
                m_x=a(j-1:j);
                m_y=b(j-1:j);
                fit_par=polyfit(m_x,m_y,1);% 計算一次函式的斜率 與 b
                record(j-length(a_mm)+15+1)= x(i)*fit_par(1)+fit_par(2);
        end  
    end
    a_mm(i)=min(record);
end
a_mm(length(a_mm)-3:length(a_mm))=0;


plot(x,y,'r-');
hold on;
% plot(a,b,'k');
plot(x,a_m,'b');
hold on;

plot(x,a_mm,'b');
hold on;