
%EXAMPLE CALL: findOptimalParameters(.5:.5:10, 0:.01:.5, 0:.01:.25, .001, 2*pi/3, 0.01, 10, .02, "oval")

function findOptimalParameters(freqs, lengths, widths, thick, swath, turnTime, wingDens, loadMass, shape)

    lenCombos = repelem(lengths, 1, length(widths) * length(freqs));
    widCombos = repelem(widths, 1, length(freqs));
    widCombos = repmat(widCombos, [1 length(lengths)]);
    freCombos = repmat(freqs, [1 (length(lengths) * length(widths))]);

    forces = calcNetForce(freCombos, lenCombos, widCombos, thick, swath, turnTime, wingDens, loadMass, shape);

    powers = calcPower(freCombos, lenCombos, widCombos, thick, swath, turnTime, shape);
    
    normDif = calcNormalizedDifference(forces, powers);
    
    optimalIndex = normDif == max(normDif);

    fprintf('Optimal parameters at:\n%0.4f %s\n%0.4f %s\n%0.4f %s\n%0.4f %s\n%s %0.2f %s\n\n',...
        swath, 'Radians Stroke Angle,',...
        thick, 'm Wing Thickness,',...
        wingDens, 'kg/m^3 wing Density,',...
        turnTime, 's proation/supination rotation time,',...
        'and', loadMass, 'kg Load Mass');

    fprintf('Optimal Width: %0.2f m\nOptimal Length: %0.2f m\nOptimal Frequency: %0.1f Hz\n\n',...
        widCombos(optimalIndex), lenCombos(optimalIndex), freCombos(optimalIndex));
    
    fprintf('Upward Force Generated: %0.2f N\nPower Required: %0.2f W\n\n\n',...
        forces(optimalIndex), powers(optimalIndex));
    
    if(shape == "rectangle")
        coef = 1;
    elseif(shape == "oval")
        coef = pi/4;
    end

    subplot(2,2,1);

    scatter3(coef * lenCombos .* widCombos, freCombos, forces, [], forces, 'filled');
    title('Net Upward Force per Cycle');

    xlabel('Area (m^2)');
    ylabel('Frequency (Hz)');
    zlabel('Force (N)');

    subplot(2,2,2);

    scatter3(coef * lenCombos .* widCombos, freCombos, powers, [], powers, 'filled');
    title('Power per Cycle');
    xlabel('Area (m^2)');
    ylabel('Frequency (Hz)');
    zlabel('Power (W)');

    subplot(2,4,6:7);

    scatter3(coef * lenCombos .* widCombos, freCombos, normDif, [], normDif, 'filled');
    title('Force-Power Normalized Difference');
    xlabel('Area (m^2)');
    ylabel('Frequency (Hz)');
    zlabel('Difference');
end

function netForces =...
      calcNetForce(freCombos, lenCombos, widCombos, thick, swath, turnTime, wingDens, loadMass, shape)
  
    G = 9.81; %m/s^2
    AIR_DENS = 1.225; %kg/m^3
    WING_COEF_DRAG = 2; %flat wing perpendicular to flow

    if (turnTime / 2.0 >= 1.0 / min(freCombos))
        disp(['Too slow of a turn time for the given frequencies.\n',...
            'The maximum frequency at the given turn time is ',...
            sprintf('%0.2f', 1.0  / (2 * turnTime))]);
    end

    swathInteg = (96*pi^4*sin(swath))/(swath^4-20*pi^2*swath^2+64*pi^4);
    coef = swathInteg * 4 * swath^2 * AIR_DENS * WING_COEF_DRAG;

    if(shape == "rectangle")
        upDownDif = (lenCombos .^ 3 / 3) .* (widCombos - thick);
        coef = 1;
    elseif(shape == "oval")
        upDownDif = (lenCombos .^ 3 / 192) .* (15 * pi .* widCombos - 64 * thick);
        coef = pi/4;
    end

    drag = coef * (freCombos ./ (1 - 2 * turnTime * freCombos)) .^ 2;
    weight = 2 * coef * turnTime * 1.0 ./ freCombos * G .* (loadMass + 2 * wingDens * thick * lenCombos .* widCombos);
    netForces = drag .* upDownDif - weight;

end

function powers = calcPower(freCombos, lenCombos, widCombos, thick, swath, turnTime, shape)
    AIR_DENS = 1.225; %kg/m^3
    WING_COEF_DRAG = 2; %flat wing perpendicular to flow
    
    if(turnTime <= 0)
        disp('Turn time cannot be zero or less than zero.');
    end
    
    coef = 5 * swath / 2.0 * 8 * swath^3 * AIR_DENS * WING_COEF_DRAG;

    if(shape == "rectangle")
        upDownCombo = (lenCombos .^ 3 / 3) .* (widCombos + thick);
        rotInteg = lenCombos .* widCombos .^ 4 / 128;
    elseif(shape == "oval")
        upDownCombo = (lenCombos .^ 3 / 192) .* (15 * pi .* widCombos + 64 * thick);
        rotInteg = lenCombos .* widCombos .^ 4 / 120;
    end
    
    powRot = 2 * AIR_DENS * WING_COEF_DRAG * (pi / turnTime) ^ 3 * rotInteg;
    powStroke = coef * (freCombos ./ (1 - 2 * turnTime * freCombos)) .^ 3;
    powers = powRot + powStroke .* upDownCombo;
end

function normDif = calcNormalizedDifference(forces, powers)
    normForces = (forces - min(forces)) / (max(forces) - min(forces));
    normPowers = (powers - min(powers)) / (max(powers) - min(powers));
    dif = normForces - normPowers;
    normDif = (dif - min(dif)) / (max(dif) - min(dif));
end
