clear all
clc
close all

%Oteviraci napětí
%    S umělým zatížením na základě souboru info.dat
%
%     posY    - originální data posuvů ve směru Y z Ansysu
%     nody    - čísla nodů na kterých byly posuvy měřeny
%     stress  - zatěžující napětí pro daný substep
%     
%     pocSt   - počet substepů
%     St      - matice substepů
%     pccelk    - počet cyklů
%     posY2   - posuvy před a po otevření trhliny [mm] (hodnota, čas)
%     OTVNAP  - otvírací napětí [MPA] (hodnota, čas)


%%

%načtení vysledků
posY=importdata('deformace.dat');
%nody=importdata('nody.dat');
stress=importdata('zatez.dat');
info=importdata('info.dat');


%pocet stepu
pocSt=length(posY);
St=(1:1:pocSt)';

%pocet cyklu 
pccelk=info(6,1);                
pc=3*pccelk;

%zatížení na základě info.dat
w=info(1,1);
a0=info(3,1)*0.001;
da=info(5,1);
R=info(9,1);
sub=info(7,1);               
sub2=info(8,1);                     
pocsub = 2*sub+sub2;            %substepu na cyklus
Smin0=info(10,1) ;                 
Smax0=info(11,1) ;
deltaK=info(12,1);
Kmin=info(13,1);
Kmax=info(14,1);
Smink=info(19,1);
Smaxk=info(20,1);

% smer=(Smax-Smin)/sub;            %smernice rustu/poklesu napeti

zat=zeros(pocSt,1);

a=a0;
ppc=0;                  %pocet proběhlých cyklů
Smaxmat=zeros(pc,1);    %matice maximalnich napětí
Sminmat=zeros(pc,1);    %matice minimalnich napětí
%kontrolni matice
amat=zeros(pc,1);
alfmat=zeros(pc,1);

for i=1:pc
    alfa=a/(w*0.001);
    fi=sqrt(1/cos(0.5*alfa*pi));
    Smax = Kmax/(sqrt(pi*a)*fi);
    Smin = R*Smax;
    Smaxmat(i,1)=Smax;
    Sminmat(i,1)=Smin;
    amat(i,1)=a;
    alfmat(i,1)=alfa;
    
    zb=mod(i,3);
    if zb==1
%         zat(ppc*pocsub+1:ppc*pocsub+sub)=linspace(Smin,Smax,sub);
        for q=1:sub
            zat(ppc*pocsub+q)=Smin+q*((Smax-Smin)/sub);
        
        end
    elseif zb==2
        for q2=1:sub2
            zat(ppc*pocsub+sub+q2)=Smax;
        end
%         zat(ppc*pocsub+sub+1:ppc*pocsub+sub+sub2)=linspace(Smax,Smax,sub2);
            a=a+da*0.001;
    else 
        for q3=1:sub
            if i-1==0
                i2=1;
            else
                i2=i;
            end
%             zat(ppc*pocsub+sub+sub2+q3)=Smax-q3*((Smax-Smin)/sub);
            zat(ppc*pocsub+sub+sub2+q3)=Smaxmat(i2-1,1)-q3*((Smaxmat(i2-1,1)-Smin)/sub);
        end
          ppc=ppc+1;
%         zat(ppc*pocsub+sub+sub2:ppc*pocsub+2*sub+sub2)=linspace(Smax,Smin,sub);
    end
    
    
    
end

%kontrola
figure(3)
plot(St,stress,'r-');
hold on
plot(St,zat,'g-');

smer2=zeros(length(zat));
for q4=1:length(zat)-1
    smer2(q4)=zat(q4+1)-zat(q4);
end

% zat(1,1)=Smin+smer;
% 
% for i=2:pocSt
%     puc=floor((i-1)/pocsub);                  %pocet uplnynulych cyklu
%     cykl=i-puc*pocsub;                         %jakoby jediny cyklus
%     if cykl<=sub               %LS 1
%         zat(i)=zat(i-1)+smer;
%     elseif cykl<=sub+sub2            %LS 2
%         zat(i)=zat(i-1);
%     else                                    %LS3
%         zat(i)=zat(i-1)-smer;
%         
%     end
% end



