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


double max(double x, double y) {  
  //Return maximum of two doubles
  
  if (x > y) return x;
  else return y;
}


double min(double x, double y) {  
  //Return minimum of two doubles
  
  if (x < y) return x;
  else return y;
}


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

  double det;
  
  det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0);

  *a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0);
  *a /= det;

  *b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0);
  *b /= det;

  return 0;
}


void _limit(int N, double beta, double* qc, double* qv, 
	    double* qmin, double* qmax) { 

  //N are the number of elements
  int k, i, k3;
  double dq, dqa[3], phi, r;
  
  //printf("INSIDE\n");
  for (k=0; k<N; k++) {
    k3 = k*3;
    
    //Find the gradient limiter (phi) across vertices  
    phi = 1.0;
    for (i=0; i<3; i++) {    
      r = 1.0;
      
      dq = qv[k3+i] - qc[k];    //Delta between vertex and centroid values
      dqa[i] = dq;              //Save dq for use in the next loop
      
      if (dq > 0.0) r = (qmax[k] - qc[k])/dq;
      if (dq < 0.0) r = (qmin[k] - qc[k])/dq;      
  
  
      phi = min( min(r*beta, 1.0), phi);    
    }
    
    //Then update using phi limiter
    for (i=0; i<3; i++) {    
      qv[k3+i] = qc[k] + phi*dqa[i];
    }
  }
}


void print_numeric_array(PyArrayObject *x) {  
  int i, j;
  for (i=0; i<x->dimensions[0]; i++) {  
    for (j=0; j<x->dimensions[1]; j++) {
      printf("%f ", *(double*) (x->data + i*x->strides[0] + j*x->strides[1]));
    }
    printf("\n");  
  }
  printf("\n");    
}

void print_numeric_vector(PyArrayObject *x) {  
  int i;
  for (i=0; i<x->dimensions[0]; i++) {
    printf("%f ", *(double*) (x->data + i*x->strides[0]));  
  }
  printf("\n");  
}

PyArrayObject *get_consecutive_array(PyObject *O, char *name) {
  PyArrayObject *A, *B;
  

  //Get array object from attribute
  
  /*
  //FIXME: THE TEST DOESN't WORK
  printf("Err = %d\n", PyObject_HasAttrString(O, name));
  if (PyObject_HasAttrString(O, name) == 1) {
    B = (PyArrayObject*) PyObject_GetAttrString(O, name);
    if (!B) return NULL;
  } else {
    return NULL;
    } 
  */
    
  B = (PyArrayObject*) PyObject_GetAttrString(O, name);
  if (!B) return NULL;     
  
  //Convert to consecutive array
  A = (PyArrayObject*) PyArray_ContiguousFromObject((PyObject*) B, 
						    B -> descr -> type, 0, 0);
  
  Py_DECREF(B); //FIXME: Is this really needed??
  
  if (!A) return NULL;
  return A;
}