Xmipp  v3.23.11-Nereus
Classes | Macros | Functions
numerical_tools.cpp File Reference
#include <vector>
#include "numerical_tools.h"
#include "core/matrix2d.h"
#include "core/numerical_recipes.h"
#include "core/xmipp_funcs.h"
#include "core/xmipp_memory.h"
Include dependency graph for numerical_tools.cpp:

Go to the source code of this file.

Classes

struct  CDAB
 

Macros

#define Element(a, b, c)   a[b*nDim+c]
 
#define RowVector(a, b)   (&a[b*nDim])
 
#define CopyVector(a, b)   memcpy((a),(b),nDim*sizeof(double))
 

Functions

void randomPermutation (int N, MultidimArray< int > &result)
 
void powellOptimizer (Matrix1D< double > &p, int i0, int n, double(*f)(double *x, void *), void *prm, double ftol, double &fret, int &iter, const Matrix1D< double > &steps, bool show)
 
double solveNonNegative (const Matrix2D< double > &C, const Matrix1D< double > &d, Matrix1D< double > &result)
 
void solveViaCholesky (const Matrix2D< double > &A, const Matrix1D< double > &b, Matrix1D< double > &result)
 
void evaluateQuadratic (const Matrix1D< double > &x, const Matrix1D< double > &c, const Matrix2D< double > &H, double &val, Matrix1D< double > &grad)
 
void quadraticProgramming_obj32 (int nparam, int j, double *x, double *fj, void *cd)
 
void quadraticProgramming_cntr32 (int nparam, int j, double *x, double *gj, void *cd)
 
void quadraticProgramming_grob32 (int nparam, double *x, double *gradfj, void *cd)
 
void quadraticProgramming_grcn32 (int nparam, int j, double *gradgj, void *cd)
 
void quadraticProgramming (const Matrix2D< double > &C, const Matrix1D< double > &d, const Matrix2D< double > &A, const Matrix1D< double > &b, const Matrix2D< double > &Aeq, const Matrix1D< double > &beq, Matrix1D< double > &bl, Matrix1D< double > &bu, Matrix1D< double > &x)
 
void leastSquare (const Matrix2D< double > &C, const Matrix1D< double > &d, const Matrix2D< double > &A, const Matrix1D< double > &b, const Matrix2D< double > &Aeq, const Matrix1D< double > &beq, Matrix1D< double > &bl, Matrix1D< double > &bu, Matrix1D< double > &x)
 
void regularizedLeastSquare (const Matrix2D< double > &A, const Matrix1D< double > &d, double lambda, const Matrix2D< double > &G, Matrix1D< double > &x)
 
double checkRandomness (const std::string &sequence)
 
template<int L1, int L2>
double ZernikeSphericalHarmonics (int n, int m, double xr, double yr, double zr, double r)
 
double ZernikeSphericalHarmonics (int l1, int n, int l2, int m, double xr, double yr, double zr, double r)
 
void spherical_index2lnm (int idx, int &l1, int &n, int &l2, int &m, int max_l1)
 
int spherical_lnm2index (int l1, int n, int l2, int m, int max_l1)
 
template<typename T >
void solveBySVD (const Matrix2D< T > &A, const Matrix1D< T > &b, Matrix1D< double > &result, double tolerance)
 
template<typename T >
void solve (const Matrix2D< T > &A, const Matrix2D< T > &b, Matrix2D< T > &result)
 
template void solveBySVD< double > (Matrix2D< double > const &, Matrix1D< double > const &, Matrix1D< double > &, double)
 
template void solve< double > (Matrix2D< double > const &, Matrix2D< double > const &, Matrix2D< double > &)
 

Macro Definition Documentation

◆ CopyVector

#define CopyVector (   a,
  b 
)    memcpy((a),(b),nDim*sizeof(double))

Definition at line 385 of file numerical_tools.cpp.

◆ Element

#define Element (   a,
  b,
  c 
)    a[b*nDim+c]

Definition at line 383 of file numerical_tools.cpp.

◆ RowVector

#define RowVector (   a,
  b 
)    (&a[b*nDim])

Definition at line 384 of file numerical_tools.cpp.

Function Documentation

◆ checkRandomness()

double checkRandomness ( const std::string &  sequence)

Check the randomness of a sequence. The function returns a z-score of randomness. The highest the Z-score in absolute value, the less random the sequence is. The sequence is supposed to be formed by two symbols:

  • and -, 0 and 1, A and B, ...

Definition at line 806 of file numerical_tools.cpp.

807 {
808  int imax=sequence.size();
809  if (imax<=1)
810  return 0;
811  double n0=1;
812  double n1=0;
813  double R=1;
814  int current=0;
815  for (int i=1; i<imax; ++i)
816  {
817  if (sequence[i]!=sequence[i-1])
818  {
819  current=(current+1)%2;
820  R++;
821  if (current==0)
822  n0++;
823  else
824  n1++;
825  }
826  else
827  {
828  if (current==0)
829  n0++;
830  else
831  n1++;
832  }
833  }
834  double m=1+2*n0*n1/(n0+n1);
835  double s=sqrt(2*n0*n1*(2*n0*n1-n0-n1)/((n0+n1)*(n0+n1)*(n0+n1-1)));
836  double z=(R-m)/s;
837  return ABS(z);
838 }
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
i
#define i
Definition: numerical_recipes.cpp:2493
ABS
#define ABS(x)
Definition: xmipp_macros.h:142
z
double z
Definition: numerical_recipes.cpp:8490
m
int m
Definition: image_find_center.cpp:57

◆ powellOptimizer()

void powellOptimizer ( Matrix1D< double > &  p,
int  i0,
int  n,
double(*)(double *x, void *)  f,
void *  prm,
double  ftol,
double &  fret,
int &  iter,
const Matrix1D< double > &  steps,
bool  show 
)

Definition at line 44 of file numerical_tools.cpp.

48 {
49  // Adapt indexes of p
50  double *pptr = p.adaptForNumericalRecipes();
51  double *auxpptr = pptr + (i0 - 1);
52 
53  // Form direction matrix
54  std::vector<double> buffer(n*n);
55  auto *xi= buffer.data()-1;
56  for (int i = 1, ptr = 1; i <= n; i++)
57  for (int j = 1; j <= n; j++, ptr++)
58  xi[ptr] = (i == j) ? steps(i - 1) : 0;
59 
60  // Optimize
61  powell(auxpptr, xi -n, n, ftol, iter, fret, f, prm, show); // xi - n because NR works with matrices starting at [1,1]
62 }
xi
double xi
Definition: numerical_recipes.cpp:4635
buffer
HBITMAP buffer
Definition: svm-toy.cpp:37
iter
glob_prnt iter
Definition: numerical_recipes.cpp:4698
i
#define i
Definition: numerical_recipes.cpp:2493
f
double * f
Definition: numerical_recipes.cpp:4414
powell
void powell(double *p, double *xi, int n, double ftol, int &iter, double &fret, double(*func)(double *, void *), void *prm, bool show)
Definition: numerical_recipes.cpp:877
j
#define j
Definition: numerical_recipes.cpp:2493
steps
double steps
Definition: numerical_recipes.cpp:5166
prm
ProgClassifyCL2D * prm
Definition: mpi_classify_CL2D.cpp:35
Matrix1D::adaptForNumericalRecipes
T * adaptForNumericalRecipes() const
Definition: matrix1d.h:844
n
int * n
Definition: numerical_recipes.cpp:2229

◆ quadraticProgramming_cntr32()

void quadraticProgramming_cntr32 ( int  nparam,
int  j,
double *  x,
double *  gj,
void *  cd 
)

Definition at line 182 of file numerical_tools.cpp.

183 {
184  auto* in = (CDAB *)cd;
185  *gj = 0;
186  for (int k = 0; k < nparam; k++)
187  *gj += in->A(j - 1, k) * x[k];
188  *gj -= in->B(j - 1, 0);
189 }
cd
void * cd
Definition: numerical_recipes.cpp:4416
nparam
static void nparam
Definition: numerical_recipes.cpp:5153
x
doublereal * x
Definition: numerical_recipes.cpp:2230
k
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
in
int in
Definition: image_find_center.cpp:56
j
#define j
Definition: numerical_recipes.cpp:2493

◆ quadraticProgramming_grcn32()

void quadraticProgramming_grcn32 ( int  nparam,
int  j,
double *  gradgj,
void *  cd 
)

Definition at line 206 of file numerical_tools.cpp.

207 {
208  auto* in = (CDAB *)cd;
209  for (int k = 0; k < nparam; k++)
210  gradgj[k] = in->A(j - 1, k);
211 }
cd
void * cd
Definition: numerical_recipes.cpp:4416
nparam
static void nparam
Definition: numerical_recipes.cpp:5153
k
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
in
int in
Definition: image_find_center.cpp:56
j
#define j
Definition: numerical_recipes.cpp:2493

◆ quadraticProgramming_grob32()

void quadraticProgramming_grob32 ( int  nparam,
double *  x,
double *  gradfj,
void *  cd 
)

Definition at line 192 of file numerical_tools.cpp.

