%% Well Setup Function
function wellParams =  wellSetup(...
    dt, ...
    rodLengths, ...
    rodDiameters, ...
    rodMaterials, ...
    rodDampingFactors, ...
    tubingHeadPressure, ...
    stuffingBoxFriction, ...
    fluidGradient, ...
    pumpDiameter,...
    tubingOuterDiameter, ...
    tubingInnerDiameter, ...
    tubingAnchorDepth, ...
    structuralRating)

    % FUNCTION PARAMETERS
    % dt
    % rodLengths
    % rodDiameters
    % rodMaterials
    % tubingHeadPressure
    % pumpDiameter
    % tubingOuterDiameter
    % tubingInnerDiameter
    wellParams = struct();
    wellParams.dt = dt;
    wellParams.rodLengths = rodLengths;
    wellParams.rodDiameters = rodDiameters;
    wellParams.rodMaterials = rodMaterials;
    wellParams.rodDampingFactors = rodDampingFactors;
    wellParams.tubingHeadPressure = tubingHeadPressure;
    wellParams.stuffingBoxFriction = stuffingBoxFriction;
    wellParams.fluidGradient = fluidGradient;
    wellParams.pumpDiameter = pumpDiameter;
    wellParams.tubingOuterDiameter = tubingOuterDiameter;
    wellParams.tubingInnerDiameter = tubingInnerDiameter;
    wellParams.tubingAnchorDepth = tubingAnchorDepth;
    wellParams.structuralRating = structuralRating;
    
    numTapers = length(rodLengths);
    wellParams.numTapers = numTapers;

    
    
    wellParams.tubingCrossSectionalArea = (pi / 4) * ...
        (power(tubingOuterDiameter, 2) - ...
        power(tubingInnerDiameter,2));
    
    wellParams.pumpArea = power(pumpDiameter, 2)  * pi;
    
    % Calculated Taper Parameters
    area = zeros(1, numTapers+1);
    a = zeros(1, numTapers);
    pressure = zeros(1, numTapers);
    rodYMs = zeros(1, numTapers);
    rodWeightPerFoots = zeros(1, numTapers);
    force = zeros(1, numTapers);
    alpha = zeros(1, numTapers);
    stretch = zeros(1, numTapers);
    xOverA = zeros(1, numTapers);
    lagIndex = zeros(1, numTapers, 'uint16');
    factorArray = zeros(1, numTapers);
    centerPoint = zeros(1, numTapers, 'uint16');
    lengthRequired = zeros(1, numTapers, 'uint16');
    
    
    % Calculated Rod String Totals
    rodDepthTotal = 0.0;
    buoyantForceTotal = 0.0;
    rodWeightAirTotal = 0.0;
    

    for i = 1:numTapers
        area(i) = pi / 4 * power(rodDiameters(i), 2);
        
        if lower(rodMaterials(i)) == "fiberglass"
            rodYMs(i) = 7.2;
        else
            rodYMs(i) = 30.5;
        end
        
        rodWeightPerFoots(i) = lookupRodWeightPerFoot(rodYMs(i), ...
            rodDiameters(i));
    end
    wellParams.area = area;
    wellParams.rodYMs = rodYMs;
    wellParams.rodWeightPerFoots = rodWeightPerFoots;

    for i = 1:numTapers
        a(i) = 1000.0 * sqrt(32.2 * rodYMs(i) * area(i) / rodWeightPerFoots(i));
        rodDepthTotal = rodDepthTotal + rodLengths(i);

        if i > 1
            pressure(i) = pressure(i - 1) + fluidGradient * rodLengths(i);
        else
            pressure(i) = tubingHeadPressure + fluidGradient * rodLengths(i);
        end

        buoyantForceTotal = buoyantForceTotal + pressure(i) * (area(i+1) - area(i));
        rodWeightAirTotal = rodWeightAirTotal + rodWeightPerFoots(i) * rodLengths(i);
        rodWeightFluidTotal = rodWeightAirTotal + buoyantForceTotal;
    end
    
    wellParams.a = a;
    wellParams.rodDepthTotal = rodDepthTotal;
    wellParams.pressure = pressure;
    wellParams.buoyantForceTotal = buoyantForceTotal;
    wellParams.rodWeightAirTotal = rodWeightAirTotal;
    wellParams.rodWeightFluidTotal = rodWeightFluidTotal;
    

    for j = 1:numTapers
        weightData = 0.0;
        annularForceData = 0.0;

        for i = j+1:numTapers
            weightData = weightData + rodWeightPerFoots(i) * rodLengths(i);
        end

        for i = j:numTapers-1
            annularForceData = annularForceData + pressure(i) * (area(i) - area(i+1));
        end
        
        force(j) = (-area(numTapers) * pressure(numTapers)) + weightData - annularForceData;
        alpha(j) = (force(j) + rodWeightPerFoots(j) * rodLengths(j)) / (rodYMs(j) * power(10, 6) * area(j));

        if j > 1
            stretch(j) = stretch(j - 1) + alpha(j) * rodLengths(j) - ...
                (rodWeightPerFoots(j) * power(rodLengths(j), 2.0)) / ...
                (2.0 * rodYMs(j) * power(10, 6) * area(j));
        else
            stretch(j) = 0.0 + alpha(j) * rodLengths(j) - ...
                (rodWeightPerFoots(j) * power(rodLengths(j), 2.0)) / ...
                (2.0 * rodYMs(j) * power(10, 6) * area(j));
        end
    end
    wellParams.force = force;
    wellParams.alpha = alpha;
    wellParams.stretch = stretch;

    for i = 1:numTapers
        xOverA(i) = rodLengths(i) / a(i);
        lagIndex(i) = floor(rodLengths(i) / (a(i) * dt));
        factorArray(i) = (xOverA(i) - double(lagIndex(i)) * dt) / dt;
        centerPoint(i) = lagIndex(i) + 2;
        lengthRequired(i) = 2 * (lagIndex(i) + 1) + 1;
    end
    
    wellParams.xOverA = xOverA;
    wellParams.lagIndex = lagIndex;
    wellParams.factorArray = factorArray;
    wellParams.centerPoint = centerPoint;
    wellParams.lengthRequired = lengthRequired;

    frictionEstimate = rodDepthTotal * 0.10;
    wellParams.frictionEstimate = frictionEstimate;
