special: test for formfactor singularities

pub/lib/special.cpp

@@ -71,48 +71,31 @@ void NSpecial::FourierCosine()
 
 void NSpecial::Test()
 {
-    // diagonalize finite chain force matrix
+    // relative error in form factor
 
     POld fout( new COld );
 
-    int N = iask("Order N");
-    if ( N<2 )
-         return;
+    fout->name = "TEST_";
+    fout->xco = CCoord("a", "");
+    fout->yco = CCoord("FF", "");
+    fout->ZCo.push_back( CCoord("b", "") );
 
-    fout->name = "EV_" +S(N);
-    fout->xco = CCoord("n", "");
-    fout->yco = CCoord("v_n", "");
-    fout->ZCo.push_back( CCoord("p", "") );
-    fout->ZCo.push_back( CCoord("nor", "") );
-
-    for ( int p=0; p<N; ++p ){
+    static int Nj = 100;
+    static int Ni = 1000;
+    
+    for ( int j=0; j<Nj; ++j ){
+        double b = pow( 1.2, j-Nj );
         PSpec s( new CSpec );
-
-        s->z.push_back( PObjInt( new CObjInt (p) ) );
-
-        s->resize(N,false);
-        for( int n=0; n<N; ++n )
-            s->x[n] = n;
-
-        /// eigenvalue:
-        double lambda = -2 + 2*cos(p*M_PI/N);
-
-        // recursive calculation of eigenvector:
-        s->y[0] = 1;
-        s->y[1] = 1+lambda;
-        for( int n=2; n<N; ++n )
-            s->y[n] = -s->y[n-2] + (2+lambda)*s->y[n-1];
-
-        // normalize:
-        double nor = 0;
-        for( int n=0; n<N; ++n )
-            nor += s->y[n] * s->y[n];
-        printf( "p=%i => check sum = %g, norm = %g\n",
-                p, s->y[N-2]-(1+lambda)*s->y[N-1], nor );
-        s->z.push_back( PObjDbl( new CObjDbl( nor ) ) );
-        for( int n=0; n<N; ++n )
-            s->y[n] /= sqrt(nor);
-
+        s->z.push_back( PObjInt( new CObjInt (b) ) );
+        s->resize(Ni,false);
+        for ( int i=0; i<Ni; ++i ){
+            double a = b + pow( 1.2, i-Ni );
+            s->x[i] = a;
+            complex<long double> ff1 = ff_test( (long double) a, (long double) b );
+            complex<double> ff2 = ff_test( a, b );
+            complex<double> diff = (complex<double>)(ff1)-ff2;
+            s->y[i] = std::abs(diff)/std::abs(ff1);
+        }
         fout->V.push_back( move(s) );
     }
     SMem::instance()->mem_store( move(fout) );

pub/lib/special.hpp

@@ -7,12 +7,26 @@
 //! \file  special.hpp
 //! \brief NSpecial: this and that, written ad hoc.
 
-
 //  What falls out of use goes to ../../arch/one-time-routines/
 
-//! This and that, written ad hoc.
+#include <complex>
+using std::complex;
 
 namespace NSpecial {
     void FourierCosine();
     void Test();
-}
\ No newline at end of file
+
+    template<typename T>
+    complex<T> pssc( T z )
+    {
+        complex<T> ci( 0.L, 1.L );
+        complex<T> cone( 1.L, 0.L );
+        return z==0 ? -ci : (cone-std::exp(ci*z))/z;
+    }
+
+    template<typename T>
+    complex<T> ff_test( T a, T b )
+    {
+        return complex<T>(T(1.L)/b,T(0.L)) * ( pssc(a+b) - pssc(a-b) );
+    }    
+}