/*
 * Alessandro MIRONE
 * ESRF
 */

/*
  #########################
  ##Disclaimer
  ##==========
  ##
  ## This  software is provided without warranty of any kind.
  ## No liability is taken for any loss or damages, direct or
  ## indirect, that may result through the use of it. No warranty
  ## is made with respect to this documentation, or the programs
  ## and functions therein.
  ## There are no warranties that the programs or their documentation
  ## are free of error, or that it is consistent with any standard,
  ## or that it will meet the requirement for a particular application.
  
  ## Copyright
  ## =========
  ##
  ## We have adopted the GNU LIBRARY GENERAL PUBLIC LICENSE to apply to
  ## the this software.
  ## For more information on the license terms, see
  ## http://www.gnu.org/copyleft/lgpl.txt
  
*/


#include"Python.h"
#include "structmember.h"
#include<stdio.h>
#include<string.h>
#include "Numeric/arrayobject.h"
#include <math.h>

void KK_M_CINTERFACE( int n_dati, double   *dati, double *x, double EMIN, double EMAX, double  *emin, double  *emax,
		      int *imin, int *imax, double zmol);

void KK_CINTERFACE( int n_dati, double   *dati, double *x,  double * trasf , double  EMIN, double  EMAX,
		    double  *emin, double  *emax , int *imin , int *imax, double zmol);


static PyObject *ErrorObject;
#define onError(message)\
  { PyErr_SetString(ErrorObject, message); return NULL;}

static char KK_doc[] = ""
"/*  5 arguments are taken by the function :\n"
" *  data,        x,         EMIN,    EMAX, ZMOL\n"
" *    |          |           |         |    \n"
" *    |          |           |         |     \n"
" *    |          |           |         |                                 \n"
" *    |          |           |         ---> Float : indicates the range where transformation is done\n"
" *    |          |           |                                 \n"
" *    |          |           --> Float : indicates the range where transformation is done\n"
" *    |          |                                  \n"
" *    |          -> A PyArray of double. the energy points\n"
" *    |\n"
" *    -> A PyArray of double. The absorptions\n"
" *                                       --------\n"
" *  data and x are have the same dimension.\n"
" * The values returned are ( transf, emin,emax, imin, imax)\n"
" * Transf is an array of double, it is the transformation of data.\n"
" * It has the same dimension. The useful values of transf are those\n"
" * whose X is between EMIN and EMAX.\n"
" * emin is the lowest energy of useful transf. emax the highest\n"
" * imin and imax are the position in X of these points.\n"
" */\n"
"";


static PyObject *
pyKK_c (PyObject *self, PyObject *args)
{
  PyArrayObject *dataO,        *xO, *trasfO;
  double EMIN,EMAX;
  int imin,imax;
  double emin,emax;
  int dime[1];
  PyObject *res ;
  double zmol;
  
  if(!PyArg_ParseTuple(args,"OOddd:pyKK_c", (PyObject *) &dataO,(PyObject *)&xO,
		       &EMIN,&EMAX, &zmol)
     )
    return NULL;
  
  
  
  /* check the Objects */
  if(!PyArray_Check((PyObject *)dataO ))    onError("not a PyArray, argument 1");
  if(!PyArray_Check((PyObject *)xO))     onError("not a PyArray, argument 2");
  
  /* check the types */
  if( dataO->descr->type_num != PyArray_DOUBLE ) onError(" arg 1 is not double " ) ;
  if( xO->descr->type_num != PyArray_DOUBLE ) onError(" arg 2  is not double " ) ;
  
  /* check the dimensions */
  if( dataO->nd != 1 )
    onError("arg. 1 has not the right number of dimensions");
  if( xO->nd != 1 )
    onError("arg. 2 has not the right number of dimensions");
  
  if( xO->dimensions[0]!=dataO->dimensions[0])
    onError("the two arrays have not the same number of dimensions");
  
  
  /* check for contiguity */
  
  if((xO)->flags %2 == 0) onError(" All arrays have to be contiguous");
  if((dataO)->flags %2 == 0) onError(" All arrays have to be contiguous");
  
  dime[0]= xO->dimensions[0];
  trasfO = (PyArrayObject*)PyArray_FromDimsAndData( 1, dime,  PyArray_DOUBLE,NULL);
  Py_BEGIN_ALLOW_THREADS
  KK_CINTERFACE( xO->dimensions[0] , (double *) dataO->data , (double *) xO->data,  (double *)trasfO->data , 
		 EMIN, EMAX,
		 &emin, &emax,	
		 &imin, &imax, zmol);
  Py_END_ALLOW_THREADS
  
  
  res = PyTuple_New(5);
  PyTuple_SetItem( res , 0 , (PyObject*)trasfO ) ;
  PyTuple_SetItem( res , 1 , PyFloat_FromDouble(emin) ) ;
  PyTuple_SetItem( res , 2 , PyFloat_FromDouble(emax) ) ;
  PyTuple_SetItem( res , 3 , PyInt_FromLong(imin) ) ;
  PyTuple_SetItem( res , 4 , PyInt_FromLong(imax) ) ;

  return res;
  

}



