function [T,S,P,z_ac] = RotorSail_Calc_wrapper(AWA,AWS,SR,D,De,height,rho,mu)
% Wrapper function for the Rotor-sail calculation of the thrust force (T), 
% Side force (S), Power required to spin the rotor (P) and the vertical 
% distance from the deck to the aerodynamic centre of lift (z_ac).
% Currently a Surrogate model is used based on the data from 161-Swanson, 2014-Marco,
% 2016-Marco, 2022-Kwon, 2023-Chen and 2023-Schot is used, with additional fixes for SR=0.

De_D = De/D;        % Endplate diameter ratio [-]
Re = (AWS*D)/mu;    % Reynolds number [-]
AR = height/D;      % Rotor Aspect ratio [-]
S = D * height;     % Rotor frontal area [m^2]

% Initialize arrays
CL = zeros(length(AWS),length(SR));
CD = zeros(length(AWS),length(SR));
Cf = zeros(length(AWS),length(SR));
Lift = zeros(length(AWS),length(SR));
Drag = zeros(length(AWS),length(SR));
P = zeros(length(AWS),length(SR));

tic
disp("Beginning calculation of specified Rotor-sail.");
for j = 1:length(AWS)
    disp("Calculating Aerodynamic Forces for AWS = " + AWS(j) + " [m/s]");
    for k = 1:length(SR)
        [CL(j,k),CD(j,k),Cf(j,k)] = RotorSail_SurrogateModel(AR,De_D,Re(j),SR(k));
        if SR(k) == 0
            CL(j,k) = 0;    % Fixing zero-lift for SR=0
        end
        Lift(j,k) = 0.5 * rho * AWS(j)^2 * S * CL(j,k);
        Drag(j,k) = 0.5 * rho * AWS(j)^2 * S * CD(j,k);
        P(j,k) = 0.5 * rho * (SR(k) * AWS(j))^3 * S * pi * Cf(j,k);
    end
end

% Calculating the Thrust and Side forces:
T = zeros(length(AWA),length(AWS),length(SR));
S = zeros(length(AWA),length(AWS),length(SR));
for i = 1:length(AWA)
    for j = 1:length(AWS)
        for k = 1:length(SR)
            T(i,j,k) = Lift(j,k) * sind(AWA(i)) - Drag(j,k) * cosd(AWA(i));
            S(i,j,k) = -Lift(j,k) * cosd(AWA(i)) - Drag(j,k) * sind(AWA(i));
        end
    end
end
z_ac = height/2 * ones(length(AWS),length(SR));    % For now
toc
disp("Calculations of " + length(AWA)*length(AWS)*length(SR) + ...
    " calculation points complete!");
end