| 49 | 
  | 
RealSphericalHarmonic::RealSphericalHarmonic() { | 
| 50 | 
  | 
} | 
| 51 | 
  | 
 | 
| 52 | 
< | 
double RealSphericalHarmonic::getValueAt(double costheta, double phi) { | 
| 52 | 
> | 
RealType RealSphericalHarmonic::getValueAt(RealType costheta, RealType phi) { | 
| 53 | 
  | 
   | 
| 54 | 
< | 
  double p, phase; | 
| 54 | 
> | 
  RealType p, phase; | 
| 55 | 
  | 
   | 
| 56 | 
  | 
  // associated Legendre polynomial | 
| 57 | 
  | 
  p = LegendreP(L,M,costheta); | 
| 58 | 
  | 
  | 
| 59 | 
  | 
  if (functionType == RSH_SIN) { | 
| 60 | 
< | 
    phase = sin((double)M * phi); | 
| 60 | 
> | 
    phase = sin((RealType)M * phi); | 
| 61 | 
  | 
  } else { | 
| 62 | 
< | 
    phase = cos((double)M * phi); | 
| 62 | 
> | 
    phase = cos((RealType)M * phi); | 
| 63 | 
  | 
  } | 
| 64 | 
  | 
   | 
| 65 | 
  | 
  return coefficient*p*phase; | 
| 68 | 
  | 
 | 
| 69 | 
  | 
//---------------------------------------------------------------------------// | 
| 70 | 
  | 
// | 
| 71 | 
< | 
// double LegendreP (int l, int m, double x); | 
| 71 | 
> | 
// RealType LegendreP (int l, int m, RealType x); | 
| 72 | 
  | 
// | 
| 73 | 
  | 
// Computes the value of the associated Legendre polynomial P_lm (x) | 
| 74 | 
  | 
// of order l at a given point. | 
| 81 | 
  | 
//   value of the polynomial in x | 
| 82 | 
  | 
// | 
| 83 | 
  | 
//---------------------------------------------------------------------------// | 
| 84 | 
< | 
double RealSphericalHarmonic::LegendreP (int l, int m, double x) { | 
| 84 | 
> | 
RealType RealSphericalHarmonic::LegendreP (int l, int m, RealType x) { | 
| 85 | 
  | 
  // check parameters | 
| 86 | 
  | 
  if (m < 0 || m > l || fabs(x) > 1.0) { | 
| 87 | 
  | 
    printf("LegendreP got a bad argument: l = %d\tm = %d\tx = %lf\n", l, m, x); | 
| 88 | 
  | 
//    return NAN; | 
| 89 | 
< | 
        return std::numeric_limits <double>:: quiet_NaN(); | 
| 89 | 
> | 
        return std::numeric_limits <RealType>:: quiet_NaN(); | 
| 90 | 
  | 
  } | 
| 91 | 
  | 
   | 
| 92 | 
< | 
  double pmm = 1.0; | 
| 92 | 
> | 
  RealType pmm = 1.0; | 
| 93 | 
  | 
  if (m > 0) { | 
| 94 | 
< | 
    double h = sqrt((1.0-x)*(1.0+x)), | 
| 94 | 
> | 
    RealType h = sqrt((1.0-x)*(1.0+x)), | 
| 95 | 
  | 
      f = 1.0; | 
| 96 | 
  | 
    for (int i = 1; i <= m; i++) { | 
| 97 | 
  | 
      pmm *= -f * h; | 
| 101 | 
  | 
  if (l == m) | 
| 102 | 
  | 
    return pmm; | 
| 103 | 
  | 
  else { | 
| 104 | 
< | 
    double pmmp1 = x * (2 * m + 1) * pmm; | 
| 104 | 
> | 
    RealType pmmp1 = x * (2 * m + 1) * pmm; | 
| 105 | 
  | 
    if (l == (m+1)) | 
| 106 | 
  | 
      return pmmp1; | 
| 107 | 
  | 
    else { | 
| 108 | 
< | 
      double pll = 0.0; | 
| 108 | 
> | 
      RealType pll = 0.0; | 
| 109 | 
  | 
      for (int ll = m+2; ll <= l; ll++) { | 
| 110 | 
  | 
        pll = (x * (2 * ll - 1) * pmmp1 - (ll + m - 1) * pmm) / (ll - m); | 
| 111 | 
  | 
        pmm = pmmp1; |