Wigner3jm.cpp

/*
 * Matpack Wigner3jm special function imported and modified for use in
 * OpenMD
 *
 * Matpack Library Release 1.9.0
 * Copyright (C) 1991-2003 by Berndt M. Gammel. All rights reserved.
 *
 * Permission to use, copy, and distribute Matpack in its entirety
 * and its documentation for non-commercial purpose and without fee
 * is hereby granted, provided that this license information and
 * copyright notice appear unmodified in all copies.  This software
 * is provided 'as is' without express or implied warranty.  In no
 * event will the author be held liable for any damages arising from
 * the use of this software.
 *
 * Note that distributing Matpack 'bundled' in with any product is
 * considered to be a 'commercial purpose'.
 *
 * The software may be modified for your own purposes, but modified
 * versions may not be distributed without prior consent of the
 * author.
 *
 * Read the COPYRIGHT and README files in this distribution about
 * registration and installation of Matpack.
 */
 
#include "Wigner3jm.hpp"
 
#include <cfloat>
#include <cmath>
#include <cstdio>
 
#include "utils/simError.h"
 
namespace MATPACK {
 
  //-----------------------------------------------------------------------------//
  //
  // void ThreeJSymbolM (RealType l1, RealType l2, RealType l3, RealType m1,
  //                     RealType &m2min, RealType &m2max, RealType *thrcof, int
  //                     ndim, int &errflag)
  //
  // Evaluate the Wigner 3j symbol
  //
  //       g(m2) = ( l1  l2     l3  )
  //               ( m1  m2  -m1-m2 )
  //
  // for all allowed values of m2, the other parameters being held fixed.
  //
  // Input Arguments:
  // ----------------
  //
  //   RealType l1
  //   RealType l2
  //   RealType l3
  //   RealType m1        Parameters in 3j symbol.
  //
  //   int  ndim          Declared length of thrcof in calling program.
  //
  // Output Arguments:
  // -----------------
  //
  //   RealType &m2min    Smallest allowable m2 in 3j symbol.
  //   RealType &m2max    Largest allowable m2 in 3j symbol.
  //   RealType *thrcof   Set of 3j coefficients generated by evaluating the
  //                      3j symbol for all allowed values of m2.  thrcof(i)
  //                      will contain g(m2min+i), i=0,2,...,m2max-m2min.
  //
  //   int &errflag       Error flag.
  //                      errflag=0  No errors.
  //                      errflag=1  Either l1 < abs(m1) or l1+abs(m1)
  //                      non-integer. errflag=2  abs(l1-l2)<= l3 <= l1+l2 not
  //                      satisfied. errflag=3  l1+l2+l3 not an integer.
  //                      errflag=4 m2max-m2min not an integer. errflag=5  m2max
  //                      less than m2min. errflag=6  ndim less than
  //                      m2max-m2min+1.
  // Description:
  // ------------
  //
  // Although conventionally the parameters of the vector addition
  // coefficients satisfy certain restrictions, such as being integers
  // or integers plus 1/2, the restrictions imposed on input to this
  // subroutine are somewhat weaker. See, for example, Section 27.9 of
  // Abramowitz and Stegun or Appendix C of Volume II of A. Messiah.
  //
  // The restrictions imposed by this subroutine are
  //
  //       1. l1 >= abs(m1) and l1+abs(m1) must be an integer
  //       2. abs(l1-l2) <= l3 <= l1+l2
  //       3. l1+l2+l3 must be an integer
  //       4. m2max-m2min must be an integer, where
  //          m2max=min(l2,l3-m1) and m2min=max(-l2,-l3-m1)
  //
  // If the conventional restrictions are satisfied, then these
  // restrictions are also met.
  //
  // The user should be cautious in using input parameters that do
  // not satisfy the conventional restrictions. For example, the
  // the subroutine produces values of
  //       g(m2) = (0.75 1.50   1.75  )
  //               (0.25  m2  -0.25-m2)
  // for m2=-1.5,-0.5,0.5,1.5 but none of the symmetry properties of the
  // 3j symbol, set forth on page 1056 of Messiah, is satisfied.
  //
  // The subroutine generates g(m2min), g(m2min+1), ..., g(m2max)
  // where m2min and m2max are defined above. The sequence g(m2) is
  // generated by a three-term recurrence algorithm with scaling to
  // control overflow. Both backward and forward recurrence are used to
  // maintain numerical stability. The two recurrence sequences are
  // matched at an interior point and are normalized from the unitary
  // property of 3j coefficients and Wigner's phase convention.
  //
  // The algorithm is suited to applications in which large quantum
  // numbers arise, such as in molecular dynamics.
  //
  // References:
  // -----------
  //  1. Abramowitz, M., and Stegun, I. A., Eds., Handbook
  //     of Mathematical Functions with Formulas, Graphs
  //     and Mathematical Tables, NBS Applied Mathematics
  //     Series 55, June 1964 and subsequent printings.
  //  2. Messiah, Albert., Quantum Mechanics, Volume II,
  //     North-Holland Publishing Company, 1963.
  //  3. Schulten, Klaus and Gordon, Roy G., Exact recursive
  //     evaluation of 3j and 6j coefficients for quantum-
  //     mechanical coupling of angular momenta, J Math
  //     Phys, v 16, no. 10, October 1975, pp. 1961-1970.
  //  4. Schulten, Klaus and Gordon, Roy G., Semiclassical
  //     approximations to 3j  and 6j coefficients for
  //     quantum-mechanical coupling of angular momenta,
  //     J Math Phys, v 16, no. 10, October 1975, pp. 1971-1988.
  //  5. Schulten, Klaus and Gordon, Roy G., Recursive
  //     evaluation of 3j and 6j coefficients, Computer
  //     Phys Comm, v 11, 1976, pp. 269-278.
  //  6. SLATEC library, category  C19,
  //     double precision algorithm DRC3JM.F
  //     Keywords: 3j coefficients, 3j symbols, Clebsch-Gordan coefficients,
  //               Racah coefficients, vector addition coefficients,
  //               Wigner coefficients
  //     Author:   Gordon, R. G., Harvard University
  //               Schulten, K., Max Planck Institute
  //     Revision history  (YYMMDD)
  //     750101  DATE WRITTEN
  //     880515  SLATEC prologue added by G. C. Nielson, NBS; parameters
  //             HUGE and TINY revised to depend on D1MACH.
  //     891229  Prologue description rewritten; other prologue sections
  //             revised; MMATCH (location of match point for recurrences)
  //             removed from argument list; argument IER changed to serve
  //             only as an error flag (previously, in cases without error,
  //             it returned the number of scalings); number of error codes
  //             increased to provide more precise error information;
  //             program comments revised; SLATEC error handler calls
  //             introduced to enable printing of error messages to meet
  //             SLATEC standards. These changes were done by D. W. Lozier,
  //             M. A. McClain and J. M. Smith of the National Institute
  //             of Standards and Technology, formerly NBS.
  //     910415  Mixed type expressions eliminated; variable C1 initialized;
  //             description of THRCOF expanded. These changes were done by
  //             D. W. Lozier.
  //  7. Rewritting of the SLATEX algorithm in C++ and adaption to the
  //     Matpack C++ Numerics and Graphics Library by Berndt M. Gammel
  //     in June 1997.
  //
  //-----------------------------------------------------------------------------//
 
