// Dieses Programm integriert die Sine-Gordon-Gleichung 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). Für die lineare Algebra wird
// dabei die Bibliothek GLS (Gnu Scientific-Library) verwendet, weswegen mit
// "-ldislin -lgsl -lgslcblas" gelinkt werden muss.

#include <math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

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

// Eine nützliche Definition zum Runden
#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 = 1800.0;
// Raumbereich [x_0,x_1]x[y_0,y_1]
double x_0 = -20.0;
double x_1 = 20.0;

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

// Ausbreitungsgeschwindigkeit (es muss |c|<1 gelten!)
double c = 0.5;

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

// 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). Es sind jeweils
// die Anfangsbedingungen für Kink, Antikink und Breather angegeben. Durch
// Ändern der Kommentare könnt Ihr die verschiedenen Lösungstypen
// durchprobieren.
double startwert(int i) {
  double x = x_0 + i * dx;

  double kinkpos = 0.0;

  // Kink
  return 4.0 * atan(exp((x-kinkpos)/sqrt(1-pow(c,2))));

  // Antikink
  //return 4.0 * atan(exp(-(x-kinkpos)/sqrt(1-pow(c,2))));

  // Breather
  //return 0.0;
}

double startableitung(int i) {
  double x = x_0 + i * dx;

  // Breather
  //return 4.0 * sqrt(1-pow(c,2))/cosh(x*sqrt(1-pow(c,2)));

  // Kink
  return 2.0 * c / sqrt(1-pow(c,2)) / cosh( -x / sqrt(1-pow(c,2)) );

  // Antikink
  return 2.0 * c / sqrt(1-pow(c,2)) / cosh( x / sqrt(1-pow(c,2)) );
.0;
}

int main () {
 // 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 Variablen fname und titel werden für die Ausgabe mit Dislin
  // benötigt (siehe unten).
  char fname[80];
  char titel[80];

  // 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;
  }
  
  // Die Zeitintegration erfolgt durch Lösung eines linearen
  // Gleichungssystems. Im Folgenden werden die dafür benötigten Vektoren und
  // eine Matrix definiert.
  gsl_vector *u = gsl_vector_alloc(x_steps + 1);
  gsl_vector *u_alt = gsl_vector_alloc(x_steps + 1);
  gsl_vector *gamma = gsl_vector_alloc(x_steps + 1);
  gsl_vector *beta = gsl_vector_alloc(x_steps + 1);
  gsl_matrix *A = gsl_matrix_alloc(x_steps + 1, x_steps + 1);

  int i = 0; // i 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(dt / dx, 2);//pow(c * dt / dx, 2);

  // Das Feld wird entsprechend der Anfangsbedingung mit startwert()
  // initialisiert.
  for (i = 0; i <= x_steps; i++) {
      gsl_vector_set(u, i, startwert(i));
      gsl_vector_set(beta, i, sin(gsl_vector_get(u, i)));
      gsl_vector_set(gamma, i, startableitung(i));
  }

  // Erstelle die Matrix A
  for ( i=0; i<=x_steps; i++ ) {
    gsl_matrix_set(A, i, i, 1-alpha_quadrat);

    if ( (i != 0) && (i != x_steps) ) {
      gsl_matrix_set(A, i, i-1, 0.5 * alpha_quadrat);
      gsl_matrix_set(A, i, i+1, 0.5 * alpha_quadrat);
    } else if ( i == 0 ) {
      gsl_matrix_set(A, i, i+1, alpha_quadrat);
    } else {
      gsl_matrix_set(A, i, i-1, alpha_quadrat);
    } 
  }


  //---{ Erster Zeitschritt }------------------------------------------- 
  gsl_vector_scale(gamma, dt);
  gsl_vector_scale(beta, -0.5 * pow(dt, 2));

  gsl_vector_memcpy(u_alt, u);
  gsl_blas_dgemv(CblasNoTrans, 1.0, A, u_alt, 0.0, u);

  gsl_vector_add(u, gamma);
  gsl_vector_add(u, beta);
  //--------------------------------------------------------------------

  for ( k = 1; k <= t_steps; k++ ) {

    //---{ Ausgabe mit Dislin }-----------------------------------------
    if ( k % steps_pro_frame == 0) {
      scrmod("REVERSE");
      metafl("PNG");
      filmod("DELETE");
      sprintf(fname, "SineGordon_frame%d.png", round(k / steps_pro_frame) );
      setfil(fname);
      disini();
      texmod("ON");
      
      sprintf(titel, "%d Stuetzstellen, $t=%f$",x_steps + 1, t_0 + (k + 1) * dt);
      titlin(titel, 4);
    
      // Graph für Darstellung von Kink-/Antikink und Breather-Lösungen
      //graf(x_0, x_1, x_0, (x_1 - x_0)/10.0, -7.0, 7.0, -7.0, 1.0);
      // Graph für die Darstellung von Antikink-Reflektionen
      graf(x_0, x_1, x_0, (x_1 - x_0)/10.0, -1.0, 28.0, -1.0, 1.0);
      title();

      hsymbl(6);
      for (i=0; i<x_steps; i++) {
        rlsymb(21, x_0 + i * dx, gsl_vector_get(u, i));
      }

      disfin();
    }
    //------------------------------------------------------------------


    //---{ Der eigentliche Zeitschritt }-------------------------------
    for ( i=0; i<=x_steps; i++ ) {
      gsl_vector_set(beta, i, sin(gsl_vector_get(u, i)));
    }

    gsl_vector_memcpy(gamma, u_alt);

    gsl_vector_scale(gamma, -1.0);
    gsl_vector_scale(beta, -pow(dt, 2));

    gsl_vector_memcpy(u_alt, u); 
    gsl_blas_dgemv(CblasNoTrans, 2.0, A, u_alt, 0.0, u);

    gsl_vector_add(u, gamma);
    gsl_vector_add(u, beta);
    //-----------------------------------------------------------------
  }
  
  // Zuletzt kann der reservierte Speicher wieder freigegeben werden.
  gsl_vector_free(u);
  gsl_vector_free(u_alt);
  gsl_vector_free(beta);
  gsl_vector_free(gamma);
  gsl_matrix_free(A);

  return 1;
  
}