#  Расчет возмущающих моментов гиростабилизатора
 with(linalg):
 with(plots):; 
 with(ColorTools):;


           
Digits:=4;

#   ----- Step 1 --------------------- 
print ("----- Step 1 - Vvod dannyh -------------------- ");
g:=10.;
Moc:= 2.;Mchel:=0.5;
Mos:=0.2;Mdop:=0.1;
Mpl:= Moc +Mchel+ Mos+Mdop;
Mpl:=ceil(Mpl);

#----- Step 2 ---------------------
print ("----- Step 2 ---Ves tela ------------------ ");
Ppl:=Mpl*g;
Mnr:=ceil(0.3*Mpl);
Pnr:=Mnr*g;
Psum:=ceil(Ppl+Pnr);

#---- Step 3 ----- Reakcii R, A ----------------
print ("----- Step 3 ---- Reakcii opor: R, A ----------------- ");
Pi180:=evalf(Pi/180);
nx:=3.;ny:=3.;nz:=3.;nv:=1.;
FRnr:=ceil(Psum/2*(evalf(sqrt(nx^2+nz^2))+nv));
FAnr:=ceil(Psum*(ny+nv));
bkGr:=20;
bk:=evalf(bkGr*Pi180);
APNnr:=ceil(0.1*evalf(1.58*FRnr*tan(bk)+0.5*FAnr))*10;

FRpl:=ceil(Ppl/2*sqrt(nx^2+ny^2+nz^2)+nv);
FApl:=ceil(Ppl*(evalf(sqrt(nx^2+nz^2))+nv));
APNpl:=ceil(0.1*evalf(1.58*FRpl*tan(bk)+0.5*FApl))*10;


#---- Step 4 ----- M treniya ----------------
print ("----- Step 4 ---- M treniya opor ----------------- ");
M0:=7.2;F0:=500;k1:=0.0069;k2:=0.0035;k3:=0.008;
Mtry:=2*ceil(0.1*evalf(M0+F0*(k1-k2)+(k2*FRnr+k3*APNnr)*100))*10;
Mtrx:=2*ceil(0.1*evalf(M0+F0*(k1-k2)+(k2*FRpl+k3*APNpl)*100))*10;

#---- Step 5 ---- M  nebalance  pri sborke----------------
print ("----- Step 5 --- M  nebalance  pri sborke ---------- ");
FRnr0:=ceil(Psum/2);
FAnr0:=0; 
APNnr0:=0;
FRpl0:=ceil(Ppl/2);
FApl0:=0;
APNpl0:=0;
M0:=7.2;F0:=500.;k1:=0.0069;k2:=0.0035;k3:=0.0008;
M00:=M0+F0*(k1-k2);

Mnby0:=2*ceil(M00+(k2*FRnr0)*100);
Mnbx0:=2*ceil(M00+(k2*FRpl0)*100);

#---- Step 6 ---- M  nebalance  pri ekspluatacii----------------
print ("----- Step 6 -- M  nebalance  pri ekspluatacii --------- ");
aM:=70;
bM:=60;
MnbYmin:=10000;
MnbYmax:=-10000;
MnbXmin:=10000;
MnbXmax:=-10000;
# ----Cycle--2----
for na from -aM by 1 to aM do
alf:=na*1.0*Pi180;ca:=cos(alf);sa:=sin(alf);
          for nb from -bM by 1 to bM do
          bet:=nb*1.0*Pi180; cb:=cos(bet);sb:=sin(bet);

          MnbY:=evalf(Mnby0*(nx*(sa-ca)+nz*(sa+ca))+Mnbx0*cb*(nx*(sa-ca*cb)+ny*sb+nz*(sa*cb+ca)));

          MnbX:=evalf( Mnbx0*((nx*ca-nz*sa)*(cb+sb)+ny*(sb-cb)));


          if MnbY < MnbYmin then MnbYmin:=MnbY:end if;
          if MnbX < MnbXmin then MnbXmin:=MnbX:end if;

          if MnbY > MnbYmax then MnbYmax:=MnbY:end if;
          if MnbX > MnbXmax then MnbXmax:=MnbX:end if;
          od:
od:

