#include <stdio.h>    // For standard input/output functions like printf, scanf
#include <stdlib.h>   // For general utilities like exit, abs
#include <string.h>   // For string manipulation functions like strcspn, fgets
#include <math.h>     // For mathematical functions like round, fabs, pow
#include <ctype.h>    // For character handling functions like toupper
#include <time.h>     // For time-related functions like time, localtime, strftime

// Constants
#define PUF 96 // plotter unit factor (dots per inch) - used for scaling plots

// Global Variables
double d[51][51];   // Matrix to solve for joint rotations and displacements (1-based indexing in Pascal, so size +1 for C to match 0-based)
double b[41][7];    // Beam Data: b[beam_idx][1]=Moment of Inertia (I), b[beam_idx][2]=Length, b[beam_idx][3]=FEM (Fixed End Moment) Left, b[beam_idx][4]=FEM Right, b[beam_idx][5]=Shear Left, b[beam_idx][6]=Shear Right
double c[41][3];    // Column Data: c[col_idx][1]=Moment of Inertia (I), c[col_idx][2]=Length (height)
double e[11][5];    // Wind loading data (FEMs and shears for lateral loads on columns) e[story_idx][1]=FEM Top (L), e[story_idx][2]=FEM Bottom (R), e[story_idx][3]=Shear Top, e[story_idx][4]=Shear Bottom
double ld[51][5];   // Load Data: ld[load_idx][1]=Story # (0 for lateral), ld[load_idx][2]=Intensity (W), ld[load_idx][3]=Location (A), ld[load_idx][4]=Length (C)
int mp[21][3];      // Macro Pointer: mp[macro_idx][1]=start_load_idx, mp[macro_idx][2]=end_load_idx (defines a group of loads)
double r[401];      // Array for storing calculated moments/deflections for plotting (0-based)

int s, p;           // s: # stories, p: # portals
int q;              // size of matrix to solve (number of unknown degrees of freedom)
int z;              // type of base joints: 1 (Fixed) or 2 (Pinned)
double l2, l3;      // l2: total frame width (total beam length), l3: total frame height (total column height)

int lp;             // Load Pointer: points to the next available row in the `ld` (load data) matrix
int mc;             // Macro Counter: counts the number of defined macro loads
int lc;             // Load Count: total number of individual loads applied during current analysis run
int lm[51];         // Load Member: saves the load numbers (indices from `ld`) of all loads applied in the current run
int psl;            // Plot Scale Length (feet/inch)
int psm;            // Plot Scale Moment/Shear (kip-ft/inch or kips/inch)
char proj_descrip[61]; // Project description string (60 chars + null terminator)
char load_descrip[61]; // Load designation string (60 chars + null terminator)
char current_time[31]; // Current timestamp string (30 chars + null terminator)
int option;         // User's menu choice
int shear;          // Boolean flag (0 for moment diagram, 1 for shear diagram)

// SVG plotting variables
int page_width = 1800;  // Width of the SVG canvas in pixels
int page_height = 1100; // Height of the SVG canvas in pixels

int cur_x;          // Current X-coordinate for drawing
int cur_y;          // Current Y-coordinate for drawing
int char_height;    // Character height for text (not directly used for SVG scaling in this port, but can be)
int char_width;     // Character width for text (not directly used for SVG scaling in this port)
int char_draw_angle;// Character draw angle (not directly used for SVG rotation in this port)
int draw_pen;       // Current drawing pen color/thickness (1:black, 2:red, 3:blue)

// Mouseover display variables for SVG output
double gX, gY;      // Global X and Y values for the current point, used in mouseover display
char gWhat[100];    // String describing what member the mouseover point belongs to
int mouseover_enabled = 0; // Flag to enable/disable mouseover functionality in SVG
int tab_cnt = 0;    // Counter for generating unique `id` attributes for tabulation spans for mouseover focus

// Function Prototypes (Declarations of all functions used in the program)
void MA(int x, int y);                     // Move Absolute: Sets the current drawing position (cur_x, cur_y)
void DA(int x, int y);                     // Draw Absolute: Draws a line from (cur_x, cur_y) to (x, y) and updates (cur_x, cur_y)
void LA(const char* text);                 // Label Absolute: Draws text at the current drawing position
void SI(int h, int w);                     // Set Ink (Size): Sets character height and width (char_height, char_width)
void DI(int a);                            // Set Direction (Angle): Sets character draw angle (char_draw_angle)
void SP(int pen);                          // Set Pen: Sets the drawing pen color/thickness (draw_pen)
void get_current_time_str(char* buffer, int buffer_size); // Gets the current time as a formatted string
void parse_up(char* s_in);                 // Converts a string to uppercase
double pw(double n, int p_val);            // Power function: calculates n to the power of p_val (n^p_val)
void readr(double* r_val);                 // Reads a real (double) value from input with error handling
void readi(int* r_val);                    // Reads an integer value from input with error handling
void reads(char* s_val);                   // Reads a string from input (up to 60 chars)
void signon();                             // Displays the program's sign-on message and prompts for project designation
void fatal_err(const char* s_err);         // Prints a fatal error message and exits the program
void initial_var();                        // Initializes global arrays and counters
void get_size();                           // Prompts for number of stories and portals, calculates matrix size (q)
void get_geometry();                       // Prompts for story heights and beam lengths, calculates total height (l3) and width (l2)
void get_member_sizes();                   // Prompts for moment of inertia (I) for beams and columns
void get_base_joint_type();                // Prompts for base joint type (fixed or pinned)
void get_macro_loads();                    // Allows the user to define and store various macro loads (groups of individual loads)
void apply_macro_loads();                  // Allows the user to apply previously defined macro loads or single loads to the frame
void setup_matrix();                       // Constructs the global stiffness matrix (d) based on frame geometry, member properties, and applied loads
void solve_matrix();                       // Solves the system of equations represented by the matrix (d) to find unknown displacements and rotations
void plot_beam_mom(int story, int beam, double m1, double v1); // Calculates and plots moment/shear diagrams for a specific beam
void output_beam();                        // Iterates through all beams, calculates their moments/shears, and calls plot_beam_mom
void plot_column_mom(int story, int column, double m1, double v1); // Calculates and plots moment/shear diagrams for a specific column
void output_column();                      // Iterates through all columns, calculates their moments/shears/axial forces, and calls plot_column_mom
void output_rotation();                    // Prints a table of joint rotations and story displacements
void reset_matrix();                       // Resets the `d` matrix and other load-related arrays to zero for a new analysis
void user_option();                        // Prompts the user for their next action (define more macros, reload frame, or exit)

// Implementations of the functions

// Graphics drawing functions (simplified for SVG output)
void MA(int x, int y) {
    cur_x = x; // Set current X coordinate
    cur_y = y; // Set current Y coordinate
}

