updates??

src/madness/bspline/CMakeLists.txt

@@ -1,5 +1,5 @@
 set (MADBSPLINE_SOURCES
-  bspline.cc
+  bspline.cc madbessel.cc solidharmonics.cc
   )
 
 set (MADBSPLINE_HEADERS
@@ -8,6 +8,8 @@ set (MADBSPLINE_HEADERS
   )
 
 add_mad_executable(bspline bspline.cc "MADchem;-lqd")
+add_mad_executable(solidharmonics solidharmonics.cc "MADchem;-lqd")
+add_mad_executable(madbessel madbessel.cc "MADchem;-lqd")
 ##add_dependencies(applications-madness bspline)
 
 

src/madness/bspline/bspline.cc

@@ -1074,24 +1074,24 @@ class BsplineFunction {
         }
     }
 
-    static void test_fit() {
-        const T Z = 10.0;
-        size_t nknots = 61;
-        size_t order = 20;
-        T xlo = 0.1; // a for KnotsRational //.01 for 103 //1e-7;
-        T xhi = 26.0;
+    // static void test_fit() {
+    //     const T Z = 10.0;
+    //     size_t nknots = 61;
+    //     size_t order = 20;
+    //     T xlo = 0.1; // a for KnotsRational //.01 for 103 //1e-7;
+    //     T xhi = 26.0;
 
-        bdataT::init(order, nknots, xlo, xhi);
+    //     bdataT::init(order, nknots, xlo, xhi);
 
-        print("knots", bdataT::knots());
+    //     print("knots", bdataT::knots());
 
-        auto A = bdataT::basis_at_GL_points(); 
-        auto [X, W] = bdataT::quadrature();
-        size_t ltest = 1;
-        ?????????????????????????????
+    //     auto A = bdataT::basis_at_GL_points(); 
+    //     auto [X, W] = bdataT::quadrature();
+    //     size_t ltest = 1;
+    //     ?????????????????????????????
 
 
-    }
+    // }
 
     static void test() {
         const T Z = 10.0;

src/madness/bspline/myqd.h

@@ -14,6 +14,8 @@
 // Inject math routines into std name space to enable templated types
 namespace std {
     qd_real abs(const qd_real& x) {return ::abs(x);}
+    qd_real exp(const qd_real& x) {return ::exp(x);}
+    qd_real ldexp(const qd_real &a, int exp) {return ::ldexp(a,exp);}
     qd_real cos(const qd_real& x) {return ::cos(x);}
     qd_real sin(const qd_real& x) {return ::sin(x);}
     qd_real pow(const qd_real& x, int n) {return ::pow(x,n);}
@@ -22,6 +24,8 @@ namespace std {
     qd_real sqrt(const qd_real& x) {return ::sqrt(x);}
 
     dd_real abs(const dd_real& x) {return ::abs(x);}
+    dd_real exp(const dd_real& x) {return ::exp(x);}
+    dd_real ldexp(const dd_real &a, int exp) {return ::ldexp(a,exp);}
     dd_real cos(const dd_real& x) {return ::cos(x);}
     dd_real sin(const dd_real& x) {return ::sin(x);}
     dd_real pow(const dd_real& x, int n) {return ::pow(x,n);}
@@ -30,16 +34,87 @@ namespace std {
     dd_real sqrt(const dd_real& x) {return ::sqrt(x);}
 }
 
+template <typename T>
+T myread(const std::string& s) {
+    T x = 0;
+    int sx=1, sexpnt=1, expnt=0;
+    int n = 0;
+    int point = 0;
+    bool doingexpt = false;
+    for (char c : s) {
+        if (doingexpt) {
+            switch (c) {
+            case '-':
+                sexpnt = -1;
+                break;
+
+            case '+':
+                break;
+
+            default:
+                expnt = 10*expnt + int(c-'0');
+            }
+
+        }
+        else {
+              switch (c) {
+              case '-':
+                  sx = -1;
+                  break;
+
+              case '.':
+                  point = n;
+                  break;
+
+              case 'e':
+              case 'E':
+                  doingexpt = true;
+                  break;
+
+              default:
+                  x = T(10)*x + int(c-'0');
+                  n++;
+              }
+        }
+    }
+    //std::cout << n << " " << expnt << " " << point << std::endl;
+    n += -sexpnt*expnt - point;
+    n = -n;
+
+    while (abs(n)) {
+        int m;
+        if (n<0) {
+            m = std::max(n, -20);
+        }
+        else {
+            m = std::min(n, 20);
+        }
+        n -= m;
+        x *= pow(T(10), m);
+    }
+
+    if (sx<0) {
+        x = -x;
+    }
+
+    return x;
+}
+
+template <typename T>
+T myread(const char* s) {
+    return myread<T>(std::string(s));
+}
+
 template <typename T> T from_str(const char* s);
 
 template <>
 qd_real from_str<qd_real>(const char* s) {
-    return qd_real(s);
+    return myread<qd_real>(s);
 }
 
 template <>
 dd_real from_str<dd_real>(const char* s) {
-    return dd_real(s);
+    return myread<dd_real>(s);
 }
 
 template <>
@@ -60,12 +135,16 @@ float from_str<float>(const char* s) {
     return d;
 }
 
+template <typename T> T from_str(const std::string& s) {
+    return from_str<T>(s.c_str());
+}
+
 std::string to_str(const qd_real& t) {return t.to_string();}
 std::string to_str(const dd_real& t) {return t.to_string();}
-std::string to_str(double t) {char buf[256]; sprintf(buf,"%.17f",t); return buf;}
+std::string to_str(double t) {char buf[256]; sprintf(buf,"%.17e",t); return buf;}
+std::string to_str(float t) {char buf[256]; sprintf(buf,"%.8e",t); return buf;}
 
 double to_double(double d) {return d;}
-
 double to_double(float f) {return f;}
 
 // Needed since double(qd_real) not implemented and I fear what would break if I added it.

src/madness/bspline/solidharmonics.cc

@@ -546,7 +546,7 @@ auto threeJ() {
                                 //std::cout << std::setw(3) << l1 << " " << std::setw(3) << m1 << " " << std::setw(3) << l2 << " " << std::setw(3) << m2 << "   --> " << std::setw(3) << l3 << " " << std::setw(3) << m3 << " " << std::setw(15) << j3value << " " << test << std::endl;
                                 j3.push_back({l3,m3,j3value});
 
-                                if (l2==1 and m2==0) {
+                                if ((l2==1 and m2==1) or (l1==1 and m2==1)) {
                                     T unnorm = j3value*Nnorm.value(l3,m3)/(Nnorm.value(l1,m1)*Nnorm.value(l2,m2));
                                     std::cout << std::setw(3) << l1 << " " << std::setw(3) << m1 << " " << std::setw(3) << l2 << " " << std::setw(3) << m2 << "   --> " << std::setw(3) << l3 << " " << std::setw(3) << m3 << " " << std::setw(15) << j3value << " " << unnorm << std::endl;
                                 }