%vyber posuvů pred a po otevreni trhliny - posY2
posY2=zeros(pccelk-1,2);
m=1;
for k=1:pocSt-1;
    if posY(k,1)<0 && posY(k+1,1)>=0
        posY2(m,:)=[posY(k,1),k];
        posY2(m+1,:)=[posY(k+1,1),k+1];
        m=m+2;
    else
    end   
end
 

 %Oteviraci napeti 
OTVNAP=zeros(pccelk-1,2);
t=zeros(pccelk-1,2);
j=1;
 for n=1:pccelk-1;
     t(n,:)= [abs(posY2(j,1))/(posY2(j+1,1)-posY2(j,1)),posY2(j,2)];
      OTVNAP(n,:)=[zat(t(n,2),1)+(t(n,1))*(zat(t(n,2)+1,1)...
          -zat(t(n,2),1)),t(n,2)+t(n,1)] ;     
       j=j+2;
 end

%vykresleni mist odecteni "casu" otevreni
figure(1)
plot([0;St],[0;posY],'b-',posY2(:,2),posY2(:,1),'bx')
hold on
plot(t(:,1)+t(:,2),0,'ko','linewidth',2)
hold on
plot([0,pocSt],[0,0],'c--')
grid on
title('Posuv Y')
xlabel('substep'),ylabel('posuv Y [mm]') 
legend('posuv Y','posuvy kolem ','čas otevírání')
print('PosuvY','-dpng')

%vykresleni odecitani napeti
figure(2) 
plot([0;St],[0;zat],'r-')
hold on
plot(OTVNAP(:,2),OTVNAP(:,1),'ko','linewidth',2)
grid on
title('Otvírací napětí')
xlabel('substep'),ylabel('napětí [MPa]')
legend('zatěžující napětí','otvírací napětí')
print('OtviraciNapeti6','-dpng')

%výpis
disp('Otvírací napětí jsou:')
disp('Otv. nap. - čas')
disp(OTVNAP)

%poměry f a u
% f = Sotviraci/Smax v danem cyklu
% u = (1-f)/(1-R)
f=zeros(pccelk-1,1);
for q5=2:pccelk
    f(q5-1,1)=OTVNAP(q5-1,1)/Smaxmat((q5-1)*3+1,1);
%     f2(q5-1,1)=OTVNAP(q5-1,1)/(Smaxmat((q5-1)*3+1,1)-Sminmat((q5-1)*3+1,1)); %pomer k celemu rozkmitu napeti
end
f=[0;f];

u=zeros(pccelk,1);
for q6=2:length(f)
    u(q6)=(1-f(q6))/(1-R);
end

%průměr posledních 
prumcast=0.6;               %cast od ktere se to bude prumerovat
prumod=round(prumcast*pccelk);
prumfmat=f(prumod:end,1);   %cast f k prumerovani
prumumat=u(prumod:end,1);   %cast u k prumerovani
prumf=mean(prumfmat);
prumu=mean(prumumat);

figure(4)
plot(f,'rx')
hold on
plot(u,'bo')
hold on
plot([prumod,pccelk],[prumf,prumf],'r-')
hold on
plot([prumod,pccelk],[prumu,prumu],'b-')
grid on
axis([0 pccelk 0 1.1*max([f;u])])
title('Pomery f a u')
xlabel('zatezovy cyklus'),ylabel('u,f [-]')
legend('f = Sotv/Smax','u = (1-f)/(1-R)','prumer f','prumer u','Location','SouthEast')
print('OtviraciPomery6','-dpng')

%zápis OTVNAP do souboru
dlmwrite('OTVNAP6.txt',OTVNAP,'newline','pc','delimiter','\t')
dlmwrite('pomerOTVNAP6.txt',[prumf;prumu],'newline','pc','delimiter','\t')




