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
%     Zcoord  - souradnice uzlu po tloustce
%     
%     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('deformace3D.dat');
%nody=importdata('nody.dat');
stress=importdata('zatez.dat');
info=importdata('info3D.dat');
Zcoords=importdata('Zcoords.dat');


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

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

%zatížení na základě info3D.dat
w=info(1,1);
a0=info(4,1)*0.001; %[m]
da=info(6,1);
R=info(10,1);
sub=info(8,1);               
sub2=info(9,1);                     
pocsub = 2*sub+sub2;            %substepu na cyklus
Smin0=info(14,1) ;                 
Smax0=info(15,1) ;
deltaK=info(11,1);
Kmin=info(12,1);
Kmax=info(13,1);
Smink=info(16,1);
Smaxk=info(17,1);
pocuzluZ=size(Zcoords,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(size(zat,1),1);
for q4=1:length(zat)-1
    smer2(q4)=zat(q4+1)-zat(q4);
end

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


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

 %Oteviraci napeti 
OTVNAP=zeros(pccelk-1,pocuzluZ);
time=zeros(pccelk-1,pocuzluZ);
t1=zeros(pccelk-1,pocuzluZ);
for i=1:pocuzluZ
j=1;
 for n=1:pccelk-1
     t1(n,i)= abs(posY2(j,i))/(posY2(j+1,i)-posY2(j,i));
     if isnan(t1(n,i))
         t1(n,i)=0;
     else
      OTVNAP(n,i)=zat(pozposY2(j,i),1)+(t1(n,i))*(zat(pozposY2(j,i)+1,1)...
          -zat(pozposY2(j,i),1)) ; 
      time(n,i)=pozposY2(j,i)+t1(n,i);
       j=j+2;
     end
 end
end

%% GRAFY
%vykresleni mist odecteni "casu" otevreni
Zcoordsmat=repmat(Zcoords(:,2)',pocSt,1);
Stmat=repmat(St,1,pocuzluZ);
Zcoordsmat2=repmat(Zcoords(:,2)',size(posY2,1),1);
Zcoordsmat3=repmat(Zcoords(:,2)',size(time,1),1);
Zcoordsmat4=repmat(Zcoords(:,2)',pccelk,1);
casynuly=zeros(pccelk-1,pocuzluZ);
plane0=zeros(2,2,3);
plane0(:,:,1)=[0,Zcoords(pocuzluZ,2);
               0,Zcoords(pocuzluZ,2)];
plane0(:,:,2)=[ 0, 0;
                pocSt,pocSt];
            
zatmat=repmat(zat,1,pocuzluZ);

figure(1)
plot3(Zcoordsmat,Stmat,posY,'b-')
hold on
plot3(Zcoordsmat2,pozposY2,posY2,'bx')
hold on
plot3(Zcoordsmat3,time,casynuly,'ko','linewidth',2)
hold on
surf(plane0(:,:,1),plane0(:,:,2),plane0(:,:,3),'FaceAlpha',0.5)
hold on

pleg11 = plot3(0,0,0,'b-','Visible','off');
pleg12 = plot3(0,0,0,'bx','Visible','off');
pleg13 = plot3(0,0,0,'ko','linewidth',2,'Visible','off');

grid on
title('Posuv Y')
xlabel('souradnice Z [mm]'),ylabel('substep'),zlabel('posuv Y [mm]') 
legend([pleg11,pleg12,pleg13],'posuv Y','posuvy kolem ','čas otevírání','Location','northeast')
print('PosuvY','-dpng')

%vykresleni odecitani napeti
figure(2) 
plot3(Zcoordsmat,Stmat,zatmat,'r-')
hold on
plot3(Zcoordsmat3,time,OTVNAP,'ko','linewidth',2)
grid on

pleg21 = plot3(0,0,0,'r-','Visible','off');
pleg22 = plot3(0,0,0,'ko','linewidth',2,'Visible','off');

title('Otvírací napětí')
xlabel('souradnice Z [mm]'),ylabel('substep'),zlabel('napeti [MPa]')
legend([pleg21,pleg22],'zatěžující napětí','otvírací napětí','Location','northeast')
print('OtviraciNapeti6','-dpng')

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

%poměry f a u
% f = Sotviraci/Smax v danem cyklu
% u = (1-f)/(1-R)
f=zeros(pccelk-1,pocuzluZ);
for i=1:pocuzluZ
for q5=2:pccelk
    f(q5-1,i)=OTVNAP(q5-1,i)/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
end
f=[zeros(1,pocuzluZ);f];

u=zeros(pccelk,pocuzluZ);
for i=1:pocuzluZ
for q6=2:size(f,1)
    if f(q6,i)==0
        u(q6,i)=0;
    else
    u(q6,i)=(1-f(q6,i))/(1-R);
    end
end
end

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

Zcoordsmat5=repmat(Zcoords(:,2)',2,1);

figure(4)
plot3(Zcoordsmat4,repmat([1:1:size(f,1)]',1,pocuzluZ),f,'rx')
hold on
plot3(Zcoordsmat4,repmat([1:1:size(u,1)]',1,pocuzluZ),u,'b.')
hold on
plot3(Zcoordsmat5,repmat([prumod;size(f,1)],1,pocuzluZ),[prumf;prumf],'r-')
hold on
plot3(Zcoordsmat5,repmat([prumod;size(u,1)],1,pocuzluZ),[prumu;prumu],'b-')

pleg41 = plot3(0,0,0,'rx','Visible','off');
pleg42 = plot3(0,0,0,'b.','Visible','off');
pleg43 = plot3(0,0,0,'r-','Visible','off');
pleg44 = plot3(0,0,0,'b-','Visible','off');

grid on
axis([0 max(Zcoords(:,2)) 0 pccelk 0 1.1*max([max(f(:)),max(u(:))])])
title('Pomery f a u')
xlabel('souradnice Z [mm]'),ylabel('zatezovy cyklus'),zlabel('u,f [-]')
legend([pleg41,pleg42,pleg43,pleg44],'f = Sotv/Smax','u = (1-f)/(1-R)','prumer f','prumer u','Location','northeast')
print('OtviraciPomery6','-dpng')

%graf OTVNAP na tloustce
figure(5)
plot(Zcoords(:,2)',OTVNAP(size(OTVNAP,1),:),'ko')
grid on

title('Otviraci napeti pri poslednim cyklu')
xlabel('souradnice Z [mm]'),ylabel('Otviraci napeti [MPa]')
print('OtviraciNapeti6naTl','-dpng')

%graf u pomer na tloustce
figure(6)
plot(Zcoords(:,2)',prumf,'rx')
hold on
plot(Zcoords(:,2)',prumu,'bo')
grid on

title('pomery f a u')
legend('f = Sotv/Smax','u = (1-f)/(1-R)','Location','northeast')
xlabel('souradnice Z [mm]'),ylabel('u, f [-]')
print('OtviraciPomery6naTl','-dpng')


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




