function BinElaJoe = SimulacaoJoelho(MolRigJoeDir1, AngJoe, NumPtsSim)

%% PROGRAMA PARA A SIMULACAO DO MOVIMENTO DO JOELHO
    parametrosfuncao = [MolRigJoeDir1 AngJoe NumPtsSim/2];
    
  %%%%%%%%%%%%%%%%%%%%%%%%%
  %% ConfiguracaoFiguras %%
  %%%%%%%%%%%%%%%%%%%%%%%%%
  figura1 = figure(1);
    set(figura1,'position',[1550 100 360 900])
    hold on; grid on;axis equal
    Title('Animation - 135º')
    %axis([-10 20 165 210]);
  figura2 = figure(2);
    hold on; grid on;
    set(figura2,'position', [1210 760 330 240])%[810 760 530 240])
    axis([0 135 0 135]);
    Title('ComprimentoElastico em funcao do AnguloJoelho')
  figura3 = figure(3);
    hold on; grid on;
    set(figura3,'position',[1210 430 330 240])
    %axis([0 135 0 100]);
    Title('ForcaElastico em funcao do AnguloJoelho')
  figura4 = figure(4);
    hold on; grid on;
    set(figura4,'position',[1210 100 330 240])
    %axis([0 135 0 2]);
    Title('BinarioElastico em funcao do AnguloJoelho')
%%%%%%%%%%%%%%%%%%%%%%%%%
%% ParametrosEstrutura %%
%%%%%%%%%%%%%%%%%%%%%%%%%
  AT = 45;              %AlturaTornozelo
  ACa= 139;             %AlturaCanela
  ACo= 142;             %AlturaCoxa
  AJ = AT+ACa;          %AlturaJoelho
  SECa = [13 AJ-17];    %SuporteElasticosCanela
  SECo = [11.5 AJ+50];  %SuporteElasticosCoxa
  
%%%%%%%%%%%%%%%
%% Contantes %%
%%%%%%%%%%%%%%%
  ConvDeg2Rad = (2*pi)/360;          %ConversaoGrausParaRadianos
  ConvRad2Deg = 360/(2*pi);
  ConvMm2M=1/1000;
  NPC = parametrosfuncao(3);                         %NumeroPontosCircunferencia
  NPE = NPC;                         %NumeroPontosElastico
  Vermelho = [1 0 0];
  Verde    = [0 1 0];
  Azul     = [0 0 1];
  Preto    = [0 0 0];
  %TS = 2;                            %TempoSimulacao
  NPS = NPC;                         %NumeroPontosSimulacao
  AngJoe = parametrosfuncao(2);                      %AnguloMaximoJoelho (AnguloMinimo=0 -> PosicaoVertical)
  %VTS = linspace(0,TS,NPS);          %VectorTempoSimulacao(s)
  VAngJoe = linspace(0,AngJoe,NPS);  %AnguloJoelho(deg)
  ComprimentoElastico(1,1,1)=0;
  k=parametrosfuncao(1);

  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% DesenhosCircunferenciasCanela %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %Cicunferencias
  %%%%%%%%%%%%%%%
    figure(1)
  %EixoRotacaoJoelho
  %%%%%%%%%%%%%%%%%%
    RCJ1 = 3;                                           %RaioCircunferencia(mm)
    ICJ1 = 0;                                           %InicioCircunferencia(º)
    FCJ1 = 360;                                         %FimCircunferencia(º)
    CCXJ1 = 0;                                          %CentroCircunferenciaX(mm)
    CCYJ1 = AJ;                                         %CentroCircunferenciaY(mm)
    x=linspace(ICJ1*ConvDeg2Rad,FCJ1*ConvDeg2Rad,NPC); 
    y=linspace(ICJ1*ConvDeg2Rad,FCJ1*ConvDeg2Rad,NPC);
    X=RCJ1*sin(x)+CCXJ1;
    Y=RCJ1*cos(y)+CCYJ1;
    fill(X,Y,'b')
  %EixoRotacaoElastico1
  %%%%%%%%%%%%%%%%%%%%%
    %RC1 = 10;
    RC1 = 5;                                            %RaioCircunferencia(mm)
    IC1 = 0;                                            %InicioCircunferencia(º)
    FC1 = 360;                                          %FimCircunferencia(º)
    CCX1 = 0+15;                                        %CentroCircunferenciaX(mm)
    CCY1 = AJ+10;                                       %CentroCircunferenciaY(mm)
    x=linspace(IC1*ConvDeg2Rad,FC1*ConvDeg2Rad,NPC);
    y=linspace(IC1*ConvDeg2Rad,FC1*ConvDeg2Rad,NPC);
    X=RC1*sin(x)+CCX1;
    Y=RC1*cos(y)+CCY1;
      TGX1=X;
      TGY1=Y;
    fill(X,Y,'b')
  %EixoRotacaoElastico2
  %%%%%%%%%%%%%%%%%%%%%
    %RC2 = 10;
    RC2 = 3;                                            %RaioCircunferencia(mm)
    IC2 = 0;                                            %InicioCircunferencia(º)
    FC2 = 360;                                          %FimCircunferencia(º)
    CCX2 = 0+6;                                         %CentroCircunferenciaX(mm)
    CCY2 = AJ+15;                                       %CentroCircunferenciaY(mm)
    x=linspace(IC2*ConvDeg2Rad,FC2*ConvDeg2Rad,NPC);
    y=linspace(IC2*ConvDeg2Rad,FC2*ConvDeg2Rad,NPC);
    X=RC2*sin(x)+CCX2;
    Y=RC2*cos(y)+CCY2;
      TGX2=X;
      TGY2=Y;
    fill(X,Y,'b')
  %EixoRotacaoElastico3
  %%%%%%%%%%%%%%%%%%%%%
    %RC3 = 10;
    RC3 = 3;                                            %RaioCircunferencia(mm)
    IC3 = 0;                                            %InicioCircunferencia(º)
    FC3 = -360;                                         %FimCircunferencia(º)
    CCX3 = 0-6;                                         %CentroCircunferenciaX(mm)
    CCY3 = AJ+15;                                       %CentroCircunferenciaY(mm)
    x=linspace(IC3*ConvDeg2Rad,FC3*ConvDeg2Rad,NPC);
    y=linspace(IC3*ConvDeg2Rad,FC3*ConvDeg2Rad,NPC);
    X=RC3*sin(x)+CCX3;
    Y=RC3*cos(y)+CCY3;
      TGX3=X;
      TGY3=Y;
    fill(X,Y,'b')

