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

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

#define DDATA(p) ((double *) (((PyArrayObject *)(p))->data))
#define ZDATA(p) ((complex<double> *) (((PyArrayObject *)(p))->data))
 
#define swap(a,b,c) c=a;a=b;b=c;


#define complex std::complex



#define Complex complex<double>
#define fComplex complex<float>

extern "C"  // some Lapack declarations ( just a fraction of this is really used ...)
{
void sgetrf_( int *m, int *n, float *A, int *lda, int *ipiv, int*info );
void sgetri_( int *n, float *A, int *lda, int *ipiv, float* work, int *lwork,
              int *info );
void dgetrf_( int *m, int *n, double *A, int *lda, int *ipiv, int*info );
void dgetri_( int *n, double *A, int *lda, int *ipiv, double* work, int *lwork,
              int *info );
void cgetrf_( int *m, int *n, fComplex *A, int *lda, int *ipiv, 
	      int *info );
void cgetri_( int *n, fComplex *A, int *lda, int *ipiv, 
	      fComplex *work, int *lwork, int *info );
void zgetrf_( int *m, int *n, Complex *A, int *lda, int *ipiv, 
	      int *info );
void zgetri_( int *n, Complex *A, int *lda, int *ipiv, Complex *work, 
	      int *lwork, int *info );

void dgeev_( char * jobvl, char* jobvr, int *N, double *A, int *lda,
	     double *wr, double *wi, double *vl, int *ldvl, double *vr,
	     int *ldvr, double *work, int *lwork, int *info );
void dsyev_( char *jobz, char *uplo, int *N, double *a, int *lda, double *w,
             double *work, int *lwork, int *info );

void sgeev_( char * jobvl, char* jobvr, int *N, float *A, int *lda,
	     float *wr, float *wi, float *vl, int *ldvl, float *vr,
	     int *ldvr, float *work, int *lwork, int *info );
void ssyev_( char *jobz, char *uplo, int *N, float *a, int *lda, float *w,
             float *woArk, int *lwork, int *info );

void cgeev_( char * jobvl, char* jobvr, int *N, fComplex *A, int *lda,
	     fComplex *w, fComplex *vl, int *ldvl, fComplex *vr,
	     int *ldvr, fComplex *work, int *lwork, float *rwork, 
	     int *info );
void cgeevx_( char *balanc, char *jobvl, char *jobvr, char *sense, int
	      *N, fComplex *A, int *lda, fComplex *w, fComplex *vl,
	      int *ldvl, fComplex *vr, int *ldvr, int *ilo, int *ihi,
	      float *scale, float *abnrm, float *rconde, float
	      *rcondv, fComplex *work, int *lwork, float *rwork,  
	     int *info );
void sgeevx_( char *balanc, char *jobvl, char *jobvr, char *sense, int
	      *N, float *A, int *lda, float *wr, float *wi, float *vl,
	      int *ldvl, float *vr, int *ldvr, int *ilo, int *ihi,
	      float *scale, float *abnrm, float *rconde, float
	      *rcondv, float *work, int *lwork, int *iwork,  
	     int *info );

void zgeevx_( char *balanc, char *jobvl, char *jobvr, char *sense, int
	      *N, Complex *A, int *lda, Complex *w, Complex *vl,
	      int *ldvl, Complex *vr, int *ldvr, int *ilo, int *ihi,
	      double *scale, double *abnrm, double *rconde, double
	      *rcondv, Complex *work, int *lwork, double *rwork,  
	     int *info );
void dgeevx_( char *balanc, char *jobvl, char *jobvr, char *sense, int
	      *N, double *A, int *lda, double *wr, double *wi, double *vl,
	      int *ldvl, double *vr, int *ldvr, int *ilo, int *ihi,
	      double *scale, double *abnrm, double *rconde, double
	      *rcondv, double *work, int *lwork, int *iwork,  
	     int *info );

void cheev_( char *jobz, char *uplo, int *N, fComplex *a, int *lda, 
	     float *w, fComplex *work, int *lwork, float *rwork,
	     int *info );
void zheev_( char *jobz, char *uplo, int *N, Complex *a, int *lda, 
	     double *w, Complex *work, int *lwork, double *rwork,
	     int *info );

void sgelss_( int *m, int *n, int *nrhs, float *A, int *lda, float *b, 
	    int *ldb, float *s, float *rcond, int *rank, float *work,
	    int *lwork, int *info );

void cgelss_( int *m, int *n, int *nrhs, fComplex *A, int *lda, fComplex *b,
	    int *ldb, float *s, float *rcond, int *rank, fComplex *work,
	    int *lwork, float *rwork, int *info );

void dgelss_( int *m, int *n, int *nrhs, double *A, int *lda, double *b, 
	    int *ldb, double *s, double *rcond, int *rank, double *work, 
	    int *lwork, int *info );

void zgelss_( int *m, int *n, int *nrhs, Complex *A, int *lda, Complex *b, 
	    int *ldb, double *s, double *rcond, int *rank, Complex *work, 
	    int *lwork, double *rwork, int *info );

void sgesvd_(char *jobu, char *jobvt, int *m,int *n,  float *A, int *lda,
	     float *s, float * U, int *ldu, float *VT, int *ldvt, float *work, 
	     int *lwork, int *info);
void cgesvd_(char *jobu, char *jobvt, int *m,int *n,  fComplex *A, int *lda,
	     float *s, fComplex * U, int *ldu, fComplex *VT, int *ldvt,
	     fComplex *work, int *lwork, float *rwork, int *info);
void dgesvd_(char *jobu, char *jobvt, int *m,int *n,  double *A, int *lda,
	     double *s, double * U, int *ldu, double *VT, int *ldvt, 
	     double *work, int *lwork, int *info);
void zgesvd_(char *jobu, char *jobvt, int *m,int *n,  Complex *A, int *lda,
	     double *s, Complex * U, int *ldu, Complex *VT, int *ldvt,
	     Complex *work, int *lwork, double *rwork, int *info);

// BLAS stuff we use
void sgemm_( char *transa, char *transb, int *m, int *n, int *k, 
	     float *alpha, float *a, int *Ida, float *b, int *Idb, 
	     float *beta, float* c, int *Idc );
void dgemm_( char *transa, char *transb, int *m, int *n, int *k, 
	     double *alpha, double *a, int *Ida, double *b, int *Idb, 
	     double *beta, double* c, int *Idc );
void cgemm_( char *transa, char *transb, int *m, int *n, int *k, 
	     fComplex *alpha, fComplex *a, int *Ida, fComplex *b, int *Idb, 
	     fComplex *beta, fComplex* c, int *Idc );
void zgemm_( char *transa, char *transb, int *m, int *n, int *k, 
	     Complex *alpha, Complex *a, int *Ida, Complex *b, int *Idb, 
	     Complex *beta, Complex* c, int *Idc );


void saxpy_(int *n, float *a, float *x, int *inx, float *y, int *incy);
void daxpy_(int *n, double *a, double  *x, int *inx, double  *y, int *incy);
void caxpy_(int *n, fComplex *a, fComplex *x, int *inx,fComplex *y, int *incy);
void zaxpy_(int *n,Complex *a,Complex *x, int *inx,Complex *y, int *incy);

void dcopy_(int *n,double *dx,int *incx,double * dy, int * incy); 
void scopy_(int *n,float *dx,int *incx,float * dy, int * incy); 
 
int zgeev_(char  * ,char  *,int *,complex<double> * , int *,complex<double> * ,complex<double> * ,int *,
	     complex<double> * , int *,  complex<double> * ,int *,double * ,int *);

}