void DA(int x, int y) {
    const char* color = "black";
    int thickness = 1;

    // Determine line color and thickness based on `draw_pen`
    switch (draw_pen) {
        case 1: color = "black"; thickness = 1; break;
        case 2: color = "red"; thickness = 1; break;
        case 3: color = "blue"; thickness = 2; break;
    }

    // If mouseover is enabled, create an invisible circle for interaction
    // and a corresponding span for tabulation data.
    if (mouseover_enabled) {
        tab_cnt++; // Increment counter for unique ID
        // SVG circle for mouse interaction. Y-coordinate is inverted for SVG (origin top-left)
        fprintf(stdout, "<circle cx=\"%d\" cy=\"%d\" r=\"5\" fill=\"transparent\" fill-opacity=\"0\" onmouseover=\"showTab(%d);\" />\n",
                x, page_height - y,
                tab_cnt);
        // Span element to display tabulation data, linked by ID to the circle
        fprintf(stdout, "<span id=\"tab%d\">%3.1f  %7.3f</span><br>\n", tab_cnt, gX, gY);
    }

    // Draw the actual line in SVG. Y-coordinates are inverted because SVG has origin at top-left.
    fprintf(stdout, "<line x1=\"%d\" y1=\"%d\" x2=\"%d\" y2=\"%d\"  style=\"stroke:%s; stroke-width:%dpx;\"/>\n",
            cur_x, page_height - cur_y, // Start point (X, inverted Y)
            x, page_height - y,         // End point (X, inverted Y)
            color,
            thickness);
    cur_x = x; // Update current X coordinate
    cur_y = y; // Update current Y coordinate
}

void LA(const char* text) {
    // Draw text in SVG at the current position. Y-coordinate is inverted.
    fprintf(stdout, "<text x=\"%d\" y=\"%d\">%s</text>\n", cur_x, page_height - cur_y, text);
}

void SI(int h, int w) { // Set character height and width
    char_height = h;
    char_width = w;
    // These values are stored but not directly applied as font-size in this SVG output,
    // as SVG text styling can be more flexible.
}

void DI(int a) { // Set character draw angle
    char_draw_angle = a;
    // This value is stored but not directly applied for text rotation in this SVG output.
}

void SP(int pen) {
    draw_pen = pen; // Set the current drawing pen (color/thickness)
}

// Utility functions
void get_current_time_str(char* buffer, int buffer_size) {
    time_t rawtime;
    struct tm *info;
    time(&rawtime); // Get current calendar time
    info = localtime(&rawtime); // Convert to local time
    strftime(buffer, buffer_size, "%Y-%m-%d %H:%M:%S", info); // Format time into string
}

void parse_up(char* s_in) {
    for (int i = 0; s_in[i] != '\0'; i++) {
        s_in[i] = toupper((unsigned char)s_in[i]); // Convert each character to uppercase
    }
}

double pw(double n, int p_val) { // Calculates n to the power of p_val
    double a = 1.0;
    for (int i = 1; i <= p_val; i++) {
        a = a * n; // Multiply n, p_val times
    }
    return a;
}

void readr(double* r_val) { // Reads a real (double) value from standard input
    int result;
    do {
        result = scanf("%lf", r_val); // Attempt to read a double
        if (result != 1) { // If reading failed
            while (getchar() != '\n'); // Clear the input buffer until newline
            printf("Invalid input. Please enter a number.\n"); // Prompt for valid input
        }
    } while (result != 1); // Repeat until a valid double is entered
}

void readi(int* r_val) { // Reads an integer value from standard input
    int result;
    do {
        result = scanf("%d", r_val); // Attempt to read an integer
        if (result != 1) { // If reading failed
            while (getchar() != '\n'); // Clear the input buffer
            printf("Invalid input. Please enter an integer.\n"); // Prompt for valid input
        }
    } while (result != 1); // Repeat until a valid integer is entered
}

void reads(char* s_val) { // Reads a string from standard input
    while (getchar() != '\n'); // Clear any leftover newline character from previous input
    if (fgets(s_val, 61, stdin) != NULL) { // Read up to 60 characters into s_val
        s_val[strcspn(s_val, "\n")] = '\0'; // Remove trailing newline character if present
    }
}

void signon() {
    printf("RIGID FRAME ANALYSIS  By Gary Argraves - 8701.15\n"); // Program title and version
    printf("\n");
    printf("Project Designation ? ");
    reads(proj_descrip); // Prompt for and read project description
    printf("%s\n", proj_descrip);
    printf("  Date: %s\n", current_time); // Display current date
    printf("\n");
}

void fatal_err(const char* s_err) {
    printf("%s\n", s_err); // Print error message
    exit(1); // Terminate program with an error code
}

void initial_var() {
    lp = 1; // Initialize load pointer
    mc = 0; // Initialize macro counter
    lc = 0; // Initialize load count
    // Initialize all global arrays to zero
    for (int i = 0; i <= 50; i++) {
        for (int j = 0; j <= 50; j++) d[i][j] = 0;
        for (int j = 0; j <= 6; j++) b[i][j] = 0;
        for (int j = 0; j <= 2; j++) c[i][j] = 0;
        for (int j = 0; j <= 4; j++) e[i][j] = 0;
        for (int j = 0; j <= 4; j++) ld[i][j] = 0;
        if (i <= 20) {
            for (int j = 0; j <= 2; j++) mp[i][j] = 0;
        }
        if (i <= 50) lm[i] = 0;
    }
}

void get_size() {
    printf("___INPUT FRAME SIZE___\n");
    printf("%25s", "How many Stories ? ");
    readi(&s); // Read number of stories
    printf("%6d\n", s);
    printf("%25s", "How many Portals ? ");
    readi(&p); // Read number of portals
    printf("%6d\n", p);
    printf("\n");
    q = (s + 1) * (p + 1) + s; // Calculate size of the stiffness matrix (degrees of freedom)
}

void get_geometry() {
    int i, j, z2_val;
    double l;
    printf("___INPUT FRAME GEOMETRY (units=feet)___\n");
    printf("   Input story heights.\n");
    l3 = 0; // Initialize total height
    for (i = 1; i <= s; i++) {
        printf("Story #%2d ? ", i);
        readr(&l); // Read story height
        printf("%8.2f\n", l);
        l3 += l; // Accumulate total height
        z2_val = (p + 1) * (s - i) + 1; // Calculate index for column data based on story
        for (j = z2_val; j <= z2_val + p; j++) c[j][2] = l; // Store column length (height)
    }
    printf("  Total height=%7.2f\n", l3);
    printf("\n");
    printf("   Input beam lengths.\n");
    l2 = 0; // Initialize total width
    for (i = 1; i <= p; i++) {
        printf("Beam #%2d ? ", i);
        readr(&l); // Read beam length
        printf("%8.2f\n", l);
        l2 += l; // Accumulate total width
        j = i;
        do {
            b[j][2] = l; // Store beam length for all beams in this portal
            j += p;
        } while (j <= i + (s - 1) * p);
    }
    printf("  Total length=%7.2f\n", l2);
    psl = 5; // Initial plot scale length
    psm = 50; // Initial plot scale moment
    // Adjust plot scales dynamically if the total length is too large for default scale
    while (l2 / psl > 9) {
        psl += 10;
        psm += 100;
    }

    printf("\n");
    printf("Page width (pixels) =%5d ? ", page_width);
    readi(&page_width); // Allow user to set SVG page width
    printf("\n");
    printf("Page height (pixels) =%5d ? ", page_height);
    readi(&page_height); // Allow user to set SVG page height

    printf("\n");
    printf("Plot Scale Length (ft/in) =%3d ? ", psl);
    readi(&psl); // Allow user to set plot scale length
    printf("\n");
    printf("Plot Scale Moment (kip-ft/inch) =%4d ? ", psm);
    readi(&psm); // Allow user to set plot scale moment
    printf("\n\n");
}

