clc;
close all;
clear;

spot=input('Which spot on bonnet would you like to investigate? (use simple quotation marks) ');
label=['SK461_SPX_A12_Sc' spot '_ch101G_rt_A.crv'];

%% IMPORT DATA

%Acceleration
[MA]=importdata(label,',',5);
[T]=importdata(label,',',5);
[XACC]=importdata(label,',',1457);
[YACC]=importdata(label,',',1663);
[ZACC]=importdata(label,',',1869);
[XRACC]=importdata(label,',',2075);
[YRACC]=importdata(label,',',2281);
[ZRACC]=importdata(label,',',2487);

%HIC
[HIC]=importdata(label,',',211);

%Displacement and angles
    %GLOBAL COORDINATE SYSTEM
% [XCOR]=importdata(label,',',221);
% [YCOR]=importdata(label,',',427);
% [ZCOR]=importdata(label,',',633);
[XANG]=importdata(label,',',3311);
[YANG]=importdata(label,',',3517);
[ZANG]=importdata(label,',',3723);
[XD]=importdata(label,',',3929);
[YD]=importdata(label,',',4135);
[ZD]=importdata(label,',',4341);
    %MAGNITUDE OF DISPLACEMENT
[MD]=importdata(label,',',4547);

%ENERGY
    %GLOBAL
[KINen]=importdata(label,',',5372);
[INTen]=importdata(label,',',5578);
[TOTen]=importdata(label,',',5784);
[HOUen]=importdata(label,',',5990);
[CONen]=importdata(label,',',6402);


T=MA.data(:,1)/1000; %in seconds
MAMETA=MA.data(:,2)*1000; 
MA=MA.data(:,2)*1000/9.81; %in multiples of 'g' constant
XACC=(XACC.data(:,2)*1000);
YACC=(YACC.data(:,2)*1000);
ZACC=(ZACC.data(:,2)*1000);
XRACC=(XRACC.data(:,2));
YRACC=(YRACC.data(:,2));
ZRACC=(ZRACC.data(:,2));
Dx=(XD.data(:,2)/1000); %in meter
Dy=(YD.data(:,2)/1000); %in meter
Dz=(ZD.data(:,2)/1000); %in meter
% DGx=XCOR.data(:,2)/1000; %in meter
% DGy=YCOR.data(:,2)/1000; %in meter
% DGz=ZCOR.data(:,2)/1000; %in meter
% DGx=DGx-DGx(1);
% DGy=DGy-DGy(1);
% DGz=DGz-DGz(1);
Dm=MD.data(:,2)/1000; %in meter

XANG=(XANG.data(:,2)); %v rad
YANG=(YANG.data(:,2)); %v rad
ZANG=(ZANG.data(:,2)); %v rad
KINen=KINen.data(:,2);
INTen=INTen.data(:,2);
TOTen=TOTen.data(:,2);
HOUen=HOUen.data(:,2);
CONen=CONen.data(:,2);

m=2.5; %kg
I=0.0036; %kgm^2
Fx=m.*XACC; %force by Newton's law in X direction
Fy=m.*YACC; %force by Newton's law in Y direction
Fz=m.*ZACC; %force by Newton's law in Z direction
F=m.*MAMETA;
% Mkx=I.*XANG; %kravina!!!!
% Mky=I.*YANG;
% Mkz=I.*ZANG;
Mkx=I.*XRACC;
Mky=I.*YRACC;
Mkz=I.*ZRACC;


[hic,time,t1,t2]=HICcalc(T,MA);

format short g

hicmeta=HIC.data(2,2);
timemeta=(HIC.data(3,1)-HIC.data(1,1)); %in ms

hic;
timems=time*1000; %in ms




 figure(1)
 plot(T,Dx,T,Dy,T,Dz,T,0)
 legend('Dx','Dy','Dz',3)
 xlabel('Time [s]')
 ylabel('Displacement [m]')
 xlim([0 0.02])
 title('Displacements vs. time')
 
 
 
 figure(2)
 plot(T,Fx,T,Fy,T,Fz,T,0)
 legend('Fx','Fy','Fz',4)
 xlabel('Time [s]')
 ylabel('Force [N]')
 xlim([0 0.02])
 title('Forces vs. time')
 
 
 figure(3)
 p = plot(Dx,Fx,Dx,0);
 title('Relation between force Fx and trajectory Dx')
 xlabel('Trajectory Dx [m]')
 ylabel('Force Fx [N]')
 set(p,'Color','blue','LineWidth',2)
 
 figure(4)
 p = plot(Dy,Fy,Dy,0);
 title('Relation between force Fy and trajectory Dy')
 xlabel('Trajectory Dy [m]')
 ylabel('Force Fy [N]')
 set(p,'Color','green','LineWidth',2)
 
 figure(5)
 p = plot(Dz,Fz,Dz,0);
 title('Relation between force Fz and trajectory Dz')
 xlabel('Trajectory Dz [m]')
 ylabel('Force Fz [N]')
 set(p,'Color','red','LineWidth',2)
 
 figure(6)
 plot(T,XANG*180/pi,T,YANG*180/pi,T,ZANG*180/pi,T,0)
 legend('angle x','angle y','angle z',3)
 xlabel('Time [s]')
 ylabel('Angles [deg]')
 title('Angle change during impact')
 xlim([0 0.02])
 
 figure(7)
 plot(T,Mkx,T,Mky,T,Mkz,T,0)
 legend('Mkx','Mky','Mkz',3)
 xlabel('Time [s]')
 ylabel('Moment [Nm]')
 title('Moment vs. time')
 
enx=trapz(Dx,Fx);
eny=trapz(Dy,Fy);
enz=trapz(Dz,Fz);
 
 en=abs(enx+eny+enz);
 
 
 enrx=trapz(XANG,Mkx);
 enry=trapz(YANG,Mky);
 enrz=trapz(ZANG,Mkz);
 enr=enrx+enry+enrz;
 %en1=trapz(s,F); %pocet presne podle METY
 
 figure(8)
 plot(T,MA)
 xlabel('Time [s]')
 ylabel('Magnitude of acceleration [g*m/s^2]')
 title('Relation between time T and magnitude of acceleration MA')
 legend(['HIC = ' num2str(round(hic))], 1)
 xlim([0 0.02])
    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);

time=(t2-t1)*1000 ; 
 
%  figure(9)
%  plot(T,Dm)
%  xlabel('Time [s]')
%  ylabel('Magnitude of displacement [g*m/s^2]')
%  title('Relation between time T and magnitude of displacement Dm')
%  xlim([0 0.02])
 
%  figure(10)
%  plot(Dm,MA)
%  xlabel('Magnitude of displacement [g*m/s^2]')
%  ylabel('Magnitude of acceleration [g*m/s^2]')
%  title('Relation between magn. of displacement Dm and magn. of acc. MA')
%  xlim([0 0.02])
%  
 
  figure(11)
 plot(T,KINen,T,INTen,T,TOTen,T,CONen)
 legend('Kinetic energy','Internal energy','Total energy','Contact energy',3)
 xlabel('Time [s]')
 ylabel('Energy [J]')
 title('Global energy')
 xlim([0 0.02])
 ylim([0 180])
 
 
 totalenergy=KINen(end)+INTen(end)+CONen(end);
 
 energyINTCON=INTen(end)+CONen(end)
 
 energyINT=INTen(end)
 
 en
 
