#include "/usr/local/dislin/dislin.h"

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

#include "../NR/nrutil.c"

#define round(x) ((x)>=0?(long)((x)+0.5):(long)((x)-0.5))

double x_0 = -10.0;
double x_1 = 10.0;
double T = 5.0;

double h = 0.05;
double k = 0.0000125;

double sech(double x) {
  return 2.0 / (exp(x) + exp(-x));
}

int main () {

  long x_points = round((x_1 - x_0) / h);
  long t_points = round(T / k);
  
  double *u_lat = dvector(1, x_points);
  double *u_old = dvector(1, x_points);
  double *u_older = dvector(1, x_points);

  double two_k_over_h = 2.0 * (k / h);
  double k_over_h_cube = k / pow(h,3.0);
  
  long x_plus = 0;
  long x_plusplus = 0;
  long x_minus = 0;
  long x_minusminus = 0;

  int i, j;

  char title_str[80];
  char fname[20];

  metafl("PNG");
 
  for (j = 1; j<= x_points; j++) {
    u_old[j] = -2.0*pow(sech(x_0+(j-1)*h),2);
    u_older[j] = u_old[j];
  }

	disini();
	texmod("ON");

	sprintf(title_str, "%d Stuetzstellen, $t=%f$",x_points, 0.0);
	titlin(title_str, 4);

	graf(x_0,x_1,x_0,0.1*(x_1-x_0),-5.0,5.0,-5.0,1.0);
	title();

	hsymbl(5);
	for (j=1; j<=x_points; j++) {
	  rlsymb(21, x_0 + (j-1)*h, u_old[j]);
	}

	disfin();

  for (i = 1; i <= t_points; i++) {
    for (j = 1; j <= x_points; j++) {

      x_minus = j-1;
      x_minusminus = j-2;
      x_plus = j+1;
      x_plusplus = j+2;

      if ( j == 2 ) {
	x_minus = 1;
	x_minusminus = x_points;
      }

      if ( j == 1 ) {
	x_minus = x_points;
	x_minusminus = x_points - 1;
      }
      
      if ( j == x_points ) {
	x_plus = 1;
	x_plusplus = 2;
      }

      if ( j == x_points-1 ) {
	x_plus = j+1;
	x_plusplus = 1;
      }

      u_lat[j] = two_k_over_h * (u_old[x_plus] - u_old[x_minus]) * (u_old[j] + u_old[x_plus] + u_old[x_minus])-	k_over_h_cube * (u_old[x_plusplus]-2.0*u_old[x_plus]+2.0*u_old[x_minus]-u_old[x_minusminus]) + u_older[j];
      //      printf("u_lat[%d]=%f\n",j,u_lat[j]);
    }

      for (j = 1; j<= x_points; j++) {
	u_older[j] = u_old[j];
	u_old[j] = u_lat[j];
      }

      if ( i % 500 == 0 ) {
	sprintf(fname, "kdv_%d_frame%d.png",x_points,i / 500);
	setfil(fname);
	disini();
	texmod("ON");

	sprintf(title_str, "%d Stuetzstellen, $t=%f$",x_points,i*k);
	titlin(title_str, 4);

	sprintf(fname, "kdv_%d_frame%d",x_points,i / 500);

	graf(x_0,x_1,x_0,0.1*(x_1-x_0),-5.0,5.0,-5.0,1.0);
	title();

	hsymbl(5);
	for (j=1; j<=x_points; j++) {
	  rlsymb(21, x_0 + (j-1)*h, u_lat[j]);
	}

	disfin();
      }

  }

  free_dvector(u_lat, 1, x_points);
  free_dvector(u_old, 1, x_points);
  free_dvector(u_older, 1, x_points);

  return 0;

}