void get_member_sizes() { // Get Moment of Inertia (I) for beams and columns by story
    int i, j;
    printf("___INPUT MEMBER SIZES - Moment of Inertia (I) (units=inch^4)___\n");
    for (i = 1; i <= s; i++) { // Loop through stories
        printf("\n");
        printf("STORY #%2d\n", i);
        printf("------- --\n");
        for (j = 1; j <= p; j++) { // Loop through beams in the current story
            printf("Beam #%2d ? ", j);
            readr(&b[j + p * (s - i)][1]); // Read and store beam moment of inertia
            printf("%10.2f\n", b[j + p * (s - i)][1]);
        }
        for (j = 1; j <= p + 1; j++) { // Loop through columns in the current story
            printf("Column #%2d ? ", j);
            readr(&c[j + (s - i) * (p + 1)][1]); // Read and store column moment of inertia
            printf("%10.2f\n", c[j + (s - i) * (p + 1)][1]);
        }
    }
    printf("\n");
}

void get_base_joint_type() {
    printf("___BASE JOINT TYPE___\n");
    printf("1. Fixed       \n");
    printf("2. Pinned      \n");
    do {
        printf("%12s", "Enter 1 or 2 ? ");
        readi(&z); // Read base joint type (1 or 2)
        printf("%3d\n", z);
    } while (!(z == 1 || z == 2)); // Loop until valid input
    printf("\n");
}

void get_macro_loads() {
    int i;
    double w, j_loc, k_len; // Load parameters: w=intensity, j_loc=location, k_len=length
    double l4;
    printf("  ___________ MACRO LOAD DEFINITION ____________\n");
    printf("\n");
    // ASCII art representation of load types and locations
    printf("            +-A--+---C---+\n");
    printf("         W  !    #########--W\n");
    printf("          \\ ====================== -- Story 2\n");
    printf("        -+-#I      I      I      I\n");
    printf("         | #I      I      I      I\n");
    printf("      C--| #I      I      I      I\n");
    printf("         | #====================== -- Story 1\n");
    printf("        -+-#I      I      I      I\n");
    printf("      A--|  I  P1  I  P2  I  P3  I\n");
    printf("       __|__I______I______I______I___ Ground\n");
    printf("\n");
    printf("  S=0 :Lateral Load   S=-1 :Next Macro Definition\n");
    printf("  S>0 :Gravity Load   S=-2 :Quit Macro Definition\n");
    printf("\n");

    do {
        mc++; // Increment macro counter
        mp[mc][1] = lp; // Store starting load pointer for current macro
        do {
            printf("Macro #%2d  Load #%2d\n", mc, lp);
            printf("------- --  ------ --\n");
            printf("Enter Story # S=? ");
            readi(&i); // Read story number (0 for lateral, >0 for gravity, <0 for control)
            printf("%5d\n", i);
            // Error checking for story number
            if ((i > s) && (i != 0) && (i != -1) && (i != -2) ) fatal_err("*** Error: Incorrect stories number ***");
            if (i == 0) l4 = l3; else l4 = l2; // Determine reference length (total height for lateral, total width for gravity)
            if (i >= 0) { // If it's a valid load (lateral or gravity)
                printf("Intensity W=? ");
                readr(&w); // Read load intensity
                printf("%8.2f\n", w);
                printf("Location  A=? ");
                readr(&j_loc); // Read load start location
                printf("%8.2f\n", j_loc);
                printf("Length    C=? ");
                readr(&k_len); // Read load length
                printf("%8.2f\n", k_len);
                // Error checking for load location and length
                if (j_loc + k_len > l4) fatal_err("*** Error: A+C> Total Length of Frame ***");
                ld[lp][1] = i;      // Store story #
                ld[lp][2] = w;      // Store intensity
                ld[lp][3] = j_loc;  // Store location
                ld[lp][4] = k_len;  // Store length
                mp[mc][2] = lp;     // Store ending load pointer for current macro
                lp++; // Increment global load pointer
            } else if (mp[mc][2] == 0) { // If macro has no loads defined (e.g., user immediately enters -1 or -2)
                mc--; // Decrement macro counter to effectively "undo" the macro creation
            }
            printf("\n");
        } while (i >= 0); // Continue defining loads for the current macro until i < 0
    } while (i > -2); // Continue defining new macros until i = -2 (Quit Macro Definition)
}

// Helper function to process lateral loads (applies to columns)
void process_lateral_load(double w, double j_loc, double k_len) { // data, FEM's & Shears
    int h, o;
    double a, b2, c2, l_col;
    double T = 0, T2 = 0; // T: current cumulative height, T2: start height of current column

    for (h = 1; h <= s; h++) { // Loop through stories
        o = 1 + (s - h) * (p + 1); // Calculate index for the first column in the current story
        T += c[o][2]; // T is the height to the top of the current story's column
        T2 = T - c[o][2]; // T2 is the height to the base of the current story's column
        // Check if the load overlaps with the current column's height
        if ((j_loc < T) && (j_loc + k_len > T2)) {
            a = 0; // Relative start position of load on this column segment
            if (j_loc > T2) a = j_loc - T2; // Adjust `a` if load starts within the column
            b2 = c[o][2] - a; // Calculate effective length of load on this column segment
            if (j_loc + k_len < T) b2 = j_loc + k_len - a - T2; // Adjust `b2` if load ends within the column
            l_col = c[o][2]; // Length of the current column
            c2 = a + b2; // Relative end position of load on this column segment

            // Calculate Fixed End Moments (FEMs) and shears caused by the lateral load segment
            // and add them to the `e` array (wind loading)
            // e[h][1]: FEM at top left of column (story h)
            e[h][1] = e[h][1] - (w * fabs(1.0 / 12.0 / (l_col * l_col) * (6 * (l_col * l_col) * (c2 * c2 - a * a) - 8 * l_col * (pw(c2, 3) - pw(a, 3)) + 3 * (pw(c2, 4) - pw(a, 4)))));
            // e[h][2]: FEM at bottom right of column (story h)
            e[h][2] = e[h][2] + (w * fabs(1.0 / 12.0 / (l_col * l_col) * (4 * l_col * (pw(c2, 3) - pw(a, 3)) - 3 * (pw(c2, 4) - pw(a, 4)))));
            // e[h][3]: Shear at top left of column (story h)
            e[h][3] = e[h][3] + (w * b2 * (a + b2 / 2.0));
            // e[h][4]: Shear at bottom right of column (story h)
            e[h][4] = e[h][4] + (w * b2 * (l_col - a - b2 / 2.0));
        }
    }
}

