/*
 *  ====================================================================
 *  Problem : Use finite difference approximation and iterative solution
 *            technique to compute cable profile in a simple suspension
 *            bridge.
 * 
 *  Copyright (C) 1998 by Mark Austin and David Chancogne.
 *                                                                     
 *  This software is provided "as is" without express or implied warranty.
 *  Permission is granted to use this software on any computer system,
 *  and to redistribute it freely, subject to the following restrictions:
 * 
 *  1. The authors are not responsible for the consequences of use of
 *     this software, even if they arise from defects in the software.
 *  2. The origin of this software must not be misrepresented, either
 *     by explicit claim or by omission.
 *  3. Altered versions must be plainly marked as such, and must not
 *     be misrepresented as being the original software.
 *  4. This notice is to remain intact.
 * 
 *  Written By : Mark Austin                                January 1995
 *  ====================================================================
 */

#include <stdio.h>
#include <math.h>

#define max(a, b) ((a) > (b) ? (a) : (b))

enum { LeftHandElevation  = 10 };  /* Elevation of left-hand support   */
enum { RightHandElevation = 20 };  /* Elevation of right-hand support  */
enum { BridgeWidth        = 10 };  /* Span length of suspension bridge */
enum { Nodes              = 11 };  /* No of finite difference nodes    */

int main( void ) {
float faW[ Nodes ];
float fTempW, fMaxChange;
float fdX = BridgeWidth/((float) (Nodes - 1));
int ii;

    /* [a] : Initialize cable profile */

    faW[ 0  ]        = (float) LeftHandElevation;
    faW[ Nodes - 1 ] = (float) RightHandElevation;
 
    for( ii = 1; ii <= Nodes - 2; ii = ii + 1 ) {
         faW[ ii ] = 0.0;
    }

    /* [b] : Compute profile over interior points */

    fMaxChange = max( LeftHandElevation, RightHandElevation );
    while( fMaxChange > 0.001 ) {
        fMaxChange = 0.0;

        /* [b.1] : Update profile on interior points */

        for( ii = 1; ii <= Nodes - 2; ii = ii + 1 ) {
             fTempW     = 0.5*(faW[ ii-1 ] + faW[ ii+1 ] - fdX*fdX);
             fMaxChange = max( fMaxChange, fabs( faW[ii] - fTempW ) );
             faW[ ii ]  = fTempW;
        }
    }

    /* [c] : Print x-coordinate and cable profile */

    printf("Cable profile                   \n");
    printf("=============================== \n");
    printf("      X   Numerical  Analytical \n");
    printf("  Coord   Elevation    Solution \n");
    printf("=============================== \n");

    for( ii = 0; ii < Nodes; ii = ii + 1 ) {
         printf("  %5.2f",   ii*fdX );
         printf(" %11.3f",   faW[ii] );
         printf(" %11.3f\n", 0.5*ii*fdX*ii*fdX - 4.0*ii*fdX + 10.0 );
    }
}