193 {
194  auto* in = (CDAB *)cd;
195  Matrix2D<double> X(1,nparam);
196  for (int i=0; i<nparam; ++i)
197  X(0,i)=x[i];
198 
199  Matrix2D<double> gradient;
200  gradient = in->C * X + in->D;
201  for (int k = 0; k < nparam; k++)
202  gradfj[k] = gradient(k, 0);
203 }
cd
void * cd
Definition: numerical_recipes.cpp:4416
nparam
static void nparam
Definition: numerical_recipes.cpp:5153
x
doublereal * x
Definition: numerical_recipes.cpp:2230
i
#define i
Definition: numerical_recipes.cpp:2493
k
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
in
int in
Definition: image_find_center.cpp:56

◆ quadraticProgramming_obj32()

void quadraticProgramming_obj32 ( int  nparam,
int  j,
double *  x,
double *  fj,
void *  cd 
)

Definition at line 169 of file numerical_tools.cpp.

170 {
171  auto* in = (CDAB *)cd;
172  Matrix2D<double> X(nparam,1);
173  for (int i=0; i<nparam; ++i)
174  X(i,0)=x[i];
175  Matrix2D<double> result;
176  result = 0.5 * X.transpose() * in->C * X + in->D.transpose() * X;
177 
178  *fj = result(0, 0);
179 }
cd
void * cd
Definition: numerical_recipes.cpp:4416
nparam
static void nparam
Definition: numerical_recipes.cpp:5153
Matrix2D::transpose
Matrix2D< T > transpose() const
Definition: matrix2d.cpp:1314
x
doublereal * x
Definition: numerical_recipes.cpp:2230
i
#define i
Definition: numerical_recipes.cpp:2493
in
int in
Definition: image_find_center.cpp:56

◆ solve< double >()

template void solve< double > ( Matrix2D< double > const &  ,
Matrix2D< double > const &  ,
Matrix2D< double > &   
)

◆ solveBySVD< double >()

template void solveBySVD< double > ( Matrix2D< double > const &  ,
Matrix1D< double > const &  ,
Matrix1D< double > &  ,
double   
)

◆ spherical_index2lnm()

void spherical_index2lnm ( int  idx,
int &  l1,
int &  n,
int &  l2,
int &  m,
int  maxl1 
)

Index to Spherical harmonics index. Given an integer consecutive index (0,1,2,3,...) this function returns the corresponding (n,l,m) for the spherical harmonics basis.

Definition at line 2159 of file numerical_tools.cpp.

2160 {
2161  auto numR = static_cast<int>(std::floor((4+4*max_l1+std::pow(max_l1,2))/4));
2162  double aux_id = std::floor(idx-(idx/numR)*numR);
2163  l1 = static_cast<int>(std::floor((1.0+std::sqrt(1.0+4.0*aux_id))/2.0) + std::floor((2.0+2.0*std::sqrt(aux_id))/2.0) - 2.0);
2164  n = static_cast<int>(std::ceil((4.0*aux_id - l1*(l1+2.0))/2.0));
2165  l2 = static_cast<int>(std::floor(std::sqrt(std::floor(idx/numR))));
2166  m = static_cast<int>(std::floor(idx/numR)-l2*(l2+1));
2167 }
floor
__host__ __device__ float2 floor(const float2 v)
Definition: cuda_basic_math.h:175
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
m
int m
Definition: image_find_center.cpp:57
n
int * n
Definition: numerical_recipes.cpp:2229

◆ spherical_lnm2index()

int spherical_lnm2index ( int  l1,
int  n,
int  l2,
int  m,
int  maxl1 
)

Index to Spherical harmonics index. Given the corresponding (n,l,m), this function returns an integer consecutive index (0,1,2,3,...) for the spherical harmonics basis.

Definition at line 2169 of file numerical_tools.cpp.

2170 {
2171  auto numR = static_cast<int>(std::floor((4+4*max_l1+std::pow(max_l1,2))/4));
2172  int id_SH = l2*(l2+1)+m;
2173  auto id_R = static_cast<int>(std::floor((2*n + l1*(l1 + 2))/4));
2174  int id_Z = id_SH*numR+id_R;
2175  return id_Z;
2176 }
floor
__host__ __device__ float2 floor(const float2 v)
Definition: cuda_basic_math.h:175
m
int m
Definition: image_find_center.cpp:57
n
int * n
Definition: numerical_recipes.cpp:2229

◆ ZernikeSphericalHarmonics() [1/2]

template<int L1, int L2>
double ZernikeSphericalHarmonics ( int  n,
int  m,
double  xr,
double  yr,
double  zr,
double  r 
)

Definition at line 1054 of file numerical_tools.cpp.