// Helper function to process beam loads (gravity loads)
void process_beam_load(int i_story, double w, double j_loc, double k_len) { // FEM'S & shears
    int h, o;
    double a, b2, c2, l_beam;
    double T = 0, T2 = 0; // T: current cumulative width, T2: start width of current beam

    for (h = 1; h <= p; h++) { // Loop through portals (beams)
        o = h + (s - i_story) * p; // Calculate index for the current beam in the specified story
        T += b[o][2]; // T is the width to the end of the current beam
        T2 = T - b[o][2]; // T2 is the width to the start of the current beam
        // Check if the load overlaps with the current beam's length
        if ((j_loc < T) && (j_loc + k_len > T2)) {
            a = 0; // Relative start position of load on this beam segment
            if (j_loc > T2) a = j_loc - T2; // Adjust `a` if load starts within the beam
            b2 = b[o][2] - a; // Calculate effective length of load on this beam segment
            if (j_loc + k_len < T) b2 = j_loc + k_len - a - T2; // Adjust `b2` if load ends within the beam
            l_beam = b[o][2]; // Length of the current beam
            c2 = a + b2; // Relative end position of load on this beam segment

            // Calculate Fixed End Moments (FEMs) and shears caused by the beam load segment
            // and add them to the `b` array (beam data)
            // b[o][3]: FEM at left end of beam
            b[o][3] = b[o][3] - (w * fabs(1.0 / 12.0 / (l_beam * l_beam) * (6 * (l_beam * l_beam) * (c2 * c2 - a * a) - 8 * l_beam * (pw(c2, 3) - pw(a, 3)) + 3 * (pw(c2, 4) - pw(a, 4)))));
            // b[o][4]: FEM at right end of beam
            b[o][4] = b[o][4] + (w * fabs(1.0 / 12.0 / (l_beam * l_beam) * (4 * l_beam * (pw(c2, 3) - pw(a, 3)) - 3 * (pw(c2, 4) - pw(a, 4)))));
            // b[o][5]: Shear at left end of beam
            b[o][5] = b[o][5] + (w * b2 * (a + b2 / 2.0));
            // b[o][6]: Shear at right end of beam
            b[o][6] = b[o][6] + (w * b2 * (l_beam - a - b2 / 2.0));
        }
    }
}

void apply_macro_loads() {
    int u, i_load, m;
    double w_load, j_load, k_load;

    printf("___APPLY LOADS TO FRAME___\n");
    printf("\n");
    printf("  # > 0 :Macro Loading\n");
    printf("  # < 0 :Single Load\n");
    printf("  # = 0 :Quit loading\n");
    printf("\n");
    lc = 0; // Reset load count for current application
    printf("Enter Loading Designation ? ");
    reads(load_descrip); // Read and store load designation
    parse_up(load_descrip); // Convert load designation to uppercase
    printf("%s\n", load_descrip);

    do {
        printf("Enter # ? ");
        readi(&m); // Read macro number (positive for macro, negative for single load, 0 to quit)
        printf("%3d\n", m);
        if ((m > 0) && (m <= mc)) { // If it's a valid macro number
            for (u = mp[m][1]; u <= mp[m][2]; u++) { // Iterate through all loads within the macro
                lc++; // Increment load count
                lm[lc] = u; // Save the load index for later (for moment diagram generation)
                i_load = (int)ld[u][1]; // Get story number for this individual load
                w_load = ld[u][2];      // Get intensity
                j_load = ld[u][3];      // Get location
                k_load = ld[u][4];      // Get length

                if (i_load == 0) process_lateral_load(w_load, j_load, k_load); // Process as lateral load
                else process_beam_load(i_load, w_load, j_load, k_load);         // Process as beam load
            }
        } else if ((m < 0) && (abs(m) < lp)) { // If it's a valid single load number (absolute value less than total defined loads)
            u = abs(m); // Get the absolute load index
            lc++; // Increment load count
            lm[lc] = u; // Save the load index
            i_load = (int)ld[u][1];
            w_load = ld[u][2];
            j_load = ld[u][3];
            k_load = ld[u][4];
            if (i_load == 0) process_lateral_load(w_load, j_load, k_load);
            else process_beam_load(i_load, w_load, j_load, k_load);
        } else if (m != 0) { // If input is not a valid macro/load or 0
            fatal_err("*** Error: No such Macro or Load. ***");
        }
    } while (m != 0); // Continue applying loads until user enters 0
    printf("\n");
}

void setup_matrix() { // Sets up the global stiffness matrix `d`
    int z2_val, u, i, j;
    double o;

    z2_val = s + 1; // Default number of story displacement unknowns
    if (z == 1) z2_val = s; // If base joints are fixed, one less displacement unknown

    // This section populates the stiffness matrix `d` with coefficients
    // derived from the member stiffnesses and geometric properties.
    // The logic is dense as it combines terms from adjacent beams and columns
    // based on their connectivity and the degrees of freedom (rotations and displacements).
    for (u = 1; u <= z2_val; u++) { // Loop related to story displacement unknowns
        for (j = s + u * (1 + p) - p; j <= s + u * (1 + p); j++) { // Loop related to joint rotation unknowns
            if (u != 1) { // If not the first story (ground level)
                o = c[j - s - p - 1][1] / c[j - s - p - 1][2]; // Column stiffness (I/L) for column below
                if (j == s + u * (1 + p) - p) d[j][q + 1] = e[s + 2 - u][1]; // Add FEM from column below if it's the leftmost column in the story
                d[j][j] += 4 * o; // Add direct stiffness term
                d[j][j - p - 1] += 2 * o; // Add cross-stiffness term with joint below
                o = o / c[j - s - p - 1][2]; // o becomes I/L^2
                d[j][u - 1] += 6 * o; // Add stiffness term related to story displacement below
                if (u == s + 1) { // Special handling for the very top story (if fixed base)
                     goto NEXT_J; // Skips remaining calculations for this `j` and moves to the next
                }
                d[u][j] += -6 * o; // Add stiffness terms for displacement equation
                d[u][j - p - 1] += -6 * o;
                d[u][u] += 12 * o / c[j - s - p - 1][2];
                d[u][u - 1] += -12 * o / c[j - s - p - 1][2];
                d[j][u] += -6 * o;
            }
            o = c[j - s][1] / c[j - s][2]; // Column stiffness (I/L) for column above
            if (j == s + u * (1 + p) - p) d[j][q + 1] += e[s + 1 - u][2]; // Add FEM from current column if leftmost
            d[j][j] += 4 * o; // Add direct stiffness term
            if ((z != 1) || (u != s)) { // If base is not fixed OR not the very top story
                d[u][j + p + 1] += 6 * o / c[j - s][2];
                d[j][j + p + 1] += 2 * o;
            }
            o = o / c[j - s][2]; // o becomes I/L^2
            d[u][j] += 6 * o; // Add stiffness terms for displacement equation
            d[j][u] += 6 * o;
            d[u][u] += 12 * o / c[j - s][2];
            if (u != s) { // If not the highest story
                d[j][u + 1] += -6 * o;
                d[u][u + 1] += -12 * o / c[j - s][2];
            }
            if (j != s + u * (1 + p)) { // If not the rightmost column of a portal
                o = b[j - s - u + 1][1] / b[j - s - u + 1][2]; // Beam stiffness (I/L) for beam to the right
                d[j][j + 1] += 2 * o; // Add cross-stiffness term with adjacent joint on beam
                d[j][j] += 4 * o; // Add direct stiffness term
                d[j][q + 1] += b[j - s - u + 1][3]; // Add FEM from beam to the right
            }
            if (j != s + u * (1 + p) - p) { // If not the leftmost column of a portal
                o = b[j - s - u][1] / b[j - s - u][2]; // Beam stiffness (I/L) for beam to the left
                d[j][j] += 4 * o; // Add direct stiffness term
                d[j][j - 1] += 2 * o; // Add cross-stiffness term with adjacent joint on beam
                d[j][q + 1] += b[j - s - u][4]; // Add FEM from beam to the left
            }
            NEXT_J:; // Label for the `goto` statement
        } // End of inner loop (j)
    } // End of outer loop (u)

    // Add lateral loads to the right-hand side vector of the matrix
    for (i = 1; i <= s; i++) {
        u = s + 1 - i;
        d[i][q + 1] = (e[u][3] + e[u][2] + e[u][1]) / c[(i - 1) * (p + 1) + 1][2];
        if (i > 1) d[i][q + 1] += (e[u + 1][4] - e[u + 1][2] - e[u + 1][1]) / c[(i - 2) * (p + 1) + 1][2];
    }
}

