#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <ctype.h> // For toupper
#include <time.h>  // For time, localtime, strftime

// Constants
#define PUF 96 // plotter unit factor (dots per inch)

// Global Variables
double d[51][51];   // Matrix to solve (1-based, so size +1)
double b[41][7];    // Beam Data (1-based)
double c[41][3];    // Column data (1-based)
double e[11][5];    // Wind loading (1-based)
double ld[51][5];   // loads, story, load, location, length (1-based)
int mp[21][3];      // Macro pointer of loadings (1-based)
double r[401];      // array of moments, deflections (0-based)

int s, p;           // # stories, # portals
int q;              // size of matrix to solve
int z;              // type of base joints 1 or 2
double l2, l3;      // total height, total width feet

int lp;             // load pointer, points to next row in load matrix
int mc;             // macro counter, multible loads per macro
int lc;             // load count of loads applied
int lm[51];         // save load#'s of all loads applied (1-based)
int psl;            // plot scale member length
int psm;            // plot scale moment kip-ft/inch
char proj_descrip[61]; // 60 chars + null terminator
char load_descrip[61]; // 60 chars + null terminator
char current_time[31]; // 30 chars + null terminator
int option;
int shear;          // boolean, using int

int page_width = 1800;
int page_height = 1100;

int cur_x;
int cur_y;
int char_height;
int char_width;
int char_draw_angle;
int draw_pen;
// plot_str and tabulation_str will be directly printed to stdout for simplicity in C
double gX, gY;      // save for mouseover display
char gWhat[100];    // show what member the mouseover display is for
int mouseover_enabled = 0; // boolean, using int
int tab_cnt = 0;    // used to make <span id> to focus the tabulation elements on mouseover

// Function Prototypes
void MA(int x, int y);
void DA(int x, int y);
void LA(const char* text);
void SI(int h, int w);
void DI(int a);
void SP(int pen);
void get_current_time_str(char* buffer, int buffer_size);
void parse_up(char* s_in);
double pw(double n, int p_val);
void readr(double* r_val);
void readi(int* r_val);
void reads(char* s_val);
void signon();
void fatal_err(const char* s_err);
void initial_var();
void get_size();
void get_geometry();
void get_member_sizes();
void get_base_joint_type();
void get_macro_loads();
void apply_macro_loads();
void setup_matrix();
void solve_matrix();
void plot_beam_mom(int story, int beam, double m1, double v1);
void output_beam();
void plot_column_mom(int story, int column, double m1, double v1);
void output_column();
void output_rotation();
void reset_matrix();
void user_option();

// Implementations

void MA(int x, int y) {
    cur_x = x;
    cur_y = y;
}

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

    switch (draw_pen) {
        case 1: color = "black"; thickness = 1; break;
        case 2: color = "red"; thickness = 1; break;
        case 3: color = "blue"; thickness = 2; break;
    }

    // use an invisible circle for easier mouse snap ofr coordinate display
    if (mouseover_enabled) {
        tab_cnt++;
        fprintf(stdout, "<circle cx=\"%d\" cy=\"%d\" r=\"5\" fill=\"transparent\" fill-opacity=\"0\" onmouseover=\"showTab(%d);\" />\n",
                x, page_height - y,
                tab_cnt);
        fprintf(stdout, "<span id=\"tab%d\">%3.1f  %7.3f</span><br>\n", tab_cnt, gX, gY);
    }

    fprintf(stdout, "<line x1=\"%d\" y1=\"%d\" x2=\"%d\" y2=\"%d\"  style=\"stroke:%s; stroke-width:%dpx;\"/>\n",
            cur_x, page_height - cur_y, // set Y origin at bottom
            x, page_height - y,
            color,
            thickness);
    cur_x = x;
    cur_y = y;
}

void LA(const char* text) {
    fprintf(stdout, "<text x=\"%d\" y=\"%d\">%s</text>\n", cur_x, page_height - cur_y, text);
}

void SI(int h, int w) { // set char height width
    char_height = h;
    char_width = w;
    // For SVG text, font-size can be set globally or per text element.
    // In this port, we are not directly using char_height/width for scaling,
    // but the original code might have implicitly used it for cursor positioning or plotting.
    // For simplicity, we just store it.
}

