Program for alternate depths in rectangular channel

Principle. For a given discharge, two depths (one subcritical, one supercritical) have the same specific energy. For a rectangular channel with $q = Q/b$:

$$E = y + \frac{q^2}{2gy^2},\qquad y_c = \left(\frac{q^2}{g}\right)^{1/3},\qquad E_{min} = 1.5\,y_c$$

Given depth $y_1$, find $E$, then find the other root $y_2$ of $f(y) = y + \dfrac{q^2}{2gy^2} - E = 0$ by Newton-Raphson, with $f'(y) = 1 - \dfrac{q^2}{gy^3}$. The starting guess is taken on the opposite side of $y_c$ from $y_1$ so that the other root is found.

Algorithm

1. Start.
2. Read $Q$, $b$, $y_1$. Set $g = 9.81$.
3. Compute $q = Q/b$, $y_c = (q^2/g)^{1/3}$ and $E = y_1 + q^2/(2g\,y_1^2)$.
4. If $E < 1.5\,y_c$, print "flow not possible" and stop.
5. If $y_1 > y_c$, set $y_2 = 0.5\,y_c$; otherwise set $y_2 = 2\,y_c$.
6. Repeat: $y_{new} = y_2 - f(y_2)/f'(y_2)$ (if $y_{new} \le 0$, halve $y_2$ instead).
7. Stop when $|y_{new} - y_2| < 10^{-8}$; set $y_2 = y_{new}$.
8. Print $E$, $y_c$ and $y_2$.
9. Stop.

Program in C

#include <stdio.h>
#include <math.h>
#define G 9.81

int main(void)
{
    double Q, b, y1, q, E, yc, y2, f, df, yn;
    int i;

    printf("Enter Q (m3/s), b (m), y1 (m): ");
    scanf("%lf %lf %lf", &Q, &b, &y1);

    q  = Q / b;
    yc = cbrt(q * q / G);                      /* critical depth */
    E  = y1 + q * q / (2.0 * G * y1 * y1);     /* specific energy */

    if (E < 1.5 * yc - 1e-9) {
        printf("E is less than Emin = %.4f m : flow impossible\n", 1.5 * yc);
        return 0;
    }
    /* start on the opposite side of yc */
    y2 = (y1 > yc) ? 0.5 * yc : 2.0 * yc;

    for (i = 0; i < 100; i++) {
        f  = y2 + q * q / (2.0 * G * y2 * y2) - E;
        df = 1.0 - q * q / (G * y2 * y2 * y2);
        if (fabs(df) < 1e-12) break;
        yn = y2 - f / df;
        if (yn <= 0.0) yn = 0.5 * y2;          /* keep depth positive */
        if (fabs(yn - y2) < 1e-8) { y2 = yn; break; }
        y2 = yn;
    }
    printf("E = %.4f m, yc = %.4f m\n", E, yc);
    printf("Alternate depth y2 = %.4f m\n", y2);
    return 0;
}

Sample run: $Q = 5\ \text{m}^3/\text{s}$, $b = 3\ \text{m}$, $y_1 = 1.5\ \text{m}$ gives $E = 1.5629\ \text{m}$, $y_c = 0.6567\ \text{m}$ and alternate depth $y_2 = 0.3403\ \text{m}$.