/* Das Programm integriert die 2D Wellengleichung und gibt die
   Ergebnisse in Form von PNG-Grafiken aus, aus denen sich z.B. mit
   Programmen wie ffmpeg kleine Filme erzeugen lassen. Die Anzahl der
   ausgegebenen Bilder läßt sich über die Variable "fps" steuern
   (siehe unten). */

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

// Die folgende Zeile muss evtl. an Eure lokale Dislin-Installation
// angepasst werden.
#include "/usr/local/dislin/dislin.h"

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

// Zeitintervall [t_0,t_0]
double t_0 = 0.0;
double t_1 = 20.0;
// Raumbereich [x_0,x_1]x[y_0,y_1]
double x_0 = 0.0;
double x_1 = 1.0;
double y_0 = 0.0;
double y_1 = 1.0;

// Raum- und Zeitschritt. Aus diesen Größen wird weiter unten die
// Anzahl der Schritte berechnet.
double dxy = 0.01;
double dt = 0.0025;

// Die Ausbreitungsgeschwindigkeit 
double c = 1.0;

// fps ist die Anzahl der Bilder, die pro simulierter Sekunde ausgegeben werden sollen
int fps = 20;

// Die Funktionen startwert und startableitung werden zu Beginn des
// Programms verwendet, um die Anfangsbedingung festzulegen.
// In der Notation der Vorlesung entsprechen sie f(x) und g(x).
double startwert(int i, int j) {
  int n = 4;

  double y = y_0 + i * dxy;
  double x = x_0 + j * dxy;

  //  return sin(n * 2.0*M_PI / (y_1 - y_0) * (y_0 + i*dxy)) * sin((double)n*M_PI / (x_1 - x_0) * (x_0 + j*dxy));
  return 4*pow(x,2)*y*(1-x)*(1-y);
}

double startableitung(int i, int j) {
  return 0.0;
}

int main () {
  // Die Variablen fname und titel werden für die Ausgabe mit Dislin
  // benötigt (siehe unten).
  char fname[80];
  char titel[80];

  // Die Anzahl der räumlichen und zeitlichen Iterationsschritte wird
  // aus der Größe der entsprechenden Intervalle und den Schrittweiten
  // berechnet.
  int x_steps = round( (x_1 - x_0) / dxy); 
  int y_steps = round( (y_1 - y_0) / dxy);
  int t_steps = round( (t_1 - t_0) / dt);

  // Die zweifache Zeitableitung wird über den zentralen
  // Differenzenquotienten approximiert. Das heißt, dass jeweils die
  // Funktionswerte von drei Zeitpunkten benötigt werden: Der aktuelle
  // Funktionswert wird aus den Funktionswerten der beiden
  // vorangegangenen Zeitpunkte ermittelt.
  double u[y_steps + 1][x_steps + 1];
  double u_alt[y_steps + 1][x_steps + 1];
  double u_aelter[y_steps + 1][x_steps + 1];

  // Je nach Wahl der globalen Variablen fps und dt wird die Zahl der
  // zeitlichen Iterationen berechnet, die zwischen zwei Ausgaben
  // durchgeführt werden sollen.
  int steps_pro_frame = round(1.0 / fps / dt);
  if ( !steps_pro_frame ) {
    steps_pro_frame = 1;
  }
  
  int i = 0; // i ist der y-Index
  int j = 0; // j ist der x-Index
  int k = 0; // k ist der Zeitindex

  // Division und Potenzierung ist relativ teuer, sollte deshalb nicht
  // öfter als nötig berechnet werden!
  double alpha_quadrat = pow(c * dt / dxy, 2);

  // Das Feld wird entsprechend der Anfangsbedingung mit startwert()
  // initialisiert.
  for (i = 0; i <= y_steps; i++) {
    for (j = 0; j <= x_steps; j++) {
      u_aelter[i][j] = startwert(i, j);
      u[i][j] = u_aelter[i][j];
    }
  }

  int y_minus = 0;
  int y_plus  = 0;
  int x_minus = 0;
  int x_plus  = 0;

  for ( k = 0; k <= t_steps; k++ ) {
    // Die Raumindizes gehen von 1 bis xy_steps-1, da feste
    // Randbedingungen erfüllt sein sollen
    for ( i = 1; i <= y_steps-1; i++ ) {

      y_minus = i-1;
      y_plus  = i+1;

      for ( j = 1; j <= x_steps-1; j++) {
	
	x_minus = j-1;
	x_plus  = j+1;

	// Im ersten Zeitschritt müssen die u_i,j aus Zeitschritt n=-1
	// mit Hilfe der ersten Zeitableitung (Anfangsbedingung g(x))
	// konstruiert werden.
	if ( k == 0 ) {
	  u_alt[i][j] = dt * startableitung(i, j) + (1.0 - 2.0 * alpha_quadrat) * u_aelter[i][j] + 0.5 * alpha_quadrat * (u_aelter[y_plus][j] + u_aelter[y_minus][j] + u_aelter[i][x_minus] + u_aelter[i][x_plus]);
	} else {
	  u[i][j] = -u_aelter[i][j] + 2.0 * (1.0-2.0*alpha_quadrat)*u_alt[i][j] + alpha_quadrat * (u_alt[y_plus][j] + u_alt[y_minus][j] + u_alt[i][x_plus] + u_alt[i][x_minus]);
	}
      }
    }
    
    if (k != 0) {
      for (i = 0; i <= y_steps-1; i++) {
	for (j = 0; j <= x_steps-1; j++) {
	  u_aelter[i][j] = u_alt[i][j];
	  u_alt[i][j] = u[i][j];
	}
      }
    }

    printf("steps: %d\n", steps_pro_frame);
    // ---{ Ausgabe mit Dislin }-------------------------
    if ( k % steps_pro_frame == 0) {
      metafl("PNG");
      filmod("DELETE");
      sprintf(fname, "2Dwellengl_frame%d.png", round(k / steps_pro_frame) );
      setfil(fname);
      disini();
      texmod("ON");
      
      sprintf(titel, "%dx%d Stuetzstellen, $t=%f$",x_steps + 1, y_steps + 1, t_0 + (k + 1) * dt);
      titlin(titel, 4);
    
      graf3d(x_0, x_1, x_0, (x_1 - x_0)/10.0, y_0, y_1, y_0, (y_1 - y_0)/10.0, -0.2, 0.2, -0.2, 0.1);
      title();

      float darst[y_steps + 1][x_steps + 1];
      for (i = 0; i<=y_steps; i++) {
        for (j = 0; j<=x_steps; j++) {
	        darst[i][j] = (float) u[i][j];
	      }
      }

      surmat((float *) darst, y_steps + 1, x_steps + 1, 1, 1);
      disfin();
    }
    //-----------------------------------------------------

  }
  
  return 1;
  
}