% VAD Analysis Script for Doppler Lidar Data with Interactive Altitude Limit
% and Wind Direction as the "Blowing To" Direction.
%
% Assumptions:
%   - lidar.RWS: 60 x 2848 radial wind speed measurements.
%   - lidar.azimuth: 60 x 2848 azimuth angles in degrees.
%   - lidar.elevation: 60 x 2848 elevation angles in degrees.
%   - lidar.ranges: 1 x 2848 range values (m), representing altitude.
%
% For each range gate, the model is:
%   V_r = u*cos(elev)*cos(az) + v*cos(elev)*sin(az) + w*sin(elev)
% where [u; v; w] are the wind components.
% Horizontal wind speed: U = sqrt(u^2 + v^2)
% Wind direction (direction the wind is blowing):
%   wind_dir = mod(270 - atan2(v,u)*180/pi + 180, 360)

% Determine number of beams and range gates
[n_beams, n_gates] = size(lidar.RWS);

% Preallocate result arrays (1 x n_gates)
U   = zeros(1, n_gates);   % horizontal wind speed [m/s]
W   = zeros(1, n_gates);   % vertical wind speed [m/s]
Dir = zeros(1, n_gates);   % wind direction [degrees] (wind is blowing toward)

% Loop over each range gate to compute wind components via least squares
for j = 1:n_gates
    % Extract data for the j-th range gate
    V_r = lidar.RWS(:, j);          % Radial wind speeds (60 x 1)
    az  = lidar.azimuth(:, j);       % Azimuth angles (degrees)
    el  = lidar.elevation(:, j);     % Elevation angles (degrees)
    
    % Build the design matrix A for the linear model:
    % [cosd(el)*cosd(az), cosd(el)*sind(az), sind(el)]
    A = [cosd(el).*cosd(az), cosd(el).*sind(az), sind(el)];
    
    % Solve for [u; v; w] using least squares
    x = A \ V_r;
    u = x(1);
    v = x(2);
    w = x(3);
    
    % Compute horizontal wind speed
    U(j) = sqrt(u^2 + v^2);
    
    % Compute wind direction representing the direction the wind is blowing toward.
    % Add 180 degrees to the traditional meteorological "from" direction.
    Dir(j) = mod(270 - atan2(v,u)*180/pi + 180, 360);
    
    % Store vertical wind speed
    W(j) = w;
end

% Plot the vertical profiles as a function of altitude (range)
figure;

subplot(1,3,1);
plot(U, lidar.ranges, 'LineWidth', 1.5);
xlabel('Horizontal Wind Speed (m/s)');
ylabel('Altitude (m)');
set(gca, 'YDir','normal');  % Altitude increases upward
title('Horizontal Wind Speed Profile');

subplot(1,3,2);
plot(W, lidar.ranges, 'LineWidth', 1.5);
xlabel('Vertical Wind Speed (m/s)');
ylabel('Altitude (m)');
set(gca, 'YDir','normal');
title('Vertical Wind Speed Profile');

subplot(1,3,3);
plot(Dir, lidar.ranges, 'LineWidth', 1.5);
xlabel('Wind Direction (deg) [Blowing To]');
ylabel('Altitude (m)');
set(gca, 'YDir','normal');
title('Wind Direction Profile');

% Prompt the user for the maximum altitude to display
max_alt = input('Enter maximum altitude for display (in m): ');

% Update the y-axis limits for each subplot to use the user-specified maximum altitude.
subplot(1,3,1);
ylim([min(lidar.ranges) max_alt]);

subplot(1,3,2);
ylim([min(lidar.ranges) max_alt]);

subplot(1,3,3);
ylim([min(lidar.ranges) max_alt]);

% Optionally, store the computed profiles in the lidar structure
lidar.U = U;
lidar.W = W;
lidar.Dir = Dir;
