%krivky
clc;
close all;
clear;

tic;

m=2.5;


nas=input('Do you want 1 or 2 peaks ? ');

type=input('Do you want plot dependecy to the time interval y=1, n=2 ');

nod=1000;

if type==2
    tused=input('How many ms the impact last? ');
    T=linspace(0,tused*10^(-3),nod);
    lT=length(T);
    a=zeros(lT,1);    
    
    if nas==1
        a(:,1)=1;
    elseif nas==2
        a(:,1)=1;
        a(round(nod*3/7):round(nod*4/7))=0;
    end
                 
        
             
        enx=-inf;
        numofiter=5000;
        for i=1:numofiter
            A=a*i;   %% units [m/s^2]
            v0=trapz(T,A);
            v=cumtrapz(T,A)-v0;   %% units [m/s]
            sx=cumtrapz(T,v);
            
            
            q=trapz(abs(sx),A*m);
            
            if q < 138.5;
                enx = q;
                koef=i;
            else
                break
            end
        end
        
        
        
        MA=abs(koef*a);
        v0=trapz(T,MA);
        VA=cumtrapz(T,MA)-v0;
        SA=cumtrapz(T,VA);
        ENmy=trapz(abs(SA),MA*m);
        [hic,time,t1,t2]=HICcalc(T,MA./9.81);
        hic;
        
elseif type==1
    
    
    
    timebound=[5 30];
    timelin=linspace(timebound(1),timebound(2),nod);
    
    for t=1:nod
        T=linspace(0,timelin(t)*10^(-3),nod);
        lT=length(T);
        a=zeros(lT,1);
        
        
        if nas==1
            a(:,1)=1;
        elseif nas==2
            
            a(:,1)=1;
            a(round(nod*2/5):round(nod*3/5))=0;
        end
        
        enx=-inf;
        numofiter=5000;
        for i= 1:numofiter
            
            A=i*a;
            v0=trapz(T,A);
            v=cumtrapz(T,A)-v0;   %% units [m/s]
            sx=cumtrapz(T,v);
            
            q=trapz(abs(sx),A*m);
            
            if q < 138.5;
                enx = q;
                koef=i;
            else
                break
            end
        end
        
        
        
        MA=abs(koef*a);
        v0=trapz(T,MA);
        VA=cumtrapz(T,MA)-v0;
        SA=cumtrapz(T,VA);
        ENmy=trapz(abs(SA),MA*m);
        [hic,time,t1,t2]=HICcalc(T,MA./9.81);
        hic;
        
        B(t,1)=(timelin(t));
        B(t,2)=(hic);
        B(t,3)=(t1)*1000;
        B(t,4)=(t2)*1000;
        B(t,5)=(time)*1000;
        B(t,6)=(max(MA));
        B(t,7)=(max(abs(VA)));
        B(t,8)=(max(abs(SA)))*1000;
        
        
    end
        
    
      if nas==1
          peak1=B;
          save('HIC1peak.mat','peak1')
          
      elseif nas==2
          peak2=B;
          save('HIC2peak.mat','peak2')
      end
      
end
    

    
    %% PLOTTING
    if type==1
    peak1l=load('HIC1peak.mat','peak1');
    peak2l=load('HIC2peak.mat','peak2');
    
    figure
    plot(peak1l.peak1(:,1),peak1l.peak1(:,2),peak2l.peak2(:,1),peak2l.peak2(:,2))
    title('HIC value vs. time interval T for rectangular function')
    xlabel('Time interval T [s]')
    ylabel('HIC value [-]')
    legend('1 peak', '2 peak')
    
    
    figure
    plot(peak1l.peak1(:,8),peak1l.peak1(:,2),peak2l.peak2(:,8),peak2l.peak2(:,2))
    
   
    
    figure
    plot(B(:,1),B(:,2),'b')
    title('Time T vs. HIC diagram '  )
    xlabel('Time interval T [s]')
    ylabel('HIC value [-]')
    legend([ num2str(nas) ' peak'])
    end
    
    if type==2
    figure
    plot(T,abs(MA),'b')
    title(['Deceleration A vs. time t diagram for time interval T = ' num2str(tused) 'ms' ] )
    xlabel('Time t [s]')
    ylabel('Deceleration A [m/s^-2]')
    legend(['HIC = ' num2str(round(hic)) ' s = ' num2str(abs(SA(end)*1000)) ' mm'],4)
    
    k1 = find(T==t1);
    k2 = find(T==t2);
    
    hold on
    plot(t1,MA(k1),'--rs','LineWidth',2,...
        'MarkerEdgeColor','k',...
        'MarkerFaceColor','g',...
        'MarkerSize',10);
    plot(t2,MA(k2),'--rs','LineWidth',2,...
        'MarkerEdgeColor','k',...
        'MarkerFaceColor','g',...
        'MarkerSize',10);
    hold off
    end
    
    figure
    plot(T,abs(VA),'g')
    title('Velocity V vs. time t diagram')
    xlabel('Time t [s]')
    ylabel('Velocity [m/s]')
    
    figure
    plot(T,abs(SA),'r')
    title('Trajectory S vs. time t diagram')
    xlabel('Time t [s]')
    ylabel('Trajectory S [m]')
    
    
    %
    %
    %
    %
    % if type==1
    %     load('traj60.mat');
    %     load('traj65.mat');
    %     load('traj70.mat');
    % figure
    % plot(S1(:,1),S1(:,2),S2(:,1),S2(:,2),S3(:,1),S3(:,2));
    % legend('s = 60 mm','s = 65 mm','s = 70 mm',2)
    % xlabel('Time t [ms]')
    % ylabel('HIC value [-]')
    % else
    % end
    %
    %
    % figure
    % p = plot(sx,MA);
    % legend('Sin (x)',2)
    %  title('Force F vs. trajectory s')
    %  xlabel('Trajectory s [m]')
    %  ylabel('Force F [N]')
    %  set(p,'Color','blue','LineWidth',2)
    %
    %
    %  figure
    % p = plot(T,sx);
    % legend('Sin (x)',2)
    %  title('Trajectory s vs. time t' )
    %  xlabel('Time t [s]')
    %  ylabel('Trajectory s [m]')
    %  set(p,'Color','green','LineWidth',2)
    %
    %
    %  k1 = find(T==t1);
    %  k2 = find(T==t2);
    %
    %  figure
    % p=plot(T,MA);
    % legend('Sin (x)',2)
    %  title(['Force F vs. time t HIC = ' num2str(hic)])
    %  xlabel('Time t [s]')
    %  ylabel('Force F [N]')
    %  set(p,'Color','red','LineWidth',2)
    %
    %  hold on
    % plot(t1,MA(k1),'--rs','LineWidth',2,...
    %                 'MarkerEdgeColor','k',...
    %                 'MarkerFaceColor','g',...
    %                 'MarkerSize',10);
    % plot(t2,MA(k2),'--rs','LineWidth',2,...
    %                 'MarkerEdgeColor','k',...
    %                 'MarkerFaceColor','g',...
    %                 'MarkerSize',10);
    %
    %
    % time=time*1000
    %
    % krivka(:,1)=sx(:,1);
    % krivka(:,2)=MA(:,1);
    %
    %
    % casvypoctu=toc;
    %
    %
    % hic
    %
