/*
 *  saite.c :     Lösung der 1D-Wellengleichung
 *
 *                u  (x,t)  =  c^2 * u  (x,t)
 *                 tt                 xx
 *
 *                auf dem Intervall [0,L=1] mit einem
 *                expliziten Verfahren 2. Ordnung.
 *
 *                Anfangsbedingungen:
 *
 *                u(x,0)  = f(x)
 *
 *                u (x,0) = g(x)
 *                 t
 *
 *                Randbedingungen (feste Enden):
 *
 *                u(0,t) = u(L,t) = 0
 *
 *  u(x,t) ...... Auslenkung
 *  c ........... Phasengeschwindigkeit
 *  N + 1 ....... Anzahl der Gitterpunkte im Ortsgitter x[0],...,x[N]
 *  dx .......... Schrittweite im Ortsgitter
 *  dt .......... Zeitschritt
 *  tmax ........ Integrationsintervall [0,tmax]
 *  iout ........ Ausgabe alle "iout" Zeitschritte
 *  saite.dat ... Ausgabefile
 */


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

#define N 100

double c;
FILE *fp;



/*
 *  Anfangsbedingung  u(x,0) = f(x)
 */
double f( double x ) {

    return( exp(-100.0*(x-0.5)*(x-0.5)) );
}


/*
 *  Anfangsbedingung  u (x,0) = g(x)
 *                     t
 */
double g( double x ) {

    return( 0.0 );
}


double g1( double x ) {

    return( -200.0*c*(x-0.5)*exp(-100.0*(x-0.5)*(x-0.5)) );
}
/*
 *  Ausgabe
 */
void output( double x[], double t, double U[] ) {

    int j ;

    for ( j = 0; j <= N; j++ ) {
          fprintf(fp, "%8.5f %8.5f %13.5e\n", x[j], t, U[j]);
    }
    fprintf(fp, "\n\n");
}



int main( void ) {

    double  x[N+1], U[N+1], Uold[N+1], Unew[N+1];
    double  t, tmax, dt, dx, alpha, alpha2;
    int     j, it, iout;

    printf(" tmax [s] (0.006): ");
    scanf("%lf", &tmax);
    printf(" iout         (3): ");
    scanf("%d", &iout);

    c  = 100.0;
    dx = 1.0/N;
    dt = dx/c;

    alpha  = c*dt/dx;
    alpha2 = alpha*alpha;


   /*
    *  Ausgabefile
    */

    fp = fopen("saite.dat", "w");


   /*
    *  Räumliche Gitterpunkte
    */
     
    for ( j = 0; j <= N; j++ ) {
        x[j] = j * dx;
    }


   /*
    *  Anfangs- und Randbedingungen
    */

    t  = 0.0;
    it =
 0;

    for ( j = 1; j <= N-1; j++ ) {
        Uold[j] = f(x[j]);
    }
    Uold[0] = 0.0;
    Uold[N] = 0.0;

    output( x, t, Uold );
   /*   (1)
    *  u      (erster Zeitschritt)
    *   j
    */

    t  = dt;
    it =  1;

    for ( j = 1; j <= N-1; j++ ) {
        U[j] = (1.0 - alpha2) * Uold[j]   +
                0.5 * alpha2 * (Uold[j+1] + Uold[j-1])  +  dt * g1(x[j]);
    }
    U[0] = 0.0;
    U[N] = 0.0;




   /*   (n+1)
    *  u
    *   j
    */

    while ( t <= tmax ) {

       if ( it % iout == 0 ) {
            output( x, t, U );
       }

       for ( j = 1; j <= N-1; j++ ) {
           Unew[j] = 2.0 * (1.0 - alpha2) * U[j] +
                     alpha2 * (U[j+1] + U[j-1])  -  Uold[j];
       }
       Unew[0] = 0.0;
       Unew[N] = 0.0;

       for ( j = 0; j <= N; j++ ) {
           Uold[j] = U[j];
           U[j]    = Unew[j];
       }

       t  = t + dt;
       it = it + 1;
    }



    fclose( fp );

    return( 0 );
}













/*
 *  --------------------------------------------------------------------
 *
 *  Gnuplot
 *
 *  --------------------------------------------------------------------
 *  #
 *  # saite.gnu
 *  #
 *  set title  '1D-Wellengleichung'
 *  set xlabel 'x'
 *  set ylabel 't'
 *  set zlabel 'u(x,t)'
 *  set ticslevel 0
 *  set format y "%0.4f"
 *  set format z "%0.2f"
 *  #
 *  splot 'saite.dat' using 1:2:3 notitle with line 1
 *
 *  --------------------------------------------------------------------
 */