Егоркины графики по Вахитову

clear all;
clf;
close all;

X_STEPS = 1;
global plot_num = 0;
global plots = [2 4];



%pt3;
D = 275E-3;
a_d = D/2;
num_req = 3;
K_em = 17;
K_fem = 3;
L_k = 7E-3;
R_k = 16;
F_m = 50;

S_d = pi * a_d^2;
W_el = 75;
rho_v = 1.23; % кг/м^3

m_1 = 46E-3;
r_1 = 22.5;
c_1 = 2.2E-4;
c_sum = 7E-5;
V_v = 340; % м/с



l_z = 31.5;




function myfig (X1, fun, arg, X, XLabels, MT)
  global plot_num;
  global plots;
  plot_num++;
  figure (gcf, 'papertype', 'A4', 'paperunits', 'centimeters', 'papersize', [21 29.7], 'Position', [10 10 810 610]);
  %orient('Landscape');
  orient('Portrait');
  subplot(plots(1), plots(2), plot_num);
  A1 = semilogx(X1, fun);
%  MT
  title(MT);
  grid('on');
  set(gca, 'XMinorGrid', 'off');
  set(gca, 'XMinorTick', 'off');
  set(gca, 'XLabel', 'f, Hz/kHz');
  set(gca, 'YLabel', arg);
  set(gca, 'XLim', 1e3 * [X(1) X(size(X,2)) ]);
  set(gca, 'XTick', X * 1E3);
  set(gca, 'XTickLabel', XLabels);

%  legend('location', 'northwest');
%  set(gcf, 'Position', [10 10 810 610]);
endfunction



F_v = 2*V_v/(a_d*pi);


points    = 1E4; % Количество точек
%oct_start = -8; % Начало графика (октава)
oct_start = log2(F_m * 1e-3);
%oct_end   = 4; % Конец  графика (октава)
oct_end   = log2(F_v * 1e-3);

X       = 2.^[floor(oct_start):X_STEPS:ceil(oct_end)]; % отсчеты по X
XLabels = round(X*1E3);
%XLabels = [round(2.^[oct_start:X_STEPS:oct_end]*1E3)]% .1*round(10*2.^[0:X_STEPS:oct_end])]
X1      = 2.^linspace(oct_start,oct_end,points) * 1e3;

Z_m1    = r_1;
Z_m2    = i * 2*pi * m_1 * X1;
Z_m3    = i ./ (2*pi .* [c_1; c_sum] .* X1);
Z_meh = Z_m1 + Z_m2 .- Z_m3;

Z_1 = R_k;
Z_2 = i * 2*pi * L_k * X1;
Z_3 = K_em^2 ./ (Z_meh);
Z   = Z_1 + Z_2 + Z_3;
clear Z_?

M_chg = rho_v / (sqrt(2)) * (i*K_em*S_d*X1) ./ (Z_meh(2,:) .* Z(2,:));
K_omb = sqrt (num_req * (1 + (pi*X1*a_d/V_v).^3));
%K_omb_fm = sqrt (num_req * (1 + (pi*F_m*a_d/V_v).^3));
K_omb_fm = sqrt (num_req * (1 + (pi*F_m*K_fem*a_d/V_v).^3));
M_chk = M_chg .* K_omb;

%E_hchk = E_hch * K_omb_fm;
E_hchk = rho_v * K_em * S_d / (sqrt(8*R_k)*pi * m_1) * K_omb_fm;
N_hchk = 20*log10(E_hchk/2E-5);
p_lN = E_hchk * sqrt(num_req * W_el) / l_z;

N_lN = 20*log10(p_lN/2e-5);



THETA = linspace(-pi/2, pi/2, points);
ka = [.25 .5 1 2 4]';
R_hor = abs(2*besselj(1,ka*sin(THETA))./(ka*sin(THETA)));
R_vert = R_hor .* abs((sin(3*ka*num_req * sin(THETA)/2)) ./ (num_req * sin(3*ka * sin(THETA)/2)));

plot_num++;
subplot(plots(1), plots(2), plot_num)
polar(THETA, R_hor);
title('Горизонтальная плоскость колонки');
%legend (num2str(ka));

plot_num++;
subplot(plots(1), plots(2), plot_num)
polar(THETA, R_vert);
title('Вертикальная плоскость колонки');
%legend (num2str(ka));


myfig (X1, abs(Z(1,:)), '|Z|', X, XLabels, 'Зависимость модуля Z от частоты');
myfig (X1, 180*angle(Z(1,:)) / pi, '\phi(Z)', X, XLabels, 'Зависимость фазы Z от частоты');

myfig (X1, abs(M_chg), 'Па/В', X, XLabels, 'ЧХЗД динамической головки');
myfig (X1, 20*log10(abs(M_chg)/max(abs(M_chg))), 'дБ', X, XLabels, 'ЧХЗД динамической головки (дБ)');

myfig (X1, abs(M_chk), 'Па/В', X, XLabels, 'ЧХЗД колонки');
myfig (X1, 20*log10(abs(M_chk)/max(abs(M_chk))), 'дБ', X, XLabels, 'ЧХЗД колонки (дБ)');