Stochastic model of equivalent wind

Initialization

Needs["ColorHueGray2`",NotebookDirectory[]<>"ColorHueGray2.m"];
	
Get[NotebookDirectory[]<>"ColorHueGray2.m"];

<< (NotebookDirectory[] <> "ColorHueGray2.m")

j = I;

Funcion and tab definitions

Ganancia Bandpass

gananciaBandPass[σ_, Bandwidth_] = Integrate[BandPassPSD[f, σ, fblade0, Bandwidth], {f, 0, ∞}, Assumptions -> {σ > 0, fblade0 > 0, 1 > Bandwidth > 0}]

(((Sqrt[2]) Bandwidth (σ^2))/((Sqrt[-2+(Bandwidth^2)-i Bandwidth (Sqrt[4-(Bandwidth^2)])])+(Sqrt[-2+(Bandwidth^2)+i Bandwidth (Sqrt[4-(Bandwidth^2)])])))

gananciaBandPass[σ_, Bandwidth_] = Assuming[{σ > 0, fblade0 > 0, 0.5 > Bandwidth > 0}, FullSimplify@ComplexExpand[gananciaBandPass[σ, Bandwidth]]]

((Bandwidth (σ^2))/((e^((1/2) i Arg[-2+(Bandwidth^2)-i Bandwidth (Sqrt[4-(Bandwidth^2)])]))+(e^((1/2) i Arg[-2+(Bandwidth^2)+i Bandwidth (Sqrt[4-(Bandwidth^2)])]))))

gananciaBandPass[σ,0.25]

TabWindModel

