kill(all)$

load("vect")$

G1:[0,0,-m1*g]$
G2:[0,0,-m2*g]$
RA:[XA,YA,ZA]$
RB:[XB,YB,0]$
Φ1:[-m1*ddx1+m1*ω^2*x1,-2*m1*ω*dx1,0]$
Φ2:[m2*ω^2*L/2,0,0]$
Mvr:[0,0,M]$

OP:[x1,0,0]$
OC:[L/2,0,0]$
OA:[0,0,-3*L/2]$
OB:[0,0,3*L/2]$

eq1:G1+G2+RA+RB+Φ1+Φ2$
eq2:express(OP~G1+OC~G2+OA~RA+OB~RB+OP~Φ1+OC~Φ2+Mvr)$

u1:eq1[1]=0;
u2:eq1[2]=0;
u3:eq1[3]=0;
u4:eq2[1]=0;
u5:eq2[2]=0;
u6:eq2[3]=0;

RES:solve([u1,u2,u3,u4,u5,u6],[XA,YA,ZA,XB,YB,M]);

NY:[0,Ny,0]$
NZ:[0,0,Nz]$
Fupr:[c*(x1-L0),0,0]$

eq3:G2+RA+RB+NY+NZ+Fupr+Φ2$

u7:eq3[1]=0;
u8:eq3[2]=0;
u9:eq3[3]=0;

RESN:solve(subst(RES,[u8,u9]),[Ny,Nz]);

eqx:ratsimp(subst(RES,u7));

eqx1:subst([x1=x1(t),dx1=diff(x1(t),t),ddx1=diff(x1(t),t,2)],eqx);

RESX:desolve(eqx1,x1(t));

zamNU:[(at('diff(x1(t),t,1),t=0))=V0,at(x1(t),t=0)=x0];

RESX2:subst(zamNU,RESX);

zampar:[m1=1,m2=3,g=9.81,c=40,L0=2,L=4,ω=0.1,V0=0,x0=1.5];

X1(t1):=at(subst(RESX2,x1(t)),t=t1)$
dX1(t1):=at(diff(subst(RESX2,x1(t)),t),t=t1)$
ddX1(t1):=at(diff(subst(RESX2,x1(t)),t,2),t=t1)$

zamx1:[x1=X1(t),dx1=dX1(t),ddx1=ddX1(t)]$

XA1:subst(zampar,subst(zamx1,subst(RES,XA)))$
YA1:subst(zampar,subst(zamx1,subst(RES,YA)))$
ZA1:subst(zampar,subst(zamx1,subst(RES,ZA)))$
XB1:subst(zampar,subst(zamx1,subst(RES,XB)))$
YB1:subst(zampar,subst(zamx1,subst(RES,YB)))$

XA2:subst(zampar,XA1*cos(ω*t)-YA1*sin(ω*t))$
YA2:subst(zampar,XA1*sin(ω*t)+YA1*cos(ω*t))$

NY1:subst(zampar,subst(zamx1,subst(RESN,Ny)))$
NZ1:subst(zampar,subst(zamx1,subst(RESN,Nz)))$

wxplot2d ([XA1,YA1,ZA1],[t,0,15], [xlabel,"t,с"],[ylabel,"XA,YA,ZA, H"],[legend,"XA","YA","ZA"]);

wxplot2d ([XB1,YB1],[t,0,15], [xlabel,"t,с"],[ylabel,"XB,YB, H"],[legend,"XB","YB"]);

wxplot2d ([XA2,YA2],[t,0,15], [xlabel,"t,с"],[ylabel,"XA,YA, H"],[legend,"XA","YA"]);

wxplot2d ([NY1,NZ1],[t,0,15], [xlabel,"t,с"],[ylabel,"Ny,Nz, H"],[legend,"Ny","Nz"]);