Matematica.h

#ifndef MATEMATICA_H
#define MATEMATICA_H
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdarg.h>
#include <fstream>
#include <iostream>
//#include <time.h>
//#include <complex>
//to look up!
//#include <boost/math/special_functions/fpclassify.hpp>
#include <assert.h>
#ifndef __GSL__
#include "randomlib.h"
#include "MatematicaVect.h"
#include "MatematicaQuadri.h"
#include "MatematicaMatrix.h"
#include "MatematicaStruct.h"
#include "MatematicaPlane.h"
#include <mt19937ar.h>
#else
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#endif
#ifndef FFTW3_H
#include  <fftw3.h>
#endif
#define STRSIZE 512
#define STRLEN  STRSIZE
#define QUAD(A) ((A)*(A))
#define SQUARE(A) QUAD(A)
#define SQR(A) QUAD(A)
#define CUB(A) (QUAD(A)*(A))
#define CUBE(A) (QUAD(A)*(A))
#define SIGN(A) (A>0?1.:-1.)
#define NTAB 32
#define DUE_PI (6.2831853071795864769252867665590057683943)
#define DEG_RAD (0.0174532925199433)
#define RAD_DEG (57.2957795130823)
#define INV_DUE_PI (0.159154943091895)
#define PI M_PI
#define SQRT_PI 1.77245385090552
#define MIN(a,b) ( ( (a)<=(b) ) ? (a) : (b) )
#define MIN3(a,b,c) ( MIN(MIN(a,b),c) )
#define MAS(a,b) ( ( (a)>=(b) ) ? (a):(b) )
#define MAX(a,b) MAS(a,b)
#define POS(a)   ( ( a>=0 ) ? (a):(-a) )
#define ASS(a)   POS(a)
#define ABS(a)   POS(a)
#define OP_IF(c,t) ( ((c)&(t))==(t) )
#define OP_INVERT     0x00001
#define OP_ROT_90     0x00002
#define OP_ROT_180    0x00004
#define OP_MIRROR     0x00008
#define OP_ROT_270    0x0000c
#define OP_TRANSPOSE  0x0000f

using namespace std;

