From cd28c04710106e59f8a9003a28ac9aa3b869b3e5 Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (o)" <j.wuttke@fz-juelich.de>
Date: Tue, 14 Dec 2021 14:32:47 +0100
Subject: [PATCH] ff devtools moved to libformfactor
---
devtools/math/ff/plot/plotff.py | 138 ------
devtools/math/ff/plot/plotpardep.py | 172 --------
devtools/math/ff/plot/plotpardeppyramid.py | 179 --------
.../math/ff/polyhedron/setup_dodecahedron.py | 79 ----
.../math/ff/polyhedron/setup_icosahedron.py | 67 ---
devtools/math/ff/test/Makefile | 15 -
devtools/math/ff/test/ff_diff.py | 44 --
devtools/math/ff/test/limit-run | 26 --
devtools/math/ff/test/runff.cpp | 403 ------------------
9 files changed, 1123 deletions(-)
delete mode 100644 devtools/math/ff/plot/plotff.py
delete mode 100644 devtools/math/ff/plot/plotpardep.py
delete mode 100644 devtools/math/ff/plot/plotpardeppyramid.py
delete mode 100755 devtools/math/ff/polyhedron/setup_dodecahedron.py
delete mode 100755 devtools/math/ff/polyhedron/setup_icosahedron.py
delete mode 100644 devtools/math/ff/test/Makefile
delete mode 100755 devtools/math/ff/test/ff_diff.py
delete mode 100644 devtools/math/ff/test/limit-run
delete mode 100644 devtools/math/ff/test/runff.cpp
diff --git a/devtools/math/ff/plot/plotff.py b/devtools/math/ff/plot/plotff.py
deleted file mode 100644
index bb50b79f06f..00000000000
--- a/devtools/math/ff/plot/plotff.py
+++ /dev/null
@@ -1,138 +0,0 @@
-# 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()
-
-
-
diff --git a/devtools/math/ff/plot/plotpardep.py b/devtools/math/ff/plot/plotpardep.py
deleted file mode 100644
index 784864bda20..00000000000
--- a/devtools/math/ff/plot/plotpardep.py
+++ /dev/null
@@ -1,172 +0,0 @@
-# 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()
-
-
-
diff --git a/devtools/math/ff/plot/plotpardeppyramid.py b/devtools/math/ff/plot/plotpardeppyramid.py
deleted file mode 100644
index bcc91437e47..00000000000
--- a/devtools/math/ff/plot/plotpardeppyramid.py
+++ /dev/null
@@ -1,179 +0,0 @@
-# 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()
-
-
-
diff --git a/devtools/math/ff/polyhedron/setup_dodecahedron.py b/devtools/math/ff/polyhedron/setup_dodecahedron.py
deleted file mode 100755
index ff5b2710bf5..00000000000
--- a/devtools/math/ff/polyhedron/setup_dodecahedron.py
+++ /dev/null
@@ -1,79 +0,0 @@
-#!/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/devtools/math/ff/polyhedron/setup_icosahedron.py b/devtools/math/ff/polyhedron/setup_icosahedron.py
deleted file mode 100755
index 77c552279be..00000000000
--- a/devtools/math/ff/polyhedron/setup_icosahedron.py
+++ /dev/null
@@ -1,67 +0,0 @@
-#!/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) )
diff --git a/devtools/math/ff/test/Makefile b/devtools/math/ff/test/Makefile
deleted file mode 100644
index 16b5785ea9e..00000000000
--- a/devtools/math/ff/test/Makefile
+++ /dev/null
@@ -1,15 +0,0 @@
-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)
diff --git a/devtools/math/ff/test/ff_diff.py b/devtools/math/ff/test/ff_diff.py
deleted file mode 100755
index 72bd9d4da73..00000000000
--- a/devtools/math/ff/test/ff_diff.py
+++ /dev/null
@@ -1,44 +0,0 @@
-#!/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 )
diff --git a/devtools/math/ff/test/limit-run b/devtools/math/ff/test/limit-run
deleted file mode 100644
index ce9f0cc5a9d..00000000000
--- a/devtools/math/ff/test/limit-run
+++ /dev/null
@@ -1,26 +0,0 @@
-# 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
diff --git a/devtools/math/ff/test/runff.cpp b/devtools/math/ff/test/runff.cpp
deleted file mode 100644
index 03a0ac8a620..00000000000
--- a/devtools/math/ff/test/runff.cpp
+++ /dev/null
@@ -1,403 +0,0 @@
-// 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);
- }
-}
--
GitLab