within AIDAModelica;

model Wind2DModel
  
  parameter Real rho = 1.225
  "Air density [kg/m3]. Used in the aerodynamic drag formulation. Typical value at sea level under ISA conditions.";

  parameter Real Cd = 1.0
    "Equivalent aerodynamic drag coefficient of the vehicle projected in the horizontal plane. Represents the global aerodynamic resistance to wind.";
  
  parameter Real A = 0.05
    "Equivalent lateral reference area exposed to the wind [m2]. Represents the projected surface area of the vehicle that contributes to aerodynamic drag.";
  
  parameter Real eps = 1e-6
    "Numerical regularization term added inside the wind speed computation to avoid singularities when the wind magnitude approaches zero.";
  
  parameter Real Wmean = 8
    "Mean horizontal wind speed magnitude [m/s]. Defines the steady-state wind component in the environment.";
  
  parameter Real windDir = 0
    "Direction of the mean wind in the ground frame [rad]. Angle measured from the X axis toward the Y axis.";
  
  parameter Real Tg = 3
    "Duration of a gust event [s]. Controls the temporal length of the cosine-shaped gust profile.";
  
  parameter Real gustRate = 0.02
    "Average occurrence rate of gust events [1/s]. The expected number of gusts per second. Used to generate stochastic gust arrival times.";
  
  parameter Real WgustMax = 6
    "Maximum amplitude of gust velocity perturbations [m/s]. Each gust amplitude is randomly generated between 0 and this value.";
  
  parameter Real At1 = 0.6
    "Amplitude of the first sinusoidal turbulence component acting along the X direction [m/s].";
  
  parameter Real At2 = 0.25
    "Amplitude of the second sinusoidal turbulence component acting along the X direction [m/s].";
  
  parameter Real Bt1 = 0.6
    "Amplitude of the first sinusoidal turbulence component acting along the Y direction [m/s].";
  
  parameter Real Bt2 = 0.25
    "Amplitude of the second sinusoidal turbulence component acting along the Y direction [m/s].";
  
  parameter Real f1 = 0.15
    "Frequency of the first turbulence harmonic [Hz]. Represents a low-frequency atmospheric fluctuation.";
  
  parameter Real f2 = 0.7
    "Frequency of the second turbulence harmonic [Hz]. Represents a higher-frequency turbulence component.";
  
  parameter Real f3 = 0.21
    "Frequency of the first turbulence harmonic along the Y direction [Hz].";
  
  parameter Real f4 = 0.9
    "Frequency of the second turbulence harmonic along the Y direction [Hz].";
  
  parameter Real phi1 = 0
    "Phase shift of the first X-direction turbulence component [rad]. Allows decorrelation of turbulence signals.";
  
  parameter Real phi2 = 0.8
    "Phase shift of the second X-direction turbulence component [rad].";
  
  parameter Real phi3 = 1.3
    "Phase shift of the first Y-direction turbulence component [rad].";
  
  parameter Real phi4 = 2.1
    "Phase shift of the second Y-direction turbulence component [rad].";

  parameter Integer seed0 = 12345 "Initial random seed";
  
  Modelica.Blocks.Interfaces.RealOutput Fy
    "Aerodynamic Y wind force in ground frame [N]" annotation(
    Placement(transformation(origin = {110, -76}, extent = {{-10, -10}, {10, 10}}), iconTransformation(origin = {70, -110}, extent = {{-10, -10}, {10, 10}}, rotation = -90)));

  Modelica.Blocks.Interfaces.RealOutput Fx
    "Aerodynamic X wind force in ground frame [N]" annotation(
    Placement(transformation(origin = {110, 76}, extent = {{-10, -10}, {10, 10}}), iconTransformation(origin = {-70, -110}, extent = {{-10, -10}, {10, 10}}, rotation = -90)));

  Modelica.Blocks.Interfaces.RealOutput windVelX "Wind velocity component in X direction [m/s]";

  Modelica.Blocks.Interfaces.RealOutput windVelY "Wind velocity component in Y direction [m/s]";
  
protected

  Real wx_mean;
  Real wy_mean;
  Real wx_turb;
  Real wy_turb;
  Real wx_gust;
  Real wy_gust;
  Real wx;
  Real wy;
  Real windSpeed;

  Real tau;
  Real gust;

  discrete Real gustStart(start = -100);
  discrete Real nextGustTime(start = 5);
  discrete Real gustAmp(start = 0);
  discrete Real gustDirection(start = 0);

  discrete Integer seed(start = seed0);

  constant Integer a = 1664525;
  constant Integer c = 1013904223;
  constant Integer m = 2147483647;
  
algorithm
  when time >= pre(nextGustTime) then
  
    seed := mod(a*pre(seed) + c, m);
    gustStart := time;
    gustAmp := WgustMax * seed / m;

    seed := mod(a*seed + c, m);
    gustDirection := 2*Modelica.Constants.pi * seed / m;

    seed := mod(a*seed + c, m);
    nextGustTime := time - log(max(seed/m,1e-6)) / gustRate;
  end when;
equation

  wx_mean = Wmean * cos(windDir);
  wy_mean = Wmean * sin(windDir);

  wx_turb =
       At1 * sin(2*Modelica.Constants.pi*f1*time + phi1)
     + At2 * sin(2*Modelica.Constants.pi*f2*time + phi2);

  wy_turb =
       Bt1 * sin(2*Modelica.Constants.pi*f3*time + phi3)
     + Bt2 * sin(2*Modelica.Constants.pi*f4*time + phi4);

  

  tau = (time - gustStart) / Tg;

  gust =
       noEvent(if tau >= 0 and tau <= 1 then
        0.5*gustAmp*(1 - cos(2*Modelica.Constants.pi*tau))
      else
        0);

  wx_gust = gust*cos(gustDirection);
  wy_gust = gust*sin(gustDirection);

  wx = wx_mean + wx_gust + wx_turb;
  wy = wy_mean + wy_gust + wy_turb;

  windVelX = wx;
  windVelY = wy;

  windSpeed = sqrt(wx*wx + wy*wy + eps);

  Fx = 0.5*rho*Cd*A*windSpeed*wx;
  Fy = 0.5*rho*Cd*A*windSpeed*wy;
  
annotation(
    Diagram,
    Icon(graphics = {Rectangle(extent = {{-100, 100}, {100, -100}})}));
end Wind2DModel;