--- trunk/src/math/SphericalHarmonic.cpp 2006/09/20 22:16:23 1042 +++ trunk/src/math/SphericalHarmonic.cpp 2009/11/25 20:02:06 1390 @@ -6,19 +6,10 @@ * redistribute this software in source and binary code form, provided * that the following conditions are met: * - * 1. Acknowledgement of the program authors must be made in any - * publication of scientific results based in part on use of the - * program. An acceptable form of acknowledgement is citation of - * the article in which the program was described (Matthew - * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher - * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented - * Parallel Simulation Engine for Molecular Dynamics," - * J. Comput. Chem. 26, pp. 252-271 (2005)) - * - * 2. Redistributions of source code must retain the above copyright + * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * - * 3. Redistributions in binary form must reproduce the above copyright + * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the * distribution. @@ -37,14 +28,24 @@ * arising out of the use of or inability to use software, even if the * University of Notre Dame has been advised of the possibility of * such damages. + * + * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your + * research, please cite the appropriate papers when you publish your + * work. Good starting points are: + * + * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). + * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). + * [4] Vardeman & Gezelter, in progress (2009). */ #include +#include #include #include "math/SphericalHarmonic.hpp" #include "utils/simError.h" -using namespace oopse; +using namespace OpenMD; SphericalHarmonic::SphericalHarmonic() { } @@ -52,67 +53,67 @@ ComplexType SphericalHarmonic::getValueAt(RealType cos ComplexType SphericalHarmonic::getValueAt(RealType costheta, RealType phi) { RealType p; - ComplexType phase; - ComplexType I(0.0, 1.0); // associated Legendre polynomial - p = Legendre(L, M, costheta); - - phase = exp(I * (ComplexType)M * (ComplexType)phi); - - return coefficient * phase * (ComplexType)p; + p = Ptilde(L, M, costheta); + ComplexType phase(0.0, (RealType)M * phi); + + return exp(phase) * (ComplexType)p; } - -//---------------------------------------------------------------------------// // -// RealType LegendreP (int l, int m, RealType x); +// Routine to calculate the associated Legendre polynomials for m>=0 // -// Computes the value of the associated Legendre polynomial P_lm (x) -// of order l at a given point. -// -// Input: -// l = degree of the polynomial >= 0 -// m = parameter satisfying 0 <= m <= l, -// x = point in which the computation is performed, range -1 <= x <= 1. -// Returns: -// value of the polynomial in x -// -//---------------------------------------------------------------------------// -RealType SphericalHarmonic::LegendreP (int l, int m, RealType x) { - // check parameters - if (m < 0 || m > l || fabs(x) > 1.0) { - printf("LegendreP got a bad argument: l = %d\tm = %d\tx = %lf\n", l, m, x); +RealType SphericalHarmonic::LegendreP(int l,int m, RealType x) { + + RealType temp1, temp2, temp3, temp4, result; + RealType temp5; + int i, ll; + + if (fabs(x) > 1.0) { + printf("LegendreP: x out of range: l = %d\tm = %d\tx = %lf\n", l, m, x); return std::numeric_limits :: quiet_NaN(); } - RealType pmm = 1.0; - if (m > 0) { - RealType h = sqrt((1.0-x)*(1.0+x)), - f = 1.0; - for (int i = 1; i <= m; i++) { - pmm *= -f * h; - f += 2.0; - } + if (m>l) { + printf("LegendreP: m > l: l = %d\tm = %d\tx = %lf\n", l, m, x); + return std::numeric_limits :: quiet_NaN(); } - if (l == m) - return pmm; - else { - RealType pmmp1 = x * (2 * m + 1) * pmm; - if (l == (m+1)) - return pmmp1; - else { - RealType pll = 0.0; - for (int ll = m+2; ll <= l; ll++) { - pll = (x * (2 * ll - 1) * pmmp1 - (ll + m - 1) * pmm) / (ll - m); - pmm = pmmp1; - pmmp1 = pll; + + if (m<0) { + printf("LegendreP: m < 0: l = %d\tm = %d\tx = %lf\n", l, m, x); + return std::numeric_limits :: quiet_NaN(); + } else { + temp3=1.0; + + if (m>0) { + temp1=sqrt(1.0-pow(x,2)); + temp5 = 1.0; + for (i=1;i<=m;++i) { + temp3 *= -temp5*temp1; + temp5 += 2.0; } - return pll; } + if (l==m) { + result = temp3; + } else { + temp4=x*(2.*m+1.)*temp3; + if (l==(m+1)) { + result = temp4; + } else { + for (ll=(m+2);ll<=l;++ll) { + temp2 = (x*(2.*ll-1.)*temp4-(ll+m-1.)*temp3)/(RealType)(ll-m); + temp3=temp4; + temp4=temp2; + } + result = temp2; + } + } } + return result; } + // // Routine to calculate the associated Legendre polynomials for all m... // @@ -124,12 +125,27 @@ RealType SphericalHarmonic::Legendre(int l, int m, Rea } else if (m >= 0) { result = LegendreP(l,m,x); } else { + //result = mpow(-m)*LegendreP(l,-m,x); result = mpow(-m)*Fact(l+m)/Fact(l-m)*LegendreP(l, -m, x); } result *=mpow(m); return result; } // +// Routine to calculate the normalized associated Legendre polynomials... +// +RealType SphericalHarmonic::Ptilde(int l,int m, RealType x){ + + RealType result; + if (m>l || m<-l) { + result = 0.; + } else { + RealType y=(RealType)(2.*l+1.)*Fact(l-m)/Fact(l+m); + result = mpow(m) * sqrt(y) * Legendre(l,m,x) / sqrt(4.0*M_PI); + } + return result; +} +// // mpow returns (-1)**m // RealType SphericalHarmonic::mpow(int m) {