%%%%%%%%%%%%%%%%%%%%
%% DesenhosCanela %%
%%%%%%%%%%%%%%%%%%%%
  %Desenho da Canela
    figure(1)
    line([0 0],[0 AJ],'LineWidth',4)
    line([0 SECa(1)],[SECa(2) SECa(2)],'LineWidth',4)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% DesenhosCircunferenciasCoxa %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %cicunferencias
  %%%%%%%%%%%%%%%
    figure(1)
  %EixoRotacaoJoelho
  %%%%%%%%%%%%%%%%%%
    RCJ2 = 2;                                           %RaioCircunferencia(mm)
    ICJ2 = 0;                                           %InicioCircunferencia(º)
    FCJ2 = 360;                                         %FimCircunferencia(º)
    ConvDeg2Rad = 2*pi/360;                             %ConversaoGrausParaRadianos
    CCXJ2 = 0;                                          %CentroCircunferenciaX(mm)
    CCYJ2 = AJ;                                         %CentroCircunferenciaY(mm)
    x=linspace(ICJ2*ConvDeg2Rad,FCJ2*ConvDeg2Rad,NPC);
    y=linspace(ICJ2*ConvDeg2Rad,FCJ2*ConvDeg2Rad,NPC);
    X=RCJ2*sin(x)+CCXJ2;
    Y=RCJ2*cos(y)+CCYJ2;
    h(1,1)=fill(X,Y,'g');

%%%%%%%%%%%%%%%%%%
%% DesenhosCoxa %%
%%%%%%%%%%%%%%%%%%
  figure(1)
  h(1,3)=line([0 0],[AJ ACo+AJ],'LineWidth',4,'Color',Verde);
  h(1,4)=line([0 SECo(1)],[SECo(2) SECo(2)],'LineWidth',4,'Color',Verde);

  pause(1)

