clc
clear all
%% data
% pricny prurezy
Sv1=0.049593; % m^2
Ss=0.01;
Sc=0.02055;
Sv2=0.05;
Smc= 0.05;
Smv=0.05;
% rod lengths
Lmc=6.66;
Lmv=5;
Lv=4.4;
Ls=4.72;
Lc=6.4528;
E=210e9; % [Pa] young modulus
g=9.81 ; %gravitacní zrychleni
Ro=7850; %hustota
al= 47.01; %angle betwen the tilted and vertical rods
%% train forces
Mk=123e3; %[Kg] tokomotiva
Mv=63e3; %[Kg] Vagon
F_lok=(Mk*g)/(2*6);
F_vag=(Mv*g)/(2*4);
%% Forces and important matrixes
syms L1 L2 L3 L4 L5 L6 L7 L8 L9 L10 L12 L13 L14 L15 L16 L17 L18 L19 ...
    L11 L20 L22 L24 L26 L28 L30 L32 L34 L36 L21 L23 L25 L27 L29 L31 L33 L35 L37
[L1 L2 L3 L4 L5 L6 L7 L8 L9 L10 L12 L13 L14 L15 L16 L17 L18 L19 ]=deal(Lv,Lv,Lv,Lv,Lv,Lv,Lv,Lv,Lv,Lv,Lv,...
    Lv,Lv,Lv,Lv,Lv,Lv,Lv);
[L11 L20 L22 L24 L26 L28 L30 L32 L34 L36]=deal(Lc,Lc,Lc,Lc,Lc,Lc,Lc,Lc,Lc,Lc);
[L21 L23 L25 L27 L29 L31 L33 L35 L37]=deal(Ls,Ls,Ls,Ls,Ls,Ls,Ls,Ls,Ls);
Length_matrix=[L1 L2 L3 L4 L5 L6 L7 L8 L9 L10 L11 L12 L13 L14 L15 L16 L17 L18 L19 ...
    L20 L21 L22 L23 L24 L25 L26 L27 L28 L29 L30 L31 L32 L33 L34 L35 L36 L37 ];
syms S1 S2 S3 S4 S5 S6 S7 S8 S9 S10 S12 S13 S14 S15 S16 S17 S18 S19 ...
    S11 S20 S22 S24 S26 S28 S30 S32 S34 S36 S21 S23 S25 S27 S29 S31 S33 S35 S37
[S1 S2 S3 S4 S5 S6 S7 S8 S9 S10]=deal(Sv2,Sv2,Sv2,Sv2,Sv2,Sv2,Sv2,Sv2,Sv2,Sv2);
[S12 S13 S14 S15 S16 S17 S18 S19]=deal(Sv1,Sv1,Sv1,Sv1,Sv1,Sv1,Sv1,Sv1);
[S11 S20 S22 S24 S26 S28 S30 S32 S34 S36]=deal(Sc,Sc,Sc,Sc,Sc,Sc,Sc,Sc,Sc,Sc);
[S21 S23 S25 S27 S29 S31 S33 S35 S37] =deal(Ss,Ss,Ss,Ss,Ss,Ss,Ss,Ss,Ss);
CrossSection_matrix=[S1 S2 S3 S4 S5 S6 S7 S8 S9 S10 S11 S12 S13 S14 S15 S16 S17 S18 ...
    S19 S20 S21 S22 S23 S24 S25 S26 S27 S28 S29 S30 S31 S32 S33 S34 S35 S36 S37];
CrossSection_matrix1=CrossSection_matrix';

Stress_y= 210; % Mpa

Fb=0.5*g*Ro*(2*Sv2*Lv +Ss*Ls+ Smc*Lmc+Smv*Lmv); 
Fb=round(Fb);
 Fd =round(Fb+(5/3)*F_vag);
 Ff =round(Fb+F_vag);
 Fh =round(Fb+((7/6)*F_lok));
 Fj =round(Fb+((4/3)*F_lok));
Fb1=round(Fb+(1/3)*F_vag);
Fc=0.5*g*Ro*(2*Sv2*Lv +Ss*Ls +2*Sc*Lc*sind(al)+Smc*Lmc+Smv*Lmv);
Fc=round(Fc);
 Fe = round(Fc+(1/3)*F_vag);
 Fg = round(Fc+((1/3)*F_lok)+F_vag);
 Feye = round(Fc+((11/6)*F_lok));
