/*******************************************************************
 * wave1d_spec.c: sequential 1d-wave solver with initial heights.
 *
 * To execute: mpicc wave1d_spec.c ; mpiexec -n 4 ./a.out Or replace
 * "4" with however many procs you want to use.
 *
 * To verify: civl verify wave1d_spec.
 * Author: Ziqing Luo
 ********************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <math.h>

#define SQR(x) ((x)*(x))

/* Input parameters */
#ifdef _CIVL
#include <civlc.cvh>

$input int NXB = 5;               /* upper bound on nx */
$input int nx;                    /* number of discrete points including endpoints */
$assume(2 <= nx && nx <= NXB);    /* setting bounds */
$input double c;                  /* physical constant to do with string */
$assume(c > 0.0);       
$input int NSTEPSB = 5;        
$input int nsteps;                /* number of time steps */
$assume(0 < nsteps && nsteps <= NSTEPSB);
$input int wstep = 1;             /* number of time steps between printing */
$input double u_init[nx];         /* arbitrary initial position of string */
$output double output[nx];

#else

int nx, height_init, width_init;
int nsteps, wstep;
double c;

#endif

/* Global varibales */
double *u_prev, *u_curr, *u_next;
double k;

void update() {
  int i;
  double *tmp;

  for (i = 1; i < nx + 1; i++){
    u_next[i] = 2.0*u_curr[i] - u_prev[i] +
      k*(u_curr[i+1] + u_curr[i-1] -2.0*u_curr[i]);
  }
  //cycle pointers:
  tmp = u_prev;
  u_prev = u_curr;
  u_curr = u_next;
  u_next = tmp;
}

/* Initialization function, initializes all parameters and data
   array. */
void initialization() {
  int i, j;

#ifndef _CIVL

  double e = exp(1.0);

  nx = 50;
  c = 0.3;
  height_init = 10;
  width_init = 10;
  nsteps = 500;
  wstep = 5;
  // intialize the initial heights with some arbitary values:
  for(i = 1; i < nx+1; i++) {
    if(idx >= width_init)
      u_prev[i] = 0.0;
    else
      u_prev[i] = height_init * e *
	exp(-1.0/(1-SQR(2.0*(idx-width_init/2.0)/width_init)));
  }

#endif

  assert(nx >= 2);
  assert(c > 0);
  assert(nsteps >= 1);
  assert(wstep >= 1 && wstep <= nsteps);
  k = c * c;
  printf("Wave1d with nx=%d, c=%f, nsteps=%d, wstep=%d\n", nx, c, nsteps, wstep);
  u_prev = (double *)malloc((nx + 2) * sizeof(double));
  assert(u_prev);
  u_curr = (double *)malloc((nx + 2) * sizeof(double));
  assert(u_curr);
  u_next = (double *)malloc((nx + 2) * sizeof(double));
  assert(u_next);
  // sets constant boundaries
  u_prev[0] = 0;
  u_curr[0] = 0;
  u_next[0] = 0;
  u_prev[nx + 1] = 0;
  u_curr[nx + 1] = 0;
  u_next[nx + 1] = 0;
}

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

  initialization();
#ifdef _CIVL

  for (int i = 0; i < nx; i++)
    u_prev[i+1] = u_init[i];

#endif
  for (int i = 0; i < nx; i++)
    u_curr[i+1] = u_prev[i+1];
  for(iter = 0; iter < nsteps; iter++)
    update();
#ifdef _CIVL

  memcpy(output, &u_curr[1], sizeof(double) * nx);

#endif  
  free(u_curr);
  free(u_prev);
  free(u_next);
  return 0;
}