print ( "MY nb =",ceil(MnbYmin),ceil(MnbYmax));
print ( "MX nb =",ceil(MnbXmin),ceil(MnbXmax));

MnbY:=ceil(max(abs(MnbYmin),abs(MnbYmax)));

MnbX:=ceil(max(abs(MnbXmin),abs(MnbXmax)));

#----- Step 7 ----inerc M  ot vrashenia LA ----------------
print ("----- Step 7 --- inerc M  ot vrashenia LA ------------------ ");

# ocenka momentov inerzii platformy
apl := 0.220; bpl := .220; cpl := 0.8e-1;
#JplX :=Mpl*(bpl^2+cpl^2)/12;
#JplY := Mpl*(apl^2+cpl^2)/12;
JplZ := Mpl*(apl^2+bpl^2)/12;

# ocenka momentov inerzii ramy
 anr := .240; bnr := .240; cnr := 0.2e-1;
t:=0.02;
Vnr1:=anr*bnr*cnr;
Vnr2:=(anr-t)*(bnr-t)*cnr;
#Vnr:=anr*bnr*cnr-(anr-t)*(bnr-t)*(cnr-t);
Vnr:=Vnr1-Vnr2;
plotn:=Mnr/Vnr;

 JnrX:=plotn*(Vnr1*(bnr^2+cnr^2) -Vnr2*( (bnr-t)^2+cnr^2) )/12;
 JnrY:=plotn*(Vnr1*(anr^2+cnr^2) -Vnr2*( (anr-t)^2+cnr^2) )/12;
 JnrZ:=plotn*(Vnr1*(bnr^2+anr^2) -Vnr2*( (bnr-t)^2+(anr-t)^2) )/12;

#vliyanie uglovyx ckorostei vrasheniya LA vokrug prodol'noi osi X i poperechnoi osi Z  

B:=JnrY+JplZ;


w100:=100*Pi180;
wxc:=w100;wyc:=w100;wzc:=w100;
wt140:=140*Pi180;
wtxc:=wt140;wtyc:=wt140;wtzc:=wt140;

wwMYmin:=10000;
wwMYmax:=-10000;

# ----Cycle--2----

for nb from -bM by 1 to bM do
bet:=evalf(nb*Pi180); cb:=cos(bet);sb:=sin(bet);tng:=evalf(tan(bet));
A:=( JnrX+JnrY/cb^2+(B+JplZ)*tng^2-JnrZ);

wwinMY:=evalf(B*tng*(wtxc*ca-wtzc*sa+wyc*(wxc*sa+wzc*ca))+A*((wxc^2-wzc^2)*s2a +2*wxc*wzc*c2a)/2);
          for na from -aM by 1 to aM do
          alf:=evalf(na*Pi180);
          sa:=evalf(sin(alf));
          ca:=evalf(cos(alf));
          s2a:=evalf(sin(2*alf));
          c2a:=evalf(cos(2*alf));


          if wwinMY < wwMYmin then wwMYmin:=wwinMY:end if;
          if wwinMY > wwMYmax then wwMYmax:=wwinMY:end if;

          od:
od:
print ( "wwMYmin =",ceil(wwMYmin*10^4));
print ( "wwMYmax =",ceil(wwMYmax*10^4));
wwinMY:=ceil(max(abs(wwMYmin),abs(wwMYmax))*10^4);


#----- Step 8 ----inerc Moment  ot uglovyh kolebanii LA -----------
print ("----- Step 8 -- inerc Moment  ot uglovyh kolebanii LA ------ ");
gamGr:=2;
gam:=gamGr*Pi180;

f:=3;
wf:=2*Pi*f;

gamMYmin:=10000;
gamMYmax:=-10000;

nt:=30;ht:=2*Pi/(f*nt); #nt:=30;ht:=2*Pi/(f*nt);

