(* Коэффициенты лин-й модели ЛА *)
a1={0.3230,0.4970,0.5110}
a2={0.7050,1.6740,2.3810}
a3={0.7720,2.2650,3.0420}
a4={0.1022,0.1247,0.1302}
a5={0.4254,1.7379,2.5529}
a6={1.5198,0.7965,0.7391}
a7={0.0327,0.0451,0.0470}
b1={1.8500,2.5220,2.7250}
b2={2.2220,5.7440,8.3700}
b3={1.3890,3.6660,4.9240}
b4={0.1730,0.0721, 0.0598}
b6={-0.0780,0.0240,0.0590}
b7={0.0174,0.0436,0.0174}

(* Выбор режима 1..3 *)
rezh=3

(* Передаточные числа АП крена с ЖОС *)
tr=4 (* зададим время регулирования 4с *)
myu=(18-b1[[rezh]]*tr)/(b3[[rezh]]*tr)
ii=108/(b3[[rezh]]*tr*tr)
nyu=216/(b3[[rezh]]*tr*tr*tr)

(* Матрмца состояния свободного самолета *)
MatrixForm[MurKrap0={{(s+b1[[rezh]]),0,b3[[rezh]]},{-1,0,0},{myu,nyu/s,-1}}]
(* Матрица с заданным углом крена *)
MatrixForm[MurKrap={{(s+b1[[rezh]]),0,b3[[rezh]]},{-1,s,0},{myu,(ii+nyu/s),-1}}]

(* Расчет передаточной функции *)
w=Simplify[Det[MurKrap0]/Det[MurKrap]]

Wt=Det[MurKrap0]/Det[MurKrap]

(* Перейдем от s к комплексной форме *)
wk=w/.s->I*omega

f=(Abs[wk])^2

fvhod=1 (*входной спектр*)

fvihod=f*fvhod

(* СПЕКТР *)
Plot[fvihod,{omega,0,5},PlotRange->{0,1}]

p=NIntegrate[fvihod,{omega,0,Infinity}]

omega0=0;
shag=0.1;
While[shag≥10^-5,{While[(NIntegrate[fvihod,{omega,0,omega0}]-0.95*p)<0,omega0=omega0+shag],omega0=omega0-shag,shag=shag/10}]
omega0=omega0-shag*5

p95=NIntegrate[fvihod,{omega,0,omega0}]

(* ПЕРИОД КВАНТОВАНИЯ *)
Tkvant=0.1*Pi/omega0 (* где 1 - поправочный коэффициент*)

ht=InverseLaplaceTransform[Wt/s,s,t]

Plot[{ht,1,0.95,1.05},{t,0,10},PlotRange->{0,1.1}]

hf[t_]=Abs[ht]

(* Время регулирования и перерегулирования *)
k=10000;(*кол-во промежутков, где находим решение*)
dt=10/k;
hust=hf[10000](*установившееся значение*);
hmax=hf[0](*зададим начальное макс. значение*);
(*Численное решение уравнений h(t)=0.95 и h(t)=1.05, объединенных в одно |h(t)-1|-0.05=0 и определение максимума*)
For[i=1,i<k+1,{If[Sign[Abs[hf[dt*(i-1)]-hust]-0.05]≠Sign[Abs[hf[dt*(i)]-hust]-0.05],Treg=dt*(i-0.5)],If[hmax<hf[dt*i],hmax=hf[dt*i]]};i++]
sigmaper=Round[100*(hmax-hust)/hust,0.01]; (*перерегулирование*)
Print["Treg = ",Treg," с"] 
Print["Sigma = ",sigmaper," %"]

(*Кол-во периодов квантования за время регулирования*)
Nkvant=Ceiling[Treg/Tkvant]

Stextrop=3(*Степень полинома экстраполяции*);
Masotklon={}(*Массив отклонений*);
hizm=Table[0,{Stextrop}](*Массив последних измеренных (реальных) значений используемых при экстраполяции*);

For[i=Stextrop-1,i<=Nkvant,i++,
{Ttek=Tkvant*i;
hTtek=hf[Ttek];
hizm[[1]]=hf[Ttek-Tkvant*2];
hizm[[2]]=hf[Ttek-Tkvant*1];
hizm[[3]]=hf[Ttek-Tkvant*0];
(*Экстраполяция Метод Адамса*)
hextrop=hizm[[3]]+Tkvant*(23*hizm[[3]]-16*hizm[[2]]+5*hizm[[1]])/12;
Masotklon=Append[Masotklon,Abs[hf[Ttek+Tkvant]-hextrop]];
}];

(*Математическое ожидание отклонения истинного угла крена от экстраполированного значения*)
Matozh=Mean[Masotklon]

(*Среднее квадратическое отклонение отклонения*)
sigma=StandardDeviation[Masotklon]

(* ПОРОГ СРАБАТЫВАНИЯ *)
b=3 (*1..5 - доверительная граница*);
Porog=Matozh+b*sigma