//  Adaptateurs pour TMatrix TVector du package Sophya
//                             R. Ansari  1/99
// LAL (Orsay) / IN2P3-CNRS  DAPNIA/SPP (Saclay) / CEA

#include "sopnamsp.h"
#include "pitvmaad.h"
#include <math.h>

// -------------------------------------------------------
// Decodage de la presence dans la chaine opt des options d'affichage pour tableaux complexes 
//   cdreal , cdimag , cdphase , cdmod , cdmod2
PICmplxDispOption StringToCmplxDispOption( string & opt )
{
PICmplxDispOption dopt=PICDO_Module;
size_t olen = opt.length();
if (opt.find("cdmod2")<olen) dopt=PICDO_Module2;
else if (opt.find("cdmod")<olen) dopt=PICDO_Module;
else if (opt.find("cdphas")<olen) dopt=PICDO_Phase;
else if (opt.find("cdimag")<olen) dopt=PICDO_Imag;
else if (opt.find("cdreal")<olen) dopt=PICDO_Real;
else dopt=PICDO_Module;
return dopt;  
}

/* --- Methodes de calcul inline pour valeurs de retour complexes ---- */
static inline double  _z_mod_( complex< r_4 > z) 
{
  return ( sqrt((double)(z.real()*z.real()+z.imag()*z.imag() ) ) );
}
static inline double  _z_mod_( complex< r_8 > z) 
{
  return ( sqrt(z.real()*z.real()+z.imag()*z.imag() )  );
}
static inline double  _z_mod2_( complex< r_4 > z) 
{
  return ( (double)(z.real()*z.real()+z.imag()*z.imag() ) );
}
static inline double  _z_mod2_( complex< r_8 > z) 
{
  return ( z.real()*z.real()+z.imag()*z.imag() );
}
/* ----------------------------  */

/* --Methode-- */
template <class T> 
POTVectorAdapter<T>::POTVectorAdapter(TVector<T>* v, bool ad, PICmplxDispOption dopt)
        : P1DArrayAdapter(v->NElts())
{
dOpt = dopt;
aDel = ad; 
mVec = v; 
}
/* --Methode-- */
template <class T> 
POTVectorAdapter<T>::~POTVectorAdapter()
{
if (aDel) delete mVec;
}

/* --Methode-- */
template <class T> 
double POTVectorAdapter<T>::Value(int i)
{
  return((double)(*mVec)(i));
}

/* --Methode-- */
DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
double POTVectorAdapter< complex<float> >::Value(int i)
{
switch (dOpt) {
 case PICDO_Module:
   return _z_mod_( (*mVec)(i) );
   break;
 case PICDO_Real:
   return (double)(((*mVec)(i)).real());
   break;
 case PICDO_Imag:
   return (double)(((*mVec)(i)).imag());
   break;
 case PICDO_Phase:
   return atan2( (double) ((*mVec)(i).real()), (double)((*mVec)(i).imag()) );
   break;
 case PICDO_Module2:
   return _z_mod2_( (*mVec)(i) );
   break;
 }
}

/* --Methode-- */
DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
double POTVectorAdapter< complex<double> >::Value(int i)
{
switch (dOpt) {
 case PICDO_Module:
   return _z_mod_( (*mVec)(i) );
   break;
 case PICDO_Real:
   return (double)(((*mVec)(i)).real());
   break;
 case PICDO_Imag:
   return (double)(((*mVec)(i)).imag());
   break;
 case PICDO_Phase:
   return atan2( (double) ((*mVec)(i).real()), (double)((*mVec)(i).imag()) );
   break;
 case PICDO_Module2:
   return _z_mod2_( (*mVec)(i) );
   break;
 }
}

