From cbfba2294c10e3e967c7b664052620b16edbf556 Mon Sep 17 00:00:00 2001
From: Tobias WEBER <tweber@ill.fr>
Date: Mon, 11 Oct 2021 16:14:50 +0200
Subject: [PATCH] moved repository to new location
---
AUTHORS | 1 -
LICENSE | 339 ------------------
Makefile | 161 ---------
NOTES | 16 -
README.md | 18 +-
ext/setup_externals.sh | 50 ---
ext/tlibs2-mag/AUTHORS | 1 -
ext/tlibs2-mag/LICENSE | 674 ------------------------------------
ext/tlibs2-mag/mag.h | 199 -----------
src/calc/drawskx.cpp | 78 -----
src/calc/drawskx.gpl | 58 ----
src/calc/dyn.cpp | 300 ----------------
src/calc/heliphase.cpp | 30 --
src/calc/skx_gs.cpp | 185 ----------
src/calc/weight.cpp | 165 ---------
src/calc/weight_sum.cpp | 226 ------------
src/core/constants.h | 96 -----
src/core/fp.cpp | 181 ----------
src/core/fp.h | 73 ----
src/core/heli.cpp | 399 ---------------------
src/core/heli.h | 111 ------
src/core/helper.h | 353 -------------------
src/core/lapacke_switch.h | 19 -
src/core/longfluct.cpp | 61 ----
src/core/longfluct.h | 49 ---
src/core/magsys.cpp | 315 -----------------
src/core/magsys.h | 74 ----
src/core/skx.cpp | 615 --------------------------------
src/core/skx.h | 127 -------
src/core/skx_default_gs.cxx | 62 ----
src/takin/convert.cpp | 163 ---------
src/takin/dump.cpp | 48 ---
src/takin/genskx.cpp | 154 --------
src/takin/merge.cpp | 104 ------
src/takin/skxgrid.py | 259 --------------
src/takin/takin.cpp | 439 -----------------------
src/takin/takin_grid.cpp | 492 --------------------------
37 files changed, 1 insertion(+), 6694 deletions(-)
delete mode 100644 AUTHORS
delete mode 100644 LICENSE
delete mode 100644 Makefile
delete mode 100644 NOTES
delete mode 100755 ext/setup_externals.sh
delete mode 100644 ext/tlibs2-mag/AUTHORS
delete mode 100644 ext/tlibs2-mag/LICENSE
delete mode 100644 ext/tlibs2-mag/mag.h
delete mode 100644 src/calc/drawskx.cpp
delete mode 100644 src/calc/drawskx.gpl
delete mode 100644 src/calc/dyn.cpp
delete mode 100644 src/calc/heliphase.cpp
delete mode 100644 src/calc/skx_gs.cpp
delete mode 100644 src/calc/weight.cpp
delete mode 100644 src/calc/weight_sum.cpp
delete mode 100644 src/core/constants.h
delete mode 100644 src/core/fp.cpp
delete mode 100644 src/core/fp.h
delete mode 100644 src/core/heli.cpp
delete mode 100644 src/core/heli.h
delete mode 100644 src/core/helper.h
delete mode 100644 src/core/lapacke_switch.h
delete mode 100644 src/core/longfluct.cpp
delete mode 100644 src/core/longfluct.h
delete mode 100644 src/core/magsys.cpp
delete mode 100644 src/core/magsys.h
delete mode 100644 src/core/skx.cpp
delete mode 100644 src/core/skx.h
delete mode 100644 src/core/skx_default_gs.cxx
delete mode 100644 src/takin/convert.cpp
delete mode 100644 src/takin/dump.cpp
delete mode 100644 src/takin/genskx.cpp
delete mode 100644 src/takin/merge.cpp
delete mode 100644 src/takin/skxgrid.py
delete mode 100644 src/takin/takin.cpp
delete mode 100644 src/takin/takin_grid.cpp
diff --git a/AUTHORS b/AUTHORS
deleted file mode 100644
index 50b29e3..0000000
--- a/AUTHORS
+++ /dev/null
@@ -1 +0,0 @@
-Tobias Weber <tweber@ill.fr>
diff --git a/LICENSE b/LICENSE
deleted file mode 100644
index d159169..0000000
--- a/LICENSE
+++ /dev/null
@@ -1,339 +0,0 @@
- GNU GENERAL PUBLIC LICENSE
- Version 2, June 1991
-
- Copyright (C) 1989, 1991 Free Software Foundation, Inc.,
- 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
- Everyone is permitted to copy and distribute verbatim copies
- of this license document, but changing it is not allowed.
-
- Preamble
-
- The licenses for most software are designed to take away your
-freedom to share and change it. By contrast, the GNU General Public
-License is intended to guarantee your freedom to share and change free
-software--to make sure the software is free for all its users. This
-General Public License applies to most of the Free Software
-Foundation's software and to any other program whose authors commit to
-using it. (Some other Free Software Foundation software is covered by
-the GNU Lesser General Public License instead.) You can apply it to
-your programs, too.
-
- When we speak of free software, we are referring to freedom, not
-price. Our General Public Licenses are designed to make sure that you
-have the freedom to distribute copies of free software (and charge for
-this service if you wish), that you receive source code or can get it
-if you want it, that you can change the software or use pieces of it
-in new free programs; and that you know you can do these things.
-
- To protect your rights, we need to make restrictions that forbid
-anyone to deny you these rights or to ask you to surrender the rights.
-These restrictions translate to certain responsibilities for you if you
-distribute copies of the software, or if you modify it.
-
- For example, if you distribute copies of such a program, whether
-gratis or for a fee, you must give the recipients all the rights that
-you have. You must make sure that they, too, receive or can get the
-source code. And you must show them these terms so they know their
-rights.
-
- We protect your rights with two steps: (1) copyright the software, and
-(2) offer you this license which gives you legal permission to copy,
-distribute and/or modify the software.
-
- Also, for each author's protection and ours, we want to make certain
-that everyone understands that there is no warranty for this free
-software. If the software is modified by someone else and passed on, we
-want its recipients to know that what they have is not the original, so
-that any problems introduced by others will not reflect on the original
-authors' reputations.
-
- Finally, any free program is threatened constantly by software
-patents. We wish to avoid the danger that redistributors of a free
-program will individually obtain patent licenses, in effect making the
-program proprietary. To prevent this, we have made it clear that any
-patent must be licensed for everyone's free use or not licensed at all.
-
- The precise terms and conditions for copying, distribution and
-modification follow.
-
- GNU GENERAL PUBLIC LICENSE
- TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
-
- 0. This License applies to any program or other work which contains
-a notice placed by the copyright holder saying it may be distributed
-under the terms of this General Public License. The "Program", below,
-refers to any such program or work, and a "work based on the Program"
-means either the Program or any derivative work under copyright law:
-that is to say, a work containing the Program or a portion of it,
-either verbatim or with modifications and/or translated into another
-language. (Hereinafter, translation is included without limitation in
-the term "modification".) Each licensee is addressed as "you".
-
-Activities other than copying, distribution and modification are not
-covered by this License; they are outside its scope. The act of
-running the Program is not restricted, and the output from the Program
-is covered only if its contents constitute a work based on the
-Program (independent of having been made by running the Program).
-Whether that is true depends on what the Program does.
-
- 1. You may copy and distribute verbatim copies of the Program's
-source code as you receive it, in any medium, provided that you
-conspicuously and appropriately publish on each copy an appropriate
-copyright notice and disclaimer of warranty; keep intact all the
-notices that refer to this License and to the absence of any warranty;
-and give any other recipients of the Program a copy of this License
-along with the Program.
-
-You may charge a fee for the physical act of transferring a copy, and
-you may at your option offer warranty protection in exchange for a fee.
-
- 2. You may modify your copy or copies of the Program or any portion
-of it, thus forming a work based on the Program, and copy and
-distribute such modifications or work under the terms of Section 1
-above, provided that you also meet all of these conditions:
-
- a) You must cause the modified files to carry prominent notices
- stating that you changed the files and the date of any change.
-
- b) You must cause any work that you distribute or publish, that in
- whole or in part contains or is derived from the Program or any
- part thereof, to be licensed as a whole at no charge to all third
- parties under the terms of this License.
-
- c) If the modified program normally reads commands interactively
- when run, you must cause it, when started running for such
- interactive use in the most ordinary way, to print or display an
- announcement including an appropriate copyright notice and a
- notice that there is no warranty (or else, saying that you provide
- a warranty) and that users may redistribute the program under
- these conditions, and telling the user how to view a copy of this
- License. (Exception: if the Program itself is interactive but
- does not normally print such an announcement, your work based on
- the Program is not required to print an announcement.)
-
-These requirements apply to the modified work as a whole. If
-identifiable sections of that work are not derived from the Program,
-and can be reasonably considered independent and separate works in
-themselves, then this License, and its terms, do not apply to those
-sections when you distribute them as separate works. But when you
-distribute the same sections as part of a whole which is a work based
-on the Program, the distribution of the whole must be on the terms of
-this License, whose permissions for other licensees extend to the
-entire whole, and thus to each and every part regardless of who wrote it.
-
-Thus, it is not the intent of this section to claim rights or contest
-your rights to work written entirely by you; rather, the intent is to
-exercise the right to control the distribution of derivative or
-collective works based on the Program.
-
-In addition, mere aggregation of another work not based on the Program
-with the Program (or with a work based on the Program) on a volume of
-a storage or distribution medium does not bring the other work under
-the scope of this License.
-
- 3. You may copy and distribute the Program (or a work based on it,
-under Section 2) in object code or executable form under the terms of
-Sections 1 and 2 above provided that you also do one of the following:
-
- a) Accompany it with the complete corresponding machine-readable
- source code, which must be distributed under the terms of Sections
- 1 and 2 above on a medium customarily used for software interchange; or,
-
- b) Accompany it with a written offer, valid for at least three
- years, to give any third party, for a charge no more than your
- cost of physically performing source distribution, a complete
- machine-readable copy of the corresponding source code, to be
- distributed under the terms of Sections 1 and 2 above on a medium
- customarily used for software interchange; or,
-
- c) Accompany it with the information you received as to the offer
- to distribute corresponding source code. (This alternative is
- allowed only for noncommercial distribution and only if you
- received the program in object code or executable form with such
- an offer, in accord with Subsection b above.)
-
-The source code for a work means the preferred form of the work for
-making modifications to it. For an executable work, complete source
-code means all the source code for all modules it contains, plus any
-associated interface definition files, plus the scripts used to
-control compilation and installation of the executable. However, as a
-special exception, the source code distributed need not include
-anything that is normally distributed (in either source or binary
-form) with the major components (compiler, kernel, and so on) of the
-operating system on which the executable runs, unless that component
-itself accompanies the executable.
-
-If distribution of executable or object code is made by offering
-access to copy from a designated place, then offering equivalent
-access to copy the source code from the same place counts as
-distribution of the source code, even though third parties are not
-compelled to copy the source along with the object code.
-
- 4. You may not copy, modify, sublicense, or distribute the Program
-except as expressly provided under this License. Any attempt
-otherwise to copy, modify, sublicense or distribute the Program is
-void, and will automatically terminate your rights under this License.
-However, parties who have received copies, or rights, from you under
-this License will not have their licenses terminated so long as such
-parties remain in full compliance.
-
- 5. You are not required to accept this License, since you have not
-signed it. However, nothing else grants you permission to modify or
-distribute the Program or its derivative works. These actions are
-prohibited by law if you do not accept this License. Therefore, by
-modifying or distributing the Program (or any work based on the
-Program), you indicate your acceptance of this License to do so, and
-all its terms and conditions for copying, distributing or modifying
-the Program or works based on it.
-
- 6. Each time you redistribute the Program (or any work based on the
-Program), the recipient automatically receives a license from the
-original licensor to copy, distribute or modify the Program subject to
-these terms and conditions. You may not impose any further
-restrictions on the recipients' exercise of the rights granted herein.
-You are not responsible for enforcing compliance by third parties to
-this License.
-
- 7. If, as a consequence of a court judgment or allegation of patent
-infringement or for any other reason (not limited to patent issues),
-conditions are imposed on you (whether by court order, agreement or
-otherwise) that contradict the conditions of this License, they do not
-excuse you from the conditions of this License. If you cannot
-distribute so as to satisfy simultaneously your obligations under this
-License and any other pertinent obligations, then as a consequence you
-may not distribute the Program at all. For example, if a patent
-license would not permit royalty-free redistribution of the Program by
-all those who receive copies directly or indirectly through you, then
-the only way you could satisfy both it and this License would be to
-refrain entirely from distribution of the Program.
-
-If any portion of this section is held invalid or unenforceable under
-any particular circumstance, the balance of the section is intended to
-apply and the section as a whole is intended to apply in other
-circumstances.
-
-It is not the purpose of this section to induce you to infringe any
-patents or other property right claims or to contest validity of any
-such claims; this section has the sole purpose of protecting the
-integrity of the free software distribution system, which is
-implemented by public license practices. Many people have made
-generous contributions to the wide range of software distributed
-through that system in reliance on consistent application of that
-system; it is up to the author/donor to decide if he or she is willing
-to distribute software through any other system and a licensee cannot
-impose that choice.
-
-This section is intended to make thoroughly clear what is believed to
-be a consequence of the rest of this License.
-
- 8. If the distribution and/or use of the Program is restricted in
-certain countries either by patents or by copyrighted interfaces, the
-original copyright holder who places the Program under this License
-may add an explicit geographical distribution limitation excluding
-those countries, so that distribution is permitted only in or among
-countries not thus excluded. In such case, this License incorporates
-the limitation as if written in the body of this License.
-
- 9. The Free Software Foundation may publish revised and/or new versions
-of the General Public License from time to time. Such new versions will
-be similar in spirit to the present version, but may differ in detail to
-address new problems or concerns.
-
-Each version is given a distinguishing version number. If the Program
-specifies a version number of this License which applies to it and "any
-later version", you have the option of following the terms and conditions
-either of that version or of any later version published by the Free
-Software Foundation. If the Program does not specify a version number of
-this License, you may choose any version ever published by the Free Software
-Foundation.
-
- 10. If you wish to incorporate parts of the Program into other free
-programs whose distribution conditions are different, write to the author
-to ask for permission. For software which is copyrighted by the Free
-Software Foundation, write to the Free Software Foundation; we sometimes
-make exceptions for this. Our decision will be guided by the two goals
-of preserving the free status of all derivatives of our free software and
-of promoting the sharing and reuse of software generally.
-
- NO WARRANTY
-
- 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
-FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN
-OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
-PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
-OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
-MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS
-TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE
-PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
-REPAIR OR CORRECTION.
-
- 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
-WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
-REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
-INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
-OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
-TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
-YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
-PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
-POSSIBILITY OF SUCH DAMAGES.
-
- END OF TERMS AND CONDITIONS
-
- How to Apply These Terms to Your New Programs
-
- If you develop a new program, and you want it to be of the greatest
-possible use to the public, the best way to achieve this is to make it
-free software which everyone can redistribute and change under these terms.
-
- To do so, attach the following notices to the program. It is safest
-to attach them to the start of each source file to most effectively
-convey the exclusion of warranty; and each file should have at least
-the "copyright" line and a pointer to where the full notice is found.
-
- <one line to give the program's name and a brief idea of what it does.>
- Copyright (C) <year> <name of author>
-
- This program is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
-
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License along
- with this program; if not, write to the Free Software Foundation, Inc.,
- 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-Also add information on how to contact you by electronic and paper mail.
-
-If the program is interactive, make it output a short notice like this
-when it starts in an interactive mode:
-
- Gnomovision version 69, Copyright (C) year name of author
- Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
- This is free software, and you are welcome to redistribute it
- under certain conditions; type `show c' for details.
-
-The hypothetical commands `show w' and `show c' should show the appropriate
-parts of the General Public License. Of course, the commands you use may
-be called something other than `show w' and `show c'; they could even be
-mouse-clicks or menu items--whatever suits your program.
-
-You should also get your employer (if you work as a programmer) or your
-school, if any, to sign a "copyright disclaimer" for the program, if
-necessary. Here is a sample; alter the names:
-
- Yoyodyne, Inc., hereby disclaims all copyright interest in the program
- `Gnomovision' (which makes passes at compilers) written by James Hacker.
-
- <signature of Ty Coon>, 1 April 1989
- Ty Coon, President of Vice
-
-This General Public License does not permit incorporating your program into
-proprietary programs. If your program is a subroutine library, you may
-consider it more useful to permit linking proprietary applications with the
-library. If this is what you want to do, use the GNU Lesser General
-Public License instead of this License.
diff --git a/Makefile b/Makefile
deleted file mode 100644
index 18a3aa8..0000000
--- a/Makefile
+++ /dev/null
@@ -1,161 +0,0 @@
-#
-# MnSi dynamics module for Takin and helper tools
-# @author Tobias Weber <tweber@ill.fr>
-# @date 2018-2020
-# @license GPLv2 (see 'LICENSE' file)
-#
-
-mingw_build = 0
-debug_build = 0
-strip_bins = 1
-
-
-# -----------------------------------------------------------------------------
-# setup
-# -----------------------------------------------------------------------------
-ifneq ($(mingw_build), 1)
- ifeq ("$(CXX)", "")
- CXX = g++
- endif
-
- SYSINCS = -I/usr/local/include \
- -I/usr/include/lapacke -I/usr/local/opt/lapack/include \
- -I/usr/include/qt5 -I/usr/include/x86_64-linux-gnu/qt5/ \
- #-I/home/tw/build/boost_1_73_0
- LIBDIRS = -L/usr/local/opt/lapack/lib -L/usr/local/lib
-
- LIBBOOSTSYS = -lboost_system
- LIBBOOSTFILESYS = -lboost_filesystem
-
- BIN_SUFFIX =
-else
- CXX = x86_64-w64-mingw32-g++
-
- SYSINCS = -I/usr/x86_64-w64-mingw32/sys-root/mingw/include \
- -I/usr/x86_64-w64-mingw32/sys-root/mingw/include/qt5
- LIBDIRS = -L/usr/x86_64-w64-mingw32/sys-root/mingw/bin/
-
- LIBBOOSTSYS = -lboost_system-x64
- LIBBOOSTFILESYS = -lboost_filesystem-x64
-
- BIN_SUFFIX = .exe
-endif
-
-
-ifneq ($(debug_build), 1)
- OPT = -O2 #-march=native
-
- ifeq ($(strip_bins), 1)
- STRIP = strip
- else
- STRIP = echo -e "Not stripping"
- endif
-else
- OPT = -g -ggdb
- STRIP = echo -e "Not stripping"
-endif
-
-
-STD = -std=c++17
-LIBDEFS = -fPIC
-DEFS = -DDEF_SKX_ORDER=7 -DDEF_HELI_ORDER=7 \
- -DNO_MINIMISATION -DNO_REDEFINITIONS \
- -D__HACK_FULL_INST__ #-DPLUGIN_APPLI
-INCS = -Isrc -Iext -Iext/takin $(SYSINCS)
-# -----------------------------------------------------------------------------
-
-
-# -----------------------------------------------------------------------------
-# meta rules
-# -----------------------------------------------------------------------------
-.PHONY: all clean
-
-all: prepare lib/skxmod.so lib/skxmod_grid.so \
- bin/genskx bin/merge bin/convert bin/dump \
- bin/drawskx bin/dyn bin/weight \
-# bin/heliphase bin/skx_gs bin/weight_sum
-
-clean:
- find . -name "*.o" -exec rm -fv {} \;
- rm -rfv bin/
- rm -rfv lib/
-
-prepare:
- mkdir -p bin/
- mkdir -p lib/
-# -----------------------------------------------------------------------------
-
-
-# -----------------------------------------------------------------------------
-# Takin plugin modules
-# -----------------------------------------------------------------------------
-lib/skxmod.so: src/takin/takin.o src/core/skx.o src/core/fp.o src/core/heli.o src/core/longfluct.o src/core/magsys.o \
- ext/takin/tools/monteconvo/sqwbase.o ext/tlibs2/libs/log.o
- @echo "Linking Takin module $@..."
- $(CXX) $(STD) $(OPT) $(DEFS) $(LIBDIRS) $(LIBDEFS) -shared -o $@ $+ -llapacke
- $(STRIP) $@
-
-lib/skxmod_grid.so: src/takin/takin_grid.o ext/takin/tools/monteconvo/sqwbase.o ext/tlibs2/libs/log.o
- @echo "Linking Takin grid module $@..."
- $(CXX) $(STD) $(OPT) $(DEFS) $(LIBDIRS) $(LIBDEFS) -shared -o $@ $+ $(LIBBOOSTSYS) -lQt5Core
- $(STRIP) $@
-# -----------------------------------------------------------------------------
-
-
-# -----------------------------------------------------------------------------
-# tools
-# -----------------------------------------------------------------------------
-bin/genskx: src/takin/genskx.o src/core/skx.o src/core/magsys.o ext/tlibs2/libs/log.o
- $(CXX) $(STD) $(OPT) $(DEFS) $(LIBDIRS) $(LIBDEFS) -o $@ $+ $(LIBBOOSTFILESYS) -llapacke -lpthread
- $(STRIP) $@$(BIN_SUFFIX)
-
-bin/merge: src/takin/merge.o
- $(CXX) $(STD) $(OPT) $(DEFS) $(LIBDIRS) -o $@ $+ $(LIBBOOSTSYS)
- $(STRIP) $@$(BIN_SUFFIX)
-
-bin/convert: src/takin/convert.o
- $(CXX) $(STD) $(OPT) $(DEFS) $(LIBDIRS) -o $@ $+
- $(STRIP) $@$(BIN_SUFFIX)
-
-bin/dump: src/takin/dump.o
- $(CXX) $(STD) $(OPT) $(DEFS) $(LIBDIRS) -o $@ $+
- $(STRIP) $@$(BIN_SUFFIX)
-
-bin/drawskx: src/calc/drawskx.o
- $(CXX) $(STD) $(OPT) $(DEFS) $(LIBDIRS) -o $@ $+
- $(STRIP) $@$(BIN_SUFFIX)
-
-bin/dyn: src/calc/dyn.o src/core/skx.o src/core/fp.o src/core/heli.o src/core/magsys.o ext/tlibs2/libs/log.o
- $(CXX) $(STD) $(OPT) $(DEFS) $(LIBDIRS) $(LIBDEFS) -o $@ $+ -llapacke -lpthread
- $(STRIP) $@$(BIN_SUFFIX)
-
-bin/weight: src/calc/weight.o src/core/skx.o src/core/fp.o src/core/heli.o src/core/magsys.o ext/tlibs2/libs/log.o
- $(CXX) $(STD) $(OPT) $(DEFS) $(LIBDIRS) $(LIBDEFS) -o $@ $+ -llapacke -lpthread
- $(STRIP) $@$(BIN_SUFFIX)
-
-bin/weight_sum: src/calc/weight_sum.o src/core/skx.o src/core/fp.o src/core/heli.o src/core/magsys.o ext/tlibs2/libs/log.o
- $(CXX) $(STD) $(OPT) $(DEFS) $(LIBDIRS) $(LIBDEFS) -o $@ $+ -llapacke -lpthread
- $(STRIP) $@$(BIN_SUFFIX)
-# -----------------------------------------------------------------------------
-
-
-# -----------------------------------------------------------------------------
-# further tools needing specialised compilation options
-# -----------------------------------------------------------------------------
-bin/heliphase: src/calc/heliphase.cpp src/core/heli.cpp src/core/magsys.cpp ext/tlibs2/libs/log.cpp
- $(CXX) $(STD) $(OPT) $(INCS) -DDEF_HELI_ORDER=4 -DNO_REDEFINITIONS -D__HACK_FULL_INST__ $(LIBDIRS) -o $@ $+ -lMinuit2 -llapacke -lgomp
- $(STRIP) $@$(BIN_SUFFIX)
-
-bin/skx_gs: src/calc/skx_gs.cpp src/core/skx.cpp src/core/heli.cpp src/core/magsys.cpp ext/tlibs2/libs/log.cpp
- $(CXX) $(STD) $(OPT) $(INCS) -DDEF_SKX_ORDER=7 -DDEF_HELI_ORDER=7 -DNO_REDEFINITIONS -D__HACK_FULL_INST__ $(LIBDIRS) -o $@ $+ -lMinuit2 -llapacke -lgomp
- $(STRIP) $@$(BIN_SUFFIX)
-# -----------------------------------------------------------------------------
-
-
-# -----------------------------------------------------------------------------
-# general rules
-# -----------------------------------------------------------------------------
-%.o: %.cpp
- @echo "Compiling $< -> $@..."
- $(CXX) $(STD) $(OPT) $(DEFS) $(INCS) $(LIBDEFS) -c $< -o $@
-# -----------------------------------------------------------------------------
diff --git a/NOTES b/NOTES
deleted file mode 100644
index 1388aa9..0000000
--- a/NOTES
+++ /dev/null
@@ -1,16 +0,0 @@
-Manual compilation of takin external plugins:
-
- x86_64-w64-mingw32-g++ -std=c++17 -O2 -Isrc -Iext -Iext/takin -I/usr/x86_64-w64-mingw32/sys-root/mingw/include -I/usr/x86_64-w64-mingw32/sys-root/mingw/include/qt5 -L/usr/x86_64-w64-mingw32/sys-root/mingw/bin/ -DDEF_SKX_ORDER=7 -DDEF_HELI_ORDER=7 -DNO_MINIMISATION -DNO_REDEFINITIONS -D__HACK_FULL_INST__ -DPLUGIN_APPLI -o takinmod_skx.exe src/takin/takin.cpp src/core/skx.cpp src/core/heli.cpp src/core/fp.cpp src/core/magsys.cpp ext/takin/tools/monteconvo/sqwbase.cpp ext/tlibs2/libs/log.cpp ext/tlibs/log/log.cpp ext/tlibs/math/rand.cpp -llapacke -lboost_system-x64 -lboost_filesystem-x64
-
- g++ -std=c++17 -O2 -Isrc -Iext -Iext/takin -I/usr/local/include -I/usr/include/lapacke -I/usr/local/opt/lapack/include -L/usr/x86_64-w64-mingw32/sys-root/mingw/bin/ -DDEF_SKX_ORDER=7 -DDEF_HELI_ORDER=7 -DNO_MINIMISATION -DNO_REDEFINITIONS -D__HACK_FULL_INST__ -DPLUGIN_APPLI -o takinmod_skx src/takin/takin.cpp src/core/skx.cpp src/core/heli.cpp src/core/fp.cpp src/core/magsys.cpp ext/takin/tools/monteconvo/sqwbase.cpp ext/tlibs2/libs/log.cpp ext/tlibs/log/log.cpp ext/tlibs/math/rand.cpp -llapacke -lboost_system -lboost_filesystem -lrt -lpthread
-
-
-Manual compilation of Takin plugins:
-
- x86_64-w64-mingw32-g++ -std=c++17 -O2 -shared -fPIC -Isrc -Iext -Iext/takin -I/usr/x86_64-w64-mingw32/sys-root/mingw/include -I/usr/x86_64-w64-mingw32/sys-root/mingw/include/qt5 -L/usr/x86_64-w64-mingw32/sys-root/mingw/bin/ -DDEF_SKX_ORDER=7 -DDEF_HELI_ORDER=7 -DNO_MINIMISATION -DNO_REDEFINITIONS -D__HACK_FULL_INST__ -o lib/skxmod.dll src/takin/takin.cpp src/core/skx.cpp src/core/heli.cpp src/core/fp.cpp src/core/magsys.cpp ext/takin/tools/monteconvo/sqwbase.cpp ext/tlibs2/libs/log.cpp ext/tlibs/log/log.cpp ext/tlibs/math/rand.cpp -llapacke -lboost_system-x64 -lboost_filesystem-x64
-
- x86_64-w64-mingw32-g++ -std=c++17 -O2 -shared -fPIC -DNO_REDEFINITIONS -Isrc -Iext -Iext/takin -I/usr/x86_64-w64-mingw32/sys-root/mingw/include -I/usr/x86_64-w64-mingw32/sys-root/mingw/include/qt5 -I/usr/include/qt5 -I/usr/include/x86_64-linux-gnu/qt5/ -L/usr/x86_64-w64-mingw32/sys-root/mingw/bin/ -o lib/skxmod_grid.dll src/takin/takin_grid.cpp ext/takin/tools/monteconvo/sqwbase.cpp ext/tlibs2/libs/log.cpp ext/tlibs/log/log.cpp ext/tlibs/math/rand.cpp -lboost_system-x64 -lboost_filesystem-x64 -lQt5Core
-
-
-Grid test:
- g++ -std=c++17 -DDO_TEST -fPIC -o grid_tst -Isrc -Iext -Iext/takin -I/usr/include/qt5 -I/usr/include/x86_64-linux-gnu/qt5/ -DNO_REDEFINITIONS src/takin/takin_grid.cpp ext/takin/tools/monteconvo/sqwbase.cpp ext/tlibs2/libs/log.cpp -lboost_system -lQt5Core
diff --git a/README.md b/README.md
index af7a312..f7563e0 100644
--- a/README.md
+++ b/README.md
@@ -1,17 +1 @@
-This is the code for the Takin MnSi plugin module and its helper tools.
-It calculates the dispersion and the dynamical structure factor for the
-helimagnetic, the field-polarised, and the skyrmion phase of MnSi.
-
-The ext/ directory contains the source code of the external libraries on
-which this code depends.
-
-The development repository can be found here:
-https://code.ill.fr/tweber/takin-mnsi
-
-This code has been used in the following papers:
- - https://doi.org/10.1103/PhysRevB.100.060404
- - https://doi.org/10.1103/PhysRevB.97.224403
- - https://doi.org/10.1063/1.5041036
-
-This repository has moved to
-[https://code.ill.fr/scientific-software/takin/plugins/mnsi](https://code.ill.fr/scientific-software/takin/plugins/mnsi).
+## This repository has moved to [https://code.ill.fr/scientific-software/takin/plugins/mnsi](https://code.ill.fr/scientific-software/takin/plugins/mnsi).
diff --git a/ext/setup_externals.sh b/ext/setup_externals.sh
deleted file mode 100755
index 7377f5b..0000000
--- a/ext/setup_externals.sh
+++ /dev/null
@@ -1,50 +0,0 @@
-#!/bin/bash
-#
-# downloads external dependencies
-# @author Tobias Weber <tweber@ill.fr>
-# @date 9-apr-20
-# @license GPLv2
-#
-
-if [ -L takin ] || [ -d takin ]; then
- echo -e "A takin directory already exists. Skipping.";
-else
- echo -e "Cloning takin/core..."
- git clone https://code.ill.fr/scientific-software/takin/core.git
- mv -v core takin
-fi
-
-
-if [ -L tlibs ] || [ -d tlibs ]; then
- echo -e "A tlibs directory already exists. Skipping.";
-else
- echo -e "Cloning tlibs..."
- git clone https://code.ill.fr/scientific-software/takin/tlibs.git
-fi
-
-
-if [ -L takin2 ] || [ -d takin2 ]; then
- echo -e "A takin2 directory already exists. Skipping.";
-else
- echo -e "Cloning takin/mag-core repo..."
- git clone https://code.ill.fr/scientific-software/takin/mag-core.git
- mv -v mag-core takin2
-fi
-
-
-if [ -L tlibs2 ] || [ -d tlibs2 ]; then
- echo -e "A tlibs2 directory already exists. Skipping.";
-else
- echo -e "Cloning tlibs2 repo..."
- git clone https://code.ill.fr/scientific-software/takin/tlibs2.git
-fi
-
-
-if [ -L tlibs2-mag ] || [ -d tlibs2-mag ]; then
- echo -e "A tlibs2-mag directory already exists. Skipping.";
-else
- echo -e "Cloning tlibs2-mag repo..."
- git clone https://code.ill.fr/tweber/tlibs2_magnon_helpers.git
- mv -v tlibs2_magnon_helpers tlibs2-mag
- cp -v tlibs2-mag/mag.h tlibs2/libs/
-fi
diff --git a/ext/tlibs2-mag/AUTHORS b/ext/tlibs2-mag/AUTHORS
deleted file mode 100644
index 50b29e3..0000000
--- a/ext/tlibs2-mag/AUTHORS
+++ /dev/null
@@ -1 +0,0 @@
-Tobias Weber <tweber@ill.fr>
diff --git a/ext/tlibs2-mag/LICENSE b/ext/tlibs2-mag/LICENSE
deleted file mode 100644
index f17f169..0000000
--- a/ext/tlibs2-mag/LICENSE
+++ /dev/null
@@ -1,674 +0,0 @@
- GNU GENERAL PUBLIC LICENSE
- Version 3, 29 June 2007
-
- Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
- Everyone is permitted to copy and distribute verbatim copies
- of this license document, but changing it is not allowed.
-
- Preamble
-
- The GNU General Public License is a free, copyleft license for
-software and other kinds of works.
-
- The licenses for most software and other practical works are designed
-to take away your freedom to share and change the works. By contrast,
-the GNU General Public License is intended to guarantee your freedom to
-share and change all versions of a program--to make sure it remains free
-software for all its users. We, the Free Software Foundation, use the
-GNU General Public License for most of our software; it applies also to
-any other work released this way by its authors. You can apply it to
-your programs, too.
-
- When we speak of free software, we are referring to freedom, not
-price. Our General Public Licenses are designed to make sure that you
-have the freedom to distribute copies of free software (and charge for
-them if you wish), that you receive source code or can get it if you
-want it, that you can change the software or use pieces of it in new
-free programs, and that you know you can do these things.
-
- To protect your rights, we need to prevent others from denying you
-these rights or asking you to surrender the rights. Therefore, you have
-certain responsibilities if you distribute copies of the software, or if
-you modify it: responsibilities to respect the freedom of others.
-
- For example, if you distribute copies of such a program, whether
-gratis or for a fee, you must pass on to the recipients the same
-freedoms that you received. You must make sure that they, too, receive
-or can get the source code. And you must show them these terms so they
-know their rights.
-
- Developers that use the GNU GPL protect your rights with two steps:
-(1) assert copyright on the software, and (2) offer you this License
-giving you legal permission to copy, distribute and/or modify it.
-
- For the developers' and authors' protection, the GPL clearly explains
-that there is no warranty for this free software. For both users' and
-authors' sake, the GPL requires that modified versions be marked as
-changed, so that their problems will not be attributed erroneously to
-authors of previous versions.
-
- Some devices are designed to deny users access to install or run
-modified versions of the software inside them, although the manufacturer
-can do so. This is fundamentally incompatible with the aim of
-protecting users' freedom to change the software. The systematic
-pattern of such abuse occurs in the area of products for individuals to
-use, which is precisely where it is most unacceptable. Therefore, we
-have designed this version of the GPL to prohibit the practice for those
-products. If such problems arise substantially in other domains, we
-stand ready to extend this provision to those domains in future versions
-of the GPL, as needed to protect the freedom of users.
-
- Finally, every program is threatened constantly by software patents.
-States should not allow patents to restrict development and use of
-software on general-purpose computers, but in those that do, we wish to
-avoid the special danger that patents applied to a free program could
-make it effectively proprietary. To prevent this, the GPL assures that
-patents cannot be used to render the program non-free.
-
- The precise terms and conditions for copying, distribution and
-modification follow.
-
- TERMS AND CONDITIONS
-
- 0. Definitions.
-
- "This License" refers to version 3 of the GNU General Public License.
-
- "Copyright" also means copyright-like laws that apply to other kinds of
-works, such as semiconductor masks.
-
- "The Program" refers to any copyrightable work licensed under this
-License. Each licensee is addressed as "you". "Licensees" and
-"recipients" may be individuals or organizations.
-
- To "modify" a work means to copy from or adapt all or part of the work
-in a fashion requiring copyright permission, other than the making of an
-exact copy. The resulting work is called a "modified version" of the
-earlier work or a work "based on" the earlier work.
-
- A "covered work" means either the unmodified Program or a work based
-on the Program.
-
- To "propagate" a work means to do anything with it that, without
-permission, would make you directly or secondarily liable for
-infringement under applicable copyright law, except executing it on a
-computer or modifying a private copy. Propagation includes copying,
-distribution (with or without modification), making available to the
-public, and in some countries other activities as well.
-
- To "convey" a work means any kind of propagation that enables other
-parties to make or receive copies. Mere interaction with a user through
-a computer network, with no transfer of a copy, is not conveying.
-
- An interactive user interface displays "Appropriate Legal Notices"
-to the extent that it includes a convenient and prominently visible
-feature that (1) displays an appropriate copyright notice, and (2)
-tells the user that there is no warranty for the work (except to the
-extent that warranties are provided), that licensees may convey the
-work under this License, and how to view a copy of this License. If
-the interface presents a list of user commands or options, such as a
-menu, a prominent item in the list meets this criterion.
-
- 1. Source Code.
-
- The "source code" for a work means the preferred form of the work
-for making modifications to it. "Object code" means any non-source
-form of a work.
-
- A "Standard Interface" means an interface that either is an official
-standard defined by a recognized standards body, or, in the case of
-interfaces specified for a particular programming language, one that
-is widely used among developers working in that language.
-
- The "System Libraries" of an executable work include anything, other
-than the work as a whole, that (a) is included in the normal form of
-packaging a Major Component, but which is not part of that Major
-Component, and (b) serves only to enable use of the work with that
-Major Component, or to implement a Standard Interface for which an
-implementation is available to the public in source code form. A
-"Major Component", in this context, means a major essential component
-(kernel, window system, and so on) of the specific operating system
-(if any) on which the executable work runs, or a compiler used to
-produce the work, or an object code interpreter used to run it.
-
- The "Corresponding Source" for a work in object code form means all
-the source code needed to generate, install, and (for an executable
-work) run the object code and to modify the work, including scripts to
-control those activities. However, it does not include the work's
-System Libraries, or general-purpose tools or generally available free
-programs which are used unmodified in performing those activities but
-which are not part of the work. For example, Corresponding Source
-includes interface definition files associated with source files for
-the work, and the source code for shared libraries and dynamically
-linked subprograms that the work is specifically designed to require,
-such as by intimate data communication or control flow between those
-subprograms and other parts of the work.
-
- The Corresponding Source need not include anything that users
-can regenerate automatically from other parts of the Corresponding
-Source.
-
- The Corresponding Source for a work in source code form is that
-same work.
-
- 2. Basic Permissions.
-
- All rights granted under this License are granted for the term of
-copyright on the Program, and are irrevocable provided the stated
-conditions are met. This License explicitly affirms your unlimited
-permission to run the unmodified Program. The output from running a
-covered work is covered by this License only if the output, given its
-content, constitutes a covered work. This License acknowledges your
-rights of fair use or other equivalent, as provided by copyright law.
-
- You may make, run and propagate covered works that you do not
-convey, without conditions so long as your license otherwise remains
-in force. You may convey covered works to others for the sole purpose
-of having them make modifications exclusively for you, or provide you
-with facilities for running those works, provided that you comply with
-the terms of this License in conveying all material for which you do
-not control copyright. Those thus making or running the covered works
-for you must do so exclusively on your behalf, under your direction
-and control, on terms that prohibit them from making any copies of
-your copyrighted material outside their relationship with you.
-
- Conveying under any other circumstances is permitted solely under
-the conditions stated below. Sublicensing is not allowed; section 10
-makes it unnecessary.
-
- 3. Protecting Users' Legal Rights From Anti-Circumvention Law.
-
- No covered work shall be deemed part of an effective technological
-measure under any applicable law fulfilling obligations under article
-11 of the WIPO copyright treaty adopted on 20 December 1996, or
-similar laws prohibiting or restricting circumvention of such
-measures.
-
- When you convey a covered work, you waive any legal power to forbid
-circumvention of technological measures to the extent such circumvention
-is effected by exercising rights under this License with respect to
-the covered work, and you disclaim any intention to limit operation or
-modification of the work as a means of enforcing, against the work's
-users, your or third parties' legal rights to forbid circumvention of
-technological measures.
-
- 4. Conveying Verbatim Copies.
-
- You may convey verbatim copies of the Program's source code as you
-receive it, in any medium, provided that you conspicuously and
-appropriately publish on each copy an appropriate copyright notice;
-keep intact all notices stating that this License and any
-non-permissive terms added in accord with section 7 apply to the code;
-keep intact all notices of the absence of any warranty; and give all
-recipients a copy of this License along with the Program.
-
- You may charge any price or no price for each copy that you convey,
-and you may offer support or warranty protection for a fee.
-
- 5. Conveying Modified Source Versions.
-
- You may convey a work based on the Program, or the modifications to
-produce it from the Program, in the form of source code under the
-terms of section 4, provided that you also meet all of these conditions:
-
- a) The work must carry prominent notices stating that you modified
- it, and giving a relevant date.
-
- b) The work must carry prominent notices stating that it is
- released under this License and any conditions added under section
- 7. This requirement modifies the requirement in section 4 to
- "keep intact all notices".
-
- c) You must license the entire work, as a whole, under this
- License to anyone who comes into possession of a copy. This
- License will therefore apply, along with any applicable section 7
- additional terms, to the whole of the work, and all its parts,
- regardless of how they are packaged. This License gives no
- permission to license the work in any other way, but it does not
- invalidate such permission if you have separately received it.
-
- d) If the work has interactive user interfaces, each must display
- Appropriate Legal Notices; however, if the Program has interactive
- interfaces that do not display Appropriate Legal Notices, your
- work need not make them do so.
-
- A compilation of a covered work with other separate and independent
-works, which are not by their nature extensions of the covered work,
-and which are not combined with it such as to form a larger program,
-in or on a volume of a storage or distribution medium, is called an
-"aggregate" if the compilation and its resulting copyright are not
-used to limit the access or legal rights of the compilation's users
-beyond what the individual works permit. Inclusion of a covered work
-in an aggregate does not cause this License to apply to the other
-parts of the aggregate.
-
- 6. Conveying Non-Source Forms.
-
- You may convey a covered work in object code form under the terms
-of sections 4 and 5, provided that you also convey the
-machine-readable Corresponding Source under the terms of this License,
-in one of these ways:
-
- a) Convey the object code in, or embodied in, a physical product
- (including a physical distribution medium), accompanied by the
- Corresponding Source fixed on a durable physical medium
- customarily used for software interchange.
-
- b) Convey the object code in, or embodied in, a physical product
- (including a physical distribution medium), accompanied by a
- written offer, valid for at least three years and valid for as
- long as you offer spare parts or customer support for that product
- model, to give anyone who possesses the object code either (1) a
- copy of the Corresponding Source for all the software in the
- product that is covered by this License, on a durable physical
- medium customarily used for software interchange, for a price no
- more than your reasonable cost of physically performing this
- conveying of source, or (2) access to copy the
- Corresponding Source from a network server at no charge.
-
- c) Convey individual copies of the object code with a copy of the
- written offer to provide the Corresponding Source. This
- alternative is allowed only occasionally and noncommercially, and
- only if you received the object code with such an offer, in accord
- with subsection 6b.
-
- d) Convey the object code by offering access from a designated
- place (gratis or for a charge), and offer equivalent access to the
- Corresponding Source in the same way through the same place at no
- further charge. You need not require recipients to copy the
- Corresponding Source along with the object code. If the place to
- copy the object code is a network server, the Corresponding Source
- may be on a different server (operated by you or a third party)
- that supports equivalent copying facilities, provided you maintain
- clear directions next to the object code saying where to find the
- Corresponding Source. Regardless of what server hosts the
- Corresponding Source, you remain obligated to ensure that it is
- available for as long as needed to satisfy these requirements.
-
- e) Convey the object code using peer-to-peer transmission, provided
- you inform other peers where the object code and Corresponding
- Source of the work are being offered to the general public at no
- charge under subsection 6d.
-
- A separable portion of the object code, whose source code is excluded
-from the Corresponding Source as a System Library, need not be
-included in conveying the object code work.
-
- A "User Product" is either (1) a "consumer product", which means any
-tangible personal property which is normally used for personal, family,
-or household purposes, or (2) anything designed or sold for incorporation
-into a dwelling. In determining whether a product is a consumer product,
-doubtful cases shall be resolved in favor of coverage. For a particular
-product received by a particular user, "normally used" refers to a
-typical or common use of that class of product, regardless of the status
-of the particular user or of the way in which the particular user
-actually uses, or expects or is expected to use, the product. A product
-is a consumer product regardless of whether the product has substantial
-commercial, industrial or non-consumer uses, unless such uses represent
-the only significant mode of use of the product.
-
- "Installation Information" for a User Product means any methods,
-procedures, authorization keys, or other information required to install
-and execute modified versions of a covered work in that User Product from
-a modified version of its Corresponding Source. The information must
-suffice to ensure that the continued functioning of the modified object
-code is in no case prevented or interfered with solely because
-modification has been made.
-
- If you convey an object code work under this section in, or with, or
-specifically for use in, a User Product, and the conveying occurs as
-part of a transaction in which the right of possession and use of the
-User Product is transferred to the recipient in perpetuity or for a
-fixed term (regardless of how the transaction is characterized), the
-Corresponding Source conveyed under this section must be accompanied
-by the Installation Information. But this requirement does not apply
-if neither you nor any third party retains the ability to install
-modified object code on the User Product (for example, the work has
-been installed in ROM).
-
- The requirement to provide Installation Information does not include a
-requirement to continue to provide support service, warranty, or updates
-for a work that has been modified or installed by the recipient, or for
-the User Product in which it has been modified or installed. Access to a
-network may be denied when the modification itself materially and
-adversely affects the operation of the network or violates the rules and
-protocols for communication across the network.
-
- Corresponding Source conveyed, and Installation Information provided,
-in accord with this section must be in a format that is publicly
-documented (and with an implementation available to the public in
-source code form), and must require no special password or key for
-unpacking, reading or copying.
-
- 7. Additional Terms.
-
- "Additional permissions" are terms that supplement the terms of this
-License by making exceptions from one or more of its conditions.
-Additional permissions that are applicable to the entire Program shall
-be treated as though they were included in this License, to the extent
-that they are valid under applicable law. If additional permissions
-apply only to part of the Program, that part may be used separately
-under those permissions, but the entire Program remains governed by
-this License without regard to the additional permissions.
-
- When you convey a copy of a covered work, you may at your option
-remove any additional permissions from that copy, or from any part of
-it. (Additional permissions may be written to require their own
-removal in certain cases when you modify the work.) You may place
-additional permissions on material, added by you to a covered work,
-for which you have or can give appropriate copyright permission.
-
- Notwithstanding any other provision of this License, for material you
-add to a covered work, you may (if authorized by the copyright holders of
-that material) supplement the terms of this License with terms:
-
- a) Disclaiming warranty or limiting liability differently from the
- terms of sections 15 and 16 of this License; or
-
- b) Requiring preservation of specified reasonable legal notices or
- author attributions in that material or in the Appropriate Legal
- Notices displayed by works containing it; or
-
- c) Prohibiting misrepresentation of the origin of that material, or
- requiring that modified versions of such material be marked in
- reasonable ways as different from the original version; or
-
- d) Limiting the use for publicity purposes of names of licensors or
- authors of the material; or
-
- e) Declining to grant rights under trademark law for use of some
- trade names, trademarks, or service marks; or
-
- f) Requiring indemnification of licensors and authors of that
- material by anyone who conveys the material (or modified versions of
- it) with contractual assumptions of liability to the recipient, for
- any liability that these contractual assumptions directly impose on
- those licensors and authors.
-
- All other non-permissive additional terms are considered "further
-restrictions" within the meaning of section 10. If the Program as you
-received it, or any part of it, contains a notice stating that it is
-governed by this License along with a term that is a further
-restriction, you may remove that term. If a license document contains
-a further restriction but permits relicensing or conveying under this
-License, you may add to a covered work material governed by the terms
-of that license document, provided that the further restriction does
-not survive such relicensing or conveying.
-
- If you add terms to a covered work in accord with this section, you
-must place, in the relevant source files, a statement of the
-additional terms that apply to those files, or a notice indicating
-where to find the applicable terms.
-
- Additional terms, permissive or non-permissive, may be stated in the
-form of a separately written license, or stated as exceptions;
-the above requirements apply either way.
-
- 8. Termination.
-
- You may not propagate or modify a covered work except as expressly
-provided under this License. Any attempt otherwise to propagate or
-modify it is void, and will automatically terminate your rights under
-this License (including any patent licenses granted under the third
-paragraph of section 11).
-
- However, if you cease all violation of this License, then your
-license from a particular copyright holder is reinstated (a)
-provisionally, unless and until the copyright holder explicitly and
-finally terminates your license, and (b) permanently, if the copyright
-holder fails to notify you of the violation by some reasonable means
-prior to 60 days after the cessation.
-
- Moreover, your license from a particular copyright holder is
-reinstated permanently if the copyright holder notifies you of the
-violation by some reasonable means, this is the first time you have
-received notice of violation of this License (for any work) from that
-copyright holder, and you cure the violation prior to 30 days after
-your receipt of the notice.
-
- Termination of your rights under this section does not terminate the
-licenses of parties who have received copies or rights from you under
-this License. If your rights have been terminated and not permanently
-reinstated, you do not qualify to receive new licenses for the same
-material under section 10.
-
- 9. Acceptance Not Required for Having Copies.
-
- You are not required to accept this License in order to receive or
-run a copy of the Program. Ancillary propagation of a covered work
-occurring solely as a consequence of using peer-to-peer transmission
-to receive a copy likewise does not require acceptance. However,
-nothing other than this License grants you permission to propagate or
-modify any covered work. These actions infringe copyright if you do
-not accept this License. Therefore, by modifying or propagating a
-covered work, you indicate your acceptance of this License to do so.
-
- 10. Automatic Licensing of Downstream Recipients.
-
- Each time you convey a covered work, the recipient automatically
-receives a license from the original licensors, to run, modify and
-propagate that work, subject to this License. You are not responsible
-for enforcing compliance by third parties with this License.
-
- An "entity transaction" is a transaction transferring control of an
-organization, or substantially all assets of one, or subdividing an
-organization, or merging organizations. If propagation of a covered
-work results from an entity transaction, each party to that
-transaction who receives a copy of the work also receives whatever
-licenses to the work the party's predecessor in interest had or could
-give under the previous paragraph, plus a right to possession of the
-Corresponding Source of the work from the predecessor in interest, if
-the predecessor has it or can get it with reasonable efforts.
-
- You may not impose any further restrictions on the exercise of the
-rights granted or affirmed under this License. For example, you may
-not impose a license fee, royalty, or other charge for exercise of
-rights granted under this License, and you may not initiate litigation
-(including a cross-claim or counterclaim in a lawsuit) alleging that
-any patent claim is infringed by making, using, selling, offering for
-sale, or importing the Program or any portion of it.
-
- 11. Patents.
-
- A "contributor" is a copyright holder who authorizes use under this
-License of the Program or a work on which the Program is based. The
-work thus licensed is called the contributor's "contributor version".
-
- A contributor's "essential patent claims" are all patent claims
-owned or controlled by the contributor, whether already acquired or
-hereafter acquired, that would be infringed by some manner, permitted
-by this License, of making, using, or selling its contributor version,
-but do not include claims that would be infringed only as a
-consequence of further modification of the contributor version. For
-purposes of this definition, "control" includes the right to grant
-patent sublicenses in a manner consistent with the requirements of
-this License.
-
- Each contributor grants you a non-exclusive, worldwide, royalty-free
-patent license under the contributor's essential patent claims, to
-make, use, sell, offer for sale, import and otherwise run, modify and
-propagate the contents of its contributor version.
-
- In the following three paragraphs, a "patent license" is any express
-agreement or commitment, however denominated, not to enforce a patent
-(such as an express permission to practice a patent or covenant not to
-sue for patent infringement). To "grant" such a patent license to a
-party means to make such an agreement or commitment not to enforce a
-patent against the party.
-
- If you convey a covered work, knowingly relying on a patent license,
-and the Corresponding Source of the work is not available for anyone
-to copy, free of charge and under the terms of this License, through a
-publicly available network server or other readily accessible means,
-then you must either (1) cause the Corresponding Source to be so
-available, or (2) arrange to deprive yourself of the benefit of the
-patent license for this particular work, or (3) arrange, in a manner
-consistent with the requirements of this License, to extend the patent
-license to downstream recipients. "Knowingly relying" means you have
-actual knowledge that, but for the patent license, your conveying the
-covered work in a country, or your recipient's use of the covered work
-in a country, would infringe one or more identifiable patents in that
-country that you have reason to believe are valid.
-
- If, pursuant to or in connection with a single transaction or
-arrangement, you convey, or propagate by procuring conveyance of, a
-covered work, and grant a patent license to some of the parties
-receiving the covered work authorizing them to use, propagate, modify
-or convey a specific copy of the covered work, then the patent license
-you grant is automatically extended to all recipients of the covered
-work and works based on it.
-
- A patent license is "discriminatory" if it does not include within
-the scope of its coverage, prohibits the exercise of, or is
-conditioned on the non-exercise of one or more of the rights that are
-specifically granted under this License. You may not convey a covered
-work if you are a party to an arrangement with a third party that is
-in the business of distributing software, under which you make payment
-to the third party based on the extent of your activity of conveying
-the work, and under which the third party grants, to any of the
-parties who would receive the covered work from you, a discriminatory
-patent license (a) in connection with copies of the covered work
-conveyed by you (or copies made from those copies), or (b) primarily
-for and in connection with specific products or compilations that
-contain the covered work, unless you entered into that arrangement,
-or that patent license was granted, prior to 28 March 2007.
-
- Nothing in this License shall be construed as excluding or limiting
-any implied license or other defenses to infringement that may
-otherwise be available to you under applicable patent law.
-
- 12. No Surrender of Others' Freedom.
-
- If conditions are imposed on you (whether by court order, agreement or
-otherwise) that contradict the conditions of this License, they do not
-excuse you from the conditions of this License. If you cannot convey a
-covered work so as to satisfy simultaneously your obligations under this
-License and any other pertinent obligations, then as a consequence you may
-not convey it at all. For example, if you agree to terms that obligate you
-to collect a royalty for further conveying from those to whom you convey
-the Program, the only way you could satisfy both those terms and this
-License would be to refrain entirely from conveying the Program.
-
- 13. Use with the GNU Affero General Public License.
-
- Notwithstanding any other provision of this License, you have
-permission to link or combine any covered work with a work licensed
-under version 3 of the GNU Affero General Public License into a single
-combined work, and to convey the resulting work. The terms of this
-License will continue to apply to the part which is the covered work,
-but the special requirements of the GNU Affero General Public License,
-section 13, concerning interaction through a network will apply to the
-combination as such.
-
- 14. Revised Versions of this License.
-
- The Free Software Foundation may publish revised and/or new versions of
-the GNU General Public License from time to time. Such new versions will
-be similar in spirit to the present version, but may differ in detail to
-address new problems or concerns.
-
- Each version is given a distinguishing version number. If the
-Program specifies that a certain numbered version of the GNU General
-Public License "or any later version" applies to it, you have the
-option of following the terms and conditions either of that numbered
-version or of any later version published by the Free Software
-Foundation. If the Program does not specify a version number of the
-GNU General Public License, you may choose any version ever published
-by the Free Software Foundation.
-
- If the Program specifies that a proxy can decide which future
-versions of the GNU General Public License can be used, that proxy's
-public statement of acceptance of a version permanently authorizes you
-to choose that version for the Program.
-
- Later license versions may give you additional or different
-permissions. However, no additional obligations are imposed on any
-author or copyright holder as a result of your choosing to follow a
-later version.
-
- 15. Disclaimer of Warranty.
-
- THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
-APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
-HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
-OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
-THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
-PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
-IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
-ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
-
- 16. Limitation of Liability.
-
- IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
-WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
-THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
-GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
-USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
-DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
-PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
-EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
-SUCH DAMAGES.
-
- 17. Interpretation of Sections 15 and 16.
-
- If the disclaimer of warranty and limitation of liability provided
-above cannot be given local legal effect according to their terms,
-reviewing courts shall apply local law that most closely approximates
-an absolute waiver of all civil liability in connection with the
-Program, unless a warranty or assumption of liability accompanies a
-copy of the Program in return for a fee.
-
- END OF TERMS AND CONDITIONS
-
- How to Apply These Terms to Your New Programs
-
- If you develop a new program, and you want it to be of the greatest
-possible use to the public, the best way to achieve this is to make it
-free software which everyone can redistribute and change under these terms.
-
- To do so, attach the following notices to the program. It is safest
-to attach them to the start of each source file to most effectively
-state the exclusion of warranty; and each file should have at least
-the "copyright" line and a pointer to where the full notice is found.
-
- in20tools
- Copyright (C) 2018 Tobias WEBER
-
- This program is free software: you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-Also add information on how to contact you by electronic and paper mail.
-
- If the program does terminal interaction, make it output a short
-notice like this when it starts in an interactive mode:
-
- in20tools Copyright (C) 2018 Tobias WEBER
- This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
- This is free software, and you are welcome to redistribute it
- under certain conditions; type `show c' for details.
-
-The hypothetical commands `show w' and `show c' should show the appropriate
-parts of the General Public License. Of course, your program's commands
-might be different; for a GUI interface, you would use an "about box".
-
- You should also get your employer (if you work as a programmer) or school,
-if any, to sign a "copyright disclaimer" for the program, if necessary.
-For more information on this, and how to apply and follow the GNU GPL, see
-<http://www.gnu.org/licenses/>.
-
- The GNU General Public License does not permit incorporating your program
-into proprietary programs. If your program is a subroutine library, you
-may consider it more useful to permit linking proprietary applications with
-the library. If this is what you want to do, use the GNU Lesser General
-Public License instead of this License. But first, please read
-<http://www.gnu.org/philosophy/why-not-lgpl.html>.
diff --git a/ext/tlibs2-mag/mag.h b/ext/tlibs2-mag/mag.h
deleted file mode 100644
index 4eb4705..0000000
--- a/ext/tlibs2-mag/mag.h
+++ /dev/null
@@ -1,199 +0,0 @@
-/**
- * tlibs2 -- magnon helper functions in preparation for the main tlibs2 repo: https://code.ill.fr/scientific-software/takin/tlibs2
- * @author Tobias Weber <tweber@ill.fr>
- * @date 30-may-2020
- * @license GPLv2 or GPLv3, see 'LICENSE' file
- * @desc based on the mathematical formalism of the skyrmion and helimagnon models by M. Garst and J. Waizner, 2016-2020, references:
- * - https://doi.org/10.1088/1361-6463/aa7573
- * - https://doi.org/10.1038/nmat4223 (supplement)
- * - https://kups.ub.uni-koeln.de/7937/
- * @desc forked on 7-Nov-2018 from my privately developed "tlibs" project (https://github.com/t-weber/tlibs).
- */
-
-#ifndef __TLIBS2_PHYS_MAG__
-#define __TLIBS2_PHYS_MAG__
-
-#include "math17.h"
-
-
-// default optimisation flags
-#ifndef MX_IS_HERM
- #define MX_IS_HERM 1
-#endif
-#ifndef MXX_IS_DIAG
- #define MXX_IS_DIAG 1
-#endif
-#ifndef INTERACTMAT_IS_HERM
- // warning: setting this to "1" leads to false low-q eigenvectors!
- #define INTERACTMAT_IS_HERM 0
-#endif
-
-
-#if !defined(__TL2_STRCONV0) && !defined(__TL2_STRCONV)
- #define __TL2_STRCONV0(__DEF) #__DEF
- #define __TL2_STRCONV(__DEF) __TL2_STRCONV0(__DEF)
-#endif
-
-#pragma message("MX_IS_HERM: " __TL2_STRCONV(MX_IS_HERM))
-#pragma message("MXX_IS_DIAG: " __TL2_STRCONV(MXX_IS_DIAG))
-#pragma message("INTERACTMAT_IS_HERM: " __TL2_STRCONV(INTERACTMAT_IS_HERM))
-
-
-namespace tl2 {
-
-/**
- * Calculates energies and dynamical structure factors from Landau-Lifshitz (M x) and fluctuation matrices.
- */
-template<class t_mat_cplx, class t_vec_cplx, class t_cplx, class t_real>
-std::tuple<std::vector<t_cplx>, std::vector<t_vec_cplx>, std::vector<t_mat_cplx>>
-calc_dynstrucfact_landau(const t_mat_cplx& Mx, const t_mat_cplx& Fluc,
- t_real normfac=1, const t_real* mineval = nullptr, const t_real *maxeval = nullptr,
- std::size_t MxsubMatSize=3, std::size_t MxsubMatRowBegin=0, t_real eps=1e-6)
-{
- constexpr t_cplx imag = t_cplx(0,1);
-
- // calculate Mx eigenvalues
- std::vector<t_vec_cplx> Mxevecs;
- std::vector<t_cplx> Mxevals;
- {
-#if MX_IS_HERM == 1
- std::vector<t_real> _Mxevals;
- if(!eigenvecsel_herm<t_real>(-imag*Mx, Mxevecs, _Mxevals, true, -1., -2., eps))
- throw std::runtime_error("Mx eigenvector determination failed!");
- for(t_real d : _Mxevals)
- Mxevals.emplace_back(t_cplx(0., d));
-#else
- if(!eigenvec_cplx<t_real>(Mx, Mxevecs, Mxevals, true))
- throw std::runtime_error("Mx eigenvector determination failed!");
-#endif
-
- // filter eigenvalues
- auto maxelem = std::max_element(Mxevals.begin(), Mxevals.end(),
- [](const t_cplx& x, const t_cplx& y) -> bool
- { return std::abs(x.imag()) < std::abs(y.imag()); });
-
- std::vector<t_vec_cplx> Mxevecs_new;
- std::vector<t_cplx> Mxevals_new;
-
- for(std::size_t elem=0; elem<Mxevals.size(); ++elem)
- {
- // upper eigenvalue limit
- if(maxeval && std::abs(Mxevals[elem].imag()) < std::abs(*maxelem)**maxeval)
- continue;
- // lower eigenvalue limit
- if(mineval && std::abs(Mxevals[elem].imag()) < *mineval)
- continue;
- Mxevecs_new.push_back(Mxevecs[elem]);
- Mxevals_new.push_back(Mxevals[elem]);
- }
-
- Mxevecs = std::move(Mxevecs_new);
- Mxevals = std::move(Mxevals_new);
- }
-
-
- // convert to eigenvector matrix
- t_mat_cplx MxEvecs(Mx.size1(), Mxevecs.size());
- for(std::size_t idx1=0; idx1<MxEvecs.size1(); ++idx1)
- for(std::size_t idx2=0; idx2<MxEvecs.size2(); ++idx2)
- MxEvecs(idx1, idx2) = Mxevecs[idx2][idx1];
-
- t_mat_cplx MxEvecsH = hermitian(MxEvecs);
- t_mat_cplx MxEvecs3 = submatrix_wnd<t_mat_cplx>(MxEvecs, MxsubMatSize, Mxevecs.size(), MxsubMatRowBegin, 0);
- t_mat_cplx MxEvecsH3 = hermitian(MxEvecs3);
-
- // transform fluctuation matrix into Mx eigenvector system
- t_mat_cplx invsuscept = prod_mm(Fluc, MxEvecs);
- invsuscept = prod_mm(MxEvecsH, invsuscept);
-
- // transform Mx into Mx eigenvector system
- // Mxx is diagonal with this construction => directly use Mxevals
-#if MXX_IS_DIAG == 1
- t_mat_cplx Mxx = diag_matrix<t_mat_cplx>(Mxevals);
-#else
- t_mat_cplx Mxx = prod_mm(Mx, MxEvecs);
- Mxx = prod_mm(MxEvecsH, Mxx);
-#endif
- Mxx *= imag / normfac;
-
- // Landau-Lifshitz: d/dt dM = -Mx B_mean, B_mean = -chi^(-1) * dM
- // E = EVals{ i Mx chi^(-1) }
- // chi_dyn^(-1) = i*E*Mx^(-1) + chi^(-1)
- // Mx*chi_dyn^(-1) = i*E + Mx*chi^(-1)
- t_mat_cplx Interactmat = prod_mm(Mxx, invsuscept);
- std::vector<t_vec_cplx> Interactevecs;
- std::vector<t_cplx> Interactevals;
-#if INTERACTMAT_IS_HERM == 1
- std::vector<t_real> _Interactevals;
- if(!eigenvecsel_herm<t_real>(Interactmat, Interactevecs, _Interactevals, true, -1., -2., eps))
- throw std::runtime_error("Mxx eigenvector determination failed!");
- for(t_real d : _Interactevals)
- Interactevals.emplace_back(t_cplx(d, 0.));
-#else
- if(!eigenvec_cplx<t_real>(Interactmat, Interactevecs, Interactevals, true))
- throw std::runtime_error("Mxx eigenvector determination failed!");
-#endif
-
- std::vector<t_mat_cplx> Interactemats;
- Interactemats.reserve(Interactevals.size());
-
- for(std::size_t iInteract=0; iInteract<Interactevals.size(); ++iInteract)
- {
- const t_cplx& eval = Interactevals[iInteract];
- const t_vec_cplx& evec = Interactevecs[iInteract];
-
- auto evec_scale = prod_mv(Mxx, evec);
-
- auto matOuter = outer_cplx<t_vec_cplx, t_mat_cplx>(evec, evec);
- matOuter /= inner_cplx<t_vec_cplx>(evec, evec_scale);
-
- t_mat_cplx emat = prod_mm(matOuter, MxEvecsH3);
- emat = prod_mm(MxEvecs3, emat);
- Interactemats.emplace_back(std::move(emat));
- //eigs.emplace_back(Eig{.eval=eval, .evec=evec, .emat=emat});
- }
-
- return std::make_tuple(Interactevals, Interactevecs, Interactemats);
-}
-
-
-
-/**
- * Gets the dynamical structure factors from the eigenvectors calculated using calc_dynstrucfact_landau.
- */
-template<class t_mat_cplx, class t_vec_cplx, class t_cplx, class t_real>
-std::tuple<t_real, std::vector<t_real>>
-get_dynstrucfact_neutron(
- const t_cplx& eval, const t_vec_cplx& evec, const t_mat_cplx& _emat,
- const t_mat_cplx* projNeutron=nullptr, const std::vector<t_mat_cplx>* pol = nullptr)
-{
- t_real E = eval.real();
- t_mat_cplx emat = _emat;
-
- // magnetic neutron scattering orthogonal projector: projNeutron = 1 - |G><G|
- if(projNeutron)
- {
- emat = prod_mm(emat, *projNeutron);
- emat = prod_mm(*projNeutron, emat);
- }
-
- std::vector<t_real> sfacts;
-
- // unpolarised structure factor
- sfacts.push_back(std::abs(trace(emat).real()));
-
- // polarised structure factors
- if(pol)
- {
- for(const t_mat_cplx& polmat : *pol)
- {
- t_mat_cplx matSF = prod_mm(polmat, emat);
- sfacts.push_back(std::abs(trace(matSF).real()));
- }
- }
-
- return std::make_tuple(E, sfacts);
-}
-
-}
-#endif
diff --git a/src/calc/drawskx.cpp b/src/calc/drawskx.cpp
deleted file mode 100644
index de607c1..0000000
--- a/src/calc/drawskx.cpp
+++ /dev/null
@@ -1,78 +0,0 @@
-/**
- * @author Tobias Weber <tweber@ill.fr>
- * @date mid-2016 (from my phd thesis), oct-2018 (adapted for paper)
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#include <fstream>
-#include "tlibs2/libs/math17.h"
-
-using t_real = double;
-using t_vec = ublas::vector<t_real>;
-using t_mat = ublas::matrix<t_real>;
-
-
-
-t_vec helix_vec(const t_vec& vecCoord)
-{
- t_vec vecRet = ublas::zero_vector<t_real>(3);
-
- // helix angles
- for(t_real dAngle : {0., 120., 240.})
- {
- dAngle = tl2::d2r(dAngle);
- t_mat matRot = tl2::rotation_matrix_3d_z(dAngle);
- t_vec vecCoordRot = ublas::prod(ublas::trans(matRot), vecCoord);
-
- // helix position
- t_vec vec(3);
- vec[2] = std::cos(vecCoordRot[0]);
- vec[1] = std::sin(vecCoordRot[0]);
- vec[0] = 0;
-
- vecRet += ublas::prod(matRot, vec);
- }
-
- vecRet /= ublas::norm_2(vecRet);
- return vecRet;
-}
-
-
-
-void calcskx(const char* pcFile)
-{
- const int iCntX = 24, iCntY = 24;
- const t_real dXScale = 2., dYScale = 2.;
- const t_real dDirScale = 0.3;
-
- std::ofstream ofstr(pcFile);
- ofstr.precision(8);
-
- for(int iX=0; iX<iCntX; ++iX)
- {
- t_real dX = -tl2::pi<t_real> + t_real(iX)/t_real(iCntX-1) * 2.*tl2::pi<t_real>;
-
- for(int iY=0; iY<iCntY; ++iY)
- {
- t_real dY = -tl2::pi<t_real> + t_real(iY)/t_real(iCntY-1) * 2.*tl2::pi<t_real>;
-
- t_vec vecCoord(3);
- vecCoord[0] = dXScale*dX;
- vecCoord[1] = dYScale*dY;
- vecCoord[2] = 0.;
-
- t_vec vec = helix_vec(vecCoord)*dDirScale;
-
- ofstr << vecCoord[0] << "\t" << vecCoord[1] << "\t" << vecCoord[2] << "\t\t";
- ofstr << vec[0] << "\t" << vec[1] << "\t" << vec[2] << std::endl;
- }
- }
-}
-
-
-
-int main(int argc, char **argv)
-{
- calcskx("drawskx.dat");
- return 0;
-}
diff --git a/src/calc/drawskx.gpl b/src/calc/drawskx.gpl
deleted file mode 100644
index 17d2f2c..0000000
--- a/src/calc/drawskx.gpl
+++ /dev/null
@@ -1,58 +0,0 @@
-#!/usr/local/bin/gnuplot -p
-#
-# Plots skx lattice
-# @author Tobias Weber <tweber@ill.fr>
-# @date mid-16 (from my phd thesis), aug-18
-# @license GPLv2 (see 'LICENSE' file)
-#
-
-set term pdf enhanced font "Helvetica,24" size 10,10
-set output "skxlattice.pdf"
-
-show_axes = 1
-sc = 1
-xoffs = pi
-set xrange [(-4*pi-xoffs)*sc : (4*pi-xoffs)*sc]
-set yrange [-4*pi*sc : 4*pi*sc]
-set zrange [-1.5*sc : 1.5*sc]
-
-set xyplane 0
-unset key
-
-unset label 1
-unset label 2
-unset label 10
-unset obj 1
-unset obj 2
-unset obj 10
-
-unset xlabel
-unset ylabel
-unset zlabel
-
-set format x ""
-set format y ""
-set format z ""
-
-unset border
-unset tics
-
-set view 45,45, 1,1
-set pm3d depthorder hidden3d
-
-if(show_axes) {
- set arrow 1 from graph 0.55, graph 0.5 to graph 0.55, graph 0.75 ls 1 lc rgb "#000000" lw 5 filled head front
- set arrow 2 from graph 0.55, graph 0.5 to graph 0.3, graph 0.5 ls 1 lc rgb "#000000" lw 5 filled head front
- set arrow 31 from graph 0.55, graph 0.5, graph 0 to graph 0.55, graph 0.5, graph 0.25 ls 1 lc rgb "#ffffff" lw 8 filled head front
- set arrow 32 from graph 0.55, graph 0.5, graph 0 to graph 0.55, graph 0.5, graph 0.25 ls 1 lc rgb "#000000" lw 5 filled head front
-
- set label 1 at graph 0.62, graph 0.75 center "q_{⟂1}" tc rgb "#000000" front
- set label 2 at graph 0.3, graph 0.43 center "q_{⟂2}" tc rgb "#000000" front
- set label 3 at graph 0.55, graph 0.4, graph 0.25 center "q_{||}" tc rgb "#000000" front
-}
-
-colrgb(r,g,b) = ((int(r*255)&0xff)<<16) | ((int(g*255)&0xff)<<8) | (int(b*255)&0xff)
-skxcol(len, zval) = (zval < 0) ? colrgb(len*1,0,-zval*2) : colrgb(len*1,zval*2,0)
-
-splot "drawskx.dat" using ($1*sc):($2*sc):($3*sc):4:5:6:(skxcol(sqrt($4*$4+$5*$5), $6)) \
- with vectors lc rgb variable lw 2 filled head
diff --git a/src/calc/dyn.cpp b/src/calc/dyn.cpp
deleted file mode 100644
index 0551961..0000000
--- a/src/calc/dyn.cpp
+++ /dev/null
@@ -1,300 +0,0 @@
-/**
- * Calculates dispersion curves
- * @author Tobias Weber <tweber@ill.fr>
- * @date aug-18
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#include "core/heli.h"
-#include "core/skx.h"
-#include "core/fp.h"
-
-#include <fstream>
-#include <future>
-#include <memory>
-
-#ifndef __MINGW32__
- #include <pwd.h>
-#endif
-
-using t_real = double;
-using t_cplx = std::complex<t_real>;
-using t_vec = ublas::vector<t_real>;
-const auto g_j = t_cplx(0,1);
-
-#include "core/skx_default_gs.cxx"
-
-
-void calc_disp(char dyntype,
- t_real Gx, t_real Gy, t_real Gz,
- t_real Bx, t_real By, t_real Bz,
- t_real Px, t_real Py, t_real Pz,
- t_real qperpdir_x, t_real qperpdir_y, t_real qperpdir_z, t_real qperp,
- const char* pcFile, bool bSwapQParaQPerp=0,
- t_real T=-1., t_real B=-1.,
- t_real qrange = 0.125, t_real delta = 0.001)
-{
- tl2::Stopwatch<t_real> timer;
- timer.start();
-
- t_vec G = tl2::make_vec<t_vec>({ Gx, Gy, Gz });
- t_vec Pdir = tl2::make_vec<t_vec>({ Px, Py, Pz });
- t_vec Bdir = tl2::make_vec<t_vec>({ Bx, By, Bz });
- t_vec qparadir = Bdir / tl2::veclen(Bdir);
- t_vec qperpdir = tl2::make_vec<t_vec>({ qperpdir_x, qperpdir_y, qperpdir_z });
- qperpdir /= tl2::veclen(qperpdir);
-
- t_real bc2 = -1.;
- std::shared_ptr<MagDynamics<t_real, t_cplx>> dyn;
-
- if(dyntype == 's')
- {
- std::cout << "Calculating skyrmion dispersion." << std::endl;
- auto skx = std::make_shared<Skx<t_real, t_cplx, DEF_SKX_ORDER>>();
-
- std::vector<ublas::vector<t_cplx>> fourier_skx;
- fourier_skx.reserve(_skxgs_allcomps.size()/3);
-
- for(std::size_t comp=0; comp<_skxgs_allcomps.size(); comp+=3)
- fourier_skx.push_back(tl2::make_vec<ublas::vector<t_cplx>>({_skxgs_allcomps[comp], _skxgs_allcomps[comp+1], _skxgs_allcomps[comp+2]}));
-
- skx->SetFourier(fourier_skx);
-
- skx->SetT(-1000.);
- skx->SetB(25.); // BC2 = 40.3425
- skx->GenFullFourier();
- skx->SetCoords(Bdir[0],Bdir[1],Bdir[2], Pdir[0],Pdir[1],Pdir[2]);
- skx->SetG(G[0], G[1], G[2]);
-
- bc2 = skx->GetBC2();
- dyn = skx;
- }
- else if(dyntype == 'h')
- {
- std::cout << "Calculating helical dispersion." << std::endl;
- auto heli = std::make_shared<Heli<t_real, t_cplx, DEF_HELI_ORDER>>();
-
- heli->SetT(T);
- heli->SetB(B);
- heli->SetCoords(Bdir[0], Bdir[1], Bdir[2]);
- heli->SetG(G[0], G[1], G[2]);
-
- bc2 = heli->GetBC2();
- dyn = heli;
- }
- else if(dyntype == 'f')
- {
- std::cout << "Calculating field-polarised dispersion." << std::endl;
- auto fp = std::make_shared<FP<t_real, t_cplx>>();
-
- fp->SetT(T);
- fp->SetB(B);
- fp->SetCoords(Bdir[0], Bdir[1], Bdir[2]);
- fp->SetG(G[0], G[1], G[2]);
-
- bc2 = fp->GetBC2();
- dyn = fp;
- }
- else
- {
- std::cerr << "Unknown dynamics type selected." << std::endl;
- return;
- }
-
-
- auto calc_spectrum = [dyntype, &dyn, &G, T, &qparadir, &qperpdir, &qperp, bSwapQParaQPerp]
- (int thid, t_real qstart, t_real qend, t_real qdelta) -> auto
- {
- std::vector<t_real> allh, allk, alll, allqpara_kh, allqperp_kh, allqpara_rlu, allqperp_rlu;
- std::vector<std::vector<t_real>> allEs, allWsUnpol, allWsSF1, allWsSF2, allWsNSF;
- auto thisdyn = dyn->copyCastDyn();
-
- for(t_real _q=qstart; _q<qend; _q+=qdelta)
- {
- t_real qpara = _q;
- t_vec Q = G + qpara*qparadir + qperp*qperpdir;
- if(bSwapQParaQPerp)
- Q = G + qpara*qperpdir + qperp*qparadir; // swap qpara and qperp
-
- std::cout << "thread " << thid << " (" << Q[0] << " " << Q[1] << " " << Q[2] << ") ... ";
- std::cout.flush();
-
- auto [Es, wsUnpol, wsSF1, wsSF2, wsNSF] = thisdyn->GetDisp(Q[0], Q[1], Q[2]);
- allEs.emplace_back(Es);
- allWsUnpol.emplace_back(wsUnpol);
- allWsSF1.emplace_back(wsSF1);
- allWsSF2.emplace_back(wsSF2);
- allWsNSF.emplace_back(wsNSF);
-
- allh.push_back(Q[0]);
- allk.push_back(Q[1]);
- alll.push_back(Q[2]);
-
- if(dyntype == 's')
- {
- allqpara_kh.push_back(qpara / g_kh_rlu_29K<t_real>);
- allqperp_kh.push_back(qperp / g_kh_rlu_29K<t_real>);
- }
- else
- {
- allqpara_kh.push_back(qpara / g_kh_rlu<t_real>(T));
- allqperp_kh.push_back(qperp / g_kh_rlu<t_real>(T));
- }
-
- allqpara_rlu.push_back(qpara);
- allqperp_rlu.push_back(qperp);
-
- std::cout << "done." << std::endl;
- }
-
- return std::make_tuple(allh, allk, alll, allEs,
- allWsUnpol, allWsSF1, allWsSF2, allWsNSF,
- allqpara_kh, allqperp_kh, allqpara_rlu, allqperp_rlu);
- };
-
-
- t_real qstep = qrange / 2.;
-
- auto fut0 = std::async(std::launch::async, calc_spectrum, 0, -qrange, -qrange+1.*qstep, delta);
- auto fut1 = std::async(std::launch::async, calc_spectrum, 1, -qrange+1.*qstep, -qrange+2.*qstep, delta);
- auto fut2 = std::async(std::launch::async, calc_spectrum, 2, -qrange+2.*qstep, -qrange+3.*qstep, delta);
- auto fut3 = std::async(std::launch::async, calc_spectrum, 3, -qrange+3.*qstep, -qrange+4.*qstep, delta);
- auto val0 = fut0.get();
- auto val1 = fut1.get();
- auto val2 = fut2.get();
- auto val3 = fut3.get();
-
- insert_vals(val0, val1, std::make_index_sequence<std::tuple_size<decltype(val1)>::value>());
- insert_vals(val0, val2, std::make_index_sequence<std::tuple_size<decltype(val2)>::value>());
- insert_vals(val0, val3, std::make_index_sequence<std::tuple_size<decltype(val3)>::value>());
-
-
- timer.stop();
- std::ofstream ofstr(pcFile);
- ofstr.precision(8);
-
- ofstr << "#\n";
- ofstr << "# Date: " << tl2::epoch_to_str<t_real>(tl2::epoch<t_real>()) << "\n";
-#ifndef __MINGW32__
- ofstr << "# User: " << getpwuid(geteuid())->pw_name << "\n";
-#endif
- ofstr << "# Calculation time: " << timer.GetDur() << " s.\n";
- ofstr << "# Skx order: " << DEF_SKX_ORDER << "\n";
- if(dyntype == 'h')
- {
- ofstr << "# Heli order: " << DEF_HELI_ORDER << "\n";
- ofstr << "# kh_A = " << g_kh_A<t_real>(T) << "\n";
- ofstr << "# kh_rlu = " << g_kh_rlu<t_real>(T) << "\n";
- }
- ofstr << "# qparadir = " << qparadir << "\n";
- ofstr << "# qperpdir = " << qperpdir << "\n";
- ofstr << "# qperp = " << qperp << "\n";
- ofstr << "# G = " << G << "\n";
- ofstr << "# Bdir = " << Bdir << "\n";
- ofstr << "# Pdir = " << Pdir << "\n";
- ofstr << "# Bc2 = " << bc2 << "\n";
- if(bSwapQParaQPerp)
- ofstr << "# WARNING: In the following, q_para_* and q_perp_* are swapped!\n";
-
- // save commutator of projectors to determine channel mixing
- /*ofstr << "# [Nrpoj, Polproj_sf1] = " << tl2::commutator<t_mat_cplx>(fp.GetNeutronProjOp(), get_chiralpol<t_mat_cplx>(1)) << "\n";
- ofstr << "# [Nrpoj, Polproj_sf2] = " << tl2::commutator<t_mat_cplx>(fp.GetNeutronProjOp(), get_chiralpol<t_mat_cplx>(2)) << "\n";
- ofstr << "# [Nrpoj, Polproj_nsf] = " << tl2::commutator<t_mat_cplx>(fp.GetNeutronProjOp(), get_chiralpol<t_mat_cplx>(3)) << "\n";*/
-
- ofstr << "#\n";
-
- ofstr
- << "#" << std::setw(15) << "h" << " "
- << std::setw(16) << "k" << " "
- << std::setw(16) << "l" << " "
- << std::setw(16) << "E" << " "
- << std::setw(16) << "w_unpol" << " "
- << std::setw(16) << "w_sf1" << " "
- << std::setw(16) << "w_sf2" << " "
- << std::setw(16) << "w_nsf" << " "
- << std::setw(16) << "q_para_kh" << " "
- << std::setw(16) << "q_perp_kh" << " "
- << std::setw(16) << "q_para_rlu" << " "
- << std::setw(16) << "q_perp_rlu\n";
-
- for(std::size_t i=0; i<std::get<0>(val0).size(); ++i)
- {
- for(std::size_t j=0; j<std::get<3>(val0)[i].size(); ++j)
- {
- ofstr << std::setw(16) << std::get<0>(val0)[i] << " " // h
- << std::setw(16) << std::get<1>(val0)[i] << " " // k
- << std::setw(16) << std::get<2>(val0)[i] << " " // l
- << std::setw(16) << std::get<3>(val0)[i][j] << " " // E
- << std::setw(16) << std::get<4>(val0)[i][j] << " " // w_unpol
- << std::setw(16) << std::get<5>(val0)[i][j] << " " // w_sf1
- << std::setw(16) << std::get<6>(val0)[i][j] << " " // w_sf2
- << std::setw(16) << std::get<7>(val0)[i][j] << " " // w_nsf
- << std::setw(16) << std::get<8>(val0)[i] << " " // q_para_kh
- << std::setw(16) << std::get<9>(val0)[i] << " " // q_perp_kh
- << std::setw(16) << std::get<10>(val0)[i] << " " // q_para_rlu
- << std::setw(16) << std::get<11>(val0)[i] << "\n"; // q_perp_rlu
- }
- }
-
-
- timer.stop();
- std::cout << "Calculation took " << timer.GetDur() << " s." << std::endl;
- std::cout << "Wrote \"" << pcFile << "\"" << std::endl;
-}
-
-
-int main()
-{
- std::cout
- << "--------------------------------------------------------------------------------\n"
- << "\tDispersion calculation tool,\n\t\tT. Weber <tweber@ill.fr>, August 2018.\n"
- << "--------------------------------------------------------------------------------\n\n";
-
- char dyntype = 's';
- t_real Gx = 1., Gy = 1., Gz = 0.;
- t_real Bx = 1., By = 1., Bz = 0.;
- t_real Px = -1., Py = 1., Pz = 0.;
- t_real qperpx = -1., qperpy = 1., qperpz = 0.;
- t_real qperp = 0.;
- t_real B = 0.17, T = 28.5;
- t_real qrange = 0.125;
- t_real qdelta = 0.001;
- int alongqpara = 0;
-
- std::cout << "Helimagnon [h], skyrmion [s] or field-polarised [f] dynamics: ";
- std::cin >> dyntype; dyntype = std::tolower(dyntype);
- std::cout << "G = ";
- std::cin >> Gx >> Gy >> Gz;
- std::cout << "q_range = ";
- std::cin >> qrange;
- std::cout << "q_delta = ";
- std::cin >> qdelta;
- std::cout << "B = ";
- std::cin >> Bx >> By >> Bz;
- if(dyntype == 'h' || dyntype == 'f')
- {
- std::cout << "|B| = ";
- std::cin >> B;
- std::cout << "T = ";
- std::cin >> T;
- }
- else if(dyntype == 's')
- {
- std::cout << "pinning = ";
- std::cin >> Px >> Py >> Pz;
- }
- std::cout << "Query dispersion along q_para || B? [1/0]: ";
- std::cin >> alongqpara;
- std::cout << "q_perp = ";
- std::cin >> qperpx >> qperpy >> qperpz;
- std::cout << "|q_perp| = ";
- std::cin >> qperp;
-
- calc_disp(dyntype, Gx,Gy,Gz, Bx,By,Bz, Px,Py,Pz,
- qperpx,qperpy,qperpz, qperp,
- "dyn.dat", alongqpara==0,
- T, B,
- qrange, qdelta);
-
- return 0;
-}
diff --git a/src/calc/heliphase.cpp b/src/calc/heliphase.cpp
deleted file mode 100644
index a5a0a67..0000000
--- a/src/calc/heliphase.cpp
+++ /dev/null
@@ -1,30 +0,0 @@
-/**
- * Calculates the helimagnetic phase diagram
- * @author tweber@ill.fr
- * @date aug-18
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#include "core/heli.h"
-
-
-int main()
-{
- using t_real = double;
- using t_cplx = std::complex<t_real>;
- const auto j = t_cplx(0,1);
-
- Heli<t_real, t_cplx, 4> heli;
- std::vector<ublas::vector<t_cplx>> fourier{
- tl2::make_vec<ublas::vector<t_cplx>>({0, 0, 0.1}),
- // helical order => Re{M} perp. Im{M}
- tl2::make_vec<ublas::vector<t_cplx>>({1.+j, 1.-j, 0}) / std::sqrt(2),
- };
-
- heli.SetFourier(fourier);
- heli.SetT(-100);
- heli.SetB(0);
- heli.SaveStates("heli.dat", 4, 0, 0, 0, 0);
-
- return 0;
-}
diff --git a/src/calc/skx_gs.cpp b/src/calc/skx_gs.cpp
deleted file mode 100644
index d5ff446..0000000
--- a/src/calc/skx_gs.cpp
+++ /dev/null
@@ -1,185 +0,0 @@
-/**
- * Calculate the ground state fourier components and free energy
- * @author tweber@ill.fr
- * @date aug-18
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#include "core/skx.h"
-#include <fstream>
-
-#define ORDER DEF_SKX_ORDER
-
-using t_real = double;
-using t_cplx = std::complex<t_real>;
-const auto g_j = t_cplx(0,1);
-
-
-// initial values
-// x and y values are the imaginary components
-// z values are the real components
-const std::vector<t_real> _allcomps =
-{{
- // order 4
- /*0, 0, 10.638285,
- -5.8392782, 0, -5.4227098,
- 0.12089063, 0, 0.091292595,
- -0.60718579, 0, -0.62415054,
- 0.10321864, 0, 0.083943453,
- 0.22050946, 0, 0.19010913,
- 0.2190992, 0, 0.19099113,
- -0.011116841, 0, -0.011740137,
- -0.0093003923, 0, -0.0039557086,
- 0.0093273084, 0, 0.010642608,
- -0.0078914558, 0, -0.0028256613,*/
-
- // order 6
- /*0, 0, 10.634562,
- -5.8326696, 0, -5.430272,
- 0.11909555, 0, 0.094198481,
- -0.60621163, 0, -0.62680963,
- 0.09945007, 0, 0.097686206,
- 0.21366494, 0, 0.2020002,
- 0.21365584, 0, 0.20195807,
- -0.014657212, 0, -0.01423701,
- -0.0047616111, 0, -0.0036271944,
- 0.0077666973, 0, 0.0094190521,
- -0.0048230381, 0, -0.0036992597,
- -0.0010215919, 0, -0.0007992665,
- -0.0061202419, 0, -0.0056528426,
- -0.01028418, 0, -0.0090655584,
- -0.010317291, 0, -0.0090892116,
- -0.0061485945, 0, -0.0056819393,
- 0.00052841818, 0, 0.00054131353,
- 0.0010331349, 0, 0.00082690631,
- 0.00081511897, 0, 0.00057051837,
- 0.00060100143, 0, 0.00034822049,
- 0.00080777361, 0, 0.0005647798,
- 0.0010332427, 0, 0.00082431208,*/
-
- // order 8
- /*0, 0, 10.638316,
- -5.8320235, 0, -5.4295657,
- 0.11855499, 0, 0.093494094,
- -0.60754754, 0, -0.62734508,
- 0.098777502, 0, 0.097705375,
- 0.21357743, 0, 0.20213531,
- 0.21357552, 0, 0.20213106,
- -0.014925466, 0, -0.014235675,
- -0.0047648918, 0, -0.0036336754,
- 0.0080509486, 0, 0.0095209452,
- -0.0047710577, 0, -0.0036412848,
- -0.00093943165, 0, -0.00099042964,
- -0.0060178406, 0, -0.0060167827,
- -0.0097664982, 0, -0.0097078794,
- -0.0097684108, 0, -0.0097127543,
- -0.006021222, 0, -0.0060215456,
- 0.00057839383, 0, 0.00059687257,
- 0.00091094805, 0, 0.00093109298,
- 0.00065394655, 0, 0.00063099335,
- 0.00037647447, 0, 0.00034809709,
- 0.0006530733, 0, 0.0006283773,
- 0.00091048399, 0, 0.0009298855,
- -8.383803e-05, 0, -7.4635357e-05,
- 7.1262628e-05, 0, 7.9723568e-05,
- 0.0003774191, 0, 0.00035035634,
- 0.00060968318, 0, 0.00055401694,
- 0.00061024156, 0, 0.00055387028,
- 0.0003778197, 0, 0.00035080693,
- 7.1661474e-05, 0, 8.0687058e-05,
- -2.5559635e-05, 0, -2.2305987e-05,
- -6.4148615e-05, 0, -5.5596784e-05,
- -9.2720594e-05, 0, -7.508628e-05,
- -8.8700786e-05, 0, -6.8695235e-05,
- -7.855418e-05, 0, -5.9074914e-05,
- -8.8137649e-05, 0, -6.8259179e-05,
- -9.1478548e-05, 0, -7.3951933e-05,
- -6.3521218e-05, 0, -5.4970986e-05,*/
-
- // order 9
- 0, 0, 10.643255,
- -5.8312418, 0, -5.4285579,
- 0.11772572, 0, 0.092512566,
- -0.60875193, 0, -0.6283387,
- 0.099009681, 0, 0.097708898,
- 0.21345703, 0, 0.20197682,
- 0.21345707, 0, 0.20197861,
- -0.014776251, 0, -0.014232656,
- -0.0046309613, 0, -0.0035734788,
- 0.0081490635, 0, 0.0096266455,
- -0.0046284703, 0, -0.003570242,
- -0.00092182745, 0, -0.00095170469,
- -0.0060308055, 0, -0.0059615124,
- -0.0098962887, 0, -0.0095705205,
- -0.009896991, 0, -0.0095664551,
- -0.006029697, 0, -0.0059582546,
- 0.00066330707, 0, 0.00055053732,
- 0.0010092654, 0, 0.0008918821,
- 0.00071552573, 0, 0.00060852518,
- 0.00046383518, 0, 0.00034571667,
- 0.00071574625, 0, 0.0006113456,
- 0.0010099268, 0, 0.00089342825,
- 5.3331775e-05, 0, -0.00010249534,
- 0.00015539651, 0, 2.5296191e-06,
- 0.00030612867, 0, 0.00024823487,
- 0.00044277835, 0, 0.00042931965,
- 0.00044175346, 0, 0.00042949054,
- 0.00030547693, 0, 0.00024768997,
- 0.00015489058, 0, 1.3435546e-06,
- 6.9271166e-05, 0, 5.7206645e-05,
- -2.7366626e-05, 0, -2.1078224e-05,
- -3.9356e-05, 0, -5.7759202e-05,
- -3.8566738e-05, 0, -4.1978018e-05,
- -2.7839047e-05, 0, -2.8671778e-05,
- -3.9337396e-05, 0, -4.2705009e-05,
- -4.1111101e-05, 0, -5.945824e-05,
- -2.8210339e-05, 0, -2.1934814e-05,
- 5.8917296e-05, 0, 3.8878991e-05,
- 7.5328415e-05, 0, 5.9727815e-05,
- 6.4888601e-05, 0, 3.2477254e-05,
- 1.0543302e-05, 0, -5.3633989e-06,
- -1.8523783e-05, 0, -2.8534835e-05,
- -1.8660178e-05, 0, -2.8629368e-05,
- 1.0359175e-05, 0, -5.4897177e-06,
- 6.4840587e-05, 0, 3.245778e-05,
- 7.5336564e-05, 0, 5.9731305e-05,
-}};
-
-
-int main()
-{
- Skx<t_real, t_cplx, ORDER> skx;
- skx.SetDebug(1);
-
- std::vector<ublas::vector<t_cplx>> fourier{
-// tl2::make_vec<ublas::vector<t_cplx>>({ 0., 0., 10. }),
-// tl2::make_vec<ublas::vector<t_cplx>>({ -5., 0., -5. }),
- };
-
- fourier.reserve(_allcomps.size()/3);
-
- for(std::size_t comp=0; comp<_allcomps.size(); comp+=3)
- fourier.push_back(tl2::make_vec<ublas::vector<t_cplx>>({_allcomps[comp], _allcomps[comp+1], _allcomps[comp+2]}));
-
-
- skx.SetFourier(fourier);
- skx.SetT(-1000.);
- skx.SetB(25.);
-
- std::cout << "Bc2 = " << get_bc2<t_real>(-1000.) << std::endl;
-
- std::cout.precision(8);
- std::cout << "Order: " << ORDER << std::endl;
- std::cout << "F_start = " << skx.F() << std::endl;
- bool ok = skx.minimise(ORDER, 0,1,0, 1,1,1);
- std::cout << "F_min = " << skx.F() << " (ok: " << std::boolalpha << ok << ")" << std::endl;
-
- std::cout << "\nFourier components:\n";
- for(const auto& fourier : skx.GetFourier())
- {
- std::cout << fourier[0].real() << ", " << fourier[1].real() << ", " << fourier[2].real() << ", ";
- std::cout << std::endl;
- }
-
- return 0;
-}
diff --git a/src/calc/weight.cpp b/src/calc/weight.cpp
deleted file mode 100644
index d84eff6..0000000
--- a/src/calc/weight.cpp
+++ /dev/null
@@ -1,165 +0,0 @@
-/**
- * Calculates the weights/energies at a specific q
- * @author Tobias Weber <tweber@ill.fr>
- * @date sep-19
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#include "core/heli.h"
-#include "core/skx.h"
-#include "core/fp.h"
-
-#include <fstream>
-#include <memory>
-
-using t_real = double;
-using t_cplx = std::complex<t_real>;
-using t_vec = ublas::vector<t_real>;
-const auto g_j = t_cplx(0,1);
-
-#include "core/skx_default_gs.cxx"
-
-
-void calc_disp(char dyntype,
- t_real Gx, t_real Gy, t_real Gz,
- t_real Bx, t_real By, t_real Bz,
- t_real Px, t_real Py, t_real Pz,
- t_real Qx, t_real Qy, t_real Qz,
- int iProj=1, t_real T=-1., t_real B=-1.)
-{
- std::shared_ptr<MagDynamics<t_real, t_cplx>> dyn;
-
- if(dyntype == 's')
- {
- auto skx = std::make_shared<Skx<t_real, t_cplx, DEF_SKX_ORDER>>();
-
- std::vector<ublas::vector<t_cplx>> fourier_skx;
- fourier_skx.reserve(_skxgs_allcomps.size()/3);
-
- for(std::size_t comp=0; comp<_skxgs_allcomps.size(); comp+=3)
- fourier_skx.push_back(tl2::make_vec<ublas::vector<t_cplx>>({_skxgs_allcomps[comp], _skxgs_allcomps[comp+1], _skxgs_allcomps[comp+2]}));
-
- skx->SetFourier(fourier_skx);
- skx->SetProjNeutron(iProj!=0);
- skx->SetT(-1000.);
- skx->SetB(25.); // BC2 = 40.3425
- skx->GenFullFourier();
- skx->SetFilterZeroWeight(1);
- skx->SetCoords(Bx,By,Bz, Px,Py,Pz);
- skx->SetG(Gx,Gy,Gz);
-
- dyn = skx;
- }
- else if(dyntype == 'h')
- {
- auto heli = std::make_shared<Heli<t_real, t_cplx, DEF_HELI_ORDER>>();
-
- heli->SetProjNeutron(iProj!=0);
- heli->SetT(T);
- heli->SetB(B);
- heli->SetFilterZeroWeight(1);
- heli->SetCoords(Bx,By,Bz);
- heli->SetG(Gx,Gy,Gz);
-
- dyn = heli;
- }
- else if(dyntype == 'f')
- {
- auto fp = std::make_shared<FP<t_real, t_cplx>>();
-
- fp->SetT(T);
- fp->SetB(B);
- fp->SetCoords(Bx,By,Bz);
- fp->SetG(Gx,Gy,Gz);
-
- dyn = fp;
- }
- else
- {
- std::cerr << "Unknown dynamics type selected." << std::endl;
- return;
- }
-
-
- t_real Erange = -1.; // negative: disable range
- auto [Es, wsUnpol, wsSF1, wsSF2, wsNSF] = dyn->GetDisp(Qx, Qy, Qz, -Erange, Erange);
-
- std::cout << "# Magnetic phase: " << dyntype << "\n"
- << std::setw(15) << std::left << "# No." << " "
- << std::setw(15) << std::left << "E (meV)" << " "
- << std::setw(15) << std::left << "w_total" << " "
- << std::setw(15) << std::left << "w_SF1" << " "
- << std::setw(15) << std::left << "w_SF2" << " "
- << std::setw(15) << std::left << "w_NSF" << "\n";
-
- for(std::size_t i=0; i<Es.size(); ++i)
- {
- std::cout
- << std::setw(15) << std::left << (i+1) << " "
- << std::setw(15) << std::left /*<< std::scientific*/ << Es[i] << " "
- << std::setw(15) << std::left /*<< std::scientific*/ << wsUnpol[i] << " "
- << std::setw(15) << std::left /*<< std::scientific*/ << wsSF1[i] << " "
- << std::setw(15) << std::left /*<< std::scientific*/ << wsSF2[i] << " "
- << std::setw(15) << std::left /*<< std::scientific*/ << wsNSF[i] << "\n";
- }
-}
-
-
-int main()
-{
- std::cout
- << "--------------------------------------------------------------------------------\n"
- << "\tDynamical structure factor calculation tool,\n\t\tT. Weber <tweber@ill.fr>, September 2019.\n"
- << "--------------------------------------------------------------------------------\n\n";
-
- while(1)
- {
- t_real Gx = 1., Gy = 1., Gz = 0.;
- t_real Bx = 1., By = 1., Bz = 0.;
- t_real Qx = 1.068, Qy = 0.932, Qz = 0.;
- t_real Px = 1., Py = -1., Pz = 0.;
- t_real B = 0.17, T = 20.;
- int proj = 1;
- char dyntype = 'h';
-
- std::cout << "Helimagnon [h], skyrmion [s] or field-polarised [f] dynamics: ";
- std::cin >> dyntype; dyntype = std::tolower(dyntype);
- std::cout << "G = ";
- std::cin >> Gx >> Gy >> Gz;
- std::cout << "Q = ";
- std::cin >> Qx >> Qy >> Qz;
- std::cout << "B = ";
- std::cin >> Bx >> By >> Bz;
- if(dyntype == 'h' || dyntype == 'f')
- {
- std::cout << "|B| = ";
- std::cin >> B;
- std::cout << "T = ";
- std::cin >> T;
- }
- else if(dyntype == 's')
- {
- std::cout << "pinning = ";
- std::cin >> Px >> Py >> Pz;
- }
- std::cout << "Q projector [0/1]: ";
- std::cin >> proj;
-
- std::cout << "\n";
- if(dyntype == 'h' || dyntype == 'f')
- {
- std::cout << "# T = " << T << "\n";
- std::cout << "# |B| = " << B << "\n";
- }
- std::cout << "# B = (" << Bx << ", " << By << ", " << Bz << ")\n";
- std::cout << "# G = (" << Gx << ", " << Gy << ", " << Gz << ")\n";
- std::cout << "# Q = (" << Qx << ", " << Qy << ", " << Qz << ")\n";
- std::cout << "# q = (" << Qx-Gx << ", " << Qy-Gy << ", " << Qz-Gz << ")\n";
- std::cout << "# Q_proj = " << proj << "\n";
-
- calc_disp(dyntype, Gx,Gy,Gz, Bx,By,Bz, Px,Py,Pz, Qx,Qy,Qz, proj, T, B);
- std::cout << std::endl;
- }
-
- return 0;
-}
diff --git a/src/calc/weight_sum.cpp b/src/calc/weight_sum.cpp
deleted file mode 100644
index 03eb896..0000000
--- a/src/calc/weight_sum.cpp
+++ /dev/null
@@ -1,226 +0,0 @@
-/**
- * Calculates the integrated weights/energies for setup 3
- * @author Tobias Weber <tweber@ill.fr>
- * @date jun-20
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#include "core/heli.h"
-#include "core/skx.h"
-#include "tlibs2/libs/phys.h"
-
-#include <fstream>
-
-#include <boost/histogram.hpp>
-namespace hist = boost::histogram;
-
-using t_real = double;
-using t_cplx = std::complex<t_real>;
-using t_vec = ublas::vector<t_real>;
-const auto g_j = t_cplx(0,1);
-
-#include "core/skx_default_gs.cxx"
-
-
-const t_real g_T = 28.5;
-const t_real g_eps = 1e-5;
-
-#define E_BINS 200
-#define NUM_ANGLES 512
-#define COL_SIZE 15
-
-
-void calc_disp(
- t_real Gx, t_real Gy, t_real Gz,
- t_real Bx, t_real By, t_real Bz,
- t_real Px, t_real Py, t_real Pz,
- t_real q, int iProj = 1)
-{
- Skx<t_real, t_cplx, DEF_SKX_ORDER> skx;
- Heli<t_real, t_cplx, DEF_HELI_ORDER> heli;
-
- std::vector<ublas::vector<t_cplx>> fourier_skx;
- fourier_skx.reserve(_skxgs_allcomps.size()/3);
-
- for(std::size_t comp=0; comp<_skxgs_allcomps.size(); comp+=3)
- fourier_skx.push_back(tl2::make_vec<ublas::vector<t_cplx>>({_skxgs_allcomps[comp], _skxgs_allcomps[comp+1], _skxgs_allcomps[comp+2]}));
-
- skx.SetFourier(fourier_skx);
-
- skx.SetProjNeutron(iProj!=0);
- heli.SetProjNeutron(iProj!=0);
-
- skx.SetT(-1000.);
- heli.SetT(g_T);
-
- skx.SetB(25.);
- heli.SetB(0.17);
-
- skx.GenFullFourier();
-
- skx.SetFilterZeroWeight(1);
- heli.SetFilterZeroWeight(1);
-
- skx.SetWeightEps(1e-6);
- heli.SetWeightEps(1e-6);
-
- skx.SetCoords(Bx,By,Bz, Px,Py,Pz);
- heli.SetCoords(Bx,By,Bz);
-
- skx.SetG(Gx, Gy, Gz);
- heli.SetG(Gx, Gy, Gz);
-
-
- t_real Erange = 0.1;
-
- t_real angle_offs = 90/180.*M_PI;
- t_real angle_begin = -135/180.*M_PI + angle_offs;
- t_real angle_end = 135/180.*M_PI + angle_offs;
- t_real angle_delta = 2*M_PI/t_real(NUM_ANGLES);
-
- auto histWeightsNSF = hist::make_histogram(hist::axis::regular<t_real>(E_BINS, -Erange, Erange, "E"));
- auto histWeightsSF = hist::make_histogram(hist::axis::regular<t_real>(E_BINS, -Erange, Erange, "E"));
-
- auto histWeightsHeliNSF = hist::make_histogram(hist::axis::regular<t_real>(E_BINS, -Erange, Erange, "E"));
- auto histWeightsHeliSF = hist::make_histogram(hist::axis::regular<t_real>(E_BINS, -Erange, Erange, "E"));
-
-
- std::ofstream ofstr_raw("weightsum_skx.dat");
- std::ofstream ofstr_raw_heli("weightsum_heli.dat");
-
- // write file header
- for(std::ostream* ostr : {&ofstr_raw, &ofstr_raw_heli})
- {
- ostr->precision(8);
- (*ostr) << std::left << std::setw(COL_SIZE) << "# angle"
- << " " << std::left << std::setw(COL_SIZE) << "qh"
- << " " << std::left << std::setw(COL_SIZE) << "qk"
- << " " << std::left << std::setw(COL_SIZE) << "ql"
- << " " << std::left << std::setw(COL_SIZE) << "E"
- << " " << std::left << std::setw(COL_SIZE) << "wSF1"
- << " " << std::left << std::setw(COL_SIZE) << "wSF2"
- << " " << std::left << std::setw(COL_SIZE) << "wNSF"
- << "\n";
- }
-
-
- for(t_real angle=angle_begin; angle<angle_end; angle+=angle_delta)
- {
- t_real Qx = q * cos(angle) + Gx;
- t_real Qy = q * sin(angle) + Gy;
- t_real Qz = Gz;
-
- std::cout << "# angle: " << angle/M_PI*180. << ", Q = (" << Qx << ", " << Qy << ", " << Qz << ")" << std::endl;
-
- {
- auto [Es, wsUnpol, wsSF1, wsSF2, wsNSF] = skx.GetDisp(Qx, Qy, Qz, -Erange, Erange);
- for(std::size_t i=0; i<Es.size(); ++i)
- {
- histWeightsNSF(Es[i], hist::weight(wsNSF[i]*0.5));
- histWeightsSF(Es[i], hist::weight(wsSF1[i]));
-
- ofstr_raw << std::left << std::setw(COL_SIZE) << angle
- << " " << std::left << std::setw(COL_SIZE) << (Qx-Gx)
- << " " << std::left << std::setw(COL_SIZE) << (Qy-Gy)
- << " " << std::left << std::setw(COL_SIZE) << (Qz-Gz)
- << " " << std::left << std::setw(COL_SIZE) << Es[i]
- << " " << std::left << std::setw(COL_SIZE) << wsSF1[i]
- << " " << std::left << std::setw(COL_SIZE) << wsSF2[i]
- << " " << std::left << std::setw(COL_SIZE) << wsNSF[i]
- << std::endl;
- }
- }
-
- {
- auto [EsH, wsUnpolH, wsSF1H, wsSF2H, wsNSFH] = heli.GetDisp(Qx, Qy, Qz, -Erange, Erange);
- for(std::size_t i=0; i<EsH.size(); ++i)
- {
- histWeightsHeliNSF(EsH[i], hist::weight(wsNSFH[i]*0.5));
- histWeightsHeliSF(EsH[i], hist::weight(wsSF1H[i]));
-
- ofstr_raw_heli << std::left << std::setw(COL_SIZE) << angle
- << " " << std::left << std::setw(COL_SIZE) << (Qx-Gx)
- << " " << std::left << std::setw(COL_SIZE) << (Qy-Gy)
- << " " << std::left << std::setw(COL_SIZE) << (Qz-Gz)
- << " " << std::left << std::setw(COL_SIZE) << EsH[i]
- << " " << std::left << std::setw(COL_SIZE) << wsSF1H[i]
- << " " << std::left << std::setw(COL_SIZE) << wsSF2H[i]
- << " " << std::left << std::setw(COL_SIZE) << wsNSFH[i]
- << std::endl;
- }
- }
- }
-
-
- std::ofstream ofstrBinned("../data/weightbin.dat");
- ofstrBinned.precision(8);
-
- ofstrBinned << std::left << std::setw(COL_SIZE) << "# E (meV)"
- << " " << std::left << std::setw(COL_SIZE) << "bose"
- << " " << std::left << std::setw(COL_SIZE) << "weight_skx_sf"
- << " " << std::left << std::setw(COL_SIZE) << "weight_skx_nsf"
- << " " << std::left << std::setw(COL_SIZE) << "weight_skx_sum"
- << " " << std::left << std::setw(COL_SIZE) << "weight_skx_sum_bose"
- << " " << std::left << std::setw(COL_SIZE) << "weight_heli_sf"
- << " " << std::left << std::setw(COL_SIZE) << "weight_heli_nsf"
- << " " << std::left << std::setw(COL_SIZE) << "weight_heli_sum"
- << " " << std::left << std::setw(COL_SIZE) << "weight_heli_sum_bose"
- << "\n";
-
- auto iterHeliNSF = boost::histogram::indexed(histWeightsHeliNSF).begin();
- auto iterHeliSF = boost::histogram::indexed(histWeightsHeliSF).begin();
- auto iterNSF = boost::histogram::indexed(histWeightsNSF).begin();
- for(const auto& val : boost::histogram::indexed(histWeightsSF))
- {
- t_real E = val.bin().lower() + 0.5*(val.bin().upper() - val.bin().lower());
- t_real w = *val / t_real{E_BINS};
-
- t_real E_nsf = iterNSF->bin().lower() + 0.5*(iterNSF->bin().upper() - iterNSF->bin().lower());
- t_real w_nsf = **iterNSF / t_real{E_BINS};
-
- t_real E_h = iterHeliSF->bin().lower() + 0.5*(iterHeliSF->bin().upper() - iterHeliSF->bin().lower());
- t_real w_h = **iterHeliSF / t_real{E_BINS};
-
- t_real E_h_nsf = iterHeliNSF->bin().lower() + 0.5*(iterHeliNSF->bin().upper() - iterHeliNSF->bin().lower());
- t_real w_h_nsf = **iterHeliNSF / t_real{E_BINS};
-
- if(!tl2::float_equal<t_real>(E, E_h, g_eps) || !tl2::float_equal<t_real>(E, E_h_nsf, g_eps) ||
- !tl2::float_equal<t_real>(E, E_nsf, g_eps))
-
- {
- std::cerr << "Energy binning mismatch: " << E << " != " << E_h << std::endl;
- break;
- }
-
- t_real bose = tl2::bose(E, g_T);
-
- ofstrBinned << std::left << std::setw(COL_SIZE) << E
- << " " << std::left << std::setw(COL_SIZE) << bose
- << " " << std::left << std::setw(COL_SIZE) << w
- << " " << std::left << std::setw(COL_SIZE) << w_nsf
- << " " << std::left << std::setw(COL_SIZE) << (w + w_nsf)
- << " " << std::left << std::setw(COL_SIZE) << (w + w_nsf)*bose
- << " " << std::left << std::setw(COL_SIZE) << w_h
- << " " << std::left << std::setw(COL_SIZE) << w_h_nsf
- << " " << std::left << std::setw(COL_SIZE) << (w_h + w_h_nsf)
- << " " << std::left << std::setw(COL_SIZE) << (w_h + w_h_nsf)*bose
- << std::endl;
-
- ++iterHeliSF;
- ++iterHeliNSF;
- ++iterNSF;
- }
-}
-
-
-int main()
-{
- // equivalent to the used setup around (000)
- t_real Gx = 1., Gy = 1., Gz = 0.;
- t_real Bx = 0., By = 0., Bz = 1.;
- t_real q = 0.0123;
- int proj = 1;
-
- calc_disp(Gx,Gy,Gz, Bx,By,Bz, 1,0,0, q, proj);
- return 0;
-}
diff --git a/src/core/constants.h b/src/core/constants.h
deleted file mode 100644
index 84ae284..0000000
--- a/src/core/constants.h
+++ /dev/null
@@ -1,96 +0,0 @@
-/**
- * MnSi constants
- * @author Tobias Weber <tweber@ill.fr>
- * @date jul-18
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#ifndef __MNSI_CONSTS_H__
-#define __MNSI_CONSTS_H__
-
-#include "tlibs2/libs/math17.h"
-#include "tlibs2/libs/units.h"
-
-
-// constants
-template<class t_real = double> constexpr t_real g_pi = tl2::pi<t_real>;
-template<class t_real = double> constexpr t_real g_kB
- = static_cast<t_real>(tl2::kB<t_real>/tl2::meV<t_real>*tl2::kelvin<t_real>);
-template<class t_real = double> constexpr t_real g_muB
- = static_cast<t_real>(tl2::muB<t_real>/tl2::meV<t_real>*tl2::tesla<t_real>);
-template<class t_real = double> constexpr t_real g_chi = 0.34;
-template<class t_real = double> constexpr t_real g_hoc = -0.0073;
-template<class t_real = double> constexpr t_real g_a = 4.558;
-template<class t_real = double> constexpr t_real g_g = 2.;
-template<class t_real = double> constexpr t_real g_kh_A_29K = 0.039;
-template<class t_real = double> constexpr t_real g_kh_rlu_29K = g_kh_A_29K<t_real> / (t_real(2.)*g_pi<t_real> / g_a<t_real>);
-
-
-/**
- * helix pitch in inverse angstroms
- */
-template<class t_real = double>
-constexpr t_real g_kh_A(t_real T)
-{
- constexpr t_real T1 = 20.;
- constexpr t_real T2 = 29.;
-
- constexpr t_real A1 = 0.036; // at 20 K
- constexpr t_real A2 = g_kh_A_29K<t_real>; // at 29 K
-
- constexpr t_real m = (A2-A1) / (T2-T1);
- constexpr t_real t = A2 - m*T2;
-
- return m*T + t;
-}
-
-
-/**
- * helix pitch in rlu
- */
-template<class t_real = double>
-constexpr t_real g_kh_rlu(t_real T)
-{
- return g_kh_A<t_real>(T) / (t_real(2.)*g_pi<t_real> / g_a<t_real>);
-}
-
-
-/**
- * get upper critical field
- * - theoretical values calculated with "heli.dat_Bc2.gpl" script generated by heliphase.cpp
- * - experimental values from A. Bauer, 2015
- */
-template<class t_real=double>
-t_real get_bc2(t_real T, bool use_theo_units=1)
-{
- if(use_theo_units)
- {
- const t_real amp = 1.62567; // scaling
- const t_real ex = 0.46491; // critical exponent
-
- if(T >= 0.) return 0.;
- return amp * std::pow(-T, ex);
- }
- else
- {
- t_real p1[] = { 5.63617388e-01, 5.52570801e-02, 2.20736277e+01, 7.39287474e-02, -5.32767610e-04 };
- t_real p2[] = { 0.07075453, -0.08217821, 30.00000534, 9.19469462, 0.38951838 };
-
- const t_real *p = (T<=20 ? p1 : p2);
- t_real Tc = p[2];
- t_real tau = (Tc-T) / Tc;
-
- if(T >= Tc) return 0.;
- return p[0]*std::pow(tau, p[1]) * (1. + p[3]*std::pow(tau, p[4]));
- }
-}
-
-
-#define G_CHI g_chi<t_real>
-#define G_HOC g_hoc<t_real>
-#define G_G g_g<t_real>
-#define G_MUB g_muB<t_real>
-#define G_KH_RLU_29K g_kh_rlu_29K<t_real>
-
-
-#endif
diff --git a/src/core/fp.cpp b/src/core/fp.cpp
deleted file mode 100644
index 331de26..0000000
--- a/src/core/fp.cpp
+++ /dev/null
@@ -1,181 +0,0 @@
-/**
- * Field-polarised phase
- * @author Tobias Weber <tweber@ill.fr>
- * @date dec-16
- * @desc This file implements the theoretical magnon model by M. Garst and J. Waizner, see:
- * - https://doi.org/10.1088/1361-6463/aa7573
- * - https://kups.ub.uni-koeln.de/7937/
- * @desc This file is based on:
- * - the Mathematica description of the field-polarised magnon model by M. Garst, 2016.
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#include "fp.h"
-
-
-// instantiation
-template class FP<double, std::complex<double>>;
-#ifdef __HACK_FULL_INST__
- template FP<double, std::complex<double>>::FP();
- template void FP<double, std::complex<double>>::SetG(double, double, double);
- template void FP<double, std::complex<double>>::SetCoords(double, double, double);
-#endif
-
-
-
-template<class t_real, class t_cplx>
-FP<t_real, t_cplx>::FP()
-{}
-
-
-/**
- * set lattice vector and orthogonal projector
- */
-template<class t_real, class t_cplx>
-void FP<t_real, t_cplx>::SetG(t_real h, t_real k, t_real l)
-{
- m_Grlu = tl2::make_vec<t_vec>({h,k,l});
-
- t_vec G = m_Grlu / tl2::veclen(m_Grlu);
- G = tl2::quat_vec_prod(m_rotCoord, G);
- m_projNeutron = tl2::unit_m<t_mat>(3) - tl2::outer<t_vec, t_mat>(G, G);
-}
-
-
-/**
- * rotates field to internal [001] convention
- */
-template<class t_real, class t_cplx>
-void FP<t_real, t_cplx>::SetCoords(t_real Bx, t_real By, t_real Bz)
-{
- t_vec B = tl2::make_vec<t_vec>( {Bx, By, Bz} );
- t_quat quatB = tl2::rotation_quat(B, tl2::make_vec<t_vec>( {0, 0, 1} ));
-
- m_rotCoord = quatB;
-}
-
-
-/**
- * query the dispersion
- */
-template<class t_real, class t_cplx>
-std::tuple<std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>>
-FP<t_real, t_cplx>::GetDisp(t_real h, t_real k, t_real l, t_real minE, t_real maxE) const
-{
- constexpr auto imag = t_cplx(0,1);
- const auto ident2 = tl2::unit_m<t_mat_cplx>(2);
- const auto sigma = tl2::get_spin_matrices<ublas::matrix, std::vector, t_real>();
- const t_real dh0_shift = 0.007523;
- const t_real dh0 = dh0_shift + (m_B - m_Bc2) / m_Bc2; // for given g_hoc
-
-
- t_vec Qvec = tl2::make_vec<t_vec>({ h, k, l });
- t_vec qvec = tl2::quat_vec_prod(m_rotCoord, Qvec) - tl2::quat_vec_prod(m_rotCoord, m_Grlu);
- qvec /= g_kh_rlu<t_real>(m_T);
- qvec = -qvec;
-
- t_real q=0., theta=0., phi=0.;
- std::tie(q, phi, theta) = tl2::cart_to_sph(qvec[0], qvec[1], qvec[2]);
-
-
- // eigensystems
- auto get_evecs = [&theta, &qvec, &q, &dh0, &imag, &sigma, &ident2](t_real phi) -> auto
- {
- // Hamiltonian
- t_cplx qp = qvec[0] + imag*qvec[1];
- t_cplx qm = qvec[0] - imag*qvec[1];
-
- t_mat_cplx H = (q*q + g_hoc<t_real>*q*q*q*q + 1. + dh0) * ident2;
- H += 2.*sigma[2] * qvec[2];
- H += g_chi<t_real>/(2.*q*q) * tl2::make_mat<t_mat_cplx>(
- {{ qm*qp, qm*qm },
- { qp*qp, qm*qp }
- });
-
-
- // eigenvectors and -values
- std::vector<t_vec_cplx> evecs;
- std::vector<t_real> evals;
- bool ev_ok = tl2::eigenvecsel_herm<t_real>(tl2::prod_mm(sigma[2], H), evecs, evals, true);
- /*if(ev_ok)
- {
- t_cplx norm = tl2::mat_elem(evecs[0], sigma[2], evecs[1]);
- evecs[0] /= norm;
- evecs[1] /= norm;
- }*/
-
- return std::make_tuple(ev_ok, evals, evecs);
- };
-
-
- /*
- * spectral weights
- *
- * transition probability to flip spin given by propagator, i.e. Green's function:
- * G(x, x') = sum{j} |v_j'> <v_j| / E_j
- * with eigenvectors |v_j> and eigenvalues E_j.
- */
- auto get_weights = [&imag](const t_vec_cplx& evec_phi, const t_vec_cplx& evec_mphi, const t_mat_cplx& projNeutron, bool bProjNeutron) -> auto
- {
- t_mat_cplx kernel = tl2::outer<t_vec_cplx, t_mat_cplx>(evec_phi, evec_mphi);
- t_mat_cplx weight(3,3);
-
- const t_vec_cplx x = tl2::make_vec<t_vec_cplx>({ 1, 0, 0 });
- const t_vec_cplx y = tl2::make_vec<t_vec_cplx>({ 0, 1, 0 });
-
- for(int i=0; i<3; ++i)
- {
- t_vec_cplx veci = tl2::make_vec<t_vec_cplx>({ x[i]-imag*y[i], x[i]+imag*y[i] });
-
- for(int j=0; j<3; ++j)
- {
- t_vec_cplx vecj = tl2::make_vec<t_vec_cplx>({ x[j]-imag*y[j], x[j]+imag*y[j] });
- weight(i,j) = tl2::mat_elem(veci, kernel, vecj);
- }
- }
-
-
- // magnetic neutron scattering projections
- t_mat_cplx neutron_weight = weight;
- if(bProjNeutron)
- {
- neutron_weight = tl2::prod_mm(weight, projNeutron);
- neutron_weight = tl2::prod_mm(projNeutron, neutron_weight);
- }
-
-
- // polarisation projections
- t_mat_cplx matSF1 = tl2::prod_mm(get_chiralpol<t_mat_cplx>(1), neutron_weight);
- t_mat_cplx matSF2 = tl2::prod_mm(get_chiralpol<t_mat_cplx>(2), neutron_weight);
- t_mat_cplx matNSF = tl2::prod_mm(get_chiralpol<t_mat_cplx>(3), neutron_weight);
-
- constexpr const t_real weightScale = 0.5;
- t_cplx wAll = tl2::trace(neutron_weight) * weightScale;
- t_cplx wSF1 = tl2::trace(matSF1) * weightScale;
- t_cplx wSF2 = tl2::trace(matSF2) * weightScale;
- t_cplx wNSF = tl2::trace(matNSF) * weightScale;
-
- return std::make_tuple(wAll, wSF1, wSF2, wNSF);
- };
-
-
- auto [ok_phi, evals_phi, evecs_phi] = get_evecs(phi);
- auto [ok_mphi, evals_mphi, evecs_mphi] = get_evecs(-phi);
- if(!ok_phi || !ok_mphi)
- {
- std::vector<t_real> empty;
- return std::make_tuple(empty, empty, empty, empty, empty);
- }
-
- auto [wAll_p, wSF1_p, wSF2_p, wNSF_p] = get_weights(evecs_phi[0], evecs_mphi[0], m_projNeutron, m_bProjNeutron);
- auto [wAll_m, wSF1_m, wSF2_m, wNSF_m] = get_weights(-evecs_phi[1], evecs_mphi[1], m_projNeutron, m_bProjNeutron);
-
-
- std::vector<t_real> Es = { evals_phi[0]*g_g<t_real>*g_muB<t_real>*m_Bc2, evals_phi[1]*g_g<t_real>*g_muB<t_real>*m_Bc2 };
- std::vector<t_real> ws_all = { std::abs(wAll_p.real()), std::abs(wAll_m.real()) };
- std::vector<t_real> ws_SF1 = { std::abs(wSF1_p.real()), std::abs(wSF1_m.real()) };
- std::vector<t_real> ws_SF2 = { std::abs(wSF2_p.real()), std::abs(wSF2_m.real()) };
- std::vector<t_real> ws_NSF = { std::abs(wNSF_p.real()), std::abs(wNSF_m.real()) };
-
- return std::make_tuple(Es, ws_all, ws_SF1, ws_SF2, ws_NSF);
-}
diff --git a/src/core/fp.h b/src/core/fp.h
deleted file mode 100644
index a4757b8..0000000
--- a/src/core/fp.h
+++ /dev/null
@@ -1,73 +0,0 @@
-/**
- * Field-polarised phase
- * @author Tobias Weber <tweber@ill.fr>
- * @date dec-16, sep-18
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#ifndef __FP_H__
-#define __FP_H__
-
-#include <tuple>
-#include <vector>
-#include <complex>
-
-#include "helper.h"
-#include "magsys.h"
-
-#define USE_LAPACK
-#include "tlibs2/libs/math17.h"
-
-
-/**
- * FP
- */
-template<class t_real=double, class t_cplx = std::complex<t_real>>
-class FP : public MagDynamics<t_real, t_cplx>
-{
-public:
- using t_mat = ublas::matrix<t_real>;
- using t_vec = ublas::vector<t_real>;
- using t_quat = boost::math::quaternion<t_real>;
- using t_mat_cplx = ublas::matrix<t_cplx>;
- using t_vec_cplx = ublas::vector<t_cplx>;
-
-
-public:
- FP();
-
- std::shared_ptr<FP<t_real, t_cplx>> copy() const
- {
- return std::make_shared<FP<t_real, t_cplx>>(*this);
- }
-
- virtual std::shared_ptr<MagDynamics<t_real, t_cplx>> copyCastDyn() const override { return copy(); }
-
-
- virtual void SetB(t_real B) override { m_B = B; }
- virtual void SetT(t_real T) override { m_T = T; m_Bc2 = get_bc2(m_T, false); }
- t_real GetBC2() const { return m_Bc2; }
-
- void SetG(t_real h, t_real k, t_real l);
- void SetCoords(t_real Bx, t_real By, t_real Bz);
- const t_mat_cplx& GetNeutronProjOp() const { return m_projNeutron; }
- void SetProjNeutron(bool b) { m_bProjNeutron = b; }
-
- virtual std::tuple<std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>>
- GetDisp(t_real h, t_real k, t_real l, t_real minE=-1., t_real maxE=-2.) const override;
-
-
-private:
- t_real m_Bc2 = 0.55;
- t_real m_B = 0.7;
- t_real m_T = 20;
-
- t_vec m_Grlu;
- t_quat m_rotCoord;
-
- t_mat_cplx m_projNeutron = tl2::unit_m<t_mat_cplx>(3);
- bool m_bProjNeutron = true;
-};
-
-
-#endif
diff --git a/src/core/heli.cpp b/src/core/heli.cpp
deleted file mode 100644
index b3acf60..0000000
--- a/src/core/heli.cpp
+++ /dev/null
@@ -1,399 +0,0 @@
-/**
- * Free energy in the helimagnetic phase
- * @author tweber@ill.fr
- * @date mid-16, jul-18
- * @desc This file implements the theoretical helimagnon model by M. Garst and J. Waizner, see:
- * - https://doi.org/10.1088/1361-6463/aa7573
- * - https://kups.ub.uni-koeln.de/7937/
- * - https://doi.org/10.1103/PhysRevLett.115.097203
- * @desc This file is based on:
- * - the descriptions and Mathematica implementations of the different helimagnon model versions by M. Garst and J. Waizner, 2014-2018,
- * - the 2015 and 2016 Python implementations by G. Brandl and M. Kugler of the first version of the helimagnon model.
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#include <cmath>
-#include <numeric>
-#include <functional>
-#include <iostream>
-
-#include "heli.h"
-
-
-// instantiation
-#ifdef DEF_HELI_ORDER
- #pragma message("Heli Order: " __TL2_STRCONV(DEF_HELI_ORDER))
- template class Heli<double, std::complex<double>, DEF_HELI_ORDER>;
-
- #ifdef __HACK_FULL_INST__
- template Heli<double, std::complex<double>, DEF_HELI_ORDER>::Heli();
- template void Heli<double, std::complex<double>, DEF_HELI_ORDER>::SetG(double, double, double);
- template void Heli<double, std::complex<double>, DEF_HELI_ORDER>::SetCoords(double, double, double);
- #endif
-#endif
-
-
-/**
- * constructor
- */
-template<class t_real, class t_cplx, int ORDER>
-Heli<t_real, t_cplx, ORDER>::Heli()
-{
- // #1
- std::vector<int> m_idx1;
- m_idx1.resize(ORDER);
- std::iota(m_idx1.begin(), m_idx1.end(), 1);
-
- // #2
- for(int i=0; i<ORDER*2+1; ++i)
- std::copy(m_idx1.begin(), m_idx1.end(), std::back_inserter(m_idx2[0]));
-
- for(int i=-ORDER; i<=ORDER; ++i)
- for(int j=0; j<ORDER; ++j)
- m_idx2[1].push_back(i);
-
- // #3
- for(int i=0; i<ORDER*2+1; ++i)
- {
- std::copy(m_idx2[0].begin(), m_idx2[0].end(), std::back_inserter(m_idx3[0]));
- std::copy(m_idx2[1].begin(), m_idx2[1].end(), std::back_inserter(m_idx3[1]));
- }
-
- for(int i=-ORDER; i<=ORDER; ++i)
- for(int j=0; j<(2*ORDER+1)*ORDER; ++j)
- m_idx3[2].push_back(i);
-
- // reduce #2
- decltype(m_idx2) idx2;
- for(std::size_t i=0; i<m_idx2[0].size(); ++i)
- {
- int val = -m_idx2[0][i] - m_idx2[1][i];
- if(std::abs(val) > ORDER)
- { continue; }
- else
- {
- idx2[0].push_back(m_idx2[0][i]);
- idx2[1].push_back(m_idx2[1][i]);
- idx2[2].push_back(val);
- }
- }
- for(int i=0; i<3; ++i)
- m_idx2[i] = std::move(idx2[i]);
-
- // reduce #3
- decltype(m_idx3) idx3;
- for(std::size_t i=0; i<m_idx3[0].size(); ++i)
- {
- int val = -m_idx3[0][i] - m_idx3[1][i] - m_idx3[2][i];
- if(std::abs(val) > ORDER)
- { continue; }
- else
- {
- idx3[0].push_back(m_idx3[0][i]);
- idx3[1].push_back(m_idx3[1][i]);
- idx3[2].push_back(m_idx3[2][i]);
- idx3[3].push_back(val);
- }
- }
- for(int i=0; i<4; ++i)
- m_idx3[i] = std::move(idx3[i]);
-}
-
-
-/**
- * free energy
- * e.g.: M = (1+i, 1-i, 0), prop_1 = (0, 0, kh), prop_2 = (0, 0, -i*kh)
- *
- * Heisenberg term:
- * - discrete, real: -J_ij S_i * S_j
- * - continuous, real: J (del M(r))^2, Fourier: J k^2 |M(k)|^2
- *
- * DM term:
- * - discrete, real: D_ij * (S_i cross S_j)
- * - continuous, real: D M(r) (del cross M(r))
- */
-template<class t_real, class t_cplx, int ORDER>
-t_real Heli<t_real, t_cplx, ORDER>::F()
-{
- // only z component used
- auto fourier0 = m_fourier[0][2];
-
- auto fourier_no0 = m_fourier;
- fourier_no0.erase(fourier_no0.begin());
-
- // add complex conjugate
- auto fourier_full = m_fourier;
- for(std::size_t i=m_fourier.size()-1; i>=1; --i)
- fourier_full.push_back(tl2::conjugate_vec(m_fourier[i]));
-
-
- std::vector<t_cplx> fourier_full_comp[3];
- for(const auto &vecC : fourier_full)
- {
- auto tup = split_vec3d(vecC);
- fourier_full_comp[0].push_back(std::get<0>(tup));
- fourier_full_comp[1].push_back(std::get<1>(tup));
- fourier_full_comp[2].push_back(std::get<2>(tup));
- }
-
- // lengths
- std::vector<t_real> lens;
- for(std::size_t i=0; i<fourier_no0.size(); ++i)
- {
- const auto [fourier_real, fourier_imag] = tl2::split_cplx_vec<
- ublas::vector<std::complex<t_real>>, ublas::vector<t_real>>(fourier_no0[i]);
-
- t_real len = ublas::inner_prod(fourier_real, fourier_real) +
- ublas::inner_prod(fourier_imag, fourier_imag);
- lens.push_back(2. * len);
- }
-
-
- // free energy
- t_cplx cF = 0;
-
- // dip
- cF += g_chi<t_real>/3. * fourier0*fourier0;
- for(std::size_t i=0; i<fourier_no0.size(); ++i)
- cF += 2. * g_chi<t_real> * std::norm(fourier_no0[i][2]);
-
- // dmi
- for(int i=0; i<ORDER; ++i)
- {
- t_real k = t_real(i+1);
-
- cF += 8. * m_pitch * k *
- ( fourier_no0[i][0].real() * fourier_no0[i][1].imag() -
- fourier_no0[i][1].real() * fourier_no0[i][0].imag() );
- }
-
- // hoc
- for(int i=0; i<ORDER; ++i)
- {
- t_real k = t_real(i+1);
-
- cF += g_hoc<t_real> * m_pitch*m_pitch * k*k*k*k * lens[i];
- }
-
- // phi^2 & phi^4
- cF += (m_T + 1.) * fourier0*fourier0 + fourier0*fourier0*fourier0*fourier0;
- for(int i=0; i<ORDER; ++i)
- {
- t_real k = t_real(i+1);
-
- // Heisenberg
- cF += m_pitch*m_pitch * k*k * lens[i];
-
- cF += (m_T + 1. + fourier0*fourier0) * lens[i];
- }
-
- for(int i=0; i<m_idx2[0].size(); ++i)
- {
- cF += 2. * fourier0 * get_comp(fourier_full_comp[2], m_idx2[0][i]) *
- (
- get_comp(fourier_full_comp[0], m_idx2[1][i]) * get_comp(fourier_full_comp[0], m_idx2[2][i]) +
- get_comp(fourier_full_comp[1], m_idx2[1][i]) * get_comp(fourier_full_comp[1], m_idx2[2][i]) +
- get_comp(fourier_full_comp[2], m_idx2[1][i]) * get_comp(fourier_full_comp[2], m_idx2[2][i])
- );
- }
-
- // phi^4
- for(int i=0; i<m_idx3[0].size(); ++i)
- {
- cF += 2. *
- (
- get_comp(fourier_full_comp[0], m_idx3[0][i]) * get_comp(fourier_full_comp[0], m_idx3[1][i]) +
- get_comp(fourier_full_comp[1], m_idx3[0][i]) * get_comp(fourier_full_comp[1], m_idx3[1][i]) +
- get_comp(fourier_full_comp[2], m_idx3[0][i]) * get_comp(fourier_full_comp[2], m_idx3[1][i])
- ) *
- (
- get_comp(fourier_full_comp[0], m_idx3[2][i]) * get_comp(fourier_full_comp[0], m_idx3[3][i]) +
- get_comp(fourier_full_comp[1], m_idx3[2][i]) * get_comp(fourier_full_comp[1], m_idx3[3][i]) +
- get_comp(fourier_full_comp[2], m_idx3[2][i]) * get_comp(fourier_full_comp[2], m_idx3[3][i])
- );
- }
-
-
- // zee
- cF += -m_B*fourier0;
- return cF.real();
-}
-
-
-/**
- * set fourier components
- */
-template<class t_real, class t_cplx, int ORDER>
-void Heli<t_real, t_cplx, ORDER>::SetFourier(const std::vector<ublas::vector<t_cplx>> &fourier)
-{
- m_fourier = fourier;
- while(m_fourier.size() < ORDER_FOURIER+1)
- m_fourier.emplace_back(tl2::make_vec<ublas::vector<t_cplx>>({0., 0., 0.}));
-}
-
-
-
-/**
- * energies and spectral weights
- */
-template<class t_real, class t_cplx, int ORDER>
-std::tuple<std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>>
-Heli<t_real, t_cplx, ORDER>::GetSpecWeights(t_real qh, t_real qk, t_real ql, t_real minE, t_real maxE) const
-{
- constexpr int SIZE = 2*ORDER+1;
- constexpr t_cplx imag = t_cplx(0,1);
- static const std::vector<t_real> empty;
- const t_real epshkl = 1e-5;
-
- t_vec qvec = tl2::make_vec<t_vec>({ qh, qk, -ql });
-
- if(tl2::float_equal<t_real>(qvec[2], 0., epshkl))
- qvec[2] += epshkl;
-
-
- const t_real Brel = m_B / m_Bc2;
- const t_real Brel2 = std::sqrt(0.5 - 0.5*Brel*Brel);
- const t_real sqrtfac = std::sqrt(Brel*Brel + 2.*Brel2*Brel2);
- const t_real E_scale_fac = (Brel*Brel + 2.*Brel2*Brel2) / std::sqrt(Brel2*Brel2) * g_g<t_real>*g_muB<t_real>*m_Bc2 * Brel2;
-
-
- // static susceptibility
- t_mat_cplx fluct = tl2::zero_m<t_mat_cplx>(3*SIZE, 3*SIZE);
-
- constexpr t_real A = g_hoc<t_real>;
- constexpr t_real A2 = A*A;
- constexpr t_real A3 = A2*A;
- static const/*expr*/ t_real _c1 = (-2.*imag*std::pow(2., 2./3.) * std::pow(3., 5./6.) * A *
- std::pow(-9.*A2 + std::sqrt(t_cplx(3.*(A3*(2.+27.*A)))), 1./3.) /
- (std::pow(2., 1./3.) * (3.+imag*std::sqrt(3.)) * A + std::pow(3., 1./6.) * (std::sqrt(3.)-imag) *
- std::pow(-9.*A2 + std::sqrt(t_cplx(3.*(A3*(2.+27.*A)))), 2./3.))).real();
- static const/*expr*/ t_real _c2 = (std::pow(std::pow(2., 1./3.) * (std::sqrt(3.)-3.*imag) * A -
- imag*std::pow(3.,1./6.)*(std::sqrt(3.)-imag) *
- std::pow(-9.*A2 + std::sqrt(t_cplx(3.*(A3*(2.+27.*A)))), 2./3.), 2.) /
- (24.*std::pow(2., 1./3.) * A*std::pow(-27.*A2+3.*std::sqrt(t_cplx(3.*(A3*(2.+27.*A)))), 2./3.))).real();
- static const/*expr*/ t_real _c3 = ((324.*A3 - 9.*imag*(std::sqrt(3.)-imag)*A2*std::pow(-54.*A2+6.*std::sqrt(t_cplx(3.*(A3*(2.+27.*A)))), 1./3.) +
- (3.*imag + std::sqrt(3)) * std::sqrt(t_cplx(A3*(2.+27.*A))) * std::pow(-54.*A2+6.*std::sqrt(t_cplx(3.*(A3*(2.+27.*A)))), 1./3.) -
- A*(36.*std::sqrt(t_cplx(3.*(A3*(2.+27.*A)))) -3.*imag*std::pow(3., 1./6.) * std::pow(-18.*A2 + 2.*std::sqrt(t_cplx(3.*(A3*(2.+27.*A)))), 2./3.) +
- std::pow(-54.*A2 + 6.*std::sqrt(t_cplx(3.*(A3*(2.+27.*A)))), 2./3.))) /
- (2.*std::pow(2., 1./3.) * std::pow(3., 1./6.) * std::pow(-9.*A2+std::sqrt(t_cplx(3.*(A3*(2.+27.*A)))), 2./3.) *
- (std::pow(2., 1./3.) * (-3.*imag+std::sqrt(3.))*A - imag*std::pow(3., 1./6) *
- (-imag+std::sqrt(3.)) * std::pow(-9.*A2+std::sqrt(t_cplx(3.*(A3*(2.+27.*A)))), 2./3.)))).real();
- const t_real _c4 = _c3 + 88./3. * Brel2*Brel2;
- const t_real _c5 = _c3 + 88./3. * Brel*Brel;
-
- for(int pos=0; pos<int(SIZE); ++pos)
- {
- t_real qz = qvec[2] - t_real(ORDER) + t_real(pos);
- const t_real q2 = qvec[0]*qvec[0] + qvec[1]*qvec[1] + qz*qz;
-
- fluct(3*pos + 0, 3*pos + 0) = (g_chi<t_real>*(qvec[0]*qvec[0] + qvec[1]*qvec[1]) + 2.*q2*(2.*_c1*qz + q2 + _c2*q2*q2 + _c4)) / (2.*q2);
- fluct(3*pos + 0, 3*pos + 1) = (g_chi<t_real>*(qvec[0] - imag*qvec[1])*(qvec[0] - imag*qvec[1])) / (2.*q2);
- fluct(3*pos + 0, 3*pos + 2) = ((qvec[0] - imag*qvec[1]) * (g_chi<t_real>*qz - 2*_c1*q2) * std::sqrt(2)) / (2.*q2);
-
- fluct(3*pos + 1, 3*pos + 0) = std::conj(fluct(3*pos + 0, 3*pos + 1));
- fluct(3*pos + 1, 3*pos + 1) = (g_chi<t_real>*(qvec[0]*qvec[0] + qvec[1]*qvec[1]) + 2.*q2*(-2.*_c1*qz + q2 + _c2*q2*q2 + _c4)) / (2.*q2);
- fluct(3*pos + 1, 3*pos + 2) = ((qvec[0] + imag*qvec[1]) * (g_chi<t_real>*qz + 2.*_c1*q2) * std::sqrt(2.)) / (2.*q2);
-
- fluct(3*pos + 2, 3*pos + 0) = std::conj(fluct(3*pos + 0, 3*pos + 2));
- fluct(3*pos + 2, 3*pos + 1) = std::conj(fluct(3*pos + 1, 3*pos + 2));
- fluct(3*pos + 2, 3*pos + 2) = (2.*g_chi<t_real>*qz*qz + 2.*q2*(q2 + _c2*q2*q2 + _c5)) / (2.*q2);
-
- if(pos > 1)
- fluct(3*pos + 1, 3*(pos-2) + 0) = Brel2 * Brel2 * 88./3.;
- if(pos > 0)
- fluct(3*pos + 1, 3*(pos-1) + 2) = fluct(3*pos + 2, 3*(pos-1) + 0) = Brel * Brel2 * 88./3.;
- if(pos < int(SIZE)-1)
- fluct(3*pos + 0, 3*(pos+1) + 2) = fluct(3*pos + 2, 3*(pos+1) + 1) = Brel * Brel2 * 88./3.;
- if(pos < int(SIZE)-2)
- fluct(3*pos + 0, 3*(pos+2) + 1) = Brel2 * Brel2 * 88./3.;
- }
-
-
- // Mx
- t_mat_cplx Mx = tl2::zero_m<t_mat_cplx>(3*SIZE, 3*SIZE);
- for(int pos=0; pos<int(SIZE); ++pos)
- {
- if(pos > 0)
- {
- Mx(3*pos + 2, 3*(pos-1) + 0) = imag*Brel2;
- Mx(3*pos + 1, 3*(pos-1) + 2) = -imag*Brel2;
- }
-
- Mx(3*pos + 0, 3*pos + 0) = -imag*Brel;
- Mx(3*pos + 1, 3*pos + 1) = imag*Brel;
-
- if(pos < int(SIZE)-1)
- {
- Mx(3*pos + 0, 3*(pos+1) + 2) = imag*Brel2;
- Mx(3*pos + 2, 3*(pos+1) + 1) = -imag*Brel2;
- }
- }
-
-
- // energies and weights
- return calc_weights<t_mat_cplx, t_vec_cplx, t_cplx, t_real>(
- Mx, fluct,
- m_bProjNeutron, m_projNeutron, m_polMat,
- sqrtfac, E_scale_fac,
- minE, maxE,
- m_eveps, /*m_evlimit*/ -1., m_weighteps,
- m_filterzeroweight, m_onlymode,
- 3*ORDER);
-}
-
-
-/**
- * set lattice vector and orthogonal projector
- */
-template<class t_real, class t_cplx, int ORDER>
-void Heli<t_real, t_cplx, ORDER>::SetG(t_real h, t_real k, t_real l)
-{
- bool bInChiralBase = true;
- m_Grlu = tl2::make_vec<t_vec>({h,k,l});
-
- t_vec _G = m_Grlu / tl2::veclen(m_Grlu);
- _G = tl2::quat_vec_prod(m_rotCoord, _G);
- t_vec_cplx G = _G;
-
- if(bInChiralBase)
- {
- t_mat_cplx chiral = get_chiralbasismat<t_mat_cplx, t_vec_cplx>();
- G = tl2::prod_mv<t_vec_cplx, t_mat_cplx>(chiral, G);
- }
-
- m_projNeutron = tl2::unit_m<t_mat_cplx>(3);
- m_projNeutron -= tl2::outer_cplx<t_vec_cplx, t_mat_cplx>(G, G);
-
- if(bInChiralBase)
- {
- m_projNeutron = tl2::conjugate_mat(m_projNeutron);
- }
-}
-
-
-/**
- * rotates field to internal [001] convention
- */
-template<class t_real, class t_cplx, int ORDER>
-void Heli<t_real, t_cplx, ORDER>::SetCoords(t_real Bx, t_real By, t_real Bz)
-{
- t_vec B = tl2::make_vec<t_vec>( {Bx, By, Bz} );
- t_quat quatB = tl2::rotation_quat(B, tl2::make_vec<t_vec>( {0, 0, 1} ));
-
- m_rotCoord = quatB;
-}
-
-
-/**
- * dispersion
- */
-template<class t_real, class t_cplx, int ORDER>
-std::tuple<std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>>
-Heli<t_real, t_cplx, ORDER>::GetDisp(t_real h, t_real k, t_real l, t_real minE, t_real maxE) const
-{
- t_vec Qvec = tl2::make_vec<t_vec>({ h, k, l });
- t_vec qvec = tl2::quat_vec_prod(m_rotCoord, Qvec) - tl2::quat_vec_prod(m_rotCoord, m_Grlu);
-
- qvec /= g_kh_rlu<t_real>(m_T);
- return GetSpecWeights(qvec[0], qvec[1], qvec[2], minE, maxE);
-}
diff --git a/src/core/heli.h b/src/core/heli.h
deleted file mode 100644
index 971c522..0000000
--- a/src/core/heli.h
+++ /dev/null
@@ -1,111 +0,0 @@
-/**
- * Free energy in the helimagnetic phase
- * @author tweber@ill.fr
- * @date mid-16, jul-18
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#ifndef __HELI_H__
-#define __HELI_H__
-
-#include <tuple>
-#include <vector>
-#include <complex>
-
-#include "helper.h"
-#include "magsys.h"
-
-#define USE_LAPACK
-#include "tlibs2/libs/math17.h"
-
-
-#ifndef DEF_HELI_ORDER
- #define DEF_HELI_ORDER 7
-#endif
-
-
-/**
- * helix
- */
-template<class t_real=double, class t_cplx = std::complex<t_real>, int ORDER=4>
-class Heli : public MagSystem<t_real, t_cplx, ORDER>, public MagDynamics<t_real, t_cplx>
-{
-public:
- static constexpr int ORDER_FOURIER = ORDER;
- using t_mat = ublas::matrix<t_real>;
- using t_vec = ublas::vector<t_real>;
- using t_quat = boost::math::quaternion<t_real>;
- using t_mat_cplx = ublas::matrix<t_cplx>;
- using t_vec_cplx = ublas::vector<t_cplx>;
-
-
-public:
- Heli();
-
- virtual std::shared_ptr<Heli<t_real, t_cplx, ORDER>> copy() const
- {
- return std::make_shared<Heli<t_real, t_cplx, ORDER>>(*this);
- }
-
- virtual std::shared_ptr<MagSystem<t_real, t_cplx, ORDER_FOURIER>> copyCastSys() const override { return copy(); }
- virtual std::shared_ptr<MagDynamics<t_real, t_cplx>> copyCastDyn() const override { return copy(); }
-
- virtual t_real F() override;
-
- virtual void SetFourier(const std::vector<ublas::vector<t_cplx>> &fourier) override;
- virtual const std::vector<ublas::vector<t_cplx>> &GetFourier() const override { return m_fourier; }
-
- // careful: free energy still uses theory units for T and B, dynamics calculation already uses real units!
- virtual void SetB(t_real B) override { m_B = B; }
- virtual void SetT(t_real T) override { m_T = T; m_Bc2 = get_bc2(m_T, false); }
- t_real GetBC2() const { return m_Bc2; }
-
- using MagSystem<t_real, t_cplx, ORDER_FOURIER>::minimise;
-
-
- void SetG(t_real h, t_real k, t_real l);
- void SetCoords(t_real Bx, t_real By, t_real Bz);
- void SetFilterZeroWeight(bool b) { m_filterzeroweight = b; }
- void SetWeightEps(t_real eps) { m_weighteps = eps; }
- void SetOnlyMode(int iMode) { m_onlymode = iMode; }
- void SetProjNeutron(bool b) { m_bProjNeutron = b; }
-
- std::tuple<std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>>
- GetSpecWeights(t_real qh, t_real qk, t_real ql, t_real minE=-1., t_real maxE=-2.) const;
-
- virtual std::tuple<std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>>
- GetDisp(t_real h, t_real k, t_real l, t_real minE=-1., t_real maxE=-2.) const override;
-
-
-private:
- t_real m_pitch = 1;
- t_real m_B = 0, m_T = -100;
-
- std::vector<ublas::vector<t_cplx>> m_fourier;
-
-
- int m_onlymode = -1;
- bool m_filterzeroweight = false;
- t_real m_eveps = 1e-6;
- t_real m_weighteps = 1e-6;
-
- t_real m_Bc2 = 0;
- t_vec m_Grlu;
- t_quat m_rotCoord;
-
- std::vector<t_mat_cplx> m_polMat =
- {{
- get_polmat<t_mat_cplx>(2), // SF1
- get_polmat<t_mat_cplx>(1), // SF2
- get_polmat<t_mat_cplx>(3) // NSF
- }};
- t_mat_cplx m_projNeutron = tl2::unit_m<t_mat_cplx>(3);
- bool m_bProjNeutron = true;
-
-
-private:
- std::vector<int> m_idx2[3], m_idx3[4];
-};
-
-
-#endif
diff --git a/src/core/helper.h b/src/core/helper.h
deleted file mode 100644
index e8559d5..0000000
--- a/src/core/helper.h
+++ /dev/null
@@ -1,353 +0,0 @@
-/**
- * Helper functions
- * @author Tobias Weber <tweber@ill.fr>
- * @date jul-18
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#ifndef __HELPER_H__
-#define __HELPER_H__
-
-#define USE_LAPACK
-#include "tlibs2/libs/math17.h"
-#include "tlibs2/libs/mag.h"
-namespace ublas = boost::numeric::ublas;
-
-
-/**
- * splits a vector into its x,y,z components
- */
-template<class T = std::complex<double>>
-std::tuple<T, T, T> split_vec3d(const ublas::vector<T>& vec)
-{
- return std::make_tuple(vec[0], vec[1], vec[2]);
-}
-
-
-
-/**
- * array indexing including negative values
- */
-template<class t_arr>
-typename t_arr::value_type& get_comp(t_arr &arr, int idx)
-{
- if(idx >= 0) return arr[idx];
-
- return arr[arr.size() + idx];
-}
-
-
-
-/**
- * matrix indexing including negative values
- */
-template<class t_mat>
-typename t_mat::value_type& get_comp(t_mat &mat, int idx1, int idx2)
-{
- if(idx1 < 0) idx1 = int(mat.size1()) + idx1;
- if(idx2 < 0) idx2 = int(mat.size2()) + idx2;
-
- return mat(idx1, idx2);
-}
-
-
-/**
- * matrix indexing including negative values
- */
-template<class t_mat>
-const typename t_mat::value_type& get_comp(const t_mat &mat, int idx1, int idx2)
-{
- if(idx1 < 0) idx1 = int(mat.size1()) + idx1;
- if(idx2 < 0) idx2 = int(mat.size2()) + idx2;
-
- return mat(idx1, idx2);
-}
-
-
-/**
- * 4-dim tensor indexing including negative values
- */
-template<class t_arr>
-typename t_arr::value_type& get_comp(t_arr& arr, int SIZE, int idx1, int idx2, int idx3, int idx4)
-{
- if(idx1 < 0) idx1 = SIZE + idx1;
- if(idx2 < 0) idx2 = SIZE + idx2;
- if(idx3 < 0) idx3 = SIZE + idx3;
- if(idx4 < 0) idx4 = SIZE + idx4;
-
- return arr[idx1*SIZE*SIZE*SIZE +
- idx2*SIZE*SIZE +
- idx3*SIZE +
- idx4];
-}
-
-
-/**
- * 5-dim tensor indexing including negative values
- */
-template<class t_arr>
-typename t_arr::value_type& get_comp(t_arr& arr, int SIZE, int idx1, int idx2, int idx3, int idx4, int idx5)
-{
- if(idx1 < 0) idx1 = SIZE + idx1;
- if(idx2 < 0) idx2 = SIZE + idx2;
- if(idx3 < 0) idx3 = SIZE + idx3;
- if(idx4 < 0) idx4 = SIZE + idx4;
- if(idx5 < 0) idx5 = SIZE + idx5;
-
- return arr[idx1*SIZE*SIZE*SIZE*SIZE +
- idx2*SIZE*SIZE*SIZE +
- idx3*SIZE*SIZE +
- idx4*SIZE +
- idx5];
-}
-
-
-/**
- * 6-dim tensor indexing including negative values
- */
-template<class t_arr>
-typename t_arr::value_type& get_comp(t_arr& arr, int SIZE, int idx1, int idx2, int idx3, int idx4, int idx5, int idx6)
-{
- if(idx1 < 0) idx1 = SIZE + idx1;
- if(idx2 < 0) idx2 = SIZE + idx2;
- if(idx3 < 0) idx3 = SIZE + idx3;
- if(idx4 < 0) idx4 = SIZE + idx4;
- if(idx5 < 0) idx5 = SIZE + idx5;
- if(idx6 < 0) idx6 = SIZE + idx6;
-
- return arr[idx1*SIZE*SIZE*SIZE*SIZE*SIZE +
- idx2*SIZE*SIZE*SIZE*SIZE +
- idx3*SIZE*SIZE*SIZE +
- idx4*SIZE*SIZE +
- idx5*SIZE +
- idx6];
-}
-
-
-template<class t_arr>
-const typename t_arr::value_type& get_virt_comp(const t_arr& arr, int ORGSIZE, int VIRTSIZE, int ORDER, int idx1, int idx2, int idx3, int idx4)
-{
- if(idx1 < 0) idx1 = VIRTSIZE + idx1;
- if(idx2 < 0) idx2 = VIRTSIZE + idx2;
- if(idx3 < 0) idx3 = VIRTSIZE + idx3;
- if(idx4 < 0) idx4 = VIRTSIZE + idx4;
-
- static const typename t_arr::value_type zero{};
- if((idx1>=ORDER && idx1<VIRTSIZE-ORDER-1) || (idx2>=ORDER && idx2<VIRTSIZE-ORDER-1)
- || (idx3>=ORDER && idx3<VIRTSIZE-ORDER-1) || (idx4>=ORDER && idx4<VIRTSIZE-ORDER-1))
- return zero;
-
- bool bLower1 = (idx1<ORDER);
- bool bLower2 = (idx2<ORDER);
- bool bLower3 = (idx3<ORDER);
- bool bLower4 = (idx4<ORDER);
-
- if(bLower1 && idx1 >= ORGSIZE) return zero;
- if(bLower2 && idx2 >= ORGSIZE) return zero;
- if(bLower3 && idx3 >= ORGSIZE) return zero;
- if(bLower4 && idx4 >= ORGSIZE) return zero;
-
- if(!bLower1)
- {
- idx1 -= VIRTSIZE-ORGSIZE;
- if(idx1 < 0) return zero;
- }
- if(!bLower2)
- {
- idx2 -= VIRTSIZE-ORGSIZE;
- if(idx2 < 0) return zero;
- }
- if(!bLower3)
- {
- idx3 -= VIRTSIZE-ORGSIZE;
- if(idx3 < 0) return zero;
- }
- if(!bLower4)
- {
- idx4 -= VIRTSIZE-ORGSIZE;
- if(idx4 < 0) return zero;
- }
-
- return arr[idx1*ORGSIZE*ORGSIZE*ORGSIZE + idx2*ORGSIZE*ORGSIZE + idx3*ORGSIZE + idx4];
-}
-
-
-/**
- * get chiral basis vectors
- */
-template<class t_vec = ublas::vector<std::complex<double>>>
-t_vec get_chiralbasis(int which)
-{
- using t_cplx = typename t_vec::value_type;
- constexpr t_cplx j = t_cplx(0,1);
-
- t_cplx n = std::sqrt<typename t_cplx::value_type>(2.);
-
- if(which == 1) // SF-+
- return tl2::make_vec<t_vec>( {1, +j, 0} ) / n;
- else if(which == 2) // SF+-
- return tl2::make_vec<t_vec>( {1, -j, 0} ) / n;
- else if(which == 3) // NSF
- return tl2::make_vec<t_vec>( {0, 0, 1} );
-
- // none
- return tl2::make_vec<t_vec>( {0, 0, 0} );
-}
-
-
-/**
- * get chiral basis matrix
- */
-template<class t_mat = ublas::matrix<std::complex<double>>,
- class t_vec = ublas::vector<std::complex<double>>>
-t_mat get_chiralbasismat()
-{
- t_vec vec1 = get_chiralbasis<t_vec>(1);
- t_vec vec2 = get_chiralbasis<t_vec>(2);
- t_vec vec3 = get_chiralbasis<t_vec>(3);
-
- return tl2::row_matrix<t_mat, t_vec>({vec1, vec2, vec3});
-}
-
-
-/**
- * get polarisation matrix
- */
-template<class t_mat = ublas::matrix<std::complex<double>>>
-t_mat get_polmat(int which)
-{
- using t_cplx = typename t_mat::value_type;
- constexpr t_cplx j = t_cplx(0,1);
-
- // all channels
- if(which < 1 || which > 3) return tl2::unit_m<t_mat>(3);
-
- t_mat mat = tl2::zero_m<t_mat>(3, 3);
- mat(which-1, which-1) = 1; // SF-+, SF+-, NSF
-
- return mat;
-}
-
-
-/**
- * get polarisation matrix in chiral basis
- */
-template<class t_mat = ublas::matrix<std::complex<double>>>
-t_mat get_chiralpol(int which)
-{
- using t_cplx = typename t_mat::value_type;
- constexpr t_cplx j = t_cplx(0,1);
-
- if(which == 1) // SF-+
- return tl2::make_mat<t_mat>({{0.5, +0.5*j, 0}, {-0.5*j, 0.5, 0}, {0, 0, 0}});
- else if(which == 2) // SF+-
- return tl2::make_mat<t_mat>({{0.5, -0.5*j, 0}, {+0.5*j, 0.5, 0}, {0, 0, 0}});
- else if(which == 3) // NSF
- return tl2::make_mat<t_mat>({{0, 0, 0}, {0, 0, 0}, {0, 0, 1}});
-
- // all channels
- return tl2::unit_m<t_mat>(3);
-}
-
-
-
-/**
- * calculate energies and weights from Landau-Lifshitz Mcross and fluctuation matrices
- * Mcross rotates magnetisation
- * - has 3 eigenvalues: 0, 1, -1:
- * - 0: parallel magnetisation, long. fluctuation, not used
- * - 1, -1: trans. fluctuations
- */
-template<class t_mat_cplx, class t_vec_cplx, class t_cplx, class t_real>
-std::tuple<std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>>
-calc_weights(const t_mat_cplx& Mx, const t_mat_cplx& Fluc,
- bool bProjNeutron, const t_mat_cplx& projNeutron, const std::vector<t_mat_cplx>& polMat,
- t_real normfac=1, t_real E_scale_fac=1, t_real minE=-1, t_real maxE=-2,
- t_real eveps = 1e-6, t_real evlimit = 0.9995, t_real weighteps = 1e-6,
- bool bfilterzeroweight=0, int onlymode=-1, std::size_t MxsubMatRowBegin=0)
-{
- try
- {
- // energies and spectral weights
- struct EW
- {
- t_real E;
- t_real wUnpol;
- t_real wSF1, wSF2, wNSF;
- };
-
- auto eigs = tl2::calc_dynstrucfact_landau<t_mat_cplx, t_vec_cplx, t_cplx, t_real>(
- Mx, Fluc, normfac, &eveps, &evlimit, 3, MxsubMatRowBegin, eveps);
-
- std::size_t numEVals = std::get<0>(eigs).size();
- std::vector<EW> EWs;
- EWs.reserve(numEVals);
-
- int iCurMode = 0;
- for(std::size_t ieval=0; ieval<numEVals; ++ieval)
- {
- auto E_weight = tl2::get_dynstrucfact_neutron<t_mat_cplx, t_vec_cplx, t_cplx, t_real>(
- std::get<0>(eigs)[ieval], std::get<1>(eigs)[ieval], std::get<2>(eigs)[ieval],
- bProjNeutron ? &projNeutron : nullptr, &polMat);
-
- // filter energies if requested
- t_real E_meV = std::get<0>(E_weight) * E_scale_fac;
- if(maxE >= minE && (E_meV < minE || E_meV > maxE))
- continue;
-
- // only consider one of the modes, if desired
- if(onlymode >=0 && onlymode != iCurMode++)
- continue;
-
-
- EW ew;
- ew.E = E_meV;
- ew.wUnpol = std::get<1>(E_weight)[0];
- ew.wSF1 = std::get<1>(E_weight)[1];
- ew.wSF2 = std::get<1>(E_weight)[2];
- ew.wNSF = std::get<1>(E_weight)[3];
- EWs.emplace_back(ew);
- }
-
-
- std::stable_sort(EWs.begin(), EWs.end(), [](const EW& ew1, const EW& ew2) -> bool
- { return ew1.E < ew2.E; });
-
- std::vector<t_real> Es, wsUnpol, wsSF1, wsSF2, wsNSF;
- for(const EW& ew : EWs)
- {
- if(bfilterzeroweight && ew.wUnpol < weighteps)
- continue;
-
- Es.push_back(ew.E);
- wsUnpol.push_back(ew.wUnpol);
- wsSF1.push_back(ew.wSF1);
- wsSF2.push_back(ew.wSF2);
- wsNSF.push_back(ew.wNSF);
- }
-
- return std::make_tuple(Es, wsUnpol, wsSF1, wsSF2, wsNSF);
- }
- catch(const std::exception& ex)
- {
- tl2::log_err(ex.what());
-
- static const std::vector<t_real> empty;
- return std::make_tuple(empty, empty, empty, empty, empty);
- }
-}
-
-
-
-/**
- * appends multiple containers to one another
- */
-template<class t_cont, std::size_t... seq>
-void insert_vals(t_cont& val0, const t_cont& val1, std::index_sequence<seq...>)
-{
- ( std::get<seq>(val0).insert(std::get<seq>(val0).end(),
- std::get<seq>(val1).begin(), std::get<seq>(val1).end()), ... );
-}
-
-
-#endif
diff --git a/src/core/lapacke_switch.h b/src/core/lapacke_switch.h
deleted file mode 100644
index eef628b..0000000
--- a/src/core/lapacke_switch.h
+++ /dev/null
@@ -1,19 +0,0 @@
-/**
- * Lapacke chooser
- * @author tweber@ill.fr
- * @date sep-18
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#ifndef __LAPACKE_SWITCH_H__
-#define __LAPACKE_SWITCH_H__
-
-#if __has_include(<mkl_lapacke.h>)
- #include <mkl_lapacke.h>
- #pragma message "Using MKL Lapacke."
-#else
- #include <lapacke.h>
- #pragma message "Using Vanilla Lapacke."
-#endif
-
-#endif
diff --git a/src/core/longfluct.cpp b/src/core/longfluct.cpp
deleted file mode 100644
index 2521f7f..0000000
--- a/src/core/longfluct.cpp
+++ /dev/null
@@ -1,61 +0,0 @@
-/**
- * Longitudinal fluctuations
- * @author Tobias Weber <tweber@ill.fr>
- * @date 13-jul-20
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#include "longfluct.h"
-#include "constants.h"
-
-
-Longfluct::Longfluct()
-{
- SetPinning(1,-1,0, 0,0,1);
-}
-
-
-/**
- * generate skx satellites depending on the pinning and the up (i.e. B) vector
- */
-void Longfluct::SetPinning(t_real Px, t_real Py, t_real Pz,
- t_real upx, t_real upy, t_real upz)
-{
- m_up = tl2::make_vec<t_vec>({upx, upy, upz});
- m_up /= tl2::veclen(m_up);
-
- t_vec sat = tl2::make_vec<t_vec>({Px, Py, Pz});
- sat = sat / tl2::veclen(sat) * g_kh_rlu_29K<t_real>;
-
- t_mat rot = tl2::rotation_matrix(m_up, tl2::d2r(60.));
-
- m_sats.clear();
- m_sats.push_back(sat);
- for(int i=1; i<6; ++i)
- m_sats.emplace_back(tl2::prod_mv(rot, m_sats[i-1]));
-}
-
-
-/**
- * calculates the dynamical structure factor for the longitudinal fluctuations
- * using the formula from M. Garst, 10/jul/2020.
- */
-Longfluct::t_real Longfluct::S_para(const Longfluct::t_vec& q, Longfluct::t_real E) const
-{
- if(tl2::float_equal<t_real>(m_A, 0))
- return 0.;
-
- constexpr t_real Ec2 = 0.04;
- constexpr t_real kh = g_kh_rlu_29K<t_real>;
-
- t_real S = 0.;
- for(const t_vec& sat : m_sats)
- {
- S += g_kB<t_real> * m_T * m_A /
- (
- (m_inv_correl*m_inv_correl + tl2::inner(q-sat, q-sat)) / (kh*kh) +
- std::pow(m_Gamma*E / Ec2, 2.)
- );
- }
- return S;
-}
diff --git a/src/core/longfluct.h b/src/core/longfluct.h
deleted file mode 100644
index 50ab06c..0000000
--- a/src/core/longfluct.h
+++ /dev/null
@@ -1,49 +0,0 @@
-/**
- * Longitudinal fluctuations
- * @author Tobias Weber <tweber@ill.fr>
- * @date 13-jul-20
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#define USE_LAPACK
-#include "tlibs2/libs/math17.h"
-
-#include <vector>
-
-
-class Longfluct
-{
-public:
- using t_real = double;
- using t_vec = ublas::vector<t_real>;
- using t_mat = ublas::matrix<t_real>;
-
-public:
- Longfluct();
-
- void SetPinning(t_real Px=1, t_real Py=-1, t_real Pz=0,
- t_real upx=0, t_real upy=0, t_real upz=1);
-
- t_real S_para(const t_vec& q, t_real E) const;
-
- void SetT(t_real T) { m_T = T; }
- void SetGamma(t_real G) { m_Gamma = G; }
- void SetInvCorrel(t_real k) { m_inv_correl = k; }
- void SetA(t_real A) { m_A = A; }
-
- t_real GetT() const { return m_T; }
- t_real GetGamma() const { return m_Gamma; }
- t_real GetInvCorrel() const { return m_inv_correl; }
- t_real GetA() const { return m_A; }
-
-private:
- // skx satellites
- std::vector<t_vec> m_sats;
-
- t_vec m_up;
-
- t_real m_T = 29.;
- t_real m_Gamma = 0.5;
- t_real m_inv_correl = 0.1;
- t_real m_A = 0.;
-};
diff --git a/src/core/magsys.cpp b/src/core/magsys.cpp
deleted file mode 100644
index 2141809..0000000
--- a/src/core/magsys.cpp
+++ /dev/null
@@ -1,315 +0,0 @@
-/**
- * Magnetic system interface
- * @author Tobias Weber <tweber@ill.fr>
- * @date jul-18
- * @desc This file implements the theoretical skyrmion model by M. Garst and J. Waizner, see:
- * - https://doi.org/10.1088/1361-6463/aa7573
- * - https://kups.ub.uni-koeln.de/7937/
- * @desc This file is based on:
- * - the descriptions and Mathematica implementations of the different skyrmion model versions by M. Garst and J. Waizner, 2016-2020,
- * - the 2016 Python implementations by M. Kugler and G. Brandl of the first version of the skyrmion model.
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#include "magsys.h"
-
-#define USE_LAPACK
-#include "tlibs2/libs/math17.h"
-#include "tlibs2/libs/str.h"
-#include "tlibs2/libs/log.h"
-
-#ifndef NO_MINIMISATION
- #include "tlibs2/libs/fit.h"
-#endif
-
-#include <fstream>
-
-
-// instantiation
-#ifdef DEF_SKX_ORDER
- template class MagSystem<double, std::complex<double>, (DEF_SKX_ORDER+1)*DEF_SKX_ORDER/2>;
-#endif
-#ifdef DEF_HELI_ORDER
- template class MagSystem<double, std::complex<double>, DEF_HELI_ORDER>;
-#endif
-
-
-#ifndef NO_MINIMISATION
-
-
-/**
- * minimisation of free energy
- */
-template<class t_real, class t_cplx, int ORDER_FOURIER>
-bool MagSystem<t_real, t_cplx, ORDER_FOURIER>::minimise(
- int iMaxOrder, bool bFixXR, bool bFixYR, bool bFixZR, bool bFixXI, bool bFixYI, bool bFixZI)
-{
- auto& sys = *this;
- const auto& fourier = sys.GetFourier();
-
- std::vector<tl2::t_real_min> vars, errs;
- std::vector<std::string> names;
- std::vector<bool> fixed;
- int fourier_idx = 0;
- t_real err_def = 5.;
-
- for(const auto &vec : fourier)
- {
- vars.push_back(vec[0].real()); vars.push_back(vec[0].imag());
- vars.push_back(vec[1].real()); vars.push_back(vec[1].imag());
- vars.push_back(vec[2].real()); vars.push_back(vec[2].imag());
-
- errs.push_back(err_def); errs.push_back(err_def);
- errs.push_back(err_def); errs.push_back(err_def);
- errs.push_back(err_def); errs.push_back(err_def);
-
- names.push_back(std::string("m") + tl2::var_to_str(fourier_idx) + std::string("_x_real"));
- names.push_back(std::string("m") + tl2::var_to_str(fourier_idx) + std::string("_x_imag"));
- names.push_back(std::string("m") + tl2::var_to_str(fourier_idx) + std::string("_y_real"));
- names.push_back(std::string("m") + tl2::var_to_str(fourier_idx) + std::string("_y_imag"));
- names.push_back(std::string("m") + tl2::var_to_str(fourier_idx) + std::string("_z_real"));
- names.push_back(std::string("m") + tl2::var_to_str(fourier_idx) + std::string("_z_imag"));
-
- /*if(fourier_idx > iMaxOrder)
- {
- fixed.push_back(1); fixed.push_back(1);
- fixed.push_back(1); fixed.push_back(1);
- fixed.push_back(1); fixed.push_back(1);
- }
- else*/ if(fourier_idx == 0)
- {
- fixed.push_back(1); fixed.push_back(1);
- fixed.push_back(1); fixed.push_back(1);
- fixed.push_back(bFixZR); fixed.push_back(1);
- }
- else
- {
- fixed.push_back(bFixXR); fixed.push_back(bFixXI);
- fixed.push_back(bFixYR); fixed.push_back(bFixYI);
- fixed.push_back(bFixZR); fixed.push_back(bFixZI);
- }
-
- ++fourier_idx;
- }
-
-
- // convert minimiser args to complex vector
- auto vars_to_fourier = [](const std::vector<tl2::t_real_min>& args)
- -> std::vector<ublas::vector<t_cplx>>
- {
- constexpr auto imag = t_cplx(0,1);
-
- std::vector<ublas::vector<t_cplx>> _fourier;
- for(std::size_t i=0; i<args.size(); i+=6)
- {
- auto m = tl2::make_vec<ublas::vector<t_cplx>>(
- { args[i]+imag*args[i+1], args[i+2]+imag*args[i+3], args[i+4]+imag*args[i+5] });
- _fourier.emplace_back(std::move(m));
- }
- return _fourier;
- };
-
-
- // lambda function to minimise
- auto func = [this, &vars_to_fourier](const auto& ...pack) -> tl2::t_real_min
- {
- const std::vector<tl2::t_real_min>& args = {{ pack... }};
-
- auto newsys = this->copyCastSys();
- newsys->SetFourier(vars_to_fourier(args));
-
- return newsys->F();
- };
-
- constexpr std::size_t NUM_MINPARAMS = 6*(ORDER_FOURIER+1);
- if(vars.size() != NUM_MINPARAMS)
- {
- tl2::log_err("Fourier component number mismatch: ", vars.size(), " != ", NUM_MINPARAMS, ".");
- return false;
- }
-
- bool ok = tl2::minimise<NUM_MINPARAMS>(func, names, vars, errs, &fixed, m_debug);
- sys.SetFourier(vars_to_fourier(vars));
- return ok;
-}
-
-
-
-/**
- * calculates and saves states for phase diagram
- */
-template<class t_real, class t_cplx, int ORDER_FOURIER>
-bool MagSystem<t_real, t_cplx, ORDER_FOURIER>::SaveStates(
- const char *file, int iMaxOrder,
- bool bFixXR, bool bFixYR, bool bFixZR, bool bFixXI, bool bFixYI, bool bFixZI) const
-{
- bool calc_angles = 1;
-
- std::ofstream ofstr(file);
- std::shared_ptr<std::ofstream> ofstrGpl, ofstrGplBc2;
-
- if(calc_angles)
- {
- ofstrGpl.reset(new std::ofstream(std::string(file)+"_ang.gpl"));
- ofstrGplBc2.reset(new std::ofstream(std::string(file)+"_Bc2.gpl"));
-
- ofstrGpl->precision(8);
- ofstrGplBc2->precision(8);
-
- (*ofstrGpl) << "#!/usr/local/bin/gnuplot -p\n\n";
- (*ofstrGplBc2) << "#!/usr/local/bin/gnuplot -p\n\n";
- }
-
-
- if(!ofstr)
- return false;
- ofstr.precision(8);
-
- ofstr << "<states>\n";
- ofstr << "\t<order> " << ORDER_FOURIER << " </order>\n";
- ofstr << "\t<chi> " << g_chi<t_real> << " </chi>\n";
- ofstr << "\t<hoc> " << g_hoc<t_real> << " </hoc>\n";
-
-
- auto mag = this->copyCastSys();
-
- std::vector<t_real> Ts, Bs;
- for(t_real T = -20000; T < 0; T += 250.)
- Ts.push_back(T);
- for(t_real B = 0; B < 200; B += 2.5)
- Bs.push_back(B);
-
-
- if(ofstrGpl)
- {
- (*ofstrGpl) << "set xlabel \"T\"\nset ylabel \"B\"\n";
- (*ofstrGpl) << "set autosc xfix\nset autosc yfix\nset cbrange [0 : pi/2]\n\n";
- (*ofstrGpl) << "plot \"-\" mat nonuni w ima\n";
- (*ofstrGpl) << std::setw(16) << "0" << " ";
- for(t_real T : Ts)
- (*ofstrGpl) << std::setw(16) << T << " ";
- (*ofstrGpl) << "\n";
- }
-
-
- std::vector<std::vector<t_real>> Bc2s;
- Bc2s.resize(Ts.size());
-
- std::size_t iState = 0;
- for(t_real B : Bs)
- {
- if(ofstrGpl) (*ofstrGpl) << std::setw(16) << B << " ";
-
- for(std::size_t T_idx=0; T_idx<Ts.size(); ++T_idx)
- {
- t_real T = Ts[T_idx];
-
- mag->SetT(T);
- mag->SetB(B);
-
- bool ok = mag->minimise(iMaxOrder, bFixXR, bFixYR, bFixZR, bFixXI, bFixYI, bFixZI);
- if(T_idx == 0) // run twice to get better initial values
- ok = mag->minimise(iMaxOrder, bFixXR, bFixYR, bFixZR, bFixXI, bFixYI, bFixZI);
- const auto& fourier = mag->GetFourier();
- auto f = mag->F();
-
- const std::string labState = "state_" + tl2::var_to_str(iState);
- ofstr << "\t<" << labState << ">\n";
-
- t_real ang = std::atan2(std::norm(fourier[1][0])*std::sqrt(t_real(2)), std::norm(fourier[0][2]));
- if(ok && ang < tl2::pi<t_real>/4. && ang > 0.01)
- {
- t_real Bc2 = B / std::cos(ang);
- Bc2s[T_idx].push_back(Bc2);
-
- ofstr << "\t\t<ang> " << ang << " </ang>\n";
- ofstr << "\t\t<Bc2> " << Bc2 << " </Bc2>\n";
- }
-
- ofstr << "\t\t<ok> " << int(ok) << " </ok>\n";
- ofstr << "\t\t<B> " << B << " </B>\n";
- ofstr << "\t\t<T> " << T << " </T>\n";
- ofstr << "\t\t<F> " << f << " </F>\n";
- for(std::size_t i=0; i<fourier.size(); ++i)
- {
- const std::string labM = "M_" + tl2::var_to_str(i);
- ofstr << "\t\t<" << labM << "> "
- << fourier[i][0] << ";\t" << fourier[i][1] << ";\t" << fourier[i][2]
- << " </" << labM << ">\n";
- }
-
-
- if(ofstrGpl) (*ofstrGpl) << std::setw(16) << ang << " ";
-
- std::cout << "\rB = " << B << ", T = " << T << " ";
- std::cout.flush();
-
- ofstr << "\t</" << labState << ">\n";
- ++iState;
- }
-
- if(ofstrGpl) (*ofstrGpl) << "\n";
- }
-
- if(ofstrGpl) (*ofstrGpl) << "e\n";
- std::cout << std::endl;
- ofstr << "</states>\n";
-
-
- if(ofstrGplBc2)
- {
- std::ostringstream ostrData;
- ostrData.precision(ofstrGplBc2->precision());
- for(std::size_t T_idx=0; T_idx<Ts.size(); ++T_idx)
- {
- t_real T = Ts[T_idx];
-
- t_real mean = tl2::mean_value(Bc2s[T_idx]);
- t_real dev = tl2::std_dev(Bc2s[T_idx]);
-
- ostrData << std::setw(16) << T << " ";
- ostrData << std::setw(16) << mean << " ";
- ostrData << std::setw(16) << dev << "\n";
- }
- ostrData << "e\n";
-
- // fit
- (*ofstrGplBc2) << "amp = 1\nex = 0.5\n";
- (*ofstrGplBc2) << "func(T) = amp * (-T)**ex\n\n";
-
- (*ofstrGplBc2) << "fit func(x) \"-\" u 1:2:3 yerr via amp, ex\n";
- (*ofstrGplBc2) << ostrData.str() << "\n";
-
- // plot
- (*ofstrGplBc2) << "set xlabel \"T\"\nset ylabel \"Bc2\"\n";
- (*ofstrGplBc2) << "set autosc xfix\nset autosc yfix\n\n";
-
- (*ofstrGplBc2) << "plot \"-\" u 1:2:3 w yerrorbars pt 7, func(x) w lines lw 2\n";
- (*ofstrGplBc2) << ostrData.str();
- }
-
- return true;
-}
-
-
-#else
-
-
-template<class t_real, class t_cplx, int ORDER_FOURIER>
-bool MagSystem<t_real, t_cplx, ORDER_FOURIER>::minimise(
- int iMaxOrder, bool bFixXR, bool bFixYR, bool bFixZR, bool bFixXI, bool bFixYI, bool bFixZI)
-{
- return 0;
-}
-
-
-template<class t_real, class t_cplx, int ORDER_FOURIER>
-bool MagSystem<t_real, t_cplx, ORDER_FOURIER>::SaveStates(
- const char *file, int iMaxOrder,
- bool bFixXR, bool bFixYR, bool bFixZR, bool bFixXI, bool bFixYI, bool bFixZI) const
-{
- return 0;
-}
-
-
-#endif
diff --git a/src/core/magsys.h b/src/core/magsys.h
deleted file mode 100644
index 4973777..0000000
--- a/src/core/magsys.h
+++ /dev/null
@@ -1,74 +0,0 @@
-/**
- * Magnetic system interface
- * @author Tobias Weber <tweber@ill.fr>
- * @date jul-18
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#ifndef __MAGSYS_H__
-#define __MAGSYS_H__
-
-
-#include <vector>
-#include <complex>
-#include <memory>
-#include "helper.h"
-#include "constants.h"
-
-
-template<class t_real = double>
-class HasTandB
-{
-public:
- virtual void SetB(t_real B) = 0;
- virtual void SetT(t_real T) = 0;
-};
-
-
-/**
- * interface for a magnetic system having a free energy and fourier components defined
- */
-template<class t_real = double, class t_cplx = std::complex<t_real>, int ORDER_FOURIER=2>
-class MagSystem : public HasTandB<t_real>
-{
-public:
- using value_type = t_real;
- using value_type_c = t_cplx;
-
-public:
- virtual t_real F() = 0;
-
- virtual void SetFourier(const std::vector<ublas::vector<t_cplx>> &fourier) = 0;
- virtual const std::vector<ublas::vector<t_cplx>> &GetFourier() const = 0;
-
- virtual bool minimise(int iMaxOrder = 9999,
- bool bFixXR=0, bool bFixYR=0, bool bFixZR=0, bool bFixXI=0, bool bFixYI=0, bool bFixZI=0);
-
- virtual bool SaveStates(const char *file, int iMaxOrder = 9999,
- bool bFixXR=0, bool bFixYR=0, bool bFixZR=0, bool bFixXI=0, bool bFixYI=0, bool bFixZI=0) const;
-
- void SetDebug(bool b) { m_debug = b; }
-
- virtual std::shared_ptr<MagSystem<t_real, t_cplx, ORDER_FOURIER>> copyCastSys() const = 0;
-
-protected:
- bool m_debug = false;
-};
-
-
-/**
- * interface for a magnon dynamics
- */
-template<class t_real = double, class t_cplx = std::complex<t_real>>
-class MagDynamics : public HasTandB<t_real>
-{
-public:
- // unpol, SF-+, SF+-, NSF
- virtual std::tuple<std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>>
- GetDisp(t_real h, t_real k, t_real l, t_real minE=-1., t_real maxE=-2.) const = 0;
-
- virtual std::shared_ptr<MagDynamics<t_real, t_cplx>> copyCastDyn() const = 0;
-};
-
-
-#endif
diff --git a/src/core/skx.cpp b/src/core/skx.cpp
deleted file mode 100644
index 4169f54..0000000
--- a/src/core/skx.cpp
+++ /dev/null
@@ -1,615 +0,0 @@
-/**
- * Free energy in the skyrmion phase
- * @author tweber@ill.fr
- * @date jul-18
- * @desc This file implements the theoretical skyrmion model by M. Garst and J. Waizner, see:
- * - https://doi.org/10.1088/1361-6463/aa7573
- * - https://kups.ub.uni-koeln.de/7937/
- * @desc This file is based on:
- * - the descriptions and Mathematica implementations of the different skyrmion model versions by M. Garst and J. Waizner, 2016-2020,
- * - the 2016 Python implementations by M. Kugler and G. Brandl of the first version of the skyrmion model.
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#include <cmath>
-#include <numeric>
-#include <functional>
-#include <iostream>
-
-#include "skx.h"
-
-
-// instantiation
-#ifdef DEF_SKX_ORDER
- #pragma message("Skx Order: " __TL2_STRCONV(DEF_SKX_ORDER))
- template class Skx<double, std::complex<double>, DEF_SKX_ORDER>;
-
- #ifdef __HACK_FULL_INST__
- template Skx<double, std::complex<double>, DEF_SKX_ORDER>::Skx();
- template void Skx<double, std::complex<double>, DEF_SKX_ORDER>::SetG(double, double, double);
- template void Skx<double, std::complex<double>, DEF_SKX_ORDER>::SetCoords(double, double, double, double,double,double);
- #endif
-#endif
-
-
-/**
- * constructor
- */
-template<class t_real, class t_cplx, int ORDER>
-Skx<t_real, t_cplx, ORDER>::Skx()
-{
- // B matrix
- tl2::inverse(m_Bmat, m_Binv);
-
- // rotation in lab
- t_mat curRot = tl2::unit_m<t_mat>(2);
- const auto rot60 = tl2::rotation_matrix_2d<t_mat>(tl2::d2r<t_real>(60));
-
- //m_rot.push_back(curRot);
- for(int i=0; i<6; ++i)
- {
- curRot = tl2::prod_mm(curRot, rot60);
- m_rot.push_back(curRot);
- }
-
- // rotation in rlu
- for(const auto& rot : m_rot)
- {
- auto M = tl2::transform<t_mat>(rot, m_Bmat, 0);
- m_rotRlu.emplace_back(M);
- }
-
-
- // all peaks
- for(t_real h=-ORDER; h<ORDER+1; ++h)
- for(t_real k=-ORDER; k<ORDER+1; ++k)
- if(std::abs(h-k) <= t_real(ORDER))
- m_allpeaks.emplace_back(tl2::make_vec<t_vec>({h,k}));
-
- for(t_real h=0; h<ORDER+1; ++h)
- for(t_real k=0; k<ORDER; ++k)
- if(h>k)
- {
- t_vec vec_rlu = tl2::make_vec<t_vec>({h,k});
- t_vec vec_lab = tl2::prod_mv(m_Bmat, vec_rlu);
- vec_lab /= tl2::veclen(vec_lab);
-
- m_peaks_60.emplace_back(std::move(vec_rlu));
- m_peaks_60_lab.emplace_back(std::move(vec_lab));
- }
-
-
- // the top and bottom 180 degrees
- std::vector<t_vec> peaks_180_t, peaks_180_b;
-
- for(const auto& rot : m_rotRlu)
- {
- std::vector<t_vec> _peaks;
- for(const auto &vec : m_peaks_60)
- {
- t_vec pk_rlu = tl2::prod_mv(rot, vec);
- t_vec pk = tl2::prod_mv(m_Bmat, pk_rlu);
-
- if(tl2::float_equal<t_real>(pk[1], 0.))
- {
- if(pk[0] < 0.)
- peaks_180_b.push_back(pk_rlu);
- else
- peaks_180_t.push_back(pk_rlu);
- }
- else
- {
- if(pk[1] < 0.)
- peaks_180_b.push_back(pk_rlu);
- else
- peaks_180_t.push_back(pk_rlu);
- }
-
- _peaks.emplace_back(pk_rlu);
- }
- m_peaks_360.emplace_back(_peaks);
- }
-
-
- // #1
- for(const auto& vec : peaks_180_t)
- m_idx1.emplace_back(std::make_pair(int(std::round(vec[0])), int(std::round(vec[1]))));
-
- // #2
- for(std::size_t i=0; i<m_allpeaks.size(); ++i)
- {
- std::copy(m_idx1.begin(), m_idx1.end(), std::back_inserter(m_idx2[0]));
- for(std::size_t j=0; j<peaks_180_t.size(); ++j)
- m_idx2[1].emplace_back(std::make_pair(int(std::round(m_allpeaks[i][0])), int(std::round(m_allpeaks[i][1]))));
- }
-
- // #3
- for(std::size_t i=0; i<m_allpeaks.size(); ++i)
- {
- std::copy(m_idx2[0].begin(), m_idx2[0].end(), std::back_inserter(m_idx3[0]));
- std::copy(m_idx2[1].begin(), m_idx2[1].end(), std::back_inserter(m_idx3[1]));
-
- for(std::size_t j=0; j<peaks_180_t.size(); ++j)
- for(std::size_t k=0; k<m_allpeaks.size(); ++k)
- m_idx3[2].emplace_back(std::make_pair(int(std::round(m_allpeaks[i][0])), int(std::round(m_allpeaks[i][1]))));
- }
-
-
- // reduce #2
- decltype(m_idx2) idx2;
- for(std::size_t i=0; i<m_idx2[0].size(); ++i)
- {
- auto val1 = -m_idx2[0][i].first - m_idx2[1][i].first;
- auto val2 = -m_idx2[0][i].second - m_idx2[1][i].second;
- if(std::abs(val1) > ORDER || std::abs(val2) > ORDER || std::abs(val1-val2) > ORDER)
- { continue; }
- else
- {
- idx2[0].push_back(m_idx2[0][i]);
- idx2[1].push_back(m_idx2[1][i]);
- idx2[2].emplace_back(std::make_pair(val1, val2));
- }
- }
- for(int i=0; i<3; ++i)
- m_idx2[i] = std::move(idx2[i]);
-
- // reduce #3
- decltype(m_idx3) idx3;
- for(std::size_t i=0; i<m_idx3[0].size(); ++i)
- {
- auto val1 = -m_idx3[0][i].first - m_idx3[1][i].first - m_idx3[2][i].first;
- auto val2 = -m_idx3[0][i].second - m_idx3[1][i].second - m_idx3[2][i].second;
- if(std::abs(val1) > ORDER || std::abs(val2) > ORDER || std::abs(val1-val2) > ORDER)
- { continue; }
- else
- {
- idx3[0].push_back(m_idx3[0][i]);
- idx3[1].push_back(m_idx3[1][i]);
- idx3[2].push_back(m_idx3[2][i]);
- idx3[3].emplace_back(std::make_pair(val1, val2));
- }
- }
- for(int i=0; i<4; ++i)
- m_idx3[i] = std::move(idx3[i]);
-}
-
-
-template<class t_real, class t_cplx, int ORDER>
-ublas::matrix<typename Skx<t_real, t_cplx, ORDER>::t_vec_cplx>
-Skx<t_real, t_cplx, ORDER>::GetFullFourier() const
-{
- constexpr auto imag = t_cplx(0,1);
-
- ublas::matrix<t_vec_cplx> M(2*ORDER+1, 2*ORDER+1);
- for(std::size_t i=0; i<M.size1(); ++i)
- for(std::size_t j=0; j<M.size2(); ++j)
- M(i,j) = tl2::make_vec<t_vec_cplx>({0,0,0});
- M(0,0) = m_fourier[0];
-
- // generate all skx fourier components
- for(std::size_t ipk=0; ipk<m_peaks_360.size(); ++ipk) // 6
- {
- const auto& rot = m_rot[ipk];
- for(std::size_t ihx=0; ihx<m_peaks_360[ipk].size(); ++ihx) // 10
- {
- int idx1 = int(std::round(m_peaks_360[ipk][ihx][0]));
- int idx2 = int(std::round(m_peaks_360[ipk][ihx][1]));
-
- const auto& vecPk = m_peaks_60_lab[ihx];
- auto fourier = tl2::make_vec<t_vec_cplx>
- ({ -vecPk[1]*m_fourier[ihx+1][0],
- +vecPk[0]*m_fourier[ihx+1][0] });
- fourier = tl2::prod_mv(rot, fourier);
-
- get_comp(M, idx1, idx2) = tl2::make_vec<t_vec_cplx>
- ({ imag*fourier[0], imag*fourier[1], m_fourier[ihx+1][2] });
- }
- }
-
- return M;
-}
-
-
-template<class t_real, class t_cplx, int ORDER>
-void Skx<t_real, t_cplx, ORDER>::GenFullFourier()
-{
- m_M = GetFullFourier();
-}
-
-
-/**
- * free energy
- */
-template<class t_real, class t_cplx, int ORDER>
-t_real Skx<t_real, t_cplx, ORDER>::F()
-{
- constexpr auto imag = t_cplx(0,1);
-
- GenFullFourier();
-
- // free energy
- t_cplx cF = 0;
- auto fourier0 = tl2::veclen(m_M(0,0));
-
- // dip
- cF += g_chi<t_real>/3. * fourier0*fourier0;
-
-
- for(const auto& pair : m_idx1)
- {
- t_vec pos_rlu = tl2::make_vec<t_vec>({t_real(pair.first), t_real(pair.second)});
- t_vec pos_lab = tl2::prod_mv(m_Bmat, pos_rlu);
-
- const auto& m = get_comp(m_M, pair.first, pair.second);
- const auto len = tl2::inner(m, m);
-
- // dmi
- cF += 8. * m_pitch * imag*(pos_lab[0]*m[1]*m[2] - pos_lab[1]*m[2]*m[0]);
-
- // hoc
- cF += 2. * g_hoc<t_real> * m_pitch*m_pitch * len
- * std::pow(pos_rlu[0]*pos_rlu[0] + pos_rlu[1]*pos_rlu[1] - pos_rlu[0]*pos_rlu[1], 2.);
-
- // phi^2 & phi^4
- cF += 2. * m_pitch*m_pitch * len
- * (pos_rlu[0]*pos_rlu[0] + pos_rlu[1]*pos_rlu[1] - pos_rlu[0]*pos_rlu[1]);
- cF += (m_T + 1. + fourier0*fourier0) * 2.*len;
- }
-
-
- // phi^2 & phi^4
- cF += (m_T + 1.) * fourier0*fourier0 + fourier0*fourier0*fourier0*fourier0;
-
-
- for(std::size_t i=0; i<m_idx2[0].size(); ++i)
- {
- const auto& m1 = get_comp(m_M, m_idx2[0][i].first, m_idx2[0][i].second);
- const auto& m2 = get_comp(m_M, m_idx2[1][i].first, m_idx2[1][i].second);
- const auto& m3 = get_comp(m_M, m_idx2[2][i].first, m_idx2[2][i].second);
-
- // the x and y components are purely imaginary, z purely real
- cF += 2. * fourier0 * m1[2] * tl2::inner(m2, m3);
- }
-
-
- for(std::size_t i=0; i<m_idx3[0].size(); ++i)
- {
- const auto& m1 = get_comp(m_M, m_idx3[0][i].first, m_idx3[0][i].second);
- const auto& m2 = get_comp(m_M, m_idx3[1][i].first, m_idx3[1][i].second);
- const auto& m3 = get_comp(m_M, m_idx3[2][i].first, m_idx3[2][i].second);
- const auto& m4 = get_comp(m_M, m_idx3[3][i].first, m_idx3[3][i].second);
-
- // the x and y components are purely imaginary, z purely real
- cF += 2. * tl2::inner(m1, m2) * tl2::inner(m3, m4);
- }
-
-
- // zee
- cF += -m_B*fourier0;
-
- return cF.real();
-}
-
-
-
-/**
- * set fourier components
- */
-template<class t_real, class t_cplx, int ORDER>
-void Skx<t_real, t_cplx, ORDER>::SetFourier(const std::vector<t_vec_cplx> &fourier)
-{
- m_fourier = fourier;
- while(m_fourier.size() < ORDER_FOURIER+1)
- m_fourier.emplace_back(tl2::make_vec<t_vec_cplx>({0., 0., 0.}));
- while(m_fourier.size() > ORDER_FOURIER+1)
- m_fourier.pop_back();
-}
-
-
-
-/**
- * cross-product and fluctuation matrices
- */
-template<class t_real, class t_cplx, int ORDER>
-std::tuple<typename Skx<t_real, t_cplx, ORDER>::t_mat_cplx, typename Skx<t_real, t_cplx, ORDER>::t_mat_cplx>
-Skx<t_real, t_cplx, ORDER>::GetMCrossMFluct(
- int iGh, int iGk, t_real qh, t_real qk, t_real ql) const
-{
- constexpr int SIZE = 2*ORDER+1;
- constexpr auto imag = t_cplx(0,1);
-
- const int CRYSSIZE = ORDER + std::max(std::abs(iGh), std::abs(iGk));
- const int VIRTSIZE = 2*CRYSSIZE + 1;
-
- t_vec q_lab = tl2::prod_mv(m_Bmat, tl2::make_vec<t_vec>({ qh, qk }));
-
-
-
- // ------------------------------------------------------------------------
- // indices
-
- // 2
- std::vector<std::pair<int, int>> idx2[3];
- for(std::size_t i=0; i<m_allpeaks.size(); ++i)
- {
- for(std::size_t j=0; j<m_allpeaks.size(); ++j)
- {
- idx2[0].emplace_back(std::make_pair(int(std::round(m_allpeaks[j][0])), int(std::round(m_allpeaks[j][1]))));
- idx2[1].emplace_back(std::make_pair(int(std::round(m_allpeaks[i][0])), int(std::round(m_allpeaks[i][1]))));
- }
- }
-
- // 3
- std::vector<std::pair<int, int>> idx3[4];
- for(std::size_t i=0; i<m_allpeaks.size(); ++i)
- {
- std::copy(idx2[0].begin(), idx2[0].end(), std::back_inserter(idx3[0]));
- std::copy(idx2[1].begin(), idx2[1].end(), std::back_inserter(idx3[1]));
-
- for(std::size_t j=0; j<m_allpeaks.size(); ++j)
- for(std::size_t k=0; k<m_allpeaks.size(); ++k)
- idx3[2].emplace_back(std::make_pair(int(std::round(m_allpeaks[i][0])), int(std::round(m_allpeaks[i][1]))));
- }
-
-
- // reduce 2
- decltype(idx2) _idx2;
- for(std::size_t i=0; i<idx2[0].size(); ++i)
- {
- auto val1 = idx2[0][i].first - idx2[1][i].first;
- auto val2 = idx2[0][i].second - idx2[1][i].second;
- if(std::abs(val1) > ORDER || std::abs(val2) > ORDER || std::abs(val1-val2) > ORDER)
- {
- continue;
- }
- else
- {
- _idx2[0].push_back(idx2[0][i]);
- _idx2[1].push_back(idx2[1][i]);
- _idx2[2].emplace_back(std::make_pair(val1, val2));
- }
- }
- for(int i=0; i<3; ++i)
- idx2[i] = std::move(_idx2[i]);
-
- // reduce 3
- decltype(idx3) _idx3;
- for(std::size_t i=0; i<idx3[0].size(); ++i)
- {
- auto val1 = idx3[0][i].first - idx3[1][i].first - idx3[2][i].first;
- auto val2 = idx3[0][i].second - idx3[1][i].second - idx3[2][i].second;
- if(std::abs(val1) > ORDER || std::abs(val2) > ORDER || std::abs(val1-val2) > ORDER)
- {
- continue;
- }
- else
- {
- _idx3[0].push_back(idx3[0][i]);
- _idx3[1].push_back(idx3[1][i]);
- _idx3[2].push_back(idx3[2][i]);
- _idx3[3].emplace_back(std::make_pair(val1, val2));
- }
- }
- for(int i=0; i<4; ++i)
- idx3[i] = std::move(_idx3[i]);
- // ------------------------------------------------------------------------
-
-
-
- // ------------------------------------------------------------------------
- // Mx_ij = eps_ikj M_k
- auto Mx = std::make_unique<std::array<t_mat_cplx, SIZE*SIZE*SIZE*SIZE>>();
-
- for(std::size_t i=0; i<idx2[2].size(); ++i)
- {
- const auto& vecM = get_comp(m_M, idx2[2][i].first, idx2[2][i].second);
- auto skew = tl2::skew<t_mat_cplx>(vecM);
-
- get_comp(*Mx, SIZE, idx2[0][i].first, idx2[0][i].second, idx2[1][i].first, idx2[1][i].second) = skew;
- }
- // ------------------------------------------------------------------------
-
-
-
- // ------------------------------------------------------------------------
- // fluct. F_q1q2_ij = 0.5 * d^2F/(dM_-q_1_q dM_q2_j)
- auto Fluc = std::make_unique<std::array<t_mat_cplx, SIZE*SIZE*SIZE*SIZE>>();
-
- for(std::size_t i=0; i<idx3[3].size(); ++i)
- {
- const auto& vecM1 = get_comp(m_M, idx3[2][i].first, idx3[2][i].second);
- const auto& vecM2 = get_comp(m_M, idx3[3][i].first, idx3[3][i].second);
-
- t_mat_cplx mat(3,3);
- for(int _midx1=0; _midx1<3; ++_midx1)
- for(int _midx2=0; _midx2<3; ++_midx2)
- mat(_midx1,_midx2) = vecM1[_midx1]*vecM2[_midx2];
- mat *= 4.;
-
- // delta term
- for(int d=0; d<3; ++d)
- mat(d,d) += 2.*tl2::inner(vecM1, vecM2);
-
- auto& oldmat = get_comp(*Fluc, SIZE, idx3[0][i].first, idx3[0][i].second, idx3[1][i].first, idx3[1][i].second);
- if(!oldmat.size1()) oldmat = 2.*mat; else oldmat += 2.*mat;
- }
-
-
- for(const auto &pk_rlu : m_allpeaks)
- {
- t_vec pk_lab = tl2::prod_mv(m_Bmat, pk_rlu);
-
- t_vec pos = pk_lab + q_lab;
- pos.resize(3, true);
- pos[2] = ql;
-
- t_real dipole = g_chi<t_real> / tl2::inner(pos, pos);
-
- t_mat_cplx mat(3,3);
- mat(0,0) = dipole * pos[0]*pos[0] + (1. + m_T) + tl2::inner(pos, pos);
- mat(0,1) = dipole * pos[0]*pos[1] - 2.*imag*pos[2];
- mat(0,2) = dipole * pos[0]*pos[2] + 2.*imag*pos[1];
- mat(1,0) = dipole * pos[1]*pos[0] + 2.*imag*pos[2];
- mat(1,1) = dipole * pos[1]*pos[1] + (1. + m_T) + tl2::inner(pos, pos);
- mat(1,2) = dipole * pos[1]*pos[2] - 2.*imag*pos[0];
- mat(2,0) = dipole * pos[2]*pos[0] - 2.*imag*pos[1];
- mat(2,1) = dipole * pos[2]*pos[1] + 2.*imag*pos[0];
- mat(2,2) = dipole * pos[2]*pos[2] + (1. + m_T) + tl2::inner(pos, pos);
-
- int idx1 = int(std::round(pk_rlu[0]));
- int idx2 = int(std::round(pk_rlu[1]));
- auto& oldmat = get_comp(*Fluc, SIZE, idx1, idx2, idx1, idx2);
- if(!oldmat.size1()) oldmat = 2.*mat; else oldmat += 2.*mat;
- }
- // ------------------------------------------------------------------------
-
-
- auto mk_2dim = [VIRTSIZE, CRYSSIZE, SIZE, iGh, iGk](const /*auto&*/ decltype(*Mx)& arr) -> t_mat_cplx
- {
- std::vector<int> pks1(VIRTSIZE);
- std::iota(pks1.begin(), pks1.begin()+CRYSSIZE+1, 0);
- std::iota(pks1.begin()+CRYSSIZE+1, pks1.end(), -CRYSSIZE);
-
- std::vector<std::pair<int, int>> pks2;
- pks2.reserve(VIRTSIZE*VIRTSIZE);
- for(int k : pks1)
- for(int h : pks1)
- if(std::abs(h-k) <= CRYSSIZE)
- pks2.emplace_back(std::make_pair(h, k));
-
- auto mat = tl2::zero_matrix<t_mat_cplx>(3*pks2.size(), 3*pks2.size());
- for(std::size_t idx1=0; idx1<pks2.size(); ++idx1)
- {
- const auto& hk1 = pks2[idx1];
- for(std::size_t idx2=0; idx2<pks2.size(); ++idx2)
- {
- const auto& hk2 = pks2[idx2];
-
- const auto& comp = get_virt_comp(arr, SIZE, VIRTSIZE, ORDER,
- hk1.first+iGh, hk1.second+iGk, hk2.first+iGh, hk2.second+iGk);
-
- if(comp.size1() && comp.size2())
- {
- for(std::size_t x1=0; x1<3; ++x1)
- for(std::size_t x2=0; x2<3; ++x2)
- mat(idx1*3+x1, idx2*3+x2) = comp(x1,x2);
- }
- }
- }
- return mat;
- };
-
-
- auto Mx2d = mk_2dim(*Mx);
- auto Fluc2d = mk_2dim(*Fluc);
- return std::make_tuple(Mx2d, Fluc2d);
-}
-
-
-/**
- * energies and spectral weights
- */
-template<class t_real, class t_cplx, int ORDER>
-std::tuple<std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>>
-Skx<t_real, t_cplx, ORDER>::GetSpecWeights(int iGhmag, int iGkmag, t_real qh, t_real qk, t_real ql, t_real minE, t_real maxE) const
-{
- const t_real epshkl = 1e-5;
- if(tl2::float_equal<t_real>(qh, 0., epshkl) && tl2::float_equal<t_real>(qk, 0., epshkl) && tl2::float_equal<t_real>(ql, 0., epshkl))
- { qh += epshkl; qk += epshkl; ql += epshkl; }
-
- qh = qh; qk = qk; ql = -ql;
-
- t_real sqrtfac = std::sqrt(t_real(-0.5 - m_T*0.5));
- const t_real E_scale_fac_heli = 0.0387; // calculated with heli.cpp
- t_real E_scale_fac = E_scale_fac_heli / (sqrtfac/10.);
-
- t_mat_cplx Mx2d, Fluc2d;
- std::tie(Mx2d, Fluc2d) = GetMCrossMFluct(iGhmag, iGkmag, qh, qk, ql);
-
- // energies and weights
- return calc_weights<t_mat_cplx, t_vec_cplx, t_cplx, t_real>(
- Mx2d, Fluc2d,
- m_bProjNeutron, m_projNeutron, m_polMat,
- sqrtfac, E_scale_fac,
- minE, maxE,
- m_eveps, m_evlimit, m_weighteps,
- m_filterzeroweight, /*m_onlymode*/-1);
-}
-
-
-/**
- * rotates field and pinning to internal conventions
- * ([001] and [100], respectively])
- */
-template<class t_real, class t_cplx, int ORDER>
-void Skx<t_real, t_cplx, ORDER>::SetCoords(
- t_real Bx, t_real By, t_real Bz,
- t_real Pinx, t_real Piny, t_real Pinz)
-{
- t_vec B = tl2::make_vec<t_vec>( {Bx, By, Bz} );
- t_vec _Pin = tl2::make_vec<t_vec>( {Pinx, Piny, Pinz} );
-
- t_quat quatB = tl2::rotation_quat(B, tl2::make_vec<t_vec>( {0, 0, 1} ));
-
- t_vec Pin = tl2::quat_vec_prod(quatB, _Pin);
- t_quat quatPin = tl2::rotation_quat(Pin, tl2::make_vec<t_vec>( {1, 0, 0} ));
-
- m_rotCoord = quatPin * quatB;
-}
-
-
-template<class t_real, class t_cplx, int ORDER>
-void Skx<t_real, t_cplx, ORDER>::SetG(t_real h, t_real k, t_real l)
-{
- m_Grlu = tl2::make_vec<t_vec>({h,k,l});
-
- t_vec G = m_Grlu / tl2::veclen(m_Grlu);
- G = tl2::quat_vec_prod(m_rotCoord, G);
- m_projNeutron = tl2::unit_m<t_mat>(3) - tl2::outer<t_vec, t_mat>(G, G);
-}
-
-
-/**
- * query the dispersion
- */
-template<class t_real, class t_cplx, int ORDER>
-std::tuple<std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>>
-Skx<t_real, t_cplx, ORDER>::GetDisp(t_real h, t_real k, t_real l, t_real minE, t_real maxE) const
-{
- t_vec Qrlu = tl2::make_vec<t_vec>( {h, k, l} );
- t_vec qrlu = Qrlu - m_Grlu;
-
- t_vec qkh = qrlu / G_KH_RLU_29K;
-
- qkh = tl2::quat_vec_prod(m_rotCoord, qkh);
- t_real _l = qkh[2];
- qkh.resize(2, true);
-
- t_vec Qmagrlu = tl2::prod_mv(m_Binv, qkh);
- t_vec Gmagrlu = tl2::make_vec<t_vec>({
- std::round(Qmagrlu[0]), std::round(Qmagrlu[1]) });
-
-
- static const std::vector<t_vec> sats =
- {
- tl2::make_vec<t_vec>({0, 0}),
- tl2::make_vec<t_vec>({0, -1}), tl2::make_vec<t_vec>({0, +1}),
- tl2::make_vec<t_vec>({-1, 0}), tl2::make_vec<t_vec>({+1, 0}),
- tl2::make_vec<t_vec>({-1, -1}), tl2::make_vec<t_vec>({+1, +1}),
- };
-
- auto iterClosest = std::min_element(sats.begin(), sats.end(),
- [&Gmagrlu, &Qmagrlu](const t_vec& sat1, const t_vec& sat2) -> bool
- {
- t_vec qmag1 = Qmagrlu - (Gmagrlu+sat1);
- t_vec qmag2 = Qmagrlu - (Gmagrlu+sat2);
-
- return (qmag1[0]*qmag1[0] + qmag1[1]*qmag1[1] - qmag1[0]*qmag1[1]) <
- (qmag2[0]*qmag2[0] + qmag2[1]*qmag2[1] - qmag2[0]*qmag2[1]);
- });
-
- Gmagrlu += *iterClosest;
- t_vec qmagrlu = Qmagrlu - Gmagrlu;
-
- return GetSpecWeights(int(Gmagrlu[0]), int(Gmagrlu[1]), qmagrlu[0], qmagrlu[1], _l, minE, maxE);
-}
diff --git a/src/core/skx.h b/src/core/skx.h
deleted file mode 100644
index a50be35..0000000
--- a/src/core/skx.h
+++ /dev/null
@@ -1,127 +0,0 @@
-/**
- * Free energy in the skyrmion phase
- * @author tweber@ill.fr
- * @date jul-18
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#ifndef __SKX_H__
-#define __SKX_H__
-
-#include "magsys.h"
-#define USE_LAPACK
-#include "tlibs2/libs/math17.h"
-
-
-#ifndef DEF_SKX_ORDER
- #define DEF_SKX_ORDER 7
-#endif
-
-
-/**
- * skyrmion
- */
-template<class t_real=double, class t_cplx = std::complex<t_real>, int ORDER=4>
-class Skx : public MagSystem<t_real, t_cplx, (ORDER+1)*ORDER/2>, public MagDynamics<t_real, t_cplx>
-{
-public:
- static constexpr int ORDER_FOURIER = (ORDER+1)*ORDER/2;
-
- using t_mat = ublas::matrix<t_real>;
- using t_vec = ublas::vector<t_real>;
- using t_quat = boost::math::quaternion<t_real>;
- using t_mat_cplx = ublas::matrix<t_cplx>;
- using t_vec_cplx = ublas::vector<t_cplx>;
-
-
-public:
- Skx();
-
- virtual t_real F() override;
-
- virtual void SetFourier(const std::vector<t_vec_cplx> &fourier) override;
- virtual const std::vector<t_vec_cplx> &GetFourier() const override { return m_fourier; }
-
- virtual void SetB(t_real B) override { m_B = B; }
- virtual void SetT(t_real T) override { m_T = T; m_Bc2 = get_bc2(m_T); }
- t_real GetBC2() const { return m_Bc2; }
-
- using MagSystem<t_real, t_cplx, ORDER_FOURIER>::minimise;
-
-
- std::shared_ptr<Skx<t_real, t_cplx, ORDER>> copy() const
- {
- return std::make_shared<Skx<t_real, t_cplx, ORDER>>(*this);
- }
-
- virtual std::shared_ptr<MagSystem<t_real, t_cplx, ORDER_FOURIER>> copyCastSys() const override { return copy(); }
- virtual std::shared_ptr<MagDynamics<t_real, t_cplx>> copyCastDyn() const override { return copy(); }
-
-
- void GenFullFourier();
- std::tuple<t_mat_cplx, t_mat_cplx> GetMCrossMFluct(
- int iGhmag, int iGkmag, t_real qh, t_real qk, t_real ql) const;
- std::tuple<std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>>
- GetSpecWeights(int iGhmag, int iGkmag, t_real qh, t_real qk, t_real ql, t_real minE=-1., t_real maxE=-2.) const;
-
- void SetCoords(t_real Bx, t_real By, t_real Bz,
- t_real Pinx, t_real Piny, t_real Pinz);
- void SetG(t_real h, t_real k, t_real l);
- void SetProjNeutron(bool b) { m_bProjNeutron = b; }
- void SetFilterZeroWeight(bool b) { m_filterzeroweight = b; }
- void SetWeightEps(t_real eps) { m_weighteps = eps; }
- const t_mat_cplx& GetNeutronProjOp() const { return m_projNeutron; }
-
- virtual std::tuple<std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>, std::vector<t_real>>
- GetDisp(t_real h, t_real k, t_real l, t_real minE=-1., t_real maxE=-2.) const override;
-
-
-protected:
- ublas::matrix<t_vec_cplx> GetFullFourier() const;
-
-
-private:
- bool m_filterzeroweight = false;
- t_real m_eveps = 1e-6;
- t_real m_weighteps = 1e-6;
- t_real m_evlimit = 0.9995;
-
- // B matrix with 120 deg between (100) and (010), Q_lab = B*Q_rlu
- t_mat m_Bmat = tl2::make_mat<t_mat>(
- {{1, std::cos(tl2::d2r<t_real>(120))}, {0, std::sin(tl2::d2r<t_real>(120))}});
- t_mat m_Binv;
-
- std::vector<t_mat> m_rot, m_rotRlu;
-
- t_real m_pitch = 1;
- t_real m_B = 0;
- t_real m_T = -100;
- t_real m_Bc2 = 0;
-
- std::vector<t_vec_cplx> m_fourier;
-
- // all six sixth
- std::vector<std::vector<t_vec>> m_peaks_360;
- // one sixth of the magnetic lattice
- std::vector<t_vec> m_peaks_60, m_peaks_60_lab;
- std::vector<t_vec> m_allpeaks;
-
- ublas::matrix<t_vec_cplx> m_M;
- t_vec m_Grlu;
- t_quat m_rotCoord;
-
- std::vector<t_mat_cplx> m_polMat =
- {{
- get_chiralpol<t_mat_cplx>(1), // SF1
- get_chiralpol<t_mat_cplx>(2), // SF2
- get_chiralpol<t_mat_cplx>(3) // NSF
- }};
-
- t_mat_cplx m_projNeutron = tl2::unit_m<t_mat_cplx>(3);
- bool m_bProjNeutron = true;
-
- std::vector<std::pair<int, int>> m_idx1, m_idx2[3], m_idx3[4];
-};
-
-
-#endif
diff --git a/src/core/skx_default_gs.cxx b/src/core/skx_default_gs.cxx
deleted file mode 100644
index d3a7054..0000000
--- a/src/core/skx_default_gs.cxx
+++ /dev/null
@@ -1,62 +0,0 @@
-/**
- * Default ground state for skx
- * @author tweber@ill.fr
- * @date aug-18
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-/**
- * output from skx_gs.cpp
- * x and y values are the imaginary components
- * z values are the real components
- */
-static const std::vector<t_real> _skxgs_allcomps =
-{{
- // order 9
- 0, 0, 10.643255,
- -5.8312418, 0, -5.4285579,
- 0.11772572, 0, 0.092512566,
- -0.60875193, 0, -0.6283387,
- 0.099009681, 0, 0.097708898,
- 0.21345703, 0, 0.20197682,
- 0.21345707, 0, 0.20197861,
- -0.014776251, 0, -0.014232656,
- -0.0046309613, 0, -0.0035734788,
- 0.0081490635, 0, 0.0096266455,
- -0.0046284703, 0, -0.003570242,
- -0.00092182745, 0, -0.00095170469,
- -0.0060308055, 0, -0.0059615124,
- -0.0098962887, 0, -0.0095705205,
- -0.009896991, 0, -0.0095664551,
- -0.006029697, 0, -0.0059582546,
- 0.00066330707, 0, 0.00055053732,
- 0.0010092654, 0, 0.0008918821,
- 0.00071552573, 0, 0.00060852518,
- 0.00046383518, 0, 0.00034571667,
- 0.00071574625, 0, 0.0006113456,
- 0.0010099268, 0, 0.00089342825,
- 5.3331775e-05, 0, -0.00010249534,
- 0.00015539651, 0, 2.5296191e-06,
- 0.00030612867, 0, 0.00024823487,
- 0.00044277835, 0, 0.00042931965,
- 0.00044175346, 0, 0.00042949054,
- 0.00030547693, 0, 0.00024768997,
- 0.00015489058, 0, 1.3435546e-06,
- 6.9271166e-05, 0, 5.7206645e-05,
- -2.7366626e-05, 0, -2.1078224e-05,
- -3.9356e-05, 0, -5.7759202e-05,
- -3.8566738e-05, 0, -4.1978018e-05,
- -2.7839047e-05, 0, -2.8671778e-05,
- -3.9337396e-05, 0, -4.2705009e-05,
- -4.1111101e-05, 0, -5.945824e-05,
- -2.8210339e-05, 0, -2.1934814e-05,
- 5.8917296e-05, 0, 3.8878991e-05,
- 7.5328415e-05, 0, 5.9727815e-05,
- 6.4888601e-05, 0, 3.2477254e-05,
- 1.0543302e-05, 0, -5.3633989e-06,
- -1.8523783e-05, 0, -2.8534835e-05,
- -1.8660178e-05, 0, -2.8629368e-05,
- 1.0359175e-05, 0, -5.4897177e-06,
- 6.4840587e-05, 0, 3.245778e-05,
- 7.5336564e-05, 0, 5.9731305e-05,
-}};
diff --git a/src/takin/convert.cpp b/src/takin/convert.cpp
deleted file mode 100644
index a30cb31..0000000
--- a/src/takin/convert.cpp
+++ /dev/null
@@ -1,163 +0,0 @@
-/**
- * converts version 1 grid data to version 2
- * @author Tobias Weber <tweber@ill.fr>
- * @date 05-jan-20
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#include <iostream>
-#include <fstream>
-#include <string>
-#include <unordered_map>
-
-using t_float_dst = double;
-
-
-int main()
-{
- const double eps = 1e-8;
- const double maxE = 1.5;
-
- // --------------------------------------------------------------------
- // modify these for each data file to process
- std::string filenameIdx = "/home/tw/tmp/skx/skxdyn.idx";
- std::string filenameDat = "/home/tw/tmp/skx/skxdyn.bin";
- std::string filenameNewDat = "/home/tw/tmp/skx/skxdyn.sqw";
-
- t_float_dst hstep = 0.0006;
- t_float_dst hmin = -0.03;
- t_float_dst hmax = 0.03 + hstep;
-
- t_float_dst kstep = 0.0006;
- t_float_dst kmin = -0.03;
- t_float_dst kmax = 0.03 + kstep;
-
- t_float_dst lstep = 0.0006;
- t_float_dst lmin = -0.03;
- t_float_dst lmax = 0.03 + lstep;
- // --------------------------------------------------------------------
-
-
- std::cout << "Converting data file ..." << std::endl;
-
- std::ifstream ifDat(filenameDat);
- std::ofstream ofDat(filenameNewDat);
-
- std::unordered_map<std::size_t, std::size_t> mapIndices;
-
- std::size_t removedBranches = 0;
-
-
- // write a dummy index file offset at the beginning (to be filled in later)
- std::size_t idx_offs = 0;
- ofDat.write((char*)&idx_offs, sizeof(idx_offs));
-
-
- // write data dimensions
- ofDat.write((char*)&hmin, sizeof(hmin));
- ofDat.write((char*)&hmax, sizeof(hmax));
- ofDat.write((char*)&hstep, sizeof(hstep));
-
- ofDat.write((char*)&kmin, sizeof(kmin));
- ofDat.write((char*)&kmax, sizeof(kmax));
- ofDat.write((char*)&kstep, sizeof(kstep));
-
- ofDat.write((char*)&lmin, sizeof(lmin));
- ofDat.write((char*)&lmax, sizeof(lmax));
- ofDat.write((char*)&lstep, sizeof(lstep));
-
-
- // header
- ofDat << "takin_grid_data_ver2|title:skxdyn_reseda|author:tweber@ill.fr|date:25/mar/2020";
-
-
- std::size_t iLoop = 0;
- while(1)
- {
- std::size_t idx = ifDat.tellg();
- std::size_t idxnew = ofDat.tellp();
-
- unsigned int numBranches = 0;
- ifDat.read((char*)&numBranches, sizeof(numBranches));
- if(ifDat.gcount() != sizeof(numBranches) || ifDat.eof())
- break;
-
- mapIndices.insert(std::make_pair(idx, idxnew));
-
- // write placeholder
- unsigned int numNewBranches = 0;
- ofDat.write((char*)&numNewBranches, sizeof(numNewBranches));
-
-
- for(unsigned int branch=0; branch<numBranches; ++branch)
- {
- double vals[4] = { 0, 0, 0, 0 };
- ifDat.read((char*)vals, sizeof(vals));
-
-
- t_float_dst w = t_float_dst(std::abs(vals[1])+std::abs(vals[2])+std::abs(vals[3]));
-
- if(w >= eps && std::abs(vals[0]) <= maxE)
- {
- t_float_dst E = t_float_dst(vals[0]);
- if(std::abs(E) < eps)
- E = 0.;
- if(std::abs(w) < eps)
- w = 0.;
-
- t_float_dst newvals[2] = { E, w };
- ofDat.write((char*)newvals, sizeof(newvals));
-
- ++numNewBranches;
- }
- else
- {
- ++removedBranches;
- }
- }
-
-
- // seek back and write real number of branches
- std::size_t lastIdx = ofDat.tellp();
- ofDat.seekp(idxnew, std::ios_base::beg);
- ofDat.write((char*)&numNewBranches, sizeof(numNewBranches));
- ofDat.seekp(lastIdx, std::ios_base::beg);
-
- ++iLoop;
- }
-
- ifDat.close();
-
- // update index at beginning
- idx_offs = ofDat.tellp();
- ofDat.seekp(0, std::ios_base::beg);
- ofDat.write((char*)&idx_offs, sizeof(idx_offs));
- ofDat.seekp(idx_offs, std::ios_base::beg);
-
- std::cout << removedBranches << " branches removed (w<eps || E<maxE)." << std::endl;
- std::cout << "\nConverting index file ..." << std::endl;
-
- std::ifstream ifIdx(filenameIdx);
- std::ofstream &ofIdx = ofDat;
-
- while(1)
- {
- std::size_t idx = 0;
-
- ifIdx.read((char*)&idx, sizeof(idx));
- if(ifIdx.gcount() != sizeof(idx) || ifIdx.eof())
- break;
-
- auto iter = mapIndices.find(idx);
- if(iter == mapIndices.end())
- {
- std::cerr << "Error: Index " << std::hex << idx << " was not found." << std::endl;
- continue;
- }
-
- std::size_t newidx = iter->second;
- ofIdx.write((char*)&newidx, sizeof(newidx));
- }
-
- return 0;
-}
diff --git a/src/takin/dump.cpp b/src/takin/dump.cpp
deleted file mode 100644
index 7cbc704..0000000
--- a/src/takin/dump.cpp
+++ /dev/null
@@ -1,48 +0,0 @@
-/**
- * dump calculated bin files for quick checking
- * @author Tobias Weber <tweber@ill.fr>
- * @date mar-20
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#include <fstream>
-#include <iostream>
-
-
-int main(int argc, char** argv)
-{
- if(argc<=1)
- {
- std::cerr << "Give a bin file" << std::endl;
- return -1;
- }
-
- std::ifstream ifstr(argv[1], std::ios_base::binary);
- if(!ifstr)
- {
- std::cerr << "Cannot read " << argv[1] << std::endl;
- return -1;
- }
-
-
- while(!ifstr.eof())
- {
- unsigned int num = 0;
- ifstr.read((char*)&num, sizeof(num));
- std::cout << "Number of branches: " << num << std::endl;
-
- for(unsigned int i=0; i<num; ++i)
- {
- double E=0, w1=0, w2=0, w3=0;
- ifstr.read((char*)&E, sizeof(E));
- ifstr.read((char*)&w1, sizeof(w1));
- ifstr.read((char*)&w2, sizeof(w2));
- ifstr.read((char*)&w3, sizeof(w3));
-
- std::cout << "#" << i << ": " << "E=" << E << "ws=[" << w1 << ", " << w2 << ", " << w3 << "]" << std::endl;
- }
- }
-
-
- return 0;
-}
diff --git a/src/takin/genskx.cpp b/src/takin/genskx.cpp
deleted file mode 100644
index 28d37da..0000000
--- a/src/takin/genskx.cpp
+++ /dev/null
@@ -1,154 +0,0 @@
-/**
- * Pre-calculates the dispersion and stores it in a database
- * @author Tobias Weber <tweber@ill.fr>
- * @date oct-18
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#include "core/skx.h"
-#include "tlibs2/libs/str.h"
-#include "tlibs2/libs/file.h"
-#include <fstream>
-
-using t_real = double;
-using t_cplx = std::complex<t_real>;
-using t_vec = ublas::vector<t_real>;
-const auto g_j = t_cplx(0,1);
-
-#include "core/skx_default_gs.cxx"
-
-
-// ranges
-const t_real qrange = 0.03;
-const t_real qstep = 0.0006;
-const t_real Erange = 0.5;
-
-
-void calc_disp(
- t_real qh,
- t_real Bx, t_real By, t_real Bz,
- t_real Px, t_real Py, t_real Pz,
- const std::string& filename)
-{
- Skx<t_real, t_cplx, DEF_SKX_ORDER> skx;
-
- std::vector<ublas::vector<t_cplx>> fourier_skx;
- fourier_skx.reserve(_skxgs_allcomps.size()/3);
-
- for(std::size_t comp=0; comp<_skxgs_allcomps.size(); comp+=3)
- fourier_skx.push_back(tl2::make_vec<ublas::vector<t_cplx>>({_skxgs_allcomps[comp], _skxgs_allcomps[comp+1], _skxgs_allcomps[comp+2]}));
-
- skx.SetFourier(fourier_skx);
- skx.SetProjNeutron(false);
- skx.SetT(-1000.);
- skx.SetB(25.); // BC2 = 40.3425
- skx.GenFullFourier();
- skx.SetFilterZeroWeight(true);
-
- //t_vec G = tl2::make_vec<t_vec>({ 1, 1, 0 });
- t_vec G = tl2::make_vec<t_vec>({ 0, 0, 0 });
- t_vec Pdir = tl2::make_vec<t_vec>({ Px, Py, Pz });
- t_vec Bdir = tl2::make_vec<t_vec>({ Bx, By, Bz });
-
- skx.SetCoords(Bdir[0],Bdir[1],Bdir[2], Pdir[0],Pdir[1],Pdir[2]);
- skx.SetG(G[0], G[1], G[2]);
-
-
-
- // outputs
- std::ofstream ofstrLog(filename + ".log");
-
- if(!ofstrLog)
- {
- std::cerr << "Could not write to output files \"" << filename << ".*\"" << std::endl;
- return;
- }
-
- ofstrLog << "h = " << qh << std::endl;
- ofstrLog << "k_start = " << -qrange << ", k_end = " << qrange << ", k_step = " << qstep << std::endl;
- ofstrLog << "l_start = " << -qrange << ", l_end = " << qrange << ", l_step = " << qstep << std::endl;
-
-
- int k_idx = 0;
- int skip_files = 0;
- for(t_real qk=-qrange; qk<qrange; qk+=qstep)
- {
- std::string thisfile = filename + "_" + tl2::var_to_str(k_idx);
- if(tl2::file_exists<char>((thisfile + ".idx").c_str()) || skip_files > 0)
- {
- std::cout << "Skipping " << thisfile << std::endl;
- ++k_idx;
- --skip_files;
- continue;
- }
-
- std::ofstream ofstrIdx(thisfile + ".idx", std::ios_base::binary); // format: [h1, k1, l1], idx1, [h2, k2, l2], idx2, ...
- std::ofstream ofstrBin(thisfile + ".bin", std::ios_base::binary); // format: num_Es, E1, w1_SF1, w1_SF2, w1_NSF, E2, ...
-
-
- for(t_real ql=-qrange; ql<qrange; ql+=qstep)
- {
- std::cout << "\rCalculating " << thisfile << " (" << qh << " " << qk << " " << ql << ") ... ";
- std::cout.flush();
-
- t_vec Q = G;
- Q[0] += qh; Q[1] += qk; Q[2] += ql;
-
- std::vector<t_real> Es, wsSF1, wsSF2, wsNSF;
- std::tie(Es, std::ignore, wsSF1, wsSF2, wsNSF) = skx.GetDisp(Q[0], Q[1], Q[2], -Erange, Erange);
-
- std::size_t offset = (std::size_t)ofstrBin.tellp();
- ofstrIdx.write((char*)&offset, sizeof(offset));
-
- unsigned int numEs = (unsigned int)Es.size();
- ofstrBin.write((char*)&numEs, sizeof(numEs));
-
- for(std::size_t idxE=0; idxE<Es.size(); ++idxE)
- {
- ofstrBin.write((char*)&Es[idxE], sizeof(t_real));
- ofstrBin.write((char*)&wsSF1[idxE], sizeof(t_real));
- ofstrBin.write((char*)&wsSF2[idxE], sizeof(t_real));
- ofstrBin.write((char*)&wsNSF[idxE], sizeof(t_real));
- }
-
-
- ofstrLog << "q = " << qh << " " << qk << " " << ql << ", "
- << "Q = " << Q[0] << " " << Q[1] << " " << Q[2] << ": "
- << Es.size() << " branches." << std::endl;
-
- ofstrIdx.flush();
- ofstrBin.flush();
- }
-
- std::cout << "done" << std::endl;
- ++k_idx;
- }
-}
-
-
-
-int main(int argc, char** argv)
-{
- std::cout
- << "--------------------------------------------------------------------------------\n"
- << "\tS(q,E) grid calculation tool,\n\t\tby T. Weber <tweber@ill.fr>, October 2018.\n"
- << "--------------------------------------------------------------------------------\n\n";
-
- if(argc<=1)
- {
- std::cerr << "Please give an index." << std::endl;
- return -1;
- }
-
- int idx = tl2::str_to_var<int>(std::string(argv[1]));
- //t_real qh = tl2::str_to_var<t_real>(std::string(argv[2]));
- t_real qh = -qrange + t_real(idx)*qstep;
-
- std::cout << "qh = " << qh << std::endl;
- std::cout << "idx = " << idx << std::endl;
-
- std::string filename = "skxdyn_" + tl2::var_to_str(idx);
- calc_disp(qh, 0,0,1, 1,-1,0, filename);
-
- return 0;
-}
diff --git a/src/takin/merge.cpp b/src/takin/merge.cpp
deleted file mode 100644
index ed6be61..0000000
--- a/src/takin/merge.cpp
+++ /dev/null
@@ -1,104 +0,0 @@
-/**
- * Merges the data files generated by genskx
- * @author Tobias Weber <tweber@ill.fr>
- * @date mar-20
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#include <iostream>
-#include <fstream>
-#include <string>
-
-#include "tlibs2/libs/file.h"
-#include "tlibs2/libs/str.h"
-
-
-int main()
-{
- // --------------------------------------------------------------------
- // modify these for each data file to process
- std::string filename = "skxdyn";
- constexpr std::size_t totaldirs = 100;
- constexpr std::size_t totalfiles = 100;
- constexpr std::size_t expectedidxsize = 808;
- // --------------------------------------------------------------------
-
-
- std::ofstream ofstrMergedIdx(filename + ".idx");
- std::ofstream ofstrMergedBin(filename + ".bin");
-
- std::size_t offs = 0;
- for(std::size_t diridx=0; diridx<=totaldirs; ++diridx)
- {
- for(std::size_t fileidx=0; fileidx<=totalfiles; ++fileidx)
- {
- std::string filename_idx = tl2::var_to_str(diridx) + "/" + filename + "_" + tl2::var_to_str(diridx) + "_" + tl2::var_to_str(fileidx) + ".idx";
- std::string filename_bin = tl2::var_to_str(diridx) + "/" + filename + "_" + tl2::var_to_str(diridx) + "_" + tl2::var_to_str(fileidx) + ".bin";
-
- std::cout << "Merging index file \"" << filename_idx << "\"..." << std::endl;
-
-
- // correct indices
- std::ifstream ifstrIdx(filename_idx);
- if(!ifstrIdx)
- {
- std::cerr << "Could not open \"" << filename_idx << "\"" << std::endl;
- return -1;
- }
-
- if(tl2::get_file_size(ifstrIdx) != expectedidxsize)
- {
- std::cerr << "Corrupted data. Dir: " << diridx << ", file: " << fileidx << std::endl;
- return -1;
- }
-
- // write new indices
- std::size_t sizeWritten=0;
- while(true)
- {
- std::size_t curidx = 0;
- ifstrIdx.read((char*)&curidx, sizeof(curidx));
- if(ifstrIdx.gcount() != sizeof(curidx) || ifstrIdx.eof())
- break;
-
- curidx += offs;
- ofstrMergedIdx.write((char*)&curidx, sizeof(curidx));
-
- sizeWritten += sizeof(curidx);
- }
- if(sizeWritten != expectedidxsize)
- {
- std::cerr << "Invalid number of indices written." << std::endl;
- return -1;
- }
-
-
- // add offset (sizes of all previous binary files) to indices
- std::ifstream ifstrBin(filename_bin);
- if(!ifstrBin)
- {
- std::cerr << "Could not open \"" << filename_bin << "\"" << std::endl;
- return -1;
- }
-
- std::size_t sizeBin = tl2::get_file_size(ifstrBin);
- offs += sizeBin;
-
-
- std::cout << "Merging binary file \"" << filename_bin << "\"..." << std::endl;
- char* data = new char[sizeBin];
- memset(data, 0, sizeBin);
- ifstrBin.read(data, sizeBin);
- if(ifstrBin.gcount() != sizeBin)
- {
- std::cerr << "Could not merge binary file." << std::endl;
- return -1;
- }
-
- ofstrMergedBin.write(data, sizeBin);
- delete[] data;
- }
- }
-
- return 0;
-}
diff --git a/src/takin/skxgrid.py b/src/takin/skxgrid.py
deleted file mode 100644
index f2b46fd..0000000
--- a/src/takin/skxgrid.py
+++ /dev/null
@@ -1,259 +0,0 @@
-#
-# skx grid
-# @author tweber@ill.fr
-# @date 4-feb-2020
-# @license GPLv2 (see 'LICENSE' file)
-#
-
-import numpy as np
-
-
-# -----------------------------------------------------------------------------
-# Test files and configuration
-
-use_thales = True
-use_reseda = False
-
-if use_thales:
- idxfile = "/home/tw/tmp/skx/skxdyn_asym.idx"
- datafile = "/home/tw/tmp/skx/skxdyn_asym.bin"
-
- hstep = 0.001
- hmin = -0.096
- hmax = 0.096
-
- kstep = 0.001
- kmin = -0.096
- kmax = 0.096
-
- lstep = 0.001
- lmin = -0.096
- lmax = 0.096
-
- xlim1 = -0.09
- xlim2 = 0.09
- ylim1 = -1.
- ylim2 = 1.
-
-
-if use_reseda:
- idxfile = "/home/tw/tmp/skx/skxdyn_reseda.idx"
- datafile = "/home/tw/tmp/skx/skxdyn_reseda.bin"
-
- hstep = 0.0006
- hmin = -0.03
- hmax = 0.03 + hstep
-
- kstep = 0.0006
- kmin = -0.03
- kmax = 0.03 + kstep
-
- lstep = 0.0006
- lmin = -0.03
- lmax = 0.03 + lstep
-
- xlim1 = -0.03
- xlim2 = 0.03
- ylim1 = -0.3
- ylim2 = 0.3
-
-
-numvals = 4 # [E, w1, w2, w3]
-
-
-# longitudinal
-#plot_hklbegin = np.array([-0.09, -0.09, -0.])
-#plot_hklend = np.array([0.09, 0.09, 0.])
-#plot_dir = 0
-
-# transversal
-plot_hklbegin = np.array([-0.09, 0.09, 0.])
-plot_hklend = np.array([0.09, -0.09, 0.])
-plot_dir = 0
-
-#plot_hklbegin = np.array([0., -0.03, 0.])
-#plot_hklend = np.array([0., 0.03, 0.])
-#plot_dir = 1
-
-#plot_hklbegin = np.array([0.03, -0.03, -0.])
-#plot_hklend = np.array([-0.03, 0.03, 0.])
-#plot_dir = 0
-
-plot_hklsteps = 512
-# -----------------------------------------------------------------------------
-
-
-
-# using q, not Q
-def hkl_to_idx(hkl):
- [h, k, l] = hkl
-
- # clamp values to boundaries
- if h < hmin:
- h = hmin
- if k < kmin:
- k = kmin
- if l < lmin:
- l = lmin
- if h >= hmax:
- h = hmax - hstep
- if k >= kmax:
- k = kmax - kstep
- if l >= lmax:
- l = lmax - lstep
-
- # max dimensions
- iHSize = int(round((hmax-hmin) / hstep))
- iKSize = int(round((kmax-kmin) / kstep))
- iLSize = int(round((lmax-lmin) / lstep))
-
- # position indices
- iH = int(round(((h - hmin) / hstep)))
- iK = int(round(((k - kmin) / kstep)))
- iL = int(round(((l - lmin) / lstep)))
-
- # clamp again
- if iH >= iHSize:
- iH = iHSize-1
- if iK >= iKSize:
- iK = iKSize-1
- if iL >= iLSize:
- iL = iLSize-1
- if iH < 0:
- iH = 0
- if iK < 0:
- iK = 0
- if iL < 0:
- iL = 0
-
- return int(iH*iKSize*iLSize + iK*iLSize + iL)
-
-
-
-def getE(idxfilehandle, datafilehandle, hkl):
- hklidx = hkl_to_idx(hkl)
- #print("Index into index file: %d." % hklidx)
-
- idx = np.memmap(idxfilehandle, dtype="uint64", mode="r", offset=int(hklidx*8))[0]
- #print("Index into data file: %d." % idx)
-
- data = np.memmap(datafilehandle, dtype="uint32", mode="r", offset=int(idx))
- numbranches = data[0]
- #print("Number of dispersion branches: %d." % numbranches)
-
- values = np.ndarray.view(data[1 : 1+numbranches*numvals*2], dtype="float64")
- values.shape = (numbranches, numvals)
- return values
-
-
-
-
-# plot an example dispersion
-def plot_disp(idxfilehandle, datafilehandle, hklbegin, hklend, hklsteps):
- import matplotlib.pyplot as plot
- symscale = 2.
- eps = 1e-8
-
- qs_h = []
- qs_k = []
- qs_l = []
- Es = []
- wsSF1 = []
- wsSF2 = []
- wsNSF = []
-
- for step in range(0, hklsteps+1):
- hkl = hklbegin + (hklend-hklbegin)/hklsteps*step
- branches = getE(_idxfilehandle, _datafilehandle, hkl)
-
- newEs = [ branch[0] for branch in branches ]
- Es.extend(newEs)
- wsSF1.extend([ branch[1] if branch[1] > eps else 0. for branch in branches ])
- wsSF2.extend([ branch[2] if branch[2] > eps else 0. for branch in branches ])
- wsNSF.extend([ branch[3] if branch[3] > eps else 0. for branch in branches ])
- qs_h.extend([hkl[0]] * len(newEs))
- qs_k.extend([hkl[1]] * len(newEs))
- qs_l.extend([hkl[2]] * len(newEs))
-
- # convert to np
- qs_h = np.array(qs_h)
- qs_k = np.array(qs_k)
- qs_l = np.array(qs_l)
- Es = np.array(Es)
- wsSF1 = np.abs(np.array(wsSF1))
- wsSF2 = np.abs(np.array(wsSF2))
- wsNSF = np.abs(np.array(wsNSF))
- wsTotal = wsNSF + wsSF1 + wsSF2
-
- fig = plot.figure()
-
- pltSF1 = fig.add_subplot(221)
- pltSF1.set_xlabel("q (rlu)")
- pltSF1.set_ylabel("E (meV)")
- pltSF1.set_xlim(xlim1, xlim2)
- pltSF1.set_ylim(ylim1, ylim2)
- pltSF1.set_title("SF1")
- if plot_dir == 0:
- pltSF1.scatter(qs_h, Es, marker=".", s=wsSF1*symscale)
- elif plot_dir == 1:
- pltSF1.scatter(qs_k, Es, marker=".", s=wsSF1*symscale)
- else:
- pltSF1.scatter(qs_l, Es, marker=".", s=wsSF1*symscale)
-
- pltSF2 = fig.add_subplot(222)
- pltSF2.set_xlabel("q (rlu)")
- pltSF2.set_ylabel("E (meV)")
- pltSF2.set_xlim(xlim1, xlim2)
- pltSF2.set_ylim(ylim1, ylim2)
- pltSF2.set_title("SF2")
- if plot_dir == 0:
- pltSF2.scatter(qs_h, Es, marker=".", s=wsSF2*symscale)
- elif plot_dir == 1:
- pltSF2.scatter(qs_k, Es, marker=".", s=wsSF2*symscale)
- else:
- pltSF2.scatter(qs_l, Es, marker=".", s=wsSF2*symscale)
-
- pltNSF = fig.add_subplot(223)
- pltNSF.set_xlabel("q (rlu)")
- pltNSF.set_ylabel("E (meV)")
- pltNSF.set_xlim(xlim1, xlim2)
- pltNSF.set_ylim(ylim1, ylim2)
- pltNSF.set_title("NSF")
- if plot_dir == 0:
- pltNSF.scatter(qs_h, Es, marker=".", s=wsNSF*symscale)
- elif plot_dir == 1:
- pltNSF.scatter(qs_k, Es, marker=".", s=wsNSF*symscale)
- else:
- pltNSF.scatter(qs_l, Es, marker=".", s=wsNSF*symscale)
-
- pltTotal = fig.add_subplot(224)
- pltTotal.set_xlabel("q (rlu)")
- pltTotal.set_ylabel("E (meV)")
- pltTotal.set_xlim(xlim1, xlim2)
- pltTotal.set_ylim(ylim1, ylim2)
- pltTotal.set_title("Total")
- if plot_dir == 0:
- pltTotal.scatter(qs_h, Es, marker=".", s=wsTotal*symscale)
- elif plot_dir == 1:
- pltTotal.scatter(qs_k, Es, marker=".", s=wsTotal*symscale)
- else:
- pltTotal.scatter(qs_l, Es, marker=".", s=wsTotal*symscale)
-
- plot.tight_layout()
- plot.show()
-
-
-
-
-# open index and data file for mapping
-try:
- _idxfilehandle = open(idxfile, "rb")
- _datafilehandle = open(datafile, "rb")
-except err as IOError:
- print(err)
- exit(-1)
-
-#print(getE(_idxfilehandle, _datafilehandle, [-0.05, -0.05, 0.]))
-#exit()
-
-plot_disp(_idxfilehandle, _datafilehandle, plot_hklbegin, plot_hklend, plot_hklsteps)
diff --git a/src/takin/takin.cpp b/src/takin/takin.cpp
deleted file mode 100644
index 047ca6c..0000000
--- a/src/takin/takin.cpp
+++ /dev/null
@@ -1,439 +0,0 @@
-/**
- * Takin module
- * @author Tobias Weber <tweber@ill.fr>
- * @date sep-18
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#ifdef PLUGIN_APPLI
- #include "takin/tools/monteconvo/sqw_proc.h"
- #include "takin/tools/monteconvo/sqw_proc_impl.h"
-#else
- #include <boost/dll/alias.hpp>
-#endif
-
-#include "takin/tools/monteconvo/sqwbase.h"
-#include "takin/libs/version.h"
-
-#include "core/skx.h"
-#include "core/fp.h"
-#include "core/heli.h"
-#include "core/longfluct.h"
-
-#include "tlibs2/libs/phys.h"
-#include "tlibs2/libs/str.h"
-
-
-#ifndef DEF_SKX_ORDER
- #define SKX_ORDER 4
-#else
- #define SKX_ORDER DEF_SKX_ORDER
-#endif
-
-#ifndef DEF_HELI_ORDER
- #define HELI_ORDER 4
-#else
- #define HELI_ORDER DEF_HELI_ORDER
-#endif
-
-
-class SqwMod : public SqwBase
-{
-public:
- using SqwBase::t_var;
- using t_real = t_real_reso;
- using t_vec = ublas::vector<t_real>;
- using t_cplx = std::complex<t_real>;
-
-protected:
- Skx<t_real, t_cplx, SKX_ORDER> m_skx;
- FP<t_real, t_cplx> m_fp;
- Heli<t_real, t_cplx, HELI_ORDER> m_heli;
- Longfluct m_lf;
-
- t_real m_dErange = 2.; // E range around queried E to take into account
- t_real m_dSigma = t_real(0.05);
- t_real m_dS0 = t_real(1.);
- t_real m_dIncAmp = t_real(0.);
- t_real m_dIncSigma = t_real(0.05);
- t_real m_dT = 29.;
- t_real m_dB = 0.35;
- t_real m_dcut = 0.02;
-
- int m_iOnlyMode = -1;
- int m_iPolChan = 0;
- int m_iwhich_disp = 0; // 0: skx, 1: fp, 2: heli
- int m_iProjNeutron = 1;
- int m_ionlylf = 0;
-
- t_vec m_vecG = tl2::make_vec<t_vec>({1,1,0});
- t_vec m_vecB = tl2::make_vec<t_vec>({1,1,0});
- t_vec m_vecPin = tl2::make_vec<t_vec>({1,-1,0});
-
-public:
- SqwMod();
- SqwMod(const std::string& strCfgFile);
- virtual ~SqwMod();
-
- virtual std::tuple<std::vector<t_real>, std::vector<t_real>> disp(t_real dh, t_real dk, t_real dl) const override;
- virtual t_real operator()(t_real dh, t_real dk, t_real dl, t_real dE) const override;
-
- virtual std::vector<t_var> GetVars() const override;
- virtual void SetVars(const std::vector<t_var>&) override;
- virtual bool SetVarIfAvail(const std::string& strKey, const std::string& strNewVal) override;
-
- virtual SqwBase* shallow_copy() const override;
-};
-
-
-using t_real = typename SqwMod::t_real;
-#include "core/skx_default_gs.cxx"
-
-
-// ----------------------------------------------------------------------------
-SqwMod::SqwMod()
-{
- tl2::log_info("--------------------------------------------------------------------------------");
- tl2::log_info("This is the Takin MnSi magnon dynamics module,");
- tl2::log_info("\tby T. Weber <tweber@ill.fr>, September 2018.");
- tl2::log_info("--------------------------------------------------------------------------------");
-
- m_skx.SetT(-1000.);
- m_skx.SetB(25.);
- m_fp.SetT(29.);
- m_fp.SetB(0.35);
- m_heli.SetT(29.);
- m_heli.SetB(0.35);
- m_lf.SetT(29.);
-
- std::vector<ublas::vector<t_cplx>> fourier_skx;
- fourier_skx.reserve(_skxgs_allcomps.size()/3);
-
- for(std::size_t comp=0; comp<_skxgs_allcomps.size(); comp+=3)
- fourier_skx.push_back(tl2::make_vec<ublas::vector<t_cplx>>({_skxgs_allcomps[comp], _skxgs_allcomps[comp+1], _skxgs_allcomps[comp+2]}));
-
- m_skx.SetFourier(fourier_skx);
- m_skx.GenFullFourier();
-
- m_skx.SetCoords(m_vecB[0],m_vecB[1],m_vecB[2], m_vecPin[0],m_vecPin[1],m_vecPin[2]);
- m_skx.SetG(m_vecG[0], m_vecG[1], m_vecG[2]);
- m_skx.SetFilterZeroWeight(true);
-
- m_fp.SetCoords(m_vecB[0],m_vecB[1],m_vecB[2]);
- m_fp.SetG(m_vecG[0], m_vecG[1], m_vecG[2]);
-
- m_heli.SetCoords(m_vecB[0],m_vecB[1],m_vecB[2]);
- m_heli.SetG(m_vecG[0], m_vecG[1], m_vecG[2]);
- m_heli.SetFilterZeroWeight(true);
-
- m_lf.SetPinning(m_vecPin[0],m_vecPin[1],m_vecPin[2], m_vecB[0],m_vecB[1],m_vecB[2]);
-
- SqwBase::m_bOk = 1;
-}
-
-SqwMod::SqwMod(const std::string& strCfgFile) : SqwMod() { SqwBase::m_bOk = 1; }
-SqwMod::~SqwMod() {}
-// ----------------------------------------------------------------------------
-
-
-
-// ----------------------------------------------------------------------------
-// S(Q,E)
-
-std::tuple<std::vector<t_real>, std::vector<t_real>>
- SqwMod::disp(t_real dh, t_real dk, t_real dl) const
-{
- std::vector<t_real> vecE, vecW[4];
- if(m_iwhich_disp == 0)
- std::tie(vecE, vecW[0], vecW[1], vecW[2], vecW[3]) = m_skx.GetDisp(dh, dk, dl);
- else if(m_iwhich_disp == 1)
- std::tie(vecE, vecW[0], vecW[1], vecW[2], vecW[3]) = m_fp.GetDisp(dh, dk, dl);
- else if(m_iwhich_disp == 2)
- std::tie(vecE, vecW[0], vecW[1], vecW[2], vecW[3]) = m_heli.GetDisp(dh, dk, dl);
-
- return std::make_tuple(vecE, vecW[m_iPolChan]);
-}
-
-
-t_real SqwMod::operator()(t_real dh, t_real dk, t_real dl, t_real dE) const
-{
- t_vec vecq = tl2::make_vec<t_vec>({dh, dk, dl}) - m_vecG;
-
- // longitudinal fluctuations
- t_real dS_lf = m_lf.S_para(vecq, dE);
- if(m_ionlylf)
- return dS_lf;
-
- std::vector<t_real> vecE, vecW[4];
- // only calculate dispersion if global weight factor is not 0
- if(!tl2::float_equal(m_dS0, t_real(0)))
- {
- if(m_iwhich_disp == 0)
- std::tie(vecE, vecW[0], vecW[1], vecW[2], vecW[3]) = m_skx.GetDisp(dh, dk, dl, dE-m_dErange, dE+m_dErange);
- else if(m_iwhich_disp == 1)
- std::tie(vecE, vecW[0], vecW[1], vecW[2], vecW[3]) = m_fp.GetDisp(dh, dk, dl, dE-m_dErange, dE+m_dErange);
- else if(m_iwhich_disp == 2)
- std::tie(vecE, vecW[0], vecW[1], vecW[2], vecW[3]) = m_heli.GetDisp(dh, dk, dl, dE-m_dErange, dE+m_dErange);
- }
-
- // incoherent
- t_real dInc=0, dS_p=0, dS_m=0;
- if(!tl2::float_equal(m_dIncAmp, t_real(0)))
- dInc = tl2::gauss_model(dE, t_real(0), m_dIncSigma, m_dIncAmp, t_real(0));
-
- t_real dS = 0;
- for(std::size_t iE=0; iE<vecE.size(); ++iE)
- {
- if(!tl2::float_equal(vecW[m_iPolChan][iE], t_real(0)))
- dS += tl2::gauss_model(dE, vecE[iE], m_dSigma, vecW[m_iPolChan][iE], t_real(0));
- }
-
- return m_dS0*dS * tl2::bose_cutoff(dE, m_dT, m_dcut) + dInc;
-}
-// ----------------------------------------------------------------------------
-
-
-
-// ----------------------------------------------------------------------------
-// getters & setters
-
-std::vector<SqwMod::t_var> SqwMod::GetVars() const
-{
- std::vector<t_var> vecVars;
-
- vecVars.push_back(SqwBase::t_var{"proj_neutron", "real", tl2::var_to_str(m_iProjNeutron)});
- vecVars.push_back(SqwBase::t_var{"which_disp", "real", tl2::var_to_str(m_iwhich_disp)});
- vecVars.push_back(SqwBase::t_var{"E_range", "real", tl2::var_to_str(m_dErange)});
- vecVars.push_back(SqwBase::t_var{"bose_cutoff", "real", tl2::var_to_str(m_dcut)});
- vecVars.push_back(SqwBase::t_var{"only_mode", "int", tl2::var_to_str(m_iOnlyMode)});
- vecVars.push_back(SqwBase::t_var{"sigma", "real", tl2::var_to_str(m_dSigma)});
- vecVars.push_back(SqwBase::t_var{"inc_amp", "real", tl2::var_to_str(m_dIncAmp)});
- vecVars.push_back(SqwBase::t_var{"inc_sigma", "real", tl2::var_to_str(m_dIncSigma)});
- vecVars.push_back(SqwBase::t_var{"S0", "real", tl2::var_to_str(m_dS0)});
- vecVars.push_back(SqwBase::t_var{"T", "real", tl2::var_to_str(m_dT)});
- vecVars.push_back(SqwBase::t_var{"B", "real", tl2::var_to_str(m_dB)});
- vecVars.push_back(SqwBase::t_var{"pol_chan", "int", tl2::var_to_str(m_iPolChan)});
- vecVars.push_back(SqwBase::t_var{"G", "vector", vec_to_str(m_vecG)});
- vecVars.push_back(SqwBase::t_var{"B_dir", "vector", vec_to_str(m_vecB)});
- vecVars.push_back(SqwBase::t_var{"Pin_dir", "vector", vec_to_str(m_vecPin)});
-
- vecVars.push_back(SqwBase::t_var{"lf_only", "int", tl2::var_to_str(m_ionlylf)});
- vecVars.push_back(SqwBase::t_var{"lf_invcorrel", "real", tl2::var_to_str(m_lf.GetInvCorrel())});
- vecVars.push_back(SqwBase::t_var{"lf_A", "real", tl2::var_to_str(m_lf.GetA())});
- vecVars.push_back(SqwBase::t_var{"lf_gamma", "real", tl2::var_to_str(m_lf.GetGamma())});
-
- return vecVars;
-}
-
-
-void SqwMod::SetVars(const std::vector<SqwMod::t_var>& vecVars)
-{
- if(!vecVars.size()) return;
-
- for(const SqwBase::t_var& var : vecVars)
- {
- const std::string& strVar = std::get<0>(var);
- const std::string& strVal = std::get<2>(var);
-
- if(strVar == "proj_neutron") m_iProjNeutron = tl2::str_to_var<int>(strVal);
- else if(strVar == "which_disp") m_iwhich_disp = tl2::str_to_var<int>(strVal);
- else if(strVar == "E_range") m_dErange = tl2::str_to_var<t_real>(strVal);
- else if(strVar == "bose_cutoff") m_dcut = tl2::str_to_var<t_real>(strVal);
- else if(strVar == "sigma") m_dSigma = tl2::str_to_var<t_real>(strVal);
- else if(strVar == "inc_amp") m_dIncAmp = tl2::str_to_var<decltype(m_dIncAmp)>(strVal);
- else if(strVar == "inc_sigma") m_dIncSigma = tl2::str_to_var<decltype(m_dIncSigma)>(strVal);
- else if(strVar == "S0") m_dS0 = tl2::str_to_var<decltype(m_dS0)>(strVal);
- else if(strVar == "only_mode")
- {
- m_iOnlyMode = tl2::str_to_var<int>(strVal);
- m_heli.SetOnlyMode(m_iOnlyMode);
- }
- else if(strVar == "T")
- {
- m_dT = tl2::str_to_var<decltype(m_dT)>(strVal);
-
- m_fp.SetT(m_dT);
- m_heli.SetT(m_dT);
- // fixed (and theo units) for skx!
-
- m_lf.SetT(m_dT);
- }
- else if(strVar == "B")
- {
- m_dB = tl2::str_to_var<decltype(m_dB)>(strVal);
-
- m_fp.SetB(m_dB);
- m_heli.SetB(m_dB);
- // fixed (and theo units) for skx!
- }
- else if(strVar == "pol_chan") m_iPolChan = tl2::str_to_var<decltype(m_iPolChan)>(strVal);
- else if(strVar == "G" || strVal == "proj_neutron")
- {
- m_vecG = str_to_vec<decltype(m_vecG)>(strVal);
-
- m_skx.SetProjNeutron(m_iProjNeutron!=0);
- m_skx.SetG(m_vecG[0], m_vecG[1], m_vecG[2]);
-
- m_fp.SetProjNeutron(m_iProjNeutron!=0);
- m_fp.SetG(m_vecG[0], m_vecG[1], m_vecG[2]);
-
- m_heli.SetProjNeutron(m_iProjNeutron!=0);
- m_heli.SetG(m_vecG[0], m_vecG[1], m_vecG[2]);
- }
- else if(strVar == "B_dir")
- {
- m_vecB = str_to_vec<decltype(m_vecB)>(strVal);
-
- m_skx.SetCoords(m_vecB[0],m_vecB[1],m_vecB[2], m_vecPin[0],m_vecPin[1],m_vecPin[2]);
- m_skx.SetG(m_vecG[0], m_vecG[1], m_vecG[2]);
-
- m_fp.SetCoords(m_vecB[0],m_vecB[1],m_vecB[2]);
- m_fp.SetG(m_vecG[0], m_vecG[1], m_vecG[2]);
-
- m_heli.SetCoords(m_vecB[0],m_vecB[1],m_vecB[2]);
- m_heli.SetG(m_vecG[0], m_vecG[1], m_vecG[2]);
-
- m_lf.SetPinning(m_vecPin[0],m_vecPin[1],m_vecPin[2],
- m_vecB[0],m_vecB[1],m_vecB[2]);
- }
- else if(strVar == "Pin_dir")
- {
- m_vecPin = str_to_vec<decltype(m_vecPin)>(strVal);
-
- m_skx.SetCoords(m_vecB[0],m_vecB[1],m_vecB[2], m_vecPin[0],m_vecPin[1],m_vecPin[2]);
- m_skx.SetG(m_vecG[0], m_vecG[1], m_vecG[2]);
-
- m_lf.SetPinning(m_vecPin[0],m_vecPin[1],m_vecPin[2],
- m_vecB[0],m_vecB[1],m_vecB[2]);
- }
-
- else if(strVar == "lf_only")
- {
- m_ionlylf = tl2::str_to_var<decltype(m_ionlylf)>(strVal);
- }
- else if(strVar == "lf_invcorrel")
- {
- t_real dInvCorrel = tl2::str_to_var<t_real>(strVal);
- m_lf.SetInvCorrel(dInvCorrel);
- }
- else if(strVar == "lf_A")
- {
- t_real dA = tl2::str_to_var<t_real>(strVal);
- m_lf.SetA(dA);
- }
- else if(strVar == "lf_gamma")
- {
- t_real dG = tl2::str_to_var<t_real>(strVal);
- m_lf.SetGamma(dG);
- }
- }
-}
-
-bool SqwMod::SetVarIfAvail(const std::string& strKey, const std::string& strNewVal)
-{
- return SqwBase::SetVarIfAvail(strKey, strNewVal);
-}
-// ----------------------------------------------------------------------------
-
-
-
-// ----------------------------------------------------------------------------
-// copy
-
-SqwBase* SqwMod::shallow_copy() const
-{
- SqwMod *pMod = new SqwMod();
-
- pMod->m_iProjNeutron = this->m_iProjNeutron;
- pMod->m_iwhich_disp = this->m_iwhich_disp;
- pMod->m_dSigma = this->m_dSigma;
- pMod->m_dIncAmp = this->m_dIncAmp;
- pMod->m_dIncSigma = this->m_dIncSigma;
- pMod->m_dS0 = this->m_dS0;
- pMod->m_dT = this->m_dT;
- pMod->m_dB = this->m_dB;
- pMod->m_skx = this->m_skx;
- pMod->m_fp = this->m_fp;
- pMod->m_heli = this->m_heli;
- pMod->m_lf = this->m_lf;
- pMod->m_ionlylf = this->m_ionlylf;
- pMod->m_iPolChan = this->m_iPolChan;
- pMod->m_vecG = this->m_vecG;
- pMod->m_vecB = this->m_vecB;
- pMod->m_vecPin = this->m_vecPin;
- pMod->m_dErange = this->m_dErange;
- pMod->m_dcut = this->m_dcut;
-
- return pMod;
-}
-// ----------------------------------------------------------------------------
-
-
-
-// ----------------------------------------------------------------------------
-// takin interface
-
-static const char* pcModIdent = "skxmod";
-static const char* pcModName = "MnSi Magnon Dynamics";
-
-
-#ifndef PLUGIN_APPLI
-
-
-// exported symbols
-std::tuple<std::string, std::string, std::string> sqw_info()
-{
- return std::make_tuple(TAKIN_VER, pcModIdent, pcModName);
-}
-
-std::shared_ptr<SqwBase> sqw_construct(const std::string& strCfgFile)
-{
- return std::make_shared<SqwMod>(strCfgFile);
-}
-
-
-// exports from so file
-#ifndef __MINGW32__
- BOOST_DLL_ALIAS(sqw_info, takin_sqw_info);
- BOOST_DLL_ALIAS(sqw_construct, takin_sqw);
-#else
- // hack because BOOST_DLL_ALIAS does not seem to work with Mingw
-
- extern "C" __declspec(dllexport)
- std::tuple<std::string, std::string, std::string> takin_sqw_info()
- {
- return sqw_info();
- }
-
- extern "C" __declspec(dllexport)
- std::shared_ptr<SqwBase> takin_sqw(const std::string& strCfgFile)
- {
- return sqw_construct(strCfgFile);
- }
-#endif
-
-
-#else
-
-
-int main(int argc, char** argv)
-{
- if(argc <= 2)
- {
- std::cout << "#\n# This is a Takin plugin module.\n#\n";
- std::cout << "module_ident: " << pcModIdent << "\n";
- std::cout << "module_name: " << pcModName << "\n";
- std::cout << "required_takin_version: " << TAKIN_VER << "\n";
- std::cout.flush();
- return 0;
- }
-
- const char* pcCfgFile = argv[1];
- const char* pcSharedMem = argv[2];
- SqwProc<SqwMod> proc(pcCfgFile, SqwProcStartMode::START_CHILD, pcSharedMem);
- return 0;
-}
-
-
-#endif
-// ----------------------------------------------------------------------------
diff --git a/src/takin/takin_grid.cpp b/src/takin/takin_grid.cpp
deleted file mode 100644
index d4834b0..0000000
--- a/src/takin/takin_grid.cpp
+++ /dev/null
@@ -1,492 +0,0 @@
-/**
- * Takin module for precalculated grid data
- * @author Tobias Weber <tweber@ill.fr>
- * @date nov-18
- * @license GPLv2 (see 'LICENSE' file)
- */
-
-#include "tools/monteconvo/sqwbase.h"
-#include "libs/version.h"
-
-#include "tlibs2/libs/math17.h"
-#include "tlibs2/libs/str.h"
-#include "tlibs2/libs/phys.h"
-#include "tlibs2/libs/file.h"
-
-#include <boost/dll/alias.hpp>
-#include <QtCore/QFile>
-
-
-using t_real = t_real_reso;
-using t_vec = ublas::vector<t_real>;
-using t_mat = ublas::matrix<t_real>;
-using t_quat = boost::math::quaternion<t_real>;
-using t_cplx = std::complex<t_real>;
-
-
-class SqwMod : public SqwBase
-{
- public:
- using SqwBase::t_var;
-
- protected:
- t_real m_dSigma = t_real(0.05);
- t_real m_dS0 = t_real(1.);
- t_real m_dIncAmp = t_real(0.);
- t_real m_dIncSigma = t_real(0.05);
- t_real m_dT = 29.;
- t_real m_dRotZ = 0.;
- t_real m_dRotX = 0.;
- t_real m_dcut = 0.02;
- int m_iPolChan = 0;
- bool m_bFlipB_or_q = false;
- t_real m_dEFactor = 0.24;
- t_real m_dEFactorGrid = 0.24;
-
- t_vec m_vecG = tl2::make_vec<t_vec>({1,1,0});
-
- t_vec m_vecRotFrom = tl2::make_vec<t_vec>({1,0,0});
- t_vec m_vecRotTo = tl2::make_vec<t_vec>({1,0,0});
-
- // grid descriptors
- std::string m_strIndexFile, m_strDataFile;
-
- t_real m_hmin=0., m_hmax=0., m_hstep=0.;
- t_real m_kmin=0., m_kmax=0., m_kstep=0.;
- t_real m_lmin=0., m_lmax=0., m_lstep=0.;
-
-
- public:
- SqwMod();
- SqwMod(const std::string& strCfgFile);
-
- virtual std::tuple<std::vector<t_real>, std::vector<t_real>>
- disp(t_real dh, t_real dk, t_real dl) const override;
- virtual t_real operator()(t_real dh, t_real dk, t_real dl, t_real dE) const override;
-
- virtual std::vector<t_var> GetVars() const override;
- virtual void SetVars(const std::vector<t_var>&) override;
- virtual bool SetVarIfAvail(const std::string& strKey, const std::string& strNewVal) override;
-
- virtual SqwBase* shallow_copy() const override;
-};
-
-
-
-// ----------------------------------------------------------------------------
-// constructors
-
-SqwMod::SqwMod()
-{
- SqwBase::m_bOk = 1;
-}
-
-
-SqwMod::SqwMod(const std::string& strCfgFile) : SqwMod()
-{
- tl2::log_info("--------------------------------------------------------------------------------");
- tl2::log_info("This is the Takin skx grid module,");
- tl2::log_info("\tby T. Weber <tweber@ill.fr>, November 2018.");
- tl2::log_info("--------------------------------------------------------------------------------");
-
- tl2::log_info("Grid description file: \"", strCfgFile, "\".");
-
- tl2::Prop<std::string> prop;
- if(prop.Load(strCfgFile.c_str(), tl2::PropType::INFO))
- {
- m_strIndexFile = prop.Query<std::string>("files/index");
- m_strDataFile = prop.Query<std::string>("files/data");
-
- m_hmin = prop.QueryAndParse<t_real>("dims/hmin");
- m_hmax = prop.QueryAndParse<t_real>("dims/hmax");
- m_hstep = prop.QueryAndParse<t_real>("dims/hstep");
- m_kmin = prop.QueryAndParse<t_real>("dims/kmin");
- m_kmax = prop.QueryAndParse<t_real>("dims/kmax");
- m_kstep = prop.QueryAndParse<t_real>("dims/kstep");
- m_lmin = prop.QueryAndParse<t_real>("dims/lmin");
- m_lmax = prop.QueryAndParse<t_real>("dims/lmax");
- m_lstep = prop.QueryAndParse<t_real>("dims/lstep");
-
- SqwBase::m_bOk = 1;
-
- tl2::log_info("Index file: ", m_strIndexFile, ", data file: ", m_strDataFile, ".");
- tl2::log_info("h range: ", m_hmin, " .. ", m_hmax, ", step: ", m_hstep);
- tl2::log_info("k range: ", m_kmin, " .. ", m_kmax, ", step: ", m_kstep);
- tl2::log_info("l range: ", m_lmin, " .. ", m_lmax, ", step: ", m_lstep);
- }
- else
- {
- tl2::log_err("Grid description file \"", strCfgFile, "\" could not be loaded.");
- SqwBase::m_bOk = 0;
- }
-}
-
-
-
-// ----------------------------------------------------------------------------
-// dispersion, spectral weight and structure factor
-
-std::tuple<std::vector<t_real>, std::vector<t_real>>
- SqwMod::disp(t_real dh, t_real dk, t_real dl) const
-{
- if(!tl2::vec_equal(m_vecRotFrom, m_vecRotTo, 1e-6))
- {
- t_quat quat = tl2::rotation_quat(m_vecRotFrom, m_vecRotTo);
- t_vec vecPos = tl2::make_vec<t_vec>({dh, dk, dl}) - m_vecG;
- vecPos = tl2::quat_vec_prod(quat, vecPos) + m_vecG;
- dh = vecPos[0];
- dk = vecPos[1];
- dl = vecPos[2];
- }
-
- if(!tl2::float_equal(m_dRotZ, 0., 1e-6))
- {
- t_mat rot = tl2::rotation_matrix_3d_z<t_mat>(m_dRotZ/180.*tl2::pi<t_real>);
- t_vec vecPos = tl2::make_vec<t_vec>({dh, dk, dl}) - m_vecG;
- vecPos = tl2::prod_mv(rot, vecPos) + m_vecG;
-
- dh = vecPos[0];
- dk = vecPos[1];
- dl = vecPos[2];
- }
-
- if(!tl2::float_equal(m_dRotX, 0., 1e-6))
- {
- t_mat rot = tl2::rotation_matrix_3d_x<t_mat>(m_dRotZ/180.*tl2::pi<t_real>);
- t_vec vecPos = tl2::make_vec<t_vec>({dh, dk, dl}) - m_vecG;
- vecPos = tl2::prod_mv(rot, vecPos) + m_vecG;
-
- dh = vecPos[0];
- dk = vecPos[1];
- dl = vecPos[2];
- }
-
-
- /**
- * calculate file index based on coordinates
- */
- auto hkl_to_idx = [this](t_real h, t_real k, t_real l) -> std::size_t
- {
- // clamp values to boundaries
- if(h < m_hmin) h = m_hmin;
- if(k < m_kmin) k = m_kmin;
- if(l < m_lmin) l = m_lmin;
- if(h >= m_hmax) h = m_hmax - m_hstep;
- if(k >= m_kmax) k = m_kmax - m_kstep;
- if(l >= m_lmax) l = m_lmax - m_lstep;
-
- // max dimensions
- std::size_t iHSize = std::size_t(((m_hmax-m_hmin) / m_hstep));
- std::size_t iKSize = std::size_t(((m_kmax-m_kmin) / m_kstep));
- std::size_t iLSize = std::size_t(((m_lmax-m_lmin) / m_lstep));
-
- // position indices
- std::size_t iH = std::size_t(std::round(((h - m_hmin) / m_hstep)));
- std::size_t iK = std::size_t(std::round(((k - m_kmin) / m_kstep)));
- std::size_t iL = std::size_t(std::round(((l - m_lmin) / m_lstep)));
-
- // clamp again
- if(iH >= iHSize) iH = iHSize-1;
- if(iK >= iKSize) iK = iKSize-1;
- if(iL >= iLSize) iL = iLSize-1;
-
- return iH*iKSize*iLSize + iK*iLSize + iL;
- };
-
-
- // ------------------------------------------------------------------------
- // the index file holds the offsets into the data file
- t_real dqh = dh-m_vecG[0], dqk = dk-m_vecG[1], dql = dl-m_vecG[2];
- if(m_bFlipB_or_q)
- {
- dqh = -dqh;
- dqk = -dqk;
- dql = -dql;
- }
- std::size_t idx_file_offs = hkl_to_idx(dqh, dqk, dql);
-
- QFile fileIdx(m_strIndexFile.c_str());
- if(!fileIdx.exists())
- {
- tl2::log_err("Index file \"", m_strIndexFile, "\" does not exist.");
- return std::make_tuple(std::vector<t_real>(), std::vector<t_real>());
- }
-
- if(!fileIdx.open(QIODevice::ReadOnly))
- {
- tl2::log_err("Index file \"", m_strIndexFile, "\" cannot be opened.");
- return std::make_tuple(std::vector<t_real>(), std::vector<t_real>());
- }
-
- const void *pMemIdx = fileIdx.map(idx_file_offs*sizeof(std::size_t), sizeof(std::size_t));
- if(!pMemIdx)
- {
- tl2::log_err("Index file \"", m_strIndexFile, "\" cannot be mapped.");
- return std::make_tuple(std::vector<t_real>(), std::vector<t_real>());
- }
-
- std::size_t dat_file_offs = *((std::size_t*)pMemIdx);
- fileIdx.unmap((unsigned char*)pMemIdx);
- fileIdx.close();
- // ------------------------------------------------------------------------
-
-
- // ------------------------------------------------------------------------
- // the data file holds the energies and spectral weights of the dispersion branches
-
- QFile fileDat(m_strDataFile.c_str());
- if(!fileDat.exists())
- {
- tl2::log_err("Data file \"", m_strDataFile, "\" does not exist.");
- return std::make_tuple(std::vector<t_real>(), std::vector<t_real>());
- }
-
- if(!fileDat.open(QIODevice::ReadOnly))
- {
- tl2::log_err("Data file \"", m_strDataFile, "\" cannot be opened.");
- return std::make_tuple(std::vector<t_real>(), std::vector<t_real>());
- }
-
- const void *pMemDat = fileDat.map(dat_file_offs, sizeof(std::size_t));
- if(!pMemDat)
- {
- tl2::log_err("Data file \"", m_strDataFile, "\" cannot be mapped (1).");
- return std::make_tuple(std::vector<t_real>(), std::vector<t_real>());
- }
-
- // number of dispersion branches and weights
- unsigned int iNumBranches = *((unsigned int*)pMemDat);
- fileDat.unmap((unsigned char*)pMemDat);
-
-
- // map actual (E, w) data
- const t_real *pBranches = (t_real*)fileDat.map(dat_file_offs+sizeof(iNumBranches), iNumBranches*sizeof(t_real)*4);
- if(!pMemDat)
- {
- tl2::log_err("Data file \"", m_strDataFile, "\" cannot be mapped (2).");
- return std::make_tuple(std::vector<t_real>(), std::vector<t_real>());
- }
-
-
- std::vector<t_real> vecE, vecw, vecwSF1, vecwSF2, vecwNSF;
- for(unsigned int iBranch=0; iBranch<iNumBranches; ++iBranch)
- {
- vecE.push_back(pBranches[iBranch*4 + 0] * m_dEFactor / m_dEFactorGrid); // energy
- vecwSF1.push_back(pBranches[iBranch*4 + 1]); // weight SF1
- vecwSF2.push_back(pBranches[iBranch*4 + 2]); // weight SF2
- vecwNSF.push_back(pBranches[iBranch*4 + 3]); // weight NSF
- vecw.push_back(pBranches[iBranch*4 + 1] + pBranches[iBranch*4 + 2] + pBranches[iBranch*4 + 3]); // full weight
- }
-
- fileDat.unmap((unsigned char*)pBranches);
- fileDat.close();
- // ------------------------------------------------------------------------
-
-
- if(m_iPolChan == 1) return std::make_tuple(vecE, vecwSF1);
- else if(m_iPolChan == 2) return std::make_tuple(vecE, vecwSF2);
- else if(m_iPolChan == 3) return std::make_tuple(vecE, vecwNSF);
- return std::make_tuple(vecE, vecw);
-}
-
-
-/**
- * S(q,E)
- */
-t_real SqwMod::operator()(t_real dh, t_real dk, t_real dl, t_real dE) const
-{
- t_real dInc=0, dS_p=0, dS_m=0;
- if(!tl2::float_equal(m_dIncAmp, t_real(0)))
- dInc = tl2::gauss_model(dE, t_real(0), m_dIncSigma, m_dIncAmp, t_real(0));
-
- t_real dS = 0;
- if(!tl2::float_equal(m_dS0, t_real(0)))
- {
- std::vector<t_real> vecE, vecW;
- std::tie(vecE, vecW) = disp(dh, dk, dl);
-
- for(std::size_t iE=0; iE<vecE.size(); ++iE)
- dS += tl2::gauss_model(dE, vecE[iE], m_dSigma, vecW[iE], t_real(0));
- }
-
- return m_dS0*dS * tl2::bose_cutoff(dE, m_dT, m_dcut) + dInc;
-}
-
-
-
-// ----------------------------------------------------------------------------
-// get & set variables
-
-std::vector<SqwMod::t_var> SqwMod::GetVars() const
-{
- std::vector<t_var> vecVars;
-
- vecVars.push_back(SqwBase::t_var{"bose_cutoff", "real", tl2::var_to_str(m_dcut)});
- vecVars.push_back(SqwBase::t_var{"sigma", "real", tl2::var_to_str(m_dSigma)});
- vecVars.push_back(SqwBase::t_var{"inc_amp", "real", tl2::var_to_str(m_dIncAmp)});
- vecVars.push_back(SqwBase::t_var{"inc_sigma", "real", tl2::var_to_str(m_dIncSigma)});
- vecVars.push_back(SqwBase::t_var{"S0", "real", tl2::var_to_str(m_dS0)});
- vecVars.push_back(SqwBase::t_var{"T", "real", tl2::var_to_str(m_dT)});
- vecVars.push_back(SqwBase::t_var{"rot_z", "real", tl2::var_to_str(m_dRotZ)});
- vecVars.push_back(SqwBase::t_var{"rot_x", "real", tl2::var_to_str(m_dRotX)});
- vecVars.push_back(SqwBase::t_var{"pol_chan", "int", tl2::var_to_str(m_iPolChan)});
- vecVars.push_back(SqwBase::t_var{"G", "vector", vec_to_str(m_vecG)});
- vecVars.push_back(SqwBase::t_var{"rot_from", "vector", vec_to_str(m_vecRotFrom)});
- vecVars.push_back(SqwBase::t_var{"rot_to", "vector", vec_to_str(m_vecRotTo)});
- vecVars.push_back(SqwBase::t_var{"skx_E_factor", "real", tl2::var_to_str(m_dEFactor)});
- vecVars.push_back(SqwBase::t_var{"flip_B_or_q", "int", tl2::var_to_str(int(m_bFlipB_or_q))});
-
- return vecVars;
-}
-
-
-void SqwMod::SetVars(const std::vector<SqwMod::t_var>& vecVars)
-{
- if(!vecVars.size()) return;
-
- for(const SqwBase::t_var& var : vecVars)
- {
- const std::string& strVar = std::get<0>(var);
- const std::string& strVal = std::get<2>(var);
-
- if(strVar == "bose_cutoff") m_dcut = tl2::str_to_var<t_real>(strVal);
- else if(strVar == "sigma") m_dSigma = tl2::str_to_var<t_real>(strVal);
- else if(strVar == "inc_amp") m_dIncAmp = tl2::str_to_var<decltype(m_dIncAmp)>(strVal);
- else if(strVar == "inc_sigma") m_dIncSigma = tl2::str_to_var<decltype(m_dIncSigma)>(strVal);
- else if(strVar == "S0") m_dS0 = tl2::str_to_var<decltype(m_dS0)>(strVal);
- else if(strVar == "T") m_dT = tl2::str_to_var<decltype(m_dT)>(strVal);
- else if(strVar == "rot_z") m_dRotZ = tl2::str_to_var<decltype(m_dRotZ)>(strVal);
- else if(strVar == "rot_x") m_dRotZ = tl2::str_to_var<decltype(m_dRotX)>(strVal);
- else if(strVar == "pol_chan") m_iPolChan = tl2::str_to_var<decltype(m_iPolChan)>(strVal);
- else if(strVar == "G") m_vecG = str_to_vec<decltype(m_vecG)>(strVal);
- else if(strVar == "rot_from") m_vecRotFrom = str_to_vec<decltype(m_vecRotFrom)>(strVal);
- else if(strVar == "rot_to") m_vecRotTo = str_to_vec<decltype(m_vecRotTo)>(strVal);
- else if(strVar == "skx_E_factor") m_dEFactor = tl2::str_to_var<decltype(m_dT)>(strVal);
- else if(strVar == "flip_B_or_q") m_bFlipB_or_q = (tl2::str_to_var<int>(strVal) != 0);
- }
-}
-
-
-bool SqwMod::SetVarIfAvail(const std::string& strKey, const std::string& strNewVal)
-{
- return SqwBase::SetVarIfAvail(strKey, strNewVal);
-}
-
-
-
-// ----------------------------------------------------------------------------
-// copy
-
-SqwBase* SqwMod::shallow_copy() const
-{
- SqwMod *pMod = new SqwMod();
-
- pMod->m_dSigma = this->m_dSigma;
- pMod->m_dIncAmp = this->m_dIncAmp;
- pMod->m_dIncSigma = this->m_dIncSigma;
- pMod->m_dS0 = this->m_dS0;
- pMod->m_dT = this->m_dT;
- pMod->m_dRotZ = this->m_dRotZ;
- pMod->m_dRotX = this->m_dRotX;
- pMod->m_iPolChan = this->m_iPolChan;
- pMod->m_vecG = this->m_vecG;
- pMod->m_vecRotFrom = this->m_vecRotFrom;
- pMod->m_vecRotTo = this->m_vecRotTo;
- pMod->m_dcut = this->m_dcut;
- pMod->m_dEFactor = this->m_dEFactor;
- pMod->m_dEFactorGrid = this->m_dEFactorGrid;
- pMod->m_bFlipB_or_q = this->m_bFlipB_or_q;
-
- pMod->m_strIndexFile = this->m_strIndexFile;
- pMod->m_strDataFile = this->m_strDataFile;
-
- pMod->m_hmin = this->m_hmin;
- pMod->m_hmax = this->m_hmax;
- pMod->m_hstep = this->m_hstep;
- pMod->m_kmin = this->m_kmin;
- pMod->m_kmax = this->m_kmax;
- pMod->m_kstep = this->m_kstep;
- pMod->m_lmin = this->m_lmin;
- pMod->m_lmax = this->m_lmax;
- pMod->m_lstep = this->m_lstep;
-
- return pMod;
-}
-
-
-
-#ifdef DO_TEST
-// ----------------------------------------------------------------------------
-// test querying of data
-
-int main(int argc, char **argv)
-{
- if(argc <= 1)
- {
- std::cerr << "Please specify a config file." << std::endl;
- return -1;
- }
-
-
- SqwMod mod(argv[1]);
- mod.SetVarIfAvail("pol_chan", "0");
-
- while(1)
- {
- t_real h, k, l;
- std::cout << "Enter hkl: ";
- std::cin >> h >> k >> l;
-
- std::vector<t_real> vecE, vecW;
- std::tie(vecE, vecW) = mod.disp(h, k, l);
-
- for(std::size_t i=0; i<vecE.size(); ++i)
- std::cout << "(" << i+1 << ") " << "E = " << vecE[i] << ", weight = " << vecW[i] << "\n";
- std::cout << std::endl;
- }
-}
-
-// ----------------------------------------------------------------------------
-#else
-// ----------------------------------------------------------------------------
-// SO interface
-
-static const char* pcModIdent = "skxmod_grid";
-static const char* pcModName = "MnSi Magnon Dynamics (Grid)";
-
-
-std::tuple<std::string, std::string, std::string> sqw_info()
-{
- return std::make_tuple(TAKIN_VER, pcModIdent, pcModName);
-}
-
-
-std::shared_ptr<SqwBase> sqw_construct(const std::string& strCfgFile)
-{
- return std::make_shared<SqwMod>(strCfgFile);
-}
-
-
-// exports from so file
-#ifndef __MINGW32__
- BOOST_DLL_ALIAS(sqw_info, takin_sqw_info);
- BOOST_DLL_ALIAS(sqw_construct, takin_sqw);
-#else
- // hack because BOOST_DLL_ALIAS does not seem to work with Mingw
-
- extern "C" __declspec(dllexport)
- std::tuple<std::string, std::string, std::string> takin_sqw_info()
- {
- return sqw_info();
- }
-
- extern "C" __declspec(dllexport)
- std::shared_ptr<SqwBase> takin_sqw(const std::string& strCfgFile)
- {
- return sqw_construct(strCfgFile);
- }
-#endif
-// ----------------------------------------------------------------------------
-
-
-#endif
--
GitLab