classdef Wind < Config
% Wind: Represents Skyward's custom wind model.
%
%   Constructor:
%       - Wind: Creates an instance of the Wind class.
%           Loaded config: windConfig.m > windCustom
%           Loaded data: -
%           Arguments:
%               - mission: Mission, mission object
%               - varIn: (optional) config source. Alternative to config.m
%               file

    properties(Dependent)
        altitudes               (1, :)

        magnitudeDistribution   (1, :) ... 
            {mustBeMember(magnitudeDistribution, {'u', 'g'})}   % [-] Magnitude distrubution: uniform (default), gaussian
        azimuthDistribution     (1, :) ...
            {mustBeMember(azimuthDistribution, {'u', 'g'})}     % [-] Azimuth   distrubution: uniform (default), gaussian
        elevationDistribution   (1, :) ... 
            {mustBeMember(elevationDistribution, {'u', 'g'})}   % [-] Elevation distrubution: uniform (default), gaussian
        
        magnitudeParameters     (2, :)                          % [m/s] Distribution parameters: [min; max], [mu; sigma]
        azimuthParameters       (2, :)                          % [deg] Distribution parameters: [min; max], [mu; sigma]
        elevationParameters     (2, :)                          % [deg] Distribution parameters: [min; max], [mu; sigma]
    end
    
    properties(Access = private)
        isReady                 (1, 1)  logical                   % true if all properties are set
        dataInerpolant          (1, 1)  griddedInterpolant      % [-] Interpolates [magnitude, azimuth, elevation]

        ALTITUDES
        
        MAGNITUDE_DISTRIBUTION  = ''
        AZIMUTH_DISTRIBUTION    = ''
        ELEVATION_DISTRIBUTION  = ''

        MAGNITUDE_PARAMETERS    = []
        AZIMUTH_PARAMETERS      = []
        ELEVATION_PARAMETERS    = []
    end

    properties(SetAccess = private)
        azimuth
        elevation
        magnitude
    end

    properties(Constant, Access = protected)
        configName = 'windConfig.m'
        variableName = 'windCustom'
    end

    methods % Getters
        function value = get.altitudes(obj), value = obj.ALTITUDES; end

        function value = get.magnitudeDistribution(obj), value = obj.MAGNITUDE_DISTRIBUTION; end
        function value = get.azimuthDistribution(obj), value = obj.AZIMUTH_DISTRIBUTION; end
        function value = get.elevationDistribution(obj), value = obj.ELEVATION_DISTRIBUTION; end
            
        function value = get.magnitudeParameters(obj), value = obj.MAGNITUDE_PARAMETERS; end
        function value = get.azimuthParameters(obj), value = obj.AZIMUTH_PARAMETERS; end
        function value = get.elevationParameters(obj), value = obj.ELEVATION_PARAMETERS; end
    end

    methods % Setters
        function obj = set.altitudes(obj, value)
            obj.ALTITUDES = value;
            obj = obj.updateAll();
        end

        function obj = set.magnitudeDistribution(obj, value)
            obj.MAGNITUDE_DISTRIBUTION = value;
            obj.isReady = obj.checkProperties();
            obj.magnitude = obj.updateDistribution(value, obj.MAGNITUDE_PARAMETERS);
            obj = obj.updateInterpolant();
        end

        function obj = set.magnitudeParameters(obj, value)
            obj.MAGNITUDE_PARAMETERS = value;
            obj.isReady = obj.checkProperties();
            obj.magnitude = obj.updateDistribution(obj.MAGNITUDE_DISTRIBUTION, value);
            obj = obj.updateInterpolant();
        end

        function obj = set.azimuthDistribution(obj, value)
            obj.AZIMUTH_DISTRIBUTION = value;
            obj.isReady = obj.checkProperties();
            obj.azimuth = obj.updateDistribution(value, obj.AZIMUTH_PARAMETERS);
            obj = obj.updateInterpolant();
        end

        function obj = set.azimuthParameters(obj, value)
            obj.AZIMUTH_PARAMETERS = value;
            obj.isReady = obj.checkProperties();
            obj.azimuth = obj.updateDistribution(obj.AZIMUTH_DISTRIBUTION, value);
            obj = obj.updateInterpolant();
        end

        function obj = set.elevationDistribution(obj, value)
            obj.ELEVATION_DISTRIBUTION = value;
            obj.isReady = obj.checkProperties();
            obj.elevation = obj.updateDistribution(value, obj.ELEVATION_PARAMETERS);
            obj = obj.updateInterpolant();
        end

        function obj = set.elevationParameters(obj, value)
            obj.ELEVATION_PARAMETERS = value;
            obj.isReady = obj.checkProperties();
            obj.elevation = obj.updateDistribution(obj.ELEVATION_DISTRIBUTION, value);
            obj = obj.updateInterpolant();
        end
    end

    methods
        function obj = Wind(mission, varIn)
            arguments(Input)
                mission Mission = Mission()
                varIn = []
            end
            obj@Config(mission, varIn);
        end

        function [uw, vw, ww] = getVels(obj, z)
            if ~obj.isReady
                error(['Parameters and distributions must be the same ' ...
                    'size as altitudes or scalar']);
            end

            if isscalar(obj.altitudes)
                uw = round( - obj.magnitude * cos(obj.azimuth) * cos(obj.elevation), 6);
                vw = round( - obj.magnitude * sin(obj.azimuth) * cos(obj.elevation), 6);
                ww = round( - obj.magnitude * (-sin(obj.elevation)), 6);
            else
                data = obj.dataInerpolant(z);
                mag = data(1); az = data(2); el = data(3);

                uw = round( - mag * cos(az) * cos(el), 6);
                vw = round( - mag * sin(az) * cos(el), 6);
                ww = round( - mag * (-sin(el)), 6);
            end
        end

        function obj = updateAll(obj)
            obj.magnitude = obj.updateDistribution(obj.MAGNITUDE_DISTRIBUTION, obj.MAGNITUDE_PARAMETERS);
            obj.azimuth   = obj.updateDistribution(obj.AZIMUTH_DISTRIBUTION, obj.AZIMUTH_PARAMETERS);
            obj.elevation = obj.updateDistribution(obj.ELEVATION_DISTRIBUTION, obj.ELEVATION_PARAMETERS);

            obj.isReady = obj.checkProperties();
            obj = obj.updateInterpolant();
        end

        function obj = updateInterpolant(obj)
            % UPDATEINTERPOLANT: Interpolates [magnitude; azimuth; elevation] based on altitude
            %
            %   Returns:
            %       - interpolant: griddedInterpolant for [magnitude; azimuth; elevation]

            if ~obj.isReady, return; end
            if isscalar(obj.ALTITUDES), return; end

            % Prepare data for interpolation
            data = [obj.magnitude; obj.azimuth; obj.elevation];
            gridVec = obj.ALTITUDES;

            % Create gridded interpolant
            obj.dataInerpolant = griddedInterpolant(gridVec, data', 'linear', 'nearest');
        end
    end

    methods(Access = private)
        function ready = checkProperties(obj)
            % Check if STATIC, DYNAMIC, GEOMETRY, and STATE are non-empty

            ready = ...
                length(obj.magnitude) == length(obj.ALTITUDES) && ...
                length(obj.azimuth) == length(obj.ALTITUDES) && ...     
                length(obj.elevation) == length(obj.ALTITUDES);
        end

        function vec = updateDistribution(obj, dist, parameters)
            s = length(obj.ALTITUDES);
            vec = [];
 
            dist = string(dist);
            lDist = length(dist);
            lPara = width(parameters);

            if lDist == 0, dist = "u"; lDist = 1; end
            if lPara == 0, parameters = [0; 0]; lPara = 1; end

            if lDist == 1, dist = repmat(dist, 1, s); lDist = s; end
            if lPara == 1, parameters = repmat(parameters, 1, s); lPara = s; end

            if lDist ~= s, return; end
            if lPara ~= s, return; end
            
            isUniform = strcmp(dist, "u");

            uniformDist = @(val1, val2) rand(1)*(val2 - val1) + val1;
            gaussianDist = @(mean, sigma) normrnd(mean, sigma, 1);

            vec = nan(1,s);
            for i = 1:s
                if isUniform(i)
                    vec(i) = uniformDist(parameters(1,i), parameters(2,i));
                else
                    vec(i) = gaussianDist(parameters(1,i), parameters(2,i));
                end
            end
        end
    end
end