close all
clc
% Hodnoty pro výpoèet:

p=3;             %poèet polù
Q1=36;           %drážky statoru
Q2=33;           %drážky rotoru
f1=50.02;           %frekvence nap. napìtí
N1=87;           %poèet závitù fáze statoru
N2=0.5;          %poèet závitù rotoru
m1=3;            %poèet fází statoru
m2=Q2;           %poèet fází rotoru m2=Q2
I1=3.52;          %proud fází statoru (zapojení Y)
y=5;             %krok vinutí
q1=2 ;           %poèet drážek na pol a fázi
tp=6 ;           %polová rozteè
td1=9.468e-3;    %vzdálenost os sousedních zubù statoru 9.456e-3
td2=10.329e-3;   %vzdálenost os sousedních zubù rotoru 10.258e-3
airgap=0.25e-3;  %vzduchová mezera
o1=3.05e-3;      %otevøení drážky statoru
o2=0;            %otevøení drážky rotoru
beta1=0.41;      %souèinitel carterova koeficientu
beta2=0;         %souèinitel carterova koeficientu
U1n=230;         %napájecí fázové napìtí
kv2=1;           %èinitel vinutí rotoru
n=0;             %otáèky
alfa=0:0.01:2*pi;
alfar=0:0.01:2*pi;
cas=0.005;
p0=2e-5;         %vztažný akustický tlak
D=108.5e-3;      %vnitøní prùmìr statoru
lfe=90e-3;       %délka statorového svazku
ki=0.95;         %koeficient proudu rotoru -závislé na úèiníku(viz. BP obrázek 4-2)

%výpoèet skluzu
ns=f1*60/p
s=(ns-n)/ns

w=2*pi*f1
%výpoèet èinitele vinutí statoru
a=2*p*pi/Q1
kr=(sin(q1*a/2))/(q1*sin(a/2))
ky=sin(pi*y/(2*tp))
kv1=kr*ky

%výpoèet proudu v tyèi klece rotoru
pi2=6*N1*kv1/Q2
I2=ki*pi2*I1

%MMN statorového vinutí
S=0
for v=p;
    S=S+(1/v)*kv1*sin(w*cas+(v*alfa));
end
S1=0
for c=1:50
    v=p*(2*c*m1+1);
    S1=S1+(1/v)*kv1*sin(w*cas+(v*alfa));
end
S2=0
for c=1:50
    v=p*(2*c*m1+1);
    S2=S2+(1/v)*kv1*sin(w*cas-(v*alfa));
end
S3=0
for c=1:50
    v2=p*(2*c*m1-1);
    S3=S3+(1/v2)*kv1*sin(w*cas+(v2*alfa));
end
S4=0
for c=1:50
    v2=p*(2*c*m1-1);
    S4=S4+(1/v2)*kv1*sin(w*cas-(v2*alfa));
end

Fs=(3*sqrt(2)*N1*I1/pi)*(S+S1+S2+S3+S4)
  

%MMN rotorového vinutí
wr=2*pi*f1*(1-s);
r=0
for vr=p;
    r=r+(1/vr)*sin(wr*cas+vr*alfa);
end
r1=0
for cr=1:50
    vr=cr*Q2+p;
    r1=r1+(1/vr)*sin(wr*cas+vr*alfa);
end
r2=0
for cr=1:50
    vr=cr*Q2+p;
    r2=r2+(1/vr)*sin(wr*cas-vr*alfa);
end
r3=0
for cr=1:50
    v2r=cr*Q2-p;
    r3=r3+(1/v2r)*sin(wr*cas+v2r*alfa);
end
r4=0
for cr=1:50
    v2r=cr*Q2-p;
    r4=r4+(1/v2r)*sin(wr*cas-v2r*alfa);
end

Fr=(Q2*I2*sqrt(2)/(2*p*pi))*(r+r1+r2+r3+r4);


%magnetická vodivost vzduchové mezery-
    %carterùv koeficient - stator
    gama1=(4/pi)*((o1/(2*airgap))*atan(o1/(2*airgap))-log(sqrt(1+(o1/(2*airgap))^2)))
    kc1=td1/(td1-gama1*airgap)
    
    %carterùv koeficient - rotor
    gama2=(4/pi)*((o2/(2*airgap))*atan(o2/(2*airgap))-log(sqrt(1+(o2/(2*airgap))^2)))
    kc2=td2/(td2-gama2*airgap)
    %carteruv koeficient
    kc12=kc1*kc2
    
    
   %mag. vodivost statrou
    a01=1/(kc1*airgap);

SUMA1=0
for vv=1:100
    SUMA1=SUMA1+((beta1*4/(airgap*vv*pi))*(0.5+((vv*o1/td1)^2/(0.78-2*(vv*o1/td1)^2)))*sin(1.6*pi*vv*o1/td1))*cos(vv*Q1*alfar);
end

A1=a01-SUMA1;
    
    %mag. vodivost rotoru
b01=1/(kc2*airgap);

SUMA2=0
for vv=1:100
    SUMA2=SUMA2+((beta2*4/(airgap*vv*pi))*(0.5+((vv*o2/td2)^2/(0.78-2*(vv*o2/td2)^2)))*sin(1.6*pi*vv*o2/td2))*cos(vv*Q2*alfar);
end

A2=b01-SUMA2;

Ss=pi*D*lfe;       
delta=(1./A1)+(1./A2)-airgap;
A12=Ss*4e-7*pi./(delta);


%výpoèet radiálního tlaku pr
pr=(((Fs+Fr).^2).*(A12).^2)/(2*4*pi*1e-7);
pr1=pr.*(54.25/1113.5);
%výpoèet akustické hladiny radiálního tlaku
Lpr=20*log10(pr/p0);
Lpr1=20*log10(pr1/p0);


figure(1)
subplot(5,1,1); plot(alfa,Fs)
xlabel('alfa[rad]')
ylabel('Fs[A]')
grid on
subplot(5,1,2); plot(alfa,Fr) 
xlabel('alfa[rad]')
ylabel('Fr[A]')
grid on
subplot(5,1,3); plot(alfa,A12) 
xlabel('alfa[rad]')
ylabel('A12[H]')
grid on
subplot(5,1,4); plot(alfa,pr) 
xlabel('alfa[rad]')
ylabel('pr[Pa]')
grid on
subplot(5,1,5); plot(alfa,Lpr)
xlabel('alfa[rad]')
ylabel('Lpr[dB]')
grid on

figure(2)
plot(alfa,Fs)
xlabel('alfa[rad]')
ylabel('Fs[A]')
grid on

figure(3)
plot(alfa,Fr)
xlabel('alfa[rad]')
ylabel('Fr[A]')
grid on

figure(4)
plot(alfa,A12)
xlabel('alfa[rad]')
ylabel('A12[H]')
grid on

figure(5)
plot(alfa,pr)
xlabel('alfa[rad]')
ylabel('pr[Pa]')
grid on

figure(6)
plot(alfa,Lpr)
xlabel('alfa[rad]')
ylabel('Lpr[dB]')
grid on

%výpoèet frekvence mag.hluku(max.složky)
nn=940;              %jmenovité otáèky
sn=(ns-nn)/ns        %skluz pøi jmenovitých otáèkách
fc=f1*(Q1/(2*p)+Q2/(2*p))*(1-sn)



