| 1 |
/* |
| 2 |
* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. |
| 3 |
* |
| 4 |
* The University of Notre Dame grants you ("Licensee") a |
| 5 |
* non-exclusive, royalty free, license to use, modify and |
| 6 |
* redistribute this software in source and binary code form, provided |
| 7 |
* that the following conditions are met: |
| 8 |
* |
| 9 |
* 1. Acknowledgement of the program authors must be made in any |
| 10 |
* publication of scientific results based in part on use of the |
| 11 |
* program. An acceptable form of acknowledgement is citation of |
| 12 |
* the article in which the program was described (Matthew |
| 13 |
* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher |
| 14 |
* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented |
| 15 |
* Parallel Simulation Engine for Molecular Dynamics," |
| 16 |
* J. Comput. Chem. 26, pp. 252-271 (2005)) |
| 17 |
* |
| 18 |
* 2. Redistributions of source code must retain the above copyright |
| 19 |
* notice, this list of conditions and the following disclaimer. |
| 20 |
* |
| 21 |
* 3. Redistributions in binary form must reproduce the above copyright |
| 22 |
* notice, this list of conditions and the following disclaimer in the |
| 23 |
* documentation and/or other materials provided with the |
| 24 |
* distribution. |
| 25 |
* |
| 26 |
* This software is provided "AS IS," without a warranty of any |
| 27 |
* kind. All express or implied conditions, representations and |
| 28 |
* warranties, including any implied warranty of merchantability, |
| 29 |
* fitness for a particular purpose or non-infringement, are hereby |
| 30 |
* excluded. The University of Notre Dame and its licensors shall not |
| 31 |
* be liable for any damages suffered by licensee as a result of |
| 32 |
* using, modifying or distributing the software or its |
| 33 |
* derivatives. In no event will the University of Notre Dame or its |
| 34 |
* licensors be liable for any lost revenue, profit or data, or for |
| 35 |
* direct, indirect, special, consequential, incidental or punitive |
| 36 |
* damages, however caused and regardless of the theory of liability, |
| 37 |
* arising out of the use of or inability to use software, even if the |
| 38 |
* University of Notre Dame has been advised of the possibility of |
| 39 |
* such damages. |
| 40 |
*/ |
| 41 |
|
| 42 |
#include <stdio.h> |
| 43 |
#include <cmath> |
| 44 |
#include <limits> |
| 45 |
#include "math/SphericalHarmonic.hpp" |
| 46 |
#include "utils/simError.h" |
| 47 |
|
| 48 |
using namespace oopse; |
| 49 |
|
| 50 |
SphericalHarmonic::SphericalHarmonic() { |
| 51 |
} |
| 52 |
|
| 53 |
ComplexType SphericalHarmonic::getValueAt(RealType costheta, RealType phi) { |
| 54 |
|
| 55 |
RealType p; |
| 56 |
|
| 57 |
// associated Legendre polynomial |
| 58 |
p = Ptilde(L, M, costheta); |
| 59 |
ComplexType phase(0.0, (RealType)M * phi); |
| 60 |
|
| 61 |
return exp(phase) * (ComplexType)p; |
| 62 |
|
| 63 |
} |
| 64 |
// |
| 65 |
// Routine to calculate the associated Legendre polynomials for m>=0 |
| 66 |
// |
| 67 |
RealType SphericalHarmonic::LegendreP(int l,int m, RealType x) { |
| 68 |
|
| 69 |
RealType temp1, temp2, temp3, temp4, result; |
| 70 |
RealType temp5; |
| 71 |
int i, ll; |
| 72 |
|
| 73 |
if (fabs(x) > 1.0) { |
| 74 |
printf("LegendreP: x out of range: l = %d\tm = %d\tx = %lf\n", l, m, x); |
| 75 |
return std::numeric_limits <RealType>:: quiet_NaN(); |
| 76 |
} |
| 77 |
|
| 78 |
if (m>l) { |
| 79 |
printf("LegendreP: m > l: l = %d\tm = %d\tx = %lf\n", l, m, x); |
| 80 |
return std::numeric_limits <RealType>:: quiet_NaN(); |
| 81 |
} |
| 82 |
|
| 83 |
if (m<0) { |
| 84 |
printf("LegendreP: m < 0: l = %d\tm = %d\tx = %lf\n", l, m, x); |
| 85 |
return std::numeric_limits <RealType>:: quiet_NaN(); |
| 86 |
} else { |
| 87 |
temp3=1.0; |
| 88 |
|
| 89 |
if (m>0) { |
| 90 |
temp1=sqrt(1.0-pow(x,2)); |
| 91 |
temp5 = 1.0; |
| 92 |
for (i=1;i<=m;++i) { |
| 93 |
temp3 *= -temp5*temp1; |
| 94 |
temp5 += 2.0; |
| 95 |
} |
| 96 |
} |
| 97 |
if (l==m) { |
| 98 |
result = temp3; |
| 99 |
} else { |
| 100 |
temp4=x*(2.*m+1.)*temp3; |
| 101 |
if (l==(m+1)) { |
| 102 |
result = temp4; |
| 103 |
} else { |
| 104 |
for (ll=(m+2);ll<=l;++ll) { |
| 105 |
temp2 = (x*(2.*ll-1.)*temp4-(ll+m-1.)*temp3)/(RealType)(ll-m); |
| 106 |
temp3=temp4; |
| 107 |
temp4=temp2; |
| 108 |
} |
| 109 |
result = temp2; |
| 110 |
} |
| 111 |
} |
| 112 |
} |
| 113 |
return result; |
| 114 |
} |
| 115 |
|
| 116 |
|
| 117 |
// |
| 118 |
// Routine to calculate the associated Legendre polynomials for all m... |
| 119 |
// |
| 120 |
RealType SphericalHarmonic::Legendre(int l, int m, RealType x) { |
| 121 |
RealType result; |
| 122 |
if ( m>l || m <-l ) { |
| 123 |
printf("Legendre got a bad argument: l = %d\tm = %d\tx = %lf\n", l, m, x); |
| 124 |
return std::numeric_limits <RealType>:: quiet_NaN(); |
| 125 |
} else if (m >= 0) { |
| 126 |
result = LegendreP(l,m,x); |
| 127 |
} else { |
| 128 |
//result = mpow(-m)*LegendreP(l,-m,x); |
| 129 |
result = mpow(-m)*Fact(l+m)/Fact(l-m)*LegendreP(l, -m, x); |
| 130 |
} |
| 131 |
result *=mpow(m); |
| 132 |
return result; |
| 133 |
} |
| 134 |
// |
| 135 |
// Routine to calculate the normalized associated Legendre polynomials... |
| 136 |
// |
| 137 |
RealType SphericalHarmonic::Ptilde(int l,int m, RealType x){ |
| 138 |
|
| 139 |
RealType result; |
| 140 |
if (m>l || m<-l) { |
| 141 |
result = 0.; |
| 142 |
} else { |
| 143 |
RealType y=(RealType)(2.*l+1.)*Fact(l-m)/Fact(l+m); |
| 144 |
result = mpow(m) * sqrt(y) * Legendre(l,m,x) / sqrt(4.0*M_PI); |
| 145 |
} |
| 146 |
return result; |
| 147 |
} |
| 148 |
// |
| 149 |
// mpow returns (-1)**m |
| 150 |
// |
| 151 |
RealType SphericalHarmonic::mpow(int m) { |
| 152 |
int result; |
| 153 |
if (m<0) m=-m; |
| 154 |
if (m & 0x1) result = -1; |
| 155 |
else result = 1; |
| 156 |
return result; |
| 157 |
} |
| 158 |
// |
| 159 |
// factorial_list is a lookup table for n! |
| 160 |
// |
| 161 |
static RealType factorial_list[171]= |
| 162 |
{1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880.,3.6288e6,3.99168e7,4.790016e8,6.2270208e9, |
| 163 |
8.71782912e10,1.307674368e12,2.0922789888e13,3.55687428096e14,6.402373705728e15, |
| 164 |
1.21645100408832e17, |
| 165 |
2.43290200817664e18,5.109094217170944e19,1.1240007277776077e21,2.585201673888498e22, |
| 166 |
6.204484017332394e23,1.5511210043330986e25,4.0329146112660565e26,1.0888869450418352e28, |
| 167 |
3.0488834461171387e29,8.841761993739702e30,2.6525285981219107e32,8.222838654177922e33, |
| 168 |
2.631308369336935e35,8.683317618811886e36,2.9523279903960416e38,1.0333147966386145e40, |
| 169 |
3.7199332678990125e41,1.3763753091226346e43,5.230226174666011e44,2.0397882081197444e46, |
| 170 |
8.159152832478977e47,3.345252661316381e49,1.40500611775288e51,6.041526306337383e52, |
| 171 |
2.658271574788449e54,1.1962222086548019e56,5.502622159812089e57,2.5862324151116818e59, |
| 172 |
1.2413915592536073e61,6.082818640342675e62,3.0414093201713376e64,1.5511187532873822e66, |
| 173 |
8.065817517094388e67,4.2748832840600255e69,2.308436973392414e71,1.2696403353658276e73, |
| 174 |
7.109985878048635e74,4.0526919504877214e76,2.3505613312828785e78,1.3868311854568984e80, |
| 175 |
8.32098711274139e81,5.075802138772248e83,3.146997326038794e85,1.98260831540444e87, |
| 176 |
1.2688693218588417e89,8.247650592082472e90,5.443449390774431e92,3.647111091818868e94, |
| 177 |
2.4800355424368305e96,1.711224524281413e98,1.1978571669969892e100,8.504785885678623e101, |
| 178 |
6.1234458376886085e103,4.4701154615126844e105,3.307885441519386e107,2.48091408113954e109, |
| 179 |
1.8854947016660504e111,1.4518309202828587e113,1.1324281178206297e115,8.946182130782976e116, |
| 180 |
7.156945704626381e118,5.797126020747368e120,4.753643337012842e122,3.945523969720659e124, |
| 181 |
3.314240134565353e126,2.81710411438055e128,2.4227095383672734e130,2.107757298379528e132, |
| 182 |
1.8548264225739844e134,1.650795516090846e136,1.4857159644817615e138,1.352001527678403e140, |
| 183 |
1.2438414054641308e142,1.1567725070816416e144,1.087366156656743e146,1.032997848823906e148, |
| 184 |
9.916779348709496e149,9.619275968248212e151,9.426890448883248e153,9.332621544394415e155, |
| 185 |
9.332621544394415e157,9.42594775983836e159,9.614466715035127e161,9.90290071648618e163, |
| 186 |
1.0299016745145628e166,1.081396758240291e168,1.1462805637347084e170,1.226520203196138e172, |
| 187 |
1.324641819451829e174,1.4438595832024937e176,1.588245541522743e178,1.7629525510902446e180, |
| 188 |
1.974506857221074e182,2.2311927486598138e184,2.5435597334721877e186,2.925093693493016e188, |
| 189 |
3.393108684451898e190,3.969937160808721e192,4.684525849754291e194,5.574585761207606e196, |
| 190 |
6.689502913449127e198,8.094298525273444e200,9.875044200833601e202,1.214630436702533e205, |
| 191 |
1.506141741511141e207,1.882677176888926e209,2.372173242880047e211,3.0126600184576594e213, |
| 192 |
3.856204823625804e215,4.974504222477287e217,6.466855489220474e219,8.47158069087882e221, |
| 193 |
1.1182486511960043e224,1.4872707060906857e226,1.9929427461615188e228,2.6904727073180504e230, |
| 194 |
3.659042881952549e232,5.012888748274992e234,6.917786472619489e236,9.615723196941089e238, |
| 195 |
1.3462012475717526e241,1.898143759076171e243,2.695364137888163e245,3.854370717180073e247, |
| 196 |
5.5502938327393044e249,8.047926057471992e251,1.1749972043909107e254,1.727245890454639e256, |
| 197 |
2.5563239178728654e258,3.80892263763057e260,5.713383956445855e262,8.62720977423324e264, |
| 198 |
1.3113358856834524e267,2.0063439050956823e269,3.0897696138473508e271,4.789142901463394e273, |
| 199 |
7.471062926282894e275,1.1729568794264145e278,1.853271869493735e280,2.9467022724950384e282, |
| 200 |
4.7147236359920616e284,7.590705053947219e286,1.2296942187394494e289,2.0044015765453026e291, |
| 201 |
3.287218585534296e293,5.423910666131589e295,9.003691705778438e297,1.503616514864999e300, |
| 202 |
2.5260757449731984e302,4.269068009004705e304,7.257415615307999e306}; |
| 203 |
|
| 204 |
// |
| 205 |
// Routine to return the factorial of j |
| 206 |
// |
| 207 |
RealType SphericalHarmonic::Fact(int j) { |
| 208 |
if (j <= 170 && j>=0) return factorial_list[j]; |
| 209 |
|
| 210 |
sprintf( painCave.errMsg, |
| 211 |
"Fact(j) for j >= 171\n"); |
| 212 |
painCave.isFatal = 0; |
| 213 |
simError(); |
| 214 |
return 0.; |
| 215 |
} |
