From 226a9625ffa917675fedd2d9d4af4230cd1b5725 Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (h)" <j.wuttke@fz-juelich.de>
Date: Sat, 16 Apr 2016 15:21:01 +0200
Subject: [PATCH] rm unused cbesselJ01. Import setup scripts for
dodeca/icosahedron.
---
dev-tools/README | 3 +-
dev-tools/math/cbesselJ01.cpp | 173 ---------------------------
dev-tools/math/setup_dodecahedron.py | 79 ++++++++++++
dev-tools/math/setup_icosahedron.py | 67 +++++++++++
4 files changed, 147 insertions(+), 175 deletions(-)
delete mode 100644 dev-tools/math/cbesselJ01.cpp
create mode 100755 dev-tools/math/setup_dodecahedron.py
create mode 100755 dev-tools/math/setup_icosahedron.py
diff --git a/dev-tools/README b/dev-tools/README
index 270f4a4fd5e..68607cb7499 100644
--- a/dev-tools/README
+++ b/dev-tools/README
@@ -3,8 +3,7 @@ Collection of various tools for code development.
Directory content:
git-utils - number of lines of code
log - results of performance tests
-math - reference code; code to compute coefficients
+math - code used to derive or/and to test math code
profiling - scripts to run code profiling using valgrind/gperftool
-python-bindings - collection of scripts for automatic generation of libBornAgainCore and libBornAgainFit boost-python bindings
user-api - collection of tools to generate User API description
diff --git a/dev-tools/math/cbesselJ01.cpp b/dev-tools/math/cbesselJ01.cpp
deleted file mode 100644
index 1a0f31e0dda..00000000000
--- a/dev-tools/math/cbesselJ01.cpp
+++ /dev/null
@@ -1,173 +0,0 @@
-// original source of complex Bessel code, modified to get it running,
-// plus a short main routine for interactive use.
-// It's almost C code, just one detail makes it C++.
-// JWu jan16
-
-// cbessjy.cpp -- complex Bessel functions.
-// Algorithms and coefficient values from "Computation of Special
-// Functions", Zhang and Jin, John Wiley and Sons, 1996.
-//
-// (C) 2003, C. Bond. All rights reserved.
-//
-#include <math.h>
-#include <complex.h>
-#include <stdio.h>
-#include <stdlib.h>
-static double eps = 2.2e-16;
-static double complex cii = I;
-static double complex cone = 1;
-static double complex czero = 0;
-
-int cbessj01(double complex z,double complex &cj0,double complex &cj1)
-{
- double complex z1,z2,cr,cp,cs,cp0,cq0,cp1,cq1,ct1,ct2,cu;
- double a0,w0,w1;
- int k,kz;
-
- static double a[] = {
- -7.03125e-2,
- 0.112152099609375,
- -0.5725014209747314,
- 6.074042001273483,
- -1.100171402692467e2,
- 3.038090510922384e3,
- -1.188384262567832e5,
- 6.252951493434797e6,
- -4.259392165047669e8,
- 3.646840080706556e10,
- -3.833534661393944e12,
- 4.854014686852901e14,
- -7.286857349377656e16,
- 1.279721941975975e19};
- static double b[] = {
- 7.32421875e-2,
- -0.2271080017089844,
- 1.727727502584457,
- -2.438052969955606e1,
- 5.513358961220206e2,
- -1.825775547429318e4,
- 8.328593040162893e5,
- -5.006958953198893e7,
- 3.836255180230433e9,
- -3.649010818849833e11,
- 4.218971570284096e13,
- -5.827244631566907e15,
- 9.476288099260110e17,
- -1.792162323051699e20};
- static double a1[] = {
- 0.1171875,
- -0.1441955566406250,
- 0.6765925884246826,
- -6.883914268109947,
- 1.215978918765359e2,
- -3.302272294480852e3,
- 1.276412726461746e5,
- -6.656367718817688e6,
- 4.502786003050393e8,
- -3.833857520742790e10,
- 4.011838599133198e12,
- -5.060568503314727e14,
- 7.572616461117958e16,
- -1.326257285320556e19};
- static double b1[] = {
- -0.1025390625,
- 0.2775764465332031,
- -1.993531733751297,
- 2.724882731126854e1,
- -6.038440767050702e2,
- 1.971837591223663e4,
- -8.902978767070678e5,
- 5.310411010968522e7,
- -4.043620325107754e9,
- 3.827011346598605e11,
- -4.406481417852278e13,
- 6.065091351222699e15,
- -9.833883876590679e17,
- 1.855045211579828e20};
-
- a0 = cabs(z);
- z2 = z*z;
- z1 = z;
- if (a0 == 0.0) {
- cj0 = cone;
- cj1 = czero;
- return 0;
- }
- if (creal(z) < 0.0) z1 = -z;
- if (a0 <= 12.0) {
- cj0 = cone;
- cr = cone;
- for (k=1;k<=40;k++) {
- cr *= -0.25*z2/(double)(k*k);
- cj0 += cr;
- if (cabs(cr) < cabs(cj0)*eps) break;
- }
- cj1 = cone;
- cr = cone;
- for (k=1;k<=40;k++) {
- cr *= -0.25*z2/(k*(k+1.0));
- cj1 += cr;
- if (cabs(cr) < cabs(cj1)*eps) break;
- }
- cj1 *= 0.5*z1;
- }
- else {
- if (a0 >= 50.0) kz = 8; // can be changed to 10
- else if (a0 >= 35.0) kz = 10; // " " " 12
- else kz = 12; // " " " 14
- ct1 = z1 - M_PI_4;
- cp0 = cone;
- for (k=0;k<kz;k++) {
- cp0 += a[k]*cpow(z1,-2.0*k-2.0);
- }
- cq0 = -0.125/z1;
- for (k=0;k<kz;k++) {
- cq0 += b[k]*cpow(z1,-2.0*k-3.0);
- }
- cu = csqrt(M_2_PI/z1);
- cj0 = cu*(cp0*ccos(ct1)-cq0*csin(ct1));
- ct2 = z1 - 0.75*M_PI;
- cp1 = cone;
- for (k=0;k<kz;k++) {
- cp1 += a1[k]*cpow(z1,-2.0*k-2.0);
- }
- cq1 = 0.375/z1;
- for (k=0;k<kz;k++) {
- cq1 += b1[k]*cpow(z1,-2.0*k-3.0);
- }
- cj1 = cu*(cp1*ccos(ct2)-cq1*csin(ct2));
- }
- if (creal(z) < 0.0) {
- cj1 = -cj1;
- }
- return 0;
-}
-
-
-int main (int argc, char *argv[]) {
- if( argc!=4 ) {
- printf("Usage: cbessy <nu> <Re z> <Im z>\n");
- return -1;
- }
- int nu = atoi( argv[1] );
- double zr = atof( argv[2] );
- double zi = atof( argv[3] );
- double complex z = zr + I*zi;
- double complex cj0;
- double complex cj1;
- cbessj01(z, cj0, cj1);
- double complex res;
- if( nu==0 )
- res = cj0;
- else if( nu==1 )
- res = cj1;
- else {
- printf( "Only nu=0,1 are supported\n" );
- return -1;
- }
- printf( " res = MathFunctions::crbond_bessel_J1(complex_t(%g,%g));\n"
- " EXPECT_NEAR( res.real(), %24.17g, eps*std::abs(res) );\n"
- " EXPECT_NEAR( res.imag(), %24.17g, eps*std::abs(res) );\n",
- zr, zi, creal(res), cimag(res) );
- return 0;
-}
diff --git a/dev-tools/math/setup_dodecahedron.py b/dev-tools/math/setup_dodecahedron.py
new file mode 100755
index 00000000000..ff5b2710bf5
--- /dev/null
+++ b/dev-tools/math/setup_dodecahedron.py
@@ -0,0 +1,79 @@
+#!/usr/bin/env python3
+
+# setup_dodecahedron.py:
+# Compute vertex coordinates of a regular dodecahedron.
+# Auxiliary program, used during implementation of dodecahedron form factor in BornAgain.
+# Joachim Wuttke, Forschungszentrum Jülich, 2016
+
+from math import sqrt,sin,cos,pi
+
+def distance(s,t):
+ return sqrt((s[0]-t[0])**2+(s[1]-t[1])**2+(s[2]-t[2])**2)
+
+q=sqrt(5)
+# v=pi/5
+# vv=2*v
+c=(1+q)/4 # cos(v)
+s=sqrt((5-q)/8) # sin(v)
+cc=(q-1)/4 # cos(vv)
+ss=sqrt((5+q)/8) # sin(vv)
+
+def plane_rotate(k):
+ t=p[k]
+ ci = cc
+ si = ss
+ x = ci*t[0]-si*t[1]
+ y = si*t[0]+ci*t[1]
+ p[k+1] = [x, y,t[2]]
+ p[k+4] = [x,-y,t[2]]
+ ci = -c
+ si = s
+ x = ci*t[0]-si*t[1]
+ y = si*t[0]+ci*t[1]
+ p[k+2] = [x, y,t[2]]
+ p[k+3] = [x,-y,t[2]]
+
+
+# a=pow((7*q-15)/5, 1/3.) # edge for V=1
+a=1.
+
+rb=a/(2*s)
+rd=a*c/s
+R2a=(9+3*q)/8
+b=a*sqrt(R2a-1/(4*s**2))
+d=a*sqrt(R2a-(c/s)**2)
+
+h=rb*c
+vol = 12 * 5 * h*a*b/6
+print("a=%g, R/a=%g, rb/a=%g, rd/a=%g, b/a=%g, d/a=%g" % (a, sqrt(R2a), rb/a, rd/a, b/a, d/a))
+print("h=%g volume=%20.16g" % (h, vol) )
+
+p=[None for i in range(20)]
+p[ 0] = [+rb,0.,-b]
+p[ 5] = [+rd,0.,-d]
+p[10] = [-rd,0.,+d]
+p[15] = [-rb,0.,+b]
+for m in range(4):
+ plane_rotate(m*5)
+
+for i in range(len(p)):
+ print(" {%20.16g*a,%20.16g*a,%20.16g*a}," % (p[i][0], p[i][1], p[i][2]))
+
+for i in range(20):
+ print("p%2i(%20.16g,%20.16g,%20.16g) at R=%20.16g" %
+ (i, p[i][0], p[i][1], p[i][2], distance(p[i],[0,0,0])))
+
+distcoll = {}
+for i in range(20):
+ for j in range(i+1,20):
+ dist = distance(p[i],p[j])
+ print("p%2i(%9.4g,%9.4g,%9.4g) - p%2i(%9.4g,%9.4g,%9.4g) = %20.17g" %
+ (i, p[i][0], p[i][1], p[i][2], j, p[j][0], p[j][1], p[j][2], dist))
+ dist=round(dist*1e14)/10**14
+ if not dist in distcoll:
+ distcoll[dist] = 1
+ else:
+ distcoll[dist] += 1
+
+for dist in sorted(distcoll.keys()):
+ print( "%2i times %20.16g" % (distcoll[dist],dist) )
diff --git a/dev-tools/math/setup_icosahedron.py b/dev-tools/math/setup_icosahedron.py
new file mode 100755
index 00000000000..77c552279be
--- /dev/null
+++ b/dev-tools/math/setup_icosahedron.py
@@ -0,0 +1,67 @@
+#!/usr/bin/env python3
+
+# setup_icosahedron.py:
+# Compute vertex coordinates of a regular icosahedron.
+# Auxiliary program, used during implementation of icosahedron form factor in BornAgain.
+# Joachim Wuttke, Forschungszentrum Jülich, 2016
+
+from math import sqrt,sin,cos,pi
+
+def distance(s,t):
+ return sqrt((s[0]-t[0])**2+(s[1]-t[1])**2+(s[2]-t[2])**2)
+
+c=sqrt(3)/2
+q=sqrt(5)
+ss=c
+cc=.5
+
+# a=pow(3*(3-q)/5, 1/3.) # edge for V=1
+a=1.
+
+R2=(5+q)/8*a**2
+rb=a/sqrt(3)
+b=a*sqrt((7+3*q)/24) # sqrt(R2-rb**2)
+rd2=a**2*(3+q)/6 # (4*R2-a**2)/2/(1+cc)
+rd=sqrt(rd2)
+d=a*sqrt((3-q)/24) #sqrt(R2-rd2)
+
+h=rb+sqrt(rb**2-(0.5*a)**2)
+vol = 20 * h*a*b/6
+print("a=%g, R=%g, rb=%g, rd=%g, b=%g, d=%g" % (a, sqrt(R2), rb, rd, b, d))
+print("check P14=%g" % (sqrt( (rd*cc-rb)**2 + (rd*ss)**2 + (d-b)**2 )))
+print("h=%g volume=%20.16g" % (h, vol) )
+
+p=[None for i in range(12)]
+p[ 0] = [ rb, 0, -b]
+p[ 1] = [-rb*cos(pi/3),+rb*sin(pi/3),-b]
+p[ 2] = [-rb*cos(pi/3),-rb*sin(pi/3),-b]
+p[ 3] = [-rd, 0, -d]
+p[ 4] = [ rd*cos(pi/3),+rd*sin(pi/3),-d]
+p[ 5] = [ rd*cos(pi/3),-rd*sin(pi/3),-d]
+p[ 6] = [ rd, 0, +d]
+p[ 7] = [-rd*cos(pi/3),+rd*sin(pi/3),+d]
+p[ 8] = [-rd*cos(pi/3),-rd*sin(pi/3),+d]
+p[ 9] = [-rb, 0, +b]
+p[10] = [ rb*cos(pi/3),+rb*sin(pi/3),+b]
+p[11] = [ rb*cos(pi/3),-rb*sin(pi/3),+b]
+
+for i in range(len(p)):
+ print(" {%20.16g*a,%20.16g*a,%20.16g*a}," % (p[i][0], p[i][1], p[i][2]))
+
+for i in range(len(p)):
+ print("p%2i(%20.16g,%20.16g,%20.16g) at R=%20.16g" %
+ (i, p[i][0], p[i][1], p[i][2], distance(p[i],[0,0,0])))
+
+distcoll = {}
+for i in range(len(p)):
+ for j in range(i+1,len(p)):
+ dist = distance(p[i],p[j])
+ print("p%2i-%2i = %20.17g" % (i, j, dist))
+ dist=round(dist*1e15)/10**15
+ if not dist in distcoll:
+ distcoll[dist] = 1
+ else:
+ distcoll[dist] += 1
+
+for dist in sorted(distcoll.keys()):
+ print( "%2i times %20.16g" % (distcoll[dist],dist) )
--
GitLab