// Python - C extension for finite_volumes util module.
//
// To compile (Python2.3):
//  gcc -c util_ext.c -I/usr/include/python2.3 -o util_ext.o -Wall -O 
//  gcc -shared util_ext.o  -o util_ext.so	
//
// See the module util.py
//
//
// Ole Nielsen, GA 2004
	
#include "Python.h"
#include "Numeric/arrayobject.h"
#include "math.h"


int _gradient(double x0, double y0, 
	      double x1, double y1, 
	      double x2, double y2, 
	      double *q0, double *q1, double *q2, 
	      double *a, double *b, 
	      int N) {

  double det;
  int i;
  
  det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0);
  
  for (i=0; i<N; i++) {
    a[i] = (y2-y0)*(q1[i]-q0[i]) - (y1-y0)*(q2[i]-q0[i]);
    a[i] /= det;

    b[i] = (x1-x0)*(q2[i]-q0[i]) - (x2-x0)*(q1[i]-q0[i]);
    b[i] /= det;
  }  

  return 0;
}



// Computational function for rotation
int _rotate(double *q, double *r, double *normal, int direction) {
  /*Rotate the momentum component q (q[1], q[2])
    from x,y coordinates to coordinates based on normal vector.
    
    To rotate in opposite direction, call rotate with direction=1
    Result is returned in r 
    
    q and r are assumed to have three components each*/


  double n0, n1, q1, q2;
  
  //Determine normal and direction
  n0 = normal[0];
  if (direction == 1) {
    n1 = normal[1];    
  } else {
    n1 = -normal[1];
  }

  //Shorthands
  q1 = q[1];  //uh momentum
  q2 = q[2];  //vh momentum

  //Rotate
  r[0] = q[0];
  r[1] = n0*q1 + n1*q2;
  r[2] = n0*q2 - n1*q1;

  return 0;
}  



// Gateway to Python
PyObject *gradient(PyObject *self, PyObject *args) {
  //
  // a,b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
  //
  
  double x0, y0, x1, y1, x2, y2;
  PyObject *Q0, *Q1, *Q2, *result;
  PyArrayObject *q0, *q1, *q2, *a, *b;
  int dimensions[1];
  int N;

  
  // Convert Python arguments to C  
  if (!PyArg_ParseTuple(args, "ddddddOOO", &x0, &y0, &x1, &y1, &x2, &y2,
			&Q0, &Q1, &Q2))  
    return NULL;

  //Input check
  if PyArray_Check(Q0) {
    q0 = (PyArrayObject *) 
      PyArray_ContiguousFromObject(Q0, PyArray_DOUBLE, 0, 0);
  } else {
    PyErr_SetString(PyExc_TypeError, "q0 must be a numeric array");
    return NULL;
  }
    
  if PyArray_Check(Q1) {
    q1 = (PyArrayObject *) 
      PyArray_ContiguousFromObject(Q1, PyArray_DOUBLE, 0, 0);
  } else {
    PyErr_SetString(PyExc_TypeError, "q1 must be a numeric array");
    return NULL;
  }
  
  if PyArray_Check(Q2) {
    q2 = (PyArrayObject *) 
      PyArray_ContiguousFromObject(Q2, PyArray_DOUBLE, 0, 0);
  } else {
    PyErr_SetString(PyExc_TypeError, "q2 must be a numeric array");
    return NULL;
  }  
    
  //Allocate space for return vectors a and b  
  
  N = q0 -> dimensions[0];
  dimensions[0] = N;
  a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
  b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);  

  
  // Call underlying routine
  _gradient(x0, y0, x1, y1, x2, y2, 
	    (double *) q0 -> data, 
	    (double *) q1 -> data, 
	    (double *) q2 -> data,
	    (double *) a -> data, 
	    (double *) b-> data, N);

  Py_DECREF(q0);
  Py_DECREF(q1);
  Py_DECREF(q2);    
	      
  // Return arrays a and b
  result = Py_BuildValue("OO", a, b);
  Py_DECREF(a);
  Py_DECREF(b);  
  return result;
}


// Gateway to Python
PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
  //
  // r = rotate(q, normal, direction=0)
  //
  // Where q is assumed to be a Float numeric array of length 3 and
  // normal a Float numeric array of length 2.

  //FIXME: Input checks and unit tests
  
  PyObject *Q, *Normal;
  PyArrayObject *q, *r, *normal;
  
  static char *argnames[] = {"q", "normal", "direction", NULL};
  int dimensions[1];
  int N, direction=0;

  // Convert Python arguments to C  
  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames, 
				   &Q, &Normal, &direction))  
    return NULL;  

  //Input checks (convert sequenced into numeric arrays)
  
  q = (PyArrayObject *) 
    PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
  normal = (PyArrayObject *) 
    PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
  
  //Allocate space for return vector r (don't DECREF) 
  N = q -> dimensions[0];
  dimensions[0] = N;
  r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);

  //Rotate
  _rotate((double *) q -> data,
	  (double *) r -> data,
	  (double *) normal -> data,	  
	  direction);	  

  //Build result and DECREF r - returning r directly causes memory leak
  //result = Py_BuildValue("O", r);	  
  //Py_DECREF(r);
  
  //return result;
  return PyArray_Return(r);
}    


// Method table for python module
//static struct PyMethodDef MethodTable[] = {
//  {"gradient", gradient, METH_VARARGS},  
//  {NULL, NULL}
//};



// Method table for python module
static struct PyMethodDef MethodTable[] = {
  /* The cast of the function is necessary since PyCFunction values
   * only take two PyObject* parameters, and rotate() takes
   * three.
   */
  
  {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},  
  {"gradient", gradient, METH_VARARGS, "Print out"},    
  {NULL, NULL, 0, NULL}   /* sentinel */
};



// Module initialisation   
void initutil_ext(void){
  Py_InitModule("util_ext", MethodTable);
  
  import_array();     //Necessary for handling of NumPY structures  
}