rhcalc.js

// Event handler for the routh hurwitz button
jQuery(document).ready(function ($) {
    $('#compute-rh').on('click', function() {
        RH();
    });

    // Hide result elements
    $("#rh-matrix-card").hide();
    $("#result-card").hide();
});

// When we change the order of the system
function systemOrderChangeEvent() {
    var $sysOrderInput = $("#sys-order");
    var $sysOrderText = $("#sys-order-text")
    var order = $sysOrderInput.val();

    if (order > 0 && order <= 11) {

        switch (order % 10) {
            case 1:
                $sysOrderText.html('-st');
                break;
            case 2:
                $sysOrderText.html('-nd');
                break;
            case 3:
                $sysOrderText.html('-rd');
                break;
            default:
                $sysOrderText.html('-th');
                break;
        }

        $("#denom_input").html(regenDenomInputHTML(order));
    }
}

// Regenerate the input boxes based on system order entered
function regenDenomInputHTML(order) {
    htmlstr = '';
    for (var i = order; i > 0; i--) {
        htmlstr += '<input id="denom_coeff_' + i + '" class="coeff-input" type="number" value="0">';
        htmlstr += '<span>s</span>';
        if (order > 1) {
            htmlstr += '<sup class="sup-dense">' + i + '</sup>';
        }
        htmlstr += '+';
    }
    // Constant term
    htmlstr += '<input id="denom_coeff_0" class="coeff-input" type="number" value="1">';
    return htmlstr;
}

// Run the main routh hurwitz function
function RH() {
    // Get order
    var order = parseInt($("#sys-order").val());

    // Check that the order is correct
    if (order < 1 || order > 11) {
        console.error("Error: wrong order");
        $("#stable_result").removeClass('alert-success').addClass('alert-danger')
            .html('<strong>Error:</strong> System order must be between 1 and 11.');
        $("#result-card").show();
        return;
    }

    // Obtain all coefficients and put into an array
    var coeff = getCoefficients(order);
    if (coeff == undefined) return;

    // Check if all coefficients are zero
    if (checkAllZero(coeff)) {
        console.warn("All denominators have 0 as coefficients; system is naturally unstable");
        $("#stable_result").removeClass('alert-success').addClass('alert-warning')
            .html('<strong>Warning:</strong> All coefficients are zero. System is naturally unstable.');
        $("#result-card").show();
        return;
    }

    // Check if leading coefficient is zero
    if (coeff[0] == 0) {
        console.warn("Leading coefficient is zero");
        $("#stable_result").removeClass('alert-success').addClass('alert-warning')
            .html('<strong>Warning:</strong> Leading coefficient is zero. This may indicate an improper transfer function.');
        $("#result-card").show();
        return;
    }

    // do some more RH stuff
    var matrix = computeRH(coeff);
    var isStable = checkStability(matrix);
    outputStability(isStable, matrix);
    outputMatrix(matrix);
}

// Get coefficients from HTML and return as an array
function getCoefficients(order) {
    var coeff = [];
    for (var i = order; i >= 0; i--) {
        var input = $('#denom_coeff_' + i);
        if (input.length == 0) {
            console.error("Error: coefficient input not found for s^" + i);
            $("#stable_result").removeClass('alert-success').addClass('alert-danger')
                .html('<strong>Error:</strong> Coefficient input not found. Please refresh and try again.');
            $("#result-card").show();
            return undefined;
        }
        
        var ci = parseFloat(input.val());

        // check that every entry is valid
        if (!isNaN(ci)) {
            coeff.push(ci);
        } else {
            console.error("Error: coefficient is not a valid number for s^" + i);
            $("#stable_result").removeClass('alert-success').addClass('alert-danger')
                .html('<strong>Error:</strong> All coefficients must be valid numbers.');
            $("#result-card").show();
            return undefined;
        }
    }

    return coeff;
}

// Returns true if all of the coefficients are 0
function checkAllZero(coeff) {
    for (var i = 0, l = coeff.length; i < l; i++) {
        if (coeff[i] != 0) {
            return false;
        }
    }
    return true;
}