void solve_matrix() { // Solves the system of linear equations using Gaussian elimination
    int i, j, l;
    double m;
    int matrix_err, go, restart;

    matrix_err = 0; // Flag for matrix errors
    do {
        restart = 0; // Flag to indicate if the outer loop needs to restart
        i = 1;
        do {
            if ((z == 1) && (i == q - p)) break; // Special case for fixed base, exit loop early
            if (d[i][i] == 0) { // If diagonal element is zero (pivot is zero)
                if (i == q) matrix_err = 1; // If it's the last row and pivot is zero, matrix error
                if (!matrix_err) {
                    j = q;
                    go = 0;
                    restart = 0;
                    do { // Search for a non-zero row to swap with current row
                        if ((j != i) && (d[j][i] != 0)) {
                            for (l = 1; l <= q + 1; l++) d[i][l] = d[i][l] + d[j][l]; // Add rows to get a non-zero pivot
                            if (j > i) go = 1; // Found a row below to swap/add from
                            else {
                                restart = 1; // Need to restart the entire process if swapped from above
                                i = q; // Set i to q to exit current loop
                            }
                        }
                        j--;
                        if ((!go) && (!restart) && (j == 0)) matrix_err = 1; // No suitable row found, matrix error
                    } while (!((j == 0) || (go) || (restart))); // Loop until row found, or end of search, or restart needed
                }
                if (matrix_err) fatal_err("*** Matrix solution Error ***"); // Terminate on matrix error
            }
            if (!restart) {
                m = d[i][i]; // Get pivot element
                for (l = 1; l <= q + 1; l++) d[i][l] = d[i][l] / m; // Normalize current row (make pivot 1)
                for (j = 1; j <= q; j++) { // Eliminate elements in other rows in the current column
                    if (j != i) {
                        m = d[j][i]; // Element to eliminate
                        for (l = 1; l <= q + 1; l++) d[j][l] = d[j][l] - m * d[i][l]; // Row operation
                    }
                }
            }
            i++;
        } while (i <= q); // Continue for all rows
    } while (restart); // Restart if necessary
}

void plot_beam_mom(int story, int beam, double m1, double v1) {
    double w, a, c2, l_beam; // load intensity, start load location, length of load, length of beam
    double sb, eb;          // start position of beam, end position of beam (cumulative width)
    double sh;              // story height (cumulative height)
    int last_lc;            // last load count index (to iterate through applied loads)
    int i, k;               // general loop counters
    int ox, oy;             // Origin X, Y for plotting the beam on the SVG canvas
    double x;               // distance along the member for plotting points

    // Inner helper function to find load segment relevant to the current beam section
    auto void find_load(double* w_out, double* a_out, double* c_out) {
        double sl, el; // start load, end load positions
        *w_out = 0; // Initialize load intensity to zero
        while ((last_lc <= lc) && (*w_out == 0)) { // Loop through applied loads until load found or all checked
            if ((int)ld[lm[last_lc]][1] == story) { // If the load is on the current story
                sl = ld[lm[last_lc]][3]; // Load start location
                el = sl + ld[lm[last_lc]][4]; // Load end location
                // Check if load overlaps with the current beam's segment
                if ((sl < sb) && (el > sb)) {
                    *w_out = ld[lm[last_lc]][2]; // Set intensity
                    *a_out = 0; // Relative start on beam
                    *c_out = el - sb; // Relative length on beam
                } else if ((sl >= sb) && (sl < eb)) {
                    *w_out = ld[lm[last_lc]][2];
                    *a_out = sl - sb;
                    *c_out = el - sl;
                }
                // Ensure the load segment does not exceed beam length
                if ((*w_out != 0) && (*a_out + *c_out > l_beam)) *c_out = l_beam - *a_out;
            }
            last_lc++; // Move to next load
        }
    }

    // Inner helper function to add moments due to distributed loads to the moment diagram array `r`
    auto void tab_moments() {
        double x_inner, mx;
        x_inner = 0.5; // Start at 0.5 ft intervals
        k = 1;
        do {
            mx = 0;
            // Calculate moment due to distributed load segment
            if ((x_inner > a) && (x_inner <= a + c2)) mx = w / 2.0 * (x_inner - a) * (x_inner - a);
            else if (x_inner > a + c2) mx = c2 * w * (x_inner - (a + c2 / 2.0));
            r[k] = r[k] + mx; // Add to existing moment at this point
            x_inner += 0.5;
            k++;
        } while (x_inner < l_beam);
        r[k] = r[k] + c2 * w * (l_beam - (a + c2 / 2.0)); // Add moment at beam end
    }

    // Inner helper function to add shears due to distributed loads to the shear diagram array `r`
    auto void tab_shears() {
        double x_inner, vx;
        x_inner = 0.5; // Start at 0.5 ft intervals
        k = 1;
        do {
            vx = 0;
            // Calculate shear due to distributed load segment
            if ((x_inner > a) && (x_inner <= a + c2)) vx = (x_inner - a) * w;
            else if (x_inner > a + c2) vx = c2 * w;
            r[k] = r[k] - vx; // Subtract from existing shear at this point
            x_inner += 0.5;
            k++;
        } while (x_inner < l_beam);
        r[k] = r[k] - c2 * w; // Subtract shear at beam end
    }

    // Calculate cumulative start and end positions of the beam within the frame's total width
    sb = 0;
    for (i = 1; i <= beam - 1; i++) sb += b[p * (s - story) + i][2];
    eb = sb + b[p * (s - story) + beam][2];
    // Calculate cumulative story height
    sh = 0;
    for (i = 1; i <= story; i++) sh += c[(p + 1) * (s - i) + 1][2];
    l_beam = eb - sb; // Length of the current beam
    last_lc = 1; // Reset load count index for `find_load`
    for (i = 0; i <= 200; i++) r[i] = 0; // Initialize moment/shear array for plotting

    // Populate `r` array with initial moment/shear values from end moments/shears (m1, v1)
    i = 0;
    x = 0;
    do {
        if (shear) r[i] = v1; // If plotting shear, `r` gets initial shear
        else r[i] = m1 - v1 * x; // If plotting moment, `r` gets moment equation (m1 - V*x)
        i++;
        x += 0.5; // Step along the beam
    } while (x < l_beam);
    if (shear) r[i] = v1;
    else r[i] = m1 - v1 * l_beam; // Value at the end of the beam

    // Add moments/shears due to distributed loads along the beam
    do {
        find_load(&w, &a, &c2); // Find any load segments on this beam
        if (w != 0) { // If a load segment is found
            if (shear) tab_shears(); // Add shear effects
            else tab_moments(); // Add moment effects
        }
    } while (last_lc <= lc); // Continue until all applied loads checked

    k = i; // Total number of points calculated for the diagram

    // Calculate plotting origin (ox, oy) for the beam based on scales and layout
    ox = (int)round((sb / psl + 0.20 * (beam) + 0.50) * PUF);
    oy = (int)round((sh / psl + 0.20 * (story)) * PUF);

    // Draw the beam itself (black line)
    mouseover_enabled = 0; // Disable mouseover for the beam outline
    SP(1); // Set pen to black
    MA(ox, oy); // Move to beam start
    DA((int)round(l_beam / psl * PUF) + ox, oy); // Draw beam line

    // Draw the zero line for the moment/shear diagram (black line)
    SP(1); // Set pen to black
    MA(ox, oy); // Move to diagram start (beam start)
    DA(ox, (int)round(r[0] / psm * PUF) + oy); // Draw line from beam start to first diagram point

    // Draw the shear or moment diagram (red for shear, blue for moment)
    mouseover_enabled = 1; // Enable mouseover for the diagram lines
    snprintf(gWhat, sizeof(gWhat), "story=%d beam=%d", story, beam); // Set member description for mouseover
    // Print header for the tabulation section
    if (shear) fprintf(stdout, "<br>Shear<br>%s<br>\n", gWhat);
    else fprintf(stdout, "<br>Moment<br>%s<br>\n", gWhat);

    if (shear) SP(2); // Set pen to red for shear
    else SP(3); // Set pen to blue for moment
    x = 0;
    for (i = 0; i <= k - 1; i++) { // Loop through calculated points
        gX = x; // Store X value for mouseover
        gY = r[i]; // Store Y value (moment/shear) for mouseover
        DA((int)round(x / psl * PUF) + ox, (int)round(r[i] / psm * PUF) + oy); // Draw segment of the diagram
        x += 0.5; // Move to next point
    }

    gX = l_beam; // Store X value at beam end
    gY = r[k]; // Store Y value at beam end
    DA((int)round(l_beam / psl * PUF) + ox, (int)round(r[k] / psm * PUF) + oy); // Draw last segment

    mouseover_enabled = 0; // Disable mouseover after drawing diagram
    SP(1); // Set pen back to black
    DA((int)round(l_beam / psl * PUF) + ox, oy); // Draw closing line for the diagram
}