static char KK_M_doc[] = ""
"/*  5 arguments are taken by the function :\n"
" *  Mdata,        x,         EMIN,    EMAX, ZMOL\n"
" *    |          |           |         |    \n"
" *    |          |           |         |     \n"
" *    |          |           |         |                                 \n"
" *    |          |           |         ---> Float : indicates the range where transformation is done\n"
" *    |          |           |                                 \n"
" *    |          |           --> Float : indicates the range where transformation is done\n"
" *    |          |                                  \n"
" *    |          -> A PyArray of double. the energy points\n"
" *    |\n"
" *    -> A PyArray of double. The absorptions\n"
" *                                       --------\n"
" *  Mdata[0] and x are have the same dimension N. Mdata is (N+1)XN\n"
" * The values returned are (  emin,emax, imin, imax). Mdata is filled.\n"
" * \n"
" *  The useful values of transf Mdata are those\n"
" * whose X is between EMIN and EMAX.\n"
" * emin is the lowest energy of useful transf. emax the highest\n"
" * imin and imax are the position in X of these points.\n"
" */\n"
"";


static PyObject *
pyKK_M_c (PyObject *self, PyObject *args)
{
  PyArrayObject *dataO,        *xO;
  double EMIN,EMAX;
  int imin,imax;
  double emin,emax;
  PyObject *res ;
  double zmol;
  
  if(!PyArg_ParseTuple(args,"OOddd:pyKK_M_c", (PyObject *) &dataO,(PyObject *)&xO,
		       &EMIN,&EMAX, &zmol)
     )
    return NULL;
  
  
  
  /* check the Objects */
  if(!PyArray_Check((PyObject *)dataO ))    onError("not a PyArray, argument 1");
  if(!PyArray_Check((PyObject *)xO))     onError("not a PyArray, argument 2");
  
  /* check the types */
  if( dataO->descr->type_num != PyArray_DOUBLE ) onError(" arg 1 is not double " ) ;
  if( xO->descr->type_num != PyArray_DOUBLE ) onError(" arg 2  is not double " ) ;
  
  /* check the dimensions */
  if( dataO->nd != 2 )
    onError("arg. 1 has not the right number of dimensions");
  if( xO->nd != 1 )
    onError("arg. 2 has not the right number of dimensions");
  
  if( xO->dimensions[0]!=dataO->dimensions[1])
    onError("the two arrays have not the same X dimension");
  
  
  /* check for contiguity */
  
  if((xO)->flags %2 == 0) onError(" All arrays have to be contiguous");
  if((dataO)->flags %2 == 0) onError(" All arrays have to be contiguous");
  
  Py_BEGIN_ALLOW_THREADS
  KK_M_CINTERFACE( xO->dimensions[0] , (double *) dataO->data ,(double *) xO->data, 
		 EMIN, EMAX,
		 &emin, &emax,	
		 &imin, &imax, zmol);
  Py_END_ALLOW_THREADS
  
  
  res = PyTuple_New(4);
  PyTuple_SetItem( res , 0 , PyFloat_FromDouble(emin) ) ;
  PyTuple_SetItem( res , 1 , PyFloat_FromDouble(emax) ) ;
  PyTuple_SetItem( res , 2 , PyInt_FromLong(imin) ) ;
  PyTuple_SetItem( res , 3 , PyInt_FromLong(imax) ) ;

  return res;
  

}








static PyMethodDef KKpy_c_functions[] = {
  {"KKpy_c", pyKK_c, METH_VARARGS ,KK_doc },
  {"KKpy_M_c", pyKK_M_c, METH_VARARGS ,KK_M_doc },
  { NULL, NULL}
};


/*  extern "C" { */
void 
initKKpy_c(void );
/*  } */

