include/EvtGenBase/EvtComplex.hh

0001
0002 /***********************************************************************
0003 * Copyright 1998-2020 CERN for the benefit of the EvtGen authors       *
0004 *                                                                      *
0005 * This file is part of EvtGen.                                         *
0006 *                                                                      *
0007 * EvtGen is free software: you can redistribute it and/or modify       *
0008 * it under the terms of the GNU General Public License as published by *
0009 * the Free Software Foundation, either version 3 of the License, or    *
0010 * (at your option) any later version.                                  *
0011 *                                                                      *
0012 * EvtGen is distributed in the hope that it will be useful,            *
0013 * but WITHOUT ANY WARRANTY; without even the implied warranty of       *
0014 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        *
0015 * GNU General Public License for more details.                         *
0016 *                                                                      *
0017 * You should have received a copy of the GNU General Public License    *
0018 * along with EvtGen.  If not, see <https://www.gnu.org/licenses/>.     *
0019 ***********************************************************************/
0020
0021 #ifndef EVTCOMPLEX_HH
0022 #define EVTCOMPLEX_HH
0023
0024 #include "EvtGenBase/EvtConst.hh"
0025
0026 #include <iostream>
0027 #include <math.h>
0028
0029 class EvtComplex {
0030     inline friend EvtComplex operator*( double d, const EvtComplex& c );
0031     inline friend EvtComplex operator*( const EvtComplex& c, double d );
0032     inline friend EvtComplex operator/( const EvtComplex& c, double d );
0033     inline friend EvtComplex operator/( double d, const EvtComplex& c );
0034     inline friend EvtComplex operator*( const EvtComplex& c1,
0035                                         const EvtComplex& c2 );
0036     inline friend EvtComplex operator/( const EvtComplex& c1,
0037                                         const EvtComplex& c2 );
0038     inline friend EvtComplex operator+( const EvtComplex& c1,
0039                                         const EvtComplex& c2 );
0040     inline friend EvtComplex operator-( const EvtComplex& c1,
0041                                         const EvtComplex& c2 );
0042     inline friend EvtComplex operator-( const EvtComplex& c );
0043     inline friend EvtComplex conj( const EvtComplex& c );
0044     inline friend double abs( const EvtComplex& c );
0045     inline friend double abs2( const EvtComplex& c );
0046     inline friend double arg( const EvtComplex& c );
0047     inline friend double real( const EvtComplex& c );
0048     inline friend double imag( const EvtComplex& c );
0049     inline friend EvtComplex exp( const EvtComplex& c );
0050     friend std::ostream& operator<<( std::ostream& s, const EvtComplex& c );
0051
0052   public:
0053     EvtComplex() : _rpart( 0.0 ), _ipart( 0.0 ) {}
0054     EvtComplex( double rpart, double ipart = 0.0 ) :
0055         _rpart( rpart ), _ipart( ipart )
0056     {
0057     }
0058     EvtComplex( const EvtComplex& c ) : _rpart( c._rpart ), _ipart( c._ipart )
0059     {
0060     }
0061     inline EvtComplex& operator*=( double d );
0062     inline EvtComplex& operator/=( double d );
0063     EvtComplex& operator*=( EvtComplex c );
0064     EvtComplex& operator/=( EvtComplex c );
0065     inline EvtComplex& operator=( const EvtComplex& c );
0066     inline EvtComplex& operator+=( const EvtComplex& c );
0067     inline EvtComplex& operator-=( const EvtComplex& c );
0068     inline EvtComplex& operator+=( double d );
0069     inline EvtComplex& operator-=( double d );
0070     inline int operator==( const EvtComplex c );
0071     inline int operator!=( const EvtComplex c );
0072
0073   private:
0074     double _rpart, _ipart;
0075 };
0076
0077 typedef EvtComplex* EvtComplexPtr;
0078 typedef EvtComplexPtr* EvtComplexPtrPtr;
0079 typedef EvtComplexPtrPtr* EvtComplexPtrPtrPtr;
0080
0081 EvtComplex& EvtComplex::operator=( const EvtComplex& c )
0082 {
0083     _rpart = c._rpart;
0084     _ipart = c._ipart;
0085
0086     return *this;
0087 }
0088
0089 EvtComplex& EvtComplex::operator+=( const EvtComplex& c )
0090 {
0091     _rpart += c._rpart;
0092     _ipart += c._ipart;
0093
0094     return *this;
0095 }
0096
0097 EvtComplex& EvtComplex::operator-=( const EvtComplex& c )
0098 {
0099     _rpart -= c._rpart;
0100     _ipart -= c._ipart;
0101
0102     return *this;
0103 }
0104
0105 EvtComplex& EvtComplex::operator+=( double d )
0106 {
0107     _rpart += d;
0108
0109     return *this;
0110 }
0111
0112 EvtComplex& EvtComplex::operator-=( double d )
0113 {
0114     _rpart -= d;
0115
0116     return *this;
0117 }
0118
0119 EvtComplex operator*( double d, const EvtComplex& c )
0120 {
0121     return EvtComplex( c._rpart * d, c._ipart * d );
0122 }
0123
0124 EvtComplex operator*( const EvtComplex& c, double d )
0125 {
0126     return EvtComplex( c._rpart * d, c._ipart * d );
0127 }
0128
0129 EvtComplex operator/( const EvtComplex& c, double d )
0130 {
0131     return EvtComplex( c._rpart / d, c._ipart / d );
0132 }
0133
0134 EvtComplex& EvtComplex::operator*=( double d )
0135 {
0136     _rpart *= d;
0137     _ipart *= d;
0138
0139     return *this;
0140 }
0141
0142 EvtComplex& EvtComplex::operator/=( double d )
0143 {
0144     _rpart /= d;
0145     _ipart /= d;
0146
0147     return *this;
0148 }
0149
0150 EvtComplex operator/( double d, const EvtComplex& c )
0151 {
0152     double Num = d / ( c._rpart * c._rpart + c._ipart * c._ipart );
0153
0154     return EvtComplex( Num * c._rpart, -Num * c._ipart );
0155 }
0156
0157 EvtComplex operator/( const EvtComplex& c1, const EvtComplex& c2 )
0158 {
0159     double inv = 1.0 / ( c2._rpart * c2._rpart + c2._ipart * c2._ipart );
0160
0161     return EvtComplex( inv * ( c1._rpart * c2._rpart + c1._ipart * c2._ipart ),
0162                        inv * ( c1._ipart * c2._rpart - c1._rpart * c2._ipart ) );
0163 }
0164
0165 EvtComplex operator*( const EvtComplex& c1, const EvtComplex& c2 )
0166 {
0167     return EvtComplex( c1._rpart * c2._rpart - c1._ipart * c2._ipart,
0168                        c1._rpart * c2._ipart + c1._ipart * c2._rpart );
0169 }
0170
0171 EvtComplex operator-( const EvtComplex& c1, const EvtComplex& c2 )
0172 {
0173     return EvtComplex( c1._rpart - c2._rpart, c1._ipart - c2._ipart );
0174 }
0175
0176 EvtComplex operator+( const EvtComplex& c1, const EvtComplex& c2 )
0177 {
0178     return EvtComplex( c1._rpart + c2._rpart, c1._ipart + c2._ipart );
0179 }
0180
0181 int EvtComplex::operator==( const EvtComplex c )
0182 {
0183     return _rpart == c._rpart && _ipart == c._ipart;
0184 }
0185
0186 int EvtComplex::operator!=( const EvtComplex c )
0187 {
0188     return _rpart != c._rpart || _ipart != c._ipart;
0189 }
0190
0191 EvtComplex operator-( const EvtComplex& c )
0192 {
0193     return EvtComplex( -c._rpart, -c._ipart );
0194 }
0195
0196 EvtComplex conj( const EvtComplex& c )
0197 {
0198     return EvtComplex( c._rpart, -c._ipart );
0199 }
0200
0201 double abs( const EvtComplex& c )
0202 {
0203     double c2 = c._rpart * c._rpart + c._ipart * c._ipart;
0204     if ( c2 <= 0.0 )
0205         return 0.0;
0206     return sqrt( c2 );
0207 }
0208
0209 double abs2( const EvtComplex& c )
0210 {
0211     return c._rpart * c._rpart + c._ipart * c._ipart;
0212 }
0213
0214 double arg( const EvtComplex& c )
0215 {
0216     if ( ( c._rpart == 0 ) && ( c._ipart == 0 ) ) {
0217         return 0.0;
0218     }
0219     if ( c._rpart == 0 ) {
0220         if ( c._ipart > 0 ) {
0221             return EvtConst::pi / 2;
0222         } else {
0223             return -EvtConst::pi / 2;
0224         }
0225     } else {
0226         return atan2( c._ipart, c._rpart );
0227     }
0228 }
0229
0230 double real( const EvtComplex& c )
0231 {
0232     return c._rpart;
0233 }
0234
0235 double imag( const EvtComplex& c )
0236 {
0237     return c._ipart;
0238 }
0239
0240 EvtComplex exp( const EvtComplex& c )
0241 {
0242     return exp( c._rpart ) * EvtComplex( cos( c._ipart ), sin( c._ipart ) );
0243 }
0244
0245 #endif