1055 {
1056  // General variables
1057  double r2=r*r;
1058  double xr2=xr*xr;
1059  double yr2=yr*yr;
1060  double zr2=zr*zr;
1061 
1062  //Variables needed for l>=5
1063  double tht=0.0;
1064  double phi=0.0;
1065  double costh=0.0;
1066  double sinth=0.0;
1067  double costh2=0.0;
1068  double sinth2=0.0;
1069  double cosph=0.0;
1070  double cosph2=0.0;
1071  if (L2>=5)
1072  {
1073  tht = atan2(yr,xr);
1074  phi = atan2(zr,sqrt(xr2 + yr2));
1075  sinth = sin(abs(m)*phi); costh = cos(tht); cosph = cos(abs(m)*phi);
1076  sinth2 = sinth*sinth; costh2 = costh*costh; cosph2 = cosph * cosph;
1077  }
1078 
1079  // Zernike polynomial
1080  double R=0.0;
1081 
1082  switch (L1)
1083  {
1084  case 0:
1085  R = std::sqrt(3.0);
1086  break;
1087  case 1:
1088  R = std::sqrt(5.0)*r;
1089  break;
1090  case 2:
1091  switch (n)
1092  {
1093  case 0:
1094  R = -0.5*std::sqrt(7.0)*(2.5*(1-2*r2)+0.5);
1095  break;
1096  case 2:
1097  R = std::sqrt(7.0)*r2;
1098  break;
1099  } break;
1100  case 3:
1101  switch (n)
1102  {
1103  case 1:
1104  R = -1.5*r*(3.5*(1-2*r2)+1.5);
1105  break;
1106  case 3:
1107  R = 3.0*r2*r;
1108  } break;
1109  case 4:
1110  switch (n)
1111  {
1112  case 0:
1113  R = std::sqrt(11.0)*((63.0*r2*r2/8.0)-(35.0*r2/4.0)+(15.0/8.0));
1114  break;
1115  case 2:
1116  R = -0.5*std::sqrt(11.0)*r2*(4.5*(1-2*r2)+2.5);
1117  break;
1118  case 4:
1119  R = std::sqrt(11.0)*r2*r2;
1120  break;
1121  } break;
1122  case 5:
1123  switch (n)
1124  {
1125  case 1:
1126  R = std::sqrt(13.0)*r*((99.0*r2*r2/8.0)-(63.0*r2/4.0)+(35.0/8.0));
1127  break;
1128  case 3:
1129  R = -0.5*std::sqrt(13.0)*r2*r*(5.5*(1-2*r2)+3.5);
1130  break;
1131  case 5:
1132  R = std::sqrt(13.0)*r2*r2*r;
1133  break;
1134  } break;
1135  case 6:
1136  switch (n)
1137  {
1138  case 0:
1139  R = 103.8 * r2 * r2 * r2 - 167.7 * r2 * r2 + 76.25 * r2 - 8.472;
1140  break;
1141  case 2:
1142  R = 69.23 * r2 * r2 * r2 - 95.86 * r2 * r2 + 30.5 * r2;
1143  break;
1144  case 4:
1145  R = 25.17 * r2 * r2 * r2 - 21.3 * r2 * r2;
1146  break;
1147  case 6:
1148  R = 3.873 * r2 * r2 * r2;
1149  break;
1150  } break;
1151  case 7:
1152  switch (n)
1153  {
1154  case 1:
1155  R = 184.3 * r2 * r2 * r2 * r - 331.7 * r2 * r2 * r + 178.6 * r2 * r - 27.06 * r;
1156  break;
1157  case 3:
1158  R = 100.5 * r2 * r2 * r2 * r - 147.4 * r2 * r2 * r + 51.02 * r2 * r;
1159  break;
1160  case 5:
1161  R = 30.92 * r2 * r2 * r2 * r - 26.8 * r2 * r2 * r;
1162  break;
1163  case 7:
1164  R = 4.123 * r2 * r2 * r2 * r;
1165  break;
1166  } break;
1167  case 8:
1168  switch (n)
1169  {
1170  case 0:
1171  R = 413.9*r2*r2*r2*r2 - 876.5*r2*r2*r2 + 613.6*r2*r2 - 157.3*r2 + 10.73;
1172  break;
1173  case 2:
1174  R = 301.0*r2*r2*r2*r2 - 584.4*r2*r2*r2 + 350.6*r2*r2 - 62.93*r2;
1175  break;
1176  case 4:
1177  R = 138.9*r2*r2*r2*r2 - 212.5*r2*r2*r2 + 77.92*r2*r2;
1178  break;
1179  case 6:
1180  R = 37.05*r2*r2*r2*r2 - 32.69*r2*r2*r2;
1181  break;
1182  case 8:
1183  R = 4.359*r2*r2*r2*r2;
1184  break;
1185  } break;
1186  case 9:
1187  switch (n)
1188  {
1189  case 1:
1190  R = 751.6*r2*r2*r2*r2*r - 1741.0*r2*r2*r2*r + 1382.0*r2*r2*r - 430.0*r2*r + 41.35*r;
1191  break;
1192  case 3:
1193  R = 462.6*r2*r2*r2*r2*r - 949.5*r2*r2*r2*r + 614.4*r2*r2*r - 122.9*r2*r;
1194  break;
1195  case 5:
1196  R = 185.0*r2*r2*r2*r2*r - 292.1*r2*r2*r2*r + 111.7*r2*r2*r;
1197  break;
1198  case 7:
1199  R = 43.53*r2*r2*r2*r2*r - 38.95*r2*r2*r2*r;
1200  break;
1201  case 9:
1202  R = 4.583*r2*r2*r2*r2*r;
1203  break;
1204  } break;
1205  case 10:
1206  switch (n)
1207  {
1208  case 0:
1209  R = 1652.0*r2*r2*r2*r2*r2 - 4326.0*r2*r2*r2*r2 + 4099.0*r2*r2*r2 - 1688.0*r2*r2 + 281.3*r2 - 12.98;
1210  break;
1211  case 2:
1212  R = 1271.0*r2*r2*r2*r2*r2 - 3147.0*r2*r2*r2*r2 + 2732.0*r2*r2*r2 - 964.4*r2*r2 + 112.5*r2;
1213  break;
1214  case 4:
1215  R = 677.7*r2*r2*r2*r2*r2 - 1452.0*r2*r2*r2*r2 + 993.6*r2*r2*r2 - 214.3*r2*r2;
1216  break;
1217  case 6:
1218  R = 239.2*r2*r2*r2*r2*r2 - 387.3*r2*r2*r2*r2 + 152.9*r2*r2*r2;
1219  break;
1220  case 8:
1221  R = 50.36*r2*r2*r2*r2*r2 - 45.56*r2*r2*r2*r2;
1222  break;
1223  case 10:
1224  R = 4.796*r2*r2*r2*r2*r2;
1225  break;
1226  } break;
1227  case 11:
1228  switch (n)
1229  {
1230  case 1:
1231  R = r*-5.865234375E+1+(r*r*r)*8.7978515625E+2-(r*r*r*r*r)*4.2732421875E+3+(r*r*r*r*r*r*r)*9.0212890625E+3-(r*r*r*r*r*r*r*r*r)*8.61123046875E+3+pow(r,1.1E+1)*3.04705078125E+3;
1232  break;
1233  case 3:
1234  R = (r*r*r)*2.513671875E+2-(r*r*r*r*r)*1.89921875E+3+(r*r*r*r*r*r*r)*4.920703125E+3-(r*r*r*r*r*r*r*r*r)*5.29921875E+3+pow(r,1.1E+1)*2.0313671875E+3;
1235  break;
1236  case 5:
1237  R = (r*r*r*r*r)*-3.453125E+2+(r*r*r*r*r*r*r)*1.5140625E+3-(r*r*r*r*r*r*r*r*r)*2.1196875E+3+pow(r,1.1E+1)*9.559375E+2;
1238  break;
1239  case 7:
1240  R = (r*r*r*r*r*r*r)*2.01875E+2-(r*r*r*r*r*r*r*r*r)*4.9875E+2+pow(r,1.1E+1)*3.01875E+2;
1241  break;
1242  case 9:
1243  R = (r*r*r*r*r*r*r*r*r)*-5.25E+1+pow(r,1.1E+1)*5.75E+1;
1244  break;
1245  case 11:
1246  R = pow(r,1.1E+1)*5.0;
1247  break;
1248  } break;
1249  case 12:
1250  switch (n)
1251  {
1252  case 0:
1253  R = (r*r)*-4.57149777110666E+2+(r*r*r*r)*3.885773105442524E+3-(r*r*r*r*r*r)*1.40627979054153E+4+(r*r*r*r*r*r*r*r)*2.460989633446932E+4-pow(r,1.0E+1)*2.05828223888278E+4+pow(r,1.2E+1)*6.597058457955718E+3+1.523832590368693E+1;
1254  break;
1255  case 2:
1256  R = (r*r)*-1.828599108443595E+2+(r*r*r*r)*2.220441774539649E+3-(r*r*r*r*r*r)*9.375198603600264E+3+(r*r*r*r*r*r*r*r)*1.789810642504692E+4-pow(r,1.0E+1)*1.583294029909372E+4+pow(r,1.2E+1)*5.277646766364574E+3;
1257  break;
1258  case 4:
1259  R = (r*r*r*r)*4.934315054528415E+2-(r*r*r*r*r*r)*3.409163128584623E+3+(r*r*r*r*r*r*r*r)*8.260664503872395E+3-pow(r,1.0E+1)*8.444234826177359E+3+pow(r,1.2E+1)*3.104498097866774E+3;
1260  break;
1261  case 6:
1262  R = (r*r*r*r*r*r)*-5.244866351671517E+2+(r*r*r*r*r*r*r*r)*2.202843867704272E+3-pow(r,1.0E+1)*2.98031817394495E+3+pow(r,1.2E+1)*1.307157093837857E+3;
1263  break;
1264  case 8:
1265  R = (r*r*r*r*r*r*r*r)*2.591581020820886E+2-pow(r,1.0E+1)*6.274354050420225E+2+pow(r,1.2E+1)*3.734734553815797E+2;
1266  break;
1267  case 10:
1268  R = pow(r,1.0E+1)*-5.975575286115054E+1+pow(r,1.2E+1)*6.49519052838441E+1;
1269  break;
1270  case 12:
1271  R = pow(r,1.2E+1)*5.19615242270811;
1272  break;
1273  } break;
1274  case 13:
1275  switch (n)
1276  {
1277  case 1:
1278  R = r*7.896313435467891E+1-(r*r*r)*1.610847940832376E+3+(r*r*r*r*r)*1.093075388422608E+4-(r*r*r*r*r*r*r)*3.400678986203671E+4+(r*r*r*r*r*r*r*r*r)*5.332882955634594E+4-pow(r,1.1E+1)*4.102217658185959E+4+pow(r,1.3E+1)*1.230665297454596E+4;
1279  break;
1280  case 3:
1281  R = (r*r*r)*-4.602422688100487E+2+(r*r*r*r*r)*4.858112837433815E+3-(r*r*r*r*r*r*r)*1.854915810656548E+4+(r*r*r*r*r*r*r*r*r)*3.281774126553535E+4-pow(r,1.1E+1)*2.734811772125959E+4+pow(r,1.3E+1)*8.687049158513546E+3;
1282  break;
1283  case 5:
1284  R = (r*r*r*r*r)*8.832932431697845E+2-(r*r*r*r*r*r*r)*5.707433263555169E+3+(r*r*r*r*r*r*r*r*r)*1.312709650617838E+4-pow(r,1.1E+1)*1.286970245704055E+4+pow(r,1.3E+1)*4.572131136059761E+3;
1285  break;
1286  case 7:
1287  R = (r*r*r*r*r*r*r)*-7.60991101808846E+2+(r*r*r*r*r*r*r*r*r)*3.088728589691222E+3-pow(r,1.1E+1)*4.064116565383971E+3+pow(r,1.3E+1)*1.741764242306352E+3;
1288  break;
1289  case 9:
1290  R = (r*r*r*r*r*r*r*r*r)*3.251293252306059E+2-pow(r,1.1E+1)*7.741174410264939E+2+pow(r,1.3E+1)*4.54373280601576E+2;
1291  break;
1292  case 11:
1293  R = pow(r,1.1E+1)*-6.731456008902751E+1+pow(r,1.3E+1)*7.269972489634529E+1;
1294  break;
1295  case 13:
1296  R = pow(r,1.3E+1)*5.385164807128604;
1297  break;
1298  } break;
1299  case 14:
1300  switch (n)
1301  {
1302  case 0:
1303  R = (r*r)*6.939451623205096E+2-(r*r*r*r)*7.910974850460887E+3+(r*r*r*r*r*r)*3.955487425231934E+4-(r*r*r*r*r*r*r*r)*1.010846786448956E+5+pow(r,1.0E+1)*1.378427436065674E+5-pow(r,1.2E+1)*9.542959172773361E+4+pow(r,1.4E+1)*2.63567443819046E+4-1.749441585683962E+1;
1304  break;
1305  case 2:
1306  R = (r*r)*2.775780649287626E+2-(r*r*r*r)*4.520557057410479E+3+(r*r*r*r*r*r)*2.636991616821289E+4-(r*r*r*r*r*r*r*r)*7.351612992358208E+4+pow(r,1.0E+1)*1.060328796973228E+5-pow(r,1.2E+1)*7.634367338204384E+4+pow(r,1.4E+1)*2.170555419689417E+4;
1307  break;
1308  case 4:
1309  R = (r*r*r*r)*-1.004568234980106E+3+(r*r*r*r*r*r)*9.589060424804688E+3-(r*r*r*r*r*r*r*r)*3.393052150321007E+4+pow(r,1.0E+1)*5.655086917197704E+4-pow(r,1.2E+1)*4.490804316592216E+4+pow(r,1.4E+1)*1.370877107170224E+4;
1310  break;
1311  case 6:
1312  R = (r*r*r*r*r*r)*1.475240065354854E+3-(r*r*r*r*r*r*r*r)*9.04813906750083E+3+pow(r,1.0E+1)*1.99591302959919E+4-pow(r,1.2E+1)*1.8908649754107E+4+pow(r,1.4E+1)*6.527986224621534E+3;
1313  break;
1314  case 8:
1315  R = (r*r*r*r*r*r*r*r)*-1.064486949119717E+3+pow(r,1.0E+1)*4.201922167569399E+3-pow(r,1.2E+1)*5.402471358314157E+3+pow(r,1.4E+1)*2.270603904217482E+3;
1316  break;
1317  case 10:
1318  R = pow(r,1.0E+1)*4.001830635787919E+2-pow(r,1.2E+1)*9.395602362267673E+2+pow(r,1.4E+1)*5.44944937011227E+2;
1319  break;
1320  case 12:
1321  R = pow(r,1.2E+1)*-7.516481889830902E+1+pow(r,1.4E+1)*8.073258326109499E+1;
1322  break;
1323  case 14:
1324  R = pow(r,1.4E+1)*5.567764362829621;
1325  break;
1326  } break;
1327  case 15:
1328  switch (n)
1329  {
1330  case 1:
1331  R = r*-1.022829477079213E+2+(r*r*r)*2.720726409032941E+3-(r*r*r*r*r)*2.448653768128157E+4+(r*r*r*r*r*r*r)*1.042945123462677E+5-(r*r*r*r*r*r*r*r*r)*2.370329826049805E+5+pow(r,1.1E+1)*2.953795629386902E+5-pow(r,1.3E+1)*1.903557183384895E+5+pow(r,1.5E+1)*4.958846444106102E+4;
1332  break;
1333  case 3:
1334  R = (r*r*r)*7.773504025805742E+2-(r*r*r*r*r)*1.088290563613176E+4+(r*r*r*r*r*r*r)*5.688791582524776E+4-(r*r*r*r*r*r*r*r*r)*1.458664508337975E+5+pow(r,1.1E+1)*1.969197086257935E+5-pow(r,1.3E+1)*1.343687423563004E+5+pow(r,1.5E+1)*3.653886853551865E+4;
1335  break;
1336  case 5:
1337  R = (r*r*r*r*r)*-1.978710115659982E+3+(r*r*r*r*r*r*r)*1.750397410005331E+4-(r*r*r*r*r*r*r*r*r)*5.834658033359051E+4+pow(r,1.1E+1)*9.266809817695618E+4-pow(r,1.3E+1)*7.072039071393013E+4+pow(r,1.5E+1)*2.08793534488678E+4;
1338  break;
1339  case 7:
1340  R = (r*r*r*r*r*r*r)*2.333863213345408E+3-(r*r*r*r*r*r*r*r*r)*1.372860713732243E+4+pow(r,1.1E+1)*2.926360995060205E+4-pow(r,1.3E+1)*2.694110122436285E+4+pow(r,1.5E+1)*9.077979760378599E+3;
1341  break;
1342  case 9:
1343  R = (r*r*r*r*r*r*r*r*r)*-1.445116540770978E+3+pow(r,1.1E+1)*5.57402094297111E+3-pow(r,1.3E+1)*7.028113362878561E+3+pow(r,1.5E+1)*2.90495352332294E+3;
1344  break;
1345  case 11:
1346  R = pow(r,1.1E+1)*4.846974733015522E+2-pow(r,1.3E+1)*1.124498138058931E+3+pow(r,1.5E+1)*6.455452274046838E+2;
1347  break;
1348  case 13:
1349  R = pow(r,1.3E+1)*-8.329615837475285E+1+pow(r,1.5E+1)*8.904072102135979E+1;
1350  break;
1351  case 15:
1352  R = pow(r,1.5E+1)*5.744562646534177;
1353  break;
1354  } break;
1355  default: break;
1356  }
1357 
1358  // Spherical harmonic
1359  double Y=0.0;
1360 
1361  switch (L2)
1362  {
1363  case 0:
1364  Y = (1.0/2.0)*sqrt(1.0/PI);
1365  break;
1366  case 1:
1367  switch (m)
1368  {
1369  case -1:
1370  Y = sqrt(3.0/(4.0*PI))*yr;
1371  break;
1372  case 0:
1373  Y = sqrt(3.0/(4.0*PI))*zr;
1374  break;
1375  case 1:
1376  Y = sqrt(3.0/(4.0*PI))*xr;
1377  break;
1378  } break;
1379  case 2:
1380  switch (m)
1381  {
1382  case -2:
1383  Y = sqrt(15.0/(4.0*PI))*xr*yr;
1384  break;
1385  case -1:
1386  Y = sqrt(15.0/(4.0*PI))*zr*yr;
1387  break;
1388  case 0:
1389  Y = sqrt(5.0/(16.0*PI))*(-xr2-yr2+2.0*zr2);
1390  break;
1391  case 1:
1392  Y = sqrt(15.0/(4.0*PI))*xr*zr;
1393  break;
1394  case 2:
1395  Y = sqrt(15.0/(16.0*PI))*(xr2-yr2);
1396  break;
1397  } break;
1398  case 3:
1399  switch (m)
1400  {
1401  case -3:
1402  Y = sqrt(35.0/(16.0*2.0*PI))*yr*(3.0*xr2-yr2);
1403  break;
1404  case -2:
1405  Y = sqrt(105.0/(4.0*PI))*zr*yr*xr;
1406  break;
1407  case -1:
1408  Y = sqrt(21.0/(16.0*2.0*PI))*yr*(4.0*zr2-xr2-yr2);
1409  break;
1410  case 0:
1411  Y = sqrt(7.0/(16.0*PI))*zr*(2.0*zr2-3.0*xr2-3.0*yr2);
1412  break;
1413  case 1:
1414  Y = sqrt(21.0/(16.0*2.0*PI))*xr*(4.0*zr2-xr2-yr2);
1415  break;
1416  case 2:
1417  Y = sqrt(105.0/(16.0*PI))*zr*(xr2-yr2);
1418  break;
1419  case 3:
1420  Y = sqrt(35.0/(16.0*2.0*PI))*xr*(xr2-3.0*yr2);
1421  break;
1422  } break;
1423  case 4:
1424  switch (m)
1425  {
1426  case -4:
1427  Y = sqrt((35.0*9.0)/(16.0*PI))*yr*xr*(xr2-yr2);
1428  break;
1429  case -3:
1430  Y = sqrt((9.0*35.0)/(16.0*2.0*PI))*yr*zr*(3.0*xr2-yr2);
1431  break;
1432  case -2:
1433  Y = sqrt((9.0*5.0)/(16.0*PI))*yr*xr*(7.0*zr2-(xr2+yr2+zr2));
1434  break;
1435  case -1:
1436  Y = sqrt((9.0*5.0)/(16.0*2.0*PI))*yr*zr*(7.0*zr2-3.0*(xr2+yr2+zr2));
1437  break;
1438  case 0:
1439  Y = sqrt(9.0/(16.0*16.0*PI))*(35.0*zr2*zr2-30.0*zr2+3.0);
1440  break;
1441  case 1:
1442  Y = sqrt((9.0*5.0)/(16.0*2.0*PI))*xr*zr*(7.0*zr2-3.0*(xr2+yr2+zr2));
1443  break;
1444  case 2:
1445  Y = sqrt((9.0*5.0)/(8.0*8.0*PI))*(xr2-yr2)*(7.0*zr2-(xr2+yr2+zr2));
1446  break;
1447  case 3:
1448  Y = sqrt((9.0*35.0)/(16.0*2.0*PI))*xr*zr*(xr2-3.0*yr2);
1449  break;
1450  case 4:
1451  Y = sqrt((9.0*35.0)/(16.0*16.0*PI))*(xr2*(xr2-3.0*yr2)-yr2*(3.0*xr2-yr2));
1452  break;
1453  } break;
1454  case 5:
1455  switch (m)
1456  {
1457  case -5:
1458  Y = (3.0/16.0)*sqrt(77.0/(2.0*PI))*sinth2*sinth2*sinth*sin(5.0*phi);
1459  break;
1460  case -4:
1461  Y = (3.0/8.0)*sqrt(385.0/(2.0*PI))*sinth2*sinth2*sin(4.0*phi);
1462  break;
1463  case -3:
1464  Y = (1.0/16.0)*sqrt(385.0/(2.0*PI))*sinth2*sinth*(9.0*costh2-1.0)*sin(3.0*phi);
1465  break;
1466  case -2:
1467  Y = (1.0/4.0)*sqrt(1155.0/(4.0*PI))*sinth2*(3.0*costh2*costh-costh)*sin(2.0*phi);
1468  break;
1469  case -1:
1470  Y = (1.0/8.0)*sqrt(165.0/(4.0*PI))*sinth*(21.0*costh2*costh2-14.0*costh2+1)*sin(phi);
1471  break;
1472  case 0:
1473  Y = (1.0/16.0)*sqrt(11.0/PI)*(63.0*costh2*costh2*costh-70.0*costh2*costh+15.0*costh);
1474  break;
1475  case 1:
1476  Y = (1.0/8.0)*sqrt(165.0/(4.0*PI))*sinth*(21.0*costh2*costh2-14.0*costh2+1)*cos(phi);
1477  break;
1478  case 2:
1479  Y = (1.0/4.0)*sqrt(1155.0/(4.0*PI))*sinth2*(3.0*costh2*costh-costh)*cos(2.0*phi);
1480  break;
1481  case 3:
1482  Y = (1.0/16.0)*sqrt(385.0/(2.0*PI))*sinth2*sinth*(9.0*costh2-1.0)*cos(3.0*phi);
1483  break;
1484  case 4:
1485  Y = (3.0/8.0)*sqrt(385.0/(2.0*PI))*sinth2*sinth2*cos(4.0*phi);
1486  break;
1487  case 5:
1488  Y = (3.0/16.0)*sqrt(77.0/(2.0*PI))*sinth2*sinth2*sinth*cos(5.0*phi);
1489  break;
1490  }break;
1491  case 6:
1492  switch (m)
1493  {
1494  case -6:
1495  Y = -0.6832*sinth*pow(costh2 - 1.0, 3);
1496  break;
1497  case -5:
1498  Y = 2.367*costh*sinth*pow(1.0 - 1.0*costh2, 2.5);
1499  break;
1500  case -4:
1501  Y = 0.001068*sinth*(5198.0*costh2 - 472.5)*pow(costh2 - 1.0, 2);
1502  break;
1503  case -3:
1504  Y = -0.005849*sinth*pow(1.0 - 1.0*costh2, 1.5)*(- 1732.0*costh2*costh + 472.5*costh);
1505  break;
1506  case -2:
1507  Y = -0.03509*sinth*(costh2 - 1.0)*(433.1*costh2*costh2 - 236.2*costh2 + 13.12);
1508  break;
1509  case -1:
1510  Y = 0.222*sinth*pow(1.0 - 1.0*costh2, 0.5)*(86.62*costh2*costh2*costh - 78.75*costh2*costh + 13.12*costh);
1511  break;
1512  case 0:
1513  Y = 14.68*costh2*costh2*costh2 - 20.02*costh2*costh2 + 6.675*costh2 - 0.3178;
1514  break;
1515  case 1:
1516  Y = 0.222*cosph*pow(1.0 - 1.0*costh2, 0.5)*(86.62*costh2*costh2*costh - 78.75*costh2*costh + 13.12*costh);
1517  break;
1518  case 2:
1519  Y = -0.03509*cosph*(costh2 - 1.0)*(433.1*costh2*costh2 - 236.2*costh2 + 13.12);
1520  break;
1521  case 3:
1522  Y = -0.005849*cosph*pow(1.0 - 1.0*costh2, 1.5)*(-1732.0*costh2*costh + 472.5*costh);
1523  break;
1524  case 4:
1525  Y = 0.001068*cosph*(5198.0*costh2 - 472.5)*pow(costh2 - 1.0, 2);
1526  break;
1527  case 5:
1528  Y = 2.367*costh*cosph*pow(1.0 - 1.0*costh2, 2.5);
1529  break;
1530  case 6:
1531  Y = -0.6832*cosph*pow(costh2 - 1.0, 3);
1532  break;
1533  }break;
1534  case 7:
1535  switch (m)
1536  {
1537  case -7:
1538  Y = 0.7072*sinth*pow(1.0 - 1.0*costh2, 3.5);
1539  break;
1540  case -6:
1541  Y = -2.646*costh*sinth*pow(costh2 - 1.0, 3);
1542  break;
1543  case -5:
1544  Y = 9.984e-5*sinth*pow(1.0 - 1.0*costh2, 2.5)*(67570.0*costh2 - 5198.0);
1545  break;
1546  case -4:
1547  Y = -0.000599*sinth*pow(costh2 - 1.0, 2)*(-22520.0*costh2*costh + 5198.0*costh);
1548  break;
1549  case -3:
1550  Y = 0.003974*sinth*pow(1.0 - 1.0*costh2, 1.5)*(5631.0*costh2*costh2 - 2599.0*costh2 + 118.1);
1551  break;
1552  case -2:
1553  Y = -0.0281*sinth*(costh2 - 1.0)*(1126.0*costh2*costh2*costh - 866.2*costh2*costh + 118.1*costh);
1554  break;
1555  case -1:
1556  Y = 0.2065*sinth*pow(1.0 - 1.0*costh2, 0.5)*(187.7*costh2*costh2*costh2 - 216.6*costh2*costh2 + 59.06*costh2 - 2.188);
1557  break;
1558  case 0:
1559  Y = 29.29*costh2*costh2*costh2*costh - 47.32*costh2*costh2*costh + 21.51*costh2*costh - 2.39*costh;
1560  break;
1561  case 1:
1562  Y = 0.2065*cosph*pow(1.0 - 1.0*costh2, 0.5)*(187.7*costh2*costh2*costh2 - 216.6*costh2*costh2 + 59.06*costh2 - 2.188);
1563  break;
1564  case 2:
1565  Y = -0.0281*cosph*(costh2 - 1.0)*(1126.0*costh2*costh2*costh - 866.2*costh2*costh + 118.1*costh);
1566  break;
1567  case 3:
1568  Y = 0.003974*cosph*pow(1.0 - 1.0*costh2, 1.5)*(5631.0*costh2*costh2 - 2599.0*costh2 + 118.1);
1569  break;
1570  case 4:
1571  Y = -0.000599*cosph*pow(costh2 - 1.0, 2)*(- 22520.0*costh2*costh + 5198.0*costh);
1572  break;
1573  case 5:
1574  Y = 9.984e-5*cosph*pow(1.0 - 1.0*costh2, 2.5)*(67570.0*costh2 - 5198.0);
1575  break;
1576  case 6:
1577  Y = -2.646*cosph*costh*pow(costh2 - 1.0, 3);
1578  break;
1579  case 7:
1580  Y = 0.7072*cosph*pow(1.0 - 1.0*costh2, 3.5);
1581  break;
1582  }break;
1583  case 8:
1584  switch (m)
1585  {
1586  case -8:
1587  Y = sinth*pow(costh2-1.0,4.0)*7.289266601746931E-1;
1588  break;
1589  case -7:
1590  Y = costh*sinth*pow(costh2*-1.0+1.0,7.0/2.0)*2.915706640698772;
1591  break;
1592  case -6:
1593  Y = sinth*(costh2*1.0135125E+6-6.75675E+4)*pow(costh2-1.0,3.0)*-7.878532816224526E-6;
1594  break;
1595  case -5:
1596  Y = sinth*pow(costh2*-1.0+1.0,5.0/2.0)*(costh*6.75675E+4-(costh2*costh)*3.378375E+5)*-5.105872826582925E-56;
1597  break;
1598  case -4:
1599  Y = sinth*pow(costh2-1.0,2.0)*(costh2*-3.378375E+4+(costh2*costh2)*8.4459375E+4+1.299375E+3)*3.681897256448963E-4;
1600  break;
1601  case -3:
1602  Y = sinth*pow(costh2*-1.0+1.0,3.0/2.0)*(costh*1.299375E+3-(costh*costh2)*1.126125E+4+(costh*costh2*costh2)*1.6891875E+4)*2.851985351334463E-3;
1603  break;
1604  case -2:
1605  Y = sinth*(costh2-1.0)*(costh2*6.496875E+2-(costh2*costh2)*2.8153125E+3+(costh2*costh2*costh2)*2.8153125E+3-1.96875E+1)*-2.316963852365461E-2;
1606  break;
1607  case -1:
1608  Y = sinth*sqrt(costh2*-1.0+1.0)*(costh*1.96875E+1-(costh*costh2)*2.165625E+2+(costh*costh2*costh2)*5.630625E+2-(costh*costh2*costh2*costh2)*4.021875E+2)*-1.938511038201796E-1;;
1609  break;
1610  case 0:
1611  Y = costh2*-1.144933081936324E+1+(costh2*costh2)*6.297131950652692E+1-(costh2*costh2*costh2)*1.091502871445846E+2+(costh2*costh2*costh2*costh2)*5.847336811327841E+1+3.180369672045344E-1;
1612  break;
1613  case 1:
1614  Y = cosph*sqrt(costh2*-1.0+1.0)*(costh*1.96875E+1-(costh*costh2)*2.165625E+2+(costh*costh2*costh2)*5.630625E+2-(costh*costh2*costh2*costh2)*4.021875E+2)*-1.938511038201796E-1;;
1615  break;
1616  case 2:
1617  Y = cosph*(costh2-1.0)*(costh2*6.496875E+2-(costh2*costh2)*2.8153125E+3+(costh2*costh2*costh2)*2.8153125E+3-1.96875E+1)*-2.316963852365461E-2;
1618  break;
1619  case 3:
1620  Y = cosph*pow(costh2*-1.0+1.0,3.0/2.0)*(costh*1.299375E+3-(costh*costh2)*1.126125E+4+(costh*costh2*costh2)*1.6891875E+4)*2.851985351334463E-3;
1621  break;
1622  case 4:
1623  Y = cosph*pow(costh2-1.0,2.0)*(costh2*-3.378375E+4+(costh2*costh2)*8.4459375E+4+1.299375E+3)*3.681897256448963E-4;
1624  break;
1625  case 5:
1626  Y = cosph*pow(costh2*-1.0+1.0,5.0/2.0)*(costh*6.75675E+4-(costh*costh2)*3.378375E+5)*-5.105872826582925E-5;
1627  break;
1628  case 6:
1629  Y = cosph*(costh2*1.0135125E+6-6.75675E+4)*pow(costh2-1.0,3.0)*-7.878532816224526E-6;
1630  break;
1631  case 7:
1632  Y = cosph*costh*pow(costh2*-1.0+1.0,7.0/2.0)*2.915706640698772;
1633  break;
1634  case 8:
1635  Y = cosph*pow(costh2-1.0,4.0)*7.289266601746931E-1;
1636  break;
1637  }break;
1638  case 9:
1639  switch (m)
1640  {
1641  case -9:
1642  Y = sinth*pow(costh2*-1.0+1.0,9.0/2.0)*7.489009518540115E-1;
1643  break;
1644  case -8:
1645  Y = costh*sinth*pow(costh2-1.0,4.0)*3.17731764895143;
1646  break;
1647  case -7:
1648  Y = sinth*pow(costh2*-1.0+1.0,7.0/2.0)*(costh2*1.72297125E+7-1.0135125E+6)*5.376406125665728E-7;
1649  break;
1650  case -6:
1651  Y = sinth*(costh*1.0135125E+6-(costh*costh2)*5.7432375E+6)*pow(costh2-1.0,3.0)*3.724883428715686E-6;
1652  break;
1653  case -5:
1654  Y = sinth*pow(costh2*-1.0+1.0,5.0/2.0)*(costh2*-5.0675625E+5+(costh2*costh2)*1.435809375E+6+1.6891875E+4)*2.885282297193648E-5;
1655  break;
1656  case -4:
1657  Y = sinth*pow(costh2-1.0,2.0)*(costh*1.6891875E+4-(costh*costh2)*1.6891875E+5+(costh*costh2*costh2)*2.87161875E+5)*2.414000363328839E-4;
1658  break;
1659  case -3:
1660  Y = sinth*pow(costh2*-1.0+1.0,3.0/2.0)*((costh2)*8.4459375E+3-(costh2*costh2)*4.22296875E+4+(costh2*costh2*costh2)*4.78603125E+4-2.165625E+2)*2.131987394015766E-3;
1661  break;
1662  case -2:
1663  Y = sinth*(costh2-1.0)*(costh*2.165625E+2-(costh*costh2)*2.8153125E+3+(costh*costh2*costh2)*8.4459375E+3-(costh*costh2*costh2*costh2)*6.8371875E+3)*1.953998722751749E-2;
1664  break;
1665  case -1:
1666  Y = sinth*sqrt(costh2*-1.0+1.0)*(costh2*-1.0828125E+2+(costh2*costh2)*7.03828125E+2-(costh2*costh2*costh2)*1.40765625E+3+(costh2*costh2*costh2*costh2)*8.546484375E+2+2.4609375)*1.833013280775049E-1;
1667  break;
1668  case 0:
1669  Y = costh*3.026024588281871-(costh*costh2)*4.438169396144804E+1+(costh*costh2*costh2)*1.730886064497754E+2-(costh*costh2*costh2*costh2)*2.472694377852604E+2+(costh*costh2*costh2*costh2*costh2)*1.167661233986728E+2;
1670  break;
1671  case 1:
1672  Y = cosph*sqrt(costh2*-1.0+1.0)*(costh2*-1.0828125E+2+(costh2*costh2)*7.03828125E+2-(costh2*costh2*costh2)*1.40765625E+3+(costh2*costh2*costh2*costh2)*8.546484375E+2+2.4609375)*1.833013280775049E-1;
1673  break;
1674  case 2:
1675  Y = cosph*(costh2-1.0)*(costh*2.165625E+2-(costh*costh2)*2.8153125E+3+(costh*costh2*costh2)*8.4459375E+3-(costh*costh2*costh2*costh2)*6.8371875E+3)*1.953998722751749E-2;
1676  break;
1677  case 3:
1678  Y = cosph*pow(costh2*-1.0+1.0,3.0/2.0)*(costh2*8.4459375E+3-(costh2*costh2)*4.22296875E+4+(costh2*costh2*costh2)*4.78603125E+4-2.165625E+2)*2.131987394015766E-3;
1679  break;
1680  case 4:
1681  Y = cosph*pow(costh2-1.0,2.0)*(costh*1.6891875E+4-(costh*costh2)*1.6891875E+5+(costh*costh2*costh2)*2.87161875E+5)*2.414000363328839E-4;
1682  break;
1683  case 5:
1684  Y = cosph*pow(costh2*-1.0+1.0,5.0/2.0)*(costh2*-5.0675625E+5+(costh2*costh2)*1.435809375E+6+1.6891875E+4)*2.885282297193648E-5;
1685  break;
1686  case 6:
1687  Y = cosph*(costh*1.0135125E+6-(costh*costh2)*5.7432375E+6)*pow(costh2-1.0,3.0)*3.724883428715686E-6;
1688  break;
1689  case 7:
1690  Y = cosph*pow(costh2*-1.0+1.0,7.0/2.0)*(costh2*1.72297125E+7-1.0135125E+6)*5.376406125665728E-7;
1691  break;
1692  case 8:
1693  Y = cosph*costh*pow(costh2-1.0,4.0)*3.17731764895143;
1694  break;
1695  case 9:
1696  Y = cosph*pow(costh2*-1.0+1.0,9.0/2.0)*7.489009518540115E-1;
1697  break;
1698  }break;
1699  default: break;
1700  }
1701 
1702  return R*Y;
1703 }
L2
#define L2
Definition: cuda_volume_deform_sph_defines.h:17
sqrt
void sqrt(Image< double > &op)
Definition: image_operate.cpp:210
abs
void abs(Image< double > &op)
Definition: image_operate.cpp:219
m
int m
Definition: image_find_center.cpp:57
r2
float r2
Definition: image_find_center.cpp:54
PI
#define PI
Definition: tools.h:43
n
int * n
Definition: numerical_recipes.cpp:2229
L1
#define L1
Definition: cuda_volume_deform_sph_defines.h:14