void DI(int a) { // set char draw angle
    char_draw_angle = a;
    // Not directly used in SVG output for text rotation in this simple port,
    // but can be added if needed for more complex rendering.
}

void SP(int pen) {
    draw_pen = pen;
}

void get_current_time_str(char* buffer, int buffer_size) {
    time_t rawtime;
    struct tm *info;
    time(&rawtime);
    info = localtime(&rawtime);
    strftime(buffer, buffer_size, "%Y-%m-%d %H:%M:%S", info); // Example format
}

void parse_up(char* s_in) {
    for (int i = 0; s_in[i] != '\0'; i++) {
        s_in[i] = toupper((unsigned char)s_in[i]);
    }
}

double pw(double n, int p_val) {
    double a = 1.0;
    for (int i = 1; i <= p_val; i++) {
        a = a * n;
    }
    return a;
}

void readr(double* r_val) {
    int result;
    do {
        result = scanf("%lf", r_val);
        if (result != 1) {
            // Clear input buffer on error
            while (getchar() != '\n');
            printf("Invalid input. Please enter a number.\n");
        }
    } while (result != 1);
}

void readi(int* r_val) {
    int result;
    do {
        result = scanf("%d", r_val);
        if (result != 1) {
            // Clear input buffer on error
            while (getchar() != '\n');
            printf("Invalid input. Please enter an integer.\n");
        }
    } while (result != 1);
}

void reads(char* s_val) {
    // Read up to 60 characters, including spaces, until newline
    // clear stdin buffer before reading string
    while (getchar() != '\n');
    if (fgets(s_val, 61, stdin) != NULL) {
        // Remove trailing newline character if present
        s_val[strcspn(s_val, "\n")] = '\0';
    }
}

void signon() {
    printf("RIGID FRAME ANALYSIS  By Gary Argraves - 8701.15\n");
    printf("\n");
    printf("Project Designation ? ");
    reads(proj_descrip);
    printf("%s\n", proj_descrip);
    printf("  Date: %s\n", current_time);
    printf("\n");
}

void fatal_err(const char* s_err) {
    printf("%s\n", s_err);
    exit(1);
}

