clc;
format long
f_c=4.8e-5
f_t=4.1666667e-3
N_P=5;

Mx(1)=0;
Mf(1)=0;
for b=1:N_E
    for z=1:N_T(b)
        a=T(N(N_T_R(b)+z)+1);
        for j=1:N_Q(N_T_R(b)+z)
            T(N(N_T_R(b)+z)+N_Q(N_T_R(b)+z)-j+1)=T(N(N_T_R(b)+z)+N_Q(N_T_R(b)+z)-j+1)-a;      
        end
        for i=1:N_P
            for p=1:N_P   
                A(i,p)=0;
            end
            B(i)=0;
        end
        for i=1:N_P
            for j=1:N_Q(N_T_R(b)+z)    
                B(i)=B(i)+X(N(N_T_R(b)+z)+j)*T(N(N_T_R(b)+z)+j)^(i-1);
            end
            for p=1:N_P
                for j=1:N_Q(N_T_R(b)+z)    
                    A(i,p)=A(i,p)+T(N(N_T_R(b)+z)+j)^(i+p-2);
                end
            end
        end
        A;
        B=transpose(B);
        C=linsolve(A,B);
        
        for j=1:N_P
            CF(N_T_R(b)*N_P+N_P*(z-1)+j)=C(j);
        end
        inc=0.001;
        max=round(T(N(N_T_R(b)+z)+N_Q(N_T_R(b)+z))/inc);
        for j=1:max
            T_G_P(Mx(N_T_R(b)+z)+j)=j*inc;
        end
        for j=1:max
            X_G_P(Mx(N_T_R(b)+z)+j)=0;
            for i=1:N_P
                X_G_P(Mx(N_T_R(b)+z)+j)=X_G_P(Mx(N_T_R(b)+z)+j)+C(i)*T_G_P(Mx(N_T_R(b)+z)+j)^(i-1);
            end
        end
        B=transpose(B);
        y_T=0;
        for j=1:N_Q(N_T_R(b)+z)
            y_T=y_T+T(N(z)+j);        
        end
        y_m=y_T/N_Q(N_T_R(b)+z);
        S_1=0;
        for j=1:(N_Q(N_T_R(b)+z))
            x=round(T(N(N_T_R(b)+z)+j)/inc);
            if(Mx(N_T_R(b)+z)+x+1>Mx(N_T_R(b)+z)+max)
                x=max-1;
            end
            S_1=S_1+(X(N(N_T_R(b)+z)+j)-X_G_P(Mx(N_T_R(b)+z)+x+1))^2;
        end
        S_2=0;
        for j=1:N_Q(N_T_R(b)+z)
            S_2=S_2+(X(N(N_T_R(b)+z)+j)-y_m)^2;        
        end
        
        R(N_T_R(b)+z)=1-S_1/S_2;
        Mx(N_T_R(b)+z+1)=Mx(N_T_R(b)+z)+max;
        
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    tf=5;
    for i=1:N_T(b)
        if (tf>T_G_P(Mx(N_T_R(b)+i+1)))
            tf=T_G_P(Mx(N_T_R(b)+i+1));
        end
    end 
        for i=1:N_P
            for p=1:N_P   
                A(i,p)=0;
            end
            B(i)=0;
        end
        for i=1:N_P
            for j=(N(N_T_R(b)+1)+1):1:N(N_T_R(b+1)+1)
                if T(j)<tf
                    B(i)=B(i)+X(j)*T(j)^(i-1);
                end
            end
            for p=1:N_P
                for j=(N(N_T_R(b)+1)+1):1:N(N_T_R(b+1)+1)
                    if T(j)<tf
                        A(i,p)=A(i,p)+T(j)^(i+p-2);
                    end
                end
            end
        end
        A;
        B=transpose(B);
        C=linsolve(A,B); 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        ver=0;
        if C(N_P)>0
            N_P=N_P-1;
            ver=1;
            B=transpose(B);
            for i=1:N_P
            for j=(N(N_T_R(b)+1)+1):1:N(N_T_R(b+1)+1)
                if T(j)<tf
                    B(i)=B(i)+X(j)*T(j)^(i-1);
                end
            end
            for p=1:N_P
                for j=(N(N_T_R(b)+1)+1):1:N(N_T_R(b+1)+1)
                    if T(j)<tf
                        A(i,p)=A(i,p)+T(j)^(i+p-2);
                    end
                end
            end
        end
        A;
        B=transpose(B);
        C=linsolve(A,B);            
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        for j=1:N_P
            CM_F((b-1)*N_P+j)=C(j);            
        end
        for j=1:N_P
            CV((b-1)*N_P+j)=C(j)*(j-1);            
        end
        inc=0.002;
        max=round(10/inc);
        count=0;
        for j=1:max
            T_G_F(Mf(b)+j)=j*inc;
            X_G_F(Mf(b)+j)=0;
            for i=1:N_P
                X_G_F(Mf(b)+j)=X_G_F(Mf(b)+j)+C(i)*T_G_F(Mf(b)+j)^(i-1);
            end        
            V_G_F(Mf(b)+j)=0;
            for i=1:(N_P)
                V_G_F(Mf(b)+j)=V_G_F(Mf(b)+j)+CV((b-1)*N_P+i)*T_G_F(Mf(b)+j)^(i-2);
            end
            if count>0.5 && X_G_F(Mf(b)+j)<X_G_F(Mf(b)+j-1)
                break;
            end
            count=count+1;
        end
        B=transpose(B);
        Mf(b+1)=Mf(b)+count;
        if ver==1
             N_P=N_P+1;
        end
        b 
