/* Das Programm integriert die Advektionsgleichung 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 = 50.0;
// Raumintervall [x_0,x_1]
double x_0 = 0.0;
double x_1 = 10.0;

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

// Die Advektionskonstante 
double c = 1.0;

// Framerate - So viele Bilder werden pro simulierter Sekunde ausgegeben
int fps = 20;


// Die Funktion startwert(i) legt die Anfangsbedingung fest, es gilt
// also u(x_i,t=0) = startwert(i).
double startwert(int i) {
  int n = 4;

  return sin(n * 2.0*M_PI / (x_1 - x_0) * (x_0 + i*dx));
}

int main () {

  // Die folgenden beiden Strings werden für den Namen der
  // Ausgabedatei und für die Titelzeile der ausgegebenen Graphen
  // reserviert.
  char fname[25];
  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) / dx); 
  int t_steps = round( (t_1 - t_0) / dt);

  // Die Funktionswerte u(x_i) an den Stützszellen werden in einem
  // Array gespeichert. Gespeichert wird jeweils der aktuelle Wert (in
  // u) und derjenige vom vorangegangenen Zeitschritt (in u_alt).
  double u[x_steps + 1];
  double u_alt[x_steps + 1];

  // Berechne, nach wie vielen Zeitschritten ein Bild ausgegeben
  // werden soll.
  int steps_pro_frame = round(1.0 / fps / dt);
  if (steps_pro_frame == 0) {
    steps_pro_frame = 1;
  }
  
  int i = 0; // Indiziert den Ort
  int j = 0; // Indiziert die Zeit

  // Über x_plus und x_minus werden weiter unten die periodischen
  // Randbedingungen realisiert.
  int x_plus = 0;
  int x_minus = 0;

  // Multiplikation ist relativ teuer, sollte deshalb nicht öfter als
  // nötig berechnet werden!
  double c_mal_dt_durch_dx = c * dt / dx;

  // Initialisiere das Gitter (Anfangsbedingung)
  for (i = 0; i < x_steps; i++) {
    u[i] = startwert(i);
  }

  // Durchlaufe alle Zeitschritte von t_0 bis t_1.
  for ( j = 0; j <= t_steps; j++ ) {

    //--- Ausgabe ---------------------------------------------
    if ( j % steps_pro_frame == 0) {
      metafl("PNG");
      filmod("DELETE");
      scrmod("REVERSE");
      // Die ausgegebenen Bilder sollen advektion_frame1.png,
      // advektion_frame2.png usw, heißen
      sprintf(fname, "advektion_frame%d.png", round(j / steps_pro_frame) );
      setfil(fname);
      disini();
      texmod("ON");
      
      // Lege die Überschrift des Graphen fest
      sprintf(titel, "%d Stuetzstellen, $t=%f$",x_steps + 1, t_0 + j * dt);
      titlin(titel, 4);
      
      graf(x_0, x_1, x_0, (x_1 - x_0)/10, -3.0, 3.0, -3.0, 1.0);
      title();
      
      // Plotte die Funktionswerte
      hsymbl(5);
      for (i = 0; i <= x_steps; i++ ) {
	rlsymb(21,x_0 + i * dx, u[i]);
      }
      
      disfin();
    }
    //--- Ausgabe Ende -------------------------------------------
    
    // Kopiere den aktuellen Zustand. Daraus wird dann der neue 
    // Zustand berechnet.
    for ( i = 0; i <= x_steps; i++ ) {
      u_alt[i] = u[i];
    }

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

      // Periodische RB 
      // Ist die Schleife am rechten Rand angelangt, so wird der linke
      // Rand als nächster rechter Nachbar behandelt
      if ( i == x_steps ) {
	x_plus = 0;
      }
      // Ist die Schleife am linken Rand angelangt, so wird der rechte
      // Rand als nächster linker Nachbar behandelt
      if ( i == 0 ) {
	x_minus = x_steps;
      }

      // Jetzt erfolgt die Zeitentwicklung. Wie in der Vorlesung
      // besprochen, wird das Upwind-Verfahren genutzt um für dt <
      // dt/|c| die Stabilität zu sichern.
      
      // Upwind-Verfahren
      // ================
      // Advektionssgeschw. c |  Örtlicher Differenzenquotient 
      //                      |
      //         < 0          |            forward
      //         > 0          |            backward
      if ( c > 0 ) {
	u[i] =  u_alt[i] - c_mal_dt_durch_dx * (u_alt[i] - u_alt[x_minus]);
      } else {
	u[i] =  u_alt[i] - c_mal_dt_durch_dx * (u_alt[x_plus] - u_alt[i]);
      }
  
    }


  }
  
  return 1;
  
}