/*
 * 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)
{


  PyArrayObject *thickness, *roughness, *wavelenght, *angle, *indexes, *Sindex;
  PyFloatObject *Srough;
  

  
  
  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++;
  }



  /* creation of the return Values */
  
  PyObject *res1,*res2;
  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> *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 */;
  
  /* 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  );
  }
  
  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]);
	}
      }
      
      // 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  );
     }
      
      if(l!=numthick-1) {
	for(int i=0; i<numwave; i++) {
	  dumRS[i]=(dumkCnext[i]-dumkC[i])/(dumkCnext[i]+dumkC[i])*fat[i];
	  dumRP[i]=(dumkCnext[i]*eps[i]-dumkC[i]*epsnext[i] )/(dumkCnext[i]*eps[i]+dumkC[i]*epsnext[i] )*fat[i];
	}   
      } else {
	for(int i=0; i<numwave; i++) {
	  dumRS[i]=(kperp[i]-dumkC[i])/(kperp[i]+dumkC[i])*fat[i];
	  dumRP[i]=(kperp[i]*eps[i]-dumkC[i] )/(kperp[i]*eps[i]+dumkC[i] )*fat[i];
	}   
      }
      
      for(int i=0; i<numwave; i++) {
	ReflS[i]=(dumRS[i]+propag[i]*ReflS[i])/(1.0+dumRS[i]*propag[i]*ReflS[i]);
	ReflP[i]=(dumRP[i]+propag[i]*ReflP[i])/(1.0+dumRP[i]*propag[i]*ReflP[i]);
      }     
      
      
    }
  }
  /* DEallocation of arrays */
  
  delete [] dumRS ;
  delete [] dumRP ;
  delete [] propag ;
  delete [] eps    ;
  delete [] epsnext    ;
  delete [] dumkC  ; 
  delete [] dumkCnext  ; 
  delete [] dumC1  ; 
  
  delete [] kpar  ; 
  delete [] kperp ; 
  delete [] K0    ; 
  delete [] fat   ; 
  
  PyObject *res = PyTuple_New(2);
  PyTuple_SetItem( res , 0 , res1 ) ;
  PyTuple_SetItem( res , 1 , res2 ) ;
  
  return res;
}


//***************************************************
//** QUESTA ROUTINE RIORDINA IN DENSO DECRESCENTE 
//** NELLA PARTE REALE UN ARRAY DI COMPLESSI
//**
void riordina (complex<double> *val ,  int n)
  {
     int i,j;
     complex<double> tmp;
     for(i=0; i<n-1; i++ )
       for ( j=i+1; j<n ; j++ )
         {
            if ( val[i].real()<val[j].real() )
                {
                  tmp = val[i] ; 
                  val[i] = val[j] ;
                  val[j] = tmp ;
                }
         }
  }
//**
//**************************************************

class Onda {
public:  
  Onda();
  ~Onda();
  void CreateModes( );
  void CreateModesScalar( );
  void CreateModesTensor( );
  void CreateModesTensor_Vide( );
  
  void mostra(int sp);
  PyObject * get_vectors(); 

  PyArrayObject* indexes;
  PyArrayObject* angles;


  complex<double> *Kx;
  complex<double> *Kz;     // componenti z dei vettori d'onda
  complex<double> *E,*B;    // vettori E e B  
                                      // (4 colonne per le polarizzazioni 
                                      //  e 3 righe per x,y,z)  
};

Onda::Onda() {
  this->Kx=NULL;
  this->Kz=NULL;
  this->E=NULL;
  this->B=NULL;
}

void Onda::mostra(int sp) {

  if(this->E) {
    for(int i=0; i<4; i++) {
      printf(" \nE vettore %d \n", i);
      for(int k=0; k<3; k++) {
	printf(" %e  %e     ", this->E[sp*3*4+i*3+k].real(),this->E[sp*3*4+i*3+k].imag());
      }
    }
  }
  if(this->B) {
    for(int i=0; i<4; i++) {
      printf(" \nB vettore %d \n", i);
      for(int k=0; k<3; k++) {
	printf(" %e  %e     ", this->B[sp*3*4+i*3+k].real(),this->B[sp*3*4+i*3+k].imag());
      }
    }
  }
}

Onda::~Onda() {
  // printf(" DISTRUGGO ONDA\n");
  complex<double> **ptr;
  complex<double> *ptrlist[] = { this->Kx,this->Kz,  this->E,  this->B  ,      NULL};
  ptr=&ptrlist[0];
  while(*ptr) {
      delete(*ptr);
    ptr++;
  }
}

extern  "C" {
  int zgeev_(char  * ,char  *,int *,complex<double> * , int *,complex<double> * ,complex<double> * ,int *,
	     complex<double> * , int *,  complex<double> * ,int *,double * ,int *);

} 

void transpose_3d(complex<double> *A) {
  complex<double> swap[9];

                 swap[1]=A[3];   swap[2]=A[6];
   swap[3]=A[1];                 swap[5]=A[7];
  swap[6]=A[2];   swap[7]=A[5];

                 A[1]=swap[1];   A[2]=swap[2];
  A[3]=swap[3];                  A[5]=swap[5];
  A[6]=swap[6];   A[7]=swap[7];

}

//**************************************************
//** DATA LA MATRICE V, QUESTA FUNZIONE TROVA
//** A PARTIRE DA QUALE PUNTO SI HA UN ZERO "NUMERICO"
//**   IL VALORE RITORNATO E LA POSIZIONE NELL ARRAY LINEARE
//**   SU CUI E ALLOCAT LA MATRICE

int staccoP( double *V, int n, double st)
{
  int i;
  for(i=1;i<n;i++)
    {
      // cout << " i  " << i << "  " << V[i] << "  " << V[i-1] << endl ;
      if( fabs(V[i]/V[i-1])<st) break;
    }
  if(i==n) i=0;
  return i;
}
//**
//***********************************************