void output_beam() { // Calculates and displays beam forces, then calls plotting for each beam
    int i, j, h;
    double o_val, m_val, m1_val; // stiffness, moment at left, moment at right

    if (!shear) { // Only print table headers if plotting moments
        printf("Beam%11sMoment-L%10sMoment-R%10sShears-L%10sShears-R%10s\n", "", "", "", "", "");
        printf("(#)%11s(Kip-ft)%10s(Kip-ft)%10s(Kips)%10s(Kips)%10s\n", "", "", "", "", "");
        printf("----%11s--------%10s--------%10s------%10s------%10s\n", "", "", "", "", "");
    }
    for (i = 0; i <= s - 1; i++) { // Loop through stories (from top to bottom)
        if (!shear) printf("Story%2d\n", s - i); // Print story number
        for (j = 1; j <= p; j++) { // Loop through beams in current story
            h = i * p + j; // Calculate beam index
            o_val = b[h][1] / b[h][2]; // Beam stiffness (I/L)
            // Calculate actual moments at beam ends using calculated rotations from solve_matrix
            m_val = o_val * (4 * d[s + i + h][q + 1] + 2 * d[s + i + h + 1][q + 1]) - b[h][3]; // Left end moment
            m1_val = o_val * (2 * d[s + i + h][q + 1] + 4 * d[s + i + h + 1][q + 1]) - b[h][4]; // Right end moment
            d[h][1] = (m_val + m1_val + b[h][6]) / b[h][2]; // Shear at left end (stored in d[h][1])
            d[h][2] = (b[h][5] - m_val - m1_val) / b[h][2]; // Shear at right end (stored in d[h][2])
            if (!shear) printf("%11d%10.3f%10.3f%10.3f%10.3f\n", j, m_val, m1_val, d[h][1], d[h][2]); // Print table row
            plot_beam_mom(s - i, j, m_val, d[h][1]); // Call plotting function for current beam
        }
    }
    if (!shear) printf("\n"); // Print newline after beam table
}

