Модель скорости ветра для анализа работы ветроустановок

z0 = 0.01; # коэффициент шероховатости земной поверхности (z0 = 0.01 - ровная поверхность, 0.1 - небольшие препятствия, 1 - городская местность)
vs = 10; # cреднее значение низкочастотной компоненты, которое сохраняется в течении 10 минут, м/c
h = 30; # высота над землей, м
Lt = 150; # длина турбулентности, м 
TF = Lt/vs; # постоянная времени фильтра, с
Tss = 600; # период дискретности для турбулентной функции (принимается равным приблизительно 10-ти минутам) 
beta = 4.20654; # бета-функция(1/2,1/3)
KF = sqrt((2*pi/beta)*(TF/Tss)); # коэффициент усиления фильтра
It = 1/(log(h/z0)); # интенсивность турбулентности
psi = It*vs; # cреднеквадратичное отклонение скорости ветра 

m1 = 0.4; # коэффициенты формирующего фильтра
m2 = 0.25; # коэффициенты формирующего фильтра

f1 = 0.001;
f2 = 0.003;
f3 = 0.006;
f4 = 0.01;
f5 = 0.03;
f6 = 0.06;
f7 = 0.1;
f8 = 0.3;
f9 = 0.6;
f10 = 1;
f11 = 3;

f=[f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11];

# Cпектральная плотность мощности колебаний 
Sv1 = (4*psi^2)*(Lt/vs)/((1+6*f1*(Lt/vs))^(5/3));                    
Sv2 = (4*psi^2)*(Lt/vs)/((1+6*f2*(Lt/vs))^(5/3));
Sv3 = (4*psi^2)*(Lt/vs)/((1+6*f3*(Lt/vs))^(5/3));
Sv4 = (4*psi^2)*(Lt/vs)/((1+6*f4*(Lt/vs))^(5/3));
Sv5 = (4*psi^2)*(Lt/vs)/((1+6*f5*(Lt/vs))^(5/3));
Sv6 = (4*psi^2)*(Lt/vs)/((1+6*f6*(Lt/vs))^(5/3));
Sv7 = (4*psi^2)*(Lt/vs)/((1+6*f7*(Lt/vs))^(5/3));
Sv8 = (4*psi^2)*(Lt/vs)/((1+6*f8*(Lt/vs))^(5/3));
Sv9 = (4*psi^2)*(Lt/vs)/((1+6*f9*(Lt/vs))^(5/3));
Sv10 = (4*psi^2)*(Lt/vs)/((1+6*f10*(Lt/vs))^(5/3));
Sv11 = (4*psi^2)*(Lt/vs)/((1+6*f11*(Lt/vs))^(5/3));

Svf1 = Sv1*f1;
Svf2 = Sv2*f2;
Svf3 = Sv3*f3;
Svf4 = Sv4*f4;
Svf5 = Sv5*f5;
Svf6 = Sv6*f6;
Svf7 = Sv7*f7;
Svf8 = Sv8*f8;
Svf9 = Sv9*f9;
Svf10 = Sv10*f10;
Svf11 = Sv11*f11;

Sv = [Sv1,Sv2,Sv3,Sv4,Sv5,Sv6,Sv7,Sv8,Sv9,Sv10,Sv11];
Svf = [Svf1,Svf2,Svf3,Svf4,Svf5,Svf6,Svf7,Svf8,Svf9,Svf10,Svf11];
  
w1 = 2*pi*f1;
w2 = 2*pi*f2;
w3 = 2*pi*f3;
w4 = 2*pi*f4;
w5 = 2*pi*f5;
w6 = 2*pi*f6;
w7 = 2*pi*f7;
w8 = 2*pi*f8;
w9 = 2*pi*f9;
w10 = 2*pi*f10;
w11 = 2*pi*f11;

# значение амплитуд косинусоид
A1 = (2/pi)*(sqrt((1/2*(Sv1+Sv2))*(w2-w1)));
A2 = (2/pi)*(sqrt((1/2*(Sv2+Sv3))*(w3-w2)));
A3 = (2/pi)*(sqrt((1/2*(Sv3+Sv4))*(w4-w3)));
A4 = (2/pi)*(sqrt((1/2*(Sv4+Sv5))*(w5-w4)));
A5 = (2/pi)*(sqrt((1/2*(Sv5+Sv6))*(w6-w5)));
A6 = (2/pi)*(sqrt((1/2*(Sv6+Sv7))*(w7-w6)));
A7 = (2/pi)*(sqrt((1/2*(Sv7+Sv8))*(w8-w7)));
A8 = (2/pi)*(sqrt((1/2*(Sv8+Sv9))*(w9-w8)));
A9 = (2/pi)*(sqrt((1/2*(Sv9+Sv10))*(w10-w9)));
A10 = (2/pi)*(sqrt((1/2*(Sv10+Sv11))*(w11-w10)));