◆ ZernikeSphericalHarmonics() [2/2]

double ZernikeSphericalHarmonics ( int  l1,
int  n,
int  l2,
int  m,
double  xr,
double  yr,
double  zr,
double  r 
)

Spherical harmonics. This function calculates R_l^n(r)Y_l^m(xr,yr,zr) where R_l^n is the Zernike polynomial, l is the degree and n goes over odd or even values depending on the value of l. For instance, l=0 (n=0), l=1 (n=1), l=2 (n=0,2), l=3 (n=1,3), l=4 (n=0,2,4), ... The Y_l^m is the real spherical harmonic. l is the degree, and m=-l,...,l. For the specific formulas see https://en.wikipedia.org/wiki/Zernike_polynomials#Radial_polynomials and https://en.wikipedia.org/wiki/Table_of_spherical_harmonics#Real_spherical_harmonics. The cartesian coordinates xr, yr and zr are supposed to be normalized between -1 and 1, that is, xr=x/r, yr=y/r, and zr=z/r. r is supposed to be between 0 and 1.

Definition at line 1705 of file numerical_tools.cpp.

1705  {
1706  switch (l1) {
1707  case 0 : {
1708  switch (l2) {
1709  case 0: return ZernikeSphericalHarmonics<0, 0>(n, m, xr, yr, zr, r);
1710  case 1: return ZernikeSphericalHarmonics<0, 1>(n, m, xr, yr, zr, r);
1711  case 2: return ZernikeSphericalHarmonics<0, 2>(n, m, xr, yr, zr, r);
1712  case 3: return ZernikeSphericalHarmonics<0, 3>(n, m, xr, yr, zr, r);
1713  case 4: return ZernikeSphericalHarmonics<0, 4>(n, m, xr, yr, zr, r);
1714  case 5: return ZernikeSphericalHarmonics<0, 5>(n, m, xr, yr, zr, r);
1715  case 6: return ZernikeSphericalHarmonics<0, 6>(n, m, xr, yr, zr, r);
1716  case 7: return ZernikeSphericalHarmonics<0, 7>(n, m, xr, yr, zr, r);
1717  case 8: return ZernikeSphericalHarmonics<0, 8>(n, m, xr, yr, zr, r);
1718  case 9: return ZernikeSphericalHarmonics<0, 9>(n, m, xr, yr, zr, r);
1719  default: break;
1720  }
1721  } break;
1722  case 1 : {
1723  switch (l2) {
1724  case 0: return ZernikeSphericalHarmonics<1, 0>(n, m, xr, yr, zr, r);
1725  case 1: return ZernikeSphericalHarmonics<1, 1>(n, m, xr, yr, zr, r);
1726  case 2: return ZernikeSphericalHarmonics<1, 2>(n, m, xr, yr, zr, r);
1727  case 3: return ZernikeSphericalHarmonics<1, 3>(n, m, xr, yr, zr, r);
1728  case 4: return ZernikeSphericalHarmonics<1, 4>(n, m, xr, yr, zr, r);
1729  case 5: return ZernikeSphericalHarmonics<1, 5>(n, m, xr, yr, zr, r);
1730  case 6: return ZernikeSphericalHarmonics<1, 6>(n, m, xr, yr, zr, r);
1731  case 7: return ZernikeSphericalHarmonics<1, 7>(n, m, xr, yr, zr, r);
1732  case 8: return ZernikeSphericalHarmonics<1, 8>(n, m, xr, yr, zr, r);
1733  case 9: return ZernikeSphericalHarmonics<1, 9>(n, m, xr, yr, zr, r);
1734  default: break;
1735  }
1736  } break;
1737  case 2 : {
1738  switch (l2) {
1739  case 0: return ZernikeSphericalHarmonics<2, 0>(n, m, xr, yr, zr, r);
1740  case 1: return ZernikeSphericalHarmonics<2, 1>(n, m, xr, yr, zr, r);
1741  case 2: return ZernikeSphericalHarmonics<2, 2>(n, m, xr, yr, zr, r);
1742  case 3: return ZernikeSphericalHarmonics<2, 3>(n, m, xr, yr, zr, r);
1743  case 4: return ZernikeSphericalHarmonics<2, 4>(n, m, xr, yr, zr, r);
1744  case 5: return ZernikeSphericalHarmonics<2, 5>(n, m, xr, yr, zr, r);
1745  case 6: return ZernikeSphericalHarmonics<2, 6>(n, m, xr, yr, zr, r);
1746  case 7: return ZernikeSphericalHarmonics<2, 7>(n, m, xr, yr, zr, r);
1747  case 8: return ZernikeSphericalHarmonics<2, 8>(n, m, xr, yr, zr, r);
1748  case 9: return ZernikeSphericalHarmonics<2, 9>(n, m, xr, yr, zr, r);
1749  default: break;
1750  }
1751  } break;
1752  case 3 : {
1753  switch (l2) {
1754  case 0: return ZernikeSphericalHarmonics<3, 0>(n, m, xr, yr, zr, r);
1755  case 1: return ZernikeSphericalHarmonics<3, 1>(n, m, xr, yr, zr, r);
1756  case 2: return ZernikeSphericalHarmonics<3, 2>(n, m, xr, yr, zr, r);
1757  case 3: return ZernikeSphericalHarmonics<3, 3>(n, m, xr, yr, zr, r);
1758  case 4: return ZernikeSphericalHarmonics<3, 4>(n, m, xr, yr, zr, r);
1759  case 5: return ZernikeSphericalHarmonics<3, 5>(n, m, xr, yr, zr, r);
1760  case 6: return ZernikeSphericalHarmonics<3, 6>(n, m, xr, yr, zr, r);
1761  case 7: return ZernikeSphericalHarmonics<3, 7>(n, m, xr, yr, zr, r);
1762  case 8: return ZernikeSphericalHarmonics<3, 8>(n, m, xr, yr, zr, r);
1763  case 9: return ZernikeSphericalHarmonics<3, 9>(n, m, xr, yr, zr, r);
1764  default: break;
1765  }
1766  } break;
1767  case 4 : {
1768  switch (l2) {
1769  case 0: return ZernikeSphericalHarmonics<4, 0>(n, m, xr, yr, zr, r);
1770  case 1: return ZernikeSphericalHarmonics<4, 1>(n, m, xr, yr, zr, r);
1771  case 2: return ZernikeSphericalHarmonics<4, 2>(n, m, xr, yr, zr, r);
1772  case 3: return ZernikeSphericalHarmonics<4, 3>(n, m, xr, yr, zr, r);
1773  case 4: return ZernikeSphericalHarmonics<4, 4>(n, m, xr, yr, zr, r);
1774  case 5: return ZernikeSphericalHarmonics<4, 5>(n, m, xr, yr, zr, r);
1775  case 6: return ZernikeSphericalHarmonics<4, 6>(n, m, xr, yr, zr, r);
1776  case 7: return ZernikeSphericalHarmonics<4, 7>(n, m, xr, yr, zr, r);
1777  case 8: return ZernikeSphericalHarmonics<4, 8>(n, m, xr, yr, zr, r);
1778  case 9: return ZernikeSphericalHarmonics<4, 9>(n, m, xr, yr, zr, r);
1779  default: break;
1780  }
1781  } break;
1782  case 5 : {
1783  switch (l2) {
1784  case 0: return ZernikeSphericalHarmonics<5, 0>(n, m, xr, yr, zr, r);
1785  case 1: return ZernikeSphericalHarmonics<5, 1>(n, m, xr, yr, zr, r);
1786  case 2: return ZernikeSphericalHarmonics<5, 2>(n, m, xr, yr, zr, r);
1787  case 3: return ZernikeSphericalHarmonics<5, 3>(n, m, xr, yr, zr, r);
1788  case 4: return ZernikeSphericalHarmonics<5, 4>(n, m, xr, yr, zr, r);
1789  case 5: return ZernikeSphericalHarmonics<5, 5>(n, m, xr, yr, zr, r);
1790  case 6: return ZernikeSphericalHarmonics<5, 6>(n, m, xr, yr, zr, r);
1791  case 7: return ZernikeSphericalHarmonics<5, 7>(n, m, xr, yr, zr, r);
1792  case 8: return ZernikeSphericalHarmonics<5, 8>(n, m, xr, yr, zr, r);
1793  case 9: return ZernikeSphericalHarmonics<5, 9>(n, m, xr, yr, zr, r);
1794  default: break;
1795  }
1796  } break;
1797  case 6 : {
1798  switch (l2) {
1799  case 0: return ZernikeSphericalHarmonics<6, 0>(n, m, xr, yr, zr, r);
1800  case 1: return ZernikeSphericalHarmonics<6, 1>(n, m, xr, yr, zr, r);
1801  case 2: return ZernikeSphericalHarmonics<6, 2>(n, m, xr, yr, zr, r);
1802  case 3: return ZernikeSphericalHarmonics<6, 3>(n, m, xr, yr, zr, r);
1803  case 4: return ZernikeSphericalHarmonics<6, 4>(n, m, xr, yr, zr, r);
1804  case 5: return ZernikeSphericalHarmonics<6, 5>(n, m, xr, yr, zr, r);
1805  case 6: return ZernikeSphericalHarmonics<6, 6>(n, m, xr, yr, zr, r);
1806  case 7: return ZernikeSphericalHarmonics<6, 7>(n, m, xr, yr, zr, r);
1807  case 8: return ZernikeSphericalHarmonics<6, 8>(n, m, xr, yr, zr, r);
1808  case 9: return ZernikeSphericalHarmonics<6, 9>(n, m, xr, yr, zr, r);
1809  default: break;
1810  }
1811  } break;
1812  case 7 : {
1813  switch (l2) {
1814  case 0: return ZernikeSphericalHarmonics<7, 0>(n, m, xr, yr, zr, r);
1815  case 1: return ZernikeSphericalHarmonics<7, 1>(n, m, xr, yr, zr, r);
1816  case 2: return ZernikeSphericalHarmonics<7, 2>(n, m, xr, yr, zr, r);
1817  case 3: return ZernikeSphericalHarmonics<7, 3>(n, m, xr, yr, zr, r);
1818  case 4: return ZernikeSphericalHarmonics<7, 4>(n, m, xr, yr, zr, r);
1819  case 5: return ZernikeSphericalHarmonics<7, 5>(n, m, xr, yr, zr, r);
1820  case 6: return ZernikeSphericalHarmonics<7, 6>(n, m, xr, yr, zr, r);
1821  case 7: return ZernikeSphericalHarmonics<7, 7>(n, m, xr, yr, zr, r);
1822  case 8: return ZernikeSphericalHarmonics<7, 8>(n, m, xr, yr, zr, r);
1823  case 9: return ZernikeSphericalHarmonics<7, 9>(n, m, xr, yr, zr, r);
1824  default: break;
1825  }
1826  } break;
1827  case 8 : {
1828  switch (l2) {
1829  case 0: return ZernikeSphericalHarmonics<8, 0>(n, m, xr, yr, zr, r);
1830  case 1: return ZernikeSphericalHarmonics<8, 1>(n, m, xr, yr, zr, r);
1831  case 2: return ZernikeSphericalHarmonics<8, 2>(n, m, xr, yr, zr, r);
1832  case 3: return ZernikeSphericalHarmonics<8, 3>(n, m, xr, yr, zr, r);
1833  case 4: return ZernikeSphericalHarmonics<8, 4>(n, m, xr, yr, zr, r);
1834  case 5: return ZernikeSphericalHarmonics<8, 5>(n, m, xr, yr, zr, r);
1835  case 6: return ZernikeSphericalHarmonics<8, 6>(n, m, xr, yr, zr, r);
1836  case 7: return ZernikeSphericalHarmonics<8, 7>(n, m, xr, yr, zr, r);
1837  case 8: return ZernikeSphericalHarmonics<8, 8>(n, m, xr, yr, zr, r);
1838  case 9: return ZernikeSphericalHarmonics<8, 9>(n, m, xr, yr, zr, r);
1839  default: break;
1840  }
1841  } break;
1842  case 9 : {
1843  switch (l2) {
1844  case 0: return ZernikeSphericalHarmonics<9, 0>(n, m, xr, yr, zr, r);
1845  case 1: return ZernikeSphericalHarmonics<9, 1>(n, m, xr, yr, zr, r);
1846  case 2: return ZernikeSphericalHarmonics<9, 2>(n, m, xr, yr, zr, r);
1847  case 3: return ZernikeSphericalHarmonics<9, 3>(n, m, xr, yr, zr, r);
1848  case 4: return ZernikeSphericalHarmonics<9, 4>(n, m, xr, yr, zr, r);
1849  case 5: return ZernikeSphericalHarmonics<9, 5>(n, m, xr, yr, zr, r);
1850  case 6: return ZernikeSphericalHarmonics<9, 6>(n, m, xr, yr, zr, r);
1851  case 7: return ZernikeSphericalHarmonics<9, 7>(n, m, xr, yr, zr, r);
1852  case 8: return ZernikeSphericalHarmonics<9, 8>(n, m, xr, yr, zr, r);
1853  case 9: return ZernikeSphericalHarmonics<9, 9>(n, m, xr, yr, zr, r);
1854  default: break;
1855  }
1856  } break;
1857  case 10 : {
1858  switch (l2) {
1859  case 0: return ZernikeSphericalHarmonics<10, 0>(n, m, xr, yr, zr, r);
1860  case 1: return ZernikeSphericalHarmonics<10, 1>(n, m, xr, yr, zr, r);
1861  case 2: return ZernikeSphericalHarmonics<10, 2>(n, m, xr, yr, zr, r);
1862  case 3: return ZernikeSphericalHarmonics<10, 3>(n, m, xr, yr, zr, r);
1863  case 4: return ZernikeSphericalHarmonics<10, 4>(n, m, xr, yr, zr, r);
1864  case 5: return ZernikeSphericalHarmonics<10, 5>(n, m, xr, yr, zr, r);
1865  case 6: return ZernikeSphericalHarmonics<10, 6>(n, m, xr, yr, zr, r);
1866  case 7: return ZernikeSphericalHarmonics<10, 7>(n, m, xr, yr, zr, r);
1867  case 8: return ZernikeSphericalHarmonics<10, 8>(n, m, xr, yr, zr, r);
1868  case 9: return ZernikeSphericalHarmonics<10, 9>(n, m, xr, yr, zr, r);
1869  default: break;
1870  }
1871  } break;
1872  case 11 : {
1873  switch (l2) {
1874  case 0: return ZernikeSphericalHarmonics<11, 0>(n, m, xr, yr, zr, r);
1875  case 1: return ZernikeSphericalHarmonics<11, 1>(n, m, xr, yr, zr, r);
1876  case 2: return ZernikeSphericalHarmonics<11, 2>(n, m, xr, yr, zr, r);
1877  case 3: return ZernikeSphericalHarmonics<11, 3>(n, m, xr, yr, zr, r);
1878  case 4: return ZernikeSphericalHarmonics<11, 4>(n, m, xr, yr, zr, r);
1879  case 5: return ZernikeSphericalHarmonics<11, 5>(n, m, xr, yr, zr, r);
1880  case 6: return ZernikeSphericalHarmonics<11, 6>(n, m, xr, yr, zr, r);
1881  case 7: return ZernikeSphericalHarmonics<11, 7>(n, m, xr, yr, zr, r);
1882  case 8: return ZernikeSphericalHarmonics<11, 8>(n, m, xr, yr, zr, r);
1883  case 9: return ZernikeSphericalHarmonics<11, 9>(n, m, xr, yr, zr, r);
1884  default: break;
1885  }
1886  } break;
1887  case 12 : {
1888  switch (l2) {
1889  case 0: return ZernikeSphericalHarmonics<12, 0>(n, m, xr, yr, zr, r);
1890  case 1: return ZernikeSphericalHarmonics<12, 1>(n, m, xr, yr, zr, r);
1891  case 2: return ZernikeSphericalHarmonics<12, 2>(n, m, xr, yr, zr, r);
1892  case 3: return ZernikeSphericalHarmonics<12, 3>(n, m, xr, yr, zr, r);
1893  case 4: return ZernikeSphericalHarmonics<12, 4>(n, m, xr, yr, zr, r);
1894  case 5: return ZernikeSphericalHarmonics<12, 5>(n, m, xr, yr, zr, r);
1895  case 6: return ZernikeSphericalHarmonics<12, 6>(n, m, xr, yr, zr, r);
1896  case 7: return ZernikeSphericalHarmonics<12, 7>(n, m, xr, yr, zr, r);
1897  case 8: return ZernikeSphericalHarmonics<12, 8>(n, m, xr, yr, zr, r);
1898  case 9: return ZernikeSphericalHarmonics<12, 9>(n, m, xr, yr, zr, r);
1899  default: break;
1900  }
1901  } break;
1902  case 13 : {
1903  switch (l2) {
1904  case 0: return ZernikeSphericalHarmonics<13, 0>(n, m, xr, yr, zr, r);
1905  case 1: return ZernikeSphericalHarmonics<13, 1>(n, m, xr, yr, zr, r);
1906  case 2: return ZernikeSphericalHarmonics<13, 2>(n, m, xr, yr, zr, r);
1907  case 3: return ZernikeSphericalHarmonics<13, 3>(n, m, xr, yr, zr, r);
1908  case 4: return ZernikeSphericalHarmonics<13, 4>(n, m, xr, yr, zr, r);
1909  case 5: return ZernikeSphericalHarmonics<13, 5>(n, m, xr, yr, zr, r);
1910  case 6: return ZernikeSphericalHarmonics<13, 6>(n, m, xr, yr, zr, r);
1911  case 7: return ZernikeSphericalHarmonics<13, 7>(n, m, xr, yr, zr, r);
1912  case 8: return ZernikeSphericalHarmonics<13, 8>(n, m, xr, yr, zr, r);
1913  case 9: return ZernikeSphericalHarmonics<13, 9>(n, m, xr, yr, zr, r);
1914  default: break;
1915  }
1916  } break;
1917  case 14 : {
1918  switch (l2) {
1919  case 0: return ZernikeSphericalHarmonics<14, 0>(n, m, xr, yr, zr, r);
1920  case 1: return ZernikeSphericalHarmonics<14, 1>(n, m, xr, yr, zr, r);
1921  case 2: return ZernikeSphericalHarmonics<14, 2>(n, m, xr, yr, zr, r);
1922  case 3: return ZernikeSphericalHarmonics<14, 3>(n, m, xr, yr, zr, r);
1923  case 4: return ZernikeSphericalHarmonics<14, 4>(n, m, xr, yr, zr, r);
1924  case 5: return ZernikeSphericalHarmonics<14, 5>(n, m, xr, yr, zr, r);
1925  case 6: return ZernikeSphericalHarmonics<14, 6>(n, m, xr, yr, zr, r);
1926  case 7: return ZernikeSphericalHarmonics<14, 7>(n, m, xr, yr, zr, r);
1927  case 8: return ZernikeSphericalHarmonics<14, 8>(n, m, xr, yr, zr, r);
1928  case 9: return ZernikeSphericalHarmonics<14, 9>(n, m, xr, yr, zr, r);
1929  default: break;
1930  }
1931  } break;
1932  case 15 : {
1933  switch (l2) {
1934  case 0: return ZernikeSphericalHarmonics<15, 0>(n, m, xr, yr, zr, r);
1935  case 1: return ZernikeSphericalHarmonics<15, 1>(n, m, xr, yr, zr, r);
1936  case 2: return ZernikeSphericalHarmonics<15, 2>(n, m, xr, yr, zr, r);
1937  case 3: return ZernikeSphericalHarmonics<15, 3>(n, m, xr, yr, zr, r);
1938  case 4: return ZernikeSphericalHarmonics<15, 4>(n, m, xr, yr, zr, r);
1939  case 5: return ZernikeSphericalHarmonics<15, 5>(n, m, xr, yr, zr, r);
1940  case 6: return ZernikeSphericalHarmonics<15, 6>(n, m, xr, yr, zr, r);
1941  case 7: return ZernikeSphericalHarmonics<15, 7>(n, m, xr, yr, zr, r);
1942  case 8: return ZernikeSphericalHarmonics<15, 8>(n, m, xr, yr, zr, r);
1943  case 9: return ZernikeSphericalHarmonics<15, 9>(n, m, xr, yr, zr, r);
1944  default: break;
1945  }
1946  } break;
1947  default: break;
1948  }
1949  REPORT_ERROR(ERR_ARG_INCORRECT, "ZernikeSphericalHarmonics not supported for l1 = " + std::to_string(l1) + " and l2 = " + std::to_string(l2));
1950 }
REPORT_ERROR
#define REPORT_ERROR(nerr, ErrormMsg)
Definition: xmipp_error.h:211
ERR_ARG_INCORRECT
Incorrect argument received.
Definition: xmipp_error.h:113
m
int m
Definition: image_find_center.cpp:57
cif::to_string
std::string to_string(bond_type bondType)
Definition: compound.cpp:43
n
int * n
Definition: numerical_recipes.cpp:2229
Generated by   doxygen 1.8.13