# ----Cycle--3----
for nb from -bM by 1 to bM do
bet:=evalf(nb*Pi180); cb:=cos(bet);sb:=sin(bet);tng:=evalf(tan(bet));
A:=evalf( JnrX+JnrY/cb^2+(B+JplZ)*tng^2-JnrZ);

       for na from -aM by 1 to aM do
       alf:=evalf(na*Pi180);
       sa:=evalf(sin(alf));
       ca:=evalf(cos(alf));
       s2a:=evalf(sin(2*alf));
       c2a:=evalf(cos(2*alf));

       t:=0;
             for it from 0 to 2*nt do
             wxc:=-gam*wf*cos(wf*t);
             wtxc:=-gam*wf^2*sin(wf*t);
             wyc:=wxc; wzc:=wxc;
             wtyc:=wtxc;wtzc:=wtxc;

             gamMY:=evalf(B*tng*(wtxc*ca-wtzc*sa+wyc*(wxc*sa+wzc*ca)+A*((wxc^2-wzc^2)*s2a +2*wxc*wzc*c2a)/2));

             if gamMY < wwMYmin then gamMYmin:=gamMY:end if;
             if gamMY > wwMYmax then gamMYmax:=gamMY:end if;
             t:=t+ht;
             od;
        od:
od:
print ( "gamMYmin =",ceil(gamMYmin*10^4));
print ( "gamMYmax =",ceil(gamMYmax*10^4));
gamMY:=ceil(max(abs(gamMYmin),abs(gamMYmax))*10^4);


#----- Step 9 ---- Mynr u  neravnozhestkogo podvesa -----------
print ("----- Step 9 -- Mynr u  neravnozhestkogo podvesa ------ ");
NBradHP:=250;
NBaxeHP:=100;
NBradPL:=270;
NBaxePL:=110;
CXHP:=1.45;
CYHP:=1.58;
CZHP:=1.39;
CqXHP:=7.80;
CqYHP:=21.6;
CqZHP:=10.4;
cHP0X:=1/(1/CXHP +1/CqXHP+ 1/NBradHP);
cHP0Z:=1/(1/CZHP +1/CqZHP+ 1/NBradHP);
cPL0Z:=NBradPL;
cPL0Y:=NBradPL;
McHPy:=ceil(evalf((nx^2+nz^2)*(Pnr^2*abs(1/cHP0X -1/cHP0Z)+
Pnr*Ppl*(1/cHP0X -1/cPL0Z)+Ppl^2/cPL0Z))/4/100);
McPLx:=ceil(evalf( (nx^2+nz^2)*Ppl^2*abs(1/cPL0Y -1/cPL0Z)) /4/100);
lambdaHP := 3;
lambdaPL := 6;
McvHPy := ceil(evalf(nv^2*lambdaHP*lambdaPL*
(Pnr^2*abs(1/cHP0X-1/cHP0Z)+Pnr*Ppl*(1/cHP0X-1/cPL0Z)+Ppl^2/cPL0Z))/4/100);


#----- Step 10 ---- Momenty v tokovyvodovah -----------
print ("----- Step 10 --- Momenty v tokovyvodovah ------------------ ");
Nknr:=58;
Nkpl:=50;
Nsh:=1;
Rk:=0.5;
Fk:=0.1;f:=0.2;
Mknr := 100*f*Fk*Rk*Nknr*Nsh;
Mkpl := 100*f*Fk*Rk*Nkpl*Nsh;



tetaGrnr := 70; 
tetanr := tetaGrnr*Pi180; 
tetaGrpl := 60; 
tetapl := tetaGrpl*Pi180;
 LtpHP := 12;
 LtpPL := 10;
 MtrwHPy := ceil(6.5*Nknr*tetanr/LtpHP);
 MtrwPLx := ceil(6.5*Nknr*tetapl/LtpPL);
 Mtoky:= min(MtrwHPy,Mknr);
 Mtokx:= min(MtrwPLx,Mkpl);

#----- Step 11 ---- Momenty privoda -----------
print ("----- Step 11 --- Momenty privoda ------------------- ");
MinPy:= 80;
MdPy:=280;

MinPx:= 23;
MdPx:=60;

#----- Step 12 ---- Summarnye Max Momenty  -----------
print ("----- Step 12 ---- Summarnye Max Momenty ------------- ");

SumMy:=Mtry+MnbY+wwinMY+gamMY+McHPy+McvHPy+Mtoky + MinPy+MdPy;
SumMx:=Mtrx+MnbX+McPLx+Mtokx + MinPx+MdPx;