end
if 0
figure(3)
title('Para os 6 primeiros tempos de capacitor testado:\nx(t) da câmera já tratado, em preto\n aproximação para essas curvas, em vermelho\n Por fim a média x(t), v(t) e v(x) em azul, verde e roxo, respectivamente.')
for b=1:6
    subplot(3,2,b);
    for z=1:N_T(b) 
        for j=1:(N_Q(N_T_R(b)+z)-1) 
        plot([T(N(N_T_R(b)+z)+j) T(N(N_T_R(b)+z)+j+1)],[X(N(N_T_R(b)+z)+j) X(N(N_T_R(b)+z)+j+1)],'-k')
        hold on
        end
        for j=1:(Mx(N_T_R(b)+z+1)-Mx(N_T_R(b)+z)-1) 
        plot([T_G_P(Mx(N_T_R(b)+z)+j) T_G_P(Mx(N_T_R(b)+z)+j+1)],[X_G_P(Mx(N_T_R(b)+z)+j) X_G_P(Mx(N_T_R(b)+z)+j+1)],'-r')
        end
    end
    for j=1:(Mf(b+1)-Mf(b)-1) 
        plot([T_G_F(Mf(b)+j) T_G_F(Mf(b)+j+1)],[X_G_F(Mf(b)+j) X_G_F(Mf(b)+j+1)],'-b')
    end
    for j=1:(Mf(b+1)-Mf(b)-1) 
        plot([T_G_F(Mf(b)+j) T_G_F(Mf(b)+j+1)],[V_G_F(Mf(b)+j) V_G_F(Mf(b)+j+1)],'-g')
    end
    for j=1:(Mf(b+1)-Mf(b)-2) 
        plot([X_G_F(Mf(b)+j) X_G_F(Mf(b)+j+1)],[V_G_F(Mf(b)+j) V_G_F(Mf(b)+j+1)],'-m')
    end
    
    hold off
    b
end
end
if 1
figure(7)
for b=1:6
    subplot(3,2,b);
    for z=1:N_T(b) 
        for j=1:(N_Q(N_T_R(b)+z)-1) 
        plot([T(N(N_T_R(b)+z)+j) T(N(N_T_R(b)+z)+j+1)],[Vel_(N(N_T_R(b)+z)+j) Vel_(N(N_T_R(b)+z)+j+1)],'-r')
        plot([T(N(N_T_R(b)+z)+j) T(N(N_T_R(b)+z)+j+1)],[Vel_k(N(N_T_R(b)+z)+j) Vel_k(N(N_T_R(b)+z)+j+1)],'-b')
        hold on 
        %M=getframe()
        end
        z
    end
    for j=1:2:(Mf(b+1)-Mf(b)-1) 
        plot([T_G_F(Mf(b)+j) T_G_F(Mf(b)+j+1)],[V_G_F(Mf(b)+j) V_G_F(Mf(b)+j+1)],'-k')
    end
    for j=1:2:(Mf(b+1)-Mf(b)-1) 
        plot([T_G_F(Mf(b)+j) T_G_F(Mf(b)+j+1)],[X_G_F(Mf(b)+j) X_G_F(Mf(b)+j+1)],'-m')
    end
    hold off
    b
   
end
end