/*
 * Alessandro MIRONE
 * April the 10th 2001, 
 * ESRF
 */

#include"Python.h"
#include "structmember.h"
#include<stdio.h>
#include<iostream.h>
#include<string.h>

#include "Numeric/arrayobject.h"
#include <math.h>
#include<complex.h>

#define DEBUG(a) printf(a); printf("\n");

#define swap(a,b,c) c=a;a=b;b=c;

/*
 * The error object to expose 
 */

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

static char multicouche_calculatescan_doc[] = "\n"
"/*  7 arguments are taken by the function :\n"
" *  thickness, roughness, wavelenght, angle, indexes, Sindex, Srough\n"
" *    |          |           |         |      |        |       |\n"
" *    |          |           |         |      |        |       |\n"
" *    |          |           |         |      |        |       |\n"
" *    |          |           |         |      |        |       =>PyDoubleObject\n"
" *    |          |           |         |      |        ->PyArray_Complex. Dim= numwave\n"
" *    |          |           |         |      --> PyArray_Complex. Dims= numthickXnumwave\n"
" *    |          |           |         |                                 ----------------\n"
" *    |          |           |         --->PyArray_double. Dim= numwave\n"
" *    |          |           |                                  ------\n"
" *    |          |           --> PyArray_double. Dim= numwave\n"
" *    |          |                                   -------\n"
" *    |          -> A PyArray of double. Same dimensions as thickness\n"
" *    |\n"
" *    -> A PyArray of double. Dimension= numthick. Angstroms\n"
" *                                       --------\n"
" * Thickness are given starting from the bottom (first layer over the substrate)\n"
" * indexes for the substrate are given separately in Sindex\n"
" * Srough is the substrate roughness.\n"
" *\n"
" * RETURN VALUE\n"
" *    A tuple containing the S reflectivity an Preflectivity given\n"
" *    as PyArray double. \n"
" */\n"
"";