void initial_var() {
    lp = 1;
    mc = 0;
    lc = 0;
    for (int i = 0; i <= 50; i++) { // Adjusted for 0-based indexing if necessary,
                                    // or stick to Pascal's 1-based logic if arrays match
        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; // e is 1..10,1..4 in Pascal, so 0-10,0-4 for C
        for (int j = 0; j <= 4; j++) ld[i][j] = 0;
        if (i <= 20) { // mp is 1..20,1..2
            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);
    printf("%6d\n", s);
    printf("%25s", "How many Portals ? ");
    readi(&p);
    printf("%6d\n", p);
    printf("\n");
    q = (s + 1) * (p + 1) + s;
}

void get_geometry() {
    int i, j, z2_val;
    double l;
    printf("___INPUT FRAME GEOMETRY (units=feet)___\n");
    printf("   Input story heights.\n");
    l3 = 0;
    for (i = 1; i <= s; i++) {
        printf("Story #%2d ? ", i);
        readr(&l);
        printf("%8.2f\n", l);
        l3 += l;
        z2_val = (p + 1) * (s - i) + 1;
        for (j = z2_val; j <= z2_val + p; j++) c[j][2] = l;
    }
    printf("  Total height=%7.2f\n", l3);
    printf("\n");
    printf("   Input beam lengths.\n");
    l2 = 0;
    for (i = 1; i <= p; i++) {
        printf("Beam #%2d ? ", i);
        readr(&l);
        printf("%8.2f\n", l);
        l2 += l;
        j = i;
        do {
            b[j][2] = l;
            j += p;
        } while (j <= i + (s - 1) * p);
    }
    printf("  Total length=%7.2f\n", l2);
    psl = 5;
    psm = 50;
    while (l2 / psl > 9) { // (*** modify scale as required ***)
        psl += 10;
        psm += 100;
    }

    printf("\n");
    printf("Page width (pixels) =%5d ? ", page_width);
    readi(&page_width);
    printf("\n");
    printf("Page height (pixels) =%5d ? ", page_height);
    readi(&page_height);

    printf("\n");
    printf("Plot Scale Length (ft/in) =%3d ? ", psl);
    readi(&psl);
    printf("\n");
    printf("Plot Scale Moment (kip-ft/inch) =%4d ? ", psm);
    readi(&psm);
    printf("\n\n");
}

void get_member_sizes() { // by story
    int i, j;
    printf("___INPUT MEMBER SIZES - Moment of Inertia (I) (units=inch^4)___\n");
    for (i = 1; i <= s; i++) {
        printf("\n");
        printf("STORY #%2d\n", i);
        printf("------- --\n");
        for (j = 1; j <= p; j++) {
            printf("Beam #%2d ? ", j);
            readr(&b[j + p * (s - i)][1]);
            printf("%10.2f\n", b[j + p * (s - i)][1]);
        }
        for (j = 1; j <= p + 1; j++) {
            printf("Column #%2d ? ", j);
            readr(&c[j + (s - i) * (p + 1)][1]);
            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);
        printf("%3d\n", z);
    } while (!(z == 1 || z == 2));
    printf("\n");
}

void get_macro_loads() {
    int i;
    double w, j_loc, k_len; // load parameters
    double l4;
    printf("  ___________ MACRO LOAD DEFINITION ____________\n");
    printf("\n");
    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");
    // if (lp < 1) lp = 1; // This line was commented out in Pascal

    do {
        mc++;
        mp[mc][1] = lp;
        do {
            printf("Macro #%2d  Load #%2d\n", mc, lp);
            printf("------- --  ------ --\n");
            printf("Enter Story # S=? ");
            readi(&i);
            printf("%5d\n", i);
            if ((i > s) && (i != 0) && (i != -1) && (i != -2) ) fatal_err("*** Error: Incorrect stories number ***"); // Added checks for -1, -2 for safety
            if (i == 0) l4 = l3; else l4 = l2; // lateral load is over total height, gravity over total width
            if (i >= 0) {
                printf("Intensity W=? ");
                readr(&w);
                printf("%8.2f\n", w);
                printf("Location  A=? ");
                readr(&j_loc);
                printf("%8.2f\n", j_loc);
                printf("Length    C=? ");
                readr(&k_len);
                printf("%8.2f\n", k_len);
                if (j_loc + k_len > l4) fatal_err("*** Error: A+C> Total Length of Frame ***");
                ld[lp][1] = i;
                ld[lp][2] = w;
                ld[lp][3] = j_loc;
                ld[lp][4] = k_len;
                mp[mc][2] = lp;
                lp++; // bump load pointer
            } else if (mp[mc][2] == 0) {
                mc--;
            }
            printf("\n");
        } while (i >= 0); // until i < 0
    } while (i > -2); // until i < -1
}

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; // Cumulative height and start height of current column

    for (h = 1; h <= s; h++) {
        o = 1 + (s - h) * (p + 1); // o is the column number (1-based for c array)
        T += c[o][2]; // T is the height to top of column
        T2 = T - c[o][2]; // T2 is the height to column base
        if ((j_loc < T) && (j_loc + k_len > T2)) { // Is the load in this section ?
            a = 0; // Yes, A is the relative start position
            if (j_loc > T2) a = j_loc - T2; // Load starts in from segment beginning
            b2 = c[o][2] - a; // B is the load length in this segment
            if (j_loc + k_len < T) b2 = j_loc + k_len - a - T2; // Load does not go to end of segment
            l_col = c[o][2]; // L is the column length
            c2 = a + b2; // End position of load segment relative to column start

            // Note: Pascal's sqr(L) is L*L, pw(C2,3) is C2^3. Use pow for clarity or direct multiplication.
            // Ensure ABS is handled for numerical stability as in Pascal.
            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] = 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] = e[h][3] + (w * b2 * (a + b2 / 2.0));
            e[h][4] = e[h][4] + (w * b2 * (l_col - a - b2 / 2.0));
        }
    }
}

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; // Cumulative width and start width of current beam

    for (h = 1; h <= p; h++) {
        o = h + (s - i_story) * p; // o is the beam number (1-based for b array)
        T += b[o][2]; // T is the width to end of beam
        T2 = T - b[o][2]; // T2 is the width to beam start
        if ((j_loc < T) && (j_loc + k_len > T2)) { // Is the load in this section ?
            a = 0;
            if (j_loc > T2) a = j_loc - T2;
            b2 = b[o][2] - a;
            if (j_loc + k_len < T) b2 = j_loc + k_len - a - T2;
            l_beam = b[o][2];
            c2 = a + b2;

            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] = 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] = b[o][5] + (w * b2 * (a + b2 / 2.0));
            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;
    printf("Enter Loading Designation ? ");
    reads(load_descrip);
    parse_up(load_descrip);
    printf("%s\n", load_descrip);

    do {
        printf("Enter # ? ");
        readi(&m);
        printf("%3d\n", m);
        if ((m > 0) && (m <= mc)) {
            for (u = mp[m][1]; u <= mp[m][2]; u++) {
                lc++;
                lm[lc] = u; // save all loads applied for moment dia
                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) && (abs(m) < lp)) { // lp points to next available row, so must be less than lp
            u = abs(m);
            lc++;
            lm[lc] = u;
            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) {
            fatal_err("*** Error: No such Macro or Load. ***");
        }
    } while (m != 0);
    printf("\n");
}