// returns a lower bound and upper bound for k
function computeRH(coeff) {
    var matrix = makeRHMatrix(coeff);
    var rows = matrix.length;
    var cols = matrix[0].length;

    // for each of the remaining rows
    for (var i = 2; i < rows; i++) {
        var divider = matrix[i - 1][0];

        // Special case: if entire previous row is zero
        if (isEntireRowZero(matrix[i - 1])) {
            console.warn("Entire row is zero - using derivative method");
            matrix[i - 1] = derivativeRow(matrix[i - 2]);
            divider = matrix[i - 1][0];
        }

        // Special case: if the first element is 0 but row is not entirely zero
        if (divider == 0 && !isEntireRowZero(matrix[i - 1])) {
            console.warn("First element is zero - using epsilon method");
            divider = 1e-10;
            matrix[i - 1][0] = divider;
        }

        // Final safety check to prevent division by zero
        if (divider == 0) {
            console.error("Critical error: divider is still zero after special case handling");
            divider = 1e-10;  // Use epsilon as last resort
        }

        var head1 = matrix[i - 2][0];
        var head2 = matrix[i - 1][0];

        // Calculate remaining elements in current row
        for (var j = 0; j < cols; j++) {
            var tail1 = matrix[i - 2][j + 1];
            var tail2 = matrix[i - 1][j + 1];
            if (tail1 === undefined) {
                tail1 = 0;
            }
            if (tail2 === undefined) {
                tail2 = 0;
            }

            var m = (-1 / divider) * determinant([head1, tail1, head2, tail2]);
            
            // Check for invalid results
            if (isNaN(m) || !isFinite(m)) {
                console.warn("Invalid calculation result at row " + i + ", col " + j + ": " + m);
                console.warn("Divider:", divider, "Determinant inputs:", [head1, tail1, head2, tail2]);
                m = 0; // Set to 0 if calculation fails
            }
            
            matrix[i][j] = m;
        }
    }

    console.log(matrix);
    return matrix;
}

// Generates a matrix with an appropriate size, then put the coeffcients inside
function makeRHMatrix(coeff) {
    var cols = Math.ceil(coeff.length / 2);
    var rows = coeff.length;

    // Create 2d array for the matrix full of 0s    
    var mat = new Array(rows);
    for (var r = 0; r < rows; r++) {
        mat[r] = new Array(cols);
        for (var c = 0; c < cols; c++) {
            mat[r][c] = 0;
        }
    }
    
    // Fill in the coefficients
    for (var i = 0; i < coeff.length; i++) {
        var row = i % 2;
        var col = Math.floor(i / 2);
        mat[row][col] = coeff[i];
    }

    return mat;
}

// Compute the determinant of matrix [[a, b], [c, d]] represented in 1D as [a, b, c, d]
function determinant(X) {
    if (X.length < 4) {
        console.error("Error: determinant takes 2x2 matrix or array of 4");
        return;
    }

    var det = X[0] * X[3] - X[1] * X[2];
    if (isNaN(det)) {
        det = 0;
    }

    return(det);
}

// Check if entire row is zero
function isEntireRowZero(row) {
    for (var i = 0; i < row.length; i++) {
        if (row[i] != 0) {
            return false;
        }
    }
    return true;
}

// Create derivative row when entire row is zero
function derivativeRow(prevRow) {
    var newRow = [];
    
    // The row represents coefficients of auxiliary polynomial
    // If we have [a, b, c, d] representing as^4 + bs^2 + cs^0 + d*0 (for even powers)
    // or [a, b, c, d] representing as^5 + bs^3 + cs^1 + d*0 (for odd powers)
    // The derivative would be [4a, 2b, 0, 0] or [5a, 3b, c, 0]
    
    for (var i = 0; i < prevRow.length; i++) {
        if (i < prevRow.length - 1) {
            // Power decreases by 2 for each column
            var power = (prevRow.length - 1 - i) * 2;
            newRow[i] = power * prevRow[i];
        } else {
            newRow[i] = 0; // Last element becomes 0 after differentiation
        }
    }
    
    return newRow;
}

function checkStability(matrix) {
    // Check if any element in first column is zero
    for (var i = 0; i < matrix.length; i++) {
        if (matrix[i][0] == 0) {
            return false;
        }
    }
    
    // Check if all elements in first column have the same sign
    var firstSign = matrix[0][0] > 0 ? 1 : -1;
    for (var i = 1; i < matrix.length; i++) {
        var currentSign = matrix[i][0] > 0 ? 1 : -1;
        if (currentSign != firstSign) {
            return false;
        }
    }
    
    return true;
}