void plot_column_mom(int story, int column, double m1, double v1) {
    double w, a, c2, l_col; // load intensity, start load location, length of load, length of column
    double sc, ec;          // start position of column, end position of column (cumulative height)
    double bw;              // beam width (cumulative) - used for X-offset
    int last_lc;            // last load count index
    int i, k;               // general loop counters
    int ox, oy;             // Origin X, Y for plotting the column
    double x;               // distance along the member

    // Inner helper function to find load segment relevant to the current column section
    auto void find_load(double* w_out, double* a_out, double* c_out) {
        double sl, el; // start load, end load positions
        *w_out = 0;
        while ((last_lc <= lc) && (*w_out == 0)) {
            // Lateral load only applies to the leftmost column for plotting purposes
            if (((int)ld[lm[last_lc]][1] == 0) && (column == 1)) {
                sl = ld[lm[last_lc]][3];
                el = sl + ld[lm[last_lc]][4];
                if ((sl < sc) && (el > sc)) {
                    *w_out = ld[lm[last_lc]][2];
                    *a_out = 0;
                    *c_out = el - sc;
                } else if ((sl >= sc) && (sl < ec)) {
                    *w_out = ld[lm[last_lc]][2];
                    *a_out = sl - sc;
                    *c_out = el - sl;
                }
                if ((*w_out != 0) && (*a_out + *c_out > l_col)) *c_out = l_col - *a_out;
            }
            last_lc++;
        }
    }

    // Inner helper function to add moments due to distributed loads to the moment diagram array `r` for columns
    auto void tab_moments() {
        double x_inner, mx;
        x_inner = 0.5;
        k = 1;
        do {
            mx = 0;
            if ((x_inner > a) && (x_inner <= a + c2)) mx = w / 2.0 * (x_inner - a) * (x_inner - a);
            else if (x_inner > a + c2) mx = c2 * w * (x_inner - (a + c2 / 2.0));
            r[k] = r[k] - mx; // Note the sign change compared to beam moments for column orientation
            x_inner += 0.5;
            k++;
        } while (x_inner < l_col);
        r[k] = r[k] - c2 * w * (l_col - (a + c2 / 2.0));
    }

    // Inner helper function to add shears due to distributed loads to the shear diagram array `r` for columns
    auto void tab_shears() {
        double x_inner, vx;
        x_inner = 0.5;
        k = 1;
        do {
            vx = 0;
            if ((x_inner > a) && (x_inner <= a + c2)) vx = w * (x_inner - a);
            else if (x_inner > a + c2) vx = w * c2;
            r[k] = r[k] + vx; // Note the sign change compared to beam shears for column orientation
            x_inner += 0.5;
            k++;
        } while (x_inner < l_col);
        r[k] = r[k] + c2 * w;
    }

    // Calculate cumulative start and end positions of the column within the frame's total height
    sc = 0;
    for (i = 1; i <= story - 1; i++) sc += c[(p + 1) * (s - i) + 1][2];
    ec = sc + c[(p + 1) * (s - story) + 1][2];
    // Calculate cumulative beam width for X-offset of column
    bw = 0;
    for (i = 1; i <= column - 1; i++) bw += b[p * (s - story) + i][2];
    l_col = ec - sc; // Length of the current column
    last_lc = 1;
    for (i = 0; i <= 200; i++) r[i] = 0; // Initialize moment/shear array for plotting
    i = 0;
    x = 0;
    do {
        if (shear) r[i] = v1; // If plotting shear, `r` gets initial shear
        else r[i] = -(m1 + v1 * x); // If plotting moment, `r` gets moment equation (-(m1 + V*x))
        i++;
        x += 0.5;
    } while (x < l_col);
    if (shear) r[i] = v1;
    else r[i] = -(m1 + v1 * l_col);

    do {
        find_load(&w, &a, &c2);
        if (w != 0) {
            if (shear) tab_shears();
            else tab_moments();
        }
    } while (last_lc <= lc);
    k = i;

    // Calculate plotting origin (ox, oy) for the column
    ox = (int)round((bw / psl + 0.10 * (2 * column - 1) + 0.50) * PUF);
    oy = (int)round((sc / psl + 0.10 * (2 * story - 1)) * PUF);

    // Draw the column itself (black line)
    mouseover_enabled = 0;
    SP(1);
    MA(ox, oy);
    DA(ox, (int)round(l_col / psl * PUF) + oy);

    // Draw the zero line for the moment/shear diagram for the column
    SP(1);
    MA(ox, oy);
    DA((int)round(r[0] / psm * PUF) + ox, oy); // Note X and Y usage for column diagram

    // Draw the shear or moment diagram (red for shear, blue for moment)
    mouseover_enabled = 1;
    snprintf(gWhat, sizeof(gWhat), "story=%d column=%d", story, column);
    if (shear) fprintf(stdout, "<br>Shear<br>%s<br>\n", gWhat);
    else fprintf(stdout, "<br>Moment<br>%s<br>\n", gWhat);

    if (shear) SP(2);
    else SP(3);
    x = 0;
    for (i = 0; i <= k - 1; i++) {
        gX = r[i]; // X value for mouseover is the moment/shear value
        gY = x;    // Y value for mouseover is the distance along the column
        DA((int)round(r[i] / psm * PUF) + ox, (int)round(x / psl * PUF) + oy); // Draw segment of the diagram
        x += 0.5;
    }

    gX = r[k];
    gY = l_col;
    DA((int)round(r[k] / psm * PUF) + ox, (int)round(l_col / psl * PUF) + oy);

    mouseover_enabled = 0;
    SP(1);
    DA(ox, (int)round(l_col / psl * PUF) + oy);

    // Add general plot information (title, project name, date, scales)
    // This is placed here to appear once, typically at the end of the last column plot.
    if ((story == 1) && (column == p + 1)) { // Only for the last column of the last story
        MA(100, page_height - 30);
        SI(35, 35); // Set text size
        if (shear) LA("RIGID FRAME ANALYSIS SHEAR DIAGRAM");
        else LA("RIGID FRAME ANALYSIS MOMENT DIAGRAM");
        SI(27, 27); // Set text size for details
        MA(100, page_height - 60);
        LA("Project Name: ");
        fprintf(stdout, "<text x=\"%d\" y=\"%d\">%s</text>\n", cur_x + 10, page_height - cur_y, proj_descrip);
        MA(100, page_height - 90);
        LA("Date: ");
        fprintf(stdout, "<text x=\"%d\" y=\"%d\">%s</text>\n", cur_x + 10, page_height - cur_y, current_time);
        MA(100, page_height - 120);
        LA("Load Designation: ");
        fprintf(stdout, "<text x=\"%d\" y=\"%d\">%s</text>\n", cur_x + 10, page_height - cur_y, load_descrip);
        MA(100, page_height - 150);
        if (shear)
            fprintf(stdout, "<text x=\"%d\" y=\"%d\">Scales: X-%3d feet/inch  Y-%4d Kips/inch ; dots per inch = %3d</text>\n", cur_x, page_height - cur_y, psl, psm, PUF);
        else
            fprintf(stdout, "<text x=\"%d\" y=\"%d\">Scales: X-%3d feet/inch  Y-%4d Kip-Feet/inch ; dots per inch = %3d</text>\n", cur_x, page_height - cur_y, psl, psm, PUF);
        MA(0, 0); // Reset drawing position
    }
}

void output_column() { // Calculates and displays column forces, then calls plotting for each column
    int i, j, h;
    double o_val, l_val, m_val, m1_val; // stiffness, length-related value, moment at top, moment at bottom

    if (!shear) { // Only print table headers if plotting moments
        printf("Column%11sMoment%10sShears%10sAxial%10s\n", "", "", "", "");
        printf("(#)%11s(Kip-ft)%10s(Kips)%10s(Kips)%10s\n", "", "", "", "");
        printf("----%11s--------%10s------%10s------%10s\n", "", "", "", "");
    }
    for (i = 0; i <= s - 1; i++) { // Loop through stories (from top to bottom)
        if (!shear) printf("Story%2d\n", s - i); // Print story number
        for (j = 1; j <= p + 1; j++) { // Loop through columns in current story
            h = i * (p + 1) + j; // Calculate column index
            o_val = c[h][1] / c[h][2]; // Column stiffness (I/L)
            l_val = 6.0 / c[h][2] * d[i + 1][q + 1]; // Term related to story displacement
            // Calculate actual moments at column ends using calculated rotations and displacements
            m_val = o_val * (4 * d[s + h][q + 1] + 2 * d[s + h + p + 1][q + 1] + l_val); // Top moment
            m1_val = o_val * (2 * d[s + h][q + 1] + 4 * d[s + h + p + 1][q + 1] + l_val); // Bottom moment
            if (i != s - 1) { // If not the lowest story
                l_val = o_val / c[h][2] * 6 * d[i + 2][q + 1]; // Term related to displacement of story below
                m1_val -= l_val; // Adjust bottom moment
                m_val -= l_val;  // Adjust top moment
            }
            if (j <= 1) { // If it's the leftmost column
                m_val -= e[s - i][2]; // Adjust top moment for wind FEM
                m1_val -= e[s - i][1]; // Adjust bottom moment for wind FEM
            }
            o_val = 0; // Initialize axial force sum
            if (j < p + 1) o_val = d[i * p + j][1]; // Add shear from beam to the right
            if (j > 1) o_val += d[i * p + j - 1][2]; // Add shear from beam to the left
            if (i > 0) o_val += d[j][3]; // Add axial force from column above
            d[j][3] = o_val; // Store calculated axial force for current column (used for next story's column axial force)
            l_val = 0;
            if (j == 1) l_val = e[s - i][3]; // Add lateral load shear at top if leftmost column
            if (!shear) printf("%11d%10.3f%10.3f%10.3f\n", j, m_val, (m_val + m1_val - l_val) / c[h][2], o_val); // Print top shear and axial
            l_val = 0;
            if (j == 1) l_val = e[s - i][4]; // Add lateral load shear at bottom if leftmost column
            if (!shear) printf("%21s%10.3f\n", "", (m_val + m1_val + l_val) / -c[h][2]); // Print bottom shear

            plot_column_mom(s - i, j, m1_val, ((m_val + m1_val + l_val) / -c[h][2])); // Call plotting for current column
        }
    }
    if (!shear) printf("\n");
}