void uSVD_3d( complex<double> *A, complex<double> * U,
	  double * S,  complex<double> * V ) 
{
   // temporaries required by lapack
   static char jobu= 'A',jobvt= 'A';
   static int info,lda= 3,ldu= 3,ldvt= 3,lwork= 32 ,m= 3,n= 3;
   static complex<double>  work[32];
   static double rwork[20];

   
   transpose_3d(A);
   
   zgesvd_( &jobu, &jobvt, &m, &n, A, &lda, S,
	    U, &ldu, V, &ldvt, work, &lwork, rwork, &info);
   if ( info != 0 )
     {
       if ( info > 0 )
	 {
	   cout << "void uSVD(double) failed to converge\n";
	   abort();
	 }
       else
	 {
	   cout << "void uSVD(double) error in parameter " << -1*info << endl;
	   abort();
	 }
     }
   
   transpose_3d(U);
      
   
}









void setS(complex<double> *S,complex<double> *die, double kx ) {

  static double  kx2,kx3,kx4,kx5,kx6;
  static complex<double> a0,a1,a2,a3,a4;
  kx2=kx*kx;
  kx3=kx2*kx;
  kx4=kx2*kx2;
  kx5=kx2*kx3;
  kx6=kx3*kx3;



    
#define d(i,j)  die[(i-1)*3+(j)-1]
 
#define Power(a,n)  a##n
    
    a4= d(3,3)  ; 
    
    a0= d(1,1)*Power(kx,4) + 
      d(1,2)*d(2,1)*Power(kx,2) - d(1,1)*d(2,2)*Power(kx,2) + 
      d(1,3)*d(3,1)*Power(kx,2) - d(1,1)*d(3,3)*Power(kx,2) - 
      d(1,3)*d(2,2)*d(3,1) + d(1,2)*d(2,3)*d(3,1)+ 
      d(1,3)*d(2,1)*d(3,2) - d(1,1)*d(2,3)*d(3,2) - 
      d(1,2)*d(2,1)*d(3,3) + d(1,1)*d(2,2)*d(3,3) ; 
    
    a3=d(1,3)*kx + d(3,1)*kx ; 
    
    a2=d(1,1)*Power(kx,2) + d(3,3)*Power(kx,2) + 
      d(1,3)*d(3,1) + d(2,3)*d(3,2) - d(1,1)*d(3,3) - 
      d(2,2)*d(3,3)   ;
    
    a1= d(1,3)*Power(kx,3) + d(3,1)*Power(kx,3) - 
      d(1,3)*d(2,2)*kx + d(1,2)*d(2,3)*kx - 
      d(2,2)*d(3,1)*kx + d(2,1)*d(3,2)*kx  ;
    
#undef Power
#undef d    
 

#define SS(i,j)  S[  (i)*4+(j) ]
    for(int i=0; i<4; i++) {
      for(int j=0; j<4; j++) {
	SS(i,j)=0.0;
      }
    }
    for(int i=0; i<3; i++)
      {
	SS(i+1,i) =  1;
      }    
    SS(0,0)= a3/a4;
    SS(0,1)=-a2/a4;
    SS(0,2)= a1/a4;
    SS(0,3)=-a0/a4;

//      cout << " Matrice SS  " << endl;
//      for(int i=0; i<4; i++) {
//        for(int j=0; j<4; j++) {
//  	cout << SS(i,j)<< " " ; 
//        }
//        cout << endl;
//      }
//      cout << endl;

#undef SS


}

void findeigenValues(complex<double> *S,complex<double> *eigvalues) {
  static char jobvl='N', jobvr='N';
  static int n=4;
  static int lda=4;
  static int info;
  static int lwork = 4*8;

  static complex<double> * work = new complex<double> [ lwork ];
  static double          *rwork = new double          [ 8 ];

  static int dum1=1, dum2=1;
  // cout << " chiamo zgeev \n";
  zgeev_(&jobvl,&jobvr,&n,S,&lda,eigvalues,NULL ,&dum1,NULL ,&dum2,work,&lwork,rwork,&info);
  // cout << "  zgeev  OK  \n";
  for ( int i=0; i<4; i++) {
    eigvalues[i] = -eigvalues[i] ;
  }
  riordina(eigvalues, 4 );

}
  
  
//****************************************************************
//** Trova le 4 diverse polarizzazioni propagantesi nel mezzo die
//** La onda che si passa in argomento deve essere inizializzata
//** precedentemente con l'angolo e l'energia desiderati
//**
void Onda::CreateModes( )
{
  if(this->indexes->nd==1) {
    // cout << " Scalar implementato !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! "<< endl;
    this->CreateModesScalar();
  } else if(this->indexes->nd==3) {
    this->CreateModesTensor();
  } else {
    cout << " Problem : neither nd==1 nor nd==3 in CreateModes\n";
    exit(0);
  }
}
PyObject * Onda::get_vectors() {
  int scanLen ;
  scanLen = this-> angles->dimensions[0];
  PyObject * res = PyList_New(0); 
  {
    int nd=1;
    int dims[]={scanLen};
    PyObject * nuovo;
    nuovo=PyArray_FromDims(nd,dims, 'D') ;
    memcpy( ((PyArrayObject*)nuovo)->data, this->  Kx,       scanLen*sizeof(complex<double>) );
    PyList_Append(  res ,  nuovo     ) ;
    Py_DECREF(nuovo);
  }
  {
    int nd=2;
    int dims[]={scanLen,4};
    PyObject * nuovo;
    nuovo=PyArray_FromDims(nd,dims, 'D' ) ;
    memcpy( ((PyArrayObject*)nuovo)->data, this->  Kz,       scanLen*4*sizeof(complex<double>) );
    PyList_Append(res ,     nuovo    ) ;
    Py_DECREF(nuovo);
  }
  {
    int nd=3;
    int dims[]={scanLen,4,3};
    PyObject * nuovo;
    nuovo=PyArray_FromDims(nd,dims, 'D' ) ;
    memcpy( ((PyArrayObject*)nuovo)->data, this->  E,       4*3*scanLen*sizeof(complex<double>) );
    PyList_Append(res ,     nuovo    ) ;
    Py_DECREF(nuovo);
  }
  {
    int nd=3;
    int dims[]={scanLen,4,3};
    PyObject * nuovo;
    nuovo=PyArray_FromDims(nd,dims, 'D' ) ;
    memcpy( ((PyArrayObject*)nuovo)->data, this->  B,       4*3*scanLen*sizeof(complex<double>) );
    PyList_Append(res ,     nuovo    ) ;
    Py_DECREF(nuovo);
  }
  return  res;
}