%%%%%%%%%%%%%%%%%%%%%%%  
%% DESENHO ELASTICOS %%
%%%%%%%%%%%%%%%%%%%%%%%
  %Canela
  %%%%%%%
    for i=1:NPE                     %NumeroPontosElastico
      a=TGX1(i)-SECa(1);
      b=TGY1(i)-SECa(2);
      L=sqrt(a^2+b^2);
      Teta=asin(a/L);
      E=0.5;%Erro
      if TGX1(i)-RC1*cos(Teta)>CCX1-E && TGX1(i)-RC1*cos(Teta)<CCX1+E ...
         && TGY1(i)+RCJ1*sin(Teta)<CCY1+E && TGY1(i)+RCJ1*sin(Teta)>CCY1-E
             line([SECa(1) TGX1(i)],[SECa(2) TGY1(i)],'LineWidth',2,'Color',Vermelho);
             break
      end
    end
  %Coxa
    for i=1:NPE
      a=TGX1(i)-SECo(1);
      b=TGY1(i)-SECo(2);
      L=sqrt(a^2+b^2);
      Teta=asin(a/L);
      E=0.1;%Erro
      if TGX1(i)-RC2*cos(Teta)>CCX1-E && TGX1(i)-RC2*cos(Teta)<CCX1+E
         %&& TGY(i)-3*sin(Teta)>45+139+10-0.1 && TGY(i)-3*sin(Teta)<45+139+10+0.1;
            h(1,5)=line([SECo(1) TGX1(i)],[SECo(2) TGY1(i)],'LineWidth',2,'Color',Vermelho);
            break
      end
    end
    
%%%%%%%%%%%%%%%%%%%%%%%%%
%% DesenhoBracoBinario %%
%%%%%%%%%%%%%%%%%%%%%%%%%
  figure(1)
  h(1,6)=line([0 SECo(1)],[AJ SECo(2)],'LineWidth',2,'Color',Preto);
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CICLO ACTUALIZACAO IMAGEM %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ComprimentoElasticoCanela = sqrt((SECa(1)-CCX1)^2+(SECa(2)-CCY1)^2);    %%Aproximacao
DistanciaBatente1e2=sqrt((CCX1-CCX2)^2+(CCY1-CCY2)^2);                  %%Aproximacao
DistanciaBatente2e3=sqrt((CCX2-CCX3)^2+(CCY2-CCY3)^2);                  %%Aproximacao
%Inicializacao de variaveis
  batente=1;
  i=0;
  c=0;
  ae=1;
  ac=1;
  f=1;
  ang=1;
  bin=1;