/* --Methode-- */
template <class T> 
POTMatrixAdapter<T>::POTMatrixAdapter(TMatrix<T>* mtx, bool ad, PICmplxDispOption dopt)
        : P2DArrayAdapter(mtx->NCols(), mtx->NRows())
{
dOpt = dopt;
aDel = ad; 
mMtx = mtx;
}
/* --Methode-- */
template <class T> 
POTMatrixAdapter<T>::~POTMatrixAdapter()
{
if (aDel)  delete mMtx;
}
/* --Methode-- */
template <class T> 
double POTMatrixAdapter<T>::Value(int ix, int iy)
{
return((double)(*mMtx)(iy, ix));
}

/* --Methode-- */
template <class T> 
double POTMatrixAdapter<T>::MeanVal(int ix1, int ix2, int jy1, int jy2)
{
int ec;
if (ix1>ix2) { ec=ix1; ix1=ix2; ix2=ec; }
if (jy1>jy2) { ec=jy1; jy1=jy2; jy2=ec; }    
double ss = 0.;
for(int j=jy1; j<=jy2; j++) 
  for(int i=ix1; i<=ix2; i++)  ss += (double)((*mMtx)(j, i));
ss /= (double)((jy2-jy1+1)*(ix2-ix1+1));
return ss;  
}   


/* --Methode-- */
DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
double POTMatrixAdapter< complex<float> >::Value(int ix, int iy)
{
switch (dOpt) {
 case PICDO_Module:
   return _z_mod_( (*mMtx)(iy, ix) );
   break;
 case PICDO_Real:
   return (double)(((*mMtx)(iy, ix)).real());
   break;
 case PICDO_Imag:
   return (double)(((*mMtx)(iy, ix)).imag());
   break;
 case PICDO_Phase:
   return atan2( (double) ((*mMtx)(iy, ix)).real(), (double) ((*mMtx)(iy, ix)).imag() );
   break;
 case PICDO_Module2:
   return _z_mod2_( (*mMtx)(iy, ix) );
   break;
 }
}

/* --Methode-- */
DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
double POTMatrixAdapter< complex<double> >::Value(int ix, int iy)
{
switch (dOpt) {
 case PICDO_Module:
   return _z_mod_( (*mMtx)(iy, ix) );
   break;
 case PICDO_Real:
   return (double)(((*mMtx)(iy, ix)).real());
   break;
 case PICDO_Imag:
   return (double)(((*mMtx)(iy, ix)).imag());
   break;
 case PICDO_Phase:
   return atan2( (double) ((*mMtx)(iy, ix)).real(), (double) ((*mMtx)(iy, ix)).imag() );
   break;
 case PICDO_Module2:
   return _z_mod2_( (*mMtx)(iy, ix) );
   break;
 }
}

/* --Methode-- */
DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
double POTMatrixAdapter< complex<float> >::MeanVal(int ix1, int ix2, int jy1, int jy2)
{
int ec;
if (ix1>ix2) { ec=ix1; ix1=ix2; ix2=ec; }
if (jy1>jy2) { ec=jy1; jy1=jy2; jy2=ec; }    
complex<float> ssz = complex<float>(0.,0.);
for(int j=jy1; j<=jy2; j++) 
  for(int i=ix1; i<=ix2; i++)  ssz += (*mMtx)(j, i);
switch (dOpt) {
 case PICDO_Module:
   return _z_mod_( ssz / complex<float>((float)((jy2-jy1+1)*(ix2-ix1+1)),(float)0.) );
   break;
 case PICDO_Real:
   return ((double)(ssz.real())/(double)((jy2-jy1+1)*(ix2-ix1+1)) ) ;
   break;
 case PICDO_Imag:
   return ((double)(ssz.imag())/(double)((jy2-jy1+1)*(ix2-ix1+1)) ) ;
   break;
 case PICDO_Phase:
   return atan2( (double) ssz.real(), (double) ssz.imag() );
   break;
 case PICDO_Module2:
   return _z_mod2_( ssz / complex<float>((float)((jy2-jy1+1)*(ix2-ix1+1)),(float)0.) );
   break;
 }
return 0.;  
}   