void 
initKKpy_c()
{
  PyObject *m, *d;
  m = Py_InitModule("KKpy_c", KKpy_c_functions);
  d = PyModule_GetDict(m);
  ErrorObject = Py_BuildValue("s","KKpy_c.error");
  PyDict_SetItemString(d,"error", ErrorObject);
  if(PyErr_Occurred())
    Py_FatalError("can't initialize module KKpy_c ");

#ifdef import_array
  import_array();
#endif
}



















/* X------------------------------------------ C_INTERFACE ---------------------------------------------------------X */




#define max(a,b) ((a) > (b) ) ? (a) : (b) 


#ifndef PI
#define PI M_PI
#endif

double KK(double *x, double * dati, int n_dati, int i);

void KK_CINTERFACE( int n_dati, double   *dati, double *x,  double * trasf , double EMIN, double EMAX, double  *emin, double  *emax,
		    int *imin, int *imax, double zmol)
{
  
/*    int swl, swh; */
/*    double sw; */


  /* la formula che si trov in Gullikson per esempio est
       lambda**2 * Re /2pi * densita * F

     dove la densita e il numero di atomi per unita di volume
     F lo scattering factor di un atomo tipo

     Quindi Re e' in unita atomiche 1/137/137 Bhor ( 0.529 Angstrom)

     La formula qui sotto e' usa angstroms angstrom e 0.6022 e' la densita di 
     atomi per Angstrom cubo per una mole per centimetro quadro

     Cioe Avogadro = 10**(8*3) *0.6022 = 6.022 * 10**(23)

  */
 

  double zfatt= 12398.52*12398.52*  0.529/137.036/137.036 *0.6022/2/PI ;
  int i;
  
/*    // ************************************************************ */
/*    // ** si assume che le energie siano decrescenti */
/*    // ** */
  
/*    if (x[0] < x[1]) // trasfomazione in energie discendenti */
/*      { */
/*        cout << " TRASFORMO IN ENERGIE DECRESCENTI \n"; */
/*        swl=0; */
/*        swh=n_dati-1; */
      
/*        while (swl<swh) */
/*  	{ */
/*  	  sw=x[swl];  x[swl]=x[swh]; x[swh]=sw; */
/*  	  sw=dati[swl];  dati[swl]=dati[swh]; dati[swh]=sw; */
/*  	  swl++; */
/*  	  swh--; */
/*  	} */
/*      } */
  
  printf("# eseguo la trasformazione di KK \n" );
     printf("   EMIN   %e  \n" , EMIN);
     printf("   EMAX   %e  \n" , EMAX);
  
     printf("  x[n_dati-1] %e \n" , x[n_dati-1] );
     printf("  x[0       ] %e \n" , x[0]        );
  
  if( x[n_dati-1]<EMAX || x[0]>EMIN)
    {
      printf("# Attenzione: occore che  EMIN > %e e che  EMAX < %e \n ", x[0],  x[n_dati-1] ); 
      exit(0);
    }
  
  //***        conv. eV2Ang**-1---el. class. rayon -----Avog
  *emin=0;
  *emax=0;
  for( i=0; i< n_dati; i++)
    {
      if( x[i]>=EMIN && x[i]<=EMAX)
        {

	  trasf[i]  = KK(x, dati , n_dati, i) ;
	  trasf[i] -= zfatt*zmol /x[i]/x[i] ;

        }
      
      if( x[i]    > EMIN  && *emin==0)   
	{               // si stabilisce la baseline

	  *emin= x[i];
	  *imin=i;

	}
      if( x[i] > EMAX  && *emax==0)
	{                 // si stabilisce la baseline

	  *emax= x[i-1];
	  *imax=i-1;

	}
    }
}


double KK(double *x, double * dati, int n_dati, int pos)
{
  
  double res=0;
  double x1,x2,y1,y2,x0,a,b;
  int i;
  
  x0= x[pos]*x[pos];
  
  for(i=0 ; i< n_dati-1; i++)
    {
      
      x1= x[i]*x[i];
      x2= x[i+1]*x[i+1];
      

      if(x1==x2) continue;
      y1= dati[i]   ;
      y2= dati[i+1]  ;
      
      b=(y2-y1)/(x2-x1)/2.;
      a=y1/2.-b*x1;
      
      if(i==(pos-1) || x2==x0 )
        {
	  res += b*(x2-x1) + (a+b*x0)*log(fabs(1./(x1-x0)));
        }
      else if(i==pos || x1==x0  )
	{
	  res += b*(x2-x1) + (a+b*x0)*log(fabs((x2-x0))); 
	}
      else if(i> pos || i < (pos-1) )
	{
	  res += b*(x2-x1) + (a+b*x0)*log(fabs((x2-x0)/(x1-x0))); 
	}

      if(!(res> -100000000000.)) {
	printf(" rpoblema con un passo di integrazione %e %e \n" ,x1,x2);
	exit(0);
      }
    }
  return res*2./PI  ;
}