void Onda::CreateModesTensor( ) {
  int scanLen ;
  scanLen = this-> angles->dimensions[0];
  
  
  
  this-> Kx = new  complex<double> [ scanLen ] ;
  // creazione di Kz , E,  B
  this-> Kz= new complex<double> [ scanLen *4];
  this-> E = new complex<double> [ scanLen *4 *3 ];
  this-> B = new complex<double> [ scanLen *4 *3 ];


  
  for(int scanPoint =0 ; scanPoint < scanLen; scanPoint++) {
    this->Kx[scanPoint]= cos(((double *) this->angles->data)[scanPoint]);
  }
  
  
  
  static complex<double > S[16] ;

  
  
  static double kx;

  complex<double>  *die;

  for(int scanPoint =0 ; scanPoint < scanLen; scanPoint++) {
    


   die = ( (complex<double> *) this->indexes->data) + scanPoint*9 ;



    kx = this->Kx[scanPoint].real();

    setS(S,die, kx );

    // LEAK no perche e statica

    static complex<double> *eigvalues =new complex<double> [4];
    findeigenValues(S,eigvalues);
    // cout << " AUTOVALORI =====\n";for ( int i=0; i<4; i++) {cout << eigvalues[i]<< endl;}


#define KZ(i) this->Kz[scanPoint*4+i]
    for ( int i=0; i<4; i++) {
      KZ(i)=eigvalues[i];
    }
#undef KZ

    int stacco0=0;
    for(int iwaves=0; iwaves<4; )
      {
	static  complex<double> Sistem[3][3];

        // cout << " Sistema " << endl;

	complex<double> K[3];

	K[0] = kx ; K[1]=0; K[2]= eigvalues[iwaves] ;
	for(int ix=0; ix<3; ix++)
	  {
	    for(int iy=0; iy<3; iy++)
	      {
		Sistem[ix][iy] = K[ix]*K[iy] ;

#define d(i,j)  die[(i)*3+(j)]
                Sistem[ix][iy]+=d(ix,iy);
#undef d


	      }
	    Sistem[ix][ix] -= kx*kx +eigvalues[iwaves]* eigvalues[iwaves] ;
	    // for(int iy=0; iy<3; iy++) cout << Sistem[ix][iy] << " ";
	    // cout << endl;

	  }
	{
	  static complex<double> U[3][3],W[3][3];
	  static double V[3];
	  uSVD_3d ( (complex<double>*)Sistem ,(complex<double>*) U , (double*)V  ,(complex<double>*) W );
	  {
	    int staccoat=staccoP(V, 3, 1.e-7), dim;
	    if(!staccoat)
	      {
		// cout << " problema con stacco\n";
		staccoat=2;  
		// (0);
	      } 
	    if(iwaves==0)
	      {
		stacco0=staccoat;
	      }
	    if(staccoat!=stacco0)
	      {
		cerr << " anomalia, stacco!=stacco0 \n";
		staccoat=stacco0;
	      }
	    
	    dim =3-staccoat;
	    // cout << " ==========stacco " << staccoat << endl;
            for(int k=0; k<3; k++) {
 
	      for(int z=staccoat; z< 3; z++) {
#define E(i,j,k)      E[  i*3*4     +(j)*3        +k     ]
		E( scanPoint , (iwaves +z-staccoat) , k ) = conj(W[k][z]);
		// cout << W[k][z] << "  " ;
#undef E
	      }
	      // cout << endl;
	    }
	    
	    //  	    onda.E.insert(uIndex(0,1,2), uIndex(i,1,i+dim-1),
	    //  			  W(uIndex(0,1,2),uIndex(staccoat,1,2))
	    //  			  );
	    iwaves += dim;
	  } 
	  
	}
      }
    for(int colo=0; colo<4; colo++)
      {
	static complex<double> E2B[3][3], K[3];
        K[1]=0;
        K[0]=kx;
        K[2]=eigvalues[colo];
        // cout <<"K[0] " << K[0] << "   K[2]  "<< K[2] << endl;
	for(int i=0; i<3; i++) for(int k=0; k<3; k++)  E2B[i][k]=0;
        E2B[0][1]=-K[2];        E2B[1][0]=+K[2];
        E2B[1][2]=-K[0];        E2B[2][1]=+K[0];
        E2B[2][0]=-K[1];        E2B[0][2]=+K[1];
	// cout << " E2B " << endl;
       {
	 for( int i=0; i<0*3; i++) {
	   for(int j=0; j<3; j++) {
	     cout << E2B[i][j] << " " ;
	   }
	   cout << endl;
	 }
       }
#define E(j,k)      E[  scanPoint*3*4     +(j)*3        +k     ]
#define B(j,k)      B[  scanPoint*3*4     +(j)*3        +k     ]
	
	for(int j=0; j<3; j++) {
	  B(colo,j)=0;
	  
	  for(int k=0; k<3; k++) {
	    B (colo,j) +=   E2B[j][k] * E(colo,k  ) ;
	  }
	}
#undef E
#undef B	
      }
    // this->mostra(scanPoint);
  }




}
  


void Onda::CreateModesTensor_Vide( ) {
  int scanLen ;
  scanLen = this-> angles->dimensions[0];
  
  this-> Kx = new  complex<double> [ scanLen ] ;
  // creazione di Kz , E,  B
  this-> Kz= new complex<double> [ scanLen *4];
  this-> E = new complex<double> [ scanLen *4 *3 ];
  this-> B = new complex<double> [ scanLen *4 *3 ];

  for(int scanPoint =0 ; scanPoint < scanLen; scanPoint++) {
    this->Kx[scanPoint]= cos(((double *) this->angles->data)[scanPoint]);
  }
  
  // static complex<double > S[16] ;

  static double kx;


  for(int scanPoint =0 ; scanPoint < scanLen; scanPoint++) {

    kx = this->Kx[scanPoint].real();

#define KZ(i) this->Kz[scanPoint*4+i]
    for ( int i=0; i<2; i++) {
      KZ(i)=+sqrt(1.0-kx*kx);
    }
    for ( int i=2; i<4; i++) {
      KZ(i)=-sqrt(1.0-kx*kx);
    }



#define E(i,j,k)      E[  i*3*4     +(j)*3        +k     ]
#define B(i,j,k)      B[  i*3*4     +(j)*3        +k     ]

    E( scanPoint , 0 , 0 ) = 0.0;
    E( scanPoint , 2 , 0 ) = 0.0;

    E( scanPoint , 1 , 0 ) = -KZ(1);
    E( scanPoint , 3 , 0 ) = -KZ(3);


    E( scanPoint , 0 , 1 ) = 1.0;
    E( scanPoint , 2 , 1 ) = 1.0;
    E( scanPoint , 1 , 1 ) = 0.0;
    E( scanPoint , 3 , 1 ) = 0.0;


    E( scanPoint , 0 , 2 ) = 0.0;
    E( scanPoint , 2 , 2 ) = 0.0;
    E( scanPoint , 1 , 2 ) = kx;
    E( scanPoint , 3 , 2 ) = kx;

    ///////////////////////////////////////////////////////////
    B( scanPoint , 0 , 0 ) = -KZ(0);
    B( scanPoint , 2 , 0 ) = -KZ(2);

    B( scanPoint , 1 , 0 ) = 0.0;
    B( scanPoint , 3 , 0 ) = 0.0;


    B( scanPoint , 0 , 1 ) = 0.0;
    B( scanPoint , 2 , 1 ) = 0.0;
    B( scanPoint , 1 , 1 ) = -1.0;
    B( scanPoint , 3 , 1 ) = -1.0;


    B( scanPoint , 0 , 2 ) = kx;
    B( scanPoint , 2 , 2 ) = kx;
    B( scanPoint , 1 , 2 ) = 0.0;
    B( scanPoint , 3 , 2 ) = 0.0;

#undef E

#undef KZ



  }
}