function outputStability(stable, matrix) {
    var txt;
    var alertClass;
    
    if (stable) {
        txt = '<strong>System is stable!</strong><br>All elements in the first column have the same sign.';
        alertClass = 'alert-success';
    } else {
        txt = '<strong>System is unstable!</strong><br>';
        alertClass = 'alert-danger';
        
        // Count sign changes in first column
        var signChanges = 0;
        if (matrix && matrix.length > 1) {
            var prevSign = matrix[0][0] > 0 ? 1 : -1;
            for (var i = 1; i < matrix.length; i++) {
                if (matrix[i][0] != 0) {
                    var currentSign = matrix[i][0] > 0 ? 1 : -1;
                    if (currentSign != prevSign) {
                        signChanges++;
                        prevSign = currentSign;
                    }
                }
            }
            txt += 'Number of poles in right half-plane: ' + signChanges;
        }
    }
    
    $("#stable_result").removeClass('alert-success alert-danger').addClass(alertClass).html(txt);
    $("#result-card").show();
}

// Generate HTML for the matrix
function outputMatrix(matrix) {
    var mat = '';
    var rowLabels = ['s^n', 's^(n-1)', 's^(n-2)', 's^(n-3)', 's^(n-4)', 's^(n-5)', 's^(n-6)', 's^(n-7)', 's^(n-8)', 's^(n-9)', 's^(n-10)', 's^0'];
    
    for (var row = 0; row < matrix.length; row++) {
        var thisRow = matrix[row];
        mat += '<tr>';
        
        // Add row label
        mat += '<td style="background-color:#e9ecef; font-weight:bold; border-right: 2px solid #dee2e6;">';
        if (row < rowLabels.length) {
            var power = matrix.length - 1 - row;
            mat += 's<sup>' + power + '</sup>';
        } else {
            mat += 'Row ' + (row + 1);
        }
        mat += '</td>';
        
        for (var col = 0; col < thisRow.length; col++) {
            mat += '<td>';
            // Format numbers to avoid very long decimals
            var value = thisRow[col];
            if (typeof value === 'number') {
                if (Math.abs(value) < 1e-10) {
                    mat += '0';
                } else if (Math.abs(value) < 0.001 || Math.abs(value) > 1000) {
                    mat += value.toExponential(3);
                } else {
                    mat += value.toFixed(6).replace(/\.?0+$/, '');
                }
            } else {
                mat += value;
            }
            mat += '</td>';
        }
        mat += '</tr>';
    }

    $("#rh-matrix").html(mat);
    $("#rh-matrix-card").show();
}

// Test example functions
function loadExample1() {
    // s^5 + s^4 - s - 1 (unstable)
    $("#sys-order").val(5);
    systemOrderChangeEvent();
    setTimeout(function() {
        $("#denom_coeff_5").val(1);
        $("#denom_coeff_4").val(1);
        $("#denom_coeff_3").val(0);
        $("#denom_coeff_2").val(0);
        $("#denom_coeff_1").val(-1);
        $("#denom_coeff_0").val(-1);
    }, 100);
}

function loadExample2() {
    // s^6 + 3s^5 + 9s^4 + 18s^3 + 22s^2 + 12s + 12 (has zero row)
    $("#sys-order").val(6);
    systemOrderChangeEvent();
    setTimeout(function() {
        $("#denom_coeff_6").val(1);
        $("#denom_coeff_5").val(3);
        $("#denom_coeff_4").val(9);
        $("#denom_coeff_3").val(18);
        $("#denom_coeff_2").val(22);
        $("#denom_coeff_1").val(12);
        $("#denom_coeff_0").val(12);
    }, 100);
}

function loadExample3() {
    // s^3 + 2s^2 + 3s + 4 (stable)
    $("#sys-order").val(3);
    systemOrderChangeEvent();
    setTimeout(function() {
        $("#denom_coeff_3").val(1);
        $("#denom_coeff_2").val(2);
        $("#denom_coeff_1").val(3);
        $("#denom_coeff_0").val(4);
    }, 100);
}

function loadExample4() {
    // s^3 + 0s^2 + 2s + 1 (zero coefficient case)
    $("#sys-order").val(3);
    systemOrderChangeEvent();
    setTimeout(function() {
        $("#denom_coeff_3").val(1);
        $("#denom_coeff_2").val(0);
        $("#denom_coeff_1").val(2);
        $("#denom_coeff_0").val(1);
    }, 100);
}