anchor.diff:          L=L+1
           RY0(L)=X
   C       XXX: should have a better initial guess for large N ~> 100 here
  -        X=X+PI+MAX((0.312+0.0852*N-0.000403*N**2)/L,0d0)
  +        X=X+PI+MAX((0.312d0+0.0852*N-0.000403*N**2)/L,0d0)
           IF (L.LT.NT) GO TO 20
   C       -- Newton method for y_{N,L}'
           IF (N.LE.20) THEN
  
positive.code:         L=L+1
          RY1(L)=X
  C       XXX: should have a better initial guess for large N ~> 100 here
          X=X+PI+MAX((0.197+0.0643*N-0.000286*N**2)/L,0d0)
          IF (L.LT.NT) GO TO 25
          RETURN
          END
  
negative.code:         L=L+1
          RJ1(L)=X
  C       XXX: should have a better initial guess for large N ~> 100 here
          X=X+PI+MAX((0.4955d0+0.0915*N-0.000435*N**2)/L, 0d0)
          IF (L.LT.NT) GO TO 15
  C       -- Newton method for y_{N,L}
          IF (N.LE.20) THEN
  