/*

Compile using:
  gcc -o cairo `pkg-config --cflags --libs gtk+-2.0 gsl` cairo.c

*/

#include <stdlib.h>
#include <cairo.h>
#include <gtk/gtk.h>
#include <math.h>
#include <gsl/gsl_linalg.h>
     
#define GETX(n) (n)
#define GETY(n) (210-n)

#define N 20

double X[N][2] = { 
    { 10, 0 }, 
    { 25, 4 }, 
    { 50, -10 }, 
    { 75, 5 }, 
    { 100, 10 }, 
    { 120, 20 }, 
    { 150, -30 }, 
    { 175, 50 }, 
    { 200, 60 },
    { 220, 30 },
    { 250, 50 }, 
    { 275, 35 }, 
    { 300, 30 },
    { 325, 15 },
    { 350, 0 }, 
    { 375, 60 }, 
    { 400, 100 },
    { 425, 110 },
    { 450, 120 },
    { 455, 190 },
};

gboolean
on_expose_event(GtkWidget *widget,
    GdkEventExpose *event,
    gpointer data)
{
  cairo_t *cr;
  gint i;

  // draw line
  cr = gdk_cairo_create (widget->window);
  cairo_set_source_rgb(cr, 0, 0, 255);
  cairo_set_line_width(cr, 1);
  for ( i = 0; i < N; i++ ) {
      cairo_line_to(cr, GETX(X[i][0]), GETY(X[i][1]));
  }
  cairo_stroke(cr);
 
  // draw circles for data points
  cairo_set_source_rgb(cr, 0, 0, 0); 
  for ( i = 0; i < N; i++ ) {
      cairo_new_sub_path(cr);
      cairo_arc(cr, GETX(X[i][0]), GETY(X[i][1]), 3, 0, 2*M_PI);
      cairo_fill(cr);
  }


  // set bezier splines drawing style
  cairo_set_source_rgb(cr, 255, 0, 0);
  cairo_set_line_width(cr, 2);

  // Pi, Qi are control points for curve (Xi, Xi+1)
  double P[N-1][2], Q[N-1][2];

  /*****************************************************************************
  
   Bezier control points system matrix
  
    P0 P1 P2 Pn-1 ... Q0 Q1 Q2 Qn-1       
   /    1              1           \   /P0  \      /  2*X1\ Pi+1 + Qi = 2*Xi+1
   | 1  2             -2 -1        |   |P1  |      |     0|
   |       1              1        |   |P2  |      |  2*X2|
   |    1  2             -2 -1     | * |Pn-1|   =  |     0| Pi + 2*Pi+1
   |          1              1     |   |Q0  |      |2*Xn-1|    - Qi+1 - 2*Qi = 0
   |       1  2             -2 -1  |   |Q1  |      |     0|
   | 1                             |   |Q2  |      |    X0| P0   = X0     
   \                             1 /   \Qn-1/      \    Xn/ Qn-1 = Xn

                A*x = b
              x = inv(A)*b

       Pi, Qi and Xi are (x,y) pairs!
  
  *****************************************************************************/
  int s, eq = 0;
  gsl_matrix *m;
  gsl_vector *bx, *by, *x;
  gsl_permutation *perm;
   
  // allocate matrix and vectors
  m  = gsl_matrix_calloc(2*(N-1), 2*(N-1)); 
  bx = gsl_vector_calloc(2*(N-1)); 
  by = gsl_vector_calloc(2*(N-1)); 
  
  // fill-in matrix
  for ( i = 0; i < N-2; i++ ) {
      gsl_matrix_set(m, eq, i+1, 1);        // Pi+1
      gsl_matrix_set(m, eq, (N-1)+i, 1);    // + Qi
      eq++;                                 // = 2Xi+1

      gsl_matrix_set(m, eq, i, 1);          // Pi
      gsl_matrix_set(m, eq, i+1, 2);        // + 2*Pi+1
      gsl_matrix_set(m, eq, (N-1)+i+1, -1); // - Qi+1
      gsl_matrix_set(m, eq, (N-1)+i, -2);   // - 2*Qi
      eq++;                                 // = 0
  }
  gsl_matrix_set(m, eq++, 0, 1);            // P0   = X0
  gsl_matrix_set(m, eq++, 2*(N-1)-1, 1);    // Qn-1 = Xn

  
  // fill-in vectors
  for ( i = 0; i < N-2; i++ ) {
      gsl_vector_set(bx, 2*i, 2*X[i+1][0]);
      gsl_vector_set(by, 2*i, 2*X[i+1][1]);
  }
  gsl_vector_set(bx, 2*(N-1)-2, X[0][0]);
  gsl_vector_set(bx, 2*(N-1)-1, X[N-1][0]);

  gsl_vector_set(by, 2*(N-1)-2, X[0][1]);
  gsl_vector_set(by, 2*(N-1)-1, X[N-1][1]);
    
  // caluclate LU decomposition, solve lin. systems...  
  perm = gsl_permutation_alloc(2*(N-1));
  gsl_linalg_LU_decomp(m, perm, &s);

  // solve for bx
  x  = gsl_vector_calloc(2*(N-1));
  gsl_linalg_LU_solve( m, perm, bx, x ); 
  // copy solution (FIXME: should be avoided!)
  for ( i = 0; i < N-1; i++ ) {
      P[i][0] = gsl_vector_get(x, i);
      Q[i][0] = gsl_vector_get(x, i+(N-1));
  }
  gsl_vector_free(x);

  // solve for by
  x  = gsl_vector_calloc(2*(N-1));
  gsl_linalg_LU_solve( m, perm, by, x ); 
  // copy solution (FIXME: should be avoided!)
  for ( i = 0; i < N-1; i++ ) {
      P[i][1] = gsl_vector_get(x, i);
      Q[i][1] = gsl_vector_get(x, i+(N-1));
  }
  gsl_vector_free(x);


  gsl_permutation_free (perm);
  
  // free matrix and vectors
  gsl_matrix_free (m);
  gsl_vector_free (bx);
  gsl_vector_free (by);

  
  for ( i = 0; i < N-1; i++ ) {

// B second derivates  
//        printf("%d: Bx''(0) = %lf\n", i+1, 6*X[i][0]-12*P[i][0]+6*Q[i][0]);
//        printf("%d: Bx''(1) = %lf\n\n", i+1, 6*P[i][0]-12*Q[i][0]+6*X[i+1][0]);

// B first derivates
//        printf("%d: Bx'(0) = %lf\n", i+1, -3*X[i][0]+3*P[i][0]);
//        printf("%d: Bx'(1) = %lf\n", i+1, -3*Q[i][0]+3*X[i+1][0]);

        cairo_move_to(cr,
            GETX(X[i][0]), GETY(X[i][1]));

        cairo_curve_to(cr,
            GETX(P[i][0]), GETY(P[i][1]),
            GETX(Q[i][0]), GETY(Q[i][1]),
            GETX(X[i+1][0]), GETY(X[i+1][1]));                  
  }
  cairo_stroke(cr);
  cairo_destroy(cr);

  return FALSE;
}


int
main (int argc, char *argv[])
{

  GtkWidget *window;
  GtkWidget *darea;

  gtk_init(&argc, &argv);

  window = gtk_window_new(GTK_WINDOW_TOPLEVEL);
  
  darea = gtk_drawing_area_new();
  gtk_container_add(GTK_CONTAINER (window), darea);

  g_signal_connect(darea, "expose-event",
      G_CALLBACK(on_expose_event), NULL);
  g_signal_connect(window, "destroy",
      G_CALLBACK(gtk_main_quit), NULL);

  gtk_window_set_position(GTK_WINDOW(window), GTK_WIN_POS_CENTER);
  gtk_window_set_default_size(GTK_WINDOW(window), 465, 260); 

  gtk_widget_show_all(window);

  gtk_main();

  return 0;
}