CAAM : Среда разработки и моделирования

00001 #ifndef BxMathH
00002 #define BxMathH
00003 
00004 #include <complex>
00005 #include <exception>
00006 #include "BxCAAMInterpretorData.h"
00007 #include "BxArray.h"
00008 #define _USE_MATH_DEFINES
00009 #include "math.h"
00010 #include "BxComplex.h"
00011 
00012 #ifndef M_PI
00013 #define M_PI       3.14159265358979323846
00014 #endif
00015 
00016 typedef std::complex<double> complex_double;
00017 typedef std::complex<long double> complex_long_double;
00018 
00019 namespace Bx{
00020 
00021   class exception : public std::exception{
00022   private:
00023     const char * mMessage;
00024   public:
00025     exception(const char * message) : std::exception(), mMessage(message) {}
00026 
00027     virtual const char * what () const//  _RWSTD_THROW_SPEC_NULL
00028     { 
00029       return mMessage;
00030     }
00031   };
00032 #if defined(__BORLANDC__)
00033     #pragma warn -8027
00034 #endif
00035   class CComplex1DArray : public CArray<complex_double>
00036   {
00037   public:
00038     CComplex1DArray(const TComplex1D & complex1D) :
00039       CArray<complex_double>(complex1D.size)
00040     {
00041       for (int i = 0; i < mSize; i++)
00042       {
00043         (*this)[i] = complex_double(complex1D.real[i], complex1D.imag[i]);
00044       }
00045     }
00046   };
00047 
00048   template <class T>
00049   class CModesArray : public CArray<T>
00050   {
00051   /*public:
00052     const PBxDouble alpha;
00053     const PBxDouble w;
00054     const PBxInteger m;
00055     const PBxInteger n;
00056   */
00057   private:
00058 
00059     virtual void updateArray(const TModes & modes){
00060       for (int i = 0; i < mSize; i++)
00061       {
00062         (*this)[i] = modes.A[i] * (cos(modes.phi[i] * M_PI) +
00063           i * sin(modes.phi[i]*M_PI));
00064       }
00065     }
00066   public:
00067     CModesArray(const TModes & modes) :
00068       CArray<T>(modes.count)//, alpha(modes.alpha), w(modes.w), m(modes.m), n(modes.n)
00069     {
00070       updateArray(modes);
00071     }
00072   };
00073 
00074   template<>
00075   void CModesArray<complex_long_double>::updateArray(const TModes & modes)
00076     {
00077       for (int i = 0; i < mSize; i++)
00078       {
00079         (*this)[i] = (long double)modes.A[i] * (cosl(modes.phi[i] * M_PI) +
00080           i * sinl(modes.phi[i]*M_PI));
00081       }
00082     }
00083 
00084   template <typename T>
00085   class Complex {
00086   public:
00087     static const std::complex<T> i;
00088   };
00089 
00090   template <typename T>
00091   const std::complex<T> Complex<T>::i(0, 1); 
00092 #if defined(__BORLANDC__)
00093     #pragma warn +8027
00094 #endif
00095 } // namespace Bx
00096 
00097   //---------------------------------------------------------------------------
00098   template <class T>
00099   std::complex<T> gammaFunctionln(const std::complex<T>& xx)
00100   {
00101     std::complex<T> x,y,tmp,ser;
00102     static T cof[6]={76.18009172947146,-86.50532032941677,
00103         24.01409824083091,-1.231739572450155,
00104         0.1208650973866179e-2,-0.5395239384953e-5};
00105     int j;
00106 
00107     y=x=xx;
00108     tmp = x + (T)5.5;
00109     tmp -= (x+(T)0.5)*log(tmp);
00110     ser=1.000000000190015;
00111     for (j=0;j<=5;j++){
00112           ser += cof[j]/y;
00113           y += 1.0;
00114         }
00115     return -tmp+log((T)2.5066282746310005*ser/x);
00116   }
00117   //---------------------------------------------------------------------------
00118   template <class T>
00119   inline std::complex<T> gammaFunction(std::complex<T> x){
00120     std::complex<T> result = exp(gammaFunctionln(x+(T)1));
00121     return result;
00122   }
00123   //---------------------------------------------------------------------------
00124   template <class T>
00125   std::complex<T> gammaFunctionlnl(const std::complex<T>& xx)
00126   {
00127     std::complex<T> x,y,tmp,ser;
00128     static T cof[6]={76.18009172947146,-86.50532032941677,
00129         24.01409824083091,-1.231739572450155,
00130         0.1208650973866179e-2,-0.5395239384953e-5};
00131     int j;
00132 
00133     y=x=xx;
00134     tmp = x + (T)5.5;
00135     tmp -= (x+(T)0.5)*logl(tmp);
00136     ser=1.000000000190015;
00137     for (j=0;j<=5;j++){
00138           ser += cof[j]/y;
00139           y += 1.0;
00140         }
00141     return -tmp+logl((T)2.5066282746310005*ser/x);
00142   }
00143   //---------------------------------------------------------------------------
00144   template <class T>
00145   inline std::complex<T> gammaFunctionl(std::complex<T> x){
00146     std::complex<T> result = expl(gammaFunctionlnl(x+(T)1));
00147     return result;
00148   }
00149   //---------------------------------------------------------------------------
00150 
00151 
00152     /* таблица коэффициентов */
00153 #define CN 8
00154 static double cof[CN]={
00155        2.5066282746310005,
00156        1.0000000000190015,
00157        76.18009172947146,
00158       -86.50532032941677,
00159        24.01409824083091,
00160       -1.231739572450155,
00161        0.1208650973866179e-2,
00162       -0.5395239384953e-5,
00163 };
00164 
00165 /* логарифм гамма-фукнции */
00166 template <class T>
00167   std::complex<T>  GammLn(std::complex<T> x) {
00168   std::complex<T> y,ser;
00169   double *co;
00170   int j;
00171   /* вычисоение последовательностей */
00172   ser=(T)cof[1]; y=x; co=cof+2;
00173   for(j=2;j<CN;j++) {
00174     y+=1.; ser+=(T)(*co)/y; co++;
00175   }
00176   /* и других частей функции */
00177   y=x+(T)5.5;
00178   y-=(x+(T)0.5)*logl(y);
00179   return(-y+logl((T)cof[0]*ser/x));
00180 }
00181 
00182 /* сама гамма-функция */
00183 template <class T>
00184   std::complex<T>  Gamm(const std::complex<T> & x) {
00185   return(expl(GammLn(x)));
00186 }
00187 
00188 /* бета-функция */
00189 template <class T>
00190 double Beta(const std::complex<T> &  x, const std::complex<T> &  y) {
00191   return(exp(GammLn(x)+GammLn(y)-GammLn(x+y)));
00192 }
00193 
00194   template <class T>
00195   T factorial(T x){
00196     T result = 1;
00197     for (T i = 2; i <= x; i++)
00198       result *= i;
00199     return result;
00200   }
00201   //---------------------------------------------------------------------------
00202   // Функция производит инверсное копирование радиальносимметричной матрицы
00203   //---------------------------------------------------------------------------
00204   template <class T>
00205   void copyRadialDevelopment(T & matrix, int size){
00206     // инверсное копирование треугольника
00207     for (int i = 1; i <= size / 2; i++)
00208       for (int j = 0; j < i; j++)
00209         matrix[i][j] = matrix[j][i];
00210     // отражающее копирование по вертикали
00211     for (int i = 1; i < size / 2; i++)
00212       for (int j = 0;j <= size / 2;j++)
00213         matrix[size - i][j] = matrix[i][j];
00214     // отражающее копирование по горизонтали
00215     for (int i = 0;i < size; i++)
00216       for (int j = 1; j < size / 2; j++)
00217         matrix[i][size - j] = matrix[i][j];
00218   }
00219 
00220   //---------------------------------------------------------------------------
00221   // Функция производит инверсное копирование эллемента радиальносимметричной матрицы
00222   //---------------------------------------------------------------------------
00223   template <class T>
00224   inline void copyRadialDevelopment(T & matrix, int size, int i, int j){
00225     matrix[j][i] = matrix[i][j];
00226     matrix[size - i][j] = matrix[i][j];
00227     matrix[size - j][i] = matrix[i][j];
00228     matrix[i][size - j] = matrix[i][j];
00229     matrix[j][size - i] = matrix[i][j];
00230     matrix[size - i][size - j] = matrix[i][j];
00231     matrix[size - j][size - i] = matrix[i][j];
00232   }
00233 
00234 #endif // BxMathH
BxMath.h
