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

#include <math.h>
#include <string.h>

#include "io.h"

/*
 * main.c : acoustic equations
 */
int
main (void)
{

    int n = 500;
    double xl = 0.0;
    double xr = 1.0;
    double dx = (xr - xl) / n;

    double dt = 5.0e-4;

    double* x = (double*) calloc (n, sizeof (double));
    double* u = (double*) calloc (4 * n, sizeof (double));
    double* b = (double*) calloc (n, sizeof (double));

    double* u12 = (double*) calloc (4 * (n + 1), sizeof (double));
    double* f12 = (double*) calloc (4 * (n + 1), sizeof (double));

    double* unew = (double*) calloc (4 * n, sizeof (double));

    // initial conditions

    for (int i = 0; i < n; i ++) {
        x[i] = 0.5 * dx + dx * i;
    }

    iocoord("coord.dat", x, n);

    double hbar    = 0.2;
    double xbar    = 0.075;
    double x0      = 0.4;

    for (int i = 0; i < n; i ++) {
        if (fabs (x[i] - x0) < xbar) {
            u[4*i] = hbar * sqrt (1.0 - pow ((x[i] - x0) / xbar, 2));
        }
    }

    for (int i = 0; i < n; i ++) {
        if (x[i] < 0.6) {
            u[4*i+2] = 1.0;
        } else {
            u[4*i+2] = 4.0;
        }
        u[4*i+3] = 1.0;
    }

    iowrite("q00.dat", u, n);

    //
    //
    //
    // -- initialized --
    //
    //
    //

    int m = 0;
    while (1) {

        printf (">>> step %i to be taken. press any key to continue.\n", m);
	printf (">>> time passed is %f s.", m * dt);
        getchar ();

        u12[0] = u[0];
        u12[1] = u[1];
        u12[2] = u[2];
        u12[3] = u[3];
        u12[4*n] = u[4*(n-1)];
        u12[4*n+1] = u[4*(n-1)+1];
        u12[4*n+2] = u[4*(n-1)+2];
        u12[4*n+3] = u[4*(n-1)+3];

        for (int i = 1; i < n; i ++) {
            double pl = u[4*(i-1)];
            double ul = u[4*(i-1)+1];
            double rl = u[4*(i-1)+2];
            double kl = u[4*(i-1)+3];

            double pr = u[4*i];
            double ur = u[4*i+1];
            double rr = u[4*i+2];
            double kr = u[4*i+3];

            double hat_a = 0.5 * (1.0 / rl + 1.0 / rr);
            double hat_k = 0.5 * (kl + kr);

            double dq1 = pr - pl;
            double dq2 = ur - ul;

            f12[4*i]   = 0.5 * hat_k * dq2;
            f12[4*i+1] = 0.5 * hat_a * dq1;
            f12[4*i+2] = 0.0;
            f12[4*i+3] = 0.0;

            for (int j = 0; j < 4; j ++) {
                u12[4*i+j] = 0.5 * (u[4*(i-1)+j] + u[4*i+j]) - dt * f12[4*i+j] / dx;
            }
        }

        memset (f12, 0.0, (4*(n+1)) * sizeof (double));

	iowrite("q12.dat", u12, n+1);

        for (int i = 1; i < n + 1; i ++) {
            double pl = u12[4*(i-1)];
            double ul = u12[4*(i-1)+1];
            double rl = u12[4*(i-1)+2];
            double kl = u12[4*(i-1)+3];

            double pr = u12[4*i];
            double ur = u12[4*i+1];
            double rr = u12[4*i+2];
            double kr = u12[4*i+3];

            double star_a = 0.5 * (1.0 / rl + 1.0 / rr);
            double star_k = 0.5 * (kr + kl);

            double dq1 = pr - pl;
            double dq2 = ur - ul;

            f12[4*i]   = 0.5 * star_k * dq2;
            f12[4*i+1] = 0.5 * star_a * dq1;
            f12[4*i+2] = 0.0;
            f12[4*i+3] = 0.0;

            for (int j = 0; j < 4; j ++) {
	      unew[4*(i-1)+j] = 0.5 * (u12[4*(i-1)+j] + u12[4*i+j]) - dt * f12[4*i+j] / dx;
	      if (j == 2)
		unew[4*(i-1)+j] = u[4*(i-1)+j];
            }
        }

	iowrite("qne.dat", unew, n);

        memcpy (u,    unew, 4*n * sizeof (double));
        memset (u12,  0.0, (4*(n+1)) * sizeof (double));
        memset (f12,  0.0, (4*(n+1)) * sizeof (double));
        memset (unew, 0.0,  4*n * sizeof (double));

        m = m + 1;

    }


    free (unew);
    free (f12);
    free (u12);
    free (b);
    free (u);
    free (x);

    return 0;
}