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

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

#include "NR/ansi/recipes/bandec.c"
#include "NR/ansi/recipes/banbks.c"

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

// Gitterparameter: a / h und k / T müssen Integer sein!

float h = 0.1;
float k = 0.1;

float a = 1.0;
float T = 1.0;

float max_temp = 5.0;

//            t ^           
//              |            |
//              |            |
//    g_1(t)    |            |  g_2(t)
//              |            |
//              |            |
//              |------------|-->
//              0    f(x)    a    x

float g_1(float t) {
  return 0.0;
}

float g_2(float t) {
  return max_temp;
}

float f(float x) {
  return 1.0;
}

int main () {

  int x_steps = (int) a / h;
  int t_steps = (int) T / k;

  float lambda = k / pow(h,2);

  if ( lambda > 0.5 ) {
    printf("Lambda is %f --> Code probably unstable!\n", lambda);
  } else {
    printf("Lambda is %f --> Code should be stable!\n", lambda);
  }

  float **lattice = matrix(0,x_steps,0,t_steps);


  int i = 0;
  int j = 0;
  
  // Randbedingungen eintragen

  for ( i = 0; i <= x_steps; i++ ) {
    lattice[i][0] = f(i * h);
  }

  for ( i = 0; i <= t_steps; i++ ) {
    lattice[0][i] = g_1(i * k); 
  }

  for ( i = 0; i <= t_steps; i++ ) {
    lattice[x_steps][i] = g_2(i * k); 
  }

  // Implizite Lösung, LGS aufstellen

  float *b = fvector(1,(x_steps - 1)*(t_steps));
  float **matr = matrix(1,(x_steps - 1)*(t_steps),1, (x_steps - 1) * (t_steps));


  float *indx = fvector(1,(x_steps - 1)*(t_steps));
  float d;

  int lat_row = 0;
  for ( i = 1; i <= (x_steps - 1)*(t_steps); i++ ) {

    matr[i][i] = -(1.0 + 2.0 * lambda);
    matr[i][i+1] = lambda;
    matr[i][i-1] = lambda;

    if ( i > 3 ) {
      matr[i][i-3] = 1.0;
    }    

    if (i <= (x_steps-1))
      b[i] -= f(i*h);

    if (i == 1 + lat_row * (x_steps-1) ) {

      matr[i][i-1] = 0.0;
      if (lat_row > 0) {
	matr[i-1][i] = 0.0;
      }
      
      lat_row++;
      b[i] -= lambda * g_1(lat_row*k);      
    }

    if (i == lat_row * (x_steps - 1))
      b[i] -= lambda * g_2(lat_row * k );

  }

  // Darstellen

  /*  for ( i = 1; i <= (x_steps - 1)*(t_steps-1); i++) {
    for ( j = 1; j <= (x_steps - 1)*(t_steps-1); j++) {
      if (matr[j][i] < 0) {
	printf("%.2f ", matr[i][j]);
      } else {
	printf(" %.2f ", matr[i][j]);
	//	printf("%d%d ",i,j );
      }
      
    }
    printf(" | %f \n", b[i]);
    }*/

  // und lösen...

  ludcmp(matr, (x_steps - 1)*(t_steps), indx, &d);
  lubksb(matr, (x_steps - 1)*(t_steps), indx, b);

  // Ergebnis in das Gitter eintragen

  float *iter = &b[1];

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

    for ( j = 1; j <= x_steps - 1; j++ ) {
      
      lattice[j][i] = *iter;
      iter++;
      //lattice[j][i] = b[(i - 1) * x_steps + j]; 

    }

  }

  metafl("CONS");
  disini();

  name("x","X");
  name("t","Y");
  name("H","Z");

  graf3(0,a,0,0.5,0,T,0,0.5,0,max_temp,0,0.5);

  color("GREEN");
  crvmat((float *) &lattice[0][0], x_steps+1, t_steps+1, 10, 10);

  disfin();
  
  /*for ( j = x_steps; j >= 0; j--) {
    for ( i = 0; i <= t_steps; i++) {
      if (lattice[j][i] < 0) {
	printf("%.2f ", lattice[j][i]);
      } else {
	printf(" %.2f ", lattice[j][i]);
	//	printf("%d%d ",i,j );
      }

    }
    printf("Line!\n");
    }*/

  free_matrix(lattice, 0,x_steps,0,t_steps);
  free_matrix(matr,1,(x_steps-1)*(t_steps-1),1,(x_steps - 1)*(t_steps -1));
}