#include <stdio.h>
#include <gsl/gsl_math.h>
#include <math.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_vector_complex_double.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>

typedef struct {
  double w;
  gsl_complex mu;
  gsl_complex eps;
  gsl_complex k;
} constants;

typedef struct {
          gsl_matrix * nodes;
          gsl_matrix * triangles;
          gsl_matrix * edges;
} mesh;

int index_find (int node_number, gsl_matrix * nodes) {
  int index = 0;
  while (gsl_matrix_get (nodes, index, 0) != node_number) index++;
  return index;
}

int
main (void)
    {
  constants c;
  c.w = 2e9;
  c.mu = gsl_complex_rect(1, 0);
  c.eps = gsl_complex_rect(1,0);
  c.k = gsl_complex_rect(50, 0);
  mesh mm;
  int i, j; 
  mm.nodes = gsl_matrix_calloc(10, 4);
  mm.triangles = gsl_matrix_calloc(10, 3);
  mm.edges = gsl_matrix_calloc(10, 6);
  
  for (i = 0; i < 10; i++)
    for (j = 0; j < 4; j++)
      gsl_matrix_set (mm.nodes, i, j, 0.23 + 100*i + j);
  
 for (i = 0; i < 10; i++)
    for (j = 0; j < 3; j++)
      gsl_matrix_set (mm.triangles, i, j, 0.33 + 100*i + j);


  for (i = 0; i < 10; i++)
    for (j = 0; j < 3; j++)
      printf("mm.triangles(%d,%d) = %g\n", i, j, 
             gsl_matrix_get (mm.triangles, i, j));

  gsl_vector_complex * A = gsl_vector_complex_calloc (10);
  gsl_vector_complex_sub (A, A);
  gsl_vector_complex_scale (A, gsl_complex_rect (1, 0));

  return 0;

}