IC_4DOF_new2.m

%% Задание безразмерных параметров

t0 = 0; tend =0.01; %с, время для решателя
%---------------------------------------

%%Обезразмерить x0
%%

%-----Размерные параметры------------------
nu=1e4; % rpm 1e4 Частота вращения двигателя В Об.сек = w=omega
p=0;  % Гц, частота внешнего воздействия, варьируем
% w=1017; % s^-1, Частота вращения ротора в 1017 рад/c примерно 10000 об/мин
w=(pi*nu)/30;
omega=w; %переобозначили
omega_square=w^2;
T0 = 2*pi / 1; % период собственных колебаний

a1=9*1e-3; %м, точка установки АМП1
a2=9*1e-3; %м, точка установки АМП2
lm=0.004; %м, длина э.двигателя
l=22*1e-3; % m,длина ротора
r=7.7*1e-3; % m,радиус ротора
h=0.2*1e-3; % m, зазор между ротором и АМП
% tolsh=0.003; %мм, толщина стенки ротора
e=4.27*1e-7; %4.27*1e-7; % m, эксцентриситет ???
delta_dim=h; %принимаем безразмерную величину зазора
h=0.2*1e-3; % m, зазор между ротором и АМП
m=(12.42*1e-3)/4;%кг,масса Э.Д. ????
cm=6000;%Н*м^-1 жесткость электродвигателя
m2=12.42*1e-3; %кг, масса ротора

gamma=0; %2.15*1e-5; %m2*e^2; %Дебаланс, тут m- масса ротора! 2.15*1e-5 рад
cw=0.91; %коэффициент сопротивления среды крови в аорте
etta=0.003;% Па*с ,динамическая вязкость крови Па*с
P_e=100;% Вт, Мощность электрическая двигателя???
R_e=40; %Ом, сопротивление обмоток ???????
Iz=3.78*1e-7; % кг*м^2, Момент инерции сечения  ротора относительно центра масс
Ix=6.9*1e-7; %kg*m^2 Момент инерции сечения  ротора относительно оси X
ksi2=0.3926991; %рад ,угол отклонения магнитной силы от вертикали/горизонтали равный 22,5°
mu0=4*pi*1e-7; %м*кг*с^(-2)*А^(-2), магнитная постоянная 
n=100; %число витков обмотки ????????
A=1e-6; %м^2, площадь сечения , занимаемя обмоткой ?????
ka=(mu0*n^2*A)*cos(ksi2)/2;%конструктивный параметр Э.Д.  0.202
kb=ka;%по условию

%Управляющие токи АМП в соответствующих направлениях
i0=0.01;%0.000001;  %A,i_0=i_max/2 ???
ia.x=0; %А, токи смещения в виде структуры  ?????
ia.y=0; %А, токи смещения ????
ib.x=0; %А, токи смещения ????
ib.y=0; %А, токи смещения????

%Внешние возмущения и моменты
Mx=0.03; %0.03;% Н*м, ????????
My=Mx;
Px=Mx/2/r; %Н,
Py=Px;

%Амплитуды внешних нагрузок
Ax=1e-3; %1e-3; %m
Ay=Ax; %m
Bx=2*Ax; %m
By=2*Ax; %m

%-----------Масштабы------------------
Si=0.2;%Масштаб токов ???????? корень из P_e/R_e
Sl=l; %  масштаб длины
Sf=h*cm; %  Масштаб силы  
Sm=2*pi*nu*r^3*l*etta/(30*h); %Масштаб моментов Ix*(h/l)*cm/m

% St=60/nu;
% St=2*pi/omega; %Масштаб времени 60/nu
St=sqrt(m/cm); %масштаб из диссера
Sw=1/St; %Масштаб угловой скорости
Sx=h; %масштаб координаты закрепления АМП 
Sbeta=h/l; %Масштаб угловой координаты 

%%---------- Начальные условия----------
x1=h/2;
y1=h/4;
alpha1=0.02; %0.02 -max угол наклона ротора
beta1=0.015; %0.015
vx1=0.1; %0.1
vy1=0.2; %0.2
om1=0.1;
om2=0.2;