Fc1=round(Fc+(5/3)*F_vag);
Fa=0.5*g*Ro*(Sv2*Lv +Sc*Lc*sind(al) + Smv*Lmv );
Fa=round(Fa);
Fk=0.5*g*Ro*(Sv2*Lv +Sc*Lc*sind(al)+ Smc*Lmc + Smv*Lmv );
Fk=round(Fk);
Fk=round(Fk+((4/3)*F_lok));
Fl=0.5*g*Ro*(Sv1*Lv +Ss*Ls +2*Sc*Lc*sind(al));
Fl=round(Fl);
 Ft = Fl ;
Fs=0.5*g*Ro*(2*Sv1*Lv +Ss*Ls);
Fs=round(Fs);
 Fq = Fs;
 Fo = Fs;
 Fm = Fs;
Fr=0.5*g*Ro*(2*Sv1*Lv +Ss*Ls +2*Sc*Lc*sind(al));
Fr=round(Fr);
 Fp = Fr;
 Fn = Fr;

    F_matrix=[Fa ;Fb1 ;Fc1 ;Fd ;Fe ;Ff ;Fg ;Fh ;Feye ;Fj ;Fk ;Fl ;Fm ;Fn ;Fo ;Fp ;Fq ;Fr ;Fs ;Ft];
format long
syms N1 N2 N3 N4 N5 N6 N7 N8 N9 N10 N11 N12 N13 N14 N15 N16 N17 N18 N19 N20...
    N21 N22 N23 N24 N25 N26 N27 N28 N29 N30 N31 N32 N33 N34 N35 N36 N37 Fy1 Fy2% Fa Fb Fc Fd Fe Ff Fg Fh Feye Fj Fk Fl Fm Fn Fo Fp Fq Fr Fs Ft
%% equations
eqn2 =Fy1 + N20*sind(al)-Fa==0;
eqn1 = N1+N20*cosd(al)==0;
eqn3 = N2-N1==0;
eqn4 = N21-Fb1==0;
eqn8 = N25-Fd==0;
eqn12 = N29-Ff==0;
eqn16 = N33-Fh==0;
eqn20 = N37-Fj==0;
eqn22 = Fy2 +N11*sind(al)-Fk==0;
eqn21 = -N10-N11*cosd(al)==0;
eqn26 = -N35-Fm==0;
eqn30 = -N31-Fo==0;
eqn34 = -N27-Fq==0;
eqn38 = -N23-Fs==0;
eqn40 = -N21-N22*sind(al)-N20*sind(al)-Ft==0;
eqn6 = N23-Fc1+N24*sind(al)+N22*sind(al)==0;
eqn36 = -N25-N26*sind(al)-N24*sind(al)-Fr==0;
eqn10 = N27-Fe+N28*sind(al)+N26*sind(al)==0;
eqn32 = -N29-N28*sind(al)-N30*sind(al)-Fp==0;
eqn14 = N31-Fg+N32*sind(al)+N30*sind(al)==0;
eqn28 = -N33-N34*sind(al)-N32*sind(al)-Fn==0;
eqn18 = N35-Feye+N36*sind(al)+N34*sind(al)==0;
eqn39 = N22*cosd(al)-N20*cosd(al)+N19==0;
eqn37 = N18-N19==0;
eqn23 = N11*cosd(al)-N36*cosd(al)-N12==0;
eqn25 = N12-N13==0;
eqn5 = N3-N2+N24*cosd(al)-N22*cosd(al)==0;
eqn7 = N4-N3==0;
eqn9 = N5-N4+N28*cosd(al)-N26*cosd(al)==0;
eqn11 = N6-N5==0;
eqn13 = N7-N6+N32*cosd(al)-N30*cosd(al)==0;
eqn15 = N8-N7==0;
eqn17 = N9-N8+N36*cosd(al)-N34*cosd(al)==0;
eqn27 = N34*cosd(al)-N32*cosd(al)-N14+N13==0;
eqn29 = N14-N15==0;
eqn31 = N30*cosd(al)-N28*cosd(al)-N16+N15==0;
eqn33 = N16-N17==0;
eqn24 = -N37-N11*sind(al)-N36*sind(al)-Fl== 0;
eqn35 = N17-N18+N26*cosd(al)-N24*cosd(al) == 0;
eqn19=N10-N9==0;
%% solution
solu=solve(eqn2,eqn1,eqn3,eqn4,eqn8,eqn12,eqn16,eqn20,eqn22,eqn21,eqn26,...
    eqn30,eqn34,eqn38,eqn40,eqn39,eqn37,eqn6,eqn23,eqn25,eqn18,eqn28,...
    eqn14,eqn32,eqn10,eqn36,eqn5,eqn7,eqn9,eqn11,eqn13,eqn15,eqn17,eqn27,...
    eqn29,eqn31,eqn33,eqn24,eqn35,eqn19);