IniWindModel[]:= Module[{},
(*S2harris[f_,U_]=Simplify[ ((4 (uf^2) L U)/(((2+(x^2)))^(5/6)))//.{ uf->((k U)/(Log[z/z0])) , k-> 0.4, z0-> 0.0025 (*surface drag coefficient*),  x-> L f/z ,L-> 1800,  z -> 10} , {f>0, U>0}];*)
PSDh1[f_, a_, b_]:=((a+2 a (b^2) (f^2) (π^2))/((a+4 (b^2) (f^2) (π^2)) (1+2 a (b^2) (f^2) (π^2))));
S2karman[f_,L2_,Uwind_,σ_]:=(((σ^2)4 L2/Uwind )/(((1+70.8(((f L2/Uwind))^2)))^(5/6)));S2kaimal[f_,L1_,Uwind_,σ_]:=(((σ^2)4 L1/Uwind )/(((1+6 f L1/Uwind))^(5/3))) ;

Davenport[f_,Λ1_,σ_]:=(σ^2)((975.8148560061259` f ((Λ1)^2))/(((49+400 (f^2) ((Λ1)^2)))^(4/3))) ;
 IEC614001[f_,Λ1_,Uwind_,σ1_]:=0.05 ((σ1)^2) (((Λ1/Uwind))^(-2/3))(f^(-5/3));
(*Blade gain (multiplicative model respect aerodynamic filter)*)PSDh2[f_,gd_,fblade_]:=((((((f^2)+(fblade^2)))^2) (((1+(gd^2)))^2))/((fblade^4)+2 (fblade^2) ((f^2)+(fblade^2)) (gd^2)+((((f^2)-(fblade^2)))^2) (gd^4)));
(*Blade gain (multiplicative model respect aerodynamic filter)*)
BandPassPSD[f_,σ_,fblade_, Bandwidth_]=((2 Bandwidth (f^2) (σ^2))/(((Bandwidth^2) (f^2)+((((f^2)-(fblade^2)))^2)) π));

(*P.Teigen,A.Naess,"Extreme Response of Floating Structures in Combined Wind and Waves".Slettringen spectrum recommended by the Norwegian Petroleum Directorate (NPD).*)
(*IntensidadTurbulenciaPorDefecto para n->0.468 => 0.019798255944295007` (U10^(3/8)) (z^(1/6))*)
(* Para las turbinas estudiadas, la intensidad es del orden del 12% *)
PSDconsorcionoruegopetroleo[f_,U10_,n_,z_]:=((320 (z/10)0.45*(U10/10)^2)/(((1+(((172 f(((z/10))^(2/3))(((U10/10))^(-3/4))))^n)))^(5/(3n))));
(* Versión escalada para fijar una desviación estándar determinada *)
PSDconsorcionoruegopetroleo[f_,U10_,n_,z_,σ_]:=(σ^2)((PSDconsorcionoruegopetroleo[f,U10,n,z])/((0.0006910360620848527` (U10^(11/4)) (z^(1/3)) Gamma[1+(1/n)] Gamma[(2/(3 n))])/(Gamma[(5/(3 n))])));
PowerCurve[v_]= (1/2) (Tanh[27.465307216702744`-1.0986122886681098` v]-Tanh[3.295836866004329`-0.3662040962227032` v]);
Λ1=42;σkaimal[v01_,Iref_]:= Iref(0.75 v01 + 5.6);Iref = 0.14;v0 = 10;fblade=1;R=41;TipoEspectro=1;GananciaPasoPala=0.5; fblade = 1; rotorfilter =2;spatialsampling = False;BandWidth=0.05;TowerShadowIntensity=0.05;
NivelGris={0,(1/2),(1/3),(5/6),(2/3),(1/6),1};
tamano = {800,400};
PuntosYrayas={{},{},{0.02},{0.04,0.03},
{0.08,0.02,0.01,0.02},{0.08,0.02},{0.01,0.02},{0.005,0.02,0.05,0.02},{0.005,0.02,0.02,0.02},{0.01},{0.06,0.02,0.01,0.02,0.01,0.02}};
EstiloGrafica ={
{GridLines-> Automatic,Frame->True, Axes-> False,AxesOrigin->{1,0},PlotStyle->Table[{Hue[(-y)/6], Thickness[0.01],Dashing[PuntosYrayas[[y+1]]]},{y,0,6}], Filling-> Axis,Background->None, ImageSize->tamano},
{GridLines-> Automatic,Frame->True, Axes-> False,AxesOrigin->{1,0},PlotStyle->Table[{Hue[(-y)/6], Thickness[0.006*If[y>0,1.5,2.5]],Dashing[PuntosYrayas[[y+1]]]},{y,0,6}], Filling-> Axis,Background->None, ImageSize->tamano},
{GridLines-> Automatic,Frame->True, Axes-> False,AxesOrigin->{1,0},PlotStyle->Table[{Hue[(-y)/6], Thickness[0.006*If[y>0,1.5,2.5]],Dashing[{}]},{y,0,6}], Filling-> Axis,Background->None, ImageSize->tamano},
{GridLines-> None,Frame->True, Axes-> False,AxesOrigin->{1,0},
PlotStyle->Table[{Hue[-y/10], Thickness[0.01],Dashing[PuntosYrayas[[y+1]]]},{y,0,10}], Filling-> Axis,Background->None},
{GridLines-> Automatic,Frame->True, Axes-> False,AxesOrigin->{1,0},PlotStyle->Table[{GrayLevel[NivelGris[[y+1]]], Thickness[0.01],Dashing[PuntosYrayas[[y+1]]]},{y,0,6}], Filling-> Axis,Background->None, ImageSize->tamano},
{GridLines-> Automatic,Frame->True, Axes-> False,AxesOrigin->{1,0},PlotStyle->Table[{GrayLevel[NivelGris[[y+1]]], Thickness[0.006*If[y>0,1.5,2.5]],Dashing[PuntosYrayas[[y+1]]]},{y,0,6}], Filling-> Axis,Background->None, ImageSize->tamano}
};
EstiloLogPlot ={Frame->True, Axes-> False,AxesOrigin->{1,0},PlotStyle->Table[{GrayLevel[NivelGris[[y+1]]], Thickness[0.006*If[y>0,1.5,2.5]],Dashing[PuntosYrayas[[y+1]]]},{y,0,6}], Filling-> None,Background->None,
FrameLabel-> {"PSD of equivalent wind (m/s)^2/Hz","frequency f [Hz]"} ,GridLines->{Automatic,Automatic},GridLinesStyle->Directive[Yellow], ImageSize->tamano};
(SetOptions[Plot, #1] &) /@ EstiloGrafica[[1]];
(SetOptions[LogLogPlot, #1] &) /@  EstiloLogPlot;
(SetOptions[ListLogLogPlot, #1] &) /@ {Frame-> True, Axes-> False,AxesOrigin->{1,0},Joined->joined, FrameLabel-> {"PSD of equivalent wind (m/s)^2/Hz","frequency f [Hz]"} ,GridLines->{Automatic,Automatic},GridLinesStyle->Directive[Yellow],PlotStyle->{Thick, Thin},Filling-> None, ImageSize->tamano};
(SetOptions[ListLogLinearPlot, #1] &) /@ {Frame-> True, Axes-> False,AxesOrigin->{1,0},Joined->joined,FrameLabel-> {"f·PSD[f] = relative variance (m/s)^2 at freq. f","frequency f [Hz]"}  ,GridLines->{Automatic,Automatic},GridLinesStyle->Directive[Yellow],PlotStyle->{Thick, Thin},Filling-> None, ImageSize->tamano};
(SetOptions[LogLinearPlot, #1] &) /@ {Frame-> True, Axes-> False,AxesOrigin->{1,0},
FrameLabel-> {"f·PSD[f] = relative variance (m/s)^2 at freq. f","frequency f [Hz]"} ,GridLines->{Automatic,Automatic},GridLinesStyle->Directive[Yellow],PlotStyle->{Thick, Thin}, Filling-> Axis, ImageSize->tamano 
};
(SetOptions[ListLogLogPlot, #1] &) /@ {Frame-> True, Axes-> False,AxesOrigin->{1,0},Joined->joined ,GridLines->{Automatic,Automatic},GridLinesStyle->Directive[Yellow],PlotStyle->{Thick, Thin},  Filling-> Axis, ImageSize->tamano,FrameLabel-> {"PSD of equivalent wind (m/s)^2/Hz","frequency f [Hz]"}
};
];

ComputeWindModel[]:= Module[{},
Gain2rotorfilter[f_] = Switch [rotorfilter,
0, (* None, unity gain *) 1,
1, (* First order *) PSDh1[f,0.55,1.3 R/v0],
2, (* Second order *) PSDh1[f,0.55,1.3 R/v0]* PSDh1[f,0.55,1.3 R/v0],
_, (* Default *) 1
];

(* The gain of spatial sampling must be divided by rotor filter since rotor filtering effect does not affect tower shadow and wind shear. In fact, this effect is more important with increasing rotor diameter (in contrast with the turbulence filtering) *)
(*Gain2spatialsampling[f_] =If[spatialsampling,PSDh2[f,Max[((GananciaPasoPala )/(π Gain2rotorfilter[fblade])),1],fblade],1];*)
Gain2spatialsampling[f_] =If[spatialsampling,BandPassPSD[f,(*(Subscript[σ, fblade])*) v0 TowerShadowIntensity,fblade, BandWidth  ],0];
espectroviento[f_]=Gain2spatialsampling[f] +Gain2rotorfilter[f]*
(Switch[TipoEspectro,
1,S2karman[f,0.7130740328122924` 8.1  Λ1 (* Lkarman *),v0,1.184220811841064` σkaimal[v0,TurbInt∞]],
2,S2kaimal[f,8.1  Λ1 (*Lkaimal*),v0,0.7495558272626599` σkaimal[v0,TurbInt∞]],
3,Davenport[f,v0],
4,IEC614001[f,Λ1,v0,(*σ1*)  σkaimal[v0,TurbInt∞=0.16 (*Iref*)]  ],
5,IEC614001[f,Λ1,v0,(*σ1*) σkaimal[v0,TurbInt∞=0.14 (*Iref*)]  ],
6,IEC614001[f,Λ1,v0,(*σ1*)  σkaimal[v0,TurbInt∞=0.12 (*Iref*)]  ],
7,PSDconsorcionoruegopetroleo[f,(*U10 =*)v0,(* n = *) 0.468,(*z = "Altura del buje" *) Hhub,(*σ*)σkaimal[v0,TurbInt∞=0.12 (*Iref*)]],
17,{ S2karman[f,0.7130740328122924` 8.1  Λ1 (* Lkarman *),v0,1.184220811841064` σkaimal[v0,TurbInt∞]],
S2kaimal[f,8.1  Λ1 (*Lkaimal*),v0,0.7495558272626599` σkaimal[v0,TurbInt∞]],
(*Davenport[f,v0],*) IEC614001[f,Λ1,v0,(*σ1*) σkaimal[v0,0.16]  ],IEC614001[f,Λ1,v0,(*σ1*) σkaimal[v0,0.14]  ], IEC614001[f,Λ1,v0,(*σ1*) σkaimal[v0,0.12]  ]},
18,1,
 _, 1 (* default value *)]
+ If[HornsRev,(* Single-sided spectrum -> 2 factor *)2(((((0.0046 v0 + 0))^2)Hhub/v0)/((((Hhub f/v0))^(5/3))(1+100.(Hhub f/v0)))),0]);
]

PlotWindModel[] :=Column[{
Switch [rotorfilter,
0, (* None, unity gain *) "",
1, (* First order *) "Rotor filter => one first order low pass filter with (Subscript[f, cut]) = 0.1224 U∞/ R = " <> ToString@N[0.1224 v0 / R]<>" Hz.",
2, (* Second order *) "Rotor filter => Two cascaded 1st order with (Subscript[f, cut]) = 0.1224 U∞/ R = " <> ToString@N[0.1224 v0 / R]<>" Hz.",
_, (* Default *) ""
],
LogLogPlot[
Gain2spatialsampling[f] +Gain2rotorfilter[f]*
(
Switch[TipoEspectro,
1,S2karman[f,0.7130740328122924` 8.1  Λ1 (* Lkarman *),v0,1.184220811841064` σkaimal[v0,TurbInt∞]],
2,S2kaimal[f,8.1  Λ1 (*Lkaimal*),v0,0.7495558272626599` σkaimal[v0,TurbInt∞]],
3,Davenport[f,v0],
4,IEC614001[f,Λ1,v0,(*σ1*)  σkaimal[v0,TurbInt∞=0.16 (*Iref*)]  ],
5,IEC614001[f,Λ1,v0,(*σ1*) σkaimal[v0,TurbInt∞=0.14 (*Iref*)]  ],
6,IEC614001[f,Λ1,v0,(*σ1*)  σkaimal[v0,TurbInt∞=0.12 (*Iref*)]  ],
7,PSDconsorcionoruegopetroleo[f,(*U10 =*)v0,(* n = *) 0.468,(*z = "Altura del buje" *) Hhub,(*σ*)σkaimal[v0,TurbInt∞=0.12 (*Iref*)]],
17,{ S2karman[f,0.7130740328122924` 8.1  Λ1 (* Lkarman *),v0,1.184220811841064` σkaimal[v0,TurbInt∞]],
S2kaimal[f,8.1  Λ1 (*Lkaimal*),v0,0.7495558272626599` σkaimal[v0,TurbInt∞]],
(*Davenport[f,v0],*) IEC614001[f,Λ1,v0,(*σ1*) σkaimal[v0,0.16]  ],IEC614001[f,Λ1,v0,(*σ1*) σkaimal[v0,0.14]  ], IEC614001[f,Λ1,v0,(*σ1*) σkaimal[v0,0.12]  ]},
18,1,
 _, 1 (* default value *)]
+ If[HornsRev,(* Single-sided spectrum -> 2 factor *)2(((((0.0046 v0 + 0))^2)Hhub/v0)/((((Hhub f/v0))^(5/3))(1+100.(Hhub f/v0)))),0])
,
{f,0.0001,10},Frame-> True, FrameLabel-> {"PSD of equivalent wind (m/s)^2/Hz","frequency f [Hz]"} ,GridLines->{Automatic,Automatic},GridLinesStyle->Directive[Yellow],PlotStyle->Thickness[0.01], PlotRange -> {(10^(-8)),1000},PlotLabel->"Spectral Density of equivalent wind (PSD)", ImageSize->tamano
](* End of LogLogPlot *),
Quiet@LogLinearPlot[f (Gain2spatialsampling[f] +Gain2rotorfilter[f]*(
Switch[TipoEspectro,
1,S2karman[f,0.7130740328122924` 8.1  Λ1 (* Lkarman *),v0,1.184220811841064` σkaimal[v0,TurbInt∞]],
2,S2kaimal[f,8.1  Λ1 (*Lkaimal*),v0,0.7495558272626599` σkaimal[v0,TurbInt∞]],
3,Davenport[f,v0],
4,IEC614001[f,Λ1,v0,(*σ1*)  σkaimal[v0,TurbInt∞=0.16 (*Iref*)]  ],
5,IEC614001[f,Λ1,v0,(*σ1*) σkaimal[v0,TurbInt∞=0.14 (*Iref*)]  ],
6,IEC614001[f,Λ1,v0,(*σ1*)  σkaimal[v0,TurbInt∞=0.12 (*Iref*)]  ],
7,PSDconsorcionoruegopetroleo[f,(*U10 =*)v0,(* n = *) 0.468,(*z = "Altura del buje" *) Hhub,(*σ*)σkaimal[v0,TurbInt∞=0.12 (*Iref*)]],
17,{ S2karman[f,0.7130740328122924` 8.1  Λ1 (* Lkarman *),v0,1.184220811841064` σkaimal[v0,TurbInt∞]],
S2kaimal[f,8.1  Λ1 (*Lkaimal*),v0,0.7495558272626599` σkaimal[v0,TurbInt∞]],
(*Davenport[f,v0],*) IEC614001[f,Λ1,v0,(*σ1*) σkaimal[v0,0.16]  ],IEC614001[f,Λ1,v0,(*σ1*) σkaimal[v0,0.14]  ], IEC614001[f,Λ1,v0,(*σ1*) σkaimal[v0,0.12]  ]},
18,1,
 _, 1 (* default value *)]
+ If[HornsRev,(* Single-sided spectrum -> 2 factor *)2(((((0.0046 v0 + 0))^2)Hhub/v0)/((((Hhub f/v0))^(5/3))(1+100.(Hhub f/v0)))),0])),{f,0.0001,100}, Frame-> True, FrameLabel-> {"f·PSD[f] = relative variance (m/s)^2 at frequency f","frequency f [Hz]"} ,GridLines->{Automatic,Automatic},GridLinesStyle->Directive[Yellow],PlotStyle->Thick, 
PlotRange ->{{0.0001,10},{0,2}},PlotLabel-> "The area beneath the f·PSD(f) plot in Log-Lin axis is the variance of the\nequivalent wind and the height is the relative contribution at each frequency f", Filling-> Axis ,Epilog->{ Inset[Style[Text["\!\(\*SubscriptBox[\(σ\), \(0.0001\\ Hz\\ ~\\ 100\\ Hz\)]\) =\n\!\(\*SqrtBox[\(shadowed\\ area\)]\)= "<>  ToString@NumberForm[√Quiet@NIntegrate[espectro[f],{f,0.0001 ,fblade,100},AccuracyGoal->6],{3,2},NumberPadding->{"","0"}]<> " m/s"],Larger,Bold],{0.65(* centered horizontally at 65% of plot*),Center}]} (* End of LogLinearPlot *)]
},Alignment->Center]

TabWindModel[]:=Manipulate[
(* Compute Spectrum *)
ComputeWindModel[];
(* Plot Graphics *)
PlotWindModel[],
{{TipoEspectro,1,"Wind spectrum type:"},{1->"Von Karman",2->"Kaimal",3->"Davenport" ,4-> "IEC 61400-1 class A (Iref = 16%)",5-> "IEC 61400-1 class B  (Iref = 14%)",6-> "IEC 61400-1 class C  (Iref = 12%)",7-> "Espectro recomendado por Consorcio Noruego del Petroleo (Slettringen)",17-> "Comparing Von Karman, Kaimal & IEC 61400",18->"Filtro equivalente del rotor"}, Alignment->Left },
{{v0,10,"Wind at hub height U∞:"Dynamic[Refresh[Plot[PowerCurve[v],{v,3,30},Frame-> False,Ticks-> None,GridLines-> None,Axes-> True,Epilog->{PointSize[Large],Point[{v0,PowerCurve[v0]}]},ImageSize->{40,40}], TrackedSymbols->{v0}]]},0,30,1,Appearance->"Labeled"},
{{Iref,14,"Intensity Turbulence at 15 m/s, \!\(\*SubscriptBox[\(I\), \(15\\ m/s\)]\)"},0.05,0.3,0.01,Appearance->"Labeled"},
{{Λ1,42,"Length (Subscript[Λ, 1]) ((Subscript[Λ, 1]) = 42 m in IEC 61400-1)"},30,70,1,Appearance->"Labeled"},
{{HornsRev,True,"Añadir oscilaciones de muy baja frecuencia (f < 1 mHz)\n observadas en P.E. Horns Rev (en Hovsore)"},{False,True},Checkbox},
{{Hhub,70,"Altura del buje"},25,100,0.5,Appearance->"Labeled",Enabled->HornsRev},
{{rotorfilter,2,"Rotor Filter"},{0-> "None",1-> "First Order (blade longitudinal average)",2-> "Second Order (area averaged, recommended)"},ControlType->Setter},{{R,41,"Rotor radius R:"},25,80,0.5,Appearance->"Labeled", Enabled->( rotorfilter>0)},
{{spatialsampling,False,"Include spatial sampling in equivalent wind (not reccomended) "},{True,False}},
{{GananciaPasoPala,0.5,"((σ^2)/(PSD ((Subscript[f, blade])))) due to spatial sampling:"},0.1,4,0.1,Appearance->"Labeled", Enabled-> spatialsampling},
{{fblade,1,"Blade frequency ( (Subscript[f, blade])= (Subscript[Ω, rotor])3/60):"},0.4,2,0.1,Appearance->"Labeled", Enabled-> spatialsampling},
{{BandWidth,0.05,"Relative bandwidth of rotor oscillations at (Subscript[f, blade]):"},0.001,0.2,0.002,Appearance->"Labeled", Enabled-> spatialsampling},
{{TowerShadowIntensity,0.05,"RMS oscillation: \!\(\*FractionBox[\(Typical\\ oscillation\\ amplitude\\ at\\ \*SubscriptBox[\(f\), \(blade\)]\), \(\*SqrtBox[\(2\)]\\ wind\\ speed\)]\)"},0,0.2,0.01,Appearance->"Labeled", Enabled-> spatialsampling},
SaveDefinitions-> True,TrackedSymbols:>Manipulate, LocalizeVariables->False,Paneled-> False,ContinuousAction->False (*Setting ContinuousAction to False updates only when the control is released:*),
SynchronousUpdating -> False (* Don't queue updates *),
ControlPlacement->{Top,Top,Top,Top,Top,Top,Top,Top,Bottom,Bottom,Bottom,Bottom,Bottom}, Alignment->Center ,(*Initialization:>( IniWindModel[];) ,*)SynchronousInitialization->True
] (* End of Manipulate *)

TabWindModel[]

TabWindResource

IniWindResource[FactorForma0_,FactorEscala0_ , v00_]:= Module[{},
pdfViento[v_]:=(E^(-((((v/FactorEscala)))^FactorForma))) (FactorEscala^(-FactorForma)) FactorForma (v^(-1+FactorForma));
cdfViento[v_]:=1-(E^(-((((v/FactorEscala)))^FactorForma)));cdfVientoEscalado[v_] :=cdfViento[v]pdfViento[FactorEscala (((((-1+FactorForma)/FactorForma)))^(1/FactorForma))] ; pdfVientoEscalado[v_] :=pdfViento[v]/pdfViento[FactorEscala (((((-1+FactorForma)/FactorForma)))^(1/FactorForma))] ;
FactorForma=FactorForma0;FactorEscala=FactorEscala0;v0=v00;
];

PlotWindResource[]:= Column[{Plot[{pdfVientoEscalado[v],cdfViento[v]  },{v,0,30},FrameLabel->{"Wind Speed (m/s)","Probability and relative probability density"},PlotLabel->"Relative Probability Density Function (PDF) and Cumulative Density Function (CDF)",ImageSize->{800,500}, Filling-> {1-> 0,2-> {2}}, Evaluate[EstiloGrafica[[2]]],
Epilog->{PointSize[Large],Point[{{v0,pdfVientoEscalado[v0]},{v0,cdfViento[v0] }}], Inset[Style[Text["Probability[wind ≤ "<>  ToString@v0<> " m/s ] = "<> ToString@NumberForm[N[100*(cdfViento[v0])],{3,2},NumberPadding->{"","0"}]<>"%"],Larger,Bold],{v0+4.5,cdfViento[v0]}],Inset[Style[Text["(Subscript[PDF, wind])[ "<> ToString@N@v0<> " m/s ] = "<>ToString@TraditionalForm@((" Probability["<>ToString@v0<>" m/s < wind ≤ "<>  ToString@v0<> " m/s + Δwind ]" )/ Δwind )<>"= \n= "<> ToString@N@pdfViento[v0]<> " (m/s)^(-1)"],Larger,Bold],{v0+7,pdfVientoEscalado[v0]}]}
](* End of plot *),
Text@TraditionalForm@Grid[{
{"(Subscript[μ, wind])"," == ",HoldForm[λ Γ(1+(1/k))]," == ", NumberForm[#,{5,4},NumberPadding->{"","0"}]&@ (vmedia =N[FactorEscala Gamma[1+(1/FactorForma)]])," m/s"},
{"(Subscript[σ, wind])"," == ",(Sqrt[HoldForm[(λ^2)Γ(1+(2/k))-((Subscript[μ, wind])^2)]]) ," == ", NumberForm[#,{5,4},NumberPadding->{"","0"}]&@(vstddev=N[(Sqrt[(FactorEscala^2)Gamma[1+(2/FactorForma)]-(vmedia^2)])])," m/s"}, 
{"wind median"," == ",HoldForm[((λ [Ln(2)])^(1/k))] ," == ", NumberForm[#,{5,4},NumberPadding->{"","0"}]&@(vmediana = N[FactorEscala ((Log[2])^(1/FactorForma))])," m/s"},
{"wind mode"," == ",HoldForm[λ (((1-(1/k)))^(1/k))] ," == ", NumberForm[#,{5,4},NumberPadding->{"","0"}]&@(vmoda = N[FactorEscala (((1-(1/FactorForma)))^(1/FactorForma))])," m/s"}
},Alignment->Center(*,ItemSize->{{10,3,12,3,10},Automatic}*)]}
,Background->None,Alignment->Center]

TabWindResource[FactorForma0_,FactorEscala0_ , v00_]:= 
Manipulate[PlotWindResource[],
{{FactorForma,FactorForma0,"Shape Factor Weibull Distr. k"},1.1,3,0.1,Appearance->"Labeled"},{{FactorEscala,FactorEscala0,"Scale Factor Weibull Distr. λ" },5,15,0.1,Appearance->"Labeled"},
(*{{v0,10,"⟨Wind⟩ [m/s]"},0,30,1,Appearance->"Labeled"},*)
{{v0,10,"Wind at hub height U∞:"Dynamic[Refresh[Plot[PowerCurve[v],{v,3,30},Frame-> False,Ticks-> None,GridLines-> None,Axes-> True,Epilog->{PointSize[Large],Point[{v0,PowerCurve[v0]}]},ImageSize->{40,40}], TrackedSymbols->{v0,FactorForma,FactorEscala}]]},0,30,1,Appearance->"Labeled"},
SaveDefinitions-> True,LocalizeVariables->False,Paneled-> False, TrackedSymbols:> Manipulate,
ContinuousAction->False (*Setting ContinuousAction to False updates only when the control is released:*),
SynchronousUpdating -> False (* Don't queue updates *),
(*Initialization:>(IniWindResource[FactorForma0,FactorEscala0, v00]),*)SynchronousInitialization->True]

TabWindResource[2,7((2 )/(Sqrt[π])),10]

TabGustStatistics

IniGustStatistics[]=Module[{},A=1;Δt =60./64.;v0 = 10; T=60;logΔt=0;joined = True; filled=True;(*espectro[f_]:=(((0.1^2)4 42/10 )/(((1+70.8(((f 42/10))^2)))^(5/6)));*) fmin = 1/600; fmax = 100;transformacion = 1;
RateString[DensityPerSecond_] =  If[DensityPerSecond> 1,  ToString@N[DensityPerSecond]<>" times per second",If[DensityPerSecond> 1/60,  ToString@N[60 DensityPerSecond]<>" times per minute",If[DensityPerSecond> 1/3600,  ToString@N[3600DensityPerSecond]<>" times per hour",  
If[DensityPerSecond> 1/(3600*24),   ToString@N[3600 24DensityPerSecond]<>" times per day",
If[DensityPerSecond> 1/(3600*24*30),   ToString@N[3600 24 30DensityPerSecond]<>" times per month",
If[DensityPerSecond> 1/(3600*24*365),   ToString@N[3600 24 365DensityPerSecond]<>" times per year",If[DensityPerSecond> 1/(3600*24*365*100),ToString@N[3600 24 365 100 DensityPerSecond]<>" times per century","--" (* Recurrence time above a century *)]]]]]]];
RecurrenceTimeString[DensityPerSecond_] =  
If[DensityPerSecond> 1/60,   ToString@N[1/(DensityPerSecond)]<>" seconds",If[DensityPerSecond> 1/(3600),   ToString@N[1/(60DensityPerSecond)]<>" minutes",  
If[DensityPerSecond> 1/(3600*24),   ToString@N[1/(3600DensityPerSecond)]<>" hours",
If[DensityPerSecond> 1/(3600*24*30),    ToString@N[1/(3600 24DensityPerSecond)]<>" days",
If[DensityPerSecond> 1/(3600*24*365),   ToString@N[1/(3600 24 30DensityPerSecond)]<>" months", 
If[DensityPerSecond> 1/(3600*24*365 100),   ToString@N[1/(3600 24 365DensityPerSecond)]<>" years",If[DensityPerSecond> 1/(3600*24*365 10000),  ToString@N[1/(3600 24 365 100 DensityPerSecond)]<>" centuries", "--" (* Recurrence time above 10000 years *)]]]]]]];

(* Wind Speed to compute mean values of gust occurences along a year *)
vmax = 35; vmin = 4;
pdfLaplace[μ_,σ_,x_]=((E^(-(((Sqrt[2]) (x-μ) Sign[x-μ])/σ)))/((Sqrt[2]) σ));
(*cdfLaplace[μ_,σ_,x_]= CDF[LaplaceDistribution[μ,σ],x];*)
cdfLaplace[μ_,σ_,x_]=(1/2) (1+(1-(E^(-(((Sqrt[2]) (x-μ) Sign[x-μ])/σ)))) Sign[x-μ]);
pdfNormal[μ_,σ_,x_]= PDF[NormalDistribution[μ,σ],x];
cdfNormal[μ_,σ_,x_]= CDF[NormalDistribution[μ,σ],x];
pdfRootNormal[μ_,σ_,x_] =(Sqrt[(Sqrt[3])((E^(-(Sqrt[3]) Abs[x-μ]/(σ^2)))/(8 (σ^2) π Abs[x-μ]))]);
cdfRootNormal[μ_,σ_,x_]=(1/2) (1+ Sign[x]Erf[(((3^(1/4)) (Sqrt[Abs[x]]))/((Sqrt[2]) σ))]);
cdfNormal[x_]=cdfNormal[0,1,x]; 
ProbabilityWindExceedsΔvGauss[Δv_,σ2_] =(1/2) (1-Erf[((Δv)/((Sqrt[2]) (Sqrt[σ2])))]);
ProbabilityWindExceedsΔvLaplace[Δv_,σ2_] =(1/2) (E^(-(((Sqrt[2]) Δv)/(Sqrt[σ2]))));
ProbabilityWindExceedsΔvSqrGauss[Δv_,σ2_] =(1/2) (1- Sign[Δv]Erf[(((3^(1/4)) (Sqrt[Abs[Δv]]))/((Sqrt[2])(Sqrt[σ2])))]);
];

UpdateModelIfNeeded[] := 
(* test constants to be defined *)(* If parameters are defined, then do computations if required*)
If[ NumericQ[ω1[[1]] σkaimal[v0,Iref]  espectro[freqs1[[2]]PSDUeq[[1]]npuntos2 Δt /T]],
(* Spectrum is defined and the model must be updated if needed *)
(* Recompute discrete spectrum if needed *)
If[ v00 ≠  v0,ComputeWindModel[]];
If[transformacion0==transformacion &&T0 == T && logΔt0 == logΔt&& espectro0[f]== espectro[f]&&A0==A ,
(* Do nothing, all variables are updated. Else, compute again discrete spectrum *),
DiscretizeSpectrum[];
MeanPeakGust[(*npuntos2,T*)] ;
ComputeGustStatistics[];
];
transformacion0=transformacion;T0 = T;logΔt0 =logΔt; espectro0[f_]= espectro[f];A0=A; v00 = v0;
(* Return true to indicate inicialization is ok *) True,

(* If wind parameters are not numeric, throw a message and return false *)
MessageDialog["The rest of wind spectra and wind parameters must be defined to be able to compute the characteristics of this gust"];
(* Return false to indicate inicialization has failed *) False
] (* End of UpdateModelIfNeeded[] *)

DiscretizeSpectrum[(*npuntos2,T, espectro*)] := Module[{},
(* Definiciones *)Δt=2^logΔt;  Δf = 1/T;
(* El número de frecuencias significativas del espectro es npuntos/2 puesto que las frecuencias por encima de fnyquist es el espectro reflejado. El mínimo número de puntos generados es 4 *)
npuntos2 = 2^Round@Min[ Max[1,Log[2,T/(2Δt)]],12 (* 2^12 = 4x1024 = 4096 *)];
(*npuntos2 está entre 2 y 4096 *)
 npuntos =2 npuntos2;
(* ajusta Δt y fnyquist después de limitar el número máximo de puntos de la muestra *)
 Δt =N[ T/npuntos]; fnyquist = npuntos2  Δf (* npuntos / (2 T) *) ;logΔt = Log[2,Δt];
(* Utiliza una semilla para cada muestra *)
freqs1 =N@Table[f,{f,Δf,fnyquist,Δf}];
ω1 = 2 π freqs1;
(* gain corresponds to the discrete value of 1-SIDED variance. Hence, the spectrum (S^2)[f] is divided by T (not 2 T), where Δf = 1/T and fnyquist = 1/(2Δt).
 (S^2)[f] is assumed to be mean-sqared units (effective values, not amplitude values). For example, Kaimal is allways expressed in effective values.
*)
PSDUeq=(espectro /@ freqs1);
(* avoid DC component *)
(* Variabilidad de la función contínua en el tiempo del viento eqivalente contínuo en el intervalo Δt de tiempo,
teniendo en cuenta 2-SIDED spectrum
(∫_t^(t+Δt))(((U[t] - ⟨U⟩ ))^2) \[DifferentialD]t == Total@Abs[ (spectrum^2)]/4 -(componentedc^2) == 
       Total@Abs[ (spectrum1^2)]/4 == Total[ ((Re[spectrum1])^2)+((Im[spectrum1])^2)]/4 ~ Total[ gain21 ] *)
(* El factor 4 T viene de;
  2 = ratio  (((amplitud/valor_rms))^2) -se puede justificar porque hay componente real e imaginaria => 1+1 = 2;
  2 = haber utilizado los coeficientes de la serie del coseno y seno. Si utilizo los coeficientes complejos de Fourier (two-sided spectrum) no hay q poner este segundo factor 2;
T => al ser una señal estocástica, la varianza de la señal es constante (la integral del cuadrado del spectro) pero magnitud de la FFT es proporcional a (Sqrt[4 T PSD+(f)]). Por tanto, hay que dividir el espectro por (Sqrt[4T]) para tener una medida espectral independiente de la longitud de la muestra.
(Sqrt[(1/(4T))PSD+(f)])
*)
(* Factor 1/2 = 1/(Sqrt[4])  due to the fact that the complex Fourier (two-sided spectrum) transform is: 
 sine:  1/2 I δ(f+f0) - 1/2 I δ(f-f0) 
or cosine: 1/2 δ(f-f0) + 1/2 δ(f-f0)
*)
(* Factor 1/(Sqrt[4]) already included in (Sqrt[gain]) *)
(*factor Δf=1/(4T) para escalar el PSD a contenido frecuencial en la banda Δf*)gain21=PSDUeq/(4 T);
(*std deviation is the square root of the variance*)
stDevSpectrum=Sqrt[gain21];
];

MeanPeakGust[(*npuntos2,T*)] := Module[{},
TransformedProcess[x_]=Switch[transformacion,
1,x,
2,v0-(Sqrt[((σinterval2)/2)]) Log[Erfc[((Abs[x-v0])/(Sqrt[2 σinterval2]))]] Sign[x-v0],
3,v0+(*Sign[x-v0] (Sqrt[Abs[x-v0]])*)(( Sign[x-v0] (((x-v0))^2))/(Sqrt[3σinterval2])) (* Transformación raiz comprobada *),
 _, x (* default value *)];
TransformedProcess[x_]=Switch[transformacion,
1,x,
2,v0-(Sqrt[((σinterval2)/2)]) Log[Erfc[((Abs[x-v0])/(Sqrt[2 σinterval2]))]] Sign[x-v0],
3,v0+(*Sign[x-v0] (Sqrt[Abs[x-v0]])*)(( Sign[x-v0] (((x-v0))^2))/(Sqrt[3σinterval2])) (* Transformación raiz comprobada *),
 _, x (* default value *)];
InverseTransformedProcess[x_]=Switch[transformacion,
1,x,
2,v0+(Sqrt[2 σinterval2]) Sign[x-v0]InverseErfc[Exp[-((Abs[x-v0])/(Sqrt[σinterval2/2]))]],
3,v0+Sign[x-v0] (Sqrt[(Sqrt[3σinterval2])Abs[x-v0]]) (* Transformación raiz comprobada *),
 _, x (* default value *)];
 WindDistributionName =Switch[ transformacion,
1,"normal",
2,"laplacian",
3,"signed square of a normal",
 _, "normal" (* default value *)];
AbsAccelDistributionName =Switch[ transformacion,
1,"half-normal",
2,"exponential",
3,"chi-square with one degree of freedom",
 _, "normal" (* default value *)];

(* COEFFICIENTS COMPUTED FROM DISCRE SYSTEM  *)
ΣS=(Total@gain21);
ωS = (ω1  gain21);
ΣωS = Total@ωS ;
ω2S = ω1 ωS ;
Σω2S =  Total@ω2S ;
(*Σω4S = ((ω1)^2) ω2S*) ;
halfcosiaverage = gain21(((InverseTransformedProcess[A+v0] -v0)/2)/(ΣS));
quarterofcosivar=gain21(1-gain21/ΣS );
halfcosistdev =(Sqrt[quarterofcosivar]);
 quarterofsinivar =gain21(1-((ω2S)/(Σω2S)));
halfsinistdev =  (Sqrt[quarterofsinivar]);

(* Average shape of the gust. 1/2 factor to transform cosine and sine 1-sided Fourier coefs into 2-sided complex Fourier coefs *)
(* Factor 1/2  due to the fact that the complex Fourier (two-sided spectrum) transform is HALF the one-sided spectrum: 
 sine:  1/2 I δ(f+f0) - 1/2 I δ(f-f0) 
or cosine: 1/2 δ(f-f0) + 1/2 δ(f-f0) 
*)
averagespectrum =(Sqrt[npuntos]) Join[{N[v0 ]},halfcosiaverage,(*Conjugate@*)Reverse@Take[halfcosiaverage,npuntos2-1]];
discreteACFwind = TransformedProcess@ RotateLeft[Re@InverseFourier[averagespectrum],npuntos2];
(* Add the point correspondant to t=T as the same point as time t=0 since time series is periodic *)
AppendTo[discreteACFwind,discreteACFwind[[1]]];
(* 1/Factor == ((( coefs.Sine-Cosine))^2) -> (((1))^2) ) ((((Sine-Cosine Amplitude/RMS))^2) -> 2, or alternatively, (Re^2)+(Im^2)-> 2) == 2
Note that only positive frequencies are in spectrum  (1 sided PSD): Factor == (1/((((1))^2) 2))-> 1/2)
*)
(* Autocorrelation function and its inverse *)
ACF=ListInterpolation[Reverse@Take[discreteACFwind,(npuntos2+1)],{{0,T/2 }},InterpolationOrder->1];
      invACF=ListInterpolation[ Range[0,T/2,Δt],{Reverse@Take[discreteACFwind,(npuntos2+1)]},InterpolationOrder->1];
LimitForGustDuration = A/E+v0;
ExpectedNormalizedGustDuration = 2 invACF[LimitForGustDuration];

timeseriestdev=√( Total@quarterofcosivar + Total@quarterofsinivar);

(* COEFFICIENTS COMPUTED FOR ESTIMATION OF NUMBER OF GUSTS OF CONTINUOUS SYSTEM *)
(* A CUT-OFF FREQUENCY IS EMPLOYED TO AVOID COUNTING BURST LEVEL CROSSINGS AS MANY GUSTS *)
(* THE CUT-OFF TIME IS THE EXPECTED NORMALIZED GUST DURATION, THIS DOES NOT MODIFY SIGNIFICANTLY THE VARIANCE OF THE SIGNAL NOR THE GUST SHAPE BUT IT SMOOTHS THE DERIVATIVE OF THE WIND, AVOIDING VERY TINY MAXIMUMS AT SOUND FREQUENCIES TO BE COUNTED AS GUSTS -fractal behaviour of the wind- *)
(* Assuming the time that the gust is above Δv/E is half period, the average gust cycle is twice the normalized gust duration *)
(* For not overfiltering shorter gust, a compromise between over or infra estimating gust is just allow gust up to half the average duration *) 
ffiltered = 1/( (*2*) ExpectedNormalizedGustDuration);

];

MeanRampGust[] := Module[{},
(* ffrequency goes from Δf to fnyquist. Spectrum is reflected around its middle (negative frequency components) for real data.*)
unfilteredComponents = Round[ffiltered / Δf];
(* Part of the spectrum up to ffiltered *)
freqs1filtered = freqs1[[1;;unfilteredComponents]];
gain21filtered = gain21[[1;;unfilteredComponents]];
ω2Sfiltered = ω2S[[1;;unfilteredComponents]];
(* In case very few elements are consideded (average gust time a fraction of the interval), the rest of the spectrum is zero padded to increase the resolution *)
If[unfilteredComponents<16,ω2Sfiltered =PadRight[ ω2Sfiltered, 16];unfilteredComponents=16;];
(* npoints in the filtered ACFslope = 2 * unfilteredComponents *)
WindSlopeSpectrum =(Sqrt[unfilteredComponents (* total number of points in spectrum *)]) Join[{0} (*dc*),ω2Sfiltered,(*Conjugate@*)Reverse@Take[ω2Sfiltered,unfilteredComponents-1]];
(* Remove frequencies above fcutoff to avoid false event counting due to signal rugosity *)
discreteACFwindSlope = Re@InverseFourier[WindSlopeSpectrum];
(* Autocorrelation function and its inverse *)
ACFslope=ListInterpolation[Take[discreteACFwindSlope,unfilteredComponents+1],{{0,T/2 }},InterpolationOrder->1];invACFslope=ListInterpolation[ Range[0,T/2,T/(2unfilteredComponents)],{Take[discreteACFwindSlope,unfilteredComponents+1]},InterpolationOrder->1];
];

ComputeGustStatistics[] := Module[{},
σinterstep2 :=Quiet@NIntegrate[espectro[f],
Evaluate@If[fnyquist<fblade,{f,fnyquist ,fblade,fmax},{f,fnyquist ,fmax}],AccuracyGoal->4];
(* Varianza entre dos puntos consecutivo ⟨(((U[t+Δt]-U[t]))^2)⟩ de viento eqivalente debido a todo el espectro del viento
σinterstep2 :=0.5 NIntegrate[espectro[f] (1-Cos[2 π f Δt]),
Evaluate@If[ Δf<fblade,{f, Δf ,fblade,fmax},{f, Δf ,fmax}],AccuracyGoal->8] ;*)
(*σinterstep2 =NIntegrate[((espectro[f])^(1/2))/2,{f,fnyquist ,fmax},AccuracyGoal-> 3] ^2;*)
σdiezminutal2=Quiet@NIntegrate[espectro[f],Evaluate@If[fmin<fblade<fnyquist,
{f,fmin ,fblade,fnyquist},{f,fmin ,fnyquist}],AccuracyGoal->4] ;
σinterval2= Quiet@NIntegrate[espectro[f],Evaluate@If[Δf<fblade<fnyquist,
{f,Δf,fblade,fnyquist},{f,Δf ,fnyquist}],AccuracyGoal->4] ;
ΣSfiltered =  NIntegrate[ espectro[f], Evaluate@If[ffiltered<fblade, {f,1/T,ffiltered}, {f,1/T,fblade,ffiltered}]];
ΣωSfiltered = NIntegrate[(2 π f) espectro[f],Evaluate@If[ffiltered<fblade, {f,1/T,ffiltered}, {f,1/T,fblade,ffiltered},AccuracyGoal->4]];
Σω2Sfiltered = NIntegrate[(((2 π f))^2) espectro[f],Evaluate@If[ffiltered<fblade, {f,1/T,ffiltered}, {f,1/T,fblade,ffiltered},AccuracyGoal->4]];
Σω3Sfiltered = NIntegrate[(((2 π f))^3) espectro[f],Evaluate@If[ffiltered<fblade, {f,1/T,ffiltered}, {f,1/T,fblade,ffiltered},AccuracyGoal->4]];
Σω4Sfiltered = NIntegrate[(((2 π f))^4) espectro[f], Evaluate@If[ffiltered<fblade, {f,1/T,ffiltered}, {f,1/T,fblade,ffiltered},AccuracyGoal->4]];

ΣSunfiltered =  NIntegrate[ espectro[f], Evaluate@If[fmax<fblade, {f,1/T,fmax}, {f,1/T,fblade,fmax},AccuracyGoal->4]];
ΣωSunfiltered= NIntegrate[(2 π f) espectro[f],Evaluate@If[ffiltered<fblade, Evaluate@If[fmax<fblade, {f,1/T,fmax}, {f,1/T,fblade,fmax},AccuracyGoal->4], {f,1/T,fblade,ffiltered}]];
Σω2Sunfiltered= NIntegrate[(((2 π f))^2) espectro[f],Evaluate@If[ffiltered<fblade, Evaluate@If[fmax<fblade, {f,1/T,fmax}, {f,1/T,fblade,fmax},AccuracyGoal->4], {f,1/T,fblade,ffiltered}]];
Σω4Sunfiltered= NIntegrate[(((2 π f))^4) espectro[f],Evaluate@If[ffiltered<fblade, Evaluate@If[fmax<fblade, {f,1/T,fmax}, {f,1/T,fblade,fmax},AccuracyGoal->4], {f,1/T,fblade,ffiltered}]];Σω4Sfiltered = NIntegrate[(((2 π f))^4) espectro[f], Evaluate@If[ffiltered<fblade, {f,1/T,ffiltered}, {f,1/T,fblade,ffiltered},AccuracyGoal->4]];
ΔUmaxSample =TransformedProcess[v0+ √ σinterval2 Log[((((T/(2 π))))^2)((Σω2Sfiltered (* = σacel2*))/(σinterval2))]]-v0;
meanAcelerationFactor=Switch[transformacion,
1,(Sqrt[2/π]),
2,1/(Sqrt[2]),
3,1,
 _,(Sqrt[2/π])(* default value *)];
pdfWind[μ_,σ_,x_]=Switch[transformacion,
1,pdfNormal[(*μ*) μ ,(* (Subscript[σ, y]) *)σ ,x],
2,pdfLaplace[(*μ*) μ ,(* (Subscript[σ, y]) *)σ ,x],
3,pdfRootNormal[(*μ*) μ ,(* (Subscript[σ, y]) *)σ ,x] (* Transformación raiz comprobada *),
 _, pdfNormal[(*μ*) μ ,(* (Subscript[σ, y]) *)σ ,x] (* default value *)];
ProbabilityWindExceedsΔv[Δv_,σ2_] =Switch[transformacion,
1,ProbabilityWindExceedsΔvGauss[Δv,σ2],
2,ProbabilityWindExceedsΔvLaplace[Δv,σ2],
3,ProbabilityWindExceedsΔvSqrGauss[Δv,σ2] (* Transformación raiz comprobada *),
 _, ProbabilityWindExceedsΔvGauss[Δv,σ2](* default value *)];
(*"Moe and Niedzwecki, \"Frequency of maxima of non-narrow banded stochastic\nprocesses\"Applied Ocean Research 27 (2005) 265-272 ((Subscript[f, cut - off]) =(1/MeanGustDuration)).\nRISO-R-1437, page 35, Eq. 4.12 ((Subscript[f, cut - off]) =\!\(\*FractionBox[\(1\), \(\(\\ \)\(MeanGustDuration\)\)]\))\n" *)(*Similar to Bierbooms, EWEA'07, page 13, Eq. D.19 except for the factor ((Σω4Sfiltered)/(2 π ΣSfiltered Σω2Sfiltered))*)
(*RiceDensityOfLocalMaximaInANormalProcess*)
PeakGustRate[ΔVc_,ΣSunfiltered_,Σω2Sfiltered_]=(1/2)meanAcelerationFactor * (* (Subscript[σ, (Overscript[y, *])]) *)(Sqrt[Σω2Sfiltered]) *pdfWind[(*μ*) 0 ,(* (Subscript[σ, y]) *)(Sqrt[ ΣSunfiltered]) ,ΔVc];

GustDensityPerSecond = PeakGustRate[A,ΣSunfiltered,Σω2Sfiltered];
DensityWindExceedsΔvPerYear[Δv_,vmin_,vmax_] := NIntegrate[pdfViento[v] PeakGustRate[A,(((((σkaimal[v,Iref])/(σkaimal[v0,Iref]))))^2)ΣSunfiltered,(((((σkaimal[v,Iref])/(σkaimal[v0,Iref]))))^2)Σω2Sfiltered],{v,vmin,vmax}, AccuracyGoal-> 5];
YearlyAveragedGustDensityPerSecond =DensityWindExceedsΔvPerYear[A,vmin,vmax];
];

TextGustStatistics[]:=Column[{Text@TraditionalForm["* The standard deviation of wind in an interval of "<> ToString@T<>" s is (Subscript[σ, T]) = "<> ToString@(Sqrt[σinterval2])<>" m/s ((Subscript[σ, T]) is the typical difference of instantaneous wind respect the interval mean wind). \n* The standard deviation of wind in an interval of 10 minutes  is (Subscript[σ, 10']) = "<>ToString@(Sqrt[ σdiezminutal2])<>" m/s.\n* The deviation from the interval average ΔU has "<> WindDistributionName<>" distribution and the probability of exceding ΔU = +"<>ToString@A<>" m/s (positive peak) is "<>ToString[100*ProbabilityWindExceedsΔv[A,(* (Subscript[σ, y]) *)(Sqrt[σinterval2])] ]<>" %."<>("\n* According to the wind spectra, the mean gust duration ⟨τ⟩ is "<> ToString@Round@ExpectedNormalizedGustDuration<>" seconds on average (mean time above ⟨U⟩+Δv/E, limit drawn in the graph ).\n* The average gust shape and duration is independent of the gust peak Δv (and the period lengh T provided T >> ⟨τ⟩), except for a small distortion from the non-linear transformation.\n* The maximum wind expected in an interval of "<> ToString@T<>" s is ⟨(Subscript[ΔU, Umax, T])⟩ = "<> ToString@ΔUmaxSample<>" m/s.\n")<>("\n* The wind U is assumed to have "<>WindDistributionName<>" distribution and hence the absolute value of the air acceleration, |(d/dt)U|, has "<>AbsAccelDistributionName<>" distribution.\n* Gust statistics are computed counting as one gust those local maxima closer than one average gust duration." )<>("* When the (equivalent) wind is ⟨U⟩ ≈ "<> ToString@v0<>" m/s (average measured on "<> ToString@T<>" s intervals), a peak gust exceding the interval mean in ΔU ≥" <> ToString[A]<>" m/s happens " <> RateString[GustDensityPerSecond] 
<>" on average ( equivalent to a gust approximately every "<>RecurrenceTimeString[GustDensityPerSecond] <>" on average).")<>("\n* Taking into account the long-term distribution of wind corresponding to a Weibull( λ = "<> ToString[NumberForm[#,{4,1}]&@N@FactorForma ]<>", k = "<>ToString[NumberForm[#,{4,1}]&@N@FactorEscala ] <> "  m/s) and a distribution of turbulence according the standard IEC 61400-1, the overall average of the number of gust excedding the interval mean in ΔU = " <> ToString@A<>" m/s is "<>RateString[YearlyAveragedGustDensityPerSecond] <>" or one gust every ")<>RecurrenceTimeString[YearlyAveragedGustDensityPerSecond] <>" on average."],
TraditionalForm@Text@TableForm[{
{RecurrenceTimeString[ProbabilityWindExceedsΔv[A/E,σinterval2] (1/ExpectedNormalizedGustDuration)]},
{RecurrenceTimeString[GustDensityPerSecond]}},TableHeadings->{{"Recurrence time ≈ \!\(\*FractionBox[\(\(Mean\)\(\\ \)\(Gust\)\(\\ \)\(Duration\)\(\\ \)\), \(Probability\\ Exceeding\\ Δv/E\)]\) = ",
"Gust Rate = RateUpCrossings[ΔVc] = (1/2)⟨|(d/dt)Δv|⟩(Subscript[PDF, Δv])[ΔVc] ≈ "},{"τ, mean time between gusts"}}, TableAlignments-> {Center,Center},TableSpacing->{3,5}],
(*"\n* The plot beneath represents the mean shape of a peak gust (in the central dash dot line) and the upper and lower dashed line are the ±(Subscript[σ, gust]) limits. A random peak gust has been added in solid thick line."*)
Quiet@ListLinePlot[discreteACFwind,DataRange->{0,T},Epilog->{{Dashed,Line[{{(T-ExpectedNormalizedGustDuration)/2,LimitForGustDuration},{(T+ExpectedNormalizedGustDuration)/2,LimitForGustDuration}}]},{Dashed,Arrow[{{(T-ExpectedNormalizedGustDuration)/2,LimitForGustDuration},{(T-ExpectedNormalizedGustDuration)/2,v0}}]},{Dashed, Arrow[{{(T+ExpectedNormalizedGustDuration)/2,LimitForGustDuration},{(T+ExpectedNormalizedGustDuration)/2,v0}}]},{Dotted,Line[{{0,v0},{T,v0}}]},Inset[Style[Text["⟨Uwind⟩"],Larger,Bold],{T/10,v0+A/25}],Inset[Style[Text["Mean gust\nduration: "<>  ToString@Round@ExpectedNormalizedGustDuration<>" s"],Larger,Bold],{T/2,v0+A/E-0.6 timeseriestdev}],
Inset[Style[Text["⟨Gust⟩"],Larger,Bold],{0.07 T,discreteACFwind[[1]] - 0.1 timeseriestdev+0.1 A}]
},
ImageSize->tamano,(*ImagePadding->{{25,25},{25,25}},*) Joined->joined,Filling->If[filled,{1->v0},None], Frame-> True, FrameLabel-> {"time t [s]","wind [m/s]"} ,GridLines->{Automatic,Automatic}, GridLinesStyle->Directive[Yellow],PlotLabel->"Average shape of a peak gust",PlotStyle->{Thick} , InterpolationOrder-> If[npuntos>60,1,2], PlotRange->{{0,T},{v0-A,v0+A}}
]

}, Alignment -> Center](*End of Column *)

TabGustStatistics[]:=Manipulate[ 
(* If wind model is numeric and updated, then show tab *)
If[UpdateModelIfNeeded[],
(* Write main properties of the process *)
TextGustStatistics[]
] (*End of if *),
{{transformacion,1,"Transform of the Gaussian process:"},{1->"None (Gaussian distribution)",2->"Laplacian Distribution",3->"±(X^2) Distribution (X is Gaussian)"}, Alignment->Left },{{A,1,"Gust peak Δv measured from interval mean  [m/s]"},0.25,5,0.25,Appearance->"Labeled"},
{{v0,10,"⟨Wind⟩ [m/s]"},0,30,0.5,Appearance->"Labeled"},
{{T,60,"Time [s]"},10,600,1,Appearance->"Labeled"},
{{logΔt,60/64,Dynamic@Refresh["Δt [s]: "<>ToString[2^logΔt],TrackedSymbols:>{logΔt}] },Log[2,0.01],Log[2,100.],0.1},

SynchronousUpdating -> False (* Don't queue updates *),
SaveDefinitions-> True,TrackedSymbols:>{TipoEspectro,Iref,Λ1,rotorfilter,R,spatialsampling,GananciaPasoPala,fblade,BandWidth,TowerShadowIntensity,espectro,transformacion, v0,T,logΔt}, LocalizeVariables->False,Paneled-> False, (*Initialization->(IniGustStatistics[];) ,*)SynchronousInitialization->True ,ControlPlacement->{Top,Top,Top,Top,Top,Bottom,Bottom,Bottom}, Alignment->Center ]

TabGustStatistics[]

TabRandomSample

Tomo n puntos :  t = 0, Δt, 2 Δt, ... , (npuntos - 1) Δt. Por tanto, el último punto es t =  (npuntos - 1) Δt.
El periodo es T = npuntos Δt (la señal se replica a la siguiente muestra de npuntos).

IniRandomSample[] := Module[{},
];

ComputeRandomSample[(*npuntos2,T, espectro*)] := Module[{},
SeedRandom[GeneratedSample];
(* npuntos/2 Componentes aleatorias del espectro *)
(* parteReal = twice the cosine amplitude at each frequency *)
(* Generate samples in a way that changing the time step, the main waveform does not change: i.e., low part of the spectrum doesn't vary *)
(* normalized Gauss Distribution *)
(* Generation of random sample in order not to vary with varying timesteps Δt*)
NormalSample =Partition[RandomReal[NormalDistribution[] (* normalized Gauss Distribution *),npuntos],2]; 
parteReal =NormalSample[[All,1]];
(* parteImaginaria = minus twice the sine amplitude at each frequency *)
parteImaginaria =NormalSample[[All,2]];

spectrum1 = stDevSpectrum( parteReal+I parteImaginaria);
(* Comprobación: Genera el espectro el método Shinozuka: el módulo corresponde a (Sqrt[2]) stDevSpectrum y argumento aleatorio 
Este espectro genera una desviación estándar constante y muy similar a la teórica (no coincide por la discretización)
spectrum1  =(Sqrt[2]) stDevSpectrum Exp[I 20 parteReal];*)
(*Complete spectrum with deterministic DC component
PrependTo[parteReal,v0  (Sqrt[npuntos]) ];
PrependTo[parteImaginaria,0.0];*)
reflectedspectrum =Join[{N[v0  ]},spectrum1, Conjugate@Reverse@Take[spectrum1,npuntos2-1]];
(* To avoid zero division, add DC gain 
gain2=Prepend[gain21,1.0];
freqs=Prepend[freqs1,0.0];*)
(* Due to Mathematica default constants in Fourier Transforms, a factor (Sqrt[npuntos]) must be added to the reflected spectrum*)
(* Si se utiliza el convenio de transformada típico de teoría de la señal para la InverseFourier, hay que multiplicar por EL NÚMERO DE PUNTOS  en vez de su raíz cuadrada, (Sqrt[npuntos])*)
 timeseriesRandomSample =  TransformedProcess[Re@InverseFourier[(Sqrt[npuntos])reflectedspectrum]];
SampleVariance = Variance@timeseriesRandomSample;
(*Since interpolation is periodic, wind at t = 0 is the same as time t = T. Otherwise, time series span is T-Δt *)
transformedSpectrum=Fourier[timeseriesRandomSample][[2;;npuntos2+1]]/(Sqrt[npuntos]);
AppendTo[timeseriesRandomSample,timeseriesRandomSample[[1]]];];

PlotRandomSample[]:=GraphicsColumn[{If[(σinterstep2>0) && (σinterval2>0) ,
Quiet@ListLinePlot[{timeseriesRandomSample,timeseriesRandomSample-2(Sqrt[σinterstep2]),
timeseriesRandomSample+2(Sqrt[σinterstep2]),{v0-(Sqrt[σinterval2]),v0-(Sqrt[σinterval2])},{v0+(Sqrt[σinterval2]),v0+(Sqrt[σinterval2])},
{v0-2(Sqrt[σinterval2]),v0-2(Sqrt[σinterval2])},{v0+2(Sqrt[σinterval2]),v0+2(Sqrt[σinterval2])},
{v0-ΔUmaxSample,v0-ΔUmaxSample},{v0+ΔUmaxSample,v0+ΔUmaxSample}},
(*AxesOrigin->{0,0},*)DataRange->{0,T},
ImageSize->tamano,(*ImagePadding->{{25,25},{25,25}},*)Joined->joined,
Filling->If[filled,{2->{1},3-> {1},4-> {5},6-> {7}, 8-> {9}},None],Frame-> True,
 FrameLabel-> {{"wind [m/s]",""},{"time t [s]",""}} ,GridLines->{Automatic,Automatic},
GridLinesStyle->Directive[Yellow],PlotLabel->"Synthetized wind series",
PlotStyle->{
(*Thick*) Thickness@Min[0.01,Max[(*Grosor como una fracción del tamaño del dibujo *)(0.4/6)(Sqrt[((σinterstep2)/( σinterval2))]),0.001]]
 ,Thin,Thin, Dashed,Dashed, Dotted,Dotted, Dotted,Dotted}, InterpolationOrder-> If[npuntos>40,1,2] ,
 PlotRange-> If[fittedrangeplot, {Min[timeseriesRandomSample],Max[timeseriesRandomSample]},
{v0-Max[3(Sqrt[ σdiezminutal2]),2], v0+Max[3(Sqrt[ σdiezminutal2]),2]}],
Epilog->{Inset[Style[Text["⟨Wind⟩ + \!\(\*SubscriptBox[\(σ\), \(wind\\ in\\ interval\\ T\)]\)"],Larger,Bold],
{(T/2),v0+(Sqrt[σinterval2])+0.15(Sqrt[σdiezminutal2])}],
Inset[Style[Text["⟨Wind⟩ - \!\(\*SubscriptBox[\(σ\), \(wind\\ in\\ interval\\ T\)]\)"],Larger,Bold],
{(T/2),v0-(Sqrt[σinterval2])-0.15(Sqrt[σdiezminutal2])}],Inset[Style[Text["⟨Wind⟩ + 2 \!\(\*SubscriptBox[\(σ\), \(wind\\ in\\ interval\\ T\)]\)"],Larger,Bold],
{(T/2),v0+2(Sqrt[σinterval2])+0.15(Sqrt[σdiezminutal2])}],Inset[Style[Text["⟨Wind⟩ - 2 \!\(\*SubscriptBox[\(σ\), \(wind\\ in\\ interval\\ T\)]\)"],Larger,Bold],
{(T/2),v0-2 (Sqrt[σinterval2])-0.15(Sqrt[σdiezminutal2])}],Inset[Style[Text["⟨(Subscript[U, max])⟩"],Larger,Bold],
{(T/5),v0+ΔUmaxSample}],Inset[Style[Text["⟨(Subscript[U, min])⟩"],Larger,Bold],
{(T/5),v0-ΔUmaxSample}]
}
],Quiet@ListLinePlot[{timeseriesRandomSample},DataRange->{0,T},
ImageSize->tamano,Joined->joined,
Frame-> True, FrameLabel-> {"wind [m/s]","time t [s]"} ,GridLines->{Automatic,Automatic},GridLinesStyle->Directive[Yellow],PlotLabel->"Synthetized wind series",PlotStyle->{Thick}, InterpolationOrder-> If[npuntos>40,1,2] , PlotRange-> If[fittedrangeplot, {Min[timeseriesRandomSample],Max[timeseriesRandomSample]},{v0-Max[3(Sqrt[ σdiezminutal2]),2], v0+Max[3(Sqrt[ σdiezminutal2]),2]}]]],
Quiet@ListLogLinearPlot[(* Use factors to transform discrete fourier complex values to PSD *) 
{
Transpose@{freqs1,aux = freqs1 PSDUeq},Transpose@{freqs1,2T  freqs1( ((Re[transformedSpectrum])^2)+ ((Im[transformedSpectrum])^2))}
(* El factor 2 viene de haber utilizado los coefs complejos de Fourier en 1 sided PSD*)
(* 1/Factor == ( (((coefs complejos / coefs.Sine-Cosine))^2) -> (((1/2))^2) ) ((((Sine-Cosine Amplitude/RMS))^2) -> 2, or alternatively, (Re^2)+(Im^2)-> 2) == 2
Note that only positive frequencies are in spectrum  (1 sided PSD): 
2/(((1/2))^2) -> Factor == 2/(((1/2))^2)==2
 *)
(* Spectrum1 are complex coefs and have 1/2 factor compared to sin~cos coefs *)
(* Abs[spectrum1(f)] ~ (1/2) (1/(Sqrt[2]))(Sqrt[((PSD(f))/T)])= (Sqrt[((PSD(f))/(8 T))]) ==> PSD(f) ~ 2 T Abs[spectrum1(f)]; 
Re[spectrum1(f)] ~ (1/2) (Sqrt[((PSD(f))/T)])= (Sqrt[((PSD(f))/(4 T))]) ==> PSD(f) ~ 4 T ((Re[spectrum1(f)])^2);
Im[spectrum1(f)] ~ (1/2) (Sqrt[((PSD(f))/T)])= (Sqrt[((PSD(f))/(4 T))]) ==> PSD(f) ~ 4 T ((Im[spectrum1(f)])^2);*)
}, Frame-> True, 
ImageSize->tamano,Joined->joined,FrameLabel-> {"f·PSD[f] = frequency power content (m/s)^2","frequency f [Hz]"} ,
GridLines->{Automatic,Automatic},GridLinesStyle->Directive[Yellow],PlotStyle->Thick,
 PlotLabel-> "Frequency content of the transformed sample vs. (Subscript[PSD, Ueq])(f)", Filling-> Axis ,
Epilog->{ Inset[Style[Text["(Subscript[(OverHat[s]), sample])=\!\(\*SqrtBox[\(sum\\ of\\ shadowed\\ area\)]\)=\n= "<> 
 ToString@NumberForm[(Sqrt[SampleVariance]),{3,2},NumberPadding->{"","0"}]<> " m/s"],
Larger,Bold],{0.65(* centered horizontally at 65% of plot*),Center}]} , PlotRange->{{ freqs1[[1]],10},{0,(*(Sqrt[SampleVariance])*)aux [[1]]  }}
 (* End of LogLinearPlot *)]
}, Alignment -> Center](*End of GraphicsColumn *)

TabRandomSample[]:=Manipulate[
(* If wind model is numeric and updated, then show tab *)
If[UpdateModelIfNeeded[],
(* Compute and plot random sample *)
ComputeRandomSample[];
PlotRandomSample[]
] (*End of if *),
{{transformacion,1,"Transform of the Gaussian process:"},{1->"None (Gaussian distribution)",2->"Laplacian Distribution",3->"±(X^2) Distribution (X is Gaussian)"}, Alignment->Left },
(*{{v0,10,"⟨Wind⟩ [m/s]"},0,30,0.5,Appearance->"Labeled"},*)
{{T,60,"Time [s]"},10,600,1,Appearance->"Labeled"},
{{logΔt,60/64,Dynamic@Refresh["Δt [s]: "<>ToString[2^logΔt],TrackedSymbols:>{logΔt}] },Log[2,0.01],Log[2,100.],0.1},
{{GeneratedSample,1,"Generated sample (generate new sample\nto update settings modified on other tabs):"},1,10,1,Appearance->"Labeled"},
{{joined,True,"joined"},{True,False}},
{{filled,True,"filling"},{True,False}},
{{fittedrangeplot,False,"Fit vertical scale to sample"},{True,False}},
ContinuousAction->False (*Setting ContinuousAction to False updates only when the control is released:*), SynchronousUpdating -> False (* Don't queue updates *),
SaveDefinitions-> True,TrackedSymbols:>{espectro,transformacion, v0,T,logΔt,GeneratedSample,joined,filled,fittedrangeplot}, LocalizeVariables->False,Paneled-> False, (*Initialization->(IniRandomSample[];)*)SynchronousInitialization->True ,ControlPlacement->{Top,Top,Top,Top,Bottom,Bottom,Bottom,Bottom}, Alignment->Center ]

TabPeakGust

IniPeakGust[]  := Module[{},
];

ComputePeakGust[]  := Module[{},
(* gain21 and other staff is assumed to be formerly computed *)
SeedRandom[GeneratedSamplePeakGust];
(* parteReal = twice the cosine amplitude at each frequency *)
(* Generate samples in a way that changing the time step, the main waveform does not change: i.e., low part of the spectrum doesn't vary *)
(* normalized Gauss Distribution *)
NormalSample =Partition[RandomReal[NormalDistribution[] (* normalized Gauss Distribution *),npuntos],2]; 
parteReal =NormalSample[[All,1]];
(* parteImaginaria = minus twice the sine amplitude at each frequency *)
parteImaginaria =NormalSample[[All,2]];
(* avoid DC component *)
(*SiRe =(Sqrt[ gain21]) parteReal / T;*)
SiRe =halfcosiaverage+halfcosistdev  parteReal ;
ΣSiRe=Total[SiRe ];  (* wind at t=0, its mean is ⟨Δv[t=0]⟩ = A *)
(*SiIm =(Sqrt[ gain21])parteImaginaria / T;*)
SiIm =halfsinistdev parteImaginaria;
(* ωiSiIm =  ω SiIm //N;*)
ΣωiSiIm = ω1. SiIm; (* derivative of wind at t=0, it has zero mean *)
(*Σωi2SiIm = (((2 π))^2)Total[  ((freqs1^2) SiIm)] //N;*)

(* solución en media A -pero no tiene el valor exacto A en t=0-,spectrum1 =0.5(halfcosiaverage + parteReal halfcosistdev+I parteImaginaria halfsinistdev);*)
(* Factor 1/2 due to the fact that the Fourier transform is:
 sine:  1/2 I δ(f+f0) - 1/2 I δ(f-f0) 
or cosine: 1/2 δ(f-f0) + 1/2 δ(f-f0) *)
(* Spectrum constrained to Δv = A and dv/dt = 0 at t=0 *)
spectrum1 =0.5(SiRe+gain21((A-ΣSiRe)/(ΣS))+I (SiIm -ωS((ΣωiSiIm)/(Σω2S))));
(* Impose second derivative negative to be a maximum. In fact, generate sample whithout impose second derivative and then reject sample if second derivative is not negative -maximum- *)
(* If high frequency is present in the spectrum, this is not significative since wind or signal is very sharp, almost fractal in nature. This can be avoided computing derivatives using a cut-off frecuency such as ffiltered *) 
(* Usually, second deivative is negative since value at t=0 is avobe mean and the signal tends to return to the mean *)
(* Factor 2 due to the fact that the Fourier transform of sine/cosine is:
 sine:  1/2 I δ(f+f0) - 1/2 I δ(f-f0) 
or cosine: 1/2 δ(f-f0) + 1/2 δ(f-f0) *)
(* If denominators are quite close to zero, instead of having a strange sample with very high variance due to numerical instability when imposing zero derivative at t=0 , just generate another random sample *)
While[((d2udt2 =-2 ((((ω1)^2)) .( Re@spectrum1)))>0) ,
NormalSample =Partition[RandomReal[NormalDistribution[] (* normalized Gauss Distribution *),npuntos],2]; 
parteReal =NormalSample[[All,1]];
(* parteImaginaria = minus twice the sine amplitude at each frequency *)
parteImaginaria =NormalSample[[All,2]];
(* avoid DC component *)
(*SiRe =(Sqrt[ gain21]) parteReal / T;*)
SiRe =halfcosiaverage+halfcosistdev  parteReal ;
ΣSiRe=Total[SiRe ];  (* wind at t=0, its mean is ⟨Δv[t=0]⟩ = A *)
(*SiIm =(Sqrt[ gain21])parteImaginaria / T;*)
SiIm =halfsinistdev parteImaginaria;
(* ωiSiIm =  ω SiIm //N;*)
ΣωiSiIm = ω1 . SiIm; (* derivative of wind at t=0, it has zero mean *)
(*Σωi2SiIm = (((2 π))^2)Total[  ((freqs1^2) SiIm)] //N;*)

(* solución en media A -pero no tiene el valor exacto A en t=0-,spectrum1 =0.5(halfcosiaverage + parteReal halfcosistdev+I parteImaginaria halfsinistdev);*)
(* Factor 1/2 due to the fact that the Fourier transform is:
 sine:  1/2 I δ(f+f0) - 1/2 I δ(f-f0) 
or cosine: 1/2 δ(f-f0) + 1/2 δ(f-f0) *)
(* Spectrum constrained to Δv = A and dv/dt = 0 at t=0 *)
spectrum1 =0.5(SiRe+gain21((A-ΣSiRe)/(ΣS))+I (SiIm -ωS((ΣωiSiIm)/(Σω2S))));
];
(* Compute upper side spectrum for the signal to be real *)
reflectedspectrum =Join[{N[v0 ]},spectrum1, Conjugate@Reverse@Take[spectrum1,npuntos2-1]];
(* To avoid zero division at f=0, add DC gain 
gain2=Prepend[gain21,1/(2T)];
freqs=Prepend[freqs1,0.0];*)
(* Normalize taking into account Sqrt(N) factor in Mathematica FFT *)
(* Due to Mathematica default constants in Fourier Transforms, a factor (Sqrt[npuntos]) must be added to the reflected spectrum*)
timeseriesPeakGust = Re@InverseFourier[(Sqrt[npuntos]) reflectedspectrum];
SampleVariance = Variance@timeseriesPeakGust;
(* shift time samples for the peack to be in the middle of the sample *) 
timeseriesPeakGust =RotateLeft[ timeseriesPeakGust,npuntos2];
(*Since interpolation is periodic, wind at t = 0 is the same as time t = T. Otherwise, time series span is T-Δt *)
AppendTo[timeseriesPeakGust,timeseriesPeakGust[[1]]];
];(* End of ComputePeakGust *)

PlotPeakGust[]  := Module[{},
Column[
{Quiet@ListLinePlot[{ discreteACFwind+timeseriestdev,discreteACFwind,discreteACFwind-timeseriestdev,timeseriesPeakGust(*,{v0+A/E,v0+A/E}*)},(*AxesOrigin->{0,0},*)DataRange->{0,T},Epilog->{{Dashed,Line[{{(T-ExpectedNormalizedGustDuration)/2,v0+A/E},{(T+ExpectedNormalizedGustDuration)/2,v0+A/E}}]},{Dashed,Arrow[{{(T-ExpectedNormalizedGustDuration)/2,v0+A/E},{(T-ExpectedNormalizedGustDuration)/2,discreteACFwind[[1]]-timeseriestdev}}]},{Dashed, Arrow[{{(T+ExpectedNormalizedGustDuration)/2,v0+A/E},{(T+ExpectedNormalizedGustDuration)/2,discreteACFwind[[1]]-timeseriestdev}}]},Inset[Style[Text["Mean gust\nduration: "<>  ToString@Round@ExpectedNormalizedGustDuration<>" s"],Larger,Bold],{T/2,v0+A/E-0.6 timeseriestdev}],
Inset[Style[Text["⟨Gust⟩ + (Subscript[σ, gust])"],Larger,Bold],{0.07 T, discreteACFwind[[1]] + 1.1 timeseriestdev+0.1 A}],
Inset[Style[Text["⟨Gust⟩"],Larger,Bold],{0.07 T,discreteACFwind[[1]] - 0.1 timeseriestdev+0.1 A}],
Inset[Style[Text["⟨Gust⟩ - (Subscript[σ, gust])"],Larger,Bold],{0.07 T, discreteACFwind[[1]] - 0.9 timeseriestdev+0.1A}]
(*Inset[Style[ Text["⟨Gust⟩ + (Subscript[σ, gust])"], Larger, Bold], {(T - ExpectedNormalizedGustDuration)/2, 
   timeseriesaverage[[npuntos2-Round[ExpectedNormalizedGustDuration / Δt]]] + 2timeseriestdev},Automatic,Automatic,{ ExpectedNormalizedGustDuration/T,(Sqrt[2]) A(1-1/E)/(timeseriesaverage[[npuntos2]]-timeseriesaverage[[1]]+2timeseriestdev) }(*{Cos[70 Degree],Sin[70 Degree]}*)],Inset[Style[ Text["⟨Gust⟩ - (Subscript[σ, gust])"], Larger, Bold], {(T - ExpectedNormalizedGustDuration)/2, 
     timeseriesaverage[[npuntos2-Round[ExpectedNormalizedGustDuration / Δt]]] - 2 timeseriestdev},Automatic,Automatic,{ ExpectedNormalizedGustDuration/T,(Sqrt[2])A(1-1/E)/(timeseriesaverage[[npuntos2]]-timeseriesaverage[[1]]+2timeseriestdev) }(*{Cos[70 Degree],Sin[70 Degree]}*)]*)},
ImageSize->tamano,(*ImagePadding->{{25,25},{25,25}},*) Joined->joined,Filling->If[filled,{1->{2},3-> {2}},None], Frame-> True, FrameLabel-> {{"wind [m/s]",""},{"time t [s]",""}} ,GridLines->{Automatic,Automatic},GridLinesStyle->Directive[Yellow],PlotLabel->"Synthetized sample of equivalent wind peak gust",PlotStyle->{Dashed,DotDashed,Dashed,Thickness[0.04/  Log[npuntos]],Dashed} , InterpolationOrder-> If[npuntos>60,1,2],

PlotRange ->{ Min[discreteACFwind[[0]]-timeseriestdev,timeseriesPeakGust],Max[discreteACFwind[[npuntos2]]+timeseriestdev,timeseriesPeakGust]}]
(*PlotRange ->{ Min[timeseriesaverage-timeseriestdev,timeseriesaverage+timeseriestdev,timeseries],Max[timeseriesaverage-timeseriestdev,timeseriesaverage+timeseriestdev,timeseries]}]*)
,(*
{(* no dibujes componente DC pues sale fuera del rango en escala logarítmica *)
Quiet@ListLogLogPlot[ (* Use factors to transform discrete fourier complex values to PSD *) 
{(* halfcosiaverage ~ (Sqrt[((PSD(f))/T)]) ==> PSD(f) ~ T (halfcosiaverage^2); 
halfcosiaverage^2+quarterofcosivar ~ ((PSD(f))/T) ==> PSD(f) ~  T(halfcosiaverage^2+quarterofcosivar);*)Transpose@{freqs1, T(halfcosiaverage^2+quarterofcosivar)},Transpose@{freqs1,T(quarterofsinivar)},
(* Spectrum1 are complex coefs and have 1/2 factor compared to sin~cos coefs *)
(* Abs[spectrum1(f)] ~ (1/2) (1/(Sqrt[2]))(Sqrt[((PSD(f))/T)]) ==> PSD(f) ~ 2 T Abs[spectrum1(f)]; 
Re[spectrum1(f)] ~ (1/2) (Sqrt[((PSD(f))/T)]) ==> PSD(f) ~ 4 T ((Re[spectrum1(f)])^2);
Im[spectrum1(f)] ~ (1/2) (Sqrt[((PSD(f))/T)]) ==> PSD(f) ~ 4 T ((Im[spectrum1(f)])^2);*)
Transpose@{freqs1,4 T ((Re[spectrum1])^2)},Transpose@{freqs1,4 T ((Im[spectrum1])^2)}},
AxesOrigin->{1,0},
ImageSize->{800,300},Joined->joined,Frame-> True, FrameLabel-> {{"PSD of equivalent wind (m/s)^2/Hz",""},{"frequency f [Hz]",""}} ,GridLines->{Automatic,Automatic},GridLinesStyle->Directive[Yellow],PlotStyle->{Thickness[0.02],Thickness[0.02], Thin, Thin},PlotRange ->{ {0.01,100},{(10^(-12)),1000}},PlotLabel->"Power Spectral Density of equivalent wind (PSD)",ImageSize->tamano]
}*)

Quiet@ListLogLinearPlot[(* Use factors to transform discrete fourier complex values to PSD *) 
{(* halfcosiaverage ~ (Sqrt[((PSD(f))/T)]) ==> PSD(f) ~ T (halfcosiaverage^2); 
halfcosiaverage^2+quarterofcosivar ~ ((PSD(f))/T) ==> PSD(f) ~  T(halfcosiaverage^2+quarterofcosivar);*)Transpose@{freqs1,aux = T freqs1(halfcosiaverage^2+quarterofcosivar)},Transpose@{freqs1,T  freqs1(quarterofsinivar)},
(* Spectrum1 are complex coefs and have 1/2 factor compared to sin~cos coefs *)
(* Abs[spectrum1(f)] ~ (1/2) (1/(Sqrt[2]))(Sqrt[((PSD(f))/T)]) ==> PSD(f) ~ 2 T Abs[spectrum1(f)]; 
Re[spectrum1(f)] ~ (1/2) (Sqrt[((PSD(f))/T)]) ==> PSD(f) ~ 4 T ((Re[spectrum1(f)])^2);
Im[spectrum1(f)] ~ (1/2) (Sqrt[((PSD(f))/T)]) ==> PSD(f) ~ 4 T ((Im[spectrum1(f)])^2);*)
Transpose@{freqs1,4 T freqs1 ((Re[spectrum1])^2)},Transpose@{freqs1,4 T freqs1 ((Im[spectrum1])^2)}}, Frame-> True, 
ImageSize->{800,300},FrameLabel-> {{"f·PSD[f] = frequency power content (m/s)^2",""},{"frequency f [Hz]",""}} ,GridLines->{Automatic,Automatic},GridLinesStyle->Directive[Yellow],PlotStyle->Thick, PlotLabel-> "Frequency content of the sample \n(real and imaginary part of the power spectrum)", Filling-> Axis ,Epilog->{ Inset[Style[Text["(Subscript[(OverHat[s]), sample])=\!\(\*SqrtBox[\(sum\\ of\\ shadowed\\ area\)]\)=\n= "<>  ToString@NumberForm[(Sqrt[Variance@timeseriesPeakGust]),{3,2},NumberPadding->{"","0"}]<> " m/s"],Larger,Bold],{0.65(* centered horizontally at 65% of plot*),Center}]} , PlotRange-> {0,(*(Sqrt[SampleVariance])*) aux [[1]] }
 (* End of LogLinearPlot *)]}
, Alignment -> Center](*End of Column *)
]; (* End of PlotPeakGust *)

TabPeakGust[]:=Manipulate[ 
(* If wind model is numeric and updated, then show tab *)
If[UpdateModelIfNeeded[],
(* Compute and plot PeakGust sample *)
ComputePeakGust[];
PlotPeakGust[]
] (*End of if *),
{{A,1,"Gust peak Δv measured from interval mean  [m/s]"},0.25,5,0.25,Appearance->"Labeled"},
(* {{B,3,"Deceleration at maximum [m/(s^2)]"},1,20,1,Appearance->"Labeled"},*)
{{GeneratedSamplePeakGust,1,"Generated sample (generate new sample\nto update settings modified on other tabs):"},1,10,1,Appearance->"Labeled"},
(*{{joined,True,"joined"},{True,False}},
{{filled,True,"filling"},{True,False}},*)
SynchronousUpdating -> False (* Don't queue updates *),
SaveDefinitions-> True,TrackedSymbols:>{A,GeneratedSamplePeakGust,espectro,transformacion, v0,T,logΔt}, LocalizeVariables->False,Paneled-> False, (*Initialization->(Δt =1.;v0 = 10; T=60;logΔt=0;joined = True; filled=True;espectro[f_]=(((0.1^2)4 42/10 )/(((1+70.8(((f 42/10))^2)))^(5/6)));(*A=1;*) fmax = 50;) ,*)SynchronousInitialization->True]

TabΔvΔtGust

IniRampGust[] := Module[{},
];

ComputeRampGust[] := Module[{},
(* COEFFICIENTS COMPUTED FROM DISCRE SYSTEM  *)
sins =  2 Sin[π freqs1 dt] //N;
auxB = gain21 sins;
auxA = auxB  sins;
sdv2= Total[auxA]; (* sdv2 = Q *)
sinaverage = (Dv/sdv2)auxB;
quarterofcosivar=gain21;
halfcosistdev =(Sqrt[quarterofcosivar]);
(* quarterofsinivar =gain21 - gain21 auxA / sdv2; realmente en esta fórmula no tengo en cuenta los términos cruzados *)
Clear[McSin];
McSin = DiagonalMatrix[gain21]-Outer[Times,auxB/sdv2,auxB];
{eivals,eivecs}=Eigensystem[McSin];
(* the eigenvalues should be positive, but in practice tiny eigenvalues could eventually turn into  negative due to roundoff errors *)
NotNegative[x_]= Max[x,0];
Λ = (Sqrt[NotNegative/@eivals]);(*W=Transpose@vecs;SqrtMcSin =W.DiagonalMatrix[Λ ] ;*)
      (* Average shape of the gust. 1/2 factor to transform cosine and sine 1-sided Fourier coefs into 2-sided complex Fourier coefs *)
timeseriestdev=(Sqrt[0.5( Total@quarterofcosivar + Total@(eivals^2))]);
averagespectrum =0.5(Sqrt[npuntos]) Join[{N[2v0 ]},I sinaverage,Reverse@Take[-I sinaverage,npuntos2-1]];
(*varespectrum =√((npuntos/8 ) Join[{0},quarterofcosivar-I quarterofsinivar,Reverse@Take[quarterofcosivar+I quarterofsinivar,npuntos2-1]]);*)
(* stdevspectrum1 =0.5(Sqrt[npuntos]) Join[{0},halfcosistdev+I halfsinistdev,Reverse@Take[halfcosistdev-I halfsinistdev,npuntos2-1]];*)
timeseriesaverage =  RotateLeft[Re@InverseFourier[averagespectrum],npuntos2];
(*timeseriestdev=RotateLeft[Abs@InverseFourier[stdevspectrum],npuntos2];*)

SeedRandom[GeneratedSampleΔvΔtGust];
(* Generate samples in a way that changing the time step, the main waveform does not change: i.e., low part of the spectrum doesn't vary *)
(* normalized Gauss Distribution *)
NormalSample =Partition[RandomReal[NormalDistribution[] (* normalized Gauss Distribution *),npuntos],2]; 
parteReal =NormalSample[[All,1]];
(* parteImaginaria = minus twice the sine amplitude at each frequency *)
parteImaginaria =NormalSample[[All,2]];
(* avoid DC component *)
(*SiRe =(Sqrt[ gain21]) parteReal / T;*)
SiRe =halfcosistdev  parteReal ;
(* Impose Δv = 0 speed at t = 0 using constrained simulation*)
SiRe -=(gain21/(Total@gain21))Total@SiRe;
(* Adjust Real part to sum 0 in order that U(t=0) == 0 exactly, not only in average
posCount=negCount=0;
Do[If[SiRe[[i]]>0, posCount +=SiRe[[i]];, negCount +=SiRe[[i]];];,{i, npuntos2}];
posFactor = ((-2 negCount )/(posCount+negCount))//N; negFactor = ((2 posCount )/(posCount+negCount))//N ;Do[If[SiRe[[i]]>0,SiRe[[i]]*= posFactor ;,SiRe[[i]]*= negFactor;];,{i,npuntos2}];
Do not use this approach since it produces distorted gusts*)
SiIm =sinaverage+(Λ  parteImaginaria).eivecs;
(* Adjust Real part to sum 0 in order that U(t=0) == 0 *)

      (*SiIm = SiIm + sinaverage(1-actualΔv / Dv); actualΔv =SiIm. sins ;*)
(* xc == x+M.GT.Inverse[Q].(Y- 2 ΣSiResinπfnΔt(*G.x*)) *)
      SiIm +=((gain21 sins (Dv-SiIm. sins))/sdv2);
(*  SiIm +=(((SiIm^2) sins (Dv-SiIm. sins))/((SiIm^2). sins2 ));*)
(* Spectrum constrained *)
spectrum1 =0.5(SiRe+I SiIm);

(* Compute upper side spectrum for the signal to be real *)
reflectedspectrum =Join[{N[v0 ]},spectrum1, Conjugate@Reverse@Take[spectrum1,npuntos2-1]];
timeseriesStepGust = Re@InverseFourier[(Sqrt[npuntos]) reflectedspectrum];
(* shift time samples for the peak to be in the middle of the sample *) 
timeseriesStepGust =RotateLeft[ timeseriesStepGust,npuntos2];
(*Since interpolation is periodic, wind at t = 0 is the same as time t = T. Otherwise, time series span is T-Δt *)
AppendTo[timeseriesStepGust,timeseriesStepGust[[1]]];
AppendTo[timeseriesaverage,timeseriesaverage[[1]]];
];

PlotRampGust[] := Column[{Text@TraditionalForm[("* The standard deviation of Δv ≡ v(t)-v(t-Δt) is (Subscript[σ, Δv]) = "<> ToString@Sqrt@sdv2<>" m/s. Notice that (Subscript[σ, Δv]) is the typical difference of instantaneous value of equivalent wind between samples separated Δt in time. Therefore, the probability of exceding Δv = +"<>ToString@Dv<>" m/s upwards in time gap Δt = "<>ToString@dt<>" s is "<>ToString@CForm@ProbabilityWindExceedsΔvGauss[Dv,sdv2] <>" (per unit).\n* The meteorological datalogger and the annemometer act as low-pass filters. Notice that the difference of averaged or filtered measures are significantly smaller than the difference of instantaneous values, Δv ≡ v(t)-v(t-Δt).\n* Estimation based on linearity, normality, ⟨Wind⟩ = "<>ToString@v0<>" m/s on average and  (Subscript[σ, 1]) = (Subscript[σ, kaimal]) = "<> ToString@σkaimal[v0,Iref]<>" m/s as defined the standard deviation of longitudinal turbulence in IEC 61400-1 for normal turbulence.")],
Quiet@ListLinePlot[{ timeseriesaverage+timeseriestdev,timeseriesaverage,timeseriesaverage-timeseriestdev,timeseriesStepGust(*,{v0+A/E,v0+A/E}*)},(*AxesOrigin->{0,0},*)DataRange->{0,T},Epilog->{
Inset[Style[Text["⟨Gust⟩ + (Subscript[σ, gust])"],Larger,Bold],{0.07 T, timeseriesaverage[[1]] + 1.1 timeseriestdev}],
Inset[Style[Text["⟨Gust⟩"],Larger,Bold],{0.07 T,timeseriesaverage[[Max[npuntos/8,1]]] }],
Inset[Style[Text["⟨Gust⟩ - (Subscript[σ, gust])"],Larger,Bold],{0.07 T, timeseriesaverage[[Max[npuntos/8,1]]] - 1.1 timeseriestdev}]
},
ImageSize->tamano,(*ImagePadding->{{25,25},{25,25}},*) Joined->joined,Filling->If[filled,{1->{2},3-> {2}},None], Frame-> True, FrameLabel-> {{"wind [m/s]",""},{"time t [s]",""}} ,GridLines->{Automatic,Automatic},GridLinesStyle->Directive[Yellow],PlotLabel->"Synthetized sample of equivalent wind peak gust",PlotStyle->{Dashed,DotDashed,Dashed,Thickness[0.04/  Log[npuntos]],Dashed} , InterpolationOrder-> If[npuntos>60,1,2],
PlotRange ->{ Min[timeseriesaverage-timeseriestdev,timeseriesaverage+timeseriestdev,timeseriesStepGust],Max[timeseriesaverage-timeseriestdev,timeseriesaverage+timeseriestdev,timeseriesStepGust]}]
,
(* no dibujes componente DC pues sale fuera del rango en escala logarítmica *)
Quiet@ListLogLogPlot[ (* Use factors to transform discrete fourier complex values to PSD *) 
{(* halfcosiaverage ~ (Sqrt[((PSD(f))/T)]) ==> PSD(f) ~ T (halfcosiaverage^2); 
halfcosiaverage^2+quarterofcosivar ~ ((PSD(f))/T) ==> PSD(f) ~  T(halfcosiaverage^2+quarterofcosivar);*)
Transpose@{freqs1, T(quarterofcosivar)},Transpose@{freqs1,T(sinaverage^2 + quarterofsinivar)},
(* Abs[spectrum1(f)] ~ (1/2)(Sqrt[((PSD(f))/T)]) ==> PSD(f) ~ 4 T Abs[spectrum1(f)]; 
Re[spectrum1(f)] ~ (1/2)(Sqrt[((PSD(f))/(2 T))]) ==> PSD(f) ~ 8 T Re[spectrum1(f)];
Im[spectrum1(f)] ~ (1/2)(Sqrt[((PSD(f))/(2 T))]) ==> PSD(f) ~ 8 T Im[spectrum1(f)];*)
Transpose@{freqs1,8 T (((Re[spectrum1]))^2)},Transpose@{freqs1,8 T((( Im[spectrum1]))^2)}},
AxesOrigin->{1,0},
ImageSize->{800,300},Joined->joined,Frame-> True, FrameLabel-> {{"PSD of equivalent wind (m/s)^2/Hz",""},{"frequency f [Hz]",""}} ,GridLines->{Automatic,Automatic},GridLinesStyle->Directive[Yellow],PlotStyle->{Thickness[0.02],Thickness[0.02], Thin, Thin},PlotRange ->{ {0.01,100},{(10^(-12)),1000}},PlotLabel->"Power Spectral Density of equivalent wind (PSD)",ImageSize->tamano]
}, Alignment -> Center
](*End of PlotRampGust Column *)

TabΔvΔtGust[]:=Manipulate[ 
(* If wind model is numeric and updated, then show tab *)
If[UpdateModelIfNeeded[],
(* Compute and plot RampGust sample *)
ComputePeakGust[];
ComputeRampGust[] ;
PlotRampGust[] 
] (*End of if *),
{{Dv,1,"Wind difference Δv  [m/s]"},0.25,5,0.25,Appearance->"Labeled"},
{{dt,6,"Time gap Δt  [s]"},1,20,1,Appearance->"Labeled"},
{{GeneratedSampleΔvΔtGust,1,"Generated sample (generate new sample\nto update settings modified on other tabs):"},1,10,1,Appearance->"Labeled"},
SynchronousUpdating -> False (* Don't queue updates *),
SaveDefinitions-> True,TrackedSymbols:>{Dv,dt,GeneratedSampleΔvΔtGust,espectro,transformacion, v0,T,logΔt}, LocalizeVariables->False,Paneled-> False, SynchronousInitialization->True,Initialization->(Dv=1; dt=6;) ]

ListPlot[(auxA / sdv2)]

TabΔvΔt

(((Sin[α]))^2)//TrigReduce

TabΔvΔt[]:=
Manipulate[If[ButtonΔvΔt=="Refresh",
ButtonΔvΔt="";
ComputePeakGust[];
ComputeRampGust[] ;
];
If[NumericQ[σkaimal[v0,Iref] timeseriesPeakGust[[2]] espectro[freqs1[[2]] gain21[[1]] fblade fmax npuntos2 Δt/T]],σΔv2[Δt_] := 2Quiet@NIntegrate[ espectro[f] (((Sin[π f Δt]))^2),{f,1/3600,fblade,fmax}];
Column[{
Quiet@LogLogPlot[(Sqrt[σΔv2[dtime] ]),{dtime,0.1,600},Frame-> True, FrameLabel-> {{"standard (equivalent) wind variation (Subscript[σ, Δv]) [m/s]",""},{"time Δt [s]",""}} ,GridLines->{Automatic,Automatic},GridLinesStyle->Directive[Yellow],PlotStyle->Thickness[0.01],PlotLabel->"Typical wind variation (Subscript[σ, Δv]) in Δt",ImageSize->{800,500} (*Evaluate[EstiloLogPlot]*)
],Text@TraditionalForm["Notes:\n* (Subscript[σ, Δv])(Δt) is the standard deviation of Δv ≡ v(t)-v(t-Δt), the difference of instantaneous value of equivalent wind between samples separated Δt in time.\n* The meteorological datalogger and the annemometer act as low-pass filters. Notice that the difference of averaged or filtered measures are significantly smaller than the difference of instantaneous values, Δv ≡ v(t)-v(t-Δt).\n* Estimation based on linearity, normality, ⟨Wind⟩ = "<>ToString@v0<>" m/s on average and  (Subscript[σ, 1]) = (Subscript[σ, kaimal]) = "<> ToString@σkaimal[v0,Iref]<>" m/s as defined the standard deviation of longitudinal turbulence in IEC 61400-1 for normal turbulence.\n* The 10 minute standard deviation of the equivalent wind is (Subscript[σ, 10  min])= "<>ToString@N[(Sqrt[σdiezminutal2])]<>" m/s. For different turbulences, (Subscript[σ, Δv]) should be scalled proportionally.\n* For time spans up to a few minutes, wind diference Δv is approximately a normal random variable with zero mean and standard deivation (Subscript[σ, Δv]) (the probability of wind changes Δv is a straightforward calculation). Some authors based on experimental measures [Larsen, 2007] recomends the use of Gumbel EV1 distribution instead of Gaussian distribution for the estimation of rare gust events.\n* For time spans larger than a few minutes, PSD at the specific site must be used instead of typical spectra (wind forecast based on non-linear models should be used for long time horizons)."]}, Alignment -> Center
] (*End of column *),"Visit gust tab to initialize estimation of the number of gusts depending on gust peak"],{{ButtonΔvΔt,"","Refresh control of  typical wind variation Δv in Δt"},{"","Refresh"}},TrackedSymbols:>Manipulate, LocalizeVariables->False,Paneled-> False]

TabΔvΔt[]

Button, Print@GraphicsColumn

Manipulate[UpdateMe = False; RefreshMe = RefreshMe+1,{{UpdateMe,False,"Push to update"},{True,False} , ControlType -> Button},
SaveDefinitions-> True,TrackedSymbols:>{UpdateMe}, LocalizeVariables->False,Paneled-> False,ContinuousAction->False (*Setting ContinuousAction to False updates only when the control is released:*),
SynchronousUpdating -> False (* Don't queue updates *),
ControlPlacement->{Top,Top,Top,Top,Top,Top}, Alignment->Center,Initialization:>(RefreshMe = 1)]

Dynamic@Refresh[ UpdateMe = False; RefreshMe = RefreshMe + 1, TrackedSymbols :> {UpdateMe}]

Button["Update occurence probability estimations", Print@GraphicsColumn[{ RefreshMe = RefreshMe + 1}]]

TabGustsPerYear and PerDay

Tab TabGustsPerYear

TabGustsPerDay[]:=Manipulate[
If[ButtonGustsPerDay=="Refresh"
,ButtonGustsPerDay="";
ComputePeakGust[];
ComputeRampGust[];
];
If[NumericQ[σkaimal[v0,Iref] timeseriesPeakGust[[2]] espectro[freqs1[[2]] gain21[[1]] npuntos2 Δt/T]],Column[{
LogLogPlot[24*3600{PeakGustDensityA[dv],PeakGustDensityB[dv],PeakGustDensityC[dv]},{dv,1 (Sqrt[σinterval2]),3(Sqrt[σinterval2])},Frame-> True, FrameLabel-> {{"Expected number of gusts exceeding Δv in a day",""},{"Gust peak Δv [m/s] respect mean wind computed in "<>ToString@T<>" s intervals",""}} ,GridLines->{Automatic,Automatic},GridLinesStyle->Directive[Yellow],PlotStyle->Thickness[0.01],PlotLabel->("Number of gust exceeding Δv in a day assuming turbulence (Subscript[σ, kaimal]) = "<> ToString@σkaimal[v0,Iref]<>" m/s and time intervals T = "<>ToString@T<>" s"),ImageSize->{800,500} (*Evaluate[EstiloLogPlot]*)
],Text@TraditionalForm["Notes:\n* Very big peak gust are very scarce and its estimates can have significative error due to assumptions (normality, linearity, turbulence spectra, etc.).\n* Ten minutes wind average is assumed to be constant and equal to ⟨Wind⟩ = "<>ToString@v0<>" m/s and the standard deviation of longitudinal turbulence\nis (Subscript[σ, 1]) = (Subscript[σ, kaimal]) = "<> ToString@σkaimal[v0,Iref]<>" m/s as defined in IEC 61400-1 for normal turbulence.\n* Standard deviation of (equivalent) wind turbulence is (Subscript[σ, 10  min])= "<>ToString@N[(Sqrt[σdiezminutal2])]<>" m/s. For other wind averages, Δv should be scalled proportionally\nto turbulence ratios.\n* The mean gust duration is "<> ToString@Round@ExpectedNormalizedGustDuration<>"s, independently of the gust amplitude and turbulence intensity."]}, Alignment -> Center
] (*End of column *) ,"Visit gust tab to initialize estimation of the number of gusts depending on gust peak"],{{ButtonGustsPerDay,"","Refresh control of peak gust DAYLY occurence:"},{"","Refresh"}},TrackedSymbols:>Manipulate, LocalizeVariables->False,Paneled-> False]

TabGustsPerDay[]

TabGustsPerYear[]:=Manipulate[If[ButtonGustsPerYear=="Refresh",ButtonGustsPerYear=""];If[NumericQ[σkaimal[v0,Iref] timeseriesPeakGust[[2]] espectro[freqs1[[2]] gain21[[1]] npuntos2 Δt/T]],
Quiet@Column[{
LogLogPlot[365*24*3600{DensityWindExceedsΔvPerYearA[dv,vmin,vmax],DensityWindExceedsΔvPerYearC[dv,vmin,vmax],DensityWindExceedsΔvPerYearC[dv,vmin,vmax]} ,{dv,1 (Sqrt[σinterval2]),3(Sqrt[σinterval2])},Frame-> True, FrameLabel-> {{"Expected number of gusts exceeding Δv, averaged in a year",""},{"Gust peak Δv [m/s] respect mean wind computed in "<>ToString@T<>" s intervals",""}} ,GridLines->{Automatic,Automatic},GridLinesStyle->Directive[Yellow],PlotStyle->Thickness[0.01],PlotLabel->("Number of gust exceeding Δv in a year assuming reference turbulence \!\(\*SubscriptBox[\(σ\), \(kaimal, \\ 15\\ m/s\)]\) = "<> ToString@σkaimal[15,Iref]<>" m/s,\n Weibull( λ = "<> ToString[NumberForm[#,{4,1}]&@N@FactorForma ]<>", k = "<>ToString[NumberForm[#,{4,1}]&@N@FactorEscala ] <> "  m/s)  and T = "<>ToString@T<>" s averaging intervals."),ImageSize->{800,500} (*Evaluate[EstiloLogPlot]*)
],Text@TraditionalForm["Notes:\n* Very big peak gust are very scarce and its estimates can have significative error due to assumptions.\n* The standard deviation of longitudinal turbulence is (Subscript[σ, kaimal]) = "<> ToString@σkaimal[15,Iref]<>" m/s at ⟨Wind⟩ = 15 m/s as defined in IEC 61400-1 for normal turbulence.\n* Standard deviation of (equivalent) wind turbulence is (Subscript[σ, 1]) = (Subscript[σ, 10  min])= "<>ToString@N[(Sqrt[σdiezminutal2])]<>" m/s. For other wind averages, Δv should be scalled proportionally\nto turbulence ratios. Some authors based on experimental measures [Larsen, 2007] recomends the use of Gumbel EV1 distribution instead of Gaussian distribution for the estimation of rare gust events.\n* The mean gust duration is "<> ToString@Round@ExpectedNormalizedGustDuration<>"s, independently of the gust amplitude and turbulence intensity.\n* The yearly distribution of wind is assumed to correspond to a Weibull( λ = "<> ToString[NumberForm[#,{4,1}]&@N@FactorForma ]<>", k = "<>ToString[NumberForm[#,{4,1}]&@N@FactorEscala ] <> "  m/s)."]}, Alignment -> Center
] (*End of column *),"Visit gust tab to initialize estimation of the number of gusts depending on gust peak"],{{ButtonGustsPerYear,"","Refresh control of peak gust YEARLY occurence:"},{"","Refresh"}},TrackedSymbols:>Manipulate, LocalizeVariables->False,Paneled-> False]

TabGustsPerYear[]

Tabbed plots

Manipulate[
TabView[{
"Equivalent Wind Model" -> TabWindModel[] ,
"Wind Resource" ->   TabWindResource[2,7((2 )/(Sqrt[π])),10],
"Gust Statistics"-> TabGustStatistics[],
"Random Sample" ->   TabRandomSample[],
"Peak gust" ->   TabPeakGust[],(**)
"Ramp gust" ->   TabΔvΔtGust[] ,
"Typical wind variation Δv in Δt" ->  TabΔvΔt[],
"Gust per day" -> TabGustsPerDay[],
"Gust per year" ->
TabGustsPerYear[](**)
}(*End of pane list *)
,1 (*default tab*)
(*, Appearance -> "Limited"*)
](* End of tab *)
,
TrackedSymbols :> {}, Paneled-> False,
SynchronousUpdating -> False,
SynchronousInitialization->True (* This makes this part to be evaluated before the tabview inside *) ,
(*TimeConstraint->1,*)(*DynamicEvaluationTimeout-> 10,*)
 SaveDefinitions -> True (*, Deployed->False*),LocalizeVariables -> False , ContinuousAction -> False, 
Initialization:>  (
(* For avoiding function shadow, this input must be run before the interactive intialization *)
(* Depending how it is statred the interactive tab, some part is executed before *)
(* This is the reason why it is replicated three times in the code *)
(* In order toproduce Mathematica Player compatible files, Gets must be used intead of Needs *)
(*Needs["ColorHueGray2`",NotebookDirectory[]<>"ColorHueGray2.m"];*)
Get[NotebookDirectory[]<>"ColorHueGray2.m"];

minfreqinspectrogram = 0.1;

(* Fija estilo gráficas *)
NivelGris={0,1/2,1/3,5/6,2/3,1/6,1};
tamano = {500,400};
tamano = Medium;
axes = Automatic;
joined = True;
PuntosYrayas={{},Tiny,{0.02},{0.04,0.03},{0.08,0.02,0.01,0.02},{0.08,0.02},{0.01,0.02},{0.005,0.02,0.05,0.02},{0.005,0.02,0.02,0.02},{0.01},{0.06,0.02,0.01,0.02,0.01,0.02}};
EstiloGrafica ={
{GridLines-> Automatic,Frame->True, Axes-> False,AxesOrigin->axes,PlotStyle->Table[{Hue[(-y)/6], Thickness[0.01],Dashing[PuntosYrayas[[y+1]]]},{y,0,6}], Filling-> Axis,Background->None, ImageSize->tamano},
{GridLines-> Automatic,Frame->True, Axes-> False,AxesOrigin->axes,PlotStyle->Table[{Hue[(-y)/6], Thickness[0.006*If[y>0,1.5,2.5]],Dashing[PuntosYrayas[[y+1]]]},{y,0,6}], Filling-> Axis,Background->None, ImageSize->tamano},
{GridLines-> Automatic,Frame->True, Axes-> False,AxesOrigin->axes,PlotStyle->Table[{Hue[(-y)/6], Thickness[0.006*If[y>0,1.5,2.5]],Dashing[{}]},{y,0,6}], Filling-> Axis,Background->None, ImageSize->tamano},
{GridLines-> None,Frame->True, Axes-> False,AxesOrigin->axes,PlotStyle->Table[{Hue[-y/10], Thickness[0.01],Dashing[PuntosYrayas[[y+1]]]},{y,0,10}], Filling-> Axis,Background->None, ImageSize->tamano},
{GridLines-> Automatic,Frame->True, Axes-> False,AxesOrigin->axes,PlotStyle->Table[{GrayLevel[NivelGris[[y+1]]], Thickness[0.01],Dashing[PuntosYrayas[[y+1]]]},{y,0,6}], Filling-> Bottom,Background->None, ImageSize->tamano},
{GridLines-> Automatic,Frame->True, Axes-> False,AxesOrigin->axes,PlotStyle->Table[{GrayLevel[NivelGris[[y+1]]], Thickness[0.006*If[y>0,1.5,2.5]],Dashing[PuntosYrayas[[y+1]]]},{y,0,6}], Filling-> Bottom,Background->None, ImageSize->tamano}
};
EstiloLogPlot ={Frame->True, Axes-> False, AxesOrigin->Automatic, PlotStyle->Table[{GrayLevel[NivelGris[[y+1]]], Thickness[0.006*If[y>0,1.5,2.5]],Dashing[PuntosYrayas[[y+1]]]},{y,0,6}], Filling-> None, Background->None, GridLines->{Automatic,Automatic}, GridLinesStyle->Directive[Yellow], ImageSize->tamano};
EstiloListLogPlot = Join[{Joined->joined}, EstiloLogPlot];

(SetOptions[Plot, #1] &) /@ EstiloGrafica[[1]];
(SetOptions[ListPlot, #1] &) /@ EstiloGrafica[[1]];
(SetOptions[ListLinePlot, #1] &) /@ EstiloGrafica[[1]];
(SetOptions[DateListPlot, #1] &) /@ EstiloGrafica[[1]];
(SetOptions[LogPlot, #1] &) /@  EstiloLogPlot;
(SetOptions[ListLogPlot, #1] &) /@ EstiloListLogPlot;
(SetOptions[LogLogPlot, #1] &) /@ EstiloLogPlot;
(SetOptions[ListLogLogPlot, #1] &) /@ EstiloListLogPlot;
(SetOptions[ListLogLinearPlot, #1] &) /@ EstiloListLogPlot;
(SetOptions[LogLinearPlot, #1] &) /@ EstiloLogPlot;
(SetOptions[ListLogLogPlot, #1] &) /@ EstiloListLogPlot;

SetOptions[Plot3D, ColorFunction->(ColorHueGray[0.9` -0.6666666666666666`#3,Min[1,0.25` +1.15` Abs[-0.4`+Cos[π (-0.15`+#3)]]],Min[1,0.45` +0.6` #3]]&), MeshFunctions->{#1&,#2&,#3&}, Mesh->{5,5,9},MeshStyle->{LightGray,LightGray,Darker@Brown}];
SetOptions[ContourPlot, ColorFunction->(RGBColor@ColorHueGray[0.2333333333333334` +0.6666666666666666` #1,Min[1,0.25` +1.15` Abs[-0.4`+Cos[2.670353755551324` -π #1]]],Min[1,1.05` -0.6` #1]]&),ContourLabels -> All ];
SetOptions[DensityPlot, ColorFunction->(RGBColor@ColorHueGray[0.2333333333333334` +0.6666666666666666` #1,Min[1,0.25` +1.15` Abs[-0.4`+Cos[2.670353755551324` -π #1]]],Min[1,1.05` -0.6` #1]]&), MeshFunctions->{#1&,#2&,#3&}, Mesh->{5,5,9}, MeshStyle->{LightGray,LightGray,Darker@Brown}];
SetOptions[ListDensityPlot, ColorFunction->(RGBColor@ColorHueGray[0.2333333333333334` +0.6666666666666666` #1,Min[1,0.25` +1.15` Abs[-0.4`+Cos[2.670353755551324` -π #1]]],Min[1,1.05` -0.6` #1]]&), MeshFunctions->{#1&,#2&,#3&}, Mesh->{5,5,9}, MeshStyle->{LightGray,LightGray,Darker@Brown}];

(SetOptions[ListLogLogPlot, #1] &) /@ {Frame-> True, Axes-> False,AxesOrigin->{1,0},Joined->joined, FrameLabel-> {{"PSD of equivalent wind (m/s)^2/Hz",""},{"frequency\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"f\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"[\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"Hz\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"]\",\nFontSlant->\"Italic\"]\)",""}} ,GridLines->{Automatic,Automatic},GridLinesStyle->Directive[Yellow],PlotStyle->{Thick, Thin},Filling-> None, ImageSize->{800,400}};
(SetOptions[ListLogLinearPlot, #1] &) /@ {Frame-> True, Axes-> False,AxesOrigin->{1,0},Joined->joined,FrameLabel-> {{"PSD of equivalent wind (m/s)^2/Hz",""},{"frequency\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"f\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"[\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"Hz\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"]\",\nFontSlant->\"Italic\"]\)",""}} ,GridLines->{Automatic,Automatic},GridLinesStyle->Directive[Yellow],PlotStyle->{Thick, Thin},Filling-> None, ImageSize->{800,400}};
(SetOptions[LogLinearPlot, #1] &) /@ {Frame-> True, Axes-> False,AxesOrigin->{1,0},
FrameLabel-> {{"\!\(\*
StyleBox[\"f\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"·\",\nFontSlant->\"Italic\"]\)PSD[f] = relative variance (m/s)^2 at freq. \!\(\*
StyleBox[\"f\",\nFontSlant->\"Italic\"]\)",""},{"frequency\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"f\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"[\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"Hz\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"]\",\nFontSlant->\"Italic\"]\)",""}} ,GridLines->{Automatic,Automatic},GridLinesStyle->Directive[Yellow],PlotStyle->{Thick, Thin}, Filling-> Axis, ImageSize->{800,400} ,PlotLabel-> "The area beneath the \!\(\*
StyleBox[\"f\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"·\",\nFontSlant->\"Italic\"]\)PSD(\!\(\*
StyleBox[\"f\",\nFontSlant->\"Italic\"]\)) plot in Log-Lin axis is the variance of the  
equivalent wind and the height is the relative contribution at each frequency f"
};
(SetOptions[ListLogLogPlot, #1] &) /@ {Frame-> True, Axes-> False,AxesOrigin->{1,0},Joined->joined,
FrameLabel-> {{"\!\(\*
StyleBox[\"f\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"·\",\nFontSlant->\"Italic\"]\)PSD[f] = relative variance (m/s)^2 at freq. \!\(\*
StyleBox[\"f\",\nFontSlant->\"Italic\"]\)",""},{"frequency\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"f\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"[\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"Hz\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"]\",\nFontSlant->\"Italic\"]\)",""}} ,GridLines->{Automatic,Automatic},GridLinesStyle->Directive[Yellow],PlotStyle->{Thick, Thin},  Filling-> Axis, ImageSize->{800,400},PlotLabel-> "The area beneath the \!\(\*
StyleBox[\"f\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"·\",\nFontSlant->\"Italic\"]\)PSD(\!\(\*
StyleBox[\"f\",\nFontSlant->\"Italic\"]\)) plot in Log-Lin axis is the variance of the  
equivalent wind and the height is the relative contribution at each frequency f"
};
(**************************************************************************************************)

AbortProtect[IniWindModel[]];
AbortProtect[IniWindResource[2 (*FactorForma0*),7((2 )/(Sqrt[π]))(*FactorEscala0*), 10 (*v0 windspeed *)]];
AbortProtect[IniGustStatistics[]];
AbortProtect[IniRandomSample[]];
AbortProtect[IniPeakGust[]] ;
AbortProtect[ComputeWindModel[]];
AbortProtect[DiscretizeSpectrum[]];
AbortProtect[MeanPeakGust[(*npuntos2,T*)]] ;
AbortProtect[ComputeGustStatistics[]];
transformacion0=transformacion;T0 = T;logΔt0 =logΔt; 
espectro0[f_]= espectro[f];
A0=A; 
v00 = v0;
)
]