diff --git a/CMakeLists.txt b/CMakeLists.txt index 887dc04..2427097 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,26 +1,28 @@ # Name of project project(wignerSymbols) set (wignerSymbols_VERSION_MAJOR 0) -set (wignerSymbols_VERSION_MINOR 1) -set (wignerSymbols_VERSION_RELEASE 3) +set (wignerSymbols_VERSION_MINOR 2) +set (wignerSymbols_VERSION_RELEASE 0) # CMake config cmake_minimum_required(VERSION 2.8) -set (CMAKE_INSTALL_PREFIX /usr/) +if (CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT) + set (CMAKE_INSTALL_PREFIX /usr) +endif() # Compiler config enable_language (Fortran) -set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -Wall -march=native") +set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -Wall -march=native") # Included files -include_directories(${CMAKE_CURRENT_SOURCE_DIR}/include) +include_directories(${CMAKE_CURRENT_SOURCE_DIR}/include INC_LIST) # Source files aux_source_directory(./src SRC_LIST) # Build a shared library -add_library(${PROJECT_NAME} SHARED - ${SRC_LIST} +add_library(${PROJECT_NAME} SHARED + ${SRC_LIST} "./src/wignerSymbols-fortran.f" "./src/machine.for" "./src/fdump.f" diff --git a/include/wignerSymbols.h b/include/wignerSymbols.h index 75e4b0a..ee23cae 100644 --- a/include/wignerSymbols.h +++ b/include/wignerSymbols.h @@ -3,7 +3,7 @@ * Lesser Public License. If a copy of the LGPL was not -/ * distributed with this file, you can obtain one at -/ * https://www.gnu.org/licenses/lgpl.html. -/ - ********************************************************/ + ********************************************************/ #ifndef WIGNER_SYMBOLS_H #define WIGNER_SYMBOLS_H diff --git a/include/wignerSymbols/commonFunctions.h b/include/wignerSymbols/commonFunctions.h index 565e70d..c8ace8e 100644 --- a/include/wignerSymbols/commonFunctions.h +++ b/include/wignerSymbols/commonFunctions.h @@ -3,24 +3,25 @@ * Lesser Public License. If a copy of the LGPL was not -/ * distributed with this file, you can obtain one at -/ * https://www.gnu.org/licenses/lgpl.html. -/ - ********************************************************/ + ******************************************************** -/** \file commonFunctions.h +/** \file commonFunctions.h * * \author Joey Dumont * * \since 2014-07-10 * - * \brief Defines some common functions to the C++ and Fortran code. + * \brief Defines some common functions to the C++ and Fortran code. * * We define some functions that will be of use in both the C++ and Fortran - * parts of this library. + * parts of this library. * */ namespace WignerSymbols { -template -double sgn(T val) + +template +double sgn(T val) { int sgn = (T(0) < val) - (val < T(0)); if (sgn == 0) @@ -28,4 +29,5 @@ double sgn(T val) else return (double)sgn; } -} + +} // namespace WignerSymbols diff --git a/include/wignerSymbols/wignerSymbols-cpp.h b/include/wignerSymbols/wignerSymbols-cpp.h index 8811958..b0fe693 100644 --- a/include/wignerSymbols/wignerSymbols-cpp.h +++ b/include/wignerSymbols/wignerSymbols-cpp.h @@ -3,28 +3,28 @@ * Lesser Public License. If a copy of the LGPL was not -/ * distributed with this file, you can obtain one at -/ * https://www.gnu.org/licenses/lgpl.html. -/ - ********************************************************/ + ********************************************************/ #ifndef WIGNER_SYMBOLS_CPP_H #define WIGNER_SYMBOLS_CPP_H -/** \file wignerSymbols-cpp.h +/** \file wignerSymbols-cpp.h * * \author Joey Dumont * * \since 2013-08-16 * - * \brief Defines utility functions for the evaluation of Wigner-3j and -6j symbols. + * \brief Defines utility functions for the evaluation of Wigner-3j and -6j symbols. * - * We compute the Wigner-3j and -6j symbols + * We compute the Wigner-3j and -6j symbols * f(L1) = ( L1 L2 L3) * (-M2-M3 M2 M3) * for all allowed values of L1, the other parameters - * being held fixed. The algorithm is based on the work + * being held fixed. The algorithm is based on the work * by Schulten and Gordon. * K. Schulten, "Exact recursive evaluation of 3j- and 6j-coefficients for quantum-mechanical coupling of angular momenta," * J. Math. Phys. 16, 1961 (1975). - * K. Schulten and R. G. Gordon, "Recursive evaluation of 3j and 6j coefficients," + * K. Schulten and R. G. Gordon, "Recursive evaluation of 3j and 6j coefficients," * Comput. Phys. Commun. 11, 269–278 (1976). */ @@ -36,8 +36,8 @@ namespace WignerSymbols { -/*! @name Evaluation of Wigner-3j and -6j symbols. - * We implement Schulten's algorithm in C++. +/*! @name Evaluation of Wigner-3j and -6j symbols. + * We implement Schulten's algorithm in C++. */ ///@{ std::vector wigner3j(double l2, double l3, @@ -64,12 +64,12 @@ double wigner6j_auxB(double l1, double l2, double l3, double l4, double l5, double l6); /*! Computes the Clebsch-Gordan coefficient by relating it to the - * Wigner 3j symbol. It sometimes eases the notation to use the + * Wigner 3j symbol. It sometimes eases the notation to use the * Clebsch-Gordan coefficients directly. */ inline double clebschGordan(double l1, double l2, double l3, double m1, double m2, double m3) { - // We simply compute it via the 3j symbol. + // We simply compute it via the 3j symbol. return (pow(-1.0,l1-l2+m3)*sqrt(2.0*l3+1.0)*wigner3j(l1,l2,l3,m1,m2,-m3)); } } diff --git a/include/wignerSymbols/wignerSymbols-fortran.h b/include/wignerSymbols/wignerSymbols-fortran.h index d25ca3f..695bee3 100644 --- a/include/wignerSymbols/wignerSymbols-fortran.h +++ b/include/wignerSymbols/wignerSymbols-fortran.h @@ -3,30 +3,30 @@ * Lesser Public License. If a copy of the LGPL was not -/ * distributed with this file, you can obtain one at -/ * https://www.gnu.org/licenses/lgpl.html. -/ - ********************************************************/ + ********************************************************/ #ifndef WIGNER_SYMBOLS_FORTRAN_H #define WIGNER_SYMBOLS_FORTRAN_H -/** \file wignerSymbols-fortran.h +/** \file wignerSymbols-fortran.h * * \author Joey Dumont * * \since 2013-08-16 * - * \brief Defines utility functions for the evaluation of Wigner-3j and -6j symbols. + * \brief Defines utility functions for the evaluation of Wigner-3j and -6j symbols. * - * We use modified SLATEC (http://netlib.org/slatec) Fortran subroutines to compute - * Wigner-3j and -6j symbols. We modified the subroutines so that they do not depend - * d1mach, r1mach or i1mach as these are obsolete routines. They have been replaced + * We use modified SLATEC (http://netlib.org/slatec) Fortran subroutines to compute + * Wigner-3j and -6j symbols. We modified the subroutines so that they do not depend + * d1mach, r1mach or i1mach as these are obsolete routines. They have been replaced * by intrinsics such as huge(), tiny(), epsilon() and spacing(), which are guaranteed * to work. The file wignerSymbols.f contain the subroutines and their - * dependencies. + * dependencies. * * We rely on the ISO C Binding to bind the Fortran subroutines to C++. This method * of working insures that proper type casts are performed, among other things. We * then use the same method as we did before (extern "C"). - * + * */ #include @@ -50,11 +50,11 @@ std::vector wigner3j_f(double l2, double l3, double m1, double m2, doubl double wigner3j_f(double l1, double l2, double l3, double m1, double m2, double m3); /*! Computes the Clebsch-Gordan coefficient by relating it to the - * Wigner 3j symbol. It sometimes eases the notation to use the + * Wigner 3j symbol. It sometimes eases the notation to use the * Clebsch-Gordan coefficients directly. */ inline double clebschGordan_f(double l1, double l2, double l3, double m1, double m2, double m3) { - // We simply compute it via the 3j symbol. + // We simply compute it via the 3j symbol. return (pow(-1.0,l1-l2+m3)*sqrt(2.0*l3+1.0)*wigner3j_f(l1,l2,l3,m1,m2,-m3)); } diff --git a/src/wignerSymbols-cpp.cpp b/src/wignerSymbols-cpp.cpp index 5bfbd8f..db5154e 100644 --- a/src/wignerSymbols-cpp.cpp +++ b/src/wignerSymbols-cpp.cpp @@ -3,7 +3,7 @@ * Lesser Public License. If a copy of the LGPL was not -/ * distributed with this file, you can obtain one at -/ * https://www.gnu.org/licenses/lgpl.html. -/ - ********************************************************/ + ********************************************************/ #include "../include/wignerSymbols/commonFunctions.h" #include "../include/wignerSymbols/wignerSymbols-cpp.h" @@ -12,7 +12,7 @@ namespace WignerSymbols { std::vector wigner3j(double l2, double l3, double m1, double m2, double m3) { - // We compute the numeric limits of double precision. + // We compute the numeric limits of double precision. double huge = sqrt(std::numeric_limits::max()/20.0); double srhuge = sqrt(huge); double tiny = std::numeric_limits::min(); @@ -29,7 +29,7 @@ std::vector wigner3j(double l2, double l3, if (!select) return std::vector(1,0.0); - // We compute the limits of l1. + // We compute the limits of l1. double l1min = std::max(std::fabs(l2-l3),std::fabs(m1)); double l1max = l2+l3; @@ -46,10 +46,10 @@ std::vector wigner3j(double l2, double l3, // Another special case where the recursion relation fails. else { - // We start with an arbitrary value. + // We start with an arbitrary value. thrcof[0] = srtiny; - // From now on, we check the variation of |alpha(l1)|. + // From now on, we check the variation of |alpha(l1)|. double alphaNew, l1(l1min); if (l1min==0.0) alphaNew = -(m3-m2+2.0*wigner3j_auxB(l1,l2,l3,m1,m2,m3))/wigner3j_auxA(1.0,l2,l3,m1,m2,m3); @@ -63,20 +63,20 @@ std::vector wigner3j(double l2, double l3, // If size > 2, we start the forward recursion. if (size>2) { - // We start with an arbitrary value. + // We start with an arbitrary value. thrcof[0] = srtiny; - - // From now on, we check the variation of |alpha(l1)|. + + // From now on, we check the variation of |alpha(l1)|. double alphaOld, alphaNew, beta, l1(l1min); if (l1min==0.0) alphaNew = -(m3-m2+2.0*wigner3j_auxB(l1,l2,l3,m1,m2,m3))/wigner3j_auxA(1.0,l2,l3,m1,m2,m3); else alphaNew = -wigner3j_auxB(l1min,l2,l3,m1,m2,m3) /(l1min*wigner3j_auxA(l1min+1.0,l2,l3,m1,m2,m3)); - + // We compute the two-term recursion. thrcof[1] = alphaNew*thrcof[0]; - + // We compute the rest of the recursion. int i = 1; bool alphaVar = false; @@ -86,17 +86,17 @@ std::vector wigner3j(double l2, double l3, i++; // Next term in recursion alphaOld = alphaNew; // Monitoring of |alpha(l1)|. l1 += 1.0; // l1 = l1+1 - - // New coefficients in recursion. + + // New coefficients in recursion. alphaNew = -wigner3j_auxB(l1,l2,l3,m1,m2,m3) /(l1*wigner3j_auxA(l1+1.0,l2,l3,m1,m2,m3)); - + beta = -(l1+1.0)*wigner3j_auxA(l1,l2,l3,m1,m2,m3) /(l1*wigner3j_auxA(l1+1.0,l2,l3,m1,m2,m3)); - - // Application of the recursion. + + // Application of the recursion. thrcof[i] = alphaNew*thrcof[i-1]+beta*thrcof[i-2]; - + // We check if we are overflowing. if (std::fabs(thrcof[i])>srhuge) { @@ -104,57 +104,57 @@ std::vector wigner3j(double l2, double l3, for (std::vector::iterator it = thrcof.begin(); it != thrcof.begin()+i; ++it) { //if (std::fabs(*it) < srtiny) *it = 0; - //else + //else *it /= srhuge; } } - + // This piece of code checks whether we have reached - // the classical region. If we have, the second if + // the classical region. If we have, the second if // sets alphaVar to true and we break this loop at the - // next iteration because we need thrcof(l1mid+1) to + // next iteration because we need thrcof(l1mid+1) to // compute the scalar lambda. if (alphaVar) break; if (std::fabs(alphaNew)-std::fabs(alphaOld)>0.0) alphaVar=true; - + } while (i<(size-1)); // Loop stops when we have computed all values. - + // If this is the case, we have stumbled upon a classical region. // We start the backwards recursion. if (i!=size-1) { // We keep the two terms around l1mid to compute the factor later. double l1midm1(thrcof[i-2]),l1mid(thrcof[i-1]),l1midp1(thrcof[i]); - - // We compute the backward recursion by providing an arbitrary + + // We compute the backward recursion by providing an arbitrary // startint value. thrcof[size-1] = srtiny; - + // We compute the two-term recursion. l1 = l1max; alphaNew = -wigner3j_auxB(l1,l2,l3,m1,m2,m3) /((l1+1.0)*wigner3j_auxA(l1,l2,l3,m1,m2,m3)); thrcof[size-2] = alphaNew*thrcof[size-1]; - + // We compute the rest of the backward recursion. int j = size-2; do { // Bookkeeping - j--; // Previous term in recursion. + j--; // Previous term in recursion. l1 -= 1.0; // l1 = l1-1 - - // New coefficients in recursion. + + // New coefficients in recursion. alphaNew = -wigner3j_auxB(l1,l2,l3,m1,m2,m3) /((l1+1.0)*wigner3j_auxA(l1,l2,l3,m1,m2,m3)); beta = -l1*wigner3j_auxA(l1+1.0,l2,l3,m1,m2,m3) /((l1+1.0)*wigner3j_auxA(l1,l2,l3,m1,m2,m3)); - - // Application of the recursion. + + // Application of the recursion. thrcof[j] = alphaNew*thrcof[j+1]+beta*thrcof[j+2]; - + // We check if we are overflowing. if (std::fabs(thrcof[j]>srhuge)) { @@ -162,17 +162,17 @@ std::vector wigner3j(double l2, double l3, for (std::vector::iterator it = thrcof.begin()+j; it != thrcof.end(); ++it) { //if (std::fabs(*it) < srtiny) *it = 0; - //else + //else *it /= srhuge; } } - + } while (j>(i-2)); // Loop stops when we are at l1=l1mid-1. - + // We now compute the scaling factor for the forward recursion. double lambda = (l1midp1*thrcof[j+2]+l1mid*thrcof[j+1]+l1midm1*thrcof[j]) /(l1midp1*l1midp1+l1mid*l1mid+l1midm1*l1midm1); - + // We scale the forward recursion. for (std::vector::iterator it = thrcof.begin(); it != thrcof.begin()+j; ++it) { @@ -202,15 +202,15 @@ std::vector wigner3j(double l2, double l3, return thrcof; } -double wigner3j(double l1, double l2, double l3, +double wigner3j(double l1, double l2, double l3, double m1, double m2, double m3) { // We enforce the selection rules. bool select(true); - select = ( - std::fabs(m1+m2+m3)<1.0e-10 + select = ( + std::fabs(m1+m2+m3)<1.0e-10 && std::floor(l1+l2+l3)==(l1+l2+l3) - && l3 >= std::fabs(l1-l2) + && l3 >= std::fabs(l1-l2) && l3 <= l1+l2 && std::fabs(m1) <= l1 && std::fabs(m2) <= l2 @@ -231,18 +231,18 @@ double wigner3j(double l1, double l2, double l3, std::vector wigner6j(double l2, double l3, double l4, double l5, double l6) { - // We compute the numeric limits of double precision. + // We compute the numeric limits of double precision. double huge = std::numeric_limits::max(); double srhuge = sqrt(huge); double tiny = std::numeric_limits::min(); double srtiny = sqrt(tiny); double eps = std::numeric_limits::epsilon(); - // We enforce the selection rules. + // We enforce the selection rules. bool select(true); // Triangle relations for the four tryads - select = ( + select = ( std::fabs(l4-l2) <= l6 && l6 <= l4+l2 && std::fabs(l4-l5) <= l3 && l3 <= l4+l5 ); @@ -273,7 +273,7 @@ std::vector wigner6j(double l2, double l3, // Otherwise, we start the forward recursion. else { - // We start with an arbitrary value. + // We start with an arbitrary value. sixcof[0] = srtiny; // From now on, we check the variation of |alpha(l1)|. @@ -291,10 +291,10 @@ std::vector wigner6j(double l2, double l3, if (size>2) { - // We start with an arbitrary value. + // We start with an arbitrary value. sixcof[0] = srtiny; - - // From now on, we check the variation of |alpha(l1)|. + + // From now on, we check the variation of |alpha(l1)|. double alphaOld, alphaNew, beta, l1(l1min); if (l1min==0) alphaNew = -(l2*(l2+1.0)+l3*(l3+1.0)+l5*(l5+1.0)+l6*(l6+1.0)-2.0*l4*(l4+1.0))/wigner6j_auxA(1.0,l2,l3,l4,l5,l6); @@ -302,10 +302,10 @@ std::vector wigner6j(double l2, double l3, else alphaNew = -wigner6j_auxB(l1,l2,l3,l4,l5,l6) /(l1min*wigner6j_auxA(l1+1.0,l2,l3,l4,l5,l6)); - + // We compute the two-term recursion. sixcof[1] = alphaNew*sixcof[0]; - + // We compute the rest of the recursion. unsigned int i = 1; bool alphaVar = false; @@ -315,17 +315,17 @@ std::vector wigner6j(double l2, double l3, i++; // Next term in recursion alphaOld = alphaNew; // Monitoring of |alpha(l1)|. l1 += 1.0; // l1 = l1+1 - - // New coefficients in recursion. + + // New coefficients in recursion. alphaNew = -wigner6j_auxB(l1,l2,l3,l4,l5,l6) /(l1*wigner6j_auxA(l1+1.0,l2,l3,l4,l5,l6)); - + beta = -(l1+1.0)*wigner6j_auxA(l1,l2,l3,l4,l5,l6) /(l1*wigner6j_auxA(l1+1.0,l2,l3,l4,l5,l6)); - - // Application of the recursion. + + // Application of the recursion. sixcof[i] = alphaNew*sixcof[i-1]+beta*sixcof[i-2]; - + // We check if we are overflowing. if (std::fabs(sixcof[i]>srhuge)) { @@ -335,53 +335,53 @@ std::vector wigner6j(double l2, double l3, *it /= srhuge; } } - + // This piece of code checks whether we have reached - // the classical region. If we have, the second if + // the classical region. If we have, the second if // sets alphaVar to true and we break this loop at the - // next iteration because we need sixcof(l1mid+1) to + // next iteration because we need sixcof(l1mid+1) to // compute the scalar. if (alphaVar) break; - + if (std::fabs(alphaNew)-std::fabs(alphaOld)>0.0) alphaVar=true; - + } while (i<(size-1)); // Loop stops when we have computed all values. - + // If this is the case, we have stumbled upon a classical region. // We start the backwards recursion. if (i!=size-1) { // We keep the two terms around l1mid to compute the factor later. double l1midm1(sixcof[i-2]),l1mid(sixcof[i-1]),l1midp1(sixcof[i]); - - // We compute the backward recursion by providing an arbitrary + + // We compute the backward recursion by providing an arbitrary // startint value. sixcof[size-1] = srtiny; - + // We compute the two-term recursion. l1 = l1max; alphaNew = -wigner6j_auxB(l1,l2,l3,l4,l5,l6) /((l1+1.0)*wigner6j_auxA(l1,l2,l3,l4,l5,l6)); sixcof[size-2] = alphaNew*sixcof[size-1]; - + // We compute the rest of the backward recursion. unsigned int j = size-2; do { // Bookkeeping - j--; // Previous term in recursion. + j--; // Previous term in recursion. l1 -= 1.0; // l1 = l1-1 - - // New coefficients in recursion. + + // New coefficients in recursion. alphaNew = -wigner6j_auxB(l1,l2,l3,l4,l5,l6) /((l1+1.0)*wigner6j_auxA(l1,l2,l3,l4,l5,l6)); beta = -l1*wigner6j_auxA(l1+1.0,l2,l3,l4,l5,l6) /((l1+1.0)*wigner6j_auxA(l1,l2,l3,l4,l5,l6)); - - // Application of the recursion. + + // Application of the recursion. sixcof[j] = alphaNew*sixcof[j+1]+beta*sixcof[j+2]; - + // We check if we are overflowing. if (std::fabs(sixcof[j]>srhuge)) { @@ -391,13 +391,13 @@ std::vector wigner6j(double l2, double l3, *it /= srhuge; } } - + } while (j>(i-2)); // Loop stops when we are at l1=l1mid-1. - + // We now compute the scaling factor for the forward recursion. double lambda = (l1midp1*sixcof[j+2]+l1mid*sixcof[j+1]+l1midm1*sixcof[j]) /(l1midp1*l1midp1+l1mid*l1mid+l1midm1*l1midm1); - + // We scale the forward recursion. for (std::vector::iterator it = sixcof.begin(); it != sixcof.begin()+j; ++it) { @@ -422,14 +422,14 @@ std::vector wigner6j(double l2, double l3, return sixcof; } -double wigner6j(double l1, double l2, double l3, +double wigner6j(double l1, double l2, double l3, double l4, double l5, double l6) { - // We enforce the selection rules. + // We enforce the selection rules. bool select(true); // Triangle relations for the four tryads - select = ( + select = ( std::fabs(l1-l2) <= l3 && l3 <= l1+l2 && std::fabs(l1-l5) <= l6 && l6 <= l1+l5 && std::fabs(l4-l2) <= l6 && l6 <= l4+l2 @@ -463,7 +463,7 @@ double wigner3j_auxA(double l1, double l2, double l3, return sqrt(T1*T2*T3); } -double wigner3j_auxB(double l1, double l2, double l3, +double wigner3j_auxB(double l1, double l2, double l3, double m1, double m2, double m3) { double T1 = -(2.0*l1+1.0); diff --git a/src/wignerSymbols-fortran.cpp b/src/wignerSymbols-fortran.cpp index b8b8cb6..34e503b 100644 --- a/src/wignerSymbols-fortran.cpp +++ b/src/wignerSymbols-fortran.cpp @@ -3,7 +3,7 @@ * Lesser Public License. If a copy of the LGPL was not -/ * distributed with this file, you can obtain one at -/ * https://www.gnu.org/licenses/lgpl.html. -/ - ********************************************************/ + ********************************************************/ #include "../include/wignerSymbols/wignerSymbols-fortran.h" @@ -13,7 +13,7 @@ namespace WignerSymbols { /*! Computes a string of Wigner-3j symbols for given l2, l3, m1, m2, m3. */ std::vector wigner3j_f(double l2, double l3, double m1, double m2, double m3) { - // We prepare the size of the resulting array. + // We prepare the size of the resulting array. int size = (int)std::ceil(l2+l3-std::max(std::fabs(l2-l3),std::fabs(m1)))+1; // We prepare the output values. @@ -33,15 +33,15 @@ std::vector wigner3j_f(double l2, double l3, double m1, double m2, doubl /*! Computes the Wigner-3j symbol for given l1,l2,l3,m1,m2,m3. We * explicitly enforce the selection rules. */ -double wigner3j_f(double l1, double l2, double l3, +double wigner3j_f(double l1, double l2, double l3, double m1, double m2, double m3) { // We enforce the selection rules. bool select(true); - select = ( - m1+m2+m3==0.0 + select = ( + m1+m2+m3==0.0 && std::floor(l1+l2+l3)==(l1+l2+l3) - && l3 >= std::fabs(l1-l2) + && l3 >= std::fabs(l1-l2) && l3 <= l1+l2 && std::fabs(m1) <= l1 && std::fabs(m2) <= l2 @@ -58,22 +58,22 @@ double wigner3j_f(double l1, double l2, double l3, double thrcof [size]; int ierr; - // External function call. + // External function call. drc3jj_wrap(l2,l3,m2,m3,&l1min,&l1max,thrcof,size,&ierr); // We fetch and return the value with the proper l1 value. - int index = (int)l1-l1min; + int index = (int)l1-l1min; return thrcof[index]; } -double wigner6j_f(double l1, double l2, double l3, +double wigner6j_f(double l1, double l2, double l3, double l4, double l5, double l6) { - // We enforce the selection rules. + // We enforce the selection rules. bool select(true); // Triangle relations for the four tryads - select = ( + select = ( std::fabs(l1-l2) <= l3 && l3 <= l1+l2 && std::fabs(l1-l5) <= l6 && l6 <= l1+l5 && std::fabs(l4-l2) <= l6 && l6 <= l4+l2 @@ -90,7 +90,7 @@ double wigner6j_f(double l1, double l2, double l3, if (!select) return 0.0; - // We compute the size of the resulting array. + // We compute the size of the resulting array. int size = (int)std::ceil(std::min(l2+l3,l5+l6)-std::max(std::fabs(l2-l3),std::fabs(l5-l6)))+1; // We prepare the output values @@ -101,7 +101,7 @@ double wigner6j_f(double l1, double l2, double l3, // External function call drc6j_wrap(l2,l3,l4,l5,l6,&l1min,&l1max,sixcof,size,&ierr); - // We fetch and return the coefficient with the proper l1 value. + // We fetch and return the coefficient with the proper l1 value. int index = (int)l1-l1min; return sixcof[index]; } diff --git a/src/wignerSymbols-fortran.f b/src/wignerSymbols-fortran.f index afdfa3e..dbdc8db 100644 --- a/src/wignerSymbols-fortran.f +++ b/src/wignerSymbols-fortran.f @@ -3,7 +3,7 @@ C Lesser Public License. If a copy of the LGPL was not -/ C distributed with this file, you can obtain one at -/ C https://www.gnu.org/licenses/lgpl.html. -/ -C ******************************************************/ +C ******************************************************/ *DECK DRC3JJ SUBROUTINE DRC3JJ (L2, L3, M2, M3, L1MIN, L1MAX, THRCOF, NDIM,