// 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"

//Shared snippets
#include "util_ext.h"



int _point_on_line(double x, double y, 
		   double x0, double y0, 
		   double x1, double y1) {
  /*Determine whether a point is on a line segment 
   
    Input: x, y, x0, x0, x1, y1: where
        point is given by x, y
	line is given by (x0, y0) and (x1, y1)
   	
  */
    
  double a0, a1, a_normal0, a_normal1, b0, b1, len_a, len_b;
  
  a0 = x - x0;
  a1 = y - y0; 
  
  a_normal0 = a1;
  a_normal1 = -a0;
           		
  b0 = x1 - x0;
  b1 = y1 - y0;
    
  if ( a_normal0*b0 + a_normal1*b1 == 0 ) {   
    //Point is somewhere on the infinite extension of the line

    len_a = sqrt(a0*a0 + a1*a1);    	       
    len_b = sqrt(b0*b0 + b1*b1);
        	       	                
    if (a0*b0 + a1*b1 >= 0 && len_a <= len_b) {
      return 1;
    } else {
      return 0;
    }
  } else {
    return 0;
  }
}



// Gateways to Python
PyObject *point_on_line(PyObject *self, PyObject *args) {
  //
  // point_on_line(x, y, x0, y0, x1, y1)
  //
  
  double x, y, x0, y0, x1, y1;
  int res;
  PyObject *result;
  
  // Convert Python arguments to C  
  if (!PyArg_ParseTuple(args, "dddddd", &x, &y, &x0, &y0, &x1, &y1))  
    return NULL;

  
  // Call underlying routine
  res = _point_on_line(x, y, x0, y0, x1, y1);

  // Return values a and b
  result = Py_BuildValue("i", res);
  return result;
}


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, q0, q1, q2, a, b;
  PyObject *result;
  
  // Convert Python arguments to C  
  if (!PyArg_ParseTuple(args, "ddddddddd", &x0, &y0, &x1, &y1, &x2, &y2,
			&q0, &q1, &q2))  
    return NULL;

  
  // Call underlying routine
  _gradient(x0, y0, x1, y1, x2, y2, 
	    q0, q1, q2, &a, &b);

  // Return values a and b
  result = Py_BuildValue("dd", a, b);
  return result;
}


// 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"},    
  {"point_on_line", point_on_line, 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  
}