Skip to content
Snippets Groups Projects
Commit cd28c047 authored by Wuttke, Joachim's avatar Wuttke, Joachim
Browse files

ff devtools moved to libformfactor

parent 4e62e705
Branches
Tags
1 merge request!552Clean up few things in directory devtools/
Pipeline #52140 passed
Stage: build
Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 1123 deletions
devtools/math/ff/plot/plotff.py deleted 100644 → 0
+ 0
− 138
View file @ 4e62e705
# script to plot the form factor
# after the ROOT window with plots appear,
# just save the picture in that format which you need
import ROOT
from libBornAgainCore import *
# define a form factor, I recommend a 10--20 nm diameter
#ff = FormFactorFullSphere(10.0*nanometer)
#ff = FormFactorSphere(10.0*nanometer, 13.0*nanometer)
ff = FormFactorPyramid(13*nanometer, 10.0*nanometer, 60*degree)
# volume
# I suggest, that the form factor evaluated for q=0 gives the volume
# please correct if not
zero = cvector_t(0,0,0)
V = abs(ff.evaluate_for_q(zero))
# number of bins for qx, qy, qz
# I recommend to put at least 400 for a better picture quality
nqy = 400
nqz = 400
nqx = 400
# minimum and maximum values for qx, qy qz
qymin = -2.0
qymax = 2.0
qzmin = -2.0
qzmax = 2.0
qxmin = -2.0
qxmax = 2.0
# first and last bin for the projection slice
# if number of bins=400, then fbin should be set to 200 and lbin to 201
#fbin=200
#lbin=201
# step for qx, qy, qz
stepqy = (qymax - qymin)/(nqy-1)
stepqz = (qzmax - qzmin)/(nqz-1)
stepqx = (qxmax - qxmin)/(nqx-1)
# define style of the diagram
ROOT.gROOT.SetStyle("Plain")
ROOT.gStyle.SetOptStat(0);
ROOT.gStyle.SetOptTitle(0);
ROOT.gStyle.SetLabelSize(0.06, "xyz");
ROOT.gStyle.SetTitleSize(0.06, "xyz");
#ROOT.gStyle.SetPadRightMargin(0.1)
ROOT.gStyle.SetPadLeftMargin(0.18)
ROOT.gStyle.SetPadBottomMargin(0.18)
t = ROOT.TText()
# create ROOT histograms
hist = ROOT.TH2D("hist","Sphere:H=R",nqy,qymin,qymax, nqz, qzmin, qzmax)
hist2 = ROOT.TH2D("hist2","Sphere:H=R",nqx,qxmin,qxmax, nqy, qymin, qymax)
# and fill them with the values
for i in range(nqy):
qy = qymin + i*stepqy
for j in range(nqz):
qz = qzmin + j*stepqz
k = cvector_t(0,qy,qz)
hist.Fill(qy,qz,(abs(ff.evaluate_for_q(k)))**2+1)
for i in range(nqy):
qy = qymin + i*stepqy
for j in range(nqx):
qx = qxmin + j*stepqx
k = cvector_t(qx,qy,0)
hist2.Fill(qx,qy,(abs(ff.evaluate_for_q(k)))**2+1)
# create a ROOT canvas and put all plots on it
c = ROOT.TCanvas("FormfactorSphere","Formfactor Sphere", 1000,500)
c.Divide(2,1)
hist.SetMinimum(1)
hist.SetMaximum(V**2)
hist2.SetMinimum(1)
hist2.SetMaximum(V**2)
hist.SetContour(50)
hist2.SetContour(50)
#hist.GetZaxis().SetTitle(" |F|^{2} ")
#hist2.GetZaxis().SetTitle(" |F|^{2} ")
c.cd(1)
ROOT.gPad.SetLogz()
hist.GetXaxis().SetTitle(" q_{y} [nm^{-1}] ")
hist.GetXaxis().CenterTitle()
hist.GetXaxis().SetTitleOffset(1.1)
hist.GetYaxis().SetTitle(" q_{z} [nm^{-1}] ")
hist.GetYaxis().CenterTitle()
hist.GetYaxis().SetTitleOffset(1.4)
hist.GetZaxis().SetTitleOffset(1.3)
hist.Draw("col")
#t.SetNDC(1)
#t.SetTextSize(8.0e-2)
#t.DrawText(0.1, 0.01, "a")
c.cd(2)
ROOT.gPad.SetLogz()
hist2.SetMinimum(1)
hist2.GetXaxis().SetTitle(" q_{x} [nm^{-1}] ")
hist2.GetXaxis().CenterTitle()
hist2.GetXaxis().SetTitleOffset(1.1)
hist2.GetYaxis().SetTitle(" q_{y} [nm^{-1}] ")
hist2.GetYaxis().CenterTitle()
hist2.GetYaxis().SetTitleOffset(1.4)
hist2.GetZaxis().SetTitleOffset(1.3)
hist2.Draw("col")
#t.DrawText(0.1, 0.004, "b")
#c.cd(3)
#ROOT.gPad.SetLogy()
#py = hist.ProjectionX("py", fbin, lbin, 'o')
#py.GetYaxis().SetTitle(" |F|^{2} ")
#py.GetYaxis().SetTitleOffset(1.3)
#py.Draw()
#t.DrawText(0.1, 0.01, "c")
#c.cd(4)
#ROOT.gPad.SetLogy()
#px = hist2.ProjectionX("px", fbin, lbin, 'o')
#px.GetYaxis().SetTitle(" |F|^{2} ")
#px.GetYaxis().SetTitleOffset(1.3)
#px.Draw()
#t.DrawText(0.1, 0.004, "d")
c.Update()
ROOT.gApplication.Run()
devtools/math/ff/plot/plotpardep.py deleted 100644 → 0
+ 0
− 172
View file @ 4e62e705
# script to plot the form factor
# after the ROOT window with plots appear,
# just save the picture in that format which you need
import ROOT
from libBornAgainCore import *
# define a form factor, I recommend a 10--20 nm diameter
#ff = FormFactorFullSphere(10.0*nanometer)
ff05 = FormFactorSphere(10.0*nanometer, 5.0*nanometer)
ff10 = FormFactorSphere(10.0*nanometer, 10.0*nanometer)
ff15 = FormFactorSphere(10.0*nanometer, 15.0*nanometer)
# zero q vector
zero = cvector_t(0,0,0)
# volume of particles
V05 = abs(ff05.evaluate_for_q(zero))
V10 = abs(ff10.evaluate_for_q(zero))
V15 = abs(ff15.evaluate_for_q(zero))
# number of bins for qx, qy, qz
# I recommend to put at least 400 for a better picture quality
nqy = 100
nqz = 100
nqx = 100
# minimum and maximum values for qx, qy qz
qymin = -2.0
qymax = 2.0
qzmin = -2.0
qzmax = 2.0
qxmin = -2.0
qxmax = 2.0
# first and last bin for the projection slice
# if number of bins=400, then fbin should be set to 200 and lbin to 201
#fbin=50
#lbin=51
# step for qx, qy, qz
stepqy = (qymax - qymin)/(nqy-1)
stepqz = (qzmax - qzmin)/(nqz-1)
stepqx = (qxmax - qxmin)/(nqx-1)
# define style of the diagram
ROOT.gROOT.SetStyle("Plain")
ROOT.gStyle.SetOptStat(0);
ROOT.gStyle.SetOptTitle(1);
ROOT.gStyle.SetLabelSize(0.06, "xyz");
ROOT.gStyle.SetTitleSize(0.06, "xyz");
ROOT.gStyle.SetTitleFontSize(0.1);
#ROOT.gStyle.SetPadRightMargin(0.19)
ROOT.gStyle.SetPadLeftMargin(0.18)
ROOT.gStyle.SetPadBottomMargin(0.18)
ROOT.gStyle.SetTitleX(0.3);
#ROOT.gStyle.SetTitleW(0.5);
#t = ROOT.TText()
# create ROOT histograms
h05 = ROOT.TH2D("h05nm","H = 5 nm",nqy,qymin,qymax, nqz, qzmin, qzmax)
h10 = ROOT.TH2D("h10nm","H = 10 nm",nqy,qymin,qymax, nqz, qzmin, qzmax)
h15 = ROOT.TH2D("h15nm","H = 15 nm",nqy,qymin,qymax, nqz, qzmin, qzmax)
# and fill them with the values
for i in range(nqy):
qy = qymin + i*stepqy
for j in range(nqz):
qz = qzmin + j*stepqz
k = cvector_t(0,qy,qz)
h05.Fill(qy,qz,(abs(ff05.evaluate_for_q(k))/V05)**2)
h10.Fill(qy,qz,(abs(ff10.evaluate_for_q(k))/V10)**2)
h15.Fill(qy,qz,(abs(ff15.evaluate_for_q(k))/V15)**2)
# create a ROOT canvas and put all plots on it
c = ROOT.TCanvas("FormfactorSphere","Formfactor Sphere", 1200,400)
c.Divide(3,1)
h05.SetMinimum(5e-8)
#h05.SetMaximum(V05**2)
h05.SetMaximum(1)
h10.SetMinimum(5e-8)
#h10.SetMaximum(V10**2)
h10.SetMaximum(1)
h15.SetMinimum(5e-8)
#h15.SetMaximum(V15**2)
h15.SetMaximum(1)
h05.SetContour(50)
h10.SetContour(50)
h15.SetContour(50)
#h05.GetZaxis().SetTitle(" |F|^{2} ")
#h10.GetZaxis().SetTitle(" |F|^{2} ")
#h15.GetZaxis().SetTitle(" |F|^{2} ")
c.cd(1)
ROOT.gPad.SetLogz()
h05.GetXaxis().SetTitle(" q_{y} [nm^{-1}] ")
h05.GetXaxis().CenterTitle()
h05.GetXaxis().SetTitleOffset(1.1)
h05.GetYaxis().SetTitle(" q_{z} [nm^{-1}] ")
h05.GetYaxis().CenterTitle()
h05.GetYaxis().SetTitleOffset(1.4)
h05.GetZaxis().SetTitleOffset(1.3)
h05.Draw("col")
#t.SetNDC(1)
#t.SetTextSize(8.0e-2)
#t.DrawText(0.1, 0.01, "H = 5 nm")
c.cd(2)
ROOT.gPad.SetLogz()
h10.GetXaxis().SetTitle(" q_{y} [nm^{-1}] ")
h10.GetXaxis().CenterTitle()
h10.GetXaxis().SetTitleOffset(1.1)
h10.GetYaxis().SetTitle(" q_{z} [nm^{-1}] ")
h10.GetYaxis().CenterTitle()
h10.GetYaxis().SetTitleOffset(1.4)
h10.GetZaxis().SetTitleOffset(1.3)
h10.Draw("col")
#t.SetNDC(1)
#t.SetTextSize(8.0e-2)
#t.DrawText(0.1, 0.01, "H = 10 nm")
c.cd(3)
ROOT.gPad.SetLogz()
h15.GetXaxis().SetTitle(" q_{y} [nm^{-1}] ")
h15.GetXaxis().CenterTitle()
h15.GetXaxis().SetTitleOffset(1.1)
h15.GetYaxis().SetTitle(" q_{z} [nm^{-1}] ")
h15.GetYaxis().CenterTitle()
h15.GetYaxis().SetTitleOffset(1.4)
h15.GetZaxis().SetTitleOffset(1.3)
h15.Draw("col")
#t.SetNDC(1)
#t.SetTextSize(8.0e-2)
#t.DrawText(0.1, 0.01, "H = 15 nm")
# draw projections
#c.cd(4)
#ROOT.gPad.SetLogy()
#p05 = h05.ProjectionX("p05", fbin, lbin, 'o')
#p05.GetYaxis().SetTitle(" |F|^{2} ")
#p05.GetYaxis().SetTitleOffset(1.3)
#p05.Draw()
#t.DrawText(0.1, 0.01, "d")
#c.cd(5)
#ROOT.gPad.SetLogy()
#p10 = h10.ProjectionX("p10", fbin, lbin, 'o')
#p10.GetYaxis().SetTitle(" |F|^{2} ")
#p10.GetYaxis().SetTitleOffset(1.3)
#p10.Draw()
#t.DrawText(0.1, 0.01, "e")
#c.cd(6)
#ROOT.gPad.SetLogy()
#p15 = h15.ProjectionX("p15", fbin, lbin, 'o')
#p15.GetYaxis().SetTitle(" |F|^{2} ")
#p15.GetYaxis().SetTitleOffset(1.3)
#p15.Draw()
#t.DrawText(0.1, 0.01, "f")
c.Update()
ROOT.gApplication.Run()
devtools/math/ff/plot/plotpardeppyramid.py deleted 100644 → 0
+ 0
− 179
View file @ 4e62e705
# script to plot the form factor
# after the ROOT window with plots appear,
# just save the picture in that format which you need
import ROOT
from libBornAgainCore import *
# define a form factor, I recommend a 10--20 nm diameter
# vary H
ff05 = FormFactorPyramid(2.5*nanometer, 5.0*nanometer, 60*degree)
ff10 = FormFactorPyramid(5.0*nanometer, 5.0*nanometer, 60*degree)
ff15 = FormFactorPyramid(7.5*nanometer, 5.0*nanometer, 60*degree)
# vary alpha
#ff05 = FormFactorPyramid(5.0*nanometer, 5.0*nanometer, 50*degree)
#ff10 = FormFactorPyramid(5.0*nanometer, 5.0*nanometer, 65*degree)
#ff15 = FormFactorPyramid(5.0*nanometer, 5.0*nanometer, 80*degree)
# zero q vector
zero = cvector_t(0,0,0)
# volume of particles
V05 = abs(ff05.evaluate_for_q(zero))
V10 = abs(ff10.evaluate_for_q(zero))
V15 = abs(ff15.evaluate_for_q(zero))
# number of bins for qx, qy, qz
# I recommend to put at least 400 for a better picture quality
nqy = 400
nqz = 400
nqx = 400
# minimum and maximum values for qx, qy qz
qymin = -2.0
qymax = 2.0
qzmin = -2.0
qzmax = 2.0
qxmin = -2.0
qxmax = 2.0
# first and last bin for the projection slice
# if number of bins=400, then fbin should be set to 200 and lbin to 201
#fbin=50
#lbin=51
# step for qx, qy, qz
stepqy = (qymax - qymin)/(nqy-1)
stepqz = (qzmax - qzmin)/(nqz-1)
stepqx = (qxmax - qxmin)/(nqx-1)
# define style of the diagram
ROOT.gROOT.SetStyle("Plain")
ROOT.gStyle.SetOptStat(0);
ROOT.gStyle.SetOptTitle(1);
ROOT.gStyle.SetLabelSize(0.06, "xyz");
ROOT.gStyle.SetTitleSize(0.06, "xyz");
ROOT.gStyle.SetTitleFontSize(0.1);
#ROOT.gStyle.SetPadRightMargin(0.19)
ROOT.gStyle.SetPadLeftMargin(0.18)
ROOT.gStyle.SetPadBottomMargin(0.18)
ROOT.gStyle.SetTitleX(0.3);
#ROOT.gStyle.SetTitleW(0.5);
#t = ROOT.TText()
# create ROOT histograms
h05 = ROOT.TH2D("h05nm","H = 2.5 nm",nqy,qymin,qymax, nqz, qzmin, qzmax)
h10 = ROOT.TH2D("h10nm","H = 5.0 nm",nqy,qymin,qymax, nqz, qzmin, qzmax)
h15 = ROOT.TH2D("h15nm","H = 7.5 nm",nqy,qymin,qymax, nqz, qzmin, qzmax)
#h05 = ROOT.TH2D("a50"," #alpha = 50^{o} ",nqy,qymin,qymax, nqz, qzmin, qzmax)
#h10 = ROOT.TH2D("a65"," #alpha = 65^{o} ",nqy,qymin,qymax, nqz, qzmin, qzmax)
#h15 = ROOT.TH2D("a80"," #alpha = 80^{o} ",nqy,qymin,qymax, nqz, qzmin, qzmax)
# and fill them with the values
for i in range(nqy):
qy = qymin + i*stepqy
for j in range(nqz):
qz = qzmin + j*stepqz
k = cvector_t(0,qy,qz)
h05.Fill(qy,qz,(abs(ff05.evaluate_for_q(k))/V05)**2)
h10.Fill(qy,qz,(abs(ff10.evaluate_for_q(k))/V10)**2)
h15.Fill(qy,qz,(abs(ff15.evaluate_for_q(k))/V15)**2)
# create a ROOT canvas and put all plots on it
c = ROOT.TCanvas("FormfactorSphere","Formfactor Sphere", 1200,400)
c.Divide(3,1)
h05.SetMinimum(1e-7)
#h05.SetMaximum(V05**2)
h05.SetMaximum(1)
h10.SetMinimum(1e-7)
#h10.SetMaximum(V10**2)
h10.SetMaximum(1)
h15.SetMinimum(1e-7)
#h15.SetMaximum(V15**2)
h15.SetMaximum(1)
h05.SetContour(50)
h10.SetContour(50)
h15.SetContour(50)
#h05.GetZaxis().SetTitle(" |F|^{2} ")
#h10.GetZaxis().SetTitle(" |F|^{2} ")
#h15.GetZaxis().SetTitle(" |F|^{2} ")
c.cd(1)
ROOT.gPad.SetLogz()
h05.GetXaxis().SetTitle(" q_{y} [nm^{-1}] ")
h05.GetXaxis().CenterTitle()
h05.GetXaxis().SetTitleOffset(1.1)
h05.GetYaxis().SetTitle(" q_{z} [nm^{-1}] ")
h05.GetYaxis().CenterTitle()
h05.GetYaxis().SetTitleOffset(1.4)
h05.GetZaxis().SetTitleOffset(1.3)
h05.Draw("col")
#t.SetNDC(1)
#t.SetTextSize(8.0e-2)
#t.DrawText(0.1, 0.01, "H = 5 nm")
c.cd(2)
ROOT.gPad.SetLogz()
h10.GetXaxis().SetTitle(" q_{y} [nm^{-1}] ")
h10.GetXaxis().CenterTitle()
h10.GetXaxis().SetTitleOffset(1.1)
h10.GetYaxis().SetTitle(" q_{z} [nm^{-1}] ")
h10.GetYaxis().CenterTitle()
h10.GetYaxis().SetTitleOffset(1.4)
h10.GetZaxis().SetTitleOffset(1.3)
h10.Draw("col")
#t.SetNDC(1)
#t.SetTextSize(8.0e-2)
#t.DrawText(0.1, 0.01, "H = 10 nm")
c.cd(3)
ROOT.gPad.SetLogz()
h15.GetXaxis().SetTitle(" q_{y} [nm^{-1}] ")
h15.GetXaxis().CenterTitle()
h15.GetXaxis().SetTitleOffset(1.1)
h15.GetYaxis().SetTitle(" q_{z} [nm^{-1}] ")
h15.GetYaxis().CenterTitle()
h15.GetYaxis().SetTitleOffset(1.4)
h15.GetZaxis().SetTitleOffset(1.3)
h15.Draw("col")
#t.SetNDC(1)
#t.SetTextSize(8.0e-2)
#t.DrawText(0.1, 0.01, "H = 15 nm")
# draw projections
#c.cd(4)
#ROOT.gPad.SetLogy()
#p05 = h05.ProjectionX("p05", fbin, lbin, 'o')
#p05.GetYaxis().SetTitle(" |F|^{2} ")
#p05.GetYaxis().SetTitleOffset(1.3)
#p05.Draw()
#t.DrawText(0.1, 0.01, "d")
#c.cd(5)
#ROOT.gPad.SetLogy()
#p10 = h10.ProjectionX("p10", fbin, lbin, 'o')
#p10.GetYaxis().SetTitle(" |F|^{2} ")
#p10.GetYaxis().SetTitleOffset(1.3)
#p10.Draw()
#t.DrawText(0.1, 0.01, "e")
#c.cd(6)
#ROOT.gPad.SetLogy()
#p15 = h15.ProjectionX("p15", fbin, lbin, 'o')
#p15.GetYaxis().SetTitle(" |F|^{2} ")
#p15.GetYaxis().SetTitleOffset(1.3)
#p15.Draw()
#t.DrawText(0.1, 0.01, "f")
c.Update()
ROOT.gApplication.Run()
devtools/math/ff/polyhedron/setup_dodecahedron.py deleted 100755 → 0
+ 0
− 79
View file @ 4e62e705
#!/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) )
devtools/math/ff/polyhedron/setup_icosahedron.py deleted 100755 → 0
+ 0
− 67
View file @ 4e62e705
#!/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) )
devtools/math/ff/test/Makefile deleted 100644 → 0
+ 0
− 15
View file @ 4e62e705
CC = g++ -ggdb -O3 -std=c++11 -Wall -I/G/ba/Core/FormFactors/ -I/G/ba/Core/Geometry -I/G/ba/Core/Tools -I/G/ba/Core/Samples -I/usr/include/eigen3
# LIBS= -lm -lgslcblas -lgsl -lyaml-cpp
BIND = ./
all: $(BIND)runff
# source from BornAgain overwrites -lBornAgainCore, cf http://stackoverflow.com/questions/36667429
SRC = runff.cpp /G/ba/Core/FormFactors/FormFactorPolyhedron.cpp
DIAG = -D POLYHEDRAL_DIAGNOSTIC=ON
$(BIND)runff: $(SRC)
$(CC) $(DIAG) -L../../../build/lib/ -l:_libBornAgainCore.so -o $(BIND)runff $(SRC)
devtools/math/ff/test/ff_diff.py deleted 100755 → 0
+ 0
− 44
View file @ 4e62e705
#!/usr/bin/env python3
import re, sys
from math import sqrt
def read_file( fn ):
fd = open( fn, 'r' )
ti = fd.read().rstrip()
fd.close
ti = re.sub( r'^\s*#.*\n', '', ti, re.M )
return ti.split("\n")
## Main
maxdiff = 0
fpair = sys.argv[1:3]
tpair = [ read_file( fn ) for fn in fpair ]
if len(tpair[0])!=len(tpair[1]):
raise Exception( "tables differ in size" )
for i in range(len(tpair[0])):
dpair = [ re.split(r'\s+', t[i]) for t in tpair]
if len(dpair[0])!=len(dpair[1]):
raise Exception( "lines %i differ in size" % i )
for j in range(5):
if dpair[0][j]!=dpair[1][j]:
raise Exception( "lines %i differ in parameter %i" % (i,j) )
xypair = [ [ float(d[8]), float(d[9]) ] for d in dpair ]
out_of_range = False
for x,y in xypair:
if abs(x)>1e166 or abs(y)>1e166:
out_of_range = True
break
if out_of_range:
print( tpair[0][i] )
print( tpair[1][i] )
print( "result too large" )
continue
apair = [ sqrt(x**2+y**2) for x,y in xypair ]
aval = (apair[0] + apair[1])/2
diff = max( abs(xypair[0][0]-xypair[1][0]), abs(xypair[0][1]-xypair[1][1]) ) / aval
maxdiff = max( maxdiff, diff )
if diff>3.e-12:
print( " ".join(dpair[0]+dpair[1][8:]+["%g" % diff]) )
print( "# maxdiff = %g" % maxdiff )
devtools/math/ff/test/limit-run deleted 100644 → 0
+ 0
− 26
View file @ 4e62e705
# no q expansion, hence sum over faces, then no expansion in qpa
runff 2 7 4 1e-99 1e-99 0 > x0
# no q expansion, hence sum over faces, then expand in qpa
runff 2 7 4 1e-99 1e99 2 > a1
runff 2 7 4 1e-99 1e99 4 > a2
runff 2 7 4 1e-99 1e99 6 > a3
runff 2 7 4 1e-99 1e99 8 > a4
runff 2 7 4 1e-99 1e99 10 > a5
runff 2 7 4 1e-99 1e99 12 > a6
runff 2 7 4 1e-99 1e99 14 > a7
runff 2 7 4 1e-99 1e99 16 > a8
runff 2 7 4 1e-99 1e99 18 > a9
runff 2 7 4 1e-99 1e99 20 > as
# q expansion (no need for qpa expansion)
runff 2 7 4 1e99 1e-99 2 > h1
runff 2 7 4 1e99 1e-99 4 > h2
runff 2 7 4 1e99 1e-99 6 > h3
runff 2 7 4 1e99 1e-99 8 > h4
runff 2 7 4 1e99 1e-99 10 > h5
runff 2 7 4 1e99 1e-99 12 > h6
runff 2 7 4 1e99 1e-99 14 > h7
runff 2 7 4 1e99 1e-99 16 > h8
runff 2 7 4 1e99 1e-99 18 > h9
runff 2 7 4 1e99 1e-99 20 > hs
devtools/math/ff/test/runff.cpp deleted 100644 → 0
+ 0
− 403
View file @ 4e62e705
// Form factor computation: test program
// JWu 2016
// Reference:
// Joachim Wuttke
// "Form factor (Fourier shape transform) of polygon and polyhedron."
// To use these tests, FormFactorPolyhedron.cpp must be compiled
// with flag -D POLYHEDRAL_DIAGNOSTIC=ON
#include <cassert>
#include <iostream>
#include <iomanip>
#include <complex>
#include <vector>
#include "ParticleShapes.h"
#include "FormFactorTriangle.h"
using std::cout;
using std::cerr;
using std::vector;
static complex_t I(0.,1.);
static double eps(2e-16);
Diagnosis diagnosis;
int nshape = 13;
extern int n_limit;
extern double q_limit_series;
extern double qpa_limit_series;
//! Returns a pointer to a particle, according to given code
IFormFactorBorn* make_particle( int ishape )
{
if ( ishape==1 ) {
return new FormFactorDodecahedron(3.);
} else if( ishape==2 ) {
return new FormFactorIcosahedron(15.);
} else if( ishape==3 ) { // true tetrahedron
double alpha = 72 * Units::degree;
return new FormFactorTetrahedron(1., tan(alpha)/2/sqrt(3), alpha);
} else if( ishape==4 ) { // tetrahedral frustum
double alpha = 72 * Units::degree;
return new FormFactorTetrahedron(1., 0.5*tan(alpha)/2/sqrt(3), alpha);
} else if( ishape==5 ) { // tetrahedral frustum, flat one
double alpha = 80 * Units::degree;
return new FormFactorTetrahedron(1., 0.1*tan(alpha)/2/sqrt(3), alpha);
} else if( ishape==6 ) { // tetrahedral frustum as in BasicTest
return new FormFactorTetrahedron(16., 4., .8);
} else if( ishape==7 ) {
double alpha = 72 * Units::degree;
return new FormFactorCone6(10., 10., alpha);
} else if( ishape==8 ) {
return new FormFactorPyramid(1.5, .24, 1.);
} else if( ishape==9 ) {
return new FormFactorAnisoPyramid(2, .4, .24, 1.);
} else if( ishape==10) {
return new FormFactorPrism3(1.2, 1.);
} else if( ishape==11) {
return new FormFactorPrism6(1., 1.);
} else if( ishape==12) {
return new FormFactorTruncatedCube(4., 1.);
} else if( ishape==13 ) {
double alpha = 73 * Units::degree;
return new FormFactorCuboctahedron(1., 1., .8, alpha);
} else if( ishape==90 ) {
return new FormFactorTriangle(1.);
} else
throw "Shape not implemented";
}
//! Print q in a form that can be easily pasted to the command line for further investigation
std::string nice_q( cvector_t q )
{
std::ostringstream ret;
double qmax = 0;
ret << std::setprecision(16);
for( int i=0; i<3; ++i )
qmax = std::max( qmax, q[i].real() );
for( int i=0; i<3; ++i )
ret << q[i].real()/qmax << " " << q[i].imag()/qmax << " ";
ret << qmax;
return ret.str();
}
//! Bisect between two q's to find possible discontinuities
void bisect(
const IFormFactorBorn* polyh, int ishape, double q0mag,
cvector_t qi, complex_t ri, const Diagnosis di,
cvector_t qf, complex_t rf, const Diagnosis df,
double& maxrelstep )
{
assert( di!=df );
if( (qi-qf).mag()<4*eps*q0mag ) {
// narrowed down to minimal step, now check for continuity
double aval = (std::abs(ri) + std::abs(rf))/2;
double step = std::abs(ri-rf);
double relstep = step/aval;
if( relstep>maxrelstep){
cout<<"ishape="<<ishape<<": relstep "<<std::setprecision(8)<<relstep<<"="<<step<<"/"<<std::setprecision(16)<<aval<<" for "<<di<<"->"<<df<<" at q between "<<nice_q(qi)<<" and "<<nice_q(qf)<<"\n";
maxrelstep = relstep;
}
return;
}
cvector_t q2 = (qi + qf)/2.;
complex_t r2 = polyh->evaluate_for_q(q2);
Diagnosis d2 = diagnosis;
if( d2!=di )
bisect( polyh, ishape, q0mag, qi, ri, di, q2, r2, d2, maxrelstep );
if( d2!=df )
bisect( polyh, ishape, q0mag, q2, r2, d2, qf, rf, df, maxrelstep );
}
//! Computes form factor, and prints result according to outfilter.
void run( const IFormFactorBorn* polyh, int ishape, cvector_t q, int outfilter )
{
complex_t ret = polyh->evaluate_for_q(q);
cout<<std::scientific<<std::setprecision(16)<<std::setfill('0');
if ( outfilter==0 )
cout<<q.mag()<<" "<<std::abs(ret)<<" "<<ret.real()<<" "<<ret.imag()<<" "<<
diagnosis.nExpandedFaces<<std::noshowpos<<" "<<diagnosis.maxOrder;
else if( outfilter==1 )
cout<<ret.real();
else if( outfilter==2 )
cout<<ret.imag();
cout<<"\n";
}
//! Compute a form factor with modified control parameter settings
complex_t ff_modified( cvector_t q, const IFormFactorBorn* polyh, bool expand_qpa, bool expand_q )
{
PolyhedralFace::setLimits( expand_qpa ? 1e99 : 1e-99, 80 );
FormFactorPolyhedron::setLimits( expand_q ? 1e99 : 1e-99, 80 );
return polyh->evaluate_for_q( q );
}
void test_matching( int ishape, const vector<vector<cvector_t>>& scans )
{
cout<<ishape<<"\n";
cerr<<"shape "<<ishape<<" ...\n";
IFormFactorBorn* polyh( make_particle( ishape ) );
int n_mag = 25;
double mag_i = 1e-3;
double mag_f = 1e1;
for( int i=1; i<n_mag; ++i ) {
double mag = mag_i*pow(mag_f/mag_i,i/(n_mag-1.));
double res = 0;
for( const vector<cvector_t>& q_scan: scans ) {
assert( q_scan.size()== 1 );
cvector_t uq = q_scan[0];
const cvector_t q = mag * uq;
complex_t ff[2];
ff[0] = ff_modified( q, polyh, false, false );
ff[1] = ff_modified( q, polyh, true, false );
double dev = std::abs(ff[0]-ff[1])*2/(std::abs(ff[0])+std::abs(ff[1]));
res = std::max(res, dev );
if( 0 && dev>.1 )
cerr<<ishape<<" "<<mag<<" "<<std::setprecision(16)<<
dev<<" "<<ff[0]<<" "<<ff[1]<<" @ "<<q<<"\n";
}
cout<<" "<<mag*polyh->getRadius()<<" "<<std::setprecision(8)<<res<<"\n";
}
cout<<"\n";
}
double test_continuity( int ishape, const vector<vector<cvector_t>>& scans )
{
IFormFactorBorn* polyh( make_particle( ishape ) );
double maxrelstep = 0;
for( size_t j=0; j<scans.size(); ++j ) {
if( !(j%100) )
fprintf( stderr, "%2i: %8li/%li\r", ishape, j, scans.size() );
const vector<cvector_t>& q_scan = scans[j];
for( size_t i=1; i<q_scan.size(); ++i ) {
complex_t last_ret = polyh->evaluate_for_q(q_scan[i-1]);
Diagnosis last_diag = diagnosis;
complex_t ret = polyh->evaluate_for_q(q_scan[i]);
Diagnosis diag = diagnosis;
if( diag!=last_diag )
bisect( polyh, ishape, q_scan[i].mag(),
q_scan[i-1], last_ret, last_diag,
q_scan[i], ret, diag,
maxrelstep );
}
}
fprintf( stderr, "\n" );
cout<<"shape "<<ishape<<" => max rel step = "<<maxrelstep<<"\n";
return maxrelstep;
}
void loop_one_shape( int ishape, const vector<vector<cvector_t>>& scans )
{
IFormFactorBorn* polyh( make_particle( ishape ) );
for( size_t j=0; j<scans.size(); ++j ) {
if( !(j%100) )
fprintf( stderr, "%2i: %8li/%li\r", ishape, j, scans.size() );
const vector<cvector_t>& q_scan = scans[j];
for( const cvector_t q: q_scan ) {
complex_t ret = polyh->evaluate_for_q(q);
cout <<
ishape<<" " <<
std::setprecision(16) <<
q[0].real()<<" "<<q[0].imag()<<" " <<
q[1].real()<<" "<<q[1].imag()<<" " <<
q[2].real()<<" "<<q[2].imag()<<" " <<
std::setprecision(17) <<
q.mag()<<" " <<
ret.real()<<" "<<ret.imag()<<" "<<std::abs(ret)<<" " <<
diagnosis.nExpandedFaces<<" "<<diagnosis.maxOrder<<"\n";
}
}
}
vector<vector<cvector_t>> create_scans( int mode )
{
static int n_dir = 7;
cvector_t dir[n_dir] = {
{ 1., 0., 0. },
{ 0., 1., 0. },
{ 0., 0., 1. },
{ 0., 1., 1. },
{ 1., 0., 1. },
{ 1., 1., 0. },
{ 1., 1., 1. }
};
static int n_absdir = 8; // absorption component
cvector_t absdir[n_absdir] = {
{ 0., 0., 0. },
{ 1e-15*I, 0., 0. },
{ 0., 1e-10*I, 0. },
{ 0., 0., 1e-5*I },
{ .1*I, 0., 0. },
{ 0., .1*I, 0. },
{ 0., 0., .1*I },
{ .1*I, .1*I, .1*I }
};
vector<double> mag;
if ( mode==0 ) {
mag.resize(1);
mag[0] = 1;
} else if ( mode==1 ) {
mag = vector<double>( { 0, 1e-12, 1e-10, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3,
.01, .06, .2, .5, 1, 2, 5, 10, 20, 50, 100 } );
} else if (mode==2 ) {
mag.resize(1001);
mag[0] = 0.;
for( size_t i=1; i<mag.size(); ++i )
mag[i] = 1e-10*pow(1e13,(i-1.)/(mag.size()-2));
}
static const int n_difmag = 17;
static double difmag[n_difmag] =
{ 0., 3e-16, 1e-15, 3e-15, 1e-14, 1e-13, 1e-12, 1e-11, 1e-9, 1e-7,
1e-5, 1e-3, .01, .03, .1, .2, .3 };
vector<vector<cvector_t>> scans;
for( int i=0; i<n_dir; ++i ) {
for( int j=0; j<n_dir; ++j ) {
if( i==j )
continue;
for( int k=0; k<n_difmag; ++k ) {
for( int l=0; l<n_absdir; ++l ) {
cvector_t uq = ( dir[i].unit() + difmag[k]*dir[j] + absdir[l] ).unit();
vector<cvector_t> scan;
for( size_t m=0; m<mag.size(); ++m )
scan.push_back( mag[m] * uq );
scans.push_back( scan );
}
}
}
}
return scans;
}
void help_and_exit()
{
cerr<<"fftest limits inmode debmsg outfilter shape qxr qxi qyr qyi qzr qzi [q]\n";
cerr<<" limits: \"def\" | qlim qpalim nlim\n";
cerr<<" inmode: get q from 0:stdin | 1:cmdline | 2:loop\n";
cerr<<" debmsg: 0..6 (only if inmode=1)\n";
cerr<<" outfilter: return 0:all | 1:real | 2:imag\n";
cerr<<"fftest limits loop shape\n";
cerr<<" limits: \"def\" | qlim qpalim nlim\n";
cerr<<" loop: 3[cont_test] | 4[matching_test] | 5[for_plot]\n";
cerr<<" shape: 0[all] | 1..[specific]\n";
exit(0);
}
#define NEXTARG if( (++arg)-argv >= argc ) help_and_exit()
#define MULTIARG(n) if( (arg+=n)-argv >= argc ) help_and_exit()
int main (int argc, const char *argv[])
{
try {
diagnosis.debmsg = 0;
diagnosis.request_convergence = false;
const char** arg = argv;
NEXTARG;
if ( !strncmp( *arg, "def", 3 ) ) {
// do nothing
} else {
double q_limit = atof( *arg );
NEXTARG;
double qpa_limit = atof( *arg );
NEXTARG;
int n_limit = atoi( *arg );
FormFactorPolyhedron::setLimits( q_limit, n_limit );
PolyhedralFace::setLimits( qpa_limit, n_limit );
}
NEXTARG;
int inmode = atoi( *arg );
if( inmode==1 ) {
NEXTARG;
diagnosis.debmsg = atoi( *arg );
}
if( inmode<=2 ) {
NEXTARG;
int outfilter = atoi( *arg );
NEXTARG;
int ishape = atoi( *arg );
IFormFactorBorn* P( make_particle( ishape ) );
MULTIARG(6);
cvector_t uq(
complex_t( atof(*(arg-5)), atof(*(arg-4)) ),
complex_t( atof(*(arg-3)), atof(*(arg-2)) ),
complex_t( atof(*(arg-1)), atof(*(arg-0)) ) );
double mag;
if( inmode==0 ) {
// for use through Frida
int nop;
std::cin >> nop;
while( std::cin >> mag )
run( P, ishape, mag*uq, outfilter );
} else if( inmode==1 ) {
NEXTARG;
mag = atof( *arg );
run( P, ishape, mag*uq, outfilter );
} else if( inmode==2 ) {
int n_mag = 2001;
double mag_i = 1e-20;
double mag_f = 1e4;
for( int i=1; i<n_mag; ++i ) {
//mag = 180.*i/(n_mag-1);
mag = mag_i*pow(mag_f/mag_i,i/(n_mag-1.));
run( P, ishape, mag*uq, outfilter );
}
}
exit(0);
}
// it's some loop
NEXTARG;
int ishapepar = atoi( *arg );
if( inmode==3 ) { // continuity test
diagnosis.request_convergence = true;
vector<vector<cvector_t>> scans = create_scans( 1 );
double totmaxrelstep = 0;
if( ishapepar==0 ) {
for( int ishape=1; ishape<=nshape; ++ishape )
totmaxrelstep = std::max( totmaxrelstep, test_continuity( ishape, scans ) );
cout<<"grand total max rel step = "<<totmaxrelstep<<"\n";
} else {
test_continuity( ishapepar, scans );
}
exit (0);
}
if( inmode==4 ) { // compare computation methods
vector<vector<cvector_t>> scans = create_scans( 0 );
if( ishapepar==0 ) {
for( int ishape=1; ishape<=nshape; ++ishape )
test_matching( ishape, scans );
} else
test_matching( ishapepar, scans );
exit (0);
}
// it's a straight loop over q
vector<vector<cvector_t>> scans = create_scans( 2 );
if( ishapepar==0 ) {
for( int ishape=1; ishape<=nshape; ++ishape )
loop_one_shape( ishape, scans );
} else
loop_one_shape( ishapepar, scans );
exit(0);
} catch( const char* ex ) {
cerr<<"F(q) failed: "<<ex<<"\n";
exit(0);
}
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment