app.py

import gradio as gr
from pint import UnitRegistry, set_application_registry
import matplotlib.pyplot as plt
import io
import base64
from sympy import symbols, Symbol, Eq, solve
from reportlab.lib.pagesizes import letter
from reportlab.pdfgen import canvas
from PIL import Image
import base64
import io
from reportlab.lib.utils import ImageReader
import os
import subprocess
import contextlib


# Initialize unit registry
u = UnitRegistry()
set_application_registry(u)
Q_ = u.Quantity

def generate_beam_diagram(n_spans, lengths, loads, moments, R_sx, R_dx,
                          cantilever_left_length=0 * u.meter, cantilever_left_load=0 * u.newton / u.meter,
                          cantilever_right_length=0 * u.meter, cantilever_right_load=0 * u.newton / u.meter,
                          unit_system='SI'):
    # Define units and their abbreviations based on the unit system
    if unit_system == 'SI':
        length_unit = u.meter
        load_unit = u.newton / u.meter
        moment_unit = u.newton * u.meter
        reaction_unit = u.newton
        length_unit_str = 'm'
        load_unit_str = 'N/m'
        moment_unit_str = 'N·m'
        reaction_unit_str = 'N'
        length_display_unit = u.meter  # For consistent plotting
    else:
        length_unit = u.foot
        load_unit = u.pound_force / u.foot
        moment_unit = u.pound_force * u.foot
        reaction_unit = u.pound_force
        length_unit_str = 'ft'
        load_unit_str = 'lb/ft'
        moment_unit_str = 'lb·ft'
        reaction_unit_str = 'lb'
        length_display_unit = u.foot  # For consistent plotting

    # Adjust lengths to include cantilevers
    lengths_display = []
    x_positions = [0]
    # Left cantilever
    if cantilever_left_length.magnitude > 0:
        l_cant_left_display = cantilever_left_length.to(length_display_unit)
        lengths_display.append(l_cant_left_display)
        x_positions.append(x_positions[-1] + l_cant_left_display.magnitude)
    # Main spans
    for l in lengths:
        l_display = l.to(length_display_unit)
        lengths_display.append(l_display)
        x_positions.append(x_positions[-1] + l_display.magnitude)
    # Right cantilever
    if cantilever_right_length.magnitude > 0:
        l_cant_right_display = cantilever_right_length.to(length_display_unit)
        lengths_display.append(l_cant_right_display)
        x_positions.append(x_positions[-1] + l_cant_right_display.magnitude)

    total_length = x_positions[-1]

    fig, ax = plt.subplots(figsize=(10, 3))
    # Draw the beam as a horizontal line
    ax.hlines(0, 0, total_length, colors='black', linewidth=2)

    # -------------------------------------------
    # Draw supports (rotated 180° so the tip is on the beam)
    # -------------------------------------------
    support_width = total_length / 50  # Adjust support width relative to total length
    n_supports = n_spans + 1
    support_positions = []
    idx = 0
    # Left support
    if cantilever_left_length.magnitude > 0:
        # Support at the end of the left cantilever (x_positions[1])
        x = x_positions[1]
        support_positions.append(x)
        support = plt.Polygon([[x - support_width, -0.2],
                               [x + support_width, -0.2],
                               [x, 0]], color='black')
        ax.add_patch(support)
        idx = 2
    else:
        # Support at start (x_positions[0])
        x = x_positions[0]
        support_positions.append(x)
        support = plt.Polygon([[x - support_width, -0.2],
                               [x + support_width, -0.2],
                               [x, 0]], color='black')
        ax.add_patch(support)
        idx = 1

    # Intermediate supports
    for i in range(1, n_supports - 1):
        x = x_positions[idx]
        support_positions.append(x)
        support = plt.Polygon([[x - support_width, -0.2],
                               [x + support_width, -0.2],
                               [x, 0]], color='black')
        ax.add_patch(support)
        idx += 1

    # Right support
    if cantilever_right_length.magnitude > 0:
        # Support before the right cantilever (x_positions[-2])
        x = x_positions[-2]
        support_positions.append(x)
        support = plt.Polygon([[x - support_width, -0.2],
                               [x + support_width, -0.2],
                               [x, 0]], color='black')
        ax.add_patch(support)
    else:
        # Support at end (x_positions[-1])
        x = x_positions[-1]
        support_positions.append(x)
        support = plt.Polygon([[x - support_width, -0.2],
                               [x + support_width, -0.2],
                               [x, 0]], color='black')
        ax.add_patch(support)

    # -------------------------------------------
    # Annotate the TOTAL reaction force (sum of left and right) at each support
    # -------------------------------------------
    for i, x in enumerate(support_positions):
        R_total = R_sx[i] + R_dx[i]
        R_total_conv = R_total.to(reaction_unit)
        R_total_num = round(R_total_conv.magnitude, 6)
        ax.text(x, 0.25, f'$R_{{{i+1}}} = {R_total_num}$ {reaction_unit_str}',
                ha='center', va='bottom', color='green', fontsize=10)

    # Annotate internal moments below the beam
    for i, x in enumerate(support_positions):
        M_value = moments[i]
        M_unit = M_value.to(moment_unit)
        M_num = round(M_unit.magnitude, 6)
        ax.text(x, -0.25, f'$M_{{{i+1}}} = {M_num}$ {moment_unit_str}',
                ha='center', va='top', color='blue', fontsize=10)

    # Remove axes and adjust limits
    ax.axis('off')
    plt.xlim(-0.05 * total_length, total_length * 1.05)
    plt.ylim(-0.4, 0.4)

    # Save plot to a buffer and encode as base64
    buf = io.BytesIO()
    plt.savefig(buf, format='png', bbox_inches='tight', dpi=150)
    plt.close(fig)
    buf.seek(0)
    img_base64 = base64.b64encode(buf.getvalue()).decode('utf-8')
    return f'<img src="data:image/png;base64,{img_base64}" alt="Beam Diagram"/>'

def calculate_reactions(n_spans, lengths, distributed_loads, moments,
                        cantilever_left_length=Q_(0.0, u.meter), cantilever_left_load=Q_(0.0, u.newton / u.meter),
                        cantilever_right_length=Q_(0.0, u.meter), cantilever_right_load=Q_(0.0, u.newton / u.meter)):
    """
    Calculate reactions.
    For a cantilever the reaction force is simply:
         R = (distributed load) * (cantilever length)
    """
    n_supports = n_spans + 1

    # Print all torques (moments) at supports
    print("\nTorques at each support:")
    for i, moment in enumerate(moments):
        print(f"Torque {i+1}: {moment}")
    print("\n")

    # Initialize reaction lists (each support will have left (R_sx) and right (R_dx) components)
    R_sx = [Q_(0.0, 'newton') for _ in range(n_supports)]
    R_dx = [Q_(0.0, 'newton') for _ in range(n_supports)]

    # ---------------------------
    # First support (m = 0)
    # ---------------------------
    print("Calculating R1:")
    # For a left cantilever, the left reaction is the cantilever force:
    if cantilever_left_length.magnitude > 0:
        R_sx[0] = cantilever_left_load * cantilever_left_length
        print(f"R1_sx = cantilever_left_load * cantilever_left_length = {R_sx[0]}")
        
    # r1_dx = (l0 * w0)/2 + (M1 - M2)/l0
    R_dx[0] = (distributed_loads[0] * lengths[0]) / 2 + (moments[0] - moments[1]) / lengths[0]
    print(f"R1_dx = ({distributed_loads[0]} * {lengths[0]})/2 + ({moments[0]} - {moments[1]})/{lengths[0]} = {R_dx[0]}")

    # ---------------------------
    # Middle supports (1 .. n-1)
    # ---------------------------
    for m in range(1, n_supports - 1):
        print(f"\nCalculating R{m+1}:")
        # Left reaction at support m: load from previous span minus the right reaction of previous support
        R_sx[m] = distributed_loads[m-1] * lengths[m-1] - R_dx[m-1]
        print(f"R{m+1}_sx = {distributed_loads[m-1]} * {lengths[m-1]} - {R_dx[m-1]} = {R_sx[m]}")
        
        # Right reaction at support m: half the load on the next span plus moment difference contribution
        R_dx[m] = (distributed_loads[m] * lengths[m]) / 2 + (moments[m] - moments[m+1]) / lengths[m]
        print(f"R{m+1}_dx = ({distributed_loads[m]} * {lengths[m]})/2 + ({moments[m]} - {moments[m+1]})/{lengths[m]} = {R_dx[m]}")

    # ---------------------------
    # Last support (m = n_supports - 1)
    # ---------------------------
    m = n_supports - 1
    print(f"\nCalculating R{m+1}:")
    # Left reaction at last support from the previous span
    R_sx[m] = distributed_loads[m-1] * lengths[m-1] - R_dx[m-1]
    print(f"R{m+1}_sx = {distributed_loads[m-1]} * {lengths[m-1]} - {R_dx[m-1]} = {R_sx[m]}")

    # For a right cantilever, the right reaction is the cantilever force;
    # otherwise it is set to zero.
    if cantilever_right_length.magnitude > 0:
        R_dx[m] = cantilever_right_load * cantilever_right_length
        print(f"R{m+1}_dx = cantilever_right_load * cantilever_right_length = {R_dx[m]}")
    else:
        R_dx[m] = Q_(0.0, 'newton')
        print(f"R{m+1}_dx is forced to 0: {R_dx[m]}")

    # Print final reactions summary
    print("\nFinal Reactions Summary:")
    for i in range(n_supports):
        print(f"R{i+1}_sx = {R_sx[i]}")
        print(f"R{i+1}_dx = {R_dx[i]}")
        print(f"R{i+1}_total = {R_sx[i] + R_dx[i]}\n")

    return R_sx, R_dx