for i=1:NPS
    %DesenhoCoxa
    %%%%%%%%%%%%
      XDataCoxa=-ACo*sin(VAngJoe(1,i)*ConvDeg2Rad);   %NovaPosicaoXExtremidadeCoxa
      set(h(1,3),'XData',[0 XDataCoxa]);
      YDataCoxa=AJ+ACo*cos(VAngJoe(1,i)*ConvDeg2Rad); %NovaPosicaoYExtremidadeCoxa
      set(h(1,3),'YData',[AJ YDataCoxa]);
    %DesenhoSuporteElasticoCoxa
    %%%%%%%%%%%%%%%%%%%%%%%%%%%
      L=sqrt(SECo(1)^2+(SECo(2)-AJ)^2);               %DistanciaEntreJoelhoEPontaSuporte
      TetaSECo=asin(SECo(1)/L)*ConvRad2Deg;           %Angulo(º)EntreCoxaEPontaSuporteElasticosJoelho(BracoBinario)
      XDataSECo1=(SECo(2)-AJ)*sin(-VAngJoe(1,i)*ConvDeg2Rad);
      XDataSECo2=L*sin((-VAngJoe(1,i)+TetaSECo)*ConvDeg2Rad);
        set(h(1,4),'XData',[XDataSECo1 XDataSECo2]);
      YDataSECo1=(SECo(2)-AJ)*cos(-VAngJoe(1,i)*ConvDeg2Rad);
      YDataSECo2=L*cos((-VAngJoe(1,i)+TetaSECo)*ConvDeg2Rad);
        set(h(1,4),'YData',[YDataSECo1+AJ YDataSECo2+AJ]);
    %DesenhoBracoBinario
    %%%%%%%%%%%%%%%%%%%%
      XDataBraco1=0;
      XDataBraco2=XDataSECo2;
      set(h(1,6),'XData',[XDataBraco1 XDataBraco2]);
      YDataBraco1=AJ;
      YDataBraco2=YDataSECo2;
      set(h(1,6),'YData',[YDataBraco1 YDataBraco2+AJ]);
      ComprimentoBracoBinario = sqrt( SECo(1)^2 + (SECo(2)-AJ)^2 );
      TetaBracoBinario = asin(SECo(1)/ComprimentoBracoBinario)*ConvRad2Deg;
    %DesenhoElasticos - ElasticoCoxa
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %Batente1 - 1ªSaliencia do Joelho
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if batente==1
          for j=1:NPS
            a = TGX1(j) - XDataSECo2;
            b = AJ+YDataSECo2-TGY1(j);
            L = sqrt(a^2+b^2);
            Teta = asin(a/L);
            E=0.2;
            if TGX1(j)-RC1*cos(Teta) > CCX1-E && TGX1(j)-RC1*cos(Teta)<CCX1+E
              AnguloElastico(1,ae)=Teta*ConvRad2Deg;    ae=ae+1;
              set(h(1,5),'XData',[XDataSECo2 TGX1(j)]);
              set(h(1,5),'YData',[YDataSECo2+AJ TGY1(j)]);
              c=c+1;
              ComprimentoElastico(c,1)= ...
                  ComprimentoElasticoCanela+sqrt((XDataSECo2-TGX1(j))^2+(YDataSECo2+AJ-TGY1(j))^2);
              AumentoComprimentoElastico(i,1) = ComprimentoElastico(i,1)-ComprimentoElastico(1,1); ac=ac+1;
              ForcaElastico(i,1)=k*AumentoComprimentoElastico(ac-1,1)*ConvMm2M;
              AnguloElastico(1,i)+TetaBracoBinario-VAngJoe(1,i);
              ForcaElasticoX(i,1)=ForcaElastico(i,1) * sin((AnguloElastico(1,i)+TetaBracoBinario-VAngJoe(1,i))*ConvDeg2Rad);
              ForcaElasticoY(i,1)=ForcaElastico(i,1) * cos((AnguloElastico(1,i)+TetaBracoBinario-VAngJoe(1,i))*ConvDeg2Rad);      f=f+1;
              Binario(bin,1)=ComprimentoBracoBinario*ConvMm2M*ForcaElasticoX(f-1,1); bin=bin+1;
              ang=ang+1;
            break %Sai do ciclo FOR --
          end     %                   |
        end       %                   |
                  %               <<--
        x1=XDataSECo2;              x2=TGX1(j);
        y1=YDataSECo2+AJ;           y2=TGY1(j);
        m = (y1-y2)/(x1-x2);        b1=y1-x1*m;     b2=y2-x2*m;
        X=linspace(x1,x2,NPC);     Y=m*X+b1;
        for j=1:1:NPS
          if sqrt((X(j)-CCX2)^2+(Y(j)-CCY2)^2) < RC2
              figure(2)
              line([VAngJoe(1,i) VAngJoe(1,i)],[0 AngJoe],'Color',Azul,'LineWidth',2)

              figure(3)
              line([VAngJoe(1,i) VAngJoe(1,i)],[0 max(ForcaElastico)],'Color',Azul,'LineWidth',2)

              figure(4)
              line([VAngJoe(1,i) VAngJoe(1,i)],[0 max(Binario)],'Color',Azul,'LineWidth',2)
            batente=2;
            break %Sai do ciclo FOR --
          end     %                   |
        end       %                   |
    end           %               <<--
    
    %Batente2
    %%%%%%%%%
    if batente==2
      for j=1:NPS
        a=TGX2(j)-XDataSECo2;
        b=AJ+YDataSECo2-TGY2(j);
        L=sqrt(a^2+b^2);
        Teta=asin(a/L);
        E=0.2;
        if TGX2(j)-RC2*cos(Teta)>CCX2-E && TGX2(j)-RC2*cos(Teta)<CCX2+E
          AnguloElastico(1,ae)=Teta*ConvRad2Deg; ae=ae+1;
          set(h(1,5),'XData',[XDataSECo2 TGX2(j)] );
          set(h(1,5),'YData',[YDataSECo2+AJ TGY2(j)] );
          ComprimentoElastico(c,1)=...
              ComprimentoElasticoCanela + DistanciaBatente1e2 + ...
              sqrt((XDataSECo2-TGX2(j))^2+(YDataSECo2+AJ-TGY2(j))^2);   c=c+1;
          AumentoComprimentoElastico(i,1)=ComprimentoElastico(i,1)-ComprimentoElastico(1,1); ac=ac+1;
          ForcaElastico(i,1)=k*AumentoComprimentoElastico(i,1)*ConvMm2M;
          AnguloElastico(1,i)+TetaBracoBinario-VAngJoe(1,i);
          ForcaElasticoX(i,1)=ForcaElastico(i,1) * sin((AnguloElastico(1,i)+TetaBracoBinario-VAngJoe(1,i))*ConvDeg2Rad);
          ForcaElasticoY(i,1)=ForcaElastico(i,1) * cos((AnguloElastico(1,i)+TetaBracoBinario-VAngJoe(1,i))*ConvDeg2Rad);      f=f+1;
          Binario(i,1)=ComprimentoBracoBinario*ConvMm2M*ForcaElasticoX(i,1); bin=bin+1;
          ang=ang+1;
          break %Sai do ciclo FOR --
        end     %                   |
      end       %                   |
                %               <<--
      x1=XDataSECo2;              x2=TGX2(j);
      y1=YDataSECo2+AJ;           y2=TGY2(j);
      m = (y1-y2) / (x1-x2);      b1=y1-x1*m;       b2=y2-x2*m;
      X=linspace(x1,x2,NPC);     Y=m*X+b1;
      for j=1:1:NPS
        if sqrt((X(j)-CCX3)^2+(Y(j)-CCY3)^2) < RC3
                figure(2)
                line([VAngJoe(1,c) VAngJoe(1,c)],[0 AngJoe],'Color',Azul,'LineWidth',2)

                figure(3)
                line([VAngJoe(1,c) VAngJoe(1,c)],[0 max(ForcaElastico)],'Color',Azul,'LineWidth',2)

                figure(4)
                line([VAngJoe(1,c) VAngJoe(1,c)],[0 max(Binario)],'Color',Azul,'LineWidth',2)
          batente=3;
          break %Sai do ciclo FOR --
        end     %                   |
      end       %                   |
    end         %               <<--
    %Batente3
    %%%%%%%%%
    if batente==3
      for j=1:NPS
        a=TGX3(j)-XDataSECo2;
        b=AJ+YDataSECo2-TGY2(j);
        L=sqrt(a^2+b^2);
        Teta=(90*ConvDeg2Rad)-asin(a/L)+(90*ConvDeg2Rad);
        E=0.1;
        if TGY3(j)-RC3*sin(Teta)>CCY3-E && TGY3(j)-RC3*sin(Teta)<CCY3+E
           %TGY3(j)-RC3*cos(Teta-pi)>CCY3-E && TGY3(j)-RC3*cos(Teta-pi)>CCY3+E
          AnguloElastico(1,ae)=Teta*ConvRad2Deg;    ae=ae+1;
          set(h(1,5),'XData',[XDataSECo2 TGX3(j)] );
          set(h(1,5),'YData',[YDataSECo2+AJ TGY3(j)] );
          ComprimentoElastico(c,1)=...
              ComprimentoElasticoCanela + DistanciaBatente1e2 + DistanciaBatente2e3 + ...
              sqrt((XDataSECo2-TGX3(j))^2+(YDataSECo2+AJ-TGY3(j))^2);   c=c+1;
          AumentoComprimentoElastico(i,1)=ComprimentoElastico(i,1)-ComprimentoElastico(1,1); ac=ac+1;
          ForcaElastico(i,1)=k*AumentoComprimentoElastico(i,1)*ConvMm2M;
          AnguloElastico(1,i)+TetaBracoBinario-VAngJoe(1,i);
          ForcaElasticoX(i,1)=ForcaElastico(i,1) * sin((AnguloElastico(1,i)+TetaBracoBinario-VAngJoe(1,i))*ConvDeg2Rad);
          ForcaElasticoY(i,1)=ForcaElastico(i,1) * cos((AnguloElastico(1,i)+TetaBracoBinario-VAngJoe(1,i))*ConvDeg2Rad);      f=f+1;
          Binario(i,1)=ComprimentoBracoBinario*ConvMm2M*ForcaElasticoX(i,1); bin=bin+1;
          ang=ang+1;
          break %Sai do ciclo FOR --
        end     %                   |
      end       %                   |
    end         %               <<--
end

figure(2)
for i=1:NPS
    ComprimentoElasticoRedimensionado(1,i)=ComprimentoElastico(i,1);
end
plot(VAngJoe,ComprimentoElasticoRedimensionado,'r');
figure(3)
plot(VAngJoe,ForcaElastico,'r')
figure(4)
plot(VAngJoe,Binario,'r')

j=NumPtsSim;
for i=1:1:NumPtsSim
  if i <= NumPtsSim/2
    Bin(i,1) = Binario(i);
  end
  if  i > NumPtsSim/2
    Bin(i,1) = Binario(j);
  end
  j=j-1;
end
BinElaJoe = Bin;