void setup_matrix() {
    int z2_val, u, i, j;
    double o;

    z2_val = s + 1;
    if (z == 1) z2_val = s;

    for (u = 1; u <= z2_val; u++) {
        for (j = s + u * (1 + p) - p; j <= s + u * (1 + p); j++) {
            if (u != 1) {
                o = c[j - s - p - 1][1] / c[j - s - p - 1][2];
                if (j == s + u * (1 + p) - p) d[j][q + 1] = e[s + 2 - u][1];
                d[j][j] += 4 * o;
                d[j][j - p - 1] += 2 * o;
                o = o / c[j - s - p - 1][2];
                d[j][u - 1] += 6 * o;
                if (u == s + 1) {
                     // GOTO NEXT_J logic implemented using a continue for inner loop
                     // Since this is the innermost loop, it effectively skips to the next j
                     goto NEXT_J; // Used goto as refactoring to avoid it might be complex and lead to less readable code
                }
                d[u][j] += -6 * o;
                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];
            if (j == s + u * (1 + p) - p) d[j][q + 1] += e[s + 1 - u][2];
            d[j][j] += 4 * o;
            if ((z != 1) || (u != s)) {
                d[u][j + p + 1] += 6 * o / c[j - s][2];
                d[j][j + p + 1] += 2 * o;
            }
            o = o / c[j - s][2];
            d[u][j] += 6 * o;
            d[j][u] += 6 * o;
            d[u][u] += 12 * o / c[j - s][2];
            if (u != s) {
                d[j][u + 1] += -6 * o;
                d[u][u + 1] += -12 * o / c[j - s][2];
            }
            if (j != s + u * (1 + p)) {
                o = b[j - s - u + 1][1] / b[j - s - u + 1][2];
                d[j][j + 1] += 2 * o;
                d[j][j] += 4 * o;
                d[j][q + 1] += b[j - s - u + 1][3];
            }
            if (j != s + u * (1 + p) - p) {
                o = b[j - s - u][1] / b[j - s - u][2];
                d[j][j] += 4 * o;
                d[j][j - 1] += 2 * o;
                d[j][q + 1] += b[j - s - u][4];
            }
            NEXT_J:; // Label for goto
        } // for j
    } // for u

    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() {
    int i, j, l;
    double m;
    int matrix_err, go, restart;

    matrix_err = 0;
    do {
        restart = 0;
        i = 1;
        do {
            if ((z == 1) && (i == q - p)) break; // Equivalent to exit for the loop
            if (d[i][i] == 0) {
                if (i == q) matrix_err = 1;
                if (!matrix_err) {
                    j = q;
                    go = 0;
                    restart = 0;
                    do {
                        if ((j != i) && (d[j][i] != 0)) {
                            for (l = 1; l <= q + 1; l++) d[i][l] = d[i][l] + d[j][l];
                            if (j > i) go = 1;
                            else {
                                restart = 1;
                                i = q; // To exit the outer do-while loop
                            }
                        }
                        j--;
                        if ((!go) && (!restart) && (j == 0)) matrix_err = 1;
                    } while (!((j == 0) || (go) || (restart)));
                }
                if (matrix_err) fatal_err("*** Matrix solution Error ***");
            }
            if (!restart) {
                m = d[i][i];
                for (l = 1; l <= q + 1; l++) d[i][l] = d[i][l] / m;
                for (j = 1; j <= q; j++) {
                    if (j != i) {
                        m = d[j][i];
                        for (l = 1; l <= q + 1; l++) d[j][l] = d[j][l] - m * d[i][l];
                    }
                }
            }
            i++;
        } while (i <= q);
    } while (restart);
}

void plot_beam_mom(int story, int beam, double m1, double v1) {
    double w, a, c2, l_beam; // load, start load, length load, length beam
    double sb, eb;          // start beam, end beam length
    double sh;              // story height
    int last_lc;            // last load count
    int i, k;               // general
    int ox, oy;
    double x;               // distance along member

    // Nested function in Pascal, implemented as local helpers in C if needed, or by passing args.
    // For simplicity, these helpers will be inlined or recreated with necessary parameters.

    auto void find_load(double* w_out, double* a_out, double* c_out) {
        double sl, el; // start load, end load
        *w_out = 0;
        while ((last_lc <= lc) && (*w_out == 0)) {
            if ((int)ld[lm[last_lc]][1] == story) {
                sl = ld[lm[last_lc]][3];
                el = sl + ld[lm[last_lc]][4];
                if ((sl < sb) && (el > sb)) {
                    *w_out = ld[lm[last_lc]][2];
                    *a_out = 0;
                    *c_out = el - sb;
                } else if ((sl >= sb) && (sl < eb)) {
                    *w_out = ld[lm[last_lc]][2];
                    *a_out = sl - sb;
                    *c_out = el - sl;
                }
                if ((*w_out != 0) && (*a_out + *c_out > l_beam)) *c_out = l_beam - *a_out;
            }
            last_lc++;
        }
    }

    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;
            x_inner += 0.5;
            k++;
        } while (x_inner < l_beam);
        r[k] = r[k] + c2 * w * (l_beam - (a + c2 / 2.0));
    }

    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 = (x_inner - a) * w;
            else if (x_inner > a + c2) vx = c2 * w;
            r[k] = r[k] - vx;
            x_inner += 0.5;
            k++;
        } while (x_inner < l_beam);
        r[k] = r[k] - c2 * w;
    }

    sb = 0;
    for (i = 1; i <= beam - 1; i++) sb += b[p * (s - story) + i][2];
    eb = sb + b[p * (s - story) + beam][2];
    sh = 0;
    for (i = 1; i <= story; i++) sh += c[(p + 1) * (s - i) + 1][2];
    l_beam = eb - sb;
    last_lc = 1;
    for (i = 0; i <= 200; i++) r[i] = 0;

    i = 0;
    x = 0;
    do {
        if (shear) r[i] = v1;
        else r[i] = m1 - v1 * x;
        i++;
        x += 0.5;
    } while (x < l_beam);
    if (shear) r[i] = v1;
    else r[i] = m1 - v1 * l_beam;

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

    ox = (int)round((sb / psl + 0.20 * (beam) + 0.50) * PUF);
    oy = (int)round((sh / psl + 0.20 * (story)) * PUF);

    // draw the beam, l is length, 1=black
    mouseover_enabled = 0; // disable onmouseover
    SP(1);
    MA(ox, oy);
    DA((int)round(l_beam / psl * PUF) + ox, oy);

    // value at 0 feet is displayed at the Y-value point, not at y=0
    SP(1);
    MA(ox, oy);
    DA(ox, (int)round(r[0] / psm * PUF) + oy);

    // draw the shear or the moment; whichever is selected, 2=red, 3=blue
    mouseover_enabled = 1; // enable onmouseover
    snprintf(gWhat, sizeof(gWhat), "story=%d beam=%d", story, beam);
    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 = x;
        gY = r[i]; // save for mouseover display
        DA((int)round(x / psl * PUF) + ox, (int)round(r[i] / psm * PUF) + oy);
        x += 0.5;
    }

    gX = l_beam;
    gY = r[k]; // save for mouseover display
    DA((int)round(l_beam / psl * PUF) + ox, (int)round(r[k] / psm * PUF) + oy);

    mouseover_enabled = 0; // disable onmouseover
    SP(1);
    DA((int)round(l_beam / psl * PUF) + ox, oy);
}