def continuous_beam_solver(unit_system, n_spans, lengths_str, loads_str,
                             cantilever_left_length_str='0', cantilever_left_load_str='0',
                             cantilever_right_length_str='0', cantilever_right_load_str='0'):
    # Parse the input strings into lists
    try:
        n_spans = int(n_spans)
        l_values = [float(val.strip()) for val in lengths_str.split(',')]
        p_values = [float(val.strip()) for val in loads_str.split(',')]
        cant_left_len_val = float(cantilever_left_length_str.strip())
        cant_left_load_val = float(cantilever_left_load_str.strip())
        cant_right_len_val = float(cantilever_right_length_str.strip())
        cant_right_load_val = float(cantilever_right_load_str.strip())
    except ValueError:
        return "Invalid input. Please ensure all inputs are numbers.", "", "", "", ""

    if len(l_values) != n_spans or len(p_values) != n_spans:
        return "The number of lengths and loads must match the number of spans.", "", "", "", ""

    n_supports = n_spans + 1  # Total number of supports

    # Define units based on the selected unit system
    if unit_system == 'SI':
        length_unit = u.meter
        load_unit = u.newton / u.meter
        moment_unit = u.newton * u.meter
        reaction_unit = u.newton
    else:
        length_unit = u.foot
        load_unit = u.pound_force / u.foot
        moment_unit = u.pound_force * u.foot
        reaction_unit = u.pound_force

    # Convert lengths and loads to quantities with units
    l = [Q_(length, length_unit) for length in l_values]
    p = [Q_(load, load_unit) for load in p_values]

    # Convert lengths and loads to SI units for calculation
    l_SI = [length.to(u.meter) for length in l]
    p_SI = [load.to(u.newton / u.meter) for load in p]

    # Process cantilever inputs
    cantilever_left_length = Q_(cant_left_len_val, length_unit)
    cantilever_left_load = Q_(cant_left_load_val, load_unit)
    cantilever_right_length = Q_(cant_right_len_val, length_unit)
    cantilever_right_load = Q_(cant_right_load_val, load_unit)

    # Convert cantilever inputs to SI units for calculations
    cantilever_left_length_SI = cantilever_left_length.to(u.meter)
    cantilever_left_load_SI = cantilever_left_load.to(u.newton / u.meter)
    cantilever_right_length_SI = cantilever_right_length.to(u.meter)
    cantilever_right_load_SI = cantilever_right_load.to(u.newton / u.meter)

    # Initialize moments at supports
    M_values_SI = []
    M_values_Imperial = []

    # Handle single-span beam separately
    if n_spans == 1:
        # Cantilever moment contributions (if any)
        M_cantilever_left = 0.0
        if cantilever_left_length_SI.magnitude > 0 and cantilever_left_load_SI.magnitude != 0:
            M_cantilever_left = (cantilever_left_load_SI * cantilever_left_length_SI**2 / 2).to(u.newton * u.meter).magnitude

        M_cantilever_right = 0.0
        if cantilever_right_length_SI.magnitude > 0 and cantilever_right_load_SI.magnitude != 0:
            M_cantilever_right = (cantilever_right_load_SI * cantilever_right_length_SI**2 / 2).to(u.newton * u.meter).magnitude

        # Assign moments for supports A and B
        M_A = M_cantilever_left
        M_B = M_cantilever_right

        M_values_SI = [Q_(M_A, u.newton * u.meter), Q_(M_B, u.newton * u.meter)]
        M_values_Imperial = [M_values_SI[0].to(u.pound_force * u.foot), M_values_SI[1].to(u.pound_force * u.foot)]

        # Calculate reactions+
        f = io.StringIO()
        with contextlib.redirect_stdout(f):
            R_sx_SI, R_dx_SI = calculate_reactions(n_spans, l_SI, p_SI, M_values_SI)
        reaction_log = f.getvalue()

        R_sx_Imperial = [R.to(u.pound_force) for R in R_sx_SI]
        R_dx_Imperial = [R.to(u.pound_force) for R in R_dx_SI]

        results_SI = ""
        results_Imperial = ""
        for i in range(n_supports):
            results_SI += f"M_{i+1} = {M_values_SI[i].magnitude:.6f} N·m\n"
            results_Imperial += f"M_{i+1} = {M_values_Imperial[i].magnitude:.6f} lb·ft\n"

        beam_diagram_SI = generate_beam_diagram(
            n_spans, l, p, M_values_SI, R_sx_SI, R_dx_SI,
            cantilever_left_length=cantilever_left_length,
            cantilever_left_load=cantilever_left_load,
            cantilever_right_length=cantilever_right_length,
            cantilever_right_load=cantilever_right_load,
            unit_system='SI'
        )
        beam_diagram_Imperial = generate_beam_diagram(
            n_spans, l, p, M_values_Imperial, R_sx_Imperial, R_dx_Imperial,
            cantilever_left_length=cantilever_left_length,
            cantilever_left_load=cantilever_left_load,
            cantilever_right_length=cantilever_right_length,
            cantilever_right_load=cantilever_right_load,
            unit_system='Imperial'
        )
        # For single-span, no complex equations need to be shown.
        equations_md = ""
        return equations_md, results_SI, results_Imperial, beam_diagram_SI, beam_diagram_Imperial, reaction_log

    # For multiple spans
    else:
        # Compute fixed-end moments due to cantilever loads
        M_cantilever_left = 0.0
        if cantilever_left_length_SI.magnitude > 0 and cantilever_left_load_SI.magnitude != 0:
            M_cantilever_left = (cantilever_left_load_SI * cantilever_left_length_SI**2 / 2).to(u.newton * u.meter).magnitude

        M_cantilever_right = 0.0
        if cantilever_right_length_SI.magnitude > 0 and cantilever_right_load_SI.magnitude != 0:
            M_cantilever_right = (cantilever_right_load_SI * cantilever_right_length_SI**2 / 2).to(u.newton * u.meter).magnitude

        # Initialize moments M_i (M_1 to M_n)
        M_symbols = []
        M_symbols.append(M_cantilever_left)  # M_1

        for i in range(1, n_supports - 1):
            M_symbols.append(symbols(f'M_{i+1}'))  # M_2 to M_{n_supports-1}

        M_symbols.append(M_cantilever_right)  # M_n

        # Set up the system of equations (for supports 2 to n-1)
        equations = []
        equations_latex = []
        for k in range(1, n_supports - 1):
            l_prev = l_SI[k - 1].magnitude
            l_curr = l_SI[k].magnitude
            p_prev = p_SI[k - 1].magnitude
            p_curr = p_SI[k].magnitude
            M_prev = M_symbols[k - 1]
            M_curr = M_symbols[k]
            M_next = M_symbols[k + 1]

            lhs = (1/24) * (l_prev**3 * p_prev + l_curr**3 * p_curr)
            rhs = (1/6) * (l_prev * M_prev + l_curr * M_next) + (1/3) * (l_prev + l_curr) * M_curr
            equation = Eq(lhs, rhs)
            equations.append(equation)
            equation_latex = f"\\frac{{1}}{{24}}(l_{{{k}}}^3 p_{{{k}}} + l_{{{k+1}}}^3 p_{{{k+1}}}) = \\frac{{1}}{{6}}(l_{{{k}}} M_{{{k}}} + l_{{{k+1}}} M_{{{k+2}}}) + \\frac{{1}}{{3}}(l_{{{k}}} + l_{{{k+1}}}) M_{{{k+1}}}"
            equations_latex.append(equation_latex)

        # Solve the system for the unknown moments
        unknown_M_symbols = [M_symbols[i] for i in range(1, n_supports - 1) if isinstance(M_symbols[i], Symbol)]
        solution = solve(equations, unknown_M_symbols, dict=True)

        if solution:
            solution = solution[0]
            results_SI = ""
            results_Imperial = ""
            M_values_SI = []
            M_values_Imperial = []
            for i in range(n_supports):
                M_i_value = M_symbols[i]
                if isinstance(M_i_value, Symbol):
                    M_i_value = float(solution.get(M_i_value, 0))
                else:
                    M_i_value = float(M_i_value)
                M_quantity_SI = Q_(M_i_value, u.newton * u.meter)
                M_values_SI.append(M_quantity_SI)
                results_SI += f"M_{i+1} = {M_quantity_SI.magnitude:.6f} N·m\n"
                M_quantity_Imperial = M_quantity_SI.to(u.pound_force * u.foot)
                M_values_Imperial.append(M_quantity_Imperial)
                results_Imperial += f"M_{i+1} = {M_quantity_Imperial.magnitude:.6f} lb·ft\n"
        else:
            return "No solution found.", "", "", "", ""

        # Calculate reactions
        f = io.StringIO()
        with contextlib.redirect_stdout(f):
            R_sx_SI, R_dx_SI = calculate_reactions(n_spans, l_SI, p_SI, M_values_SI,
                                                    cantilever_left_length=cantilever_left_length_SI,
                                                    cantilever_left_load=cantilever_left_load_SI,
                                                    cantilever_right_length=cantilever_right_length_SI,
                                                    cantilever_right_load=cantilever_right_load_SI)
        reaction_log = f.getvalue()

        R_sx_Imperial = [R.to(u.pound_force) for R in R_sx_SI]
        R_dx_Imperial = [R.to(u.pound_force) for R in R_dx_SI]

        beam_diagram_SI = generate_beam_diagram(
            n_spans, l, p, M_values_SI, R_sx_SI, R_dx_SI,
            cantilever_left_length=cantilever_left_length,
            cantilever_left_load=cantilever_left_load,
            cantilever_right_length=cantilever_right_length,
            cantilever_right_load=cantilever_right_load,
            unit_system='SI')
        beam_diagram_Imperial = generate_beam_diagram(
            n_spans, l, p, M_values_Imperial, R_sx_Imperial, R_dx_Imperial,
            cantilever_left_length=cantilever_left_length,
            cantilever_left_load=cantilever_left_load,
            cantilever_right_length=cantilever_right_length,
            cantilever_right_load=cantilever_right_load,
            unit_system='Imperial')

        equations_md = "\n\n".join(
            [f"**Equation {i+1}:**\n\n$$ {eq} $$" for i, eq in enumerate(equations_latex)]
        )
        return equations_md, results_SI, results_Imperial, beam_diagram_SI, beam_diagram_Imperial, reaction_log


        