static PyObject *
multicouche_calculatescan(PyObject *self, PyObject *args)
{
  FILE *f;
  char *name;
  PyArrayObject *thickness, *roughness, *wavelenght, *angle, *indexes, *Sindex;
  PyFloatObject *Srough;
  
  int i,j, lenght;
  char * endian;
  
  
  if(!PyArg_ParseTuple(args,"OOOOOOO:multicouche_calculatescan", (PyObject *) &thickness,(PyObject *)&roughness,
		       (PyObject *)&wavelenght,
		       (PyObject *)&angle,(PyObject *) &indexes, 
		       (PyObject *)&Sindex,(PyObject *)&Srough  )
     )
    return NULL;
  
  /* check the Objects */
  if(!PyArray_Check((PyObject *)thickness ))    onError("not a PyArray, argument 1");
  if(!PyArray_Check((PyObject *)roughness))     onError("not a PyArray, argument 2");
  if(!PyArray_Check((PyObject *)wavelenght))    onError("not a PyArray, argument 3");
  if(!PyArray_Check((PyObject *)angle))         onError("not a PyArray, argument 4");
  if(!PyArray_Check((PyObject *)indexes))       onError("not a PyArray, argument 5");
  if(!PyArray_Check((PyObject *)Sindex))        onError("not a PyArray, argument 6");
  if(!PyFloat_Check((PyObject *)Srough))        onError("not a PyFloat, argument 7");
  
  /* check the types */
  if( thickness->descr->type_num != PyArray_DOUBLE ) onError(" thickness is not double " ) ;
  if( roughness->descr->type_num != PyArray_DOUBLE ) onError(" roughness  is not double " ) ;
  if( wavelenght->descr->type_num != PyArray_DOUBLE ) onError(" wavelenght is not double " ) ;
  if( angle->descr->type_num != PyArray_DOUBLE ) onError(" angle is not double " ) ;
  if( indexes->descr->type_num != PyArray_CDOUBLE ) onError(" indexes is not cdouble " ) ;
  if( Sindex->descr->type_num != PyArray_CDOUBLE ) onError(" Sindex is not double " ) ;
  
  /* check the dimensions */
  if( thickness->nd != 1 )
    onError("The thickness array (arg. 1) has not the right number of dimensions");
  if(roughness ->nd != 1 )
    onError("The roughness array (arg. 2) has not the right number of dimensions");
  if( wavelenght->nd != 1 )
    onError("The wavelenght array (arg. 3) has not the right number of dimensions");
  if(angle->nd != 1 )
    onError("The angle array (arg. 4) has not the right number of dimensions");
  if(indexes->nd != 2 )
    onError("The indexes array (arg. 5) has not the right number of dimensions");
  if(Sindex->nd != 1 )
    onError("The Sindex array (arg. 6) has not the right number of dimensions");
  
  
  int numthick, numwave;
  numthick=thickness->dimensions[0];
  if(numthick!=roughness->dimensions[0] )
    onError("The roughness array (arg. 2) has not the right  dimension");
  
  numwave  = wavelenght->dimensions[0] ;
  
  if(numwave!=angle->dimensions[0] )
    onError("The  angles array (arg. 4) has not the right  dimension");
  
  if(numthick!=indexes->dimensions[0] )
    onError("The indexes  array (arg. 5) has not the right  first dimension");
  if(numwave!=indexes->dimensions[1] ) {
    printf("%d %d\n", numwave, indexes->dimensions[1]);
    onError("The indexes  array (arg. 5) has not the right second   dimension");
  }
  if(numwave!=Sindex->dimensions[0] )
    onError("The Sindex  array (arg. 6) has not the right  second dimension");
  
  

  /* check that everything is contiguous */
  
  PyArrayObject *toverify[]={thickness,roughness,wavelenght,angle,indexes,Sindex,NULL};
  PyArrayObject **ptr;
  ptr=toverify;
  while(*ptr) {
    if((*ptr)->flags %2 == 0) onError(" All arrays have to be contiguous");
    ptr++;
  }
  
  PyObject *res = PyTuple_New(2);


//            PyThreadState *_save;

//            _save = PyThreadState_Swap(NULL);
//            PyEval_ReleaseLock();
//  	  cout << "  Py_BEGIN_ALLOW_THREADS " << endl;

Py_BEGIN_ALLOW_THREADS

  /* creation of the return Values */
  
  PyObject *res1,*res2, *tuple;
  int d[1];
  d[0]= numwave;
  res1= (PyObject*)PyArray_FromDimsAndData( 1, d,  PyArray_CDOUBLE,NULL);
  res2= (PyObject*)PyArray_FromDimsAndData( 1, d,  PyArray_CDOUBLE,NULL);
  complex<double> *ReflS =(complex<double> *) (((PyArrayObject*) res1)->data) ;
  complex<double> *ReflP =(complex<double> *)  (((PyArrayObject*) res2)->data) ;
  
  
  /* allocation of arrays */
  
  complex<double> *dumRS = new complex<double>[numwave];
  complex<double> *dumRP = new complex<double>[numwave];
  complex<double> *dumRS_ = new complex<double>[numwave];
  complex<double> *dumRP_ = new complex<double>[numwave];
  complex<double> *propag = new complex<double>[numwave];
  complex<double> *eps    = new complex<double>[numwave] ;
  complex<double> *epsnext    = new complex<double>[numwave] ;
  complex<double> *dumkC  = new complex<double>[numwave]; /* auxiliary array complex */
  complex<double> *dumkCnext  = new complex<double>[numwave]; /* auxiliary array complex */
  complex<double> *dumC1  = new complex<double>[numwave]; /* auxiliary array complex */
  for(int k=0; k<  numwave; k++) dumC1[k]=complex<double>(1.0,0.0);

  double *kpar  = new double  [numwave]; /* auxiliary array double */;
  double *kperp = new double  [numwave]; /* auxiliary array double */;
  double *K0    = new double  [numwave]; /* auxiliary array double */;
  // double *fat   = new double  [numwave]; /* auxiliary array double */;
  complex<double> *fat   = new  complex<double> [numwave]; /* auxiliary array double */;
  complex<double> caux;

  /* useful  pointers */
  complex<double> *index, *indexnext;
  
  /* useful variable */ 
  double thick, rough;
  

  /* useful initialisations */
  for(int i=0; i<numwave; i++) {
    K0[i]= 2*M_PI/( (double*)wavelenght->data ) [i] ;
  }
  for(int i=0; i<numwave; i++) {
    kpar[i]= K0[i]*cos(( (double*)angle->data ) [i]);
  }
  for(int i=0; i<numwave; i++) {
    kperp[i]= K0[i]*sin(( (double*)angle->data ) [i]);
  }
  
  
  /* initialisation of reflectivity for the substrate */
  index     = (complex<double> *) Sindex->data;
  if(numthick>0) {
    indexnext = ((complex<double> *) (indexes->data) ) +0*numwave;
  } else {
    indexnext = dumC1;
  }
  rough = PyFloat_AsDouble( (PyObject*) (Srough) );
  
  for(int i=0; i<numwave; i++) {
    eps[i] = (index[i])*(index[i]);
    epsnext[i] = (indexnext[i])*(indexnext[i]);
    dumkC[i]=sqrt((index[i])*(index[i])*K0[i]*K0[i] -kpar[i]*kpar[i]);
    dumkCnext[i]=sqrt((indexnext[i])*(indexnext[i])*K0[i]*K0[i] -kpar[i]*kpar[i]);
  }
  
  for(int i=0; i<numwave; i++) {
    // fat[i] = exp(-rough*rough*kperp[i]*kperp[i]*2  );
    caux = 2.0*dumkC[i]*dumkCnext[i] ; 
    fat[i] = exp(-rough*rough*caux );
  }
  
  for(int i=0; i<numwave; i++) {
    ReflS[i]=(dumkCnext[i]-dumkC[i])/(dumkCnext[i]+dumkC[i])*fat[i];
    ReflP[i]=(dumkCnext[i]*eps[i]-dumkC[i]*epsnext[i] )/(kperp[i]*eps[i]+dumkC[i]*epsnext[i] )*fat[i];
  }
  
  if(1) {
    /* And now for each layer */
    for(int l=0; l<numthick; l++) {
      index = ((complex<double> *) (indexes->data) ) +l*numwave;
      if(l==numthick-1) {
        indexnext = dumC1;
      } else {
	indexnext = ((complex<double> *) (indexes->data) ) +(l+1)*numwave;
      }
      rough = *(((double *) (roughness->data) )+l)  ;
      thick = *(((double *) (thickness->data) )+l)  ;
      
      
      complex<double> * dumc;
      
      swap(eps,epsnext,dumc);
      swap(dumkC,dumkCnext,dumc);
      
      if(l!=numthick-1) {
	for(int i=0; i<numwave; i++) {
	  epsnext[i] = (indexnext[i])*(indexnext[i]);
	  dumkCnext[i]=sqrt((indexnext[i])*(indexnext[i])*K0[i]*K0[i] -kpar[i]*kpar[i]);
	}
      } else {
	for(int i=0; i<numwave; i++) {
	  dumkCnext[i]=kperp[i];
	}
      }
      
      // cout<< " l " << l  << " numthick " <<numthick << " index " << index[0]<<endl;

      for(int i=0; i<numwave; i++) {
	propag[i] = exp( +2*thick*complex<double>(0.0,1.0)*dumkC[i]     );
      }
      
      for(int i=0; i<numwave; i++) {
	// fat[i] = exp(-rough*rough*kperp[i]*kperp[i]*2  );
	caux = 2.0*dumkC[i]*dumkCnext[i] ; 
	fat[i] = exp(-rough*rough*caux );
      }
      
      if(l!=numthick-1) {
	for(int i=0; i<numwave; i++) {
	  dumRS_[i]= (dumkCnext[i]-dumkC[i])/(dumkCnext[i]+dumkC[i]);
	  // dumRP_[i]= (dumkCnext[i]*eps[i]-dumkC[i]*epsnext[i] )/(dumkCnext[i]*eps[i]+dumkC[i]*epsnext[i] );
	  dumRS[i] = dumRS_[i]*fat[i];
	  // dumRP[i] = dumRP_[i]*fat[i];
	}   
      } else {
	for(int i=0; i<numwave; i++) {
	  dumRS_[i]=(kperp[i]-dumkC[i])/(kperp[i]+dumkC[i]);
	  // dumRP_[i]=(kperp[i]*eps[i]-dumkC[i] )/(kperp[i]*eps[i]+dumkC[i] );
	  dumRS[i] = dumRS_[i]*fat[i];
	  // dumRP[i] = dumRP_[i]*fat[i];
	}   
      }
      
      for(int i=0; i<numwave; i++) {
	ReflS[i]=(dumRS[i]+propag[i]*ReflS[i]*(1.0-dumRS_[i]*dumRS_[i] +dumRS[i]*dumRS[i]) )/(1.0+dumRS[i]*propag[i]*ReflS[i]);
	// ReflP[i]=(dumRP[i]+propag[i]*ReflP[i]*(1-dumRP_[i]*dumRP_[i] +dumRP[i]*dumRP[i]) )/(1+dumRP[i]*propag[i]*ReflP[i]);
      }     
      
      
    }
  }
  /* DEallocation of arrays */
  
  delete [] dumRS ;
  delete [] dumRP ;
  delete [] dumRS_ ;
  delete [] dumRP_ ;
  delete [] propag ;
  delete [] eps    ;
  delete [] epsnext    ;
  delete [] dumkC  ; 
  delete [] dumkCnext  ; 
  delete [] dumC1  ; 
  
  delete [] kpar  ; 
  delete [] kperp ; 
  delete [] K0    ; 
  delete [] fat   ; 
  
  PyTuple_SetItem( res , 0 , res1 ) ;
  PyTuple_SetItem( res , 1 , res2 ) ;

Py_END_ALLOW_THREADS
//    PyEval_AcquireLock();
//    PyThreadState_Swap(_save);


  return res;
}


static PyMethodDef multicouche_functions[] = {
  {"PPM_calculatescan", multicouche_calculatescan, METH_VARARGS , multicouche_calculatescan_doc },
  { NULL, NULL}
};
  
  
  extern "C" {
    void 
initPPMcore();
}
void 
initPPMcore()
{
  PyObject *m, *d;
  m = Py_InitModule("PPMcore", multicouche_functions);
  d = PyModule_GetDict(m);
  ErrorObject = Py_BuildValue("s","PPMcore.error");
  PyDict_SetItemString(d,"error", ErrorObject);
  if(PyErr_Occurred())
    Py_FatalError("can't initialize module PPMcore");

#ifdef import_array
  import_array();
#endif
}





