void output_beam() {
    int i, j, h;
    double o_val, m_val, m1_val;
    if (!shear) {
        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++) {
        if (!shear) printf("Story%2d\n", s - i);
        for (j = 1; j <= p; j++) {
            h = i * p + j;
            o_val = b[h][1] / b[h][2];
            m_val = o_val * (4 * d[s + i + h][q + 1] + 2 * d[s + i + h + 1][q + 1]) - b[h][3];
            m1_val = o_val * (2 * d[s + i + h][q + 1] + 4 * d[s + i + h + 1][q + 1]) - b[h][4];
            d[h][1] = (m_val + m1_val + b[h][6]) / b[h][2]; // Shear-L
            d[h][2] = (b[h][5] - m_val - m1_val) / b[h][2]; // Shear-R
            if (!shear) printf("%11d%10.3f%10.3f%10.3f%10.3f\n", j, m_val, m1_val, d[h][1], d[h][2]);
            plot_beam_mom(s - i, j, m_val, d[h][1]);
        }
    }
    if (!shear) printf("\n");
}

void plot_column_mom(int story, int column, double m1, double v1) {
    double w, a, c2, l_col; // load, start load, length load, length column
    double sc, ec;          // start colum, end column length
    double bw;              // beam width (cumulative)
    int last_lc;            // last load count
    int i, k;               // general
    int ox, oy;
    double x;

    auto void find_load(double* w_out, double* a_out, double* c_out) {
        double sl, el; // start load, end load
        *w_out = 0;
        while ((last_lc <= lc) && (*w_out == 0)) {
            if (((int)ld[lm[last_lc]][1] == 0) && (column == 1)) { // Lateral load only applies to leftmost column for display
                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++;
        }
    }

    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 from beam moments
            x_inner += 0.5;
            k++;
        } while (x_inner < l_col);
        r[k] = r[k] - c2 * w * (l_col - (a + c2 / 2.0));
    }

    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 from beam shears
            x_inner += 0.5;
            k++;
        } while (x_inner < l_col);
        r[k] = r[k] + c2 * w;
    }

    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];
    bw = 0;
    for (i = 1; i <= column - 1; i++) bw += b[p * (s - story) + i][2];
    l_col = ec - sc;
    last_lc = 1;
    for (i = 0; i <= 200; i++) r[i] = 0;
    i = 0;
    x = 0;
    do {
        if (shear) r[i] = v1;
        else r[i] = -(m1 + v1 * x); // Note the sign change
        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;

    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, 1=black
    mouseover_enabled = 0; // disable onmouseover
    SP(1);
    MA(ox, oy);
    DA(ox, (int)round(l_col / psl * PUF) + oy);

    // value at 0 feet is displayed at the Y-value point, not at y=0
    SP(1);
    MA(ox, oy);
    DA((int)round(r[0] / psm * PUF) + ox, oy);

    // draw the shear or the moment; whichever is selected, 2=red, 3=blue
    mouseover_enabled = 1; // enable onmouseover
    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 and y roles reversed for column plot
        gY = x;    // when plotting with origin at bottom-left and rotating
        DA((int)round(r[i] / psm * PUF) + ox, (int)round(x / psl * PUF) + oy);
        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; // disable onmouseover
    SP(1);
    DA(ox, (int)round(l_col / psl * PUF) + oy);

    if ((story == 1) && (column == p + 1)) {
        MA(100, page_height - 30);
        SI(35, 35); // set character height
        if (shear) LA("RIGID FRAME ANALYSIS SHEAR DIAGRAM");
        else LA("RIGID FRAME ANALYSIS MOMENT DIAGRAM");
        SI(27, 27);
        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);
    }
}