/* --Methode-- */
DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
double POTMatrixAdapter< complex<double> >::MeanVal(int ix1, int ix2, int jy1, int jy2)
{
int ec;
if (ix1>ix2) { ec=ix1; ix1=ix2; ix2=ec; }
if (jy1>jy2) { ec=jy1; jy1=jy2; jy2=ec; }    
complex<double> ssz = complex<double>(0.,0.);
for(int j=jy1; j<=jy2; j++) 
  for(int i=ix1; i<=ix2; i++)  ssz += (*mMtx)(j, i);
switch (dOpt) {
 case PICDO_Module:
   return _z_mod_( ssz / complex<double>((double)((jy2-jy1+1)*(ix2-ix1+1)),(double)0.) );
   break;
 case PICDO_Real:
   return ((double)(ssz.real())/(double)((jy2-jy1+1)*(ix2-ix1+1)) ) ;
   break;
 case PICDO_Imag:
   return ((double)(ssz.imag())/(double)((jy2-jy1+1)*(ix2-ix1+1)) ) ;
   break;
 case PICDO_Phase:
   return atan2( (double) ssz.real(), (double) ssz.imag() );
   break;
 case PICDO_Module2:
   return _z_mod2_( ssz / complex<double>((double)((jy2-jy1+1)*(ix2-ix1+1)),(double)0.) );
   break;
 }
return 0.;  
}   





#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template POTVectorAdapter< uint_1 >
#pragma define_template POTVectorAdapter< uint_2 >
#pragma define_template POTVectorAdapter< uint_4 >
#pragma define_template POTVectorAdapter< uint_8 >
#pragma define_template POTVectorAdapter< int_1 >
#pragma define_template POTVectorAdapter< int_2 >
#pragma define_template POTVectorAdapter< int_4 >
#pragma define_template POTVectorAdapter< int_8 >
#pragma define_template POTVectorAdapter< float >
#pragma define_template POTVectorAdapter< double >
#pragma define_template POTVectorAdapter< complex<float> > 
#pragma define_template POTVectorAdapter< complex<double> > 

#pragma define_template POTMatrixAdapter< uint_1 >
#pragma define_template POTMatrixAdapter< uint_2 >
#pragma define_template POTMatrixAdapter< uint_4 >
#pragma define_template POTMatrixAdapter< uint_8 >
#pragma define_template POTMatrixAdapter< int_1 >
#pragma define_template POTMatrixAdapter< int_2 >
#pragma define_template POTMatrixAdapter< int_4 >
#pragma define_template POTMatrixAdapter< int_8 >
#pragma define_template POTMatrixAdapter< float >
#pragma define_template POTMatrixAdapter< double >
#pragma define_template POTMatrixAdapter< complex<float> > 
#pragma define_template POTMatrixAdapter< complex<double> > 
#endif

#if defined(ANSI_TEMPLATES)
template class POTVectorAdapter< uint_1 >;
template class POTVectorAdapter< uint_2 >;
template class POTVectorAdapter< uint_4 >;
template class POTVectorAdapter< uint_8 >;
template class POTVectorAdapter< int_1 >;
template class POTVectorAdapter< int_2 >;
template class POTVectorAdapter< int_4 >;
template class POTVectorAdapter< int_8 >;
template class POTVectorAdapter< float >;
template class POTVectorAdapter< double >;
template class POTVectorAdapter< complex<float> >;
template class POTVectorAdapter< complex<double> >;

template class POTMatrixAdapter< uint_1 >;
template class POTMatrixAdapter< uint_2 >;
template class POTMatrixAdapter< uint_4 >;
template class POTMatrixAdapter< uint_8 >;
template class POTMatrixAdapter< int_1 >;
template class POTMatrixAdapter< int_2 >;
template class POTMatrixAdapter< int_4 >;
template class POTMatrixAdapter< int_8 >;
template class POTMatrixAdapter< float >;
template class POTMatrixAdapter< double >;
template class POTMatrixAdapter< complex<float> >;
template class POTMatrixAdapter< complex<double> >;
#endif