clc;

% вариант
var=4;

M=[800 900 1000 1100 1200 1300 1400 1500];
m=M(var);
%m=1000;
Ppit	=2200 * 10000;   % P пит
c_sh    =1750 * 100;
A2      =0.87 / 100;
w       =20;
F0      =10000;
ymax    =3.4 / 100;
c_k     =2*10^4 * 100;
ro      =0.0008 * 1000000;
E       =0.7*10^5* 10000;

% раскомментить, если считать пример с листка
%Ppit	=2000 * 10000;   % P пит
%c_sh    =1770 * 100;

fprintf('m = %.0f kg\n', m);
disp('Поуось эллипса нагрузки по оси сил');
Fm=abs(A2*(c_sh - m*w^2));
fprintf('Fm = A2( Сш - mw^2 ) = %.1f H\n\n', Fm);

disp('Координата силы нагружения, которая соответствует максимальному значению мощности нагрузки');
Fk=3/4*F0 + sqrt( (3/4*F0)^2 - 1/2*(F0^2-Fm^2) );
%Fk=10330;
disp( 'Fk=3/4*F0 + sqrt( (3/4*F0)^2 - 1/2*(F0^2-Fm^2) )' );
fprintf('Fk = %.0f H\n', Fk);
disp(' ');

disp('Площадь поршня силового гидропривода');
Ap=3/2*Fk/Ppit;
fprintf('Ap >= 3/2*Fk/Ppit = %.2f см^2\n', Ap*10000);
%Ap=8/10000;
Ap=ceil(Ap*10000)/10000;
fprintf('Для расчета принимаем Ap=%d cм^2\n\n', Ap*10000);
Ap=pi*2^2/10000;

disp('Коэффициент жесткости "гидравлической пружины" гидроцилиндра');
V0=27.2/1000000 + 1.2*Ap*ymax;


cr=2*E*Ap^2/V0;
%cr=2.74*10^5*10000;
disp('cr=2*E*Ap^2/V0');
fprintf('cr=%f H/м^2\n', cr);
fprintf('cr = %.3f*10^3 H/cм^2\n\n', cr/1000/10000);

disp('Обобщенная жесткость гидромеханической системы привода');
%V0=27.2/1000000 + 1.2*Ap*ymax;
%V0=0.00005882;
disp('c1=(c_k*cr)/(c_k+cr)');
c1=(c_k*cr)/(c_k+cr);
%disp(sprintf('%f', c_k+cr));
fprintf('c1 = %.2f*10^4 H/cм\n\n', c1/10000/100);

disp('Частота собственных колебаний массы нагрузки');
wc=sqrt( (c1+c_sh)/m );
disp('wc = sqrt( (c1+c_sh)/m )');
fprintf('wc = %.0f 1/c\n\n', wc);

disp('Полуось эллипса по оси скоростей');
Vm= A2*(1+c_sh/c1)*(1-w^2/wc^2)*w;
disp('Vm=A2*(1+c_ш/c1)*(1-w^2/wc^2)*w');
fprintf('Vm = %.1f см/с\n\n', Vm*100);

disp('Уравнение эллипса нагрузки');
fprintf('( V/%.1f )^2 + ( (F-10000)/%.1f )^2 = 1 \n\n', Vm*100, Fm);

disp('Скорость гиропривода, при которой сила нагружения F=Fk');
Vk = Vm * sqrt( 1 - ((Fk-F0)/Fm)^2 );
%Vk=14.6 / 100;
disp('Vk = Vm * sqrt( 1 - ((Fk-F0)/Fm)^2 )');
fprintf('Vk = %.1f см/с\n\n', Vk*100);

disp('Гидралическая проводимость управляемого дросселя золотника');
Gm = Vk*Ap / sqrt( 1/2*(Ppit-Fk/Ap) );
disp('Gm = Vk*Aп / sqrt( 1/2*(Pпит-Fk/Aп) )');
fprintf('Gm = %.2f см^4/(H^1/2*с)\n\n', Gm*100000000);

disp('Площадь рабочего окна управляемого дросселя при µ=0.7 и x=xm');
u=0.7;
A0 = Gm/(u*sqrt( 2/ro ));
disp('A0 = Gm/(µ*sqrt( 2/p ))');
fprintf('A0 = %.1f мм^2\n\n', A0*1000000);

disp('Длина прямоугольного рабочего окна при xm=0.05см');
xm=0.05/100;
b = A0/xm;
disp('b = A0/xm');
fprintf('b = %.1f мм\n\n', b*1000);

disp('Скорость гиропривода на холостом ходу (F=0)');
Vxx=Gm/Ap*sqrt( 0.5*Ppit );
disp('Vxx=Gm/Aп*sqrt( 0.5*Pпит )');
fprintf('Vxx = %.1f см/c\n\n', Vxx*100);

disp('Уравнение механической характеристики гидропривода');
disp('Vp=Gm/Aп * sqrt( 0.5*(Pпит-F/Aп) )');
fprintf('Vp=%.2f * sqrt( %0.1f- %0.4f F) ) \n\n', Gm*100000000/Ap/10000*sqrt(11), 0.5*Ppit/10000/11, 1/2/Ap/10000/11);

disp('Пусковое значение движущего усилия гидропривода');
disp('Fпуск = Aп*Pпит');
Fpusk=Ap*Ppit;
fprintf('Fпуск = %0.0f H\n\n', Fpusk);

disp('Пусковое значение движущего усилия гидропривода');
disp('Nmax = 2/3*Pпит*Vk/100');
Nmax = 2/3*Fpusk*Vk;
fprintf('Nmax = %0.0f Вт\n\n', Nmax);