void output_column() {
    int i, j, h;
    double o_val, l_val, m_val, m1_val;
    if (!shear) {
        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++) {
        if (!shear) printf("Story%2d\n", s - i);
        for (j = 1; j <= p + 1; j++) {
            h = i * (p + 1) + j;
            o_val = c[h][1] / c[h][2];
            l_val = 6.0 / c[h][2] * d[i + 1][q + 1];
            m_val = o_val * (4 * d[s + h][q + 1] + 2 * d[s + h + p + 1][q + 1] + l_val);
            m1_val = o_val * (2 * d[s + h][q + 1] + 4 * d[s + h + p + 1][q + 1] + l_val);
            if (i != s - 1) {
                l_val = o_val / c[h][2] * 6 * d[i + 2][q + 1];
                m1_val -= l_val;
                m_val -= l_val;
            }
            if (j <= 1) {
                m_val -= e[s - i][2];
                m1_val -= e[s - i][1];
            }
            o_val = 0;
            if (j < p + 1) o_val = d[i * p + j][1];
            if (j > 1) o_val += d[i * p + j - 1][2];
            if (i > 0) o_val += d[j][3];
            d[j][3] = o_val; // Axial force stored in d[j][3] for column
            l_val = 0;
            if (j == 1) l_val = e[s - i][3];
            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);
            l_val = 0;
            if (j == 1) l_val = e[s - i][4];
            if (!shear) printf("%21s%10.3f\n", "", (m_val + m1_val + l_val) / -c[h][2]); // Shear at bottom

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

void output_rotation() {
    int i, j, l;

    printf("Joint%11sRotation%15sDisplacement%15s\n", "", "", "");
    printf("(#)%11s(Radians)%15s(Feet)%15s\n", "", "", "");
    printf("----%11s---------%15s------%15s\n", "", "", "");
    l = 0; // This `l` variable seems to be used as a counter for story displacements
    for (i = s + 1; i <= q; i++) {
        j = i - s;
        if ((j - 1) % (p + 1) == 0) { // Check if it's the start of a new column group (new story)
            printf("Story%2d\n", (int)(s - (j - 1) / (p + 1)));
            l++;
            if (l <= s) // Displacement for the story (stored in d[l][q+1])
                printf("%11d%15.9f%15.8f\n", j, d[i][q + 1] * 0.0048, d[l][q + 1] * 0.0048);
            else
                printf("%11d%15.9f\n", j, d[i][q + 1] * 0.0048);
        } else {
            printf("%11d%15.9f\n", j, d[i][q + 1] * 0.0048);
        }
    }
    printf("\n");
}


void reset_matrix() {
    int i, j;
    for (i = 0; i <= q + 1; i++) { // Adjusted to handle array sizes correctly
        for (j = 0; j <= q + 1; j++) d[i][j] = 0;
    }
    for (i = 0; i <= 40; i++) { // Adjusted for array size
        for (j = 0; j <= 6; j++) { // b[i][j] is 1..40, 1..6
            if (j >= 3 && j <= 6) b[i][j] = 0;
            if (i <= 10 && j >= 1 && j <= 4) e[i][j] = 0; // e[i][j-2] becomes e[i][j] if j from 1..4
        }
    }
}

void user_option() {
    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);
        printf("%3d\n", option);
    } while (!(option >= 1 && option <= 3));
    printf("\n");
}