  void Wigner3jm(RealType l1, RealType l2, RealType l3, RealType m1,
                 RealType& m2min, RealType& m2max, RealType* thrcof, int ndim,
                 int& errflag) {
    // In single precision, the largest floating point number is not
    // the same as in double precision:
#ifdef SINGLE_PRECISION
    RealType MaxFloat = FLT_MAX;
#else
    RealType MaxFloat = DBL_MAX;
#endif
 
    const RealType zero = 0.0, eps = 0.01, one = 1.0, two = 2.0;
 
    int nfin, nlim, i, n, index, lstep, nfinp1, nfinp2, nfinp3, nstep2;
    RealType oldfac, dv, newfac,
        sumbac = 0.0, thresh, a1s, sumfor, sumuni, sum1, sum2, x, y, m2, m3, x1,
        x2, x3, y1, y2, y3, cnorm, ratio, a1, c1, c2, c1old = 0.0, sign1, sign2;
 
    // Parameter adjustments
    --thrcof;
 
    errflag = 0;
 
    // "hugeRealType" is the square root of one twentieth of the
    // largest floating point number, approximately.
 
    RealType hugeRealType = sqrt(MaxFloat / 20.0), srhuge = sqrt(hugeRealType),
             tiny = one / hugeRealType, srtiny = one / srhuge;
 
    // lmatch = zero
 
    //  Check error conditions 1, 2, and 3.
    if (l1 - fabs(m1) + eps < zero ||
        fmod(l1 + fabs(m1) + eps, one) >= eps + eps) {
      errflag = 1;
 
      snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH, "%s: %s",
               "ThreeJSymbolM",
               "l1-abs(m1) less than zero or l1+abs(m1) not integer.");
      painCave.isFatal = 1;
      simError();
 
      return;
    } else if (l1 + l2 - l3 < -eps || l1 - l2 + l3 < -eps ||
               -(l1) + l2 + l3 < -eps) {
      errflag = 2;
 
      snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH, "%s: %s",
               "ThreeJSymbolM",
               "l1, l2, l3 do not satisfy triangular condition.");
      painCave.isFatal = 1;
      simError();
 
