clc
vf=0.4; %max 11
dmax=2.5;
x=1;
N_Qf=N_E;
z0=0;
x0=0;
for z=1:N_Qf
    while 1
        if(V_G_F(Mf(z)+x)>vf && V_G_F(Mf(z)+x+1)<vf)
            x
            z
            X_f(z-z0)=X_G_F(Mf(z)+x);
            T_f(z-z0)=T_G_F(Mf(z)+x);
            if X_f(z-z0)>dmax
                x0=x0+1;
            end
            x=1;            
            break;
        end
        x=x+1;
        if(x>(Mf(z+1)-Mf(z)))
            %X_f(z-z0)=0;
            %T_f(z-z0)=0;
            x=1;
            z0=z0+1;
            z0            
            break;
        end
    end
end
N_Qf=N_Qf-z0-x0;
X_f
T_f
N_P_q=4;
a=X_f(1)
b=T_C(1)
c=T_f(1)
for i=1:N_Qf
    X_f(i)=X_f(i)-a;
    T_C(i)=T_C(i)-b;   
    T_f(i)=T_f(i)-c; 
end  
      
        for i=1:N_P_q
            for p=1:N_P_q   
                Af(i,p)=0;
            end
            Bf(i)=0;
        end
        for i=1:N_P_q
            for j=1:N_Qf    
                Bf(i)=Bf(i)+T_C(j)*X_f(j)^(i-1);
            end
            for p=1:N_P_q
                for j=1:N_Qf    
                    Af(i,p)=Af(i,p)+X_f(j)^(i+p-2);
                end
            end
        end
        Af;
        Bf=transpose(Bf);
        C_c=linsolve(Af,Bf);
        inc=0.01;
        max=round(X_f(N_Qf)/inc);
        ini=round(X_f(1)/inc);
        for j=ini:max
            X_f_c(j-ini+1)=j*inc;
        end
        for j=ini:max
            T_C_c(j-ini+1)=0;
            for i=1:N_P_q
                T_C_c(j-ini+1)=T_C_c(j-ini+1)+C_c(i)*X_f_c(j-ini+1)^(i-1);
            end
        end
        Bf=transpose(Bf);
        y_T=0;
        for j=1:N_Qf
            y_T=y_T+X_f_c(j);        
        end
        y_m=y_T/N_Qf;
        S_1=0;
        for j=1:N_Qf
            x=round(X_f(j)/inc);
            if(x>max)
                x=max-1;
            end
            S_1=S_1+(T_C(j)-T_C_c(x-ini+1))^2;
        end
        S_2=0;
        for j=1:N_Qf
            S_2=S_2+(T_C(j)-y_m)^2;        
        end        
        Rc=1-S_1/S_2
        C_c
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:N_P_q
            for p=1:N_P_q   
                Af(i,p)=0;
            end
            Bf(i)=0;
        end
        for i=1:N_P_q
            for j=1:N_Qf    
                Bf(i)=Bf(i)+T_f(j)*X_f(j)^(i-1);
            end
            for p=1:N_P_q
                for j=1:N_Qf    
                    Af(i,p)=Af(i,p)+X_f(j)^(i+p-2);
                end
            end
        end
        Af;
        Bf=transpose(Bf);
        C_t=linsolve(Af,Bf);
       
        for j=ini:max
            T_f_c(j-ini+1)=0;
            for i=1:N_P_q
                T_f_c(j-ini+1)=T_f_c(j-ini+1)+C_t(i)*X_f_c(j-ini+1)^(i-1);
            end
        end
        Bf=transpose(Bf);
        y_T=0;
        for j=1:N_Qf
            y_T=y_T+X_f_c(j);        
        end
        y_m=y_T/N_Qf;
        S_1=0;
        for j=1:N_Qf
            x=round(X_f(j)/inc);
            if(x>max)
                x=max-1;
            end
            S_1=S_1+(T_f(j)-T_f_c(x-ini+1))^2;
        end
        S_2=0;
        for j=1:N_Qf
            S_2=S_2+(T_f(j)-y_m)^2;        
        end        
        Rt=1-S_1/S_2
        C_t        
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:N_Qf
    X_f(i)=X_f(i)+a;
    T_C(i)=T_C(i)+b;
    T_f(i)=T_f(i)+c;
end
for i=ini:max
    X_f_c(i-ini+1)=X_f_c(i-ini+1)+a;
    T_C_c(i-ini+1)=T_C_c(i-ini+1)+b; 
    T_f_c(i-ini+1)=T_f_c(i-ini+1)+c;
end        
    
    figure(4)
    for j=1:N_Qf-1 
    plot([X_f(j) X_f(j+1)],[T_C(j) T_C(j+1)],'-k')
    hold on
    end
    for j=ini:(max-1)
    plot([X_f_c(j-ini+1) X_f_c(j-ini+2)],[T_C_c(j-ini+1) T_C_c(j-ini+2)],'-r')    
    end    
    hold off
    figure(5)
    for j=1:N_Qf-1
    plot([X_f(j) X_f(j+1)],[T_f(j) T_f(j+1)],'-b')
    hold on
    end
    for j=ini:(max-1)
    plot([X_f_c(j-ini+1) X_f_c(j-ini+2)],[T_f_c(j-ini+1) T_f_c(j-ini+2)],'-g')
    end
    hold off
fprintf('\nEquação para força do chute:');
fprintf(' %d',C_c(1)+b);
for i=2:N_P_q
    fprintf(' %d*d^%d',C_c(i),i-1);
end
fprintf('\n');
fprintf('\nEquação para tempo do chute:');
fprintf(' %d',C_t(1)+c);
for i=2:N_P_q
    fprintf(' %d*d^%d',C_t(i),i-1);
end
fprintf('\n');
fprintf('\nNumero de pontos utilizados: %d\n',N_Qf);
fprintf('Numero de pontos com vi<vf: %d\n',z0-1);
fprintf('Numero de pontos com xd<dmax: %d\n',x0-1);