int main() {
    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");
    printf("=====================================================\n");
    get_current_time_str(current_time, sizeof(current_time));
    signon();

    initial_var();
    get_size();
    get_geometry();
    get_member_sizes();
    get_base_joint_type();
    option = 1;
    do {
        if (option == 1) get_macro_loads();
        apply_macro_loads();
        setup_matrix();
        solve_matrix();

        printf("\n");
        // For tabulation_str and plot_str, we print directly to stdout as they build up.
        // This is a simplification from dynamically building large strings.
        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);
        // Reset tab_cnt for each new plot section
        tab_cnt = 0;
        printf("</div>\n");
        printf("</td>\n");

        printf("<td>\n");
        printf("<svg width=\"%d\" height=\"%d\">\n", page_width, page_height);
        // box the page
        SP(1);
        MA(0, 0);
        DA(0, page_height);
        DA(page_width, page_height);
        DA(page_width, 0);
        DA(0, 0);
        shear = 0; // false
        output_beam();
        output_column();
        printf("</svg>\n");
        printf("</td>\n");
        printf("</tr></table>\n");


        printf("\n");
        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);
        // Reset tab_cnt for each new plot section
        tab_cnt = 0;
        printf("</div>\n");
        printf("</td>\n");

        printf("<td>\n");
        printf("<svg width=\"%d\" height=\"%d\">\n", page_width, page_height);
        // box the page
        SP(1);
        MA(0, 0);
        DA(0, page_height);
        DA(page_width, page_height);
        DA(page_width, 0);
        DA(0, 0);
        shear = 1; // true
        output_beam();
        output_column();
        printf("</svg>\n");
        printf("</td>\n");
        printf("</tr></table>\n");


        output_rotation();
        reset_matrix();
        user_option();
    } while (option != 3);
    printf("***** E N D   O F   R U N  *****\n");

    return 0;
}