% Define constants
V_ship_knots = 14; % Desired ship speed in knots
V_ship_mps = V_ship_knots * 0.514444; % Convert ship speed to m/s
max_diameter = 4; % Maximum allowable diameter in meters

% Define ranges for J and specific values for P/D
J_values = linspace(0, 1.3, 200); % Denser range for advance coefficients
P_D_values = [0.6, 0.8, 1.0, 1.2, 1.4]; % Extended P/D ratios to match reference plot
K_Q_scaling_factor = 10;

best_efficiency = 0;
optimal_P_D = 0;
optimal_J = 0;
D_optimal = 0;
P_optimal = 0;

% Coefficients from Table 6.17 for Nozzle No. 194 (example values)
coefficients = [
    0,0,0.03055,0.006735;0,1,-0.148687,0;
    0,2,0,-0.016306;0,3,-0.391137,0;
    0,4,0,-0.007244;0,5,0,0;
    0,6,0,0;1,0,0,0;
    1,1,-0.432612,0;1,2,0,-0.024012;
    1,3,0,0;
    1,4,0,0;1,5,0,0;
    1,6,0,0;2,0,0.667657,0;
    2,1,0,0;2,2,0.285076,0.005193;
    2,3,0,0;2,4,0,0;
    2,5,0,0;2,6,0,0;
    3,0,-0.172529,0.046605;3,1,0,0;
    3,2,0,0;
    3,3,0,0;3,4,0,0;
    3,5,0,0;3,6,0,0;
    4,0,0,-0.007366;4,1,0,0;
    4,2,0,0;4,3,0,0;
    4,4,0,0;
    4,5,0,0;4,6,0,0;
    5,0,0,0;5,1,0,0;
    5,2,0,0;5,3,0,0;
    5,4,0,0;5,5,0,0;
    5,6,0,0;6,0,0,-0.00173;
    6,1,-0.017283,-0.000337;6,2,0,0.000861;
    6,3,0,0;6,4,0,0;
    6,5,0,0;
    6,6,0,0;7,7,0,0;
];

% Initialize figure
figure;
hold on;

% Function to calculate K_T and K_Q
function [K_T, K_Q] = calculate_KT_KQ(J, P_D, coefficients)
    K_T = 0;
    K_Q = 0;
    for i = 1:size(coefficients, 1)
        x = coefficients(i, 1);
        y = coefficients(i, 2);
        Axy = coefficients(i, 3);
        Cxy = coefficients(i, 4);
        term_KT = Axy * (P_D)^x * J^y;
        term_KQ = Cxy * (P_D)^x * J^y;
        K_T = K_T + term_KT;
        K_Q = K_Q + term_KQ;
    end
end

% Store intersection points
intersections = [];

% Loop over each P/D value to calculate and plot K_T and K_Q
for pd = 1:length(P_D_values)
    P_D = P_D_values(pd); % Current value of P/D
    K_T_values = NaN(size(J_values));
    K_Q_values = NaN(size(J_values));
    eta_values = NaN(size(J_values)); % Efficiency values, initialized with NaN to ignore values out of range

    % Loop over each J value for the current P/D
    for j = 1:length(J_values)
        J = J_values(j); % Current value of J
        [K_T, K_Q] = calculate_KT_KQ(J, P_D, coefficients);

        % Store calculated values
        K_T_values(j) = K_T;
        K_Q_values(j) = K_Q * K_Q_scaling_factor;

        % Calculate efficiency η, avoid division by zero
        if K_Q ~= 0
            eta = (J * K_T) / (2 * pi * K_Q);
            if eta >= 0 && eta <= 1
                eta_values(j) = eta;
            end
        end
    end    

    % Find max efficiency and corresponding J for this P/D
    [max_eta, max_idx] = max(eta_values);
    if max_eta > best_efficiency
        % Calculate diameter and pitch under current P/D and J
        D_current = V_ship_mps / (J_values(max_idx));
        if D_current <= max_diameter
            % Update optimal values if within diameter limit
            best_efficiency = max_eta;
            optimal_P_D = P_D;
            optimal_J = J_values(max_idx);
            D_optimal = D_current;
            P_optimal = optimal_P_D * D_optimal;
        end
    end

    % Plot K_T, K_Q, and efficiency η vs J for the current P/D
    plot(J_values, K_T_values, ':','DisplayName', sprintf('K_T for P/D = %.1f', P_D), 'LineWidth', 2.5);
    plot(J_values, K_Q_values, '--', 'DisplayName', sprintf('K_Q for P/D = %.1f', P_D), 'LineWidth', 2.5);
    plot(J_values, eta_values, 'DisplayName', sprintf('Efficiency for P/D = %.1f', P_D), 'LineWidth', 2.5);

    % Find intersections of K_T and KTR
    KTR_values =  0.5861 * J_values.^2;
    [intersection_J, idx] = intersect_curves(J_values, K_T_values, KTR_values);
    for k = 1:length(intersection_J)
        J_intersect = intersection_J(k);
        K_T_intersect = KTR_values(idx(k));
        K_Q_intersect = K_Q_values(idx(k));
        eta_intersect = eta_values(idx(k));

        % Store the intersection data
        intersections = [intersections; P_D, J_intersect, K_T_intersect, K_Q_intersect, eta_intersect];

        % Annotate the intersection points
        plot(J_intersect, K_T_intersect, 'ro', 'MarkerSize', 8, 'DisplayName', 'Intersection Point');
        text(J_intersect, K_T_intersect, sprintf('J=%.2f', J_intersect), 'FontSize', 8, 'VerticalAlignment', 'bottom');
    end
end

% Add KTR curve
KTR_values = 0.5861 * J_values.^2;
plot(J_values, KTR_values, 'k-', 'LineWidth', 3, 'DisplayName', 'KTR');

% Display intersection points in a table
fprintf('\nIntersection Points:\n');
fprintf('P/D \t J \t K_T \t K_Q \t Efficiency\n');
for i = 1:size(intersections, 1)
    fprintf('%.2f \t %.2f \t %.4f \t %.4f \t %.4f\n', intersections(i, :));
end

% Add labels, title, and legend
xlabel('Advance Coefficient, J');
ylabel('Thrust and Torque Coefficients (K_T and K_Q)');
title('Thrust, Torque, and Efficiency vs. Advance Coefficient J');
ylim([0, 1]); % Set y-axis limits to [0, 1] for all variables
legend('show');
grid on;
hold off;

% Function to find intersections between two curves
function [intersection_J, idx] = intersect_curves(J_values, curve1, curve2)
    diff = curve1 - curve2;
    sign_change = diff(1:end-1) .* diff(2:end) < 0;
    idx = find(sign_change);
    intersection_J = (J_values(idx) + J_values(idx+1)) / 2;
end