void output_rotation() { // Displays joint rotations and story displacements
    int i, j, l;

    printf("Joint%11sRotation%15sDisplacement%15s\n", "", "", "");
    printf("(#)%11s(Radians)%15s(Feet)%15s\n", "", "", "");
    printf("----%11s---------%15s------%15s\n", "", "", "");
    l = 0; // Counter for story displacement output (starts from 1)
    for (i = s + 1; i <= q; i++) { // Iterate through the solved degrees of freedom (joint rotations and story displacements)
        j = i - s; // Adjust index to represent original joint numbers
        // Check if it's the start of a new column group (which corresponds to a new story's displacement)
        if ((j - 1) % (p + 1) == 0) {
            printf("Story%2d\n", (int)(s - (j - 1) / (p + 1))); // Print story number
            l++;
            if (l <= s) // If it's a story displacement row
                // The factor 0.0048 might be a unit conversion (e.g., radians to degrees or pixel scale)
                printf("%11d%15.9f%15.8f\n", j, d[i][q + 1] * 0.0048, d[l][q + 1] * 0.0048);
            else // Otherwise, it's just a joint rotation
                printf("%11d%15.9f\n", j, d[i][q + 1] * 0.0048);
        } else { // For other joint rotations within a story
            printf("%11d%15.9f\n", j, d[i][q + 1] * 0.0048);
        }
    }
    printf("\n");
}


void reset_matrix() { // Resets the matrix and load-related arrays
    int i, j;
    for (i = 0; i <= q + 1; i++) {
        for (j = 0; j <= q + 1; j++) d[i][j] = 0; // Clear stiffness matrix
    }
    for (i = 0; i <= 40; i++) {
        for (j = 0; j <= 6; j++) {
            if (j >= 3 && j <= 6) b[i][j] = 0; // Clear FEMs and shears for beams
            if (i <= 10 && j >= 1 && j <= 4) e[i][j] = 0; // Clear wind loading (FEMs and shears)
        }
    }
}

void user_option() { // Presents options to the user after analysis
    printf("============================================\n");
    printf("1. Define More Macros.\n");
    printf("2. Reload the Frame.\n");
    printf("3. Exit Program.\n");
    printf("\n");
    do {
        option = 0;
        printf("Enter Option ? ");
        readi(&option); // Read user's choice
        printf("%3d\n", option);
    } while (!(option >= 1 && option <= 3)); // Validate input
    printf("\n");
}

int main() {
    // JavaScript for SVG interactivity (mouseover)
    printf("<script>\n");
    printf("var last_o=null;\n");
    printf("function showTab(id_no){ let id=`tab${id_no}`; let o=document.getElementById(id); o.scrollIntoView({behavior: \"auto\", block: \"nearest\", inline: \"nearest\"}); o.style.backgroundColor=\"lightblue\"; if (last_o != null) last_o.style.backgroundColor=\"white\"; last_o=o;}\n");
    printf("</script>\n");

    printf("RIGID FRAME ANALYSIS  By Gary Argraves - Ver. 8809.08\n"); // Program title and version
    printf("=====================================================\n");
    get_current_time_str(current_time, sizeof(current_time)); // Get current timestamp
    signon(); // Display sign-on and get project designation

    initial_var();      // Initialize global variables
    get_size();         // Get frame dimensions (stories, portals)
    get_geometry();     // Get story heights, beam lengths, and plot scales
    get_member_sizes(); // Get moment of inertia for members
    get_base_joint_type(); // Get base support type

    option = 1; // Set initial option to define macros
    do {
        if (option == 1) get_macro_loads(); // If option is 1, allow defining new macro loads
        apply_macro_loads();    // Apply loads to the frame based on user selection
        setup_matrix();         // Formulate the global stiffness matrix
        solve_matrix();         // Solve the matrix for joint rotations and displacements

        printf("\n");
        // Output for Moment Diagram (SVG plot and tabulation)
        printf("<table><tr>\n");
        printf("<td>\n");
        printf("<div style=\"height:700px; width:180px; border:1px solid black; overflow-y:scroll; margin:10 10 10 10;padding:10 10 10 10;\">\n");
        printf("<b>%s</b><br>\n", load_descrip); // Load designation for tabulation
        tab_cnt = 0; // Reset tabulation counter for new plot
        printf("</div>\n");
        printf("</td>\n");

        printf("<td>\n");
        printf("<svg width=\"%d\" height=\"%d\">\n", page_width, page_height); // Start SVG block
        // Draw a bounding box for the SVG page
        SP(1); MA(0, 0); DA(0, page_height); DA(page_width, page_height); DA(page_width, 0);  DA(0, 0);
        shear = 0; // Set flag to plot moment diagrams
        output_beam();   // Generate and plot beam moments
        output_column(); // Generate and plot column moments
        printf("</svg>\n"); // End SVG block
        printf("</td>\n");
        printf("</tr></table>\n");

        printf("\n");
        // Output for Shear Diagram (SVG plot and tabulation)
        printf("<table><tr>\n");
        printf("<td>\n");
        printf("<div style=\"height:700px; width:180px; border:1px solid black; overflow-y:scroll; margin:10 10 10 10;padding:10 10 10 10;\">\n");
        printf("<b>%s</b><br>\n", load_descrip); // Load designation for tabulation
        tab_cnt = 0; // Reset tabulation counter for new plot
        printf("</div>\n");
        printf("</td>\n");

        printf("<td>\n");
        printf("<svg width=\"%d\" height=\"%d\">\n", page_width, page_height); // Start SVG block
        // Draw a bounding box for the SVG page
        SP(1); MA(0, 0); DA(0, page_height); DA(page_width, page_height); DA(page_width, 0);  DA(0, 0);
        shear = 1; // Set flag to plot shear diagrams
        output_beam();   // Generate and plot beam shears
        output_column(); // Generate and plot column shears
        printf("</svg>\n"); // End SVG block
        printf("</td>\n");
        printf("</tr></table>\n");

        output_rotation(); // Display joint rotations and story displacements
        reset_matrix();    // Clear matrix and load data for next analysis run
        user_option();     // Prompt user for next action
    } while (option != 3); // Continue until user chooses to exit (option 3)
    printf("***** E N D   O F   R U N  *****\n");

    return 0; // Indicate successful program execution
}