      return;
    } else if (fmod(l1 + l2 + l3 + eps, one) >= eps + eps) {
      errflag = 3;
 
      snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH, "%s: %s",
               "ThreeJSymbolM", "l1+l2+l3 not integer.");
      painCave.isFatal = 1;
      simError();
 
      return;
    }
 
    // limits for m2
    m2min = MpMax(-l2, -l3 - m1);
    m2max = MpMin(l2, l3 - m1);
 
    // Check error condition 4.
    if (fmod(m2max - m2min + eps, one) >= eps + eps) {
      errflag = 4;
 
      snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH, "%s: %s",
               "ThreeJSymbolM", "m2max-m2min not integer.");
      painCave.isFatal = 1;
      simError();
 
      return;
    }
    if (m2min < m2max - eps) goto L20;
    if (m2min < m2max + eps) goto L10;
 
    //  Check error condition 5.
    errflag = 5;
 
    snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH, "%s: %s",
             "ThreeJSymbolM", "m2min greater than m2max.");
    painCave.isFatal = 1;
    simError();
 
    return;
 
    // This is reached in case that m2 and m3 can take only one value.
  L10:
    // mscale = 0
    thrcof[1] = (odd(int(fabs(l2 - l3 - m1) + eps)) ? -one : one) /
                sqrt(l1 + l2 + l3 + one);
    return;
 
    // This is reached in case that M1 and M2 take more than one
    // value.
  L20:
    // mscale = 0
    nfin = int(m2max - m2min + one + eps);
    if (ndim - nfin >= 0) goto L23;
 
    // Check error condition 6.
 
    errflag = 6;
    snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH, "%s: %s",
             "ThreeJSymbolM",
             "Dimension of result array for 3j coefficients too small.");
    painCave.isFatal = 1;
    simError();
 
    return;
 
    //  Start of forward recursion from m2 = m2min
 
  L23:
    m2        = m2min;
    thrcof[1] = srtiny;
    newfac    = 0.0;
    c1        = 0.0;
    sum1      = tiny;
 
    lstep = 1;
  L30:
    ++lstep;
    m2 += one;
    m3 = -m1 - m2;
 
    oldfac = newfac;
    a1     = (l2 - m2 + one) * (l2 + m2) * (l3 + m3 + one) * (l3 - m3);
    newfac = sqrt(a1);
 
    dv = (l1 + l2 + l3 + one) * (l2 + l3 - l1) -
         (l2 - m2 + one) * (l3 + m3 + one) - (l2 + m2 - one) * (l3 - m3 - one);
 
    if (lstep - 2 > 0) c1old = fabs(c1);
 
    // L32:
    c1 = -dv / newfac;
 
    if (lstep > 2) goto L60;
 
    //  If m2 = m2min + 1, the third term in the recursion equation
    //  vanishes, hence
 
    x         = srtiny * c1;
    thrcof[2] = x;
    sum1 += tiny * c1 * c1;
    if (lstep == nfin) goto L220;
    goto L30;
 
  L60:
    c2 = -oldfac / newfac;
 
    // Recursion to the next 3j coefficient
    x             = c1 * thrcof[lstep - 1] + c2 * thrcof[lstep - 2];
    thrcof[lstep] = x;
    sumfor        = sum1;
    sum1 += x * x;
    if (lstep == nfin) goto L100;
 
    // See if last unnormalized 3j coefficient exceeds srhuge
 
    if (fabs(x) < srhuge) goto L80;
 
    // This is reached if last 3j coefficient larger than srhuge, so
    // that the recursion series thrcof(1), ... , thrcof(lstep) has to
    // be rescaled to prevent overflow
 
    // mscale = mscale + 1
    for (i = 1; i <= lstep; ++i) {
      if (fabs(thrcof[i]) < srtiny) thrcof[i] = zero;
      thrcof[i] /= srhuge;
    }
    sum1 /= hugeRealType;
    sumfor /= hugeRealType;
    x /= srhuge;
 
    // As long as abs(c1) is decreasing, the recursion proceeds
    // towards increasing 3j values and, hence, is numerically stable.
    // Once an increase of abs(c1) is detected, the recursion
    // direction is reversed.
 
  L80:
    if (c1old - fabs(c1) > 0.0) goto L30;
 
    //  Keep three 3j coefficients around mmatch for comparison later
    //  with backward recursion values.
 
  L100:
    // mmatch = m2 - 1
    nstep2 = nfin - lstep + 3;
    x1     = x;
    x2     = thrcof[lstep - 1];
    x3     = thrcof[lstep - 2];
 
    //  Starting backward recursion from m2max taking nstep2 steps, so
    //  that forwards and backwards recursion overlap at the three
    //  points m2 = mmatch+1, mmatch, mmatch-1.
 
    nfinp1       = nfin + 1;
    nfinp2       = nfin + 2;
    nfinp3       = nfin + 3;
    thrcof[nfin] = srtiny;
    sum2         = tiny;
 
    m2    = m2max + two;
    lstep = 1;
  L110:
    ++lstep;
    m2 -= one;
    m3     = -m1 - m2;
    oldfac = newfac;
    a1s = (l2 - m2 + two) * (l2 + m2 - one) * (l3 + m3 + two) * (l3 - m3 - one);
    newfac = sqrt(a1s);
    dv     = (l1 + l2 + l3 + one) * (l2 + l3 - l1) -
         (l2 - m2 + one) * (l3 + m3 + one) - (l2 + m2 - one) * (l3 - m3 - one);
    c1 = -dv / newfac;
    if (lstep > 2) goto L120;
 
    // if m2 = m2max + 1 the third term in the recursion equation
    // vanishes
 
    y                = srtiny * c1;
    thrcof[nfin - 1] = y;
    if (lstep == nstep2) goto L200;
    sumbac = sum2;
    sum2 += y * y;
    goto L110;
 
  L120:
    c2 = -oldfac / newfac;
 
    // Recursion to the next 3j coefficient
 
    y = c1 * thrcof[nfinp2 - lstep] + c2 * thrcof[nfinp3 - lstep];
 
    if (lstep == nstep2) goto L200;
 
    thrcof[nfinp1 - lstep] = y;
    sumbac                 = sum2;
    sum2 += y * y;
 
    // See if last 3j coefficient exceeds SRHUGE
 
    if (fabs(y) < srhuge) goto L110;
 
    // This is reached if last 3j coefficient larger than srhuge, so
    // that the recursion series
    // thrcof(nfin), ... , thrcof(nfin-lstep+1)
    // has to be rescaled to prevent overflow.
 
    // mscale = mscale + 1
    for (i = 1; i <= lstep; ++i) {
      index = nfin - i + 1;
      if (fabs(thrcof[index]) < srtiny) thrcof[index] = zero;
      thrcof[index] /= srhuge;
    }
    sum2 /= hugeRealType;
    sumbac /= hugeRealType;
 
    goto L110;
 
    //  The forward recursion 3j coefficients x1, x2, x3 are to be
    //  matched with the corresponding backward recursion values y1,
    //  y2, y3.
 
  L200:
    y3 = y;
    y2 = thrcof[nfinp2 - lstep];
    y1 = thrcof[nfinp3 - lstep];
 
    //  Determine now ratio such that yi = ratio * xi (i=1,2,3) holds
    //  with minimal error.
 
    ratio = (x1 * y1 + x2 * y2 + x3 * y3) / (x1 * x1 + x2 * x2 + x3 * x3);
    nlim  = nfin - nstep2 + 1;
 
    if (fabs(ratio) < one) goto L211;
    for (n = 1; n <= nlim; ++n)
      thrcof[n] = ratio * thrcof[n];
    sumuni = ratio * ratio * sumfor + sumbac;
    goto L230;
 
  L211:
    ++nlim;
    ratio = one / ratio;
    for (n = nlim; n <= nfin; ++n)
      thrcof[n] = ratio * thrcof[n];
    sumuni = sumfor + ratio * ratio * sumbac;
    goto L230;
 
  L220:
    sumuni = sum1;
 
    // Normalize 3j coefficients
 
  L230:
    cnorm = one / sqrt((l1 + l1 + one) * sumuni);
 
    // Sign convention for last 3j coefficient determines overall
    // phase
 
    sign1 = sign(thrcof[nfin]);
    sign2 = odd(int(fabs(l2 - l3 - m1) + eps)) ? -one : one;
    if (sign1 * sign2 <= 0.0)
      goto L235;
    else
      goto L236;
 
  L235:
    cnorm = -cnorm;
 
  L236:
    if (fabs(cnorm) < one) goto L250;
 
    for (n = 1; n <= nfin; ++n)
      thrcof[n] = cnorm * thrcof[n];
    return;
 
  L250:
    thresh = tiny / fabs(cnorm);
    for (n = 1; n <= nfin; ++n) {
      if (fabs(thrcof[n]) < thresh) thrcof[n] = zero;
      thrcof[n] = cnorm * thrcof[n];
    }
  }
}  // namespace MATPACK