end

%% Determines Rod Weight Per Foot given Young's Modulus and Diameter
function wtPerFt = lookupRodWeightPerFoot(i_ym, i_diam)
    YM_STEEL=30.5;
    YM_FIBERGLASS=7.2;
    wtPerFt = -1;
    if (i_ym == YM_STEEL)
        if (i_diam <= 2 && i_diam > 1.75)
            wtPerFt = 10.7;
        elseif (i_diam <= 1.75 && i_diam > 1.65)
            wtPerFt = 8.2;
        elseif (i_diam <= 1.65 && i_diam > 1.5)
            wtPerFt = 7;
        elseif (i_diam <= 1.5 && i_diam > 1.375)
            wtPerFt = 6;
        elseif (i_diam <= 1.375 && i_diam > 1.125)
            wtPerFt = 5;
        elseif (i_diam <= 1.125 && i_diam > 1)
            wtPerFt = 3.676;
        elseif (i_diam <= 1 && i_diam > 0.875)
            wtPerFt = 2.904;
        elseif (i_diam <= 0.875 && i_diam > 0.75)
            wtPerFt = 2.224;
        elseif (i_diam <= 0.75 && i_diam > 0.625)
            wtPerFt = 1.634;
        elseif (i_diam <= 0.625 && i_diam > 0.5)
            wtPerFt = 1.13;
        elseif (i_diam <= 0.5)
            wtPerFt = 0.72;
        else
            wtPerFt = 0;
        end
    elseif (i_ym == YM_FIBERGLASS)
        if (i_diam <= 1.25 && i_diam > 1.125)
            wtPerFt = 1.2879;
        elseif (i_diam <= 1.125 && i_diam > 1)
            wtPerFt = 1.09;
        elseif (i_diam <= 1 && i_diam > 0.875)
            wtPerFt = 0.8188;
        elseif (i_diam <= 0.875 && i_diam > 0.75)
            wtPerFt = 0.6108;
        elseif (i_diam <= 0.75)
            wtPerFt = 0.484;
        else
            wtPerFt = 0;
        end
    end
end