void Onda::CreateModesScalar( ) {
  int scanLen ;
  scanLen = this-> angles->dimensions[0];
  
  this-> Kx = new  complex<double> [ scanLen ] ;
  // creazione di Kz , E,  B
  this-> Kz= new complex<double> [ scanLen *4];
  this-> E = new complex<double> [ scanLen *4 *3 ];
  this-> B = new complex<double> [ scanLen *4 *3 ];

  for(int scanPoint =0 ; scanPoint < scanLen; scanPoint++) {
    this->Kx[scanPoint]= cos(((double *) this->angles->data)[scanPoint]);
  }
  
  // static complex<double > S[16] ;

  static double kx;


  for(int scanPoint =0 ; scanPoint < scanLen; scanPoint++) {

    kx = this->Kx[scanPoint].real();

    complex<double> nind,eps;
    nind =  *(((complex<double>*) (this->indexes->data))+scanPoint);

    // eps=nind*nind;
    eps=nind;

#define KZ(i) this->Kz[scanPoint*4+i]
    for ( int i=0; i<2; i++) {
      KZ(i)=+sqrt(1.0*eps-kx*kx);
    }
    for ( int i=2; i<4; i++) {
      KZ(i)=-sqrt(1.0*eps-kx*kx);
    }



#define E(i,j,k)      E[  i*3*4     +(j)*3        +k     ]
#define B(i,j,k)      B[  i*3*4     +(j)*3        +k     ]

    E( scanPoint , 0 , 0 ) = 0.0;
    E( scanPoint , 2 , 0 ) = 0.0;

    E( scanPoint , 1 , 0 ) = -KZ(1);
    E( scanPoint , 3 , 0 ) = -KZ(3);


    E( scanPoint , 0 , 1 ) = 1.0;
    E( scanPoint , 2 , 1 ) = 1.0;
    E( scanPoint , 1 , 1 ) = 0.0;
    E( scanPoint , 3 , 1 ) = 0.0;


    E( scanPoint , 0 , 2 ) = 0.0;
    E( scanPoint , 2 , 2 ) = 0.0;
    E( scanPoint , 1 , 2 ) = kx;
    E( scanPoint , 3 , 2 ) = kx;

    ///////////////////////////////////////////////////////////
    B( scanPoint , 0 , 0 ) = -KZ(0);
    B( scanPoint , 2 , 0 ) = -KZ(2);

    B( scanPoint , 1 , 0 ) = 0.0;
    B( scanPoint , 3 , 0 ) = 0.0;


    B( scanPoint , 0 , 1 ) = 0.0;
    B( scanPoint , 2 , 1 ) = 0.0;
    B( scanPoint , 1 , 1 ) = -1.0*eps;
    B( scanPoint , 3 , 1 ) = -1.0*eps;


    B( scanPoint , 0 , 2 ) = kx;
    B( scanPoint , 2 , 2 ) = kx;
    B( scanPoint , 1 , 2 ) = 0.0;
    B( scanPoint , 3 , 2 ) = 0.0;

#undef E

#undef KZ



  }
}
  














  

int ThoseArraysAreEqual( PyArrayObject *a  , PyArrayObject *b   ) {
  int size_a=1;
  if( a->nd!=b->nd) return 0;
  for (int k=0; k<a->nd; k++){
    size_a*=a->dimensions[k];
    if(a->dimensions[k]!=b->dimensions[k]) return 0;
  };
  if(  (a->descr->type  != b->descr->type) || a->descr->type !='D' ) {
    cout << " ThoseArraysAreEqual is made for Complex \n";
    exit(0);
  }
  int res= ( memcmp( a->data,  b->data, size_a*sizeof(complex<double>) )   ==0);
  return res ;
}

