ScalarFieldMP.h

00001 #ifndef gridripper_multipole_ScalarFieldMP_h
00002 #define gridripper_multipole_ScalarFieldMP_h
00003 
00004 
00005 /*
00006  * Revision may be needed for l>65000. (int range is fully exploited 
00007  * for e.g. l*l.) In this case, int operations, like l*l, 
00008  * should be replaced by GReal_t operations. All >2 powers of 
00009  * integers are currently calculated with floating point (hopefully).
00010  */
00011 
00012 
00013 #include "Bases.h"
00014 #include "Coeff.h"
00015 #include "Gaunt.h"
00016 #include <gridripper/amr1d/FieldComponents.h>
00017 #include <gridripper/tvector.h>
00018 
00019 
00020 namespace gridripper { namespace multipole {
00021 
00022 using namespace gridripper;
00023 using namespace gridripper::amr1d;
00024 
00029 const int SobolevOrderS2=1000000;
00030 
00031 
00039 GReal_t SobolevConstantS2(const int m);
00040 
00041 
00047 GReal_t tailSumS2(const int mup, const int mu, const int L);
00048 
00049 
00058 template <typename K>
00059 class ScalarFieldMP;
00061 template <>
00062 class ScalarFieldMP<GReal_t>;
00064 template <>
00065 class ScalarFieldMP<GComplex_t>;
00066 
00067 
00075 template <>
00076 class ScalarFieldMP<GReal_t> : public FieldComponents<GReal_t>
00077 {
00078     protected:
00080         int maxOrder;
00081     public:
00082         ScalarFieldMP<GReal_t>();
00083         ScalarFieldMP<GReal_t>(const ScalarFieldMP<GReal_t>&);
00088         ScalarFieldMP<GReal_t>(const int lmax);
00095         ScalarFieldMP<GReal_t>(GReal_t*, const int lmax);
00097         ScalarFieldMP<GReal_t>(FieldComponents<GReal_t>*);
00098         ~ScalarFieldMP<GReal_t>();
00100         int getMaxOrder() const;
00102         ScalarFieldMP<GReal_t>& setValue(const GReal_t);
00104         ScalarFieldMP<GReal_t>& annulate();
00106         GReal_t operator()(const int l, const int m) const;
00107         GReal_t& operator()(const int l, const int m);
00109         GReal_t eval(const GReal_t theta, const GReal_t phi) const;
00114         bool resize(int lmax);
00116         ScalarFieldMP<GReal_t>& expand(const tvector<GReal_t>& theta, const tvector<GReal_t>& phi, const tvector< tvector<GReal_t> >& f);
00117 #ifdef HAVE_AGM
00118         ScalarFieldMP<GReal_t>& expand(GReal_t (*f)(const GReal_t theta, const GReal_t phi, void *params), void *params);
00119 #endif
00120 
00121         ScalarFieldMP<GReal_t>& operator=(const ScalarFieldMP<GReal_t>&);
00122         ScalarFieldMP<GReal_t>& operator+=(const ScalarFieldMP<GReal_t>&);
00123         ScalarFieldMP<GReal_t>& operator-=(const ScalarFieldMP<GReal_t>&);
00124         ScalarFieldMP<GReal_t>& operator*=(const GReal_t);
00126         ScalarFieldMP<GReal_t>& operator*=(const ScalarFieldMP<GReal_t>&);
00127         ScalarFieldMP<GReal_t>& operator/=(const GReal_t);
00129         friend class ScalarFieldMP<GComplex_t>;
00130 };
00131 
00133 ScalarFieldMP<GReal_t> operator+(const ScalarFieldMP<GReal_t>&, const ScalarFieldMP<GReal_t>&);
00134 ScalarFieldMP<GReal_t> operator-(const ScalarFieldMP<GReal_t>&, const ScalarFieldMP<GReal_t>&);
00135 ScalarFieldMP<GReal_t> operator-(const ScalarFieldMP<GReal_t>&);
00136 ScalarFieldMP<GReal_t> operator*(const ScalarFieldMP<GReal_t>&, const GReal_t);
00137 ScalarFieldMP<GReal_t> operator*(const GReal_t, const ScalarFieldMP<GReal_t>&);
00138 ScalarFieldMP<GReal_t> operator*(const ScalarFieldMP<GReal_t>&, const ScalarFieldMP<GReal_t>&); // Found some problem on numerical verification.
00139 ScalarFieldMP<GReal_t> operator/(const ScalarFieldMP<GReal_t>&, const GReal_t);
00140 
00147 GReal_t L2S2(const ScalarFieldMP<GReal_t>&, const ScalarFieldMP<GReal_t>&);
00148 
00152 GReal_t H2S2(const ScalarFieldMP<GReal_t>&, const ScalarFieldMP<GReal_t>&, const int mu);
00153 
00158 ScalarFieldMP<GReal_t> reflect(const ScalarFieldMP<GReal_t>&);
00159 
00164 ScalarFieldMP<GReal_t> mulCosTheta(const ScalarFieldMP<GReal_t>&);
00165 
00170 ScalarFieldMP<GReal_t> mulCosThetaSqr(const ScalarFieldMP<GReal_t>&);
00171 
00176 ScalarFieldMP<GReal_t> mulSinThetaSqr(const ScalarFieldMP<GReal_t>&);
00177 
00182 ScalarFieldMP<GReal_t> partial_phi(const ScalarFieldMP<GReal_t>&);
00183 
00188 ScalarFieldMP<GReal_t> Sqr_partial_phi(const ScalarFieldMP<GReal_t>&);
00189 
00195 ScalarFieldMP<GReal_t> sinTheta_partial_theta(const ScalarFieldMP<GReal_t>&);
00196 
00206 ScalarFieldMP<GReal_t> partial_x(const ScalarFieldMP<GReal_t>&);
00207 
00215 ScalarFieldMP<GReal_t> LaplaceS2(const ScalarFieldMP<GReal_t>&);
00216 
00226 GReal_t normBoundNeumann(const GReal_t a00, const GReal_t H2S2aa);
00227 
00236 GReal_t optimalScalingNeumann(const GReal_t a00, const GReal_t H2S2aa);
00237 
00242 ScalarFieldMP<GReal_t> divideByNeumann(const ScalarFieldMP<GReal_t>& a, const ScalarFieldMP<GReal_t>& f, const GReal_t norm, const GReal_t errortolerance);
00243 
00248 GReal_t tailErrorBoundS2(const ScalarFieldMP<GReal_t>&, const int mup, const int mu, int *lmax, int *mmax);
00249 
00250 
00251 
00252 
00260 template <>
00261 class ScalarFieldMP<GComplex_t> : public FieldComponents<GComplex_t>
00262 {
00263     protected:
00265         int maxOrder;
00266     public:
00267         ScalarFieldMP<GComplex_t>();
00268         ScalarFieldMP<GComplex_t>(const ScalarFieldMP<GComplex_t>&);
00270         ScalarFieldMP<GComplex_t>(const ScalarFieldMP<GReal_t>&);
00275         ScalarFieldMP<GComplex_t>(const int lmax);
00282         ScalarFieldMP<GComplex_t>(GReal_t*, const int lmax);
00284         ScalarFieldMP<GComplex_t>(FieldComponents<GComplex_t>*);
00285         ~ScalarFieldMP<GComplex_t>();
00287         int getMaxOrder() const;
00289         ScalarFieldMP<GComplex_t>& setValue(const GComplex_t&);
00291         ScalarFieldMP<GComplex_t>& annulate();
00293         const GComplex_t& operator()(const int l, const int m) const;
00294         GComplex_t& operator()(const int l, const int m);
00296         GComplex_t eval(const GReal_t theta, const GReal_t phi) const;
00301         bool resize(int lmax);
00303         ScalarFieldMP<GComplex_t>& expand(const tvector<GReal_t>& theta, const tvector<GReal_t>& phi, const tvector< tvector<GComplex_t> >& f);
00304 #ifdef HAVE_AGM
00305         ScalarFieldMP<GComplex_t>& expand(GComplex_t (*f)(const GReal_t theta, const GReal_t phi, void *params), void *params);
00306 #endif
00307 
00308         ScalarFieldMP<GComplex_t>& operator=(const ScalarFieldMP<GComplex_t>&);
00309         ScalarFieldMP<GComplex_t>& operator=(const ScalarFieldMP<GReal_t>&);
00310         ScalarFieldMP<GComplex_t>& operator+=(const ScalarFieldMP<GComplex_t>&);
00311         ScalarFieldMP<GComplex_t>& operator+=(const ScalarFieldMP<GReal_t>&);
00312         ScalarFieldMP<GComplex_t>& operator-=(const ScalarFieldMP<GComplex_t>&);
00313         ScalarFieldMP<GComplex_t>& operator-=(const ScalarFieldMP<GReal_t>&);
00314         ScalarFieldMP<GComplex_t>& operator*=(const GReal_t);
00315         ScalarFieldMP<GComplex_t>& operator*=(const GComplex_t&);
00317         ScalarFieldMP<GComplex_t>& operator*=(const ScalarFieldMP<GComplex_t>&);
00319         ScalarFieldMP<GComplex_t>& operator*=(const ScalarFieldMP<GReal_t>&);
00320         ScalarFieldMP<GComplex_t>& operator/=(const GReal_t);
00321         ScalarFieldMP<GComplex_t>& operator/=(const GComplex_t&);
00322 
00324         friend class ScalarFieldMP<GReal_t>;
00325 
00335         static GComplex_t toComplex(const ScalarFieldMP<GReal_t>&, const int l, const int m);
00336         static GComplex_t toComplex(const Coeff<GReal_t>&, const int l, const int m);
00341         static GReal_t toReal(const ScalarFieldMP<GComplex_t>&, const int l, const int m);
00342         static GReal_t toReal(const Coeff<GComplex_t>&, const int l, const int m);
00347         static GReal_t toImag(const ScalarFieldMP<GComplex_t>&, const int l, const int m);
00348         static GReal_t toImag(const Coeff<GComplex_t>&, const int l, const int m);
00353         static GComplex_t toConj(const ScalarFieldMP<GComplex_t>&, const int l, const int m);
00354         static GComplex_t toConj(const Coeff<GComplex_t>&, const int l, const int m);
00355 };
00356 
00358 ScalarFieldMP<GComplex_t> operator+(const ScalarFieldMP<GComplex_t>&, const ScalarFieldMP<GComplex_t>&);
00359 ScalarFieldMP<GComplex_t> operator+(const ScalarFieldMP<GComplex_t>&, const ScalarFieldMP<GReal_t>&);
00360 ScalarFieldMP<GComplex_t> operator+(const ScalarFieldMP<GReal_t>&, const ScalarFieldMP<GComplex_t>&);
00361 ScalarFieldMP<GComplex_t> operator-(const ScalarFieldMP<GComplex_t>&, const ScalarFieldMP<GComplex_t>&);
00362 ScalarFieldMP<GComplex_t> operator-(const ScalarFieldMP<GComplex_t>&, const ScalarFieldMP<GReal_t>&);
00363 ScalarFieldMP<GComplex_t> operator-(const ScalarFieldMP<GReal_t>&, const ScalarFieldMP<GComplex_t>&);
00364 ScalarFieldMP<GComplex_t> operator-(const ScalarFieldMP<GComplex_t>&);
00365 ScalarFieldMP<GComplex_t> operator*(const ScalarFieldMP<GComplex_t>&, const GReal_t);
00366 ScalarFieldMP<GComplex_t> operator*(const GReal_t, const ScalarFieldMP<GComplex_t>&);
00367 ScalarFieldMP<GComplex_t> operator*(const ScalarFieldMP<GComplex_t>&, const GComplex_t&);
00368 ScalarFieldMP<GComplex_t> operator*(const GComplex_t&, const ScalarFieldMP<GComplex_t>&);
00369 ScalarFieldMP<GComplex_t> operator*(const ScalarFieldMP<GReal_t>&, const GComplex_t&);
00370 ScalarFieldMP<GComplex_t> operator*(const GComplex_t&, const ScalarFieldMP<GReal_t>&);
00371 ScalarFieldMP<GComplex_t> operator*(const ScalarFieldMP<GComplex_t>&, const ScalarFieldMP<GComplex_t>&); // Numerically verified by a C++ code.
00372 ScalarFieldMP<GComplex_t> operator*(const ScalarFieldMP<GComplex_t>&, const ScalarFieldMP<GReal_t>&);
00373 ScalarFieldMP<GComplex_t> operator*(const ScalarFieldMP<GReal_t>&, const ScalarFieldMP<GComplex_t>&);
00374 ScalarFieldMP<GComplex_t> operator/(const ScalarFieldMP<GComplex_t>&, const GReal_t);
00375 ScalarFieldMP<GComplex_t> operator/(const ScalarFieldMP<GComplex_t>&, const GComplex_t&);
00376 ScalarFieldMP<GComplex_t> operator/(const ScalarFieldMP<GReal_t>&, const GComplex_t&);
00377 
00382 ScalarFieldMP<GComplex_t> conj(const ScalarFieldMP<GComplex_t>&);
00383 
00385 ScalarFieldMP<GReal_t> real(const ScalarFieldMP<GComplex_t>&);
00386 
00388 ScalarFieldMP<GReal_t> imag(const ScalarFieldMP<GComplex_t>&);
00389 
00396 GComplex_t L2S2(const ScalarFieldMP<GComplex_t>&, const ScalarFieldMP<GReal_t>&);
00397 GComplex_t L2S2(const ScalarFieldMP<GReal_t>&, const ScalarFieldMP<GComplex_t>&);
00398 GComplex_t L2S2(const ScalarFieldMP<GComplex_t>&, const ScalarFieldMP<GComplex_t>&);
00399 
00403 GComplex_t H2S2(const ScalarFieldMP<GComplex_t>&, const ScalarFieldMP<GReal_t>&, const int mu);
00404 GComplex_t H2S2(const ScalarFieldMP<GReal_t>&, const ScalarFieldMP<GComplex_t>&, const int mu);
00405 GComplex_t H2S2(const ScalarFieldMP<GComplex_t>&, const ScalarFieldMP<GComplex_t>&, const int mu);
00406 
00416 GComplex_t L2S2limited(const ScalarFieldMP<GComplex_t>&, const ScalarFieldMP<GComplex_t>&, const GReal_t theta, tvector< tvector< tvector<GReal_t> > > *coeffstorage);
00417 
00422 ScalarFieldMP<GComplex_t> reflect(const ScalarFieldMP<GComplex_t>&);
00423 
00436 ScalarFieldMP<GComplex_t> stepping(const ScalarFieldMP<GComplex_t>&, const int upOrDown /*+1 or -1*/);
00437 
00444 ScalarFieldMP<GComplex_t> mulCosTheta(const ScalarFieldMP<GComplex_t>&);
00445 
00456 ScalarFieldMP<GComplex_t> mulCosThetaSqr(const ScalarFieldMP<GComplex_t>&);
00457 
00462 ScalarFieldMP<GComplex_t> mulSinThetaSqr(const ScalarFieldMP<GComplex_t>&);
00463 
00469 ScalarFieldMP<GComplex_t> mulExpPmIPhi(const ScalarFieldMP<GComplex_t>&, const int sign);
00470 
00475 ScalarFieldMP<GComplex_t> partial_phi(const ScalarFieldMP<GComplex_t>&);
00476 
00481 ScalarFieldMP<GComplex_t> Sqr_partial_phi(const ScalarFieldMP<GComplex_t>&);
00482 
00491 ScalarFieldMP<GComplex_t> sinTheta_partial_theta(const ScalarFieldMP<GComplex_t>&);
00492 
00500 ScalarFieldMP<GComplex_t> partial_theta(const ScalarFieldMP<GComplex_t>&);
00501 
00519 ScalarFieldMP<GComplex_t> partial_x(const ScalarFieldMP<GComplex_t>&);
00520 
00528 ScalarFieldMP<GComplex_t> LaplaceS2(const ScalarFieldMP<GComplex_t>&);
00529 
00539 GReal_t normBoundNeumann(const GComplex_t& a00, const GComplex_t& H2S2aa);
00540 
00549 GComplex_t optimalScalingNeumann(const GComplex_t& a00, const GComplex_t& H2S2aa);
00550 
00568 ScalarFieldMP<GComplex_t> divideByNeumann(const ScalarFieldMP<GComplex_t>& a, const ScalarFieldMP<GReal_t>& f, const GReal_t norm, const GReal_t errortolerance);
00569 ScalarFieldMP<GComplex_t> divideByNeumann(const ScalarFieldMP<GReal_t>& a, const ScalarFieldMP<GComplex_t>& f, const GReal_t norm, const GReal_t errortolerance);
00570 ScalarFieldMP<GComplex_t> divideByNeumann(const ScalarFieldMP<GComplex_t>& a, const ScalarFieldMP<GComplex_t>& f, const GReal_t norm, const GReal_t errortolerance);
00571 
00581 GReal_t tailErrorBoundS2(const ScalarFieldMP<GComplex_t>& a, const int mup, const int mu, int *lmax, int *mmax);
00582 
00583 #ifdef HAVE_AGM
00584 
00588 void maxEval(const unsigned int);
00589 unsigned int maxEval();
00590 void reqAbsError(const GReal_t);
00591 GReal_t reqAbsError();
00592 void reqRelError(const GReal_t);
00593 GReal_t reqRelError();
00594 #endif
00595 
00596 } } // namespace gridripper::multipole
00597 
00598 
00599 #endif /* gridripper_multipole_ScalarFieldMP_h */