Adnan Almulla

public static (StationData<string> rebarTexts, StationValues providedAreas) FindStationRebarConfiguration
(int count, StationValues requiredAreas, string designCode, double beamArea)
{
    double minTotalExcess = double.MaxValue;

    // Stores the best configuration found: (Master Sequence, Layers for Left, Layers for Middle, Layers for Right)
    (List<int> masterSequence, int layersL, int layersM, int layersR)? bestConfig = null;
    StationValues bestProvidedAreas = StationValues.Zero;
    StationValues smallestDiameters = StationValues.Zero;
    StationValues largestDiameters = StationValues.Zero;

    // Generate "Master Sequences" of max length 4.
    // These sequences represent the potential deepest stack of rebar.
    var masterSequences = GetPermutations(_rebarDiameters, _maximumLayers);

    foreach (var sequence in masterSequences)
    {
        // For each station, find the minimum number of layers from this sequence
        // needed to satisfy the requirement.
        // The layers MUST be taken from the start of the sequence (continuity).

        var (validL, layersL, providedL) = FindMinLayersForStation(sequence, count, requiredAreas.Left, designCode, beamArea);
        var (validM, layersM, providedM) = FindMinLayersForStation(sequence, count, requiredAreas.Middle, designCode, beamArea);
        var (validR, layersR, providedR) = FindMinLayersForStation(sequence, count, requiredAreas.Right, designCode, beamArea);

        if (validL && validM && validR)
        {
            double totalExcess = (providedL - requiredAreas.Left) +
                                 (providedM - requiredAreas.Middle) +
                                 (providedR - requiredAreas.Right);

            if (totalExcess < minTotalExcess)
            {
                minTotalExcess = totalExcess;
                bestConfig = (sequence, layersL, layersM, layersR);
                bestProvidedAreas = new StationValues(providedL, providedM, providedR);
            }
        }
    }

    if (bestConfig == null)
    {
        return (new StationData<string>("Unable to Satisfy Requirements", "Unable to Satisfy Requirements", 
        "Unable to Satisfy Requirements"), StationValues.Zero, StationValues.Zero, StationValues.Zero);
    }

    var (seq, l, m, r) = bestConfig.Value;

    string TextFor(int numLayers)
    {
        // Take first 'numLayers' from 'seq', sort descending to formatting convention "Big + Small"
        var usedLayers = seq.Take(numLayers).OrderByDescending(d => d);
        return string.Join(" + ", usedLayers.Select(dia => $"{count}T{dia}"));
    }

    return (
        new StationData<string>(TextFor(l), TextFor(m), TextFor(r)),
        bestProvidedAreas,
    );
}

import numpy as np
import scipy  # pyright:ignore

#! Reference youtube video: https://www.youtube.com/watch?v=5csVr2bRKN4
#! Reference documentation: https://help.emd.dk/knowledgebase/content/ReferenceManual/ExtremeWinds.pdf

# ! Value as obtained by AOR.
mean_wind_velocity = 24  # m/s
std_velocity = 4.58  # m/s
eulers_constant = 0.5772

# Find beta and mu:
scale_beta = round((std_velocity * np.sqrt(6)) / np.pi, 2)  # m/s
location_mu = round(mean_wind_velocity - 0.577 * scale_beta, 2)  # m/s


def create_equation(years):
    """
    Creates the equation to solve based on number of years
    """

    def equation(Ue):
        #! The equation in the form where we want to find the root
        #! 1/years = 1 - exp(-exp(-(Ue-location)/scale))
        #! Rearranged to: 1/years - (1 - exp(-exp(-(Ue-location)/scale))) = 0
        return 1 / years - (1 - np.exp(-np.exp(-(Ue - location_mu) / scale_beta)))

    return equation


def solve_wind_speed(years, initial_guess=30):
    """
    Solves for maximum wind speed given a number of years

    Parameters:
    years (int): Number of years to calculate for
    initial_guess (float): Initial guess for the solver

    Returns:
    float: Maximum wind speed in m/s
    """
    equation = create_equation(years)

    # Use scipy's root finder to solve the equation
    solution = scipy.optimize.root(equation, initial_guess)

    if solution.success:
        return solution.x[0]
    else:
        raise ValueError("Failed to find solution")


# Solve for both 50 and 700 year periods
years_to_solve = [50, 700]

for years in years_to_solve:
    Ue = solve_wind_speed(years)
    print(f"Maximum wind speed for {years} years: {Ue:.1f} m/s")

ratio = round(
    solve_wind_speed(years_to_solve[1]) / solve_wind_speed(years_to_solve[0]), 2
)

print(f"Ratio of 700 year to 50 year annual exceedance: {ratio}")