%% normal forces
N20=vpa(solu.N20);N1=vpa(solu.N1);N2=vpa(solu.N2);N21=vpa(solu.N21);N25=vpa(solu.N25);
N29=vpa(solu.N29);N33=vpa(solu.N33);N37=vpa(solu.N37);N11=vpa(solu.N11); N10=vpa(solu.N10);
N35=vpa(solu.N35);N31=vpa(solu.N31);N27=vpa(solu.N27);N28=vpa(solu.N28);N22=vpa(solu.N22);
N19=vpa(solu.N19);N18=vpa(solu.N18);N36=vpa(solu.N36);N12=vpa(solu.N12);N13=vpa(solu.N13);
N34=vpa(solu.N34);N32=vpa(solu.N32);N30=vpa(solu.N30);N26=vpa(solu.N26);N24=vpa(solu.N24);
N23=vpa(solu.N23);N3=vpa(solu.N3);N4=vpa(solu.N4);N5=vpa(solu.N5);N6=vpa(solu.N6);
N7=vpa(solu.N7);N8=vpa(solu.N8);N9=vpa(solu.N9);N14=vpa(solu.N14);N15=vpa(solu.N15);
N16=vpa(solu.N16);N17=vpa(solu.N17);Fy1=vpa(solu.Fy1);Fy2=vpa(solu.Fy2);


N=[N1; N2 ;N3 ;N4; N5; N6; N7 ;N8 ;N9 ;N10; N11; N12; N13; N14; N15; N16; N17; N18; N19; N20; N21; N22; N23;...
    N24; N25; N26; N27; N28; N29; N30; N31; N32; N33; N34; N35; N36; N37];
%% stresses
stressV1= (1e-6/Sv1)*[N1 N2 N3 N4 N5 N6 N7 N8 N9 N10];
stressV2= (1e-6/Sv2)*[N12 N13 N14 N15 N16 N17 N18 N19];
stressS= (1e-6/Ss)*[N21 N23 N25 N27 N29 N31 N33 N35 N37];
stressC= (1e-6/Sc)*[N20 N22 N24 N26 N28 N30 N32 N34 N36 N11];
stress = [stressV1 stressV2 stressS stressC];
stress_1=abs(stress);
stress_max=max(stress_1)
Kk=Stress_y/stress_max
%% maximum diflection
[Diff1,Diff2,Diff3,Diff4,Diff5,Diff6,Diff7,Diff8,Diff9,Diff10,Diff11]=deformation1;
 w=0;
Wa=0; Wb=0; Wc=0; Wd=0; We=0; Wf=0; Wg=0; Wh=0; Wi=0; Wj=0; Wk=0; 
for i=1:size(N) 
    Wa=Wa+(((N(i)*Length_matrix(i))/(E*CrossSection_matrix(i)))*Diff1(i));
    Wb=Wb+(((N(i)*Length_matrix(i))/(E*CrossSection_matrix(i)))*Diff2(i));
    Wc=Wc+(((N(i)*Length_matrix(i))/(E*CrossSection_matrix(i)))*Diff3(i));
    Wd=Wd+(((N(i)*Length_matrix(i))/(E*CrossSection_matrix(i)))*Diff4(i));
    We=We+(((N(i)*Length_matrix(i))/(E*CrossSection_matrix(i)))*Diff5(i));
    Wf=Wf+(((N(i)*Length_matrix(i))/(E*CrossSection_matrix(i)))*Diff6(i));
    Wg=Wg+(((N(i)*Length_matrix(i))/(E*CrossSection_matrix(i)))*Diff7(i));
    Wh=Wh+(((N(i)*Length_matrix(i))/(E*CrossSection_matrix(i)))*Diff8(i));
    Wi=Wi+(((N(i)*Length_matrix(i))/(E*CrossSection_matrix(i)))*Diff9(i));
    Wj=Wj+(((N(i)*Length_matrix(i))/(E*CrossSection_matrix(i)))*Diff10(i));
    Wk=Wk+(((N(i)*Length_matrix(i))/(E*CrossSection_matrix(i)))*Diff11(i));
end
W=1000*[Wa;Wb;Wc;Wd;We;Wf;Wg;Wh;Wi;Wj;Wk]