//---------------------Miscellaneous------------------
int NumeroStringa(const char *NumeroC);
void abort_(const char * s, ...);
void SigErr(int IfWrong,const char * s, ...);
void Message(const char * s, ...);
double rnor();
//float nfix(void);
//float efix(void);
float zigset(unsigned long);
/* double Evalx(void (*pt2Function)(double x)){ */
/*   return pt2Function(.2); */
/* } */
//-------------------Matematica------------------------
class Matematica{
 private:
  SPLINE SpMem;
  double PrecMinimo;
  int NPassi;
#ifdef __GSL__
  gsl_rng *rng;
#endif
 public:
  Matematica();
  ~Matematica();
  //-----------------------Filter----------------------------
  int PermuteRandomAll(int *Sequence,int NMass);
  int PermuteRandomAll(PERMUTE *Sequence,int NMass);
  int ApplyFilter(Matrice *Point,Matrice *Res,Matrice *Mask);
  int ApplyFilter(Matrice *Res,Matrice *Mask);
  int Transform(int *Out,int *In,int NEdge,int operation);
  void BackFold(Matrice *In,Matrice *Out,int NShift);
  void Smooth(double *st,double *sw,int NIn,int NOut);
  //------------------------Func------------------------
  typedef double(Matematica::*FUNC)(double x);
  FUNC Func;
  double Evalx(double x){return (*this.*Func)(x);}
  double ContactAngle(double x);
  double fProva(double x);
  /* /// Boh */
  /* double *fp; */
  double F(double TD,double T);
  double Df(double x,double Delta);
  //template <class Elabora> Elabora ElabFunc(double *st,double *sw,int NMass);
  //template <class Elabora> Elabora Elaboro(double *st,double *sw,int NMass){return (*this.*ElabFunc)(st,sw,NMass);};
  // NPunti = Elabora<int>(st,sw,NMass);
  typedef void(Matematica::*ELAB)(double *st,double *sw,int NMass);
  ELAB Elab;
  void ElabSt(double *st,double *sw,int NMass){return (*this.*Elab)(st,sw,NMass);};
  void Derivata(double *st,double *sw,int NMass);
  void DerO4(double *st,double *sw,int NMass);
  double Integrazione(double *Punti,double *sw,int NMass);
  void IntegraA3();
  void SquareGradient(double *st,double *sw,int NMass);
  double Integrazione(double a,double b);
  double IntegrazioneGauss(double a,double b,double Scarto);
  int Zeri(double a,double b,double *Radici,int NRadici);
  RADICE RegulaFalsi(double a,double b);
  RADICE Newton(double a);
  double Estremo(double a,double b);
  double Fattoriale(int n);
  double Gamma(int n);
  double Elevato(double x,int Volte);
  double Bessel(double Val,int Ord);
  double Neumann(double Val,int Ord);
  double Segno(int n);
  double QuasiBessel(double Val,int Ord);
  double QuasiNeumann(double Val,int Ord);
  double WeightFunction(double x,double a);
  double WeightFunction2(double x,double a);
  double LJHamaker(double r,double r_np,double theta);
  double LJHamakerCum(double Rad,double RadNpMin,double RadNpMax);
  double LJHamaker(double Rad,double r_np);
  double LJ39(double r,double r_np);
  //-----------------------Sign---------------------
  void ExecCommand(double *st,double *st1,int NMass,char *cmd);
  double ExecFormula(double x,double y,char *cmd);
  double ExecFormula(double **st,int n,char *cmd);
  double Gauss(double Media,double Scarto,double x);
  bool InizializzaGaussiano(double Scarto,int N);
  double Gaussiano(double Media,double Scarto);
  void Spettro(double *st,double *sw,int NMass);
  void Spettro2d(double *st,double *sw,int NMass);
  void Spettro2d(double *st,double **sw,int NMass);
  void SpettroDFT(double *st,double *sw,int NMass);
  void Radice(double *st,double *sw,int N);
  void Autocor(bool *st,double *sAuto,int N);
  void Autocor(double *st,double *sAuto,int NMass);
  double Norm(double *st,int NMass);
  int NormalizeArea(double *st,int NMass);
  void NormalizeVect(double *st,int NMass);
  int Normalizza(double *st,int NMass);
  int  Normalizza(double *st,double *sw,int NMass);
  void Modulo(double *st,double *sw,int NMass);
  void MediaMobile(double *st,int NMass,double *sw,int Parti);
  int MediaMobile(double *st,int NMass,double *sw,double *sErr,int Parti);
  int CorrelaDuePunti(double *st,int NMass,double *sw,int Punti);
  void Autosimilarita(double *st,int NMass,double *sw,int Valori);
  MOMENTI Distribuzione(const double *st,int NMass);
  MOMENTI Distribuzione(const double *st,int NMass,double *Intervalli,int Valori,int IfNorm);
  MOMENTI Distribuzione(const double *st,int NMass,double *Intervalli,int Valori,double *Confine,int IfNorm);
  MOMENTI DistrErr(const double *st,int NMass,double *Intervalli,double *Err,int Valori,double *Confine,int IfNorm);
  MOMENTI DistribuzioneGauss(const double *st,int NMass,double *Intervalli,double *dInt,int Valori,int IfNorm);
  MOMENTI DistribuzioneMaxwell(const double *st,int NMass,double *Intervalli,double *dInt,int Valori,int IfNorm);
  void DistrSample(double *Px,double *Py,int NMax,double **Distr,int NBin,const int NSample,int IfNorm,double *xBound);
  MOMENTI WeightAverage(const double *sx,const double *sy,int NMax);
  void WeightHisto(double **hist,double *Border,int NBin,int NHisto,double tolerance,double *OrPos,double *kSpring);
  void Sort(double *Sign,int NMass);
  void Swap(int i,int j,double *Sign);
  void Swap(double *s,int si,double *t,int ti,const int NDim);
  void Sort(int *Sign,int NMass);
  void Swap(int i,int j,int *Sign);
  void FileSin1d(char *FName);
  void FileSin2d(char *FName);
  void ConvWeight(double *st,int NMax,double *sw,int *WIndex,int NWeight);
  void FillWeightGauss(double *st,int *WIndex,int NWeight,double CutOff,double Sigma);
  //--------------------Interp--------------------------------
  double LinInterp(double Px1,double Px2,double Py1,double Py2,double x);
  RETTA InterRett(double *Px,double *Py,int NMass);
  RETTA InterExp(double *Px,double *Py,int NMass);
  MOMENTI InterGauss(double *Px,double *Py,int NMass);
  RETTA InterRett(double *Px,double *Py,double *Peso,int NMass);
  PARABOLA MinimoParabola(double a, double b,double *Px,double *Py,int NMass);
  PARABOLA MinimoParabola(double *Px,double *Py,int NMass);
  SPLINE Parab(double *P1,double *P2,double *P3,int x,int y);
  SPLINE Parab2(double *PA,double *PB,double *PC,int x,int y);
  CIRCLE Osculante(double *PA,double *PB,double *PC,int x,int y);
  SPLINE Cubica(double *PA,double *PB,double *PC,double *PD,int x,int y);
  SPLINE Forth(double *PA,double *PB,double *PC,double *PD,double *PE,int x,int y);
  SPLINE Spline3(double *P1,double *P2,double *P3,int x,int y);
  SPLINE Spline3Beg(double *P1,double *P2,double *P3,int x,int y);  
  SPLINE Spline3End(double *P1,double *P2,int x,int y);  
  SPLINE Spline4Beg(double *P1,double *P2,double *P3,double *P4,int x,int y);
  SPLINE Spline4(double *P1,double *P2,double *P3,double *P4,int x,int y);
  SPLINE Spline4(double *P1,double *P2,double *P3,int x,int y);
  SPLINE Spline4PreEnd(double *P1,double *P2,double *P3,int x,int y);
  SPLINE Spline4End(double *P1,double *P2,int x,int y);
  int Polinomio(double *P1,double *P2,int NMass,Spline *Sp);
  int DerMatrix(double *Px,double *Py,int NMass,SPLINE Wg,Spline *Sp);
  double Casuale();
  double RandDiscrProb(double *Prob,int NBin);
  double QBezier(double *P1,double *P2,double *P3,double x,int y);
  int Factorial(int times);
  double Binomial(int times,int n);
  double Blend(const double *dPoint,double x,int nPoint,int nOrder);
  double Blend(double *dPoint,size_t Incr,double x,int nPoint,int nOrder);
  int InterBSpline2D(Matrice *MaIn,Matrice *MaOut);
  int Voronoi();
  //---------------------------Algebra------------------------------
  //  int Invert(Matrice A,Matrice B);
  double Ypsilon;
  double PreFact;
};
//-----------------------Char---------------------
uint sTable(char *String);
uint sTable(char Char);
//--------------------Trials---------------------
class Funtore{
 public:
  virtual double operator()(double x)=0;
  virtual double Call(double x)=0;
};
template <class Classe> class SpecFuntore : public Funtore{
 private:
  double (Classe::*ftp)(const double x);
  Classe* pt2Object;
 public:
  SpecFuntore(Classe* _pt2Object,void (Classe::*_ftp)(double x))
    { pt2Object = _pt2Object; ftp = _ftp; };
  virtual double operator()(const double x)
  { return (*pt2Object.*ftp)(x);};
  virtual double Call(double x)
  { return (*pt2Object.*ftp)(x);};
};

#endif //MATEMATICA_H