static char multicouche_calculatescanTensor_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_calculatescanTensor(PyObject *self, PyObject *args)
{
  // static  complex<double> Sistem[9];



  PyArrayObject *thickness, *roughness, *wavelenght, *angle, /* *indexes , */*Sindex;
  PyObject * indexes;  // indexes are now a python list of pyarrayobjects
  PyFloatObject *Srough;
  

  
  if(!PyArg_ParseTuple(args,"OOOOOOO:multicouche_calculatescanTensor", (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(!PyList_Check((PyObject *)indexes))        onError("not a PyList, 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 " ) ;


  int num_thicks = PyList_Size(indexes) ;
  for(int k=0; k< num_thicks; k++) {
    if(!PyArray_Check((PyObject *)  PyList_GetItem(indexes,k )   )   ) {
      char errormessage[200];
      sprintf(errormessage,"not a PyArray, argument 5, item number %d",k);
      onError(errormessage);
    }
    if( ((PyArrayObject*) PyList_GetItem(indexes,k )) ->descr->type_num != PyArray_CDOUBLE )  {
      char errormessage[200];
      sprintf(errormessage,"not a CDOUBLE PyArray, argument 5, item number %d",k);
      onError(errormessage);
    }
  }

  if( Sindex->descr->type_num != PyArray_CDOUBLE ) onError(" Sindex is not double " ) ;
  
  /* check the Number of 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");
  */

  int nd_layer_i[num_thicks];
  for(int k=0; k< num_thicks; k++) {
    int nd =  ((PyArrayObject *)  PyList_GetItem(indexes,k )   ) ->nd;
    if(  nd!=1 && nd!=3   ) {
      char errormessage[200];
      sprintf(errormessage,"PyArray, argument 5, item number %d is neither an array neither a tensor array, nd= %d",k,nd);
      onError(errormessage);
    }
    nd_layer_i[k]=nd;
  }
  int nd_subs_i;
  {
    int nd =  ((PyArrayObject *)  Sindex  ) ->nd;
    if(  nd != 1 && nd != 3  )
      onError("The Sindex array (arg. 6) is neither an array neither a tensor array");
    nd_subs_i=nd;
  }

  /* check the 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!=PyList_Size(indexes)   )
    onError("The indexes  array (arg. 5) has not the right  lenght");
  
  for(int k=0; k< numthick; k++) {
    PyArrayObject *li;
    li = ((PyArrayObject *)  PyList_GetItem(indexes,k ))  ;
    if(   li->dimensions[0] != numwave   ) {
      onError("The indexes  array (arg. 5) has not the right second   dimension");
    }
    if( nd_layer_i[k]==3) {
      if(   li->dimensions[1] != 3 || li->dimensions[2] != 3   ) {
	onError("The indexes  array (arg. 5) has not the right third and/or  fourth where it should be a tensor");
      }
    }
  }
  
  if(numwave!=Sindex->dimensions[0] )
    onError("The Sindex  array (arg. 6) has not the right  second dimension");
  {
    if(nd_subs_i==3) {
      if(   Sindex->dimensions[1] != 3 || Sindex->dimensions[2] != 3   ) {
	onError("The Substrate Index  array (arg. 6) has not the right dimensions where it should be a tensor");
      }      
    }
  }
  
  /* check that everything is contiguous */
  
  PyArrayObject *toverify[]={thickness,roughness,wavelenght,angle,Sindex,NULL};
  PyArrayObject **ptr;
  ptr=toverify;
  while(*ptr) {
    if((*ptr)->flags %2 == 0) onError(" All arrays have to be contiguous, either solve this problem at the python level or at the c level using PyArray_CopyFromObject");
    ptr++;
  }
  for(int k=0; k< numthick; k++) {
    PyArrayObject *li;
    li = ((PyArrayObject *)  PyList_GetItem(indexes,k ))  ;
    if((li)->flags %2 == 0) onError(" All arrays in Indexes have to be contiguous, either solve this problem at the python level or at the c level using PyArray_CopyFromObject");
  }

  /* take into account how many eigenmodes have to be calculated */
  Onda * Eigenmodes= new Onda [numthick];
  Onda EigenSubstrate;
  Onda * * wheretolook = new   Onda * [numthick];
  int ntodo=0;


  for(int k=0; k< numthick; k++) {
    PyArrayObject *li;
    li = ((PyArrayObject *)  PyList_GetItem(indexes,k ))  ;
    int isnew=1;
    for(int l=0; l< ntodo; l++) {
      if(   ThoseArraysAreEqual(li  ,  Eigenmodes[l].indexes)    ) {
	isnew=0;
	wheretolook[k]= &( Eigenmodes[l]  );
	break;
      }
    }
    if(isnew){
      Eigenmodes[ntodo].indexes = li;
      Eigenmodes[ntodo].angles=angle    ;
      wheretolook[k]= &( Eigenmodes[ntodo]  );
      ntodo++;
    }
  }

  EigenSubstrate.indexes = Sindex;
  EigenSubstrate.angles=angle;

  PyObject * res = PyList_New(0); 

  EigenSubstrate.CreateModes();
  PyObject * nuovo = EigenSubstrate.get_vectors();
  PyList_Append(res ,   nuovo     ) ;
  Py_DECREF(nuovo);

  for(int i=0; i<ntodo; i++) {
    if(i==ntodo-1) {
      Eigenmodes[i].CreateModesTensor_Vide();
    } else {
      Eigenmodes[i].CreateModes();
    }
//      PyList_Append(res ,  Eigenmodes[i] . get_vectors()      ) ;
  }

  for(int k=0; k< numthick; k++) {
    PyObject * nuovo = wheretolook[k]-> get_vectors();
    PyList_Append(res ,     nuovo   ) ;
    Py_DECREF(nuovo);

  }

  /* creation of the return Values */

  delete []  Eigenmodes;
  delete [] wheretolook ;

  return res;
}


static char MAPDOT_doc[] = "\n"
"/*  2 arguments are taken by the function :\n"
" *  matA(NW,2,2), matB(NW,2,2)\n"
" *    where NW is the number of wavelenghts\n"
" */\n";
static PyObject *
MAPDOT(PyObject *self, PyObject *args)
{
  PyArrayObject *A, *B;
  int nw, nnw;
  
  if(!PyArg_ParseTuple(args,"OO:MAPDOT", (PyObject *) &A,(PyObject *)&B )
     )
    return NULL;

   PyArrayObject * tocheck[]= {A,B};


  nw=-10;
  for(int check=0; check<2; check++) {
  
    if(!PyArray_Check((PyObject *)tocheck[check] ) )    onError("argument is not a PyArray, routine MAPDOT");
    if( tocheck[check]->descr->type_num != PyArray_CDOUBLE ) onError(" argument  is not Cdouble, MAPDOT " ) ;
    if( tocheck[check]->nd != 3 ) onError("argument  has not three dimensions, MAPDOT");
    nnw = tocheck[check]->dimensions[0] ; 
    if(nw>0 && nnw!=nw) {
      onError("arguments  have not the same number of wavelenghts, MAPDOT");
    }
    nw=nnw;
    if(tocheck[check]->dimensions[1]!=2  || tocheck[check]->dimensions[2]!=2  ) {
      printf(" %d %d %d\n", tocheck[check]->dimensions[0], tocheck[check]->dimensions[1],tocheck[check]->dimensions[2]);
      onError("argument  is not NWX2X2 dimensioned , MAPDOT");
    }
  }
  
  PyArrayObject *C;
  int d[3];
  d[0]=nw ;
  d[1]=2 ;
  d[2]=2 ;
  C= (PyArrayObject*)PyArray_FromDims( 3, d,  PyArray_CDOUBLE );

  complex<double> *a,*b,*c;

  a= (complex<double> *)A->data;
  b= (complex<double> *)B->data;
  c= (complex<double> *)C->data;
  
  int ast0,ast1,ast2;
  int bst0,bst1,bst2;
  int cst0,cst1,cst2;


  ast0 = A->strides[0]/sizeof(complex<double>);
  ast1 = A->strides[1]/sizeof(complex<double>);
  ast2 = A->strides[2]/sizeof(complex<double>);

  bst0 = B->strides[0]/sizeof(complex<double>);
  bst1 = B->strides[1]/sizeof(complex<double>);
  bst2 = B->strides[2]/sizeof(complex<double>);

  cst0 = C->strides[0]/sizeof(complex<double>);
  cst1 = C->strides[1]/sizeof(complex<double>);
  cst2 = C->strides[2]/sizeof(complex<double>);


#define posA(iw,i,j) a[( (iw)*ast0 + i*ast1 + j*ast2  )]
#define posB(iw,i,j) b[( (iw)*bst0 + i*bst1 + j*bst2  )]
#define posC(iw,i,j) c[( (iw)*cst0 + i*cst1 + j*cst2  )]

  for(int iw=0; iw<nw; iw++) {
    posC(iw,0,0)  =  posA(iw,0,0)*posB(iw,0,0) + posA(iw,0,1)*posB(iw,1,0);
    posC(iw,0,1)  =  posA(iw,0,0)*posB(iw,0,1) + posA(iw,0,1)*posB(iw,1,1);
    posC(iw,1,0)  =  posA(iw,1,0)*posB(iw,0,0) + posA(iw,1,1)*posB(iw,1,0);
    posC(iw,1,1)  =  posA(iw,1,0)*posB(iw,0,1) + posA(iw,1,1)*posB(iw,1,1);
  }
#undef posA
#undef posB
#undef posC

  return (PyObject * )C ; 
}


static char MAPDOT4_doc[] = "\n"
"/*  2 arguments are taken by the function :\n"
" *  matA(NW,2,2), matB(NW,2,2)\n"
" *    where NW is the number of wavelenghts\n"
" */\n";
static PyObject *
MAPDOT4(PyObject *self, PyObject *args)
{
  PyArrayObject *A, *B;
  int nw, nnw;
  
  if(!PyArg_ParseTuple(args,"OO:MAPDOT4", (PyObject *) &A,(PyObject *)&B )
     )
    return NULL;

   PyArrayObject * tocheck[]= {A,B};
   PyArrayObject ** tochange[]= {&A,&B};
   int todestroy[] = { 0,0};


  nw=-10;
  for(int check=0; check<2; check++) {
  
    if(!PyArray_Check((PyObject *)tocheck[check] ) )    onError("argument is not a PyArray, routine MAPDOT4");
    if((tocheck[check])->flags %2 == 0) {
      todestroy[check] = 1;
      *tochange[check]   = (PyArrayObject *) PyArray_CopyFromObject((PyObject*)tocheck[check], PyArray_CDOUBLE, 3,3); ;
      tocheck[check] = *tochange[check];
      // onError(" All arrays have to be contiguous, routine MAPDOT4");
    }
    if( tocheck[check]->descr->type_num != PyArray_CDOUBLE ) onError(" argument  is not Cdouble, MAPDOT4 " ) ;
    if( tocheck[check]->nd != 3 ) onError("argument  has not three dimensions, MAPDOT4");
    nnw = tocheck[check]->dimensions[0] ; 
    if(nw>0 && nnw!=nw) {
      onError("arguments  have not the same number of wavelenghts, MAPDOT4");
    }
    nw=nnw;
    if(tocheck[check]->dimensions[1]!=4  || tocheck[check]->dimensions[2]!=4  ) {
      printf(" %d %d %d\n", tocheck[check]->dimensions[0], tocheck[check]->dimensions[1],tocheck[check]->dimensions[2]);
      onError("argument  is not NWX4X4 dimensioned , MAPDOT4");
    }
  }
  
  PyArrayObject *C;
  int d[3];
  d[0]=nw ;
  d[1]=4 ;
  d[2]=4 ;
  C= (PyArrayObject*)PyArray_FromDims( 3, d,  PyArray_CDOUBLE );

  complex<double> *a,*b,*c;

  a= (complex<double> *)A->data;
  b= (complex<double> *)B->data;
  c= (complex<double> *)C->data;
  
  int ast0,ast1,ast2;
  int bst0,bst1,bst2;
  int cst0,cst1,cst2;


  ast0 = A->strides[0]/sizeof(complex<double>);
  ast1 = A->strides[1]/sizeof(complex<double>);
  ast2 = A->strides[2]/sizeof(complex<double>);

  bst0 = B->strides[0]/sizeof(complex<double>);
  bst1 = B->strides[1]/sizeof(complex<double>);
  bst2 = B->strides[2]/sizeof(complex<double>);

  cst0 = C->strides[0]/sizeof(complex<double>);
  cst1 = C->strides[1]/sizeof(complex<double>);
  cst2 = C->strides[2]/sizeof(complex<double>);


#define posA(iw,i,j) a[( (iw)*ast0 + i*ast1 + j*ast2  )]
#define posB(iw,i,j) b[( (iw)*bst0 + i*bst1 + j*bst2  )]
#define posC(iw,i,j) c[( (iw)*cst0 + i*cst1 + j*cst2  )]
  char *transa = "N";
  char *transb = "N";
  char *transc = "N";
  int m = 4;
  int n = 4;
  int k = 4;
  int lda=4, ldb=4, ldc=4;
  Complex alpha(1.0);
  Complex beta(0.0);
  for(int iw=0; iw<nw; iw++) {

     zgemm_( transb, transa, &m, &n, &k, 
	     &alpha, &posB(iw,0,0), &lda, &posA(iw,0,0), &ldb, 
	     &beta, &posC(iw,0,0), &ldc );
  }
#undef posA
#undef posB
#undef posC


  for(int check=0; check<2; check++) {
    if(todestroy[check]) {
      Py_DECREF(tocheck[check]);
    }
  }
  return (PyObject * )C ; 
}




static char MAPDOT_V_doc[] = "\n"
"/*  2 arguments are taken by the function :\n"
" *  matA(NW,2,2), matB(NW,2)\n"
" *    where NW is the number of wavelenghts\n"
" */\n";
static PyObject *
MAPDOT_V(PyObject *self, PyObject *args)
{
  PyArrayObject *A, *B;
  int nw, nnw;
  
  if(!PyArg_ParseTuple(args,"OO:MAPDOT_V", (PyObject *) &A,(PyObject *)&B )
     )
    return NULL;

   PyArrayObject * tocheck[]= {A,B};


  nw=-10;
  for(int check=0; check<2; check++) {
  
    if(!PyArray_Check((PyObject *)tocheck[check] ) )    onError("argument is not a PyArray, routine MAPDOT");
    if( tocheck[check]->descr->type_num != PyArray_CDOUBLE ) onError(" argument  is not Cdouble, MAPDOT " ) ;

    if(check==0) {
      if( tocheck[check]->nd != 3 ) onError("argument  has not three dimensions, MAPDOT");
    } else {
      if( tocheck[check]->nd != 2 ) onError("argument  has not three dimensions, MAPDOT");
    }
    nnw = tocheck[check]->dimensions[0] ; 
    
    if(nw>0 && nnw!=nw) {
      onError("arguments  have not the same number of wavelenghts, MAPDOT");
    }
    
    nw=nnw;
    
    if(check==0) {
      if(tocheck[check]->dimensions[1]!=2  || tocheck[check]->dimensions[2]!=2  ) {
	onError("argument  is not NWX2X2 dimensioned , MAPDOT");
      }
    } else {
      if(tocheck[check]->dimensions[1]!=2   ) {
	onError("argument  is not NWX2 dimensioned , MAPDOT");
      }
    }
  }
  
  PyArrayObject *C;
  int d[2];
  d[0]=nw ;
  d[1]=2 ;

  C= (PyArrayObject*)PyArray_FromDims( 2, d,  PyArray_CDOUBLE );

  complex<double> *a,*b,*c;

  a= (complex<double> *)A->data;
  b= (complex<double> *)B->data;
  c= (complex<double> *)C->data;
  
  int ast0,ast1,ast2;
  int bst0,bst1;
  int cst0,cst1;


  ast0 = A->strides[0]/sizeof(complex<double>);
  ast1 = A->strides[1]/sizeof(complex<double>);
  ast2 = A->strides[2]/sizeof(complex<double>);

  bst0 = B->strides[0]/sizeof(complex<double>);
  bst1 = B->strides[1]/sizeof(complex<double>);

  cst0 = C->strides[0]/sizeof(complex<double>);
  cst1 = C->strides[1]/sizeof(complex<double>);


#define posA(iw,i,j) a[( (iw)*ast0 + i*ast1 + j*ast2  )]
#define posB(iw,i) b[( (iw)*bst0 + i*bst1 )]
#define posC(iw,i) c[( (iw)*cst0 + i*cst1 )]

  for(int iw=0; iw<nw; iw++) {
    posC(iw,0)  =  posA(iw,0,0)*posB(iw,0) + posA(iw,0,1)*posB(iw,1);
    posC(iw,1)  =  posA(iw,1,0)*posB(iw,0) + posA(iw,1,1)*posB(iw,1);
  }
#undef posA
#undef posB
#undef posC

  return (PyObject * )C ; 
}















static char MAPINV_doc[] = "\n"
"/*  1 arguments is taken by the function :\n"
" *  matA(NW,2,2)\n"
" *    where NW is the number of wavelenghts\n"
" */\n";

static PyObject *
MAPINV(PyObject *self, PyObject *args)
{
  PyArrayObject *A;
  int nw;
  
  if(!PyArg_ParseTuple(args,"O:MAPINV", (PyObject *) &A  )
     )
    return NULL;
  
  PyArrayObject * tocheck[]= {A};
  
  nw=-10;
  for(int check=0; check<1; check++) {
    
    if(!PyArray_Check((PyObject *)tocheck[check] ) )    onError("argument is not a PyArray, routine MAPINV");
    if( tocheck[check]->descr->type_num != PyArray_CDOUBLE ) onError(" argument  is not Cdouble, MAPINV " ) ;
    if( tocheck[check]->nd != 3 ) onError("argument  has not three dimensions, MAPINV");
    nw = tocheck[check]->dimensions[0] ; 
    if(tocheck[check]->dimensions[1]!=2  || tocheck[check]->dimensions[2]!=2  ) {
      printf(" %d %d %d\n", tocheck[check]->dimensions[0], tocheck[check]->dimensions[1],tocheck[check]->dimensions[2]);
      onError("argument  is not NWX2X2 dimensioned , MAPINV");
    }
  }
  
  PyArrayObject *C;
  int d[3];
  d[0]=nw ;
  d[1]=2 ;
  d[2]=2 ;
  C= (PyArrayObject*)PyArray_FromDimsAndData( 3, d,  PyArray_CDOUBLE,NULL);

  complex<double> *a,*c;

  a= (complex<double> *) A->data;
  c= (complex<double> *) C->data;
  
  int ast0,ast1,ast2;
  int cst0,cst1,cst2;

  ast0 = A->strides[0]/sizeof(complex<double>);
  ast1 = A->strides[1]/sizeof(complex<double>);
  ast2 = A->strides[2]/sizeof(complex<double>);

  cst0 = C->strides[0]/sizeof(complex<double>);
  cst1 = C->strides[1]/sizeof(complex<double>);
  cst2 = C->strides[2]/sizeof(complex<double>);

#define posA(iw,i,j) a[( (iw)*ast0 + i*ast1 + j*ast2  )]
#define posC(iw,i,j) c[( (iw)*cst0 + i*cst1 + j*cst2  )]
  complex<double> det;
  for(int iw=0; iw<nw; iw++) {
    det = posA(iw,0,0)*posA(iw,1,1)-posA(iw,0,1)*posA(iw,1,0);
    posC(iw,0,0)  =  posA(iw,1,1)/det  ;
    posC(iw,0,1)  = -posA(iw,0,1)/det;
    posC(iw,1,0)  = -posA(iw,1,0)/det;
    posC(iw,1,1)  =  posA(iw,0,0)/det;
  }
#undef posA
#undef posC
  return (PyObject *)C ; 
}




static char MAPINV4_doc[] = "\n"
"/*  1 arguments is taken by the function :\n"
" *  matA(NW,4,4)\n"
" *    where NW is the number of wavelenghts\n"
" */\n";


static PyObject *
MAPINV4(PyObject *self, PyObject *args)
{
  PyArrayObject *input, *C;
  int nw;
  
  if(!PyArg_ParseTuple(args,"O:MAPINV4", (PyObject *) &input  )
     )
    return NULL;
  
  PyArrayObject * tocheck[]= {input};
  
  nw=-10;
  for(int check=0; check<1; check++) {
    
    if(!PyArray_Check((PyObject *)tocheck[check] ) )    onError("argument is not a PyArray, routine MAPINV4");
    nw = tocheck[check]->dimensions[0] ; 
    if(tocheck[check]->dimensions[1]!=4  || tocheck[check]->dimensions[2]!=4  ) {
      onError("argument  is not NWX4X4 dimensioned , MAPINV");
    }
  }
  
  C = (PyArrayObject *)PyArray_CopyFromObject((PyObject *)input, PyArray_CDOUBLE, 3,3);
    
  complex<double> *c;
  
  c= (complex<double> *) C->data;
  
  int cst0,cst1,cst2;
  
  cst0 = C->strides[0]/sizeof(complex<double>);
  cst1 = C->strides[1]/sizeof(complex<double>);
  cst2 = C->strides[2]/sizeof(complex<double>);
  
#define posC(iw,i,j) c[( (iw)*cst0 + i*cst1 + j*cst2  )]
  
  int size=4;
  int ipiv[size];
  complex<double> * work = new complex<double> [size];
  int info;
  
  for(int iw=0; iw<nw; iw++) {
    zgetrf_( &size, &size, &(posC(iw,0,0)), &size, ipiv, &info );
    if ( info != 0 )      onError("Error while factorizing matrix in MAPINV4 )");
    zgetri_( &size, &(posC(iw,0,0))  , &size, ipiv, work, &size, &info );
    if ( info > 0 ) onError("Warning: Matrix is singular in MAPINV4");
    if ( info < 0 ) onError(" Error while inverting matrix  in MAPINV4");
  }
#undef posA
#undef posC
  delete [] work;
  return (PyObject *)C ; 
}

























static PyMethodDef multicouche_functions[] = {
  {"PPM_calculatescanTensor", multicouche_calculatescanTensor, METH_VARARGS , multicouche_calculatescanTensor_doc },
  {"PPM_calculatescan", multicouche_calculatescan, METH_VARARGS , multicouche_calculatescan_doc },
  {"PPM_MAPDOT", MAPDOT, METH_VARARGS , MAPDOT_doc },
  {"PPM_MAPDOT4", MAPDOT4, METH_VARARGS , MAPDOT4_doc },
  {"PPM_MAPDOT_V", MAPDOT_V, METH_VARARGS , MAPDOT_V_doc },
  {"PPM_MAPINV", MAPINV, METH_VARARGS , MAPINV_doc },
  {"PPM_MAPINV4", MAPINV4, METH_VARARGS , MAPINV4_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
}





