void  KK_M(double *x, double * dati, int n_dati, int i);

void KK_M_CINTERFACE( int n_dati, double   *dati, double *x, double EMIN, double EMAX, double  *emin, double  *emax,
		    int *imin, int *imax, double zmol)
{
  
/*    int swl, swh; */
/*    double sw; */
  double zfatt= 12398.52*12398.52*  0.529/137.036/137.036 *0.6022/2/PI ;
  int i;
  
/*    // ************************************************************ */
/*    // ** si assume che le energie siano decrescenti */
/*    // ** */
  
/*    if (x[0] < x[1]) // trasfomazione in energie discendenti */
/*      { */
/*        cout << " TRASFORMO IN ENERGIE DECRESCENTI \n"; */
/*        swl=0; */
/*        swh=n_dati-1; */
      
/*        while (swl<swh) */
/*  	{ */
/*  	  sw=x[swl];  x[swl]=x[swh]; x[swh]=sw; */
/*  	  sw=dati[swl];  dati[swl]=dati[swh]; dati[swh]=sw; */
/*  	  swl++; */
/*  	  swh--; */
/*  	} */
/*      } */
  
  printf("# eseguo la trasformazione di KK \n" );
  printf("   EMIN   %e  \n" , EMIN);
  printf("   EMAX   %e  \n" , EMAX);
  
  printf("  x[n_dati-1] %e \n" , x[n_dati-1] );
  printf("  x[0       ] %e \n" , x[0]        );
  
  if( x[n_dati-1]<EMAX || x[0]>EMIN)
    {
      printf("# Attenzione: occore che  EMIN > %e e che  EMAX < %e \n ", x[0],  x[n_dati-1] ); 
      exit(0);
    }
  
  //***        conv. eV2Ang**-1---el. class. rayon -----Avog
  *emin=0;
  *emax=0;
  for( i=0; i< n_dati; i++)
    {
      if( x[i]>=EMIN && x[i]<=EMAX)
        {

	  KK_M(dati+ i*n_dati, x, n_dati, i) ;
	  dati[ n_dati*n_dati + i ] = -zfatt*zmol /x[i]/x[i] ;

        }
      
      if( x[i]    > EMIN  && *emin==0)   
	{               // si stabilisce la baseline

	  *emin= x[i];
	  *imin=i;

	}
      if( x[i] > EMAX  && *emax==0)
	{                 // si stabilisce la baseline

	  *emax= x[i-1];
	  *imax=i-1;

	}
    }
}


void  KK_M( double * dati , double *x, int n_dati, int pos)
{
  
  double x1,x2,y1,y2,x0,a,b;
  int i,j;
  
  x0= x[pos]*x[pos];
  dati[0]=dati[n_dati-1]=0;
  for( j=1 ; j< n_dati-1; j++)
    {
      dati[j]=0;
      for( i=j-1; i<j+1; i++) {
	x1= x[i]*x[i];
	x2= x[i+1]*x[i+1];
        if(i==j) {
	  y1= 1.0  ;
	  y2= 0.0  ;
	} else {
	  y1= 0.0  ;
	  y2= 1.0  ;

	}
	
	b=(y2-y1)/(x2-x1)/2.;
	a=y1/2.-b*x1;
	
	if(i> pos || i < (pos-1) )
	  {
	    dati[j] += b*(x2-x1) + (a+b*x0)*log(fabs((x2-x0)/(x1-x0))); 
	  }
	else if(i==(pos-1) )
	  {
	    dati[j] += b*(x2-x1) + (a+b*x0)*log(fabs(1./(x1-x0)));
	  }
	else if(i==pos)
	  {
	    dati[j] += b*(x2-x1) + (a+b*x0)*log(fabs((x2-x0))); 
	  }
      
	if(!(dati[j]> -100000000000.)) exit(0);
      }
      dati[j] *=2.0/PI;

    }
}