def gradio_interface(unit_system, n_spans, lengths_str, loads_str,
                     cantilever_left_length, cantilever_left_load,
                     cantilever_right_length, cantilever_right_load):
    # Call the continuous beam solver to obtain all outputs including the reaction log.
    (equations_md, results_SI, results_Imperial,
     beam_diagram_SI, beam_diagram_Imperial, reaction_log) = continuous_beam_solver(
        unit_system, n_spans, lengths_str, loads_str,
        cantilever_left_length, cantilever_left_load,
        cantilever_right_length, cantilever_right_load)
    # Return outputs in the new order: equations, moments, diagrams, and then the reaction log.
    return equations_md, results_SI, results_Imperial, beam_diagram_SI, beam_diagram_Imperial, reaction_log


# Build the Gradio interface.
# Note that the input labels now indicate the expected units.
iface = gr.Interface(
    fn=gradio_interface,
    inputs=[
        gr.Radio(['SI', 'Imperial'], label="Unit System", value='SI'),
        gr.Number(label="Number of Spans (n)", value=3, precision=0),
        gr.Textbox(label="Lengths l_i (comma-separated) [m (SI) or ft (Imperial)]", 
                   placeholder="e.g., 7.92, 7.92, 7.92", value='7.91667,7.91667,7.91667'),
        gr.Textbox(label="Loads p_i (comma-separated) [N/m (SI) or lb/ft (Imperial)]", 
                   placeholder="e.g., 200,200,200", value='200,200,200'),
        gr.Textbox(label="Cantilever Left Length [m (SI) or ft (Imperial)]", 
                   placeholder="e.g., 6.67", value='6.66667'),
        gr.Textbox(label="Cantilever Left Load [N/m (SI) or lb/ft (Imperial)]", 
                   placeholder="e.g., 200", value='200'),
        gr.Textbox(label="Cantilever Right Length [m (SI) or ft (Imperial)]", 
                   placeholder="e.g., 6.67", value='6.66667'),
        gr.Textbox(label="Cantilever Right Load [N/m (SI) or lb/ft (Imperial)]", 
                   placeholder="e.g., 200", value='200'),
    ],
    outputs=[
        gr.Markdown(label="Equations Used"),
        gr.Textbox(label="Internal Moments at Supports (SI Units)"),
        gr.Textbox(label="Internal Moments at Supports (Imperial Units)"),
        gr.HTML(label="Beam Diagram (SI Units)"),
        gr.HTML(label="Beam Diagram (Imperial Units)"),
        gr.Textbox(label="Reactions Calculation Log"),
    ],
    title="Continuous Beam Solver with Cantilevers",
    description=(
        "Solve for internal moments at supports of a continuous beam with multiple spans, including cantilevers.\n\n"
        "**Input Units:**\n"
        "- For SI: Lengths in meters (m), Loads in Newtons per meter (N/m), Moments in N·m, Reaction forces in N.\n"
        "- For Imperial: Lengths in feet (ft), Loads in pounds per foot (lb/ft), Moments in lb·ft, Reaction forces in lb.\n\n"
        "The outputs are arranged with the equations on top, followed by the resulting internal moments, "
        "then the beam diagram, and finally the complete reaction calculation log."
    ),
    allow_flagging="never",
)

if __name__ == "__main__":
    iface.launch()