% Sx=h; %масштаб координаты закрепления АМП 
% Sbeta=1*e-2; %Масштаб угловой координаты ????????????
% 
x1=x1/Sx;
y1=y1/Sx;
alpha1=alpha1/Sbeta;
beta1=beta1/Sbeta;
om1=om1/Sw;
om2=om2/Sw;

x0 = [x1, y1, alpha1, beta1,...
    vx1, vy1 , om1, om2]; % [4 перем., 4 скор.] 8 параметров

%------------Переход к безразмерным величинам--------
r_dim=r/Sl; 
h_dim=h/Sl;
l_dim=l/Sl;
lm_dim=lm/Sl;

Px_dim=Px/Sf;
Py_dim=Py/Sf;
Mx_dim=Mx/Sm;
My_dim=My/Sm;
Ax_dim=Ax/Sf;
Ay_dim=Ay/Sf;
Bx_dim=Bx/Sm;
By_dim=By/Sm;

delta_dim=delta_dim/Sx;
e_dim=e/Sx;

% Fax_dim=Fax/Sf;
% Fay_dim=Fay/Sf;
% Fbx_dim=Fbx/Sf;
% Fby_dim=Fby/Sf;

%t_dim=t/St;
a1_dim=a1/Sl;
a2_dim=a2/Sl;

omega_dim=omega/Sw;
p_dim=p/Sw;
gamma_dim=gamma/Sbeta;

%Безразмерные токи
i0_dim=i0/Si;
iax_dim=ia.x/Si;
iay_dim=ia.y/Si;
ibx_dim=ib.x/Si;
iby_dim=ib.y/Si;

%-----------Критерии подобия------------------
Pn=(2*pi*cw*etta*r_dim^3*l^2*l_dim)/(h_dim*(m*cm)^(1/2)); %гидродинамика 1.48*1e-7
P1=Iz/Ix; %0.548
P2=(m*l^2)/Ix; %8.714
Pk=ka*P_e/(cm*R_e*delta_dim^3); %конструктивные параметры АМП 10^5

%%----------Постоянные матрицы---------------

%Коэффициенты
km=((lm_dim)^2)/(3*(a1_dim+a2_dim)^2); %Коэффициент безразмерный
Km=((a1_dim)^2+(a2_dim)^2)*(1+km)-2*a1_dim*a2_dim*(1-km);
%Матрица Жесткости
C=[ 1 0 0 (-a1_dim+a2_dim)/2;
    0 1 (a1-a2)/2 0;
    0 (P2*(a1_dim-a2_dim))/2 (P2*Km)/4 0;
    (P2*(-a1_dim+a2_dim))/2 0 0 (P2*Km)/4];
%симметричная часть
Hs=[0 0 0 Pn;
    0 0 -Pn 0;
    0 0 0 0;
    0 0 0  0];
%Гироскопическая матрица
G=[0 0 0 0;
    0 0 0 0;
    0 0 0 -P1;
    0 0  P1 0];
%
Hd=[0 0 -Pn 0;
    0 0 0 -Pn;
    0 0 0 0;
    0 0 0 0];

%Матрица перехода
Zb=[1 0 0 -a1_dim;
    0 1 a1_dim 0;
    1 0 0 a2_dim;
    0 1 -a2_dim 0];
%---------------------------------


%-----------TEST---------------------------
% a = .3*1 : 0.01 : 1.5; % к-т пропорциональности амплитуды силы тока
a = 0.3;
%nu = 0.1: 0.1: 1.4; % частота возбуждения
kp = [0.1 : 0.1 : 1.2]; % к-т относительного поджатия = kp_dim/(4*i0^2*k)
% Параметры внешнего возбуждения (ВВ), с которым согласована частота nu
T0 = 2*pi / 1; % период собственных колебаний
aext = 0.2; % доля периода, в которую подаётся ВВ
for nuNum = 1 : length(nu)
    
    Text = 2*pi / nu(nuNum); % период ВВ
    pulseT{ nuNum } = [0 : Text : tend]; % моменты времени, в которые происходит внешнее возбуждение
    width(nuNum) = aext* Text ; % ширина внешнего импульса в долях T0 или Text ??? -> стационар vs нестационар
    
    % Нулевое внешнее воздействие - закомментировать в случае наличия
    % внешней силы
    pulseT{ nuNum } = 0;
    width(nuNum) = 0;
end