/************************************************************************/
/*                                                                      */
/* CS 4331: Parallel Programming                        Fall, 2012      */
/*                                                                      */
/* This program computes a facsimile of a vibrating string.  It does	*/
/* not actually simulate a vibrating string because it does not "do the	*/
/* physics" correctly.							*/
/*                                                                      */
/* The input is described below.  The output is triples:		*/
/* 	(step, point on the string, vertical displacement) 		*/
/* that can be plotted with gnuplot (and probably other 3d renders) if	*/
/* PRINTALL is defined.  If PRINTALL is not defined, then the only	*/
/* output is the list of pairs at the last step:			*/
/*	(point on the string, vertical displacement)			*/
/*                                                                      */
/************************************************************************/

#include <stdio.h>
#include <stdlib.h>

// #define PRINTALL TRUE

int main(argc,argv)
    int argc;
    char *argv[];

{   int n,		// problem size: number of points on the string
	t,		// time (seconds)
	i, step;	// counters

    double yint = 2.0,	// initial displacement of point n/2
	   yinc,	// displacement increment (for initialization)
	   disp,	// displacement (for initialization)
	   g,		// gravity
	   dt,		// delta t
	   *Sy,		// string point y displacement vector
	   *Synew,	// new displacement vector
	   *Sv,		// string point momentum (velocity) vector
	   *Svnew;	// new momentum vector

    double signy;	// the sign of Sy[i]

    // number of points: 50 to 1000 is a reasonable range
    n = atoi(argv[1]);
    // there are about t/dt steps in the simulation
    t = atoi(argv[2]);
    dt = atof(argv[3]);
    // gravitational constant: 1 to 0.1 is a reasonable range
    g = atof(argv[4]);
    printf("#  n=%d  t=%d  dt=%lf  g=%lf\n", n, t, dt, g);

    // Allocate points of the string plus two extra, fixed endpoints
    Sy = malloc(sizeof(double)*(n+2));
    Synew = malloc(sizeof(double)*(n+2));
    Sv = malloc(sizeof(double)*(n+2));
    Svnew = malloc(sizeof(double)*(n+2));

    // The string is intially displaced a distance of yint at point n/2.
    // Other points are initially displaced proportional amounts.
    // The simulation keeps Sy[0] and Sy[n+1] fixed (at 0.0 here).
    // (This is only a sample initial state.)
    // 
    //      yint=2.0 --> /\
    //           y |    /  \
    //             |   /    \
    //           a |  /      \
    //           x | /        \
    //           i |/          \
    //           s +-----+------+
    //      x axis 0    n/2     n+1

    yinc = 2*yint/(n+1);

    if (n%2 == 1)
	Sy[n/2 + 1] = yint;

    for (i=0, disp=0.0; i <= n/2; ++i, disp+=yinc)
    {
	Sy[i] = disp;
	Sy[n-i+1] = disp;
    }

    // Initial velocity of all points is 0.0.

    for (i=0; i <= n+1; ++i)
	Sv[i] = 0.0;

    // Print the initial configuration of the string.
    for (i=0; i<=n+1; ++i)
#ifdef PRINTALL
	printf("0 %d %lf\n", i, Sy[i]);
#else
	// Do not print the step number in this case
	printf("%d %lf\n", i, Sy[i]);
#endif

    // Main loop
 
    for (step=1; step<t/dt; ++step)
    {
#ifdef PRINTALL
	// This line break tells gnuplot how to break data into blocks
	printf("\n");
#endif
	for (i=1; i<=n; ++i)
	{
		signy = (Sy[i] < 0) ? 1.0 : -1.0;

		// Compute new velocity of point i
		// The first term expresses the assumption that gravity
		// increases as the square of the distance of point Sy[i]
		// from 0.0 on the y axis.  (This is the inverse of how
		// gravity actually works.)  The remaining terms express the
		// momentum contributions of point i and its two neighbors.
		// These express the elastic strength of the string, that
		// is, how much a point's neighbors' momentums affect it.
		Svnew[i] = Sv[i]/2.0 + (Sv[i-1]+Sv[i+1])/4.0
			+ (signy*Sy[i]*Sy[i]*g);
		// Compute new vertical displacement of point i
		Synew[i] = Sy[i] + dt*Sv[i];
	}
	// update velocity and position
	for (i=1; i<=n; ++i)
	{
		Sv[i] = Svnew[i];
		Sy[i] = Synew[i];
	}
#ifdef PRINTALL
	// Print triples
	for (i=0; i<=n+1; ++i)
		printf("%d %d %lf\n", step, i, Sy[i]);
#endif
    }	// End of main loop

#ifndef PRINTALL
	// Print only pairs final pairs
	// These two blank lines are significant to gnuplot -
	// They are interpreted as an EOF between two files.
	printf("\n");
	printf("\n");
	for (i=0; i<=n+1; ++i)
		printf("%d %lf\n", i, Sy[i]);
#endif

    // Almost forgot!
    free(Sv);
    free(Sy);
    free(Svnew);
    free(Synew);
}