Compare commits

..

37 Commits

Author SHA1 Message Date
Marek Nečada 165dbe3d61 Forgotten source catch_aux.C 2022-05-29 14:37:50 +03:00
Marek Nečada b756453d12 Forgotten header lattice_types.h 2022-05-29 14:36:32 +03:00
Marek Nečada 16a0a7dc0a Dummy (zero) T-matrix generator 2022-04-04 09:43:02 +03:00
Marek Nečada 3052895ac6 VSWF fill test, fix memory issue 2022-03-30 14:03:49 +03:00
Marek Nečada 9e788b9d2f Missing include 2022-03-30 13:07:27 +03:00
Marek Nečada 923e1e0a77 Indexing.h C++ compatibility; flag C++ incompatible headers. 2022-03-30 11:15:18 +03:00
Marek Nečada 856903a98b Restructure headers to ensure ewald.h C++ compatibility 2022-03-30 10:58:07 +03:00
Marek Nečada 89f7400c71 vectors.h explicit type conversions for C++ compatibility 2022-03-30 09:59:04 +03:00
Marek Nečada 6820eb3366 extern "C" in headers 2022-03-30 09:47:06 +03:00
Marek Nečada e66f95aa4f Introduce testing witch catch2 (no runtime tests yet) 2022-03-30 09:31:25 +03:00
Marek Nečada b04a4deebe Replace "complex" with "_Complex" in headers
Needed for compatibility with C++ (in which "complex" is
a template).
2022-03-30 09:06:47 +03:00
Marek Nečada 736adb6974 Expressing VSWF in terms of SSWF 2022-03-29 14:57:02 +03:00
Marek Nečada 10d3beca1c Indexing documentation, pedantic use of qpms_y_sc_t 2022-03-25 21:53:45 +02:00
Marek Nečada f9081c5922 Remove duplicit function declaration.
qpms_scatsysw_scattered_field_basis()
2022-03-25 16:10:24 +02:00
Marek Nečada 54936ee722 analysis: fix symmetry group info 2022-01-03 23:51:13 +02:00
Marek Nečada 63f7f216b0 Analysis script also for finiterectlat constant driving results 2022-01-03 20:58:01 +02:00
Marek Nečada 0315481977 Expose more testing functions via cython api. 2021-12-21 12:42:23 +02:00
Marek Nečada dca04153f2 Fix unitialised diagonal blocks in full modeproblem matrix.
The bug in qpms_scatsysw_build_modeproblem_matrix_full()
via qpms_scatsysw_scatsyswk_build_modeproblem_matrix_full()
affected finite-size arrays: instead of filling the diagonal
interaction blocks with zeros, they were filled with undefined
values from a superfluous zgemm() call.
2021-12-21 11:29:44 +02:00
Marek Nečada d40f1fd004 analysis.py: fix duplicities in position mapping due to rounding errors 2021-12-09 14:53:57 +02:00
Marek Nečada dd0ee9e2d4 Helper class to analyse results of finiterectlat-modes.py 2021-12-09 14:22:46 +02:00
Marek Nečada 0a7e591aa3 WIP fixing some example scripts 2021-12-09 11:18:14 +02:00
Marek Nečada d600977ca5 Move old cython translation_calculator wrapper 2021-12-09 11:18:14 +02:00
Marek Nečada f283bc01d1 Attempt to fix the Ewald sum SSWF phase bug 2021-12-09 11:18:14 +02:00
Marek Nečada b088464e3e Revert scatsystem cython parallelisation.
It didn't work correctly on my laptop.
This must be reviewed later, as the waiting is annoying.
2021-12-08 14:18:32 +02:00
Marek Nečada 227aaf4021 Dirty fix of the otherwise non-literal arrays. 2021-12-08 14:09:15 +02:00
Marek Nečada 0078f4db56 Cython evaluation of "normal" and Ewald-summed SSWF 2021-12-08 14:09:15 +02:00
Marek Nečada 210ddaa527 ScatteringSystem, ... definitions to pxd file 2021-12-08 14:09:15 +02:00
Marek Nečada c4d0700aca SSWF evaluation functions for debugging Ewald sum. 2021-12-08 14:09:15 +02:00
Marek Nečada 86d5245ea1 Additional phase factors in Ewald sum for debugging 2021-12-08 14:09:15 +02:00
Marek Nečada afda889a39 WIP documentation on phase/normalisation implementation. 2021-12-08 14:09:15 +02:00
Marek Nečada 8069634930 Ewald summation consistency test for lattice basis fields.
In-plane test currently passes, out-of-plane does not.
2021-12-08 14:09:15 +02:00
Marek Nečada a7c70344c5 finiterect-constant-driving more metadata output 2021-12-08 14:09:15 +02:00
Marek Nečada c0ade8cf2b Remove bugs of removed code from BUGS.rst.
Some of the other bugs might have been already fixed
as well, this has to be reviewed, and also doxygenised.
2021-12-08 14:09:15 +02:00
Marek Nečada fa6d93d3e8 Enable separation of long- and short-range of Ewald sum in python 2021-12-08 14:09:15 +02:00
Marek Nečada 7e390180cd Fix prange race conditions. 2021-12-08 14:09:15 +02:00
Marek Nečada 1158e116d2 "Basis fields" for finite systems.
There seem to be race conditions in the prange'd cython parts.
2021-12-08 14:09:15 +02:00
Marek Nečada 414bee6c9f "Basis fields" for finite systems 2021-12-08 14:09:15 +02:00
99 changed files with 2065 additions and 1343 deletions

View File

@ -1,83 +0,0 @@
---
kind: pipeline
type: docker
name: buildqpms-alpine-preinstlibs
workspace:
path: /home/qpmsbuild/qpms
# don't run in master until the python/lapacke linking problem is resolved
trigger:
branch:
exclude:
- master
steps:
- name: chown
image: qpms/buildenv/alpine/pkgdnumlib
pull: never
commands:
- chown -R qpmsbuild.qpmsbuild .
- name: submodules
image: qpms/buildenv/alpine/pkgdnumlib
pull: never
user: qpmsbuild
commands:
- git submodule init
- git submodule update
- name: build
image: qpms/buildenv/alpine/pkgdnumlib
pull: never
user: qpmsbuild
commands:
- cmake -DCMAKE_INSTALL_PREFIX=$HOME/.local .
- make install
- export LIBRARY_PATH=$HOME/.local/lib
- python3 setup.py install --user
- cd examples/rectangular/modes
- pip3 install --user matplotlib #needed to run the example
- export LD_LIBRARY_PATH=$HOME/.local/lib
- ./01a_realfreq_svd.sh
---
kind: pipeline
type: docker
name: buildqpms-debian-preinstlibs
workspace:
path: /home/qpmsbuild/qpms
steps:
- name: chown
image: qpms/buildenv/debian/pkgdnumlib
pull: never
commands:
- chown -R qpmsbuild.qpmsbuild .
- name: submodules
image: qpms/buildenv/debian/pkgdnumlib
pull: never
user: qpmsbuild
commands:
- git submodule init
- git submodule update
- name: build
image: qpms/buildenv/debian/pkgdnumlib
pull: never
user: qpmsbuild
commands:
- cmake -DCMAKE_INSTALL_PREFIX=/home/qpmsbuild/.local .
- make install
- export LIBRARY_PATH=$HOME/.local/lib
- python3 setup.py install --user
- pip3 install --user matplotlib #needed to run the examples
- export LD_LIBRARY_PATH=$HOME/.local/lib
- cd examples/rectangular/modes
- ./01a_realfreq_svd.sh
- cd -
- cd examples/hexagonal/modes
#- ./01a_realfreq_svd.sh
#- ./01_compute_modes.sh
#- ./02b_compute_disp_0M.sh
#- ./02_compute_disp.sh
#- ./02x_compute_disp.sh

2
.gitmodules vendored
View File

@ -1,4 +1,4 @@
[submodule "camos"]
path = camos
url = https://codeberg.org/QPMS/zbessel.git
url = https://github.com/texnokrates/zbessel.git
branch = purec

View File

@ -23,28 +23,6 @@ Plane wave decompositions gives wrong value on the longitudinal part.
The implementation of the L coefficients OR the longitudinal waves
is thus probably wrong.
Scattering result asymmetry
---------------------------
The lattice scattering code (such as finitesqlatzsym-scattery.py) produces
asymmetric results where one should not get them, due to the system symmetry.
It seems that the asymmetry appears mostly in the y-direction (i.e.
for example the scattering/absorption cross section at k = (kx, ky, kz)
is not exactly the same as k = (kx, -ky, kz).
What has been checked (hopefully):
- The flip operators for electric waves
- Some weird kind of rounding or other numerical error depending on
the position ordering of the matrix (randomized indices give
the same asymmetry).
What has not been checked (so well):
- The x,y,z-flip operators for magnetic waves (i.e. are they really
supposet to bejust the
same as the operators for electric waves, just with opposite sign?)
- zplane_pq_y
- the translation operators
Singular value asymmetry
------------------------
@ -60,6 +38,3 @@ Moreover, the non-normalized legendre functions that are used in translations.c
are likely to overflow for even smaller values of l.
Therefore TODO: modify translations.c in a way that uses normalised legendre functions everywhere.
Memory management
-----------------
When compiled with optimizations, I get stack smashing error in qpms_trans_calculator_get_2DFT_longrange_AB_arrays().

View File

@ -10,7 +10,7 @@ include(GNUInstallDirs)
project (QPMS)
list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake/")
macro(use_c99)
if (CMAKE_VERSION VERSION_LESS "3.1")
@ -48,6 +48,7 @@ add_subdirectory(faddeeva)
add_subdirectory (qpms)
add_subdirectory (tests/catch EXCLUDE_FROM_ALL)
enable_testing()
add_subdirectory (tests EXCLUDE_FROM_ALL)

View File

@ -753,7 +753,7 @@ WARN_LOGFILE =
# spaces.
# Note: If this tag is empty the current directory is searched.
INPUT = qpms notes misc finite_systems.md MIRRORS.md CLIUTILS.md README.md README.Triton.md finite_systems.md lattices.md TODO.md
INPUT = qpms notes misc finite_systems.md CLIUTILS.md README.md README.Triton.md finite_systems.md lattices.md TODO.md
# This tag can be used to specify the character encoding of the source files
# that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses

View File

@ -1,14 +0,0 @@
QPMS source code mirrors
========================
QPMS source code is available at several locations; in all of the following,
upstream `master` branch is kept up-to-date. Various development branches
are not necessarily pushed everywhere (and they should be considered
unstable in the sense that rebases and forced pushes are possible).
mirror | note | provider | backend
----------------------------------------------- | ----------------------- | ------------------------------------------------- | ------
<https://repo.or.cz/qpms.git> | primary public upstream | [repo.or.cz](https://repo.or.cz/) | girocco
<https://codeberg.org/QPMS/qpms> | | [Codeberg](https://codeberg.org) | gitea
<https://git.piraattipuolue.fi/QPMS/qpms.git> | | [Pirate Party Finland](https://piraattipuolue.fi) | gitea
<https://version.aalto.fi/gitlab/qpms/qpms.git> | | [Aalto University](https://aalto.fi) | gitlab

View File

@ -1,9 +1,7 @@
[![Build Status](https://drone.perkele.eu/api/badges/QPMS/qpms/status.svg)](https://drone.perkele.eu/QPMS/qpms)
QPMS README
===========
[QPMS][homepage] (standing for QPMS Photonic Multiple Scattering)
QPMS (standing for QPMS Photonic Multiple Scattering)
is a toolkit for frequency-domain simulations of photonic systems
consisting of compact objects (particles) inside a homogeneous medium. Scattering
properties of the individual particles are described by their T-matrices
@ -45,8 +43,7 @@ Infinite systems (lattices)
Getting the code
================
The codebase is available at the main upstream public repository
<https://repo.or.cz/qpms.git> or any of the [maintained mirrors][MIRRORS].
The main upstream public repository is located at <https://repo.or.cz/qpms.git>.
Just clone the repository with `git` and proceed to the installation instructions
below.
@ -92,21 +89,14 @@ Special care might need to be taken when installing QPMS in cluster environments
Specific installation instructions for Aalto University's Triton cluster
can be found in a [separate document][TRITON-README].
Instructions for installation on Android-based devices are
in [another document][INSTALL-ANDROID].
Documentation
=============
[QPMS documentation][homepage] is a work in progress. Most of the newer code
is documented using [doxygen][] comments. Documentation generated for the
upstream version is hosted on the QPMS homepage <https://qpms.necada.org>.
To build the documentation yourself,
just run
Documentation of QPMS is a work in progress. Most of the newer code
is documented using [doxygen][] comments. To build the documentation, just run
`doxygen`
in the QPMS source root directory; the documentation will then be found in
in the root directory; the documentation will then be found in
`docs/html/index.html`.
Of course, the prerequisite of this is having doxygen installed.
@ -176,13 +166,12 @@ You are also warmly welcome to the [QPMS user chat][telegramchat]
in Telegram!
[homepage]: https://qpms.necada.org
[SCUFF-EM]: https://homerreid.github.io/scuff-em-documentation/
[OpenBLAS]: https://www.openblas.net/
[GSL]: https://www.gnu.org/software/gsl/
[cmake]: https://cmake.org
[TRITON-README]: README.Triton.md
[INSTALL-ANDROID]: notes/INSTALL_ANDROID.md
[tutorial-finite]: finite_systems.md
[tutorial-infinite]: lattices.md
[doxygen]: http://doxygen.nl/
@ -191,4 +180,3 @@ in Telegram!
[telegramchat]: https://t.me/QPMScattering
[authormail]: mailto:marek@necada.org
[cliutils]: CLIUTILS.md
[MIRRORS]: MIRRORS.md

View File

@ -1,7 +0,0 @@
#!/bin/sh
DP=Dockerfile_parts
# "Build environment" Dockerfiles
cat >Dockerfile.benv.debian.bnl ${DP}/00_common.debian ${DP}/01_numlibs.built
cat >Dockerfile.benv.debian.pnl ${DP}/00_common.debian ${DP}/01_numlibs.debian.pkgd
cat >Dockerfile.benv.alpine.bnl ${DP}/00_common.alpine ${DP}/01_numlibs.built
cat >Dockerfile.benv.alpine.pnl ${DP}/00_common.alpine ${DP}/01_numlibs.alpine.pkgd

View File

@ -1,5 +0,0 @@
#!/bin/sh
docker build -t qpms/buildenv/debian/builtnumlib -f Dockerfile.benv.debian.bnl .
docker build -t qpms/buildenv/debian/pkgdnumlib -f Dockerfile.benv.debian.pnl .
docker build -t qpms/buildenv/alpine/builtnumlib -f Dockerfile.benv.alpine.bnl .
docker build -t qpms/buildenv/alpine/pkgdnumlib -f Dockerfile.benv.alpine.pnl .

View File

@ -1,5 +0,0 @@
FROM alpine:latest AS commondeps
RUN apk update \
&& apk add cmake python3-dev py3-pip gcc g++ wget git make libc-dev bc \
&& adduser -D qpmsbuild

View File

@ -1,6 +0,0 @@
FROM debian:stable AS commondeps
RUN apt-get update \
&& apt-get -y install --no-install-recommends build-essential cmake python3 python3-pip git wget python3-dev bc \
&& apt-get clean \
&& useradd -m qpmsbuild

View File

@ -1,3 +0,0 @@
FROM commondeps AS numlibs
# openblas-dev adds gfortran :(
RUN apk add openblas-dev gsl-dev

View File

@ -1,16 +0,0 @@
FROM commondeps AS buildopenblas
USER qpmsbuild
RUN cd && git clone --depth 1 https://github.com/xianyi/OpenBLAS.git \
&& cd OpenBLAS && make \
&& make install PREFIX=$HOME/.local/ \
&& make clean \
&& cd .. && rm -rf OpenBLAS
FROM buildopenblas AS numlibs
USER qpmsbuild
RUN cd && wget https://ftp.gnu.org/gnu/gsl/gsl-latest.tar.gz \
&& tar xf gsl-latest.tar.gz \
&& cd $( tar tf gsl-latest.tar.gz | head -n 1 ) \
&& ./configure --prefix=$HOME/.local \
&& make && make install && make clean \
&& cd .. && rm -rf $OLDPWD gsl-latest.tar.gz

View File

@ -1,4 +0,0 @@
FROM commondeps AS numlibs
RUN apt-get -y install --no-install-recommends libopenblas-dev libgsl-dev liblapacke-dev \
&& apt-get clean

View File

@ -1,12 +0,0 @@
FROM numlibs AS buildqpms
USER qpmsbuild
ENV LD_LIBRARY_PATH /home/qpmsbuild/.local/lib
ENV LIBRARY_PATH /home/qpmsbuild/.local/lib
ENV C_INCLUDE_PATH /home/qpmsbuild/.local/include
RUN cd && git clone --depth 1 https://repo.or.cz/qpms.git \
&& cd qpms && git submodule init && git submodule update
RUN cd ~/qpms && cmake -DCMAKE_INSTALL_PREFIX=$HOME/.local . \
&& make \
&& make install
RUN cd ~/qpms && python3 setup.py install --user

View File

@ -1 +0,0 @@
../.drone.yml

View File

@ -1,284 +0,0 @@
# - Returns a version string from Git
#
# These functions force a re-configure on each git commit so that you can
# trust the values of the variables in your build system.
#
# get_git_head_revision(<refspecvar> <hashvar> [ALLOW_LOOKING_ABOVE_CMAKE_SOURCE_DIR])
#
# Returns the refspec and sha hash of the current head revision
#
# git_describe(<var> [<additional arguments to git describe> ...])
#
# Returns the results of git describe on the source tree, and adjusting
# the output so that it tests false if an error occurs.
#
# git_describe_working_tree(<var> [<additional arguments to git describe> ...])
#
# Returns the results of git describe on the working tree (--dirty option),
# and adjusting the output so that it tests false if an error occurs.
#
# git_get_exact_tag(<var> [<additional arguments to git describe> ...])
#
# Returns the results of git describe --exact-match on the source tree,
# and adjusting the output so that it tests false if there was no exact
# matching tag.
#
# git_local_changes(<var>)
#
# Returns either "CLEAN" or "DIRTY" with respect to uncommitted changes.
# Uses the return code of "git diff-index --quiet HEAD --".
# Does not regard untracked files.
#
# Requires CMake 2.6 or newer (uses the 'function' command)
#
# Original Author:
# 2009-2020 Ryan Pavlik <ryan.pavlik@gmail.com> <abiryan@ryand.net>
# http://academic.cleardefinition.com
#
# Copyright 2009-2013, Iowa State University.
# Copyright 2013-2020, Ryan Pavlik
# Copyright 2013-2020, Contributors
# SPDX-License-Identifier: BSL-1.0
# Distributed under the Boost Software License, Version 1.0.
# (See accompanying file LICENSE_1_0.txt or copy at
# http://www.boost.org/LICENSE_1_0.txt)
if(__get_git_revision_description)
return()
endif()
set(__get_git_revision_description YES)
# We must run the following at "include" time, not at function call time,
# to find the path to this module rather than the path to a calling list file
get_filename_component(_gitdescmoddir ${CMAKE_CURRENT_LIST_FILE} PATH)
# Function _git_find_closest_git_dir finds the next closest .git directory
# that is part of any directory in the path defined by _start_dir.
# The result is returned in the parent scope variable whose name is passed
# as variable _git_dir_var. If no .git directory can be found, the
# function returns an empty string via _git_dir_var.
#
# Example: Given a path C:/bla/foo/bar and assuming C:/bla/.git exists and
# neither foo nor bar contain a file/directory .git. This wil return
# C:/bla/.git
#
function(_git_find_closest_git_dir _start_dir _git_dir_var)
set(cur_dir "${_start_dir}")
set(git_dir "${_start_dir}/.git")
while(NOT EXISTS "${git_dir}")
# .git dir not found, search parent directories
set(git_previous_parent "${cur_dir}")
get_filename_component(cur_dir "${cur_dir}" DIRECTORY)
if(cur_dir STREQUAL git_previous_parent)
# We have reached the root directory, we are not in git
set(${_git_dir_var}
""
PARENT_SCOPE)
return()
endif()
set(git_dir "${cur_dir}/.git")
endwhile()
set(${_git_dir_var}
"${git_dir}"
PARENT_SCOPE)
endfunction()
function(get_git_head_revision _refspecvar _hashvar)
_git_find_closest_git_dir("${CMAKE_CURRENT_SOURCE_DIR}" GIT_DIR)
if("${ARGN}" STREQUAL "ALLOW_LOOKING_ABOVE_CMAKE_SOURCE_DIR")
set(ALLOW_LOOKING_ABOVE_CMAKE_SOURCE_DIR TRUE)
else()
set(ALLOW_LOOKING_ABOVE_CMAKE_SOURCE_DIR FALSE)
endif()
if(NOT "${GIT_DIR}" STREQUAL "")
file(RELATIVE_PATH _relative_to_source_dir "${CMAKE_SOURCE_DIR}"
"${GIT_DIR}")
if("${_relative_to_source_dir}" MATCHES "[.][.]" AND NOT ALLOW_LOOKING_ABOVE_CMAKE_SOURCE_DIR)
# We've gone above the CMake root dir.
set(GIT_DIR "")
endif()
endif()
if("${GIT_DIR}" STREQUAL "")
set(${_refspecvar}
"GITDIR-NOTFOUND"
PARENT_SCOPE)
set(${_hashvar}
"GITDIR-NOTFOUND"
PARENT_SCOPE)
return()
endif()
# Check if the current source dir is a git submodule or a worktree.
# In both cases .git is a file instead of a directory.
#
if(NOT IS_DIRECTORY ${GIT_DIR})
# The following git command will return a non empty string that
# points to the super project working tree if the current
# source dir is inside a git submodule.
# Otherwise the command will return an empty string.
#
execute_process(
COMMAND "${GIT_EXECUTABLE}" rev-parse
--show-superproject-working-tree
WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}"
OUTPUT_VARIABLE out
ERROR_QUIET OUTPUT_STRIP_TRAILING_WHITESPACE)
if(NOT "${out}" STREQUAL "")
# If out is empty, GIT_DIR/CMAKE_CURRENT_SOURCE_DIR is in a submodule
file(READ ${GIT_DIR} submodule)
string(REGEX REPLACE "gitdir: (.*)$" "\\1" GIT_DIR_RELATIVE
${submodule})
string(STRIP ${GIT_DIR_RELATIVE} GIT_DIR_RELATIVE)
get_filename_component(SUBMODULE_DIR ${GIT_DIR} PATH)
get_filename_component(GIT_DIR ${SUBMODULE_DIR}/${GIT_DIR_RELATIVE}
ABSOLUTE)
set(HEAD_SOURCE_FILE "${GIT_DIR}/HEAD")
else()
# GIT_DIR/CMAKE_CURRENT_SOURCE_DIR is in a worktree
file(READ ${GIT_DIR} worktree_ref)
# The .git directory contains a path to the worktree information directory
# inside the parent git repo of the worktree.
#
string(REGEX REPLACE "gitdir: (.*)$" "\\1" git_worktree_dir
${worktree_ref})
string(STRIP ${git_worktree_dir} git_worktree_dir)
_git_find_closest_git_dir("${git_worktree_dir}" GIT_DIR)
set(HEAD_SOURCE_FILE "${git_worktree_dir}/HEAD")
endif()
else()
set(HEAD_SOURCE_FILE "${GIT_DIR}/HEAD")
endif()
set(GIT_DATA "${CMAKE_CURRENT_BINARY_DIR}/CMakeFiles/git-data")
if(NOT EXISTS "${GIT_DATA}")
file(MAKE_DIRECTORY "${GIT_DATA}")
endif()
if(NOT EXISTS "${HEAD_SOURCE_FILE}")
return()
endif()
set(HEAD_FILE "${GIT_DATA}/HEAD")
configure_file("${HEAD_SOURCE_FILE}" "${HEAD_FILE}" COPYONLY)
configure_file("${_gitdescmoddir}/GetGitRevisionDescription.cmake.in"
"${GIT_DATA}/grabRef.cmake" @ONLY)
include("${GIT_DATA}/grabRef.cmake")
set(${_refspecvar}
"${HEAD_REF}"
PARENT_SCOPE)
set(${_hashvar}
"${HEAD_HASH}"
PARENT_SCOPE)
endfunction()
function(git_describe _var)
if(NOT GIT_FOUND)
find_package(Git QUIET)
endif()
get_git_head_revision(refspec hash)
if(NOT GIT_FOUND)
set(${_var}
"GIT-NOTFOUND"
PARENT_SCOPE)
return()
endif()
if(NOT hash)
set(${_var}
"HEAD-HASH-NOTFOUND"
PARENT_SCOPE)
return()
endif()
# TODO sanitize
#if((${ARGN}" MATCHES "&&") OR
# (ARGN MATCHES "||") OR
# (ARGN MATCHES "\\;"))
# message("Please report the following error to the project!")
# message(FATAL_ERROR "Looks like someone's doing something nefarious with git_describe! Passed arguments ${ARGN}")
#endif()
#message(STATUS "Arguments to execute_process: ${ARGN}")
execute_process(
COMMAND "${GIT_EXECUTABLE}" describe --tags --always ${hash} ${ARGN}
WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}"
RESULT_VARIABLE res
OUTPUT_VARIABLE out
ERROR_QUIET OUTPUT_STRIP_TRAILING_WHITESPACE)
if(NOT res EQUAL 0)
set(out "${out}-${res}-NOTFOUND")
endif()
set(${_var}
"${out}"
PARENT_SCOPE)
endfunction()
function(git_describe_working_tree _var)
if(NOT GIT_FOUND)
find_package(Git QUIET)
endif()
if(NOT GIT_FOUND)
set(${_var}
"GIT-NOTFOUND"
PARENT_SCOPE)
return()
endif()
execute_process(
COMMAND "${GIT_EXECUTABLE}" describe --dirty ${ARGN}
WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}"
RESULT_VARIABLE res
OUTPUT_VARIABLE out
ERROR_QUIET OUTPUT_STRIP_TRAILING_WHITESPACE)
if(NOT res EQUAL 0)
set(out "${out}-${res}-NOTFOUND")
endif()
set(${_var}
"${out}"
PARENT_SCOPE)
endfunction()
function(git_get_exact_tag _var)
git_describe(out --exact-match ${ARGN})
set(${_var}
"${out}"
PARENT_SCOPE)
endfunction()
function(git_local_changes _var)
if(NOT GIT_FOUND)
find_package(Git QUIET)
endif()
get_git_head_revision(refspec hash)
if(NOT GIT_FOUND)
set(${_var}
"GIT-NOTFOUND"
PARENT_SCOPE)
return()
endif()
if(NOT hash)
set(${_var}
"HEAD-HASH-NOTFOUND"
PARENT_SCOPE)
return()
endif()
execute_process(
COMMAND "${GIT_EXECUTABLE}" diff-index --quiet HEAD --
WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}"
RESULT_VARIABLE res
OUTPUT_VARIABLE out
ERROR_QUIET OUTPUT_STRIP_TRAILING_WHITESPACE)
if(res EQUAL 0)
set(${_var}
"CLEAN"
PARENT_SCOPE)
else()
set(${_var}
"DIRTY"
PARENT_SCOPE)
endif()
endfunction()

View File

@ -1,43 +0,0 @@
#
# Internal file for GetGitRevisionDescription.cmake
#
# Requires CMake 2.6 or newer (uses the 'function' command)
#
# Original Author:
# 2009-2010 Ryan Pavlik <rpavlik@iastate.edu> <abiryan@ryand.net>
# http://academic.cleardefinition.com
# Iowa State University HCI Graduate Program/VRAC
#
# Copyright 2009-2012, Iowa State University
# Copyright 2011-2015, Contributors
# Distributed under the Boost Software License, Version 1.0.
# (See accompanying file LICENSE_1_0.txt or copy at
# http://www.boost.org/LICENSE_1_0.txt)
# SPDX-License-Identifier: BSL-1.0
set(HEAD_HASH)
file(READ "@HEAD_FILE@" HEAD_CONTENTS LIMIT 1024)
string(STRIP "${HEAD_CONTENTS}" HEAD_CONTENTS)
if(HEAD_CONTENTS MATCHES "ref")
# named branch
string(REPLACE "ref: " "" HEAD_REF "${HEAD_CONTENTS}")
if(EXISTS "@GIT_DIR@/${HEAD_REF}")
configure_file("@GIT_DIR@/${HEAD_REF}" "@GIT_DATA@/head-ref" COPYONLY)
else()
configure_file("@GIT_DIR@/packed-refs" "@GIT_DATA@/packed-refs" COPYONLY)
file(READ "@GIT_DATA@/packed-refs" PACKED_REFS)
if(${PACKED_REFS} MATCHES "([0-9a-z]*) ${HEAD_REF}")
set(HEAD_HASH "${CMAKE_MATCH_1}")
endif()
endif()
else()
# detached HEAD
configure_file("@GIT_DIR@/HEAD" "@GIT_DATA@/head-ref" COPYONLY)
endif()
if(NOT HEAD_HASH)
file(READ "@GIT_DATA@/head-ref" HEAD_HASH LIMIT 1024)
string(STRIP "${HEAD_HASH}" HEAD_HASH)
endif()

View File

@ -1,23 +0,0 @@
Boost Software License - Version 1.0 - August 17th, 2003
Permission is hereby granted, free of charge, to any person or organization
obtaining a copy of the software and accompanying documentation covered by
this license (the "Software") to use, reproduce, display, distribute,
execute, and transmit the Software, and to prepare derivative works of the
Software, and to permit third-parties to whom the Software is furnished to
do so, all subject to the following:
The copyright notices in the Software and this entire statement, including
the above license grant, this restriction and the following disclaimer,
must be included in all copies of the Software, in whole or in part, and
all derivative works of the Software, unless such copies or derivative
works are solely in the form of machine-executable object code generated by
a source language processor.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
DEALINGS IN THE SOFTWARE.

View File

@ -9,10 +9,10 @@ ${MISCDIR}/lat2d_modes.py \
-b s${A1X_nm}e-9 s${A1Y_nm}e-9 \
-b s${A2X_nm}e-9 s${A2Y_nm}e-9 \
-p s${P1X_nm}e-9 s${P1Y_nm}e-9 \
-p s${P2X_nm}e-9 s${P2Y_nm}e-9 \
-L 3 -m $METAL -r ${RADIUS_nm}e-9 -H ${HEIGHT_nm}e-9 \
-k s${KPOINTX_nmi}e9 s${KPOINTY_nmi}e9 \
-k 0 0 \
-d -3 \
-t 0.01 \
-c 250 \
-P
#-k s${B1X_nmi}e9 s${B1Y_nmi}e9 \

View File

@ -1,7 +1,7 @@
#!/usr/bin/env python3
import math
from qpms.argproc import ArgParser, make_dict_action, sslice, annotate_pdf_metadata
from qpms.argproc import ArgParser, make_dict_action, sslice
figscale=3
ap = ArgParser(['rectlattice2d_finite', 'single_particle', 'single_lMax', 'single_omega'])
@ -44,7 +44,7 @@ import numpy as np
import qpms
from qpms.cybspec import BaseSpec
from qpms.cytmatrices import CTMatrix, TMatrixGenerator
from qpms.qpms_c import Particle, qpms_library_version
from qpms.qpms_c import Particle
from qpms.cymaterials import EpsMu, EpsMuGenerator, LorentzDrudeModel, lorentz_drude
from qpms.cycommon import DebugFlags, dbgmsg_enable
from qpms import FinitePointGroup, ScatteringSystem, BesselType, eV, hbar
@ -225,7 +225,7 @@ for iri in range(ss.nirreps):
scattered_full[spi, y] += scattered_ir_unpacked[spi, y]
if a.save_gradually:
iriout = outfile_tmp + ".%d" % iri
np.savez(iriout, iri=iri, meta={**vars(a), 'qpms_version' : qpms.__version__()},
np.savez(iriout, iri=iri, meta=vars(a),
omega=omega, wavenumber=wavenumber, nelem=nelem, wavevector=np.array(a.wavevector), phases=phases,
positions = ss.positions[:,:2],
scattered_ir_packed = scattered_ir[iri],
@ -251,10 +251,10 @@ if not math.isnan(a.ccd_distance):
logging.info("Far fields done")
outfile = defaultprefix + ".npz" if a.output is None else a.output
np.savez(outfile, meta={**vars(a), 'qpms_version' : qpms.__version__()},
np.savez(outfile, meta=vars(a),
omega=omega, wavenumber=wavenumber, nelem=nelem, wavevector=np.array(a.wavevector), phases=phases,
positions = ss.positions[:,:2],
scattered_ir_packed = np.array(scattered_ir, dtype=np.object),
scattered_ir_packed = scattered_ir,
scattered_full = scattered_full,
ir_contained = ir_contained,
t=t, l=l, m=m,
@ -264,6 +264,7 @@ np.savez(outfile, meta={**vars(a), 'qpms_version' : qpms.__version__()},
#ccd_size = ccd_size if not math.isnan(a.ccd_distance) else None,
ccd_points = ccd_points if not math.isnan(a.ccd_distance) else None,
ccd_fields = ccd_fields if not math.isnan(a.ccd_distance) else None,
fullvec_poffsets = ss.fullvec_poffsets,
)
logging.info("Saved to %s" % outfile)
@ -355,7 +356,6 @@ if a.plot or (a.plot_out is not None):
axes[y,gg].yaxis.set_major_formatter(plt.NullFormatter())
fig.text(0, 0, str(slicepairs[spi]), horizontalalignment='left', verticalalignment='bottom')
pp.savefig()
annotate_pdf_metadata(pp, scriptname="finiterectlat-constant-driving.py")
pp.close()
exit(0)

View File

@ -1,7 +1,7 @@
#!/usr/bin/env python3
import math
from qpms.argproc import ArgParser, annotate_pdf_metadata
from qpms.argproc import ArgParser
ap = ArgParser(['rectlattice2d_finite', 'background_analytical', 'single_particle', 'single_lMax', ])
@ -100,7 +100,7 @@ results['inside_contour'] = inside_ellipse((results['eigval'].real, results['eig
results['refractive_index_internal'] = [medium(om).n for om in results['eigval']]
outfile = defaultprefix + (('_ir%s_%s.npz' % (str(iri), irname)) if iri is not None else '.npz') if a.output is None else a.output
np.savez(outfile, meta={**vars(a), 'qpms_version' : qpms.__version__()}, **results)
np.savez(outfile, meta=vars(a), **results)
logging.info("Saved to %s" % outfile)
exit(0)
@ -110,7 +110,7 @@ if a.plot or (a.plot_out is not None):
import matplotlib
matplotlib.use('pdf')
from matplotlib import pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(sinalpha_list, σ_ext*1e12,label='$\sigma_\mathrm{ext}$')
@ -121,9 +121,7 @@ if a.plot or (a.plot_out is not None):
ax.set_ylabel('$\sigma/\mathrm{\mu m^2}$')
plotfile = defaultprefix + ".pdf" if a.plot_out is None else a.plot_out
with PdfPages(plotfile) as pdf:
pdf.savefig(fig)
annotate_pdf_metadata(pdf, scriptname='finiterectlat-modes.py')
fig.savefig(plotfile)
exit(0)

View File

@ -1,8 +1,8 @@
#!/usr/bin/env python3
from qpms.argproc import ArgParser, annotate_pdf_metadata
import math
pi = math.pi
from qpms.argproc import ArgParser
ap = ArgParser(['rectlattice2d_finite', 'single_particle', 'single_lMax', 'omega_seq_real_ng', 'planewave'])
@ -111,7 +111,7 @@ for i, omega in enumerate(ap.allomegas):
σ_scat_arr_ir[i, j, iri] = np.vdot(fi,np.dot(translation_matrix, fi)).real/wavenumber**2
if a.save_gradually:
iriout = outfile_tmp + ".%d.%d" % (i, iri)
np.savez(iriout, omegai=i, iri=iri, meta={**vars(a), 'qpms_version' : qpms.__version__()}, omega=omega, k_sph=k_sph_list, k_cart = k_cart_arr, E_cart=E_cart_list, E_sph=np.array(E_sph),
np.savez(iriout, omegai=i, iri=iri, meta=vars(a), omega=omega, k_sph=k_sph_list, k_cart = k_cart_arr, E_cart=E_cart_list, E_sph=np.array(E_sph),
wavenumber=wavenumber, σ_ext_list_ir=σ_ext_arr_ir[i,:,iri], σ_scat_list_ir=σ_scat_list_ir[i,:,iri])
logging.info("partial results saved to %s"%iriout)
@ -122,8 +122,7 @@ for i, omega in enumerate(ap.allomegas):
outfile = defaultprefix + ".npz" if a.output is None else a.output
np.savez(outfile, meta={**vars(a), 'qpms_version' : qpms.__version__()},
k_sph=k_sph_list, k_cart = k_cart_arr, E_cart=E_cart_list, E_sph=np.array(E_sph),
np.savez(outfile, meta=vars(a), k_sph=k_sph_list, k_cart = k_cart_arr, E_cart=E_cart_list, E_sph=np.array(E_sph),
σ_ext=σ_ext_arr,σ_abs=σ_abs_arr,σ_scat=σ_scat_arr,
σ_ext_ir=σ_ext_arr_ir,σ_abs_ir=σ_abs_arr_ir,σ_scat_ir=σ_scat_arr_ir, omega=ap.allomegas, wavenumbers=wavenumbers
)
@ -232,7 +231,6 @@ if a.plot or (a.plot_out is not None):
pdf.savefig(fig)
plt.close(fig)
annotate_pdf_metadata(pdf, scriptname="finiterectlat-scatter.py")
exit(0)

View File

@ -1,7 +1,7 @@
#!/usr/bin/env python3
import math
from qpms.argproc import ArgParser, annotate_pdf_metadata
from qpms.argproc import ArgParser
ap = ArgParser(['rectlattice2d', 'single_particle', 'single_lMax', 'omega_seq'])
ap.add_argument("-o", "--output", type=str, required=False, help='output path (if not provided, will be generated automatically)')
@ -93,7 +93,7 @@ for i, omega in enumerate(omegas):
SVs[iri][i] = np.linalg.svd(ImTW_packed, compute_uv = False)
outfile = defaultprefix + ".npz" if a.output is None else a.output
np.savez(outfile, meta={**vars(a), 'qpms_version' : qpms.__version__()}, omegas=omegas, wavenumbers=wavenumbers, SVs=np.concatenate(SVs, axis=-1), irrep_names=ss1.irrep_names, irrep_sizes=ss1.saecv_sizes, unitcell_area=ss.unitcell_volume
np.savez(outfile, meta=vars(a), omegas=omegas, wavenumbers=wavenumbers, SVs=np.concatenate(SVs, axis=-1), irrep_names=ss1.irrep_names, irrep_sizes=ss1.saecv_sizes, unitcell_area=ss.unitcell_volume
)
logging.info("Saved to %s" % outfile)
@ -102,7 +102,6 @@ if a.plot or (a.plot_out is not None):
import matplotlib
matplotlib.use('pdf')
from matplotlib import pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
fig = plt.figure()
ax = fig.add_subplot(111)
@ -120,9 +119,7 @@ if a.plot or (a.plot_out is not None):
ax.legend()
plotfile = defaultprefix + ".pdf" if a.plot_out is None else a.plot_out
with PdfPages(plotfile) as pdf:
pdf.savefig(fig)
annotate_pdf_metadata(pdf, scriptname='infiniterectlat-k0realfreqsvd.py')
fig.savefig(plotfile)
exit(0)

View File

@ -2,7 +2,7 @@
import math
pi = math.pi
from qpms.argproc import ArgParser, annotate_pdf_metadata
from qpms.argproc import ArgParser
ap = ArgParser(['rectlattice2d', 'single_particle', 'single_lMax', 'omega_seq_real_ng', 'planewave'])
ap.add_argument("-o", "--output", type=str, required=False, help='output path (if not provided, will be generated automatically)')
@ -106,7 +106,7 @@ with pgsl_ignore_error(15): # avoid gsl crashing on underflow
σ_abs_arr = σ_ext_arr - σ_scat_arr
outfile = defaultprefix + ".npz" if a.output is None else a.output
np.savez(outfile, meta={**vars(a), 'qpms_version' : qpms.__version__()}, dir_sph=dir_sph_list, k_cart = k_cart_arr, omega = ap.allomegas, E_cart = E_cart_list, wavenumbers= wavenumbers, σ_ext=σ_ext_arr,σ_abs=σ_abs_arr,σ_scat=σ_scat_arr, unitcell_area=ss.unitcell_volume
np.savez(outfile, meta=vars(a), dir_sph=dir_sph_list, k_cart = k_cart_arr, omega = ap.allomegas, E_cart = E_cart_list, wavenumbers= wavenumbers, σ_ext=σ_ext_arr,σ_abs=σ_abs_arr,σ_scat=σ_scat_arr, unitcell_area=ss.unitcell_volume
)
logging.info("Saved to %s" % outfile)
@ -214,6 +214,6 @@ if a.plot or (a.plot_out is not None):
pdf.savefig(fig)
plt.close(fig)
annotate_pdf_metadata(pdf)
exit(0)

View File

@ -1,7 +1,7 @@
#!/usr/bin/env python3
import math
from qpms.argproc import ArgParser, sfloat, annotate_pdf_metadata
from qpms.argproc import ArgParser, sfloat
ap = ArgParser(['const_real_background', 'lattice2d', 'multi_particle']) # TODO general analytical background
@ -112,7 +112,7 @@ res['inside_contour'] = inside_ellipse((res['eigval'].real, res['eigval'].imag),
#del res['omega'] If contour points are not needed...
#del res['ImTW'] # not if dbg=false anyway
outfile = defaultprefix + ".npz" if a.output is None else a.output
np.savez(outfile, meta={**vars(a), 'qpms_version' : qpms.__version__()}, empty_freqs=np.array(empty_freqs),
np.savez(outfile, meta=vars(a), empty_freqs=np.array(empty_freqs),
ss_positions=ss.positions, ss_fullvec_poffsets=ss.fullvec_poffsets,
ss_fullvec_psizes=ss.fullvec_psizes,
ss_bspecs_flat = np.concatenate(ss.bspecs),
@ -132,7 +132,6 @@ if a.plot or (a.plot_out is not None):
import matplotlib
matplotlib.use('pdf')
from matplotlib import pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
fig = plt.figure()
ax = fig.add_subplot(111,)
@ -156,11 +155,8 @@ if a.plot or (a.plot_out is not None):
ax.set_ylim([ymin-.1*yspan, ymax+.1*yspan])
ax.set_xlabel('$\hbar \Re \omega / \mathrm{eV}$')
ax.set_ylabel('$\hbar \Im \omega / \mathrm{meV}$')
plotfile = defaultprefix + ".pdf" if a.plot_out is None else a.plot_out
with PdfPages(plotfile) as pdf:
pdf.savefig(fig)
annotate_pdf_metadata(pdf, scriptname='lat2d_modes.py')
fig.savefig(plotfile)
exit(0)

View File

@ -1,7 +1,7 @@
#!/usr/bin/env python3
import math
from qpms.argproc import ArgParser, sfloat, annotate_pdf_metadata
from qpms.argproc import ArgParser, sfloat
ap = ArgParser(['background', 'lattice2d', 'multi_particle', 'omega_seq'])
@ -88,7 +88,7 @@ for i, omega in enumerate(omegas):
SVs[iri][i] = np.linalg.svd(ImTW_packed, compute_uv = False)
outfile = defaultprefix + ".npz" if a.output is None else a.output
np.savez(outfile, meta={**vars(a), 'qpms_version' : qpms.__version__()}, omegas=omegas, wavenumbers=wavenumbers, SVs=np.concatenate(SVs, axis=-1), irrep_names=ss1.irrep_names, irrep_sizes=ss1.saecv_sizes, unitcell_area=ss.unitcell_volume
np.savez(outfile, meta=vars(a), omegas=omegas, wavenumbers=wavenumbers, SVs=np.concatenate(SVs, axis=-1), irrep_names=ss1.irrep_names, irrep_sizes=ss1.saecv_sizes, unitcell_area=ss.unitcell_volume
)
logging.info("Saved to %s" % outfile)
@ -97,7 +97,6 @@ if a.plot or (a.plot_out is not None):
import matplotlib
matplotlib.use('pdf')
from matplotlib import pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
fig = plt.figure()
ax = fig.add_subplot(111)
@ -121,9 +120,7 @@ if a.plot or (a.plot_out is not None):
ax.legend()
plotfile = defaultprefix + ".pdf" if a.plot_out is None else a.plot_out
with PdfPages(plotfile) as pdf:
pdf.savefig(fig)
annotate_pdf_metadata(pdf, scriptname='lat2d_realfreqsvd.py')
fig.savefig(plotfile)
exit(0)

View File

@ -1,7 +1,7 @@
#!/usr/bin/env python3
import math
from qpms.argproc import ArgParser, annotate_pdf_metadata
from qpms.argproc import ArgParser
ap = ArgParser(['rectlattice2d', 'const_real_background', 'single_particle', 'single_lMax']) # const_real_background needed for calculation of the diffracted orders
ap.add_argument("-k", nargs=2, type=float, required=True, help='k vector', metavar=('K_X', 'K_Y'))
@ -116,7 +116,7 @@ res['refractive_index_internal'] = [emg(om).n for om in res['eigval']]
#del res['omega'] If contour points are not needed...
#del res['ImTW'] # not if dbg=false anyway
outfile = defaultprefix + ".npz" if a.output is None else a.output
np.savez(outfile, meta={**vars(a), 'qpms_version' : qpms.__version__()}, empty_freqs=np.array(empty_freqs),
np.savez(outfile, meta=vars(a), empty_freqs=np.array(empty_freqs),
ss_positions=ss.positions, ss_fullvec_poffsets=ss.fullvec_poffsets,
ss_fullvec_psizes=ss.fullvec_psizes,
ss_bspecs_flat = np.concatenate(ss.bspecs),
@ -133,7 +133,6 @@ if a.plot or (a.plot_out is not None):
import matplotlib
matplotlib.use('pdf')
from matplotlib import pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
fig = plt.figure()
ax = fig.add_subplot(111,)
@ -157,11 +156,8 @@ if a.plot or (a.plot_out is not None):
ax.set_ylim([ymin-.1*yspan, ymax+.1*yspan])
ax.set_xlabel('$\hbar \Re \omega / \mathrm{eV}$')
ax.set_ylabel('$\hbar \Im \omega / \mathrm{meV}$')
plotfile = defaultprefix + ".pdf" if a.plot_out is None else a.plot_out
with PdfPages(plotfile) as pdf:
pdf.savefig(fig)
annotate_pdf_metadata(pdf, scriptname="rectlat_simple_modes.py")
fig.savefig(plotfile)
exit(0)

View File

@ -1,107 +0,0 @@
# Installing QPMS on Android/AOSP (-based) systems
Yes, it is possible. Basically all you need is a device capable of running [Termux](https://termux.com/) with enough memory to build everything.
The following instructions have been tested with Termux version 0.118.0 on
[e/OS/ R development build on a Samsung Galaxy S10e](https://doc.e.foundation/devices/beyond0lte/install)
([e-1.0-r-20220526188878-dev-beyond0lte](https://images.ecloud.global/dev/beyond0lte/)).
Presumably, they should work also on stock Android as well, but who
in their right mind would run all the spyware by Google & al.?
Physical keyboard or [remote access](https://wiki.termux.com/wiki/Remote_Access) is strongly recommended. :D
## Get Termux
Just [install the Termux app from F-Droid or Github as per instructions](https://github.com/termux/termux-app#f-droid).
Open Termux; the following steps of these instructions are basically
just commands you need to type in Termux.
## Install prerequisities from termux repositories
```
pkg install python3 cmake git clang build-essential binutils
```
## Build and install GSL
```
curl -O https://www.nic.funet.fi/pub/gnu/ftp.gnu.org/pub/gnu/gsl/gsl-latest.tar.gz
tar xf gsl-latest.tar.gz
cd gsl-2.7.1
./configure --prefix=$PREFIX
make # add -j4 or so to make it faster on multicore systems
make install
cd -
```
## Build and install OpenBLAS
```
git clone https://github.com/xianyi/OpenBLAS.git
cd OpenBLAS
make
make PREFIX=$PREFIX install
cd -
```
### Workaround for "broken" setup.py script
The goal is to fix `setup.py` so that it finds the correct libraries automatically, but in the meantime, you can use this workaround to get the Python part of QPMS installed:
```
ln -s $PREFIX/lib/libopenblas.so $PREFIX/LIB/liblapacke.so
```
## Build and install Numpy
(Successful build requires the `MATHLIB` environmental variable set, otherwise linking will fail; see https://wiki.termux.com/wiki/Python.)
```
MATHLIB=m pip3 install numpy
```
### Install Sympy
```
pip3 install sympy
```
## Build and install QPMS
```
git clone https://repo.or.cz/qpms.git
cd qpms
cmake -DCMAKE_INSTALL_PREFIX=$PREFIX . # ugly, TODO make a separate build tree!
make install
python3 setup.py install
```
Hopefully, QPMS has installed successfully. At this point, you should be able
to import and use QPMS with some exceptions. First, there is some legacy code
in the `qpms.qpms_p` module (which is no longer imported automatically with
bare `import qpms`). You shouldn't need this unless you are trying to run some
historic Jupyter notebooks or other old custom scripts. It contains a scipy
dependence, and scipy is hard to get working in Android environment (as
it requires a Fortran compiler to build).
## Install matplotlib
If you try to run just `pip3 install matplotlib` in Termux, it might likely
fail when installing the `pillow` dependency.
First, according to [Termux wiki](https://wiki.termux.com/wiki/Python#Python_module_installation_tips_and_tricks),
pillow depends on `libpng` and `libjpeg-turbo`, which are fortunately
available in Termux packages.
Second, pillow instalation requires an additional environment variable
`LDFLAGS="-L/system/lib64"` to be set on 64-bit devices.
Hence:
```
pkg install libpng libjpeg-turbo
export LDFLAGS="-L/system/lib64" # on 64-bit devices
pip3 install matplotlib
```
After this step, you should be able to run the command-line scripts
from `misc/` directory and examples from `examples/` directory.

33
notes/VSWF_from_SSWF.md Normal file
View File

@ -0,0 +1,33 @@
VSWF expansions in terms of SSWF
================================
From
\cite necada_multiple-scattering_2021, eq. (2.19)
\f[
\wfkcout_{\tau lm}\left(\kappa (\vect r - \vect r_1) \right) =
\sum_{\tau'l'm'} \tropSr{\kappa(\vect r_2 - \vect r_1)}_{\tau l m;\tau'l'm} \wfkcreg_{\tau'l'm'}(\vect r -\vect r_2),
\qquad |\vect r -\vect r_2| < |\vect r_1 - \vect r_2|,
\f]
setting \f$ \vect r = \vect r_2\f$ and considering that
\f$ \wfkcreg_{\tau'l'm'}(\vect 0) \ne \vect 0 \f$ only for electric dipole waves (\f$ \tau = \mathrm{E}, l=1 \f$),
we have
\f[
\wfkcout_{\tau lm}\left(\kappa (\vect r - \vect r_1) \right) =
\sum_{m'} \tropSr{\kappa(\vect r - \vect r_1)}_{\tau l m;\mathrm{E}1m} \wfkcreg_{\mathrm{E}1m'}(\vect 0),
\qquad \vect r \ne \vect r_1 .
\f]
Combining this with
\cite necada_multiple-scattering_2021, eq. (2.25)
\f[
\tropSr{\vect d}_{\tau l m; \tau' l' m'} = \sum_{\lambda =|l-l'|+|\tau-\tau'|}^{l+l'}
C^{\lambda}_{\tau l m;\tau' l'm'} \underbrace{ \spharm{\lambda}{m-m'}(\uvec d) h_\lambda^{(1)}(d)}_{\sswfout_\lambda^{m-m'}(\vect d)},
\f]
we get
\f[
\wfkcout_{\tau lm}(\vect d) = \sum_{m'=-1}^1 \wfkcreg_{\mathrm{E}1m'}(\vect 0)
\sum_{\lambda=l-1+|\tau-\tau'|}^{l+1}
C^\lambda_{\tau l m;\mathrm{E}1m'} \sswfout_\lambda^{m-m'}(\vect d).
\f]
Note that the VSWF components in this expression are given in global "cartesian" basis,
*not* the local orthonormal basis derived from spherical coordinates.
(This is mostly desirable, because in lattices we need to work with flat coordinates anyway.)

View File

@ -1,6 +1,10 @@
VSWF conventions {#vswf_conventions}
====================================
*This page provides reference about the VSWF conventions used in the literature.
For VSWF convention specification in QPMS API, see
[SWF Conventions in QPMS](@ref swf_conventions_qpms).*
In general, the (transversal) VSWFs can be defined using (some) vector spherical harmonics
as follows: \f[
\wfm\pr{k\vect r}_{lm} = \sphbes_l(kr) \vshrot_{lm} (\uvec r),\\
@ -50,7 +54,7 @@ where the connection to negative orders is
\dlmfFer{\nu}{m}(x) = (-1)^m \frac{\Gamma\pr{\nu-m+1}}{\Gamma\pr{\nu+m+1}}\dlmfFer{\nu}{m}(x),\\
%\dlmfLeg{\nu}{m}(x) = \frac{\Gamma\pr{\nu-m+1}}{\Gamma\pr{\nu+m+1}}\dlmfLeg{\nu}{m}(x).\\
\f]
Note that there are called "Ferrers" functions in DLMF, while the "Legendre" functions have slightly
Note that they are called "Ferrers" functions in DLMF, while the "Legendre" functions have slightly
different meaning / conventions (Ferrers functions being defined for \f$ \abs{x} \le 1 \f$, whereas
Legendre for \f$ \abs{x} \ge 1 \f$. We will not use the DLMF "Legendre" functions here.

114
notes/conventions_qpms.md Normal file
View File

@ -0,0 +1,114 @@
SWF conventions in QPMS {#swf_conventions_qpms}
=================================================
*This page describes how (V)SWF conventions are specified
internally and in QPMS API. For a general overview of VSWF
conventions in the literature, see [VSWF Conventions](@ref vswf_conventions).*
Convention enumerator
---------------------
Most of the meaningful phase and normalisation conventions for spherical waves
can be specified by the enum type @ref qpms_normalisation_t.
The type can be also used to specify conventions that are currently not fully
supported in QPMS (such as those based on real spherical harmonics).
As an enum type, it does not cover all the conventions possibly imaginable,
but it does cover the most meaningful ones and most of those that can be found
in the literature.
(Most notably, it does not cover the “anti-normalisation”
that does appear in certain expressions in some literature where the spherical
harmonics contain unnormalised Legendre functions, so that the basis set of
of spherical harmonics has different norms for different signs of *m*
for the same *l*. This is a bad idea overall and an absolutely atrocious
approach for numerics. Do not use that.)
VSWF evaluation
---------------
Evaluation of VSWFs using qpms_vswf_fill(), qpms_eval_vswf(),
qpms_uvswf_fill() and other functions from vswf.h are evaluated as follows.
These fuctions take a @ref qpms_normalisation_t as an argument.
The Ferrers-Legendre functions and the π, τ functions are evaluated
by qpms_pitau_get(), which internally uses gsl_sf_legendre_deriv_array_e().
Note only the information about the Condon-Shortley
phase is passed to qpms_pitau_get() the result of this function
uses always the GSL_SF_LEGENDRE_SPHARM normalisation,
and possible normalisation and other phase factors are evaluated afterwards
using the inline functions
qpms_normalisation_factor_L_noCS(),
qpms_normalisation_factor_M_noCS(),
qpms_normalisation_factor_N_noCS().
Evaluation of vector spherical harmonics only with qpms_vecspharm_fill()
works similarly but TODO.
TODO reference to pi, tau.
VSWF translation operator evaluation
------------------------------------
In practice, translation operators are calculated by first creating
an instance of the qpms_trans_calculator structure, which contains
a table of constant normalisation factors for a given phase/normalisation
convention (it is assumed that the phase/normalisation conventions do not
change with the translation), and then calling
qpms_trans_calculator_get_AB_arrays()
(or others).
The precomputed factor table in qpms_trans_calculator_t contains a CS phase
related factor (via qpms_trans_normfac()).
Function qpms_trans_calculator_get_AB_arrays_buf() then calculates
the unnormalised (GSL_SF_LEGENDRE_NONE)
associated Legendre functions always with CS phase -1;
and the A, B arrays are filled (via qpms_trans_calculator_get_AB_arrays_precalcbuf())
by individual calls of qpms_trans_calculator_get_A_precalcbuf()
and qpms_trans_calculator_B_precalcbuf(), which basically just multiply
and sum the precalculated constant factors with the radial (Bessel),
polar (Legendre) and azimuthal (exponential/trigonometric) functions.
**This means that the normalisation and phase convention is fully embedded
in the constant factor tables, and nothing is calculated during "runtime".**
The "higher-level" qpms_trans_calculator_get_trans_array() currently
just calls qpms_trans_calculator_get_AB_arrays() and then reorders
the elements (using qpms_trans_array_from_AB()), asserting that
the normalisation conventions remain the same.
There seems to be an inconsistency between
qpms_trans_calculator_get_B_buf() and
qpms_trans_calculator_get_A_buf() on one hand, and
qpms_trans_calculator_get_AB_buf_p() and
qpms_trans_calculator_get_AB_arrays_buf() on the other.
While the latter two functions use always -1 as the CS phase,
the former two take it from the normalisation enumerator.
**Although the former two are probably used nowhere in the production,
this needs to be fixed.**
Lattice sums
------------
### Scalar SWFs
### Translation operators
Function qpms_trans_calculator_get_AB_arrays_e32_e()
first compute the scalar lattice sums (using ewald3_sigma_short(),
ewald3_sigma_long() and ewald3_sigma0() calls).
These are then transformed into the VSWF translation operator
elements in a similar manner as in
qpms_trans_calculator_get_A_precalcbuf() and
qpms_trans_calculator_get_B_precalcbuf(), although there some optical
differences (CHECK!).
### VSWFs

View File

@ -31,8 +31,15 @@ MathJax.Hub.Config({
spharm: ["{{Y_{\\mathrm{#1}}}_{#2}^{#3}}", 3, ""], // Spherical harmonics
spharmR: ["{{Y_{\\mathrm{#1}}}_{\\mathrm{#1}{#2}{#3}}", 4, ""], // Spherical harmonics
csphase: "\\mathsf{C_{CS}}", // Condon-Shortley phase
tropSrr: ["{{S^\\mathrm{#1}}\\pr{{#2} \\leftarrow {#3}}}", 3, ""], // Translation operator singular
tropRrr: ["{{R^\\mathrm{#1}}\\pr{{#2} \\leftarrow {#3}}}", 3, ""], // Translation operator regular
tropS: "{\\mathcal{S}}", // Translation operator singular
tropR: "{\\mathcal{R}}", // Translation operator regular
tropSr: ["{{\\mathcal{S}^\\mathrm{#1}}\\pr{{#2}}}", 2, ""], // Translation operator singular
tropRr: ["{{\\mathcal{R}^\\mathrm{#1}}\\pr{{#2}}}", 2, ""], // Translation operator regular
tropSrr: ["{{\\mathcal{S}^\\mathrm{#1}}\\pr{{#2} \\leftarrow {#3}}}", 3, ""], // Translation operator singular
tropRrr: ["{{\\mathcal{R}^\\mathrm{#1}}\\pr{{#2} \\leftarrow {#3}}}", 3, ""], // Translation operator regular
sswfout: "{\\psi}", // outgoing SSWF
sswfreg: "{\\phi}", // regular SSWF
// Kristensson's VSWFs, complex version (2014 notes)
wfkc: "{\\vect{y}}", // any wave

View File

@ -7,23 +7,16 @@ find_package(LAPACK REQUIRED)
# and other not very relevant warnings
set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -Wno-int-in-bool-context -Wno-comment")
# version file
include(GetGitRevisionDescription)
get_git_head_revision(GIT_REFSPEC GIT_SHA1)
configure_file("${CMAKE_CURRENT_SOURCE_DIR}/version.c.in" "${CMAKE_CURRENT_BINARY_DIR}/version.c" @ONLY)
list(APPEND SOURCES "${CMAKE_CURRENT_BINARY_DIR}/version.c" version.h)
#includes
set (DIRS ${GSL_INCLUDE_DIRS} ${GSLCBLAS_INCLUDE_DIRS})
include_directories(${DIRS})
add_library (qpms SHARED translations.c tmatrices.c vecprint.c vswf.c wigner.c ewald.c
ewaldsf.c pointgroups.c latticegens.c
lattices2d.c gaunt.c error.c legendre.c symmetries.c vecprint.c
bessel.c own_zgemm.c parsing.c scatsystem.c materials.c drudeparam_data.c
lll.c beyn.c trivialgroup.c version.c
lll.c beyn.c trivialgroup.c
translations_dbg.c scatsystem_dbg.c
)
use_c99()
@ -31,9 +24,9 @@ set(LIBS ${LIBS} ${GSL_LIBRARIES} ${GSLCBLAS_LIBRARIES})
target_link_libraries (qpms
${GSL_LIBRARIES}
${LAPACK_LIBRARIES}
${BLAS_LIBRARIES}
gsl
lapack
blas
${QPMS_AMOSLIB}
Faddeeva
)

View File

@ -1,3 +1,5 @@
from pkg_resources import get_distribution
__version__ = get_distribution('qpms').version
import os as __os
from sys import platform as __platform
@ -5,7 +7,7 @@ from sys import platform as __platform
import warnings as __warnings
try:
from .qpms_c import PointGroup, FinitePointGroup, FinitePointGroupElement, Particle, scatsystem_set_nthreads, ScatteringSystem, ScatteringMatrix, pitau, set_gsl_pythonic_error_handling, pgsl_ignore_error, gamma_inc, lll_reduce, qpms_library_version
from .qpms_c import PointGroup, FinitePointGroup, FinitePointGroupElement, Particle, scatsystem_set_nthreads, ScatteringSystem, ScatteringMatrix, pitau, set_gsl_pythonic_error_handling, pgsl_ignore_error, gamma_inc, lll_reduce
except ImportError as ex:
if __platform == "linux" or __platform == "linux2":
if 'LD_LIBRARY_PATH' not in __os.environ.keys():
@ -17,17 +19,12 @@ except ImportError as ex:
from .cyquaternions import CQuat, IRot3
from .cybspec import VSWFNorm, BaseSpec, default_bspec
from .cytmatrices import CTMatrix, TMatrixInterpolator, TMatrixGenerator
from .cytranslations import trans_calculator
#from .cytranslations import trans_calculator
from .cymaterials import MaterialInterpolator, EpsMu, LorentzDrudeModel, lorentz_drude, EpsMuGenerator
from .cycommon import dbgmsg_enable, dbgmsg_disable, dbgmsg_active, BesselType, VSWFType
from .cywaves import vswf_single
def __version__():
from pkg_resources import get_distribution
librev = qpms_library_version()
return get_distribution('qpms').version + (("lr:"+librev) if librev else "")
#from .qpms_p import * # don't import automatically (adds around 0.5 s delay and depends on scipy)
from .qpms_p import * # maybe don't import automatically in the future (adds around 0.5 s delay)
from .constants import *
# legacy code which brutally slows down the whole package init:

150
qpms/analysis.py Normal file
View File

@ -0,0 +1,150 @@
'''
Auxiliary functions to analyse results obtained using scripts in
the ../misc directory.
This needs to be kept in sync with the scripts and qpms/argproc.py.
'''
import math
import numpy as np
import qpms
from qpms.qpms_c import Particle
from qpms.cytmatrices import CTMatrix, TMatrixGenerator
from qpms.cybspec import BaseSpec
from qpms.cymaterials import EpsMu, EpsMuGenerator, LorentzDrudeModel, lorentz_drude
from qpms.symmetries import point_group_info
from qpms import FinitePointGroup, ScatteringSystem, eV, hbar
eh = eV / hbar
def cleanarray(a, atol=1e-8, copy=True, fullNaN=True):
a = np.array(a, copy=copy)
sieve = abs(a.real) < atol
a[sieve] = 1j * a[sieve].imag
sieve = abs(a.imag) < atol
a[sieve] = a[sieve].real
if fullNaN:
if np.all(a == 0):
a[...] = 0
return a
def nicerot(a, atol=1e-10, copy=True):
'''
Gives an array a "nice" phase.
'''
a = np.array(a, copy=copy)
i = np.argmax(abs(a))
a = a / a[i] * abs(a[i])
return a
def _meta_matspec2emg(matspec):
if isinstance(matspec, (float, complex)): # number is interpreted as refractive index
return EpsMuGenerator(EpsMu(matspec**2))
else:
return EpsMuGenerator(lorentz_drude[matspec])
class FiniteRectLatAnalysis:
'''Recreates the scattering system structure from finiterectlat-modes.py or
finiterectlat-constant-driving.py'''
def __init__(self, data):
meta = data['meta'][()]
scatter = 'symmetry_adapted' in meta.keys()# finiterectlat-constant-driving ; TODO unify
if scatter:
#thegroup = meta['symmetry_adapted'] # always D2h, this option does something else!
thegroup = 'D2h'
else: # finiterectlat-modes
thegroup = 'D2h' if meta['D2'] else 'D4h'
Nx, Ny = meta['size']
px, py = meta['period']
orig_x = (np.arange(Nx/2) + (0 if (Nx % 2) else .5)) * px
orig_y = (np.arange(Ny/2) + (0 if (Ny % 2) else .5)) * py
orig_xy = np.stack(np.meshgrid(orig_x, orig_y), axis = -1)
bspec = BaseSpec(lMax = meta['lMax'])
nelem = len(bspec)
bg = _meta_matspec2emg(meta['background'])
fg = _meta_matspec2emg(meta['material'])
if meta.get('height', None) is None:
tmgen = TMatrixGenerator.sphere(bg, fg, meta['radius'])
else:
tmgen = TMatrixGenerator.cylinder(bg, fg, meta['radius'], meta['height'], lMax_extend=meta['lMax_extend'])
particles = [Particle(orig_xy[i], tmgen, bspec) for i in np.ndindex(orig_xy.shape[:-1])]
sym = FinitePointGroup(point_group_info[thegroup])
omega = data['omega'][()] if scatter else meta['centre'] * eh
self.ss, self.ssw = ScatteringSystem.create(particles, bg, omega, sym=sym)
ss1, ssw1 = ScatteringSystem.create([Particle([0,0,0], tmgen, bspec)], bg, omega, sym=sym)
self.ss1, self.ssw1 = ss1, ssw1
# per-particle projectors/transformation to symmetry-adapted bases
fvcs1 = np.empty((nelem, nelem), dtype=complex)
y = 0
iris1 = []
for iri1 in range(ss1.nirreps):
for j in range(ss1.saecv_sizes[iri1]):
pvc1 = np.zeros((ss1.saecv_sizes[iri1],), dtype=complex)
pvc1[j] = 1
fvcs1[y] = ss1.unpack_vector(pvc1, iri1)
fvcs1[y] = cleanarray(nicerot(fvcs1[y], copy=False), copy=False)
y += 1
iris1.append(iri1)
# Mapping between ss particles and grid positions
positions = self.ss.positions.astype(np.float32)
# round to get rid of duplicities due to rounding errors
positions = positions.round(7-int(math.log(np.amax(abs(positions)),10)))
xpositions = np.unique(positions[:,0])
ypositions = np.unique(positions[:,1])
assert(len(xpositions) == Nx)
assert(len(ypositions == Ny))
# particle positions as integer indices
posmap = np.empty((positions.shape[0],2), dtype=int)
#invposmap = np.empty((Nx, Ny), dtype=int)
for i, pos in enumerate(positions):
posmap[i,0] = np.searchsorted(xpositions, positions[i,0])
posmap[i,1] = np.searchsorted(ypositions, positions[i,1])
#invposmap[posmap[i,0], posmap[i, 1]] = i
self.meta = meta
self.npzfile = data
self.fvcs1 = fvcs1
self.iris1 = np.array(iris1)
self.bspec = bspec
self.nelem = nelem
self.Nx, self.Ny, self.px, self.py = Nx, Ny, px, py
self.posmap = posmap
self.fullvec_poffsets = nelem * np.arange(Nx*Ny)
self.group = thegroup
if scatter:
self.typ = 'scatter'
self.fvcs1_fromdata = data['fvcs1']
self.iris1_fromdata = data['iris1']
self.positions_fromdata = data['positions']
self.driving_symmetry_adapted = data['symmetry_adapted']
else: # modes
self.typ = 'modes'
self.eigval = data['eigval']
self.neig = len(self.eigval)
self.eigvec = data['eigvec']
self.residuals = data['residuals']
self.ranktest_SV = data['ranktest_SV']
self.iri = data['iri'][()]
def fullvec2grid(self, fullvec, swapxy=False):
arr = np.empty((self.Nx, self.Ny, self.nelem), dtype=complex)
for pi, offset in enumerate(self.fullvec_poffsets):
ix, iy = self.posmap[pi]
arr[ix, iy] = fullvec[offset:offset+self.nelem]
return np.swapaxes(arr, 0, 1) if swapxy else arr
def saecv2grid(self, saecv, swapxy=False):
fullvec = self.ss.unpack_vector(saecv, iri=self.iri)
return self.fullvec2grid(fullvec, swapxy)
def eigvec_asgrid(self, i, swapxy = False):
'''
Returns i-th eigenmode as a (Nx, Ny, nelem)-shaped grid of excitation
coefficients.
'''
if self.iri is not None:
return self.saecv2grid(self.eigvec[i], swapxy)
else:
return self.fullvec2grid(self.eigvec[i], swapxy)

View File

@ -540,22 +540,3 @@ class ArgParser:
plist.extend([Particle(pos, t, bspec=bspec) for pos in poss])
return plist
#TODO perhaps move into another module
def annotate_pdf_metadata(pdfPages, scriptname=None, keywords=None, author=None, title=None, subject=None, **kwargs):
"""Adds QPMS version-related metadata to a matplotlib PdfPages object
Use before closing the PDF file.
"""
from .qpms_c import qpms_library_version
d = pdfPages.infodict()
d['Creator'] = "QPMS%s (lib rev. %s), https://qpms.necada.org" % (
"" if scriptname is None else (" "+scriptname), qpms_library_version())
if author is not None:
d['Author'] = author
if title is not None:
d['Title'] = title
if subject is not None:
d['Subject'] = subject
if keywords is not None:
d['Keywords'] = ' '.join(keywords)
d.update(kwargs)

View File

@ -3,18 +3,20 @@
*/
#ifndef BEYN_H
#define BEYN_H
#ifdef __cplusplus
extern "C" {
#endif
#include <stddef.h>
#include <complex.h>
/// User-supplied function that provides the (row-major) m × m matrix M(z) whose "roots" are to be found.
/** Pure C array version */
typedef int (*beyn_function_M_t)(complex double *target_M, size_t m, complex double z, void *params);
typedef int (*beyn_function_M_t)(_Complex double *target_M, size_t m, _Complex double z, void *params);
/// (optional) User-supplied function that, given \f$ \hat V \f$, calculates \f$ M(z)^{-1} \hat V \f$.
/** Pure C array version */
typedef int (*beyn_function_M_inv_Vhat_t)(complex double *target_M_inv_Vhat, size_t m, size_t l,
const complex double *Vhat, complex double z, void *params);
typedef int (*beyn_function_M_inv_Vhat_t)(_Complex double *target_M_inv_Vhat, size_t m, size_t l,
const _Complex double *Vhat, _Complex double z, void *params);
/// Complex plane integration contour structure.
typedef struct beyn_contour_t {
@ -26,15 +28,15 @@ typedef struct beyn_contour_t {
* It does not have to be a centre in some strictly defined sense,
* but it should be "somewhere around" where the contour is.
*/
complex double centre;
_Complex double centre;
/// Function testing that a point \a z lies inside the contour (optional).
_Bool (*inside_test)(struct beyn_contour_t *, complex double z);
complex double z_dz[][2]; ///< Pairs of contour points and derivatives in that points.
_Bool (*inside_test)(struct beyn_contour_t *, _Complex double z);
_Complex double z_dz[][2]; ///< Pairs of contour points and derivatives in that points.
} beyn_contour_t;
/// Complex plane elliptic integration contour with axes parallel to the real, imaginary axes.
/** Free using free(). */
beyn_contour_t *beyn_contour_ellipse(complex double centre, double halfax_re, double halfax_im, size_t npoints);
beyn_contour_t *beyn_contour_ellipse(_Complex double centre, double halfax_re, double halfax_im, size_t npoints);
typedef enum {
@ -47,23 +49,23 @@ typedef enum {
/// Complex plane "half-elliptic" integration contour with axes parallel to the real, imaginary axes.
/** Free using free(). */
beyn_contour_t *beyn_contour_halfellipse(complex double centre, double halfax_re, double halfax_im, size_t npoints,
beyn_contour_halfellipse_orientation or);
beyn_contour_t *beyn_contour_halfellipse(_Complex double centre, double halfax_re, double halfax_im, size_t npoints,
beyn_contour_halfellipse_orientation o);
/// Similar to halfellipse but with rounded corners.
beyn_contour_t *beyn_contour_kidney(complex double centre, double halfax_re, double halfax_im,
beyn_contour_t *beyn_contour_kidney(_Complex double centre, double halfax_re, double halfax_im,
double rounding, ///< Must be in interval [0, 0.5)
size_t n, beyn_contour_halfellipse_orientation or);
size_t n, beyn_contour_halfellipse_orientation o);
/// Beyn algorithm result structure (pure C array version).
typedef struct beyn_result_t {
size_t neig; ///< Number of eigenvalues found.
size_t vlen; ///< Vector space dimension (also the leading dimension of eigvec).
complex double *eigval;
complex double *eigval_err;
_Complex double *eigval;
_Complex double *eigval_err;
double *residuals;
complex double *eigvec; // Rows are the eigenvectors
_Complex double *eigvec; // Rows are the eigenvectors
double *ranktest_SV;
} beyn_result_t;
@ -83,4 +85,7 @@ beyn_result_t *beyn_solve(
double res_tol ///< (default: `0.0`) TODO DOC.
);
#ifdef __cplusplus
}
#endif
#endif // BEYN_H

View File

@ -1,26 +1,7 @@
"""
Constants and unit conversion, mostly for standalode usage.
Previously, QPMS used scipy.constants in the python parts.
Now the relevant ones are set explicitly here in order
to avoid dependency on the whole scipy (which needs a fortran
compiler to build.
The values are taken from scipy / "2018 CODATA recommended values".
They slightly differ from the constants in GSL that are used
in the C code. It would be desirable to use the same source,
hence the values here might be subject to change in future versions.
"""
epsilon_0 = ε_0 = 8.8541878128e-12 # ± 0.0000000013e-12 F m^-1
c = speed_of_light = 299792458.
eV = e = elementary_charge = 1.602176487e-19 # ± 0000000040e-19 C
hbar = = 1.054571800e-34 # ± 0.000000013e-34 J s
mu_0 = μ_0 = 1.25663706212e-6 # ± 0.00000000019 e-6 N A^-2
from math import pi
π = pi
# unit conversions, mostly for standalone usage
# TODO avoid importing the "heavy" qpms parts
from scipy.constants import epsilon_0 as ε_0, c, pi as π, e as eV, hbar, hbar as , mu_0 as μ_0
pi = π
μm = 1e-6
nm = 1e-9
# "SCUFF FREQUENCY UNIT"

View File

@ -20,6 +20,12 @@ class BesselType(enum.IntEnum):
HANKEL_PLUS = QPMS_HANKEL_PLUS
HANKEL_MINUS = QPMS_HANKEL_MINUS
class EwaldPart(enum.IntEnum):
LONG_RANGE = QPMS_EWALD_LONG_RANGE
SHORT_RANGE = QPMS_EWALD_SHORT_RANGE
FULL = QPMS_EWALD_FULL
ZEROTERM = QPMS_EWALD_0TERM
class PointGroupClass(enum.IntEnum):
CN = QPMS_PGS_CN
S2N = QPMS_PGS_S2N
@ -138,6 +144,46 @@ def get_y_mn_unsigned(int nmax):
i = i + 1
return(ymn_plus, ymn_minus)
def get_nelem_scalar(lMax):
# TODO DOC, exceptions etc.
#return qpms_lMax2nelem(lMax)
return lMax * (lMax + 2) + 1 # = (lMax + 1)**2
## Auxillary function for retrieving the "meshgrid-like" indices; inc. nmax
@cython.boundscheck(False)
def get_mn_y_scalar(int nmax):
"""
Auxillary function for retreiving the 'meshgrid-like' indices from the flat indexing;
inc. nmax.
('y to mn' conversion)
Parameters
----------
nmax : int
The maximum order to which the SSWFs / Legendre functions etc. will be evaluated.
Returns
-------
output : (m, n)
Tuple of two arrays of type np.array(shape=(nmax*nmax + 2*nmax + 1), dtype=np.int),
where [(m[y],n[y]) for y in range(nmax*nmax + 2*nmax + 1)] covers all possible
integer pairs n >= 0, -n <= m <= n.
"""
cdef Py_ssize_t nelems = (nmax + 1)**2
cdef np.ndarray[np.int_t,ndim=1] m_arr = np.empty([nelems], dtype=np.int)
cdef np.ndarray[np.int_t,ndim=1] n_arr = np.empty([nelems], dtype=np.int)
cdef Py_ssize_t i = 0
cdef np.int_t n, m
for n in range(0,nmax+1):
for m in range(-n,n+1):
m_arr[i] = m
n_arr[i] = n
i = i + 1
return (m_arr, n_arr)
def tlm2uvswfi(t, l, m):
''' TODO doc

View File

@ -1,12 +1,36 @@
from .qpms_cdefs cimport *
from .qpms_c cimport ScatteringSystem, _ScatteringSystemAtOmega, _ScatteringSystemAtOmegaK
from .cycommon import EwaldPart, BesselType
from libc.stdlib cimport malloc, free, calloc
import numpy as np
from cython.parallel cimport prange, parallel
from cython cimport boundscheck, wraparound
cdef extern from "ewald.h":
void ewald3_2_sigma_long_Delta(cdouble *target, double *err, int maxn, cdouble x, int xbranch, cdouble z)
void ewald3_2_sigma_long_Delta_series(cdouble *target, double *err, int maxn, cdouble x, int xbranch, cdouble z)
void ewald3_2_sigma_long_Delta_recurrent(cdouble *target, double *err, int maxn, cdouble x, int xbranch, cdouble z, bint bigimz)
int complex_gamma_inc_e(double a, cdouble x, int xbranch, qpms_csf_result *result)
extern int ewald_factor_ipow_l
extern int ewald_factor_ipow_m
cdef extern from "translations_dbg.h":
int qpms_trans_calculator_test_sswf(const qpms_trans_calculator *c,
cdouble *dest, const csph_t kdlj, const qpms_bessel_t J) nogil
int qpms_trans_calculator_test_sswf_lc3p(const qpms_trans_calculator *c,
cdouble *dest, const cdouble k, const cart3_t r, const qpms_bessel_t J) nogil
int qpms_trans_calculator_test_sswf_e32(const qpms_trans_calculator *c,
cdouble * const sigmas_total, double * serr_total,
const double eta, const cdouble k,
const cart2_t b1, const cart2_t b2,
const cart2_t beta,
const cart3_t particle_shift,
double maxR, double maxK,
const qpms_ewald_part parts) nogil
cdef extern from "scatsystem_dbg.h":
int qpms_scatsyswk_test_sswf_basis_e(cdouble *target, const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t typ, cart3_t evalpoint, qpms_ewald_part parts) nogil
def e32_Delta(int maxn, cdouble x, cdouble z, int xbranch = 0, get_err=True, method='auto'):
cdef np.ndarray[double, ndim=1] err_np
@ -29,9 +53,117 @@ def e32_Delta(int maxn, cdouble x, cdouble z, int xbranch = 0, get_err=True, met
else:
return target_np
def scatsyswk_sswf_basis(_ScatteringSystemAtOmegaK sswk, evalpos, btyp=QPMS_HANKEL_PLUS, ewaldparts=EwaldPart.FULL):
"""Evaluate Ewald-summed scalar SWFs that are used to compose the tranlation operators.
Parameters
----------
evalpos: array_like of floats and shape (..., 3)
Evaluation points in cartesian coordinates.
Returns
-------
ndarray of complex, with the shape `evalpos.shape[:-1] + (particle_count, pnelem)`
"""
if(btyp != QPMS_HANKEL_PLUS):
raise NotImplementedError("Only first kind Bessel function-based fields are supported")
cdef qpms_bessel_t btyp_c = BesselType(btyp)
cdef qpms_ewald_part ewaldparts_c = EwaldPart(ewaldparts)
evalpos = np.array(evalpos, dtype=float, copy=False)
if evalpos.shape[-1] != 3:
raise ValueError("Last dimension of evalpos has to be 3")
cdef ScatteringSystem ss = sswk.ssw_pyref.ss_pyref
cdef qpms_scatsys_t *ss_c = ss.s
cdef qpms_scatsys_at_omega_k_t *sswk_c = sswk.rawpointer()
cdef const qpms_trans_calculator *c = ss_c[0].c
cdef qpms_l_t sswf_lMax = 2*c[0].lMax+1
cdef qpms_y_t pnelem = qpms_lMax2nelem_sc(sswf_lMax)
cdef qpms_ss_pi_t particle_count = sswk.ssw_pyref.ss_pyref.particle_count
cdef np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
cdef np.ndarray[complex, ndim=3] results = np.empty((evalpos_a.shape[0], particle_count, pnelem), dtype=complex)
cdef cdouble *res
cdef cart3_t pos
cdef Py_ssize_t i, pi, j
with nogil, wraparound(False), parallel():
res = <cdouble *> malloc(particle_count * pnelem * sizeof(cdouble))
for i in prange(evalpos_a.shape[0]):
pos.x = evalpos_a[i,0]
pos.y = evalpos_a[i,1]
pos.z = evalpos_a[i,2]
qpms_scatsyswk_test_sswf_basis_e(res, sswk_c, btyp_c, pos, ewaldparts_c)
for pi in range(particle_count):
for j in range(pnelem):
results[i,pi,j] = res[pi * pnelem + j]
free(res)
return results.reshape(evalpos.shape[:-1] + (particle_count, pnelem))
def scatsysw_eval_sswf(_ScatteringSystemAtOmega ssw, evalpos, btyp=QPMS_HANKEL_PLUS):
"""Evaluate sswfs by internal qpms_trans_calculator at given locations.
Parameters
----------
ssw: _ScatteringSystemAtOmega
We take the wavenumber and qpms_trans_calculator from this, nothing else.
(In a way, this is an overkill, but there is no up-to-date standalone cython
wrapper for qpms_trans_calculator.)
evalpos: array_like of floats and shape (..., 3)
Evaluation points in cartesian coordinates.
btyp: BesselType
Kind of Bessel functions to be used
Returns
-------
ndarray of complex, with the shape `evalpos.shape[:-1] + (pnelem,)`
where `pnelem` is TODO
"""
cdef qpms_bessel_t btyp_c = BesselType(btyp)
evalpos = np.array(evalpos, dtype=float, copy=False)
if evalpos.shape[-1] != 3:
raise ValueError("Last dimension of evalpos has to be 3")
cdef ScatteringSystem ss = ssw.ss_pyref
cdef qpms_scatsys_t *ss_c = ss.s
cdef const qpms_trans_calculator *c = ss_c[0].c
cdef qpms_l_t sswf_lMax = 2*c[0].lMax+1
cdef qpms_y_t pnelem = qpms_lMax2nelem_sc(sswf_lMax)
cdef cdouble k = ssw.wavenumber
cdef np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
cdef np.ndarray[complex,ndim=2] results = np.empty((evalpos_a.shape[0], pnelem), dtype=complex)
cdef cdouble *res
cdef cart3_t pos
cdef csph_t kdlj
cdef Py_ssize_t i, j
with nogil, wraparound(False), parallel():
res = <cdouble *> malloc(pnelem * sizeof(cdouble))
for i in prange(evalpos_a.shape[0]):
pos.x = evalpos_a[i,0]
pos.y = evalpos_a[i,1]
pos.z = evalpos_a[i,2]
kdlj = cart2csph(pos)
kdlj.r *= k
qpms_trans_calculator_test_sswf(c, res, kdlj, btyp_c)
for j in range(pnelem):
results[i,j] = res[j]
free(res)
return results.reshape(evalpos.shape[:-1] + (pnelem,))
def gamma_inc(double a, cdouble x, int xbranch=0):
cdef qpms_csf_result res
complex_gamma_inc_e(a, x, xbranch, &res)
return res.val
def ipow_l_factor(n = None):
global ewald_factor_ipow_l
if n is not None:
ewald_factor_ipow_l = n
return ewald_factor_ipow_l
def ipow_m_factor(n = None):
global ewald_factor_ipow_m
if n is not None:
ewald_factor_ipow_m = n
return ewald_factor_ipow_m

View File

@ -18,7 +18,7 @@ cdef class MaterialInterpolator:
cdef readonly double omegamin
cdef readonly double omegamax
cdef inline void *rawpointer(self):
return <void *>(self.interp)
return <void *>&(self.interp)
cdef class EpsMu:
cdef public qpms_epsmu_t em

View File

@ -9,7 +9,7 @@ cimport cython
import enum
import warnings
import os
from .constants import e as eV, hbar, c
from scipy.constants import e as eV, hbar, c
from libc.stdlib cimport malloc, free, calloc, abort
class EpsMuGeneratorType(enum.Enum):

View File

@ -286,6 +286,10 @@ cdef class TMatrixGenerator:
self.holder = what
self.g.function = qpms_tmatrix_generator_interpolator
self.g.params = <void*>(<TMatrixInterpolator?>self.holder).rawpointer()
elif what == 0:
self.holder = 0
self.g.function = qpms_tmatrix_generator_zero
self.g.params = <void*> 0
elif callable(what):
warnings.warn("Custom python T-matrix generators are an experimental feature. Also expect it to be slow.")
self.holder = what
@ -403,3 +407,15 @@ cdef class TMatrixGenerator:
EpsMuGenerator(outside), EpsMuGenerator(inside),
ArcFunction(__ArcCylinder(r, h)), *args, **kwargs))
@staticmethod
def dummy():
"""Returns a dummy T-matrix generator (returns a zero T-matrix).
Returns
-------
tmgen_dummy : TMatrixGenerator
"""
return tmgen_dummy
# pre-generate a dummy TMatrixGenerator (which returns zero T-matrix)
tmgen_dummy = TMatrixGenerator(0)

View File

@ -7,11 +7,23 @@
#include <complex.h>
#include "tiny_inlines.h"
#include "qpms_error.h"
#include "lattices.h"
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_sf_expint.h>
/*
* Control additional phase factors
* that could possibly affect the resulting VSWF
* shape.
* Currently not used in ewald3_1_z_sigma_long()
*/
int ewald_factor_ipow_l = 0;
int ewald_factor_ipow_m = 0;
// parameters for the quadrature of integral in (4.6)
#ifndef INTEGRATION_WORKSPACE_LIMIT
#define INTEGRATION_WORKSPACE_LIMIT 30000
@ -125,7 +137,7 @@ qpms_ewald3_constants_t *qpms_ewald3_constants_init(const qpms_l_t lMax /*, cons
c->s1_constfacs_base = malloc(s1_constfacs_sz * sizeof(complex double));
size_t s1_constfacs_sz_cumsum = 0;
for (qpms_y_t y = 0; y < c->nelem_sc; ++y) {
for (qpms_y_sc_t y = 0; y < c->nelem_sc; ++y) {
qpms_l_t n; qpms_m_t m; qpms_y2mn_sc_p(y, &m, &n);
if ((m + n) % 2 == 0) {
c->s1_constfacs[y] = c->s1_constfacs_base + s1_constfacs_sz_cumsum;
@ -161,7 +173,7 @@ qpms_ewald3_constants_t *qpms_ewald3_constants_init(const qpms_l_t lMax /*, cons
c->S1_constfacs = malloc((1+c->nelem_sc) * sizeof(complex double *));
//determine sizes
size_t S1_constfacs_sz = 0;
for (qpms_y_t y = 0; y < c->nelem_sc; ++y) {
for (qpms_y_sc_t y = 0; y < c->nelem_sc; ++y) {
qpms_l_t n; qpms_m_t m; qpms_y2mn_sc_p(y, &m, &n);
const qpms_l_t L_M = n - abs(m);
for(qpms_l_t j = 0; j <= L_M; ++j) { // outer sum
@ -176,7 +188,7 @@ qpms_ewald3_constants_t *qpms_ewald3_constants_init(const qpms_l_t lMax /*, cons
c->S1_constfacs_base = malloc(S1_constfacs_sz * sizeof(complex double));
size_t S1_constfacs_sz_cumsum = 0; // second count
for (qpms_y_t y = 0; y < c->nelem_sc; ++y) {
for (qpms_y_sc_t y = 0; y < c->nelem_sc; ++y) {
qpms_l_t n; qpms_m_t m; qpms_y2mn_sc_p(y, &m, &n);
const complex double yfactor = -2 * ipow(n+1) * M_SQRTPI
* factorial((n-m)/2) * factorial((n+m)/2);
@ -303,7 +315,7 @@ int ewald3_21_xy_sigma_long (
const bool k_is_real = (cimag(k) == 0);
assert((latdim & LAT_XYONLY) && (latdim & SPACE3D));
assert((latdim & LAT1D) || (latdim & LAT2D));
const qpms_y_t nelem_sc = c->nelem_sc;
const qpms_y_sc_t nelem_sc = c->nelem_sc;
assert(nelem_sc > 0);
const qpms_l_t lMax = c->lMax;
@ -435,7 +447,7 @@ int ewald3_21_xy_sigma_long (
for(qpms_m_t m = -n; m <= n; ++m) {
if((particle_shift.z == 0) && ((m+n) % 2 != 0)) // odd coefficients are zero.
continue;
const qpms_y_t y = qpms_mn2y_sc(m, n);
const qpms_y_sc_t y = qpms_mn2y_sc(m, n);
size_t constidx = 0; // constants offset
const complex double e_imalpha_pq = cexp(I*m*arg_pq);
complex double jsum, jsum_c; ckahaninit(&jsum, &jsum_c);
@ -444,7 +456,7 @@ int ewald3_21_xy_sigma_long (
assert((n-abs(m))/2 == c->s1_jMaxes[y]);
for(qpms_l_t j = 0; j <= c->s1_jMaxes[y]/*(n-abs(m))/2*/; ++j) { // FIXME </<= ?
complex double summand = rbeta_pq_div_k_pow[n-2*j]
* e_imalpha_pq * c->legendre0[gsl_sf_legendre_array_index(n,abs(m))] * min1pow_m_neg(m) // This line can actually go outside j-loop
* e_imalpha_pq * c->legendre0[gsl_sf_legendre_array_index(n,abs(m))] * min1pow(m) // This line can actually go outside j-loop
* gamma_pq_powm1[2*j]// * Gamma_pq[j] bellow (GGG) after error computation
* c->s1_constfacs[y][j];
if(err) {
@ -455,6 +467,7 @@ int ewald3_21_xy_sigma_long (
ckahanadd(&jsum, &jsum_c, summand);
}
jsum *= phasefac * factor1d; // PFC
jsum *= ipow(n * ewald_factor_ipow_l) * ipow(m * ewald_factor_ipow_m);
ckahanadd(target + y, target_c + y, jsum);
#ifdef EWALD_AUTO_CUTOFF
kahanadd(&lsum, &lsum_c, cabs(jsum));
@ -483,6 +496,7 @@ int ewald3_21_xy_sigma_long (
ckahanadd(&jsum, &jsum_c, jsummand);
}
jsum *= phasefac; // factor1d not here, off-axis sums not implemented/allowed.
jsum *= ipow(n * ewald_factor_ipow_l) * ipow(m * ewald_factor_ipow_m);
ckahanadd(target + y, target_c + y, jsum);
#ifdef EWALD_AUTO_CUTOFF
kahanadd(&lsum, &lsum_c, cabs(jsum));
@ -496,7 +510,7 @@ int ewald3_21_xy_sigma_long (
#ifdef EWALD_AUTO_CUTOFF
++li;
double cursum_min = INFINITY;
for (qpms_y_t y = 0; y < nelem_sc; ++y) {
for (qpms_y_sc_t y = 0; y < nelem_sc; ++y) {
const double a = cabs(target[y]);
if (a) // Skip possible symmetry-induced zeros; TODO this might need to be done in completely different way, perhaps on per-y basis
cursum_min = MIN(cursum_min, a);
@ -516,10 +530,10 @@ ewald3_21_xy_sigma_long_end_point_loop:
free(err_c);
free(target_c);
for(qpms_y_t y = 0; y < nelem_sc; ++y) // CFC common factor from above
for(qpms_y_sc_t y = 0; y < nelem_sc; ++y) // CFC common factor from above
target[y] *= commonfac;
if(err)
for(qpms_y_t y = 0; y < nelem_sc; ++y)
for(qpms_y_sc_t y = 0; y < nelem_sc; ++y)
err[y] *= commonfac;
return 0;
}
@ -549,7 +563,7 @@ int ewald3_1_z_sigma_long (
assert(particle_shift.x == 0 && particle_shift.y == 0);
const double beta_z = beta.z;
const double particle_shift_z = particle_shift_z;
const qpms_y_t nelem_sc = c->nelem_sc;
const qpms_y_sc_t nelem_sc = c->nelem_sc;
const qpms_l_t lMax = c->lMax;
// Manual init of the ewald summation targets
@ -600,7 +614,7 @@ int ewald3_1_z_sigma_long (
// TODO optimisations: all the j-dependent powers can be done for each j only once, stored in array
// and just fetched for each n
for(qpms_l_t n = 0; n <= lMax; ++n) {
const qpms_y_t y = qpms_mn2y_sc(0, n);
const qpms_y_sc_t y = qpms_mn2y_sc(0, n);
complex double jsum, jsum_c; ckahaninit(&jsum, &jsum_c);
double jsum_err, jsum_err_c; kahaninit(&jsum_err, &jsum_err_c); // TODO do I really need to kahan sum errors?
for(qpms_l_t j = 0; j <= n/2; ++j) {
@ -634,10 +648,10 @@ int ewald3_1_z_sigma_long (
free(err_c);
free(target_c);
for(qpms_y_t y = 0; y < nelem_sc; ++y) // CFC common factor from above
for(qpms_y_sc_t y = 0; y < nelem_sc; ++y) // CFC common factor from above
target[y] *= commonfac;
if(err)
for(qpms_y_t y = 0; y < nelem_sc; ++y)
for(qpms_y_sc_t y = 0; y < nelem_sc; ++y)
err[y] *= commonfac;
return 0;
}
@ -780,7 +794,7 @@ int ewald3_sigma_short(
{
const bool k_is_real = (cimag(k) == 0); // TODO check how the compiler optimises the loops
const double kreal = creal(k);
const qpms_y_t nelem_sc = c->nelem_sc;
const qpms_y_sc_t nelem_sc = c->nelem_sc;
const qpms_l_t lMax = c->lMax;
gsl_integration_workspace *workspace =
gsl_integration_workspace_alloc(INTEGRATION_WORKSPACE_LIMIT);
@ -902,7 +916,8 @@ int ewald3_sigma_short(
if(err)
kahanadd(err + y, err_c + y, cabs(leg * (prefacn / I) * R_pq_pown
* interr[n])); // TODO include also other errors
complex double thesummand = prefacn * R_pq_pown * leg * cintres[n] * e_beta_Rpq * e_imf * min1pow_m_neg(m);
complex double thesummand = prefacn * R_pq_pown * leg * cintres[n] * e_beta_Rpq * e_imf * min1pow(m);
thesummand *= ipow(n * ewald_factor_ipow_l) * ipow(m * ewald_factor_ipow_m);
ckahanadd(target + y, target_c + y, thesummand);
#ifdef EWALD_AUTO_CUTOFF
kahanadd(&lsum, &lsum_c, cabs(thesummand));
@ -916,7 +931,7 @@ int ewald3_sigma_short(
#ifdef EWALD_AUTO_CUTOFF
++li;
double cursum_min = INFINITY;
for (qpms_y_t y = 0; y < nelem_sc; ++y) {
for (qpms_y_sc_t y = 0; y < nelem_sc; ++y) {
const double a = cabs(target[y]);
if (a) // Skip possible symmetry-induced zeros; TODO this might need to be done in completely different way, perhaps on per-y basis
cursum_min = MIN(cursum_min, a);

View File

@ -26,22 +26,18 @@
#ifndef EWALD_H
#define EWALD_H
#ifdef __cplusplus
extern "C" {
#endif
#include <gsl/gsl_sf_result.h>
#include <stdlib.h>
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_errno.h>
#include <math.h> // for inlined lilgamma
#include <complex.h>
#include "qpms_types.h"
#include "lattices.h"
typedef enum {
QPMS_EWALD_LONG_RANGE = 1,
QPMS_EWALD_SHORT_RANGE = 2,
QPMS_EWALD_0TERM = 4,
QPMS_EWALD_FULL = QPMS_EWALD_LONG_RANGE | QPMS_EWALD_SHORT_RANGE | QPMS_EWALD_0TERM,
} qpms_ewald_part;
#include "lattices_types.h"
#include <complex.h>
/// Use this handler to ignore underflows of incomplete gamma.
@ -55,11 +51,11 @@ gsl_error_handler_t IgnoreUnderflowsGSLErrorHandler;
*/
typedef struct qpms_ewald3_constants_t {
qpms_l_t lMax;
qpms_y_t nelem_sc;
qpms_y_sc_t nelem_sc;
/// The values of maximum \a j's in the long-range part summation, `[(l-|m|/2)]`.
qpms_l_t *s1_jMaxes;
/// The constant factors for the long range part of a 2D Ewald sum.
complex double **s1_constfacs; // indices [y][j] where j is same as in [1, (4.5)]
_Complex double **s1_constfacs; // indices [y][j] where j is same as in [1, (4.5)]
/* These are the actual numbers now: (in the EWALD32_CONSTANTS_AGNOSTIC version)
* for m + n EVEN:
*
@ -73,7 +69,7 @@ typedef struct qpms_ewald3_constants_t {
*
* s1_constfacs[y(m,n)][j] = 0
*/
complex double *s1_constfacs_base; ///< Internal pointer holding memory for the 2D Ewald sum constant factors.
_Complex double *s1_constfacs_base; ///< Internal pointer holding memory for the 2D Ewald sum constant factors.
// similarly for the 1D z-axis aligned case; now the indices are [n][j] (as m == 0)
/// The constant factors for the long range part of a 1D Ewald sum along the \a z axis.
/** If the summation points lie along a different direction, use the formula for
@ -86,7 +82,7 @@ typedef struct qpms_ewald3_constants_t {
// TODO indexing mechanisms
/// The constant factors for the long range part of a 2D Ewald sum.
complex double **S1_constfacs; // indices [y][j] where j is same as in [1, (4.5)]
_Complex double **S1_constfacs; // indices [y][j] where j is same as in [1, (4.5)]
/* These are the actual numbers now: (in the EWALD32_CONSTANTS_AGNOSTIC version)
* for m + n EVEN:
*
@ -100,7 +96,7 @@ typedef struct qpms_ewald3_constants_t {
*
* S1_constfacs[y(m,n)][j] = 0
*/
complex double *S1_constfacs_base; ///< Internal pointer holding memory for the 2D Ewald sum constant factors.
_Complex double *S1_constfacs_base; ///< Internal pointer holding memory for the 2D Ewald sum constant factors.
/// The constant factors for the long range part of a 1D Ewald sum along the \a z axis.
/** If the summation points lie along a different direction, use the formula for
* 2D sum with additional factor of
@ -108,7 +104,7 @@ typedef struct qpms_ewald3_constants_t {
*/
complex double **s1_constfacs_1Dz;
_Complex double **s1_constfacs_1Dz;
/* These are the actual numbers now:
* s1_constfacs_1Dz[n][j] =
*
@ -116,7 +112,7 @@ typedef struct qpms_ewald3_constants_t {
* --------------------------
* j! * 2**(2*j) * (n - 2*j)!
*/
complex double *s1_constfacs_1Dz_base; ///<Internal pointer holding memory for the 1D Ewald sum constant factors.
_Complex double *s1_constfacs_1Dz_base; ///<Internal pointer holding memory for the 1D Ewald sum constant factors.
double *legendre0; /* now with GSL_SF_LEGENDRE_NONE normalisation, because this is what is
* what the multipliers from translations.c count with.
@ -139,13 +135,13 @@ void qpms_ewald3_constants_free(qpms_ewald3_constants_t *);
/// Structure for holding complex-valued result of computation and an error estimate.
/** Similar to gsl_sf_result, but with complex val. */
typedef struct qpms_csf_result {
complex double val; ///< Calculation result.
_Complex double val; ///< Calculation result.
double err; ///< Error estimate.
} qpms_csf_result;
// [1, (A.9)]
static inline complex double lilgamma(double t) {
static inline _Complex double lilgamma(double t) {
t = fabs(t);
if (t >= 1)
return sqrt(t*t - 1);
@ -154,8 +150,8 @@ static inline complex double lilgamma(double t) {
}
// [1, (A.8)], complex version of lilgamma()
static inline complex double clilgamma(complex double z) {
complex double a1 = z - 1, a2 = z + 1;
static inline _Complex double clilgamma(_Complex double z) {
_Complex double a1 = z - 1, a2 = z + 1;
// ensure -pi/2 < arg(z + 1) < 3*pi/2
if (creal(a2) < 0 && cimag(a2) <= 0)
a2 = -csqrt(a2);
@ -180,7 +176,7 @@ static inline complex double clilgamma(complex double z) {
* even if `z1 == z2`, because `-0 == 0` according to IEEE 754.
* The side of the branch cut can be determined using `signbit(creal(z))`.
*/
int cx_gamma_inc_series_e(double a, complex double z, qpms_csf_result * result);
int cx_gamma_inc_series_e(double a, _Complex double z, qpms_csf_result * result);
/// Incomplete Gamma function as continued fractions.
/**
@ -192,7 +188,7 @@ int cx_gamma_inc_series_e(double a, complex double z, qpms_csf_result * result);
* even if `z1 == z2`, because `-0 == 0` according to IEEE 754.
* The side of the branch cut can be determined using `signbit(creal(z))`.
*/
int cx_gamma_inc_CF_e(double a, complex double z, qpms_csf_result * result);
int cx_gamma_inc_CF_e(double a, _Complex double z, qpms_csf_result * result);
/// Incomplete gamma for complex second argument.
/**
@ -209,7 +205,7 @@ int cx_gamma_inc_CF_e(double a, complex double z, qpms_csf_result * result);
* Another than principal branch can be selected using non-zero \a m
* argument.
*/
int complex_gamma_inc_e(double a, complex double x,
int complex_gamma_inc_e(double a, _Complex double x,
/// Branch index.
/** If zero, the principal value is calculated.
* Other branches might be chosen using non-zero \a m.
@ -226,7 +222,7 @@ int complex_gamma_inc_e(double a, complex double x,
/// Exponential integral for complex second argument.
/** If x is (almost) positive real, it just uses gsl_sf_expint_En_e(). */
int complex_expint_n_e(int n, complex double x, qpms_csf_result *result);
int complex_expint_n_e(int n, _Complex double x, qpms_csf_result *result);
/// Hypergeometric 2F2, used to calculate some errors.
@ -245,15 +241,15 @@ int ewald32_sr_integral(double r, double k, double n, double eta, double *result
* unsuitable especially for big values of \a maxn.
*
*/
void ewald3_2_sigma_long_Delta(complex double *target, double *target_err, int maxn, complex double x,
int xbranch, complex double z);
void ewald3_2_sigma_long_Delta(_Complex double *target, double *target_err, int maxn, _Complex double x,
int xbranch, _Complex double z);
/// The Delta_n factor from [Kambe II], Appendix 3, used in 2D-in-3D long range sum.
/** This function always uses Kambe's (corrected) recurrent formula.
* For production, use ewald3_2_sigma_long_Delta() instead.
*/
void ewald3_2_sigma_long_Delta_recurrent(complex double *target, double *target_err, int maxn, complex double x,
int xbranch, complex double z, _Bool bigimz);
void ewald3_2_sigma_long_Delta_recurrent(_Complex double *target, double *target_err, int maxn, _Complex double x,
int xbranch, _Complex double z, _Bool bigimz);
/// The Delta_n factor from [Kambe II], Appendix 3, used in 2D-in-3D long range sum.
/** This function always uses Taylor expansion in \a z.
@ -263,26 +259,26 @@ void ewald3_2_sigma_long_Delta_recurrent(complex double *target, double *target_
* parameters maxn = 40, z = 0.5, x = -3. This might be related to the exponential growth
* of the error.
*/
void ewald3_2_sigma_long_Delta_series(complex double *target, double *target_err, int maxn, complex double x,
int xbranch, complex double z);
void ewald3_2_sigma_long_Delta_series(_Complex double *target, double *target_err, int maxn, _Complex double x,
int xbranch, _Complex double z);
// General functions acc. to [2], sec. 4.6 currently valid for 2D and 1D lattices in 3D space
/// The Ewald sum "self-interaction" term that appears in the lattice sums with zero (direct-space) Bravais lattice displacement.
int ewald3_sigma0(complex double *result, ///< Pointer to save the result (single complex double).
int ewald3_sigma0(_Complex double *result, ///< Pointer to save the result (single _Complex double).
double *err, ///< Pointer to save the result error estimate (single double).
const qpms_ewald3_constants_t *c, ///< Constant factors structure initialised by qpms_ewald3_constants_init().
double eta, ///< Ewald parameter.
complex double wavenumber ///< Wavenumber of the background medium.
_Complex double wavenumber ///< Wavenumber of the background medium.
);
/// Short-range part of outgoing scalar spherical wavefunctions' lattice sum \f$ \sigma_{l,m}^\mathrm{S}(\vect k,\vect s)\f$.
int ewald3_sigma_short(
complex double *target_sigmasr_y, ///< Target array for \f$ \sigma_{l,m}^\mathrm{S} \f$, must be `c->nelem_sc` long.
_Complex double *target_sigmasr_y, ///< Target array for \f$ \sigma_{l,m}^\mathrm{S} \f$, must be `c->nelem_sc` long.
double *target_sigmasr_y_err, ///< Target array for error estimates, must be `c->nelem_sc` long or `NULL`.
const qpms_ewald3_constants_t *c, ///< Constant factors structure initialised by qpms_ewald3_constants_init().
double eta, ///< Ewald parameter.
complex double wavenumber, ///< Wavenumber of the background medium.
_Complex double wavenumber, ///< Wavenumber of the background medium.
/// Lattice dimensionality.
/** Ignored apart from asserts and possible optimisations, as the SR formula stays the same. */
LatticeDimensionality latdim,
@ -293,7 +289,7 @@ int ewald3_sigma_short(
* In such case, it is the responsibility of the caller to deallocate
* the generator.
*/
PGen *pgen_R,
struct PGen *pgen_R,
/// Indicates whether pgen_R already generates shifted points.
/** If false, the behaviour corresponds to the old ewald32_sigma_short_points_and_shift(),
* so the function assumes that the generated points correspond to the unshifted Bravais lattice,
@ -310,11 +306,11 @@ int ewald3_sigma_short(
/// Long-range part of outgoing scalar spherical wavefunctions' lattice sum \f$ \sigma_{l,m}^\mathrm{L}(\vect k,\vect s)\f$.
int ewald3_sigma_long( // calls ewald3_21_sigma_long or ewald3_3_sigma_long, depending on latdim
complex double *target_sigmalr_y, ///< Target array for \f$ \sigma_{l,m}^\mathrm{L} \f$, must be `c->nelem_sc` long.
_Complex double *target_sigmalr_y, ///< Target array for \f$ \sigma_{l,m}^\mathrm{L} \f$, must be `c->nelem_sc` long.
double *target_sigmalr_y_err, ///< Target array for error estimates, must be `c->nelem_sc` long or `NULL`.
const qpms_ewald3_constants_t *c, ///< Constant factors structure initialised by qpms_ewald3_constants_init().
double eta, ///< Ewald parameter.
complex double wavenumber, ///< Wavenumber of the background medium.
_Complex double wavenumber, ///< Wavenumber of the background medium.
double unitcell_volume, ///< Volume of the (direct lattice) unit cell (with dimension corresponding to the lattice dimensionality).
/// Lattice dimensionality.
LatticeDimensionality latdim,
@ -325,7 +321,7 @@ int ewald3_sigma_long( // calls ewald3_21_sigma_long or ewald3_3_sigma_long, dep
* In such case, it is the responsibility of the caller to deallocate
* the generator.
*/
PGen *pgen_K,
struct PGen *pgen_K,
/// Indicates whether pgen_K already generates shifted points.
/** If false, the behaviour corresponds to the old ewald32_sigma_long_points_and_shift(),
* so the function assumes that the generated points correspond to the unshifted reciprocal Bravais lattice,
@ -339,4 +335,13 @@ int ewald3_sigma_long( // calls ewald3_21_sigma_long or ewald3_3_sigma_long, dep
cart3_t particle_shift
);
// If nonzero, adds an additional factor \f$ i^{nl} \f$ to the Ewald sum result (for debugging).
extern int ewald_factor_ipow_l;
// If nonzero, adds an additional factor \f$ i^{nm} \f$ to the Ewald sum result (for debubbing).
extern int ewald_factor_ipow_m;
#ifdef __cplusplus
}
#endif
#endif //EWALD_H

View File

@ -14,6 +14,7 @@
#include <stdbool.h>
#include <Faddeeva.h>
#include "tiny_inlines.h"
#include "qpms_error.h"
// Some magic constants

View File

@ -3,6 +3,10 @@
*/
#ifndef GAUNT_H
#define GAUNT_H
#ifdef __cplusplus
extern "C" {
#endif
#include <stdlib.h>
#define _GAUNT_H_MIN(x,y) (((x) > (y)) ? (y) : (x))
@ -30,4 +34,7 @@ double const * gaunt_table_retrieve_allq(int m, int n, int mu, int nu);
int gaunt_table_or_xu_fill(double *target, int m, int n, int mu, int nu);
#endif //GAUNT_PRECOMPILED
#ifdef __cplusplus
}
#endif
#endif //GAUNT_H

View File

@ -23,6 +23,9 @@
*/
#ifndef QPMS_GROUPS_H
#define QPMS_GROUPS_H
#ifdef __cplusplus
extern "C" {
#endif
#include "qpms_types.h"
#include <assert.h>
@ -35,7 +38,7 @@ struct qpms_finite_group_irrep_t {
/** The r-th row, c-th column of the representation of the i'th element is retrieved as
* m[i * dim * dim + r * dim + c]
*/
complex double *m;
_Complex double *m;
};
/// A point group with its irreducible representations and some metadata.
@ -93,4 +96,7 @@ qpms_iri_t qpms_finite_group_find_irrep_by_name(qpms_finite_group_t *G, char *na
extern const qpms_finite_group_t QPMS_FINITE_GROUP_TRIVIAL;
extern const qpms_finite_group_t QPMS_FINITE_GROUP_TRIVIAL_G;
#ifdef __cplusplus
}
#endif
#endif // QPMS_GROUPS_H

View File

@ -1,57 +1,97 @@
/*! \file indexing.h
* \brief Various index conversion functions.
*
* Index conventions in QPMS
* -------------------------
*
* Depending on circumstances (e.g. whether scalar or vector spherical wavefunctions
* are being used), QPMS uses two mappings between the multipole degree/order index
* pair (\a l, \a m), due to historical reasons ocassionally also noted as (\a n, \a m),
* and a flat non-negative integer index.
*
* Both mappings map the (\a l, \a m) pair in lexicographic order (with
* the basic constraint \f$ |m| \le l \f$), the difference is that the "scalar"
* mapping starts from \a l = 0, whereas the "vector" mapping starts from \a l = 1,
* so that for the zeroth elements we have
* * (0, 0) 0, (1, -1) 1 in the "scalar" mapping,
* * (1, -1) 0 in the "vector" mapping.
*
* The vector SWF flat index is typically denoted by `y` and its type is \ref qpms_y_t,
* whereas * the scalar SWF flat index is typically denoted by `y_sc` and its
* type is \ref qpms_y_sc_t.
*
* Moreover, there is another mapping that encodes also the electric/magnetic/scalar
* WSWF type index, see \ref qpms_uvswfi_t.
*
*/
#ifndef QPMS_INDEXING_H
#define QPMS_INDEXING_H
#ifdef __cplusplus
extern "C" {
#endif
#include "qpms_types.h"
#include <math.h>
#include "optim.h"
static inline qpms_y_t qpms_mn2y(qpms_m_t m, qpms_l_t n) {
/// Mapping from multipole degree, order to flat index (degree ≥ 1)
static inline qpms_y_t qpms_mn2y(
qpms_m_t m, ///< multipole order
qpms_l_t n ///< multipole degree (a.k.a. \a l); must be greater than 0
) {
return n * (n + 1) + m - 1;
}
/// Mapping from flat index to multipole degree (degree ≥ 1)
static inline qpms_lm_t qpms_y2n(qpms_y_t y) {
//return (sqrt(5+y)-2)/2; // the cast will truncate the fractional part, which is what we want
return sqrt(y+1);
}
/// Mapping from flat index, multipole degree (≥ 1) to multipole order
static inline qpms_m_t qpms_yn2m(qpms_y_t y, qpms_l_t n) {
return y-qpms_mn2y(0,n);
}
/// Mapping from flat index to multipole degree, order (degree ≥ 1)
static inline void qpms_y2mn_p(qpms_y_t y, qpms_m_t *m, qpms_l_t *n){
*m=qpms_yn2m(y,*n=qpms_y2n(y));
}
/// Number of spherical multipoles with `lmax` ≥ degree ≥ 1
static inline qpms_y_t qpms_lMax2nelem(qpms_l_t lmax){
return lmax * ((qpms_y_t)lmax + 2);
}
// Scalar versions: they have a place for the 0, 0 term in the beginning
static inline qpms_y_t qpms_mn2y_sc(qpms_m_t m, qpms_l_t n) {
/// Mapping from multipole degree, order to flat index (degree ≥ 0)
static inline qpms_y_sc_t qpms_mn2y_sc(
qpms_m_t m, ///< multipole order
qpms_l_t n ///< multipole degree (a.k.a. \a l); muste be greater or equal than 0
) {
return n * (n + 1) + m;
}
static inline qpms_lm_t qpms_y2n_sc(qpms_y_t y) {
/// Mapping from flat index to multipole degree (degree ≥ 0)
static inline qpms_lm_t qpms_y2n_sc(qpms_y_sc_t y) {
//return (sqrt(5+y)-2)/2; // the cast will truncate the fractional part, which is what we want
return sqrt(y);
}
static inline qpms_m_t qpms_yn2m_sc(qpms_y_t y, qpms_l_t n) {
/// Mapping from flat index, multipole degree (≥ 0) to multipole order
static inline qpms_m_t qpms_yn2m_sc(qpms_y_sc_t y, qpms_l_t n) {
return y-qpms_mn2y_sc(0,n);
}
static inline void qpms_y2mn_sc_p(qpms_y_t y, qpms_m_t *m, qpms_l_t *n){
/// Mapping from flat index to multipole degree, order (degree ≥ 0)
static inline void qpms_y2mn_sc_p(qpms_y_sc_t y, qpms_m_t *m, qpms_l_t *n){
*m=qpms_yn2m_sc(y,*n=qpms_y2n_sc(y));
}
static inline qpms_y_t qpms_lMax2nelem_sc(qpms_l_t lmax){
return lmax * ((qpms_y_t)lmax + 2) + 1;
/// Number of spherical multipoles with `lmax` ≥ degree ≥ 0
static inline qpms_y_sc_t qpms_lMax2nelem_sc(qpms_l_t lmax){
return lmax * ((qpms_y_sc_t)lmax + 2) + 1;
}
// TODO maybe enable crashing / validity control by macro definitions...
@ -69,7 +109,7 @@ static const qpms_uvswfi_t QPMS_UI_L00 = 0;
/** Returns a non-zero value if the u value is invalid. */
static inline qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u,
qpms_vswf_type_t *t, qpms_m_t *m, qpms_l_t *n) {
*t = u & 3;
*t = (qpms_vswf_type_t)(u & 3);
qpms_y_sc_t y_sc = u / 4;
qpms_y2mn_sc_p(y_sc, m, n);
// Test validity
@ -82,7 +122,7 @@ static inline qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u,
/** Does *not* allow for longitudinal waves. */
static inline qpms_errno_t qpms_uvswfi2ty(qpms_uvswfi_t u,
qpms_vswf_type_t *t, qpms_y_t *y) {
*t = u & 3;
*t = (qpms_vswf_type_t)(u & 3);
*y = u / 4 - 1;
if (QPMS_UNLIKELY(*t == 0 || *t == 3)) return QPMS_ERROR;
if (QPMS_UNLIKELY(*y < 0)) return QPMS_ERROR;
@ -95,7 +135,7 @@ static inline qpms_errno_t qpms_uvswfi2ty(qpms_uvswfi_t u,
*/
static inline qpms_errno_t qpms_uvswfi2ty_l(qpms_uvswfi_t u,
qpms_vswf_type_t *t, qpms_y_t *y) {
*t = u & 3;
*t = (qpms_vswf_type_t)(u & 3);
*y = u / 4 - 1;
if (QPMS_UNLIKELY(*t == 3)) return QPMS_ERROR;
if (QPMS_UNLIKELY(*y < 0)) return QPMS_ERROR;
@ -110,4 +150,7 @@ static inline qpms_m_t qpms_uvswfi2m(qpms_uvswfi_t u) {
}
#ifdef __cplusplus
}
#endif
#endif //QPMS_INDEXING_H

View File

@ -3,8 +3,9 @@
*/
#ifndef KAHANSUM_H
#define KAHANSUM_H
#include <complex.h>
#ifdef __cplusplus
extern "C" {
#endif
static inline void kahaninit(double * const sum, double * const compensation) {
*sum = 0;
@ -19,16 +20,19 @@ static inline void kahanadd(double *sum, double *compensation, double input) {
}
static inline void ckahaninit(complex double * const sum, complex double * const compensation) {
static inline void ckahaninit(_Complex double * const sum, _Complex double * const compensation) {
*sum = 0;
*compensation = 0;
}
static inline void ckahanadd(complex double *sum, complex double *compensation, complex double input) {
complex double compensated_input = input - *compensation;
complex double nsum = *sum + compensated_input;
static inline void ckahanadd(_Complex double *sum, _Complex double *compensation, _Complex double input) {
_Complex double compensated_input = input - *compensation;
_Complex double nsum = *sum + compensated_input;
*compensation = (nsum - *sum) - compensated_input;
*sum = nsum;
}
#ifdef __cplusplus
}
#endif
#endif //KAHANSUM_H

View File

@ -1,9 +1,15 @@
/*! \file lattices.h
* \brief Lattice point generators and lattice vector analysis / transformation.
*
* \bug Header file not C++ compatible.
*/
#ifndef LATTICES_H
#define LATTICES_H
#ifdef __cplusplus // FIXME Not C++ compatible. Include "lattices_types.h" for minimal necessary enum decls.
extern "C" {
#endif
#include "lattices_types.h"
#include <math.h>
#include <stdbool.h>
@ -30,22 +36,6 @@
*
*/
typedef enum LatticeDimensionality {
LAT1D = 1,
LAT2D = 2,
LAT3D = 4,
SPACE1D = 8,
SPACE2D = 16,
SPACE3D = 32,
LAT_1D_IN_3D = 33,
LAT_2D_IN_3D = 34,
LAT_3D_IN_3D = 40,
// special coordinate arrangements (indicating possible optimisations)
LAT_ZONLY = 64,
LAT_XYONLY = 128,
LAT_1D_IN_3D_ZONLY = 97, // LAT1D | SPACE3D | 64
LAT_2D_IN_3D_XYONLY = 162 // LAT2D | SPACE3D | 128
} LatticeDimensionality;
inline static bool LatticeDimensionality_checkflags(
LatticeDimensionality a, LatticeDimensionality flags_a_has_to_contain) {
@ -945,4 +935,7 @@ int honeycomb_lattice_gen_extend_to_steps(honeycomb_lattice_gen_t *g, int maxste
int honeycomb_lattice_gen_extend_to_r(honeycomb_lattice_gen_t *g, double r);
void honeycomb_lattice_gen_free(honeycomb_lattice_gen_t *g);
#ifdef __cplusplus
}
#endif
#endif // LATTICES_H

37
qpms/lattices_types.h Normal file
View File

@ -0,0 +1,37 @@
/*! \file lattices_types.h
* \brief Auxiliary lattice point generator and related enum declarations.
*
* This file contains necessary declarations needed in other header files.
* In contrast to lattices.h, this one is C++ compatible.
*
* \see lattices.h
*/
#ifndef LATTICES_TYPES_H
#define LATTICES_TYPES_H
#ifdef __cplusplus
extern "C" {
#endif
typedef enum LatticeDimensionality {
LAT1D = 1,
LAT2D = 2,
LAT3D = 4,
SPACE1D = 8,
SPACE2D = 16,
SPACE3D = 32,
LAT_1D_IN_3D = 33,
LAT_2D_IN_3D = 34,
LAT_3D_IN_3D = 40,
// special coordinate arrangements (indicating possible optimisations)
LAT_ZONLY = 64,
LAT_XYONLY = 128,
LAT_1D_IN_3D_ZONLY = 97, // LAT1D | SPACE3D | 64
LAT_2D_IN_3D_XYONLY = 162 // LAT2D | SPACE3D | 128
} LatticeDimensionality;
struct PGen;
#ifdef __cplusplus
}
#endif
#endif // LATTICES_TYPES_H

View File

@ -5,6 +5,7 @@
#include <string.h>
#include "materials.h"
#include "qpms_error.h"
#include <stdbool.h>
#define SQ(x) ((x)*(x))

View File

@ -3,7 +3,12 @@
*/
#ifndef QPMS_MATERIALS_H
#define QPMS_MATERIALS_H
#ifdef __cplusplus
extern "C" {
#endif
#include "qpms_types.h"
#include <complex.h>
#include <gsl/gsl_spline.h>
#ifndef SPEED_OF_LIGHT
@ -20,36 +25,36 @@ typedef struct qpms_epsmu_generator_t {
* qpms_permittivity_interpolator_epsmu_g(),
* qpms_lorentzdrude_epsmu_g().
*/
qpms_epsmu_t (*function) (complex double omega, const void *params);
qpms_epsmu_t (*function) (_Complex double omega, const void *params);
const void *params;
} qpms_epsmu_generator_t;
/// Convenience function for generating material properties at given frequency.
static inline qpms_epsmu_t qpms_epsmu_generator_eval(
qpms_epsmu_generator_t gen, complex double omega) {
qpms_epsmu_generator_t gen, _Complex double omega) {
return gen.function(omega, gen.params);
}
/// Constant optical property "generator" for qpms_epsmu_generator_t.
qpms_epsmu_t qpms_epsmu_const_g(complex double omega, ///< Frequency ignored.
qpms_epsmu_t qpms_epsmu_const_g(_Complex double omega, ///< Frequency ignored.
const void *epsmu ///< Points to the qpms_epsmu_t to be returned.
);
/// Gets refractive index of a material from its permeability and permittivity.
/** \f[ n = \sqrt{\mu_r \varepsilon_r} \f] */
static inline complex double qpms_refindex(qpms_epsmu_t em) {
static inline _Complex double qpms_refindex(qpms_epsmu_t em) {
return csqrt(em.eps * em.mu);
}
/// Gets wave number \a k from angular frequency and material permeability and permittivity.
/** \f[ k = \frac{n\omega}{c_0} = \frac{\omega\sqrt{\mu_r \varepsilon_r}}{c_0} \f] */
static inline complex double qpms_wavenumber(complex double omega, qpms_epsmu_t em) {
static inline _Complex double qpms_wavenumber(_Complex double omega, qpms_epsmu_t em) {
return qpms_refindex(em)*omega/SPEED_OF_LIGHT;
}
/// Gets (relative) wave impedance \f$ \eta_r \f$ from material permeability and permittivity.
/** \eta_r = \sqrt{\mu_r / \varepsilon_r} \f] */
static inline complex double qpms_waveimpedance(qpms_epsmu_t em) {
static inline _Complex double qpms_waveimpedance(qpms_epsmu_t em) {
return csqrt(em.mu / em.eps);
}
@ -67,7 +72,7 @@ typedef struct qpms_ldparams_triple_t {
* \f]
*/
typedef struct qpms_ldparams_t {
complex double eps_inf; ///< Permittivity at infinity.
_Complex double eps_inf; ///< Permittivity at infinity.
double omega_p; ///< Plasma frequency.
size_t n; ///< Number of "oscillators".
qpms_ldparams_triple_t data[]; ///< "Oscillator" parameters.
@ -86,14 +91,14 @@ extern const qpms_ldparams_t *const QPMS_LDPARAMS_PT; ///< Lorentz-Drude paramet
extern const qpms_ldparams_t *const QPMS_LDPARAMS_W ; ///< Lorentz-Drude parameters from \cite rakic_optical_1998 for tungsten.
/// Lorentz-Drude permittivity.
complex double qpms_lorentzdrude_eps(complex double omega, const qpms_ldparams_t *);
_Complex double qpms_lorentzdrude_eps(_Complex double omega, const qpms_ldparams_t *);
/// Lorentz-Drude optical properties, with relative permeability set always to one.
qpms_epsmu_t qpms_lorentzdrude_epsmu(complex double omega, const qpms_ldparams_t *);
qpms_epsmu_t qpms_lorentzdrude_epsmu(_Complex double omega, const qpms_ldparams_t *);
/// Lorentz-Drude optical properties, with relative permeability set always to one, compatible with qpms_epsmu_generator_t.
qpms_epsmu_t qpms_lorentzdrude_epsmu_g(
complex double omega,
_Complex double omega,
const void *ldparams ///< Lorentz-Drude parameters, in reality const qpms_ldparams_t *.
);
@ -125,14 +130,14 @@ qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_from_yml(
);
/// Evaluates interpolated material permittivity at a given angular frequency.
complex double qpms_permittivity_interpolator_eps_at_omega(
_Complex double qpms_permittivity_interpolator_eps_at_omega(
const qpms_permittivity_interpolator_t *interp, double omega_SI);
/// Evaluates interpolated material permittivity at a given angular frequency, qpms_epsmu_generator_t compatible version.
/** Permeability is always set to one. Imaginary part of omega is discarded.
*/
qpms_epsmu_t qpms_permittivity_interpolator_epsmu_g(
complex double omega, ///< Angular frequency. The imaginary part is ignored!
_Complex double omega, ///< Angular frequency. The imaginary part is ignored!
const void * interpolator ///< Interpolator of type qpms_permittivity_interpolator_t
);
@ -148,14 +153,17 @@ double qpms_permittivity_interpolator_omega_max(
void qpms_permittivity_interpolator_free(qpms_permittivity_interpolator_t *interp);
/// Relative permittivity from the Drude model.
static inline complex double qpms_drude_epsilon(
complex double eps_inf, ///< Relative permittivity "at infinity".
complex double omega_p, ///< Plasma frequency \f$ \omega_p \f$ of the material.
complex double gamma_p, ///< Decay constant \f$ \gamma_p \f$ of the material.
complex double omega ///< Frequency \f$ \omega \f$ at which the permittivity is evaluated.
static inline _Complex double qpms_drude_epsilon(
_Complex double eps_inf, ///< Relative permittivity "at infinity".
_Complex double omega_p, ///< Plasma frequency \f$ \omega_p \f$ of the material.
_Complex double gamma_p, ///< Decay constant \f$ \gamma_p \f$ of the material.
_Complex double omega ///< Frequency \f$ \omega \f$ at which the permittivity is evaluated.
) {
return eps_inf - omega_p*omega_p/(omega*(omega+I*gamma_p));
}
#ifdef __cplusplus
}
#endif
#endif //QPMS_MATERIALS_H

View File

@ -2,14 +2,18 @@
* \brief Convention-dependent coefficients for VSWFs.
*
* See also @ref qpms_normalisation_t and @ref vswf_conventions.
*
* \bug Header file not C++ compatible.
*/
#ifndef NORMALISATION_H
#define NORMALISATION_H
#ifdef __cplusplus //FIXME not C++ compatible yet (enum bit operations)
extern "C" {
#endif
#include "qpms_types.h"
#include "qpms_error.h"
#include <math.h>
#include <complex.h>
#include "indexing.h"
#include "optim.h"
@ -36,8 +40,8 @@ static inline double qpms_normalisation_normfactor(qpms_normalisation_t norm, qp
* This version ignores the Condon-Shortley phase bit (perhaps because the Condon-Shortley
* phase is already taken into account in a `gsl_sf_legendre_*_e()` call.)
*/
static inline complex double qpms_normalisation_factor_M_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
complex double fac = qpms_normalisation_normfactor(norm, l, m);
static inline _Complex double qpms_normalisation_factor_M_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
_Complex double fac = qpms_normalisation_normfactor(norm, l, m);
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_M_MINUS)) fac *= -1;
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_M_I)) fac *= I;
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac;
@ -51,8 +55,8 @@ static inline complex double qpms_normalisation_factor_M_noCS(qpms_normalisation
* Do not use if the C.-S. has already been taken into account e.g. in
* a `gsl_sf_legendre_*_e()` call.
*/
static inline complex double qpms_normalisation_factor_M(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
complex double fac = qpms_normalisation_factor_M_noCS(norm, l, m);
static inline _Complex double qpms_normalisation_factor_M(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
_Complex double fac = qpms_normalisation_factor_M_noCS(norm, l, m);
return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac;
}
@ -62,8 +66,8 @@ static inline complex double qpms_normalisation_factor_M(qpms_normalisation_t no
* This version ignores the Condon-Shortley phase bit (perhaps because the Condon-Shortley
* phase is already taken into account in a `gsl_sf_legendre_*_e()` call.)
*/
static inline complex double qpms_normalisation_factor_N_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
complex double fac = qpms_normalisation_normfactor(norm, l, m);
static inline _Complex double qpms_normalisation_factor_N_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
_Complex double fac = qpms_normalisation_normfactor(norm, l, m);
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_N_MINUS)) fac *= -1;
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_N_I)) fac *= I;
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac;
@ -77,14 +81,14 @@ static inline complex double qpms_normalisation_factor_N_noCS(qpms_normalisation
* Do not use if the C.-S. has already been taken into account e.g. in
* a `gsl_sf_legendre_*_e()` call.
*/
static inline complex double qpms_normalisation_factor_N(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
complex double fac = qpms_normalisation_factor_N_noCS(norm, l, m);
static inline _Complex double qpms_normalisation_factor_N(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
_Complex double fac = qpms_normalisation_factor_N_noCS(norm, l, m);
return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac;
}
/// Returns the factors of a electric basis VSWF divided by the factor of a magnetic VWFS of a given convention, compared to the reference one.
static inline complex double qpms_normalisation_factor_N_M(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
static inline _Complex double qpms_normalisation_factor_N_M(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
return qpms_normalisation_factor_N_noCS(norm, l, m)
/ qpms_normalisation_factor_M_noCS(norm, l, m);
}
@ -95,8 +99,8 @@ static inline complex double qpms_normalisation_factor_N_M(qpms_normalisation_t
* This version ignores the Condon-Shortley phase bit (perhaps because the Condon-Shortley
* phase is already taken into account in a `gsl_sf_legendre_*_e()` call.)
*/
static inline complex double qpms_normalisation_factor_L_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
complex double fac = qpms_normalisation_normfactor(norm, l, m);
static inline _Complex double qpms_normalisation_factor_L_noCS(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
_Complex double fac = qpms_normalisation_normfactor(norm, l, m);
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_L_MINUS)) fac *= -1;
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_L_I)) fac *= I;
if (QPMS_UNLIKELY(norm & QPMS_NORMALISATION_INVERSE)) fac = 1/fac;
@ -109,13 +113,13 @@ static inline complex double qpms_normalisation_factor_L_noCS(qpms_normalisation
* Do not use if the C.-S. has already been taken into account e.g. in
* a `gsl_sf_legendre_*_e()` call.
*/
static inline complex double qpms_normalisation_factor_L(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
complex double fac = qpms_normalisation_factor_L_noCS(norm, l, m);
static inline _Complex double qpms_normalisation_factor_L(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
_Complex double fac = qpms_normalisation_factor_L_noCS(norm, l, m);
return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac;
}
/// Returns the factors of a basis VSWF of a given convention compared to the reference convention.
static inline complex double qpms_normalisation_factor_uvswfi(const qpms_normalisation_t norm, qpms_uvswfi_t ui) {
static inline _Complex double qpms_normalisation_factor_uvswfi(const qpms_normalisation_t norm, qpms_uvswfi_t ui) {
qpms_vswf_type_t t; qpms_m_t m; qpms_l_t l;
qpms_uvswfi2tmn(ui, &t, &m, &l);
switch(t) {
@ -156,7 +160,7 @@ static inline qpms_normalisation_t qpms_normalisation_dual(qpms_normalisation_t
* 0 \quad \mbox{if } m>0. \\
* \f]
*/
static inline complex double qpms_spharm_azimuthal_part(qpms_normalisation_t norm, qpms_m_t m, double phi) {
static inline _Complex double qpms_spharm_azimuthal_part(qpms_normalisation_t norm, qpms_m_t m, double phi) {
switch(QPMS_EXPECT(norm, QPMS_NORMALISATION_DEFAULT)
& (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) {
case 0:
@ -194,7 +198,7 @@ static inline complex double qpms_spharm_azimuthal_part(qpms_normalisation_t nor
*
*
*/
static inline complex double qpms_spharm_azimuthal_part_derivative_div_m(qpms_normalisation_t norm, qpms_m_t m, double phi) {
static inline _Complex double qpms_spharm_azimuthal_part_derivative_div_m(qpms_normalisation_t norm, qpms_m_t m, double phi) {
if(m==0) return 0;
switch(QPMS_EXPECT(norm, QPMS_NORMALISATION_DEFAULT)
& (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) {
@ -283,4 +287,7 @@ static inline double qpms_normalisation_t_factor_abssquare(qpms_normalisation_t
}
#endif
#ifdef __cplusplus
}
#endif
#endif //NORMALISATION_H

View File

@ -3,6 +3,9 @@
*/
#ifndef QPMS_PARSING_H
#define QPMS_PARSING_H
#ifdef __cplusplus
extern "C" {
#endif
#include <stddef.h>
@ -57,4 +60,7 @@ size_t qpms_parse_doubles_fromfile(
const char *filepath //< File to read from, or NULL, "", "-" to read from stdin.
);
#ifdef __cplusplus
}
#endif
#endif // QPMS_PARSING_H

View File

@ -1,6 +1,7 @@
#include "pointgroups.h"
#include <search.h>
#include <stdlib.h>
#include <stdbool.h>
#include <assert.h>
double qpms_pg_quat_cmp_atol = QPMS_QUAT_ATOL;

View File

@ -3,6 +3,9 @@
*/
#ifndef POINTGROUPS_H
#define POINTGROUPS_H
#ifdef __cplusplus
extern "C" {
#endif
#include "qpms_error.h"
#include "quaternions.h"
@ -18,9 +21,9 @@ static inline _Bool qpms_pg_is_finite_axial(qpms_pointgroup_class cls) {
case QPMS_PGS_DN:
case QPMS_PGS_DND:
case QPMS_PGS_DNH:
return true;
return 1;
default:
return false;
return 0;
}
}
@ -81,4 +84,7 @@ _Bool qpms_pg_is_subgroup(qpms_pointgroup_t a, qpms_pointgroup_t b);
#ifdef __cplusplus
}
#endif
#endif //POINTGROUPS_H

View File

@ -1,4 +1,6 @@
from .qpms_cdefs cimport qpms_finite_group_t
from .qpms_cdefs cimport qpms_finite_group_t, qpms_scatsys_t, qpms_iri_t,\
qpms_scatsys_at_omega_t, qpms_scatsys_at_omega_k_t
from .cymaterials cimport EpsMuGenerator
cdef class FinitePointGroup:
cdef readonly bint owns_data
@ -6,3 +8,30 @@ cdef class FinitePointGroup:
cdef inline qpms_finite_group_t *rawpointer(self):
return self.G
cdef class ScatteringSystem:
cdef list tmgobjs # here we keep the references to occuring TMatrixFunctions (and hence BaseSpecs and TMatrixGenerators)
#cdef list Tmatrices # Here we keep the references to occuring T-matrices
cdef EpsMuGenerator medium_holder # Here we keep the reference to medium generator
cdef qpms_scatsys_t *s
cdef FinitePointGroup sym
cdef qpms_iri_t iri_py2c(self, iri, allow_None=?)
cdef class _ScatteringSystemAtOmegaK:
cdef qpms_scatsys_at_omega_k_t sswk
cdef _ScatteringSystemAtOmega ssw_pyref
cdef inline qpms_scatsys_at_omega_k_t *rawpointer(self):
return &self.sswk
cdef class _ScatteringSystemAtOmega:
cdef qpms_scatsys_at_omega_t *ssw
cdef ScatteringSystem ss_pyref
cdef inline qpms_scatsys_at_omega_t *rawpointer(self):
return self.ssw

0
qpms/qpms_c.py Normal file
View File

View File

@ -11,17 +11,14 @@ from .qpms_cdefs cimport *
from .cyquaternions cimport IRot3, CQuat
from .cybspec cimport BaseSpec
from .cycommon cimport make_c_string
from .cycommon import string_c2py, PointGroupClass, BesselType
from .cycommon import string_c2py, PointGroupClass, BesselType, EwaldPart
from .cytmatrices cimport CTMatrix, TMatrixFunction, TMatrixGenerator, TMatrixInterpolator
from .cymaterials cimport EpsMuGenerator, EpsMu
from libc.stdlib cimport malloc, free, calloc
import warnings
def qpms_library_version():
return "git_" + qpms_version_git_sha1.decode('utf-8')
from cython.parallel cimport prange, parallel
from cython cimport boundscheck, wraparound
from cython.parallel import prange, parallel
from cython import boundscheck, wraparound
# Set custom GSL error handler. N.B. this is obviously not thread-safe.
cdef char *pgsl_err_reason
@ -372,12 +369,6 @@ cdef class ScatteringSystem:
Currently, it does not have a standard constructor. Use the
ScatteringSystem.create() method instead.
'''
cdef list tmgobjs # here we keep the references to occuring TMatrixFunctions (and hence BaseSpecs and TMatrixGenerators)
#cdef list Tmatrices # Here we keep the references to occuring T-matrices
cdef EpsMuGenerator medium_holder # Here we keep the reference to medium generator
cdef qpms_scatsys_t *s
cdef FinitePointGroup sym
cdef qpms_iri_t iri_py2c(self, iri, allow_None = True):
if iri is None and allow_None:
return QPMS_NO_IRREP
@ -702,7 +693,7 @@ cdef class ScatteringSystem:
self.s, iri, 0)
return target_np
def pack_matrix(self, fullmatrix, iri):
def pack_matrix(self, fullmatrix, iri, version='normal'):
"""Converts (projects) a matrix into an irrep subspace.
Parameters
@ -734,10 +725,14 @@ cdef class ScatteringSystem:
cdef np.ndarray[np.complex_t, ndim=2] target_np = np.empty(
(rlen, rlen), dtype=complex, order='C')
cdef cdouble[:,::1] target_view = target_np
qpms_scatsys_irrep_pack_matrix(&target_view[0][0], &fullmatrix_view[0][0],
if version == 'stupid':
qpms_scatsys_irrep_pack_matrix_stupid(&target_view[0][0], &fullmatrix_view[0][0],
self.s, iri)
else:
qpms_scatsys_irrep_pack_matrix(&target_view[0][0], &fullmatrix_view[0][0],
self.s, iri)
return target_np
def unpack_matrix(self, packedmatrix, iri):
def unpack_matrix(self, packedmatrix, iri, version='normal'):
"""Unpacks an "irrep-packed" excitation coefficient vector to full coordinates.
Parameters
@ -770,7 +765,11 @@ cdef class ScatteringSystem:
cdef np.ndarray[np.complex_t, ndim=2] target_np = np.empty(
(flen, flen), dtype=complex, order='C')
cdef cdouble[:,::1] target_view = target_np
qpms_scatsys_irrep_unpack_matrix(&target_view[0][0], &packedmatrix_view[0][0],
if version == 'stupid':
qpms_scatsys_irrep_unpack_matrix_stupid(&target_view[0][0], &packedmatrix_view[0][0],
self.s, iri, 0)
else:
qpms_scatsys_irrep_unpack_matrix(&target_view[0][0], &packedmatrix_view[0][0],
self.s, iri, 0)
return target_np
@ -822,6 +821,16 @@ cdef class ScatteringSystem:
qpms_scatsys_periodic_build_translation_matrix_full(&target_view[0][0], self.s, wavenumber, &blochvector_c, eta)
return target
def irrep_transform_matrix(self, qpms_iri_t iri):
self.check_s()
cdef size_t rlen = self.saecv_sizes[iri]
cdef size_t fullen = self.fecv_size
cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
(rlen, fullen), dtype=complex, order='C')
cdef cdouble[:,::1] target_view = target
qpms_scatsys_irrep_transform_matrix(&target_view[0][0], self.s, iri)
return target
def translation_matrix_packed(self, cdouble wavenumber, qpms_iri_t iri, J = QPMS_HANKEL_PLUS):
self.check_s()
cdef size_t rlen = self.saecv_sizes[iri]
@ -965,8 +974,8 @@ cdef class ScatteringSystem:
cdef ccart3_t res
cdef cart3_t pos
cdef Py_ssize_t i
with nogil, parallel(), boundscheck(False), wraparound(False):
for i in prange(evalpos_a.shape[0]):
with nogil:
for i in range(evalpos_a.shape[0]):
pos.x = evalpos_a[i,0]
pos.y = evalpos_a[i,1]
pos.z = evalpos_a[i,2]
@ -1006,18 +1015,13 @@ def empty_lattice_modes_xy(EpsMu epsmu, reciprocal_basis, wavevector, double max
free(omegas_c)
return omegas
cdef class _ScatteringSystemAtOmegaK:
'''
Wrapper over the C qpms_scatsys_at_omega_k_t structure
This represents an infinite periodic system with a given (fixed) frequency and Bloch vector.
'''
cdef qpms_scatsys_at_omega_k_t sswk
cdef _ScatteringSystemAtOmega ssw_pyref
cdef qpms_scatsys_at_omega_k_t *rawpointer(self):
return &self.sswk
property eta:
"""Ewald parameter η"""
def __get__(self):
@ -1026,7 +1030,8 @@ cdef class _ScatteringSystemAtOmegaK:
self.sswk.eta = eta
@boundscheck(False)
def scattered_E(self, scatcoeffvector_full, evalpos, btyp=QPMS_HANKEL_PLUS):
def scattered_E(self, scatcoeffvector_full, evalpos, btyp=QPMS_HANKEL_PLUS,
ewaldparts = EwaldPart.FULL):
"""Evaluate electric field for a given excitation coefficient vector (periodic system)
Parameters
@ -1045,6 +1050,7 @@ cdef class _ScatteringSystemAtOmegaK:
if(btyp != QPMS_HANKEL_PLUS):
raise NotImplementedError("Only first kind Bessel function-based fields are supported")
cdef qpms_bessel_t btyp_c = BesselType(btyp)
cdef qpms_ewald_part ewaldparts_c = EwaldPart(ewaldparts)
evalpos = np.array(evalpos, dtype=float, copy=False)
if evalpos.shape[-1] != 3:
raise ValueError("Last dimension of evalpos has to be 3")
@ -1055,18 +1061,19 @@ cdef class _ScatteringSystemAtOmegaK:
cdef ccart3_t res
cdef cart3_t pos
cdef Py_ssize_t i
with nogil, wraparound(False), parallel():
for i in prange(evalpos_a.shape[0]):
with nogil:
for i in range(evalpos_a.shape[0]):
pos.x = evalpos_a[i,0]
pos.y = evalpos_a[i,1]
pos.z = evalpos_a[i,2]
res = qpms_scatsyswk_scattered_E(&self.sswk, btyp_c, &scv_view[0], pos)
res = qpms_scatsyswk_scattered_E_e(&self.sswk, btyp_c, &scv_view[0], pos, ewaldparts_c)
results[i,0] = res.x
results[i,1] = res.y
results[i,2] = res.z
return results.reshape(evalpos.shape)
def scattered_field_basis(self, evalpos, btyp=QPMS_HANKEL_PLUS):
@boundscheck(False)
def scattered_field_basis(self, evalpos, btyp=QPMS_HANKEL_PLUS, ewaldparts=EwaldPart.FULL):
# TODO examples
"""Evaluate scattered field "basis" (periodic system)
@ -1088,6 +1095,7 @@ cdef class _ScatteringSystemAtOmegaK:
if(btyp != QPMS_HANKEL_PLUS):
raise NotImplementedError("Only first kind Bessel function-based fields are supported")
cdef qpms_bessel_t btyp_c = BesselType(btyp)
cdef qpms_ewald_part ewaldparts_c = EwaldPart(ewaldparts)
cdef Py_ssize_t fecv_size = self.fecv_size
evalpos = np.array(evalpos, dtype=float, copy=False)
if evalpos.shape[-1] != 3:
@ -1095,20 +1103,20 @@ cdef class _ScatteringSystemAtOmegaK:
cdef np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
cdef np.ndarray[complex, ndim=3] results = np.empty((evalpos_a.shape[0], fecv_size, 3), dtype=complex)
cdef ccart3_t *res
res = <ccart3_t *> malloc(fecv_size*sizeof(ccart3_t))
cdef cart3_t pos
cdef Py_ssize_t i, j
with nogil, wraparound(False), parallel():
for i in prange(evalpos_a.shape[0]):
with nogil:
res = <ccart3_t *> malloc(fecv_size*sizeof(ccart3_t))
for i in range(evalpos_a.shape[0]):
pos.x = evalpos_a[i,0]
pos.y = evalpos_a[i,1]
pos.z = evalpos_a[i,2]
qpms_scatsyswk_scattered_field_basis(res, &self.sswk, btyp_c, pos)
qpms_scatsyswk_scattered_field_basis_e(res, &self.sswk, btyp_c, pos, ewaldparts_c)
for j in range(fecv_size):
results[i,j,0] = res[j].x
results[i,j,1] = res[j].y
results[i,j,2] = res[j].z
free(res)
free(res)
return results.reshape(evalpos.shape[:-1] + (self.fecv_size, 3))
property fecv_size:
@ -1120,9 +1128,6 @@ cdef class _ScatteringSystemAtOmega:
that keeps the T-matrix and background data evaluated
at specific frequency.
'''
cdef qpms_scatsys_at_omega_t *ssw
cdef ScatteringSystem ss_pyref
def check(self): # cdef instead?
if not self.ssw:
raise ValueError("_ScatteringSystemAtOmega's ssw-pointer not set. You must not use the default constructor; ScatteringSystem.create() instead")
@ -1151,9 +1156,6 @@ cdef class _ScatteringSystemAtOmega:
qpms_scatsysw_apply_Tmatrices_full(&target_view[0], &a_view[0], self.ssw)
return target_np
cdef qpms_scatsys_at_omega_t *rawpointer(self):
return self.ssw
def scatter_solver(self, iri=None, k=None):
self.check()
cdef _ScatteringSystemAtOmegaK sswk # used only for periodic systems
@ -1237,7 +1239,6 @@ cdef class _ScatteringSystemAtOmega:
@boundscheck(False)
def scattered_E(self, scatcoeffvector_full, evalpos, blochvector=None, btyp=QPMS_HANKEL_PLUS, bint alt=False):
# FIXME TODO this obviously does not work for periodic systems
"""Evaluate electric field for a given excitation coefficient vector
Parameters
@ -1274,19 +1275,85 @@ cdef class _ScatteringSystemAtOmega:
cdef ccart3_t res
cdef cart3_t pos
cdef Py_ssize_t i
with wraparound(False), nogil, parallel():
for i in prange(evalpos_a.shape[0]):
for i in range(evalpos_a.shape[0]):
pos.x = evalpos_a[i,0]
pos.y = evalpos_a[i,1]
pos.z = evalpos_a[i,2]
if alt:
res = qpms_scatsysw_scattered_E__alt(self.ssw, btyp_c, &scv_view[0], pos)
else:
res = qpms_scatsysw_scattered_E(self.ssw, btyp_c, &scv_view[0], pos)
results[i,0] = res.x
results[i,1] = res.y
results[i,2] = res.z
return results.reshape(evalpos.shape)
@boundscheck(False)
def scattered_field_basis(self, evalpos, blochvector=None, particle=None, btyp=QPMS_HANKEL_PLUS):
# TODO examples
# FIXME periodic case not implemented
"""Evaluate scattered field "basis"
This function enables the evaluation of "scattered" fields
generated by the system for many different excitation
coefficients vectors, without the expensive re-evaluation of the respective
translation operators for each excitation coefficient vector.
Parameters
----------
evalpos: array_like of floats and shape (..., 3)
Evaluation points in cartesian coordinates.
blochvector: array_like or None
Bloch vector, must be supplied (non-None) for periodic systems, else None.
particle_index: int or None (default), optional
A valid particle index; if specified, only the contribution of the given particle
is evaluated. Not applicable for periodic arrays.
btyp: BesselType, optional
Kind of the waves. Defaults to BesselType.HANKEL_PLUS.
Returns
-------
ndarray of complex, with the shape `evalpos.shape[:-1] + (n, 3)`
"Basis" fields at the positions given in `evalpos`, in cartesian coordinates.
`n` here stays either for `self.fecv_size` (if `particle==None`)
or `len(self.bspec_pi(particle))`.
"""
#if(btyp != QPMS_HANKEL_PLUS): # #TODO, IIRC not supported only for periodic systems
# raise NotImplementedError("Only first kind Bessel function-based fields are supported")
cdef qpms_bessel_t btyp_c = BesselType(btyp)
cdef Py_ssize_t basissize = self.fecv_size if particle is None else len(self.bspec_pi(particle))
cdef qpms_ss_pi_t pi = particle
evalpos = np.array(evalpos, dtype=float, copy=False)
if evalpos.shape[-1] != 3:
raise ValueError("Last dimension of evalpos has to be 3")
cdef np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
cdef np.ndarray[complex, ndim=3] results = np.empty((evalpos_a.shape[0], basissize, 3), dtype=complex)
cdef ccart3_t *res
cdef cart3_t pos
cdef Py_ssize_t i, j
with nogil:
res = <ccart3_t *> malloc(basissize*sizeof(ccart3_t)) # thread-local
for i in range(evalpos_a.shape[0]):
pos.x = evalpos_a[i,0]
pos.y = evalpos_a[i,1]
pos.z = evalpos_a[i,2]
if alt:
res = qpms_scatsysw_scattered_E__alt(self.ssw, btyp_c, &scv_view[0], pos)
if particle is None:
qpms_scatsysw_scattered_field_basis(res, self.ssw, btyp_c, pos)
else:
res = qpms_scatsysw_scattered_E(self.ssw, btyp_c, &scv_view[0], pos)
results[i,0] = res.x
results[i,1] = res.y
results[i,2] = res.z
return results.reshape(evalpos.shape)
qpms_scatsysw_scattered_field_basis_pi(res, self.ssw, pi, btyp_c, pos)
for j in range(basissize):
results[i,j,0] = res[j].x
results[i,j,1] = res[j].y
results[i,j,2] = res[j].z
free(res)
return results.reshape(evalpos.shape[:-1] + (basissize, 3))
def bspec_pi(self, pi):
return self.ss_pyref.bspec_pi(pi)
property bspecs:
def __get__(self):
return self.ss_pyref.bspecs
cdef class ScatteringMatrix:

View File

@ -4,9 +4,6 @@ ctypedef double complex cdouble
from libc.stdint cimport *
cdef extern from "version.h":
const char *qpms_version_git_sha1
cdef extern from "gsl/gsl_errno.h":
ctypedef void gsl_error_handler_t (const char *reason, const char *file,
int line, int gsl_errno)
@ -150,6 +147,11 @@ cdef extern from "qpms_types.h":
cdouble mu
ctypedef enum qpms_coord_system_t:
pass
ctypedef enum qpms_ewald_part:
QPMS_EWALD_LONG_RANGE
QPMS_EWALD_SHORT_RANGE
QPMS_EWALD_FULL
QPMS_EWALD_0TERM
# maybe more if needed
cdef extern from "qpms_error.h":
@ -188,10 +190,11 @@ cdef extern from "vswf.h":
csphvec_t qpms_vswf_single_mg_csph(qpms_m_t m, qpms_l_t n, csph_t kdlj, qpms_bessel_t btyp, qpms_normalisation_t norm)
cdef extern from "indexing.h":
qpms_y_t qpms_lMax2nelem(qpms_l_t lMax)
qpms_uvswfi_t qpms_tmn2uvswfi(qpms_vswf_type_t t, qpms_m_t m, qpms_l_t n)
qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u, qpms_vswf_type_t* t, qpms_m_t* m, qpms_l_t* n)
qpms_m_t qpms_uvswfi2m(qpms_uvswfi_t u)
qpms_y_t qpms_lMax2nelem(qpms_l_t lMax) nogil
qpms_uvswfi_t qpms_tmn2uvswfi(qpms_vswf_type_t t, qpms_m_t m, qpms_l_t n) nogil
qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u, qpms_vswf_type_t* t, qpms_m_t* m, qpms_l_t* n) nogil
qpms_m_t qpms_uvswfi2m(qpms_uvswfi_t u) nogil
qpms_y_t qpms_lMax2nelem_sc(qpms_l_t lmax) nogil
# maybe more if needed
# Point generators from lattices.h
@ -254,95 +257,95 @@ cdef extern from "lattices.h":
cdef extern from "vectors.h":
cart2_t cart2_add(const cart2_t a, const cart2_t b)
cart2_t cart2_substract(const cart2_t a, const cart2_t b)
cart2_t cart2_scale(const double c, const cart2_t v)
double cart2_dot(const cart2_t a, const cart2_t b)
double cart2_normsq(const cart2_t a)
double cart2norm(const cart2_t v)
pol_t cart2pol(const cart2_t cart)
sph_t pol2sph_equator(const pol_t pol)
sph_t cart22sph(const cart2_t cart)
sph_t cart12sph_zaxis(double z)
cart3_t cart12cart3z(double z)
cart3_t cart22cart3xy(const cart2_t a)
cart2_t cart3xy2cart2(const cart3_t a)
double cart3_dot(const cart3_t a, const cart3_t b)
cart3_t cart3_vectorprod(const cart3_t a, const cart3_t b)
double cart3_tripleprod(const cart3_t a, const cart3_t b, const cart3_t c)
double cart3_normsq(const cart3_t a)
double cart3norm(const cart3_t v)
sph_t cart2sph(const cart3_t cart)
cart3_t sph2cart(const sph_t sph)
cart2_t pol2cart(const pol_t pol)
cart3_t pol2cart3_equator(const pol_t pol)
cart3_t cart3_add(const cart3_t a, const cart3_t b)
cart3_t cart3_substract(const cart3_t a, const cart3_t b)
cart3_t cart3_scale(const double c, const cart3_t v)
cart3_t cart3_divscale( const cart3_t v, const double c)
double cart3_dist(const cart3_t a, const cart3_t b)
bint cart3_isclose(const cart3_t a, const cart3_t b, double rtol, double atol)
ccart3_t ccart3_scale(const cdouble c, const ccart3_t v)
ccart3_t ccart3_add(const ccart3_t a, const ccart3_t b)
ccart3_t ccart3_substract(const ccart3_t a, const ccart3_t b)
cdouble ccart3_dotnc(const ccart3_t a, const ccart3_t b)
ccart3_t cart32ccart3(cart3_t c)
csphvec_t csphvec_add(const csphvec_t a, const csphvec_t b)
csphvec_t csphvec_substract(const csphvec_t a, const csphvec_t b)
csphvec_t csphvec_scale(cdouble c, const csphvec_t v)
cdouble csphvec_dotnc(const csphvec_t a, const csphvec_t b)
sph_t sph_scale(double c, const sph_t s)
csph_t sph_cscale(cdouble c, const sph_t s)
ccart3_t csphvec2ccart(const csphvec_t sphvec, const sph_t at)
ccart3_t csphvec2ccart_csph(const csphvec_t sphvec, const csph_t at)
csphvec_t ccart2csphvec(const ccart3_t cartvec, const sph_t at)
csph_t sph2csph(sph_t s)
sph_t csph2sph(csph_t s)
csph_t ccart2csph(const ccart3_t cart)
csph_t cart2csph(const cart3_t cart)
ccart3_t csph2ccart(const csph_t sph)
void ccart3_kahaninit(ccart3_t *sum, ccart3_t *compensation)
void csphvec_kahaninit(csphvec_t *sum, csphvec_t *compensation)
void ccart3_kahanadd(ccart3_t *sum, ccart3_t *compensation, const ccart3_t input)
void csphvec_kahanadd(csphvec_t *sum, csphvec_t *compensation, const csphvec_t input)
double csphvec_norm(const csphvec_t a)
double csphvec_reldiff_abstol(const csphvec_t a, const csphvec_t b, double tolerance)
double csphvec_reldiff(const csphvec_t a, const csphvec_t b)
sph_t anycoord2sph(anycoord_point_t p, qpms_coord_system_t t)
cart3_t anycoord2cart3(anycoord_point_t p, qpms_coord_system_t t)
double anycoord_norm(anycoord_point_t p, qpms_coord_system_t t)
pol_t anycoord2pol(anycoord_point_t p, qpms_coord_system_t t)
cart2_t anycoord2cart2(anycoord_point_t p, qpms_coord_system_t t)
double anycoord2cart1(anycoord_point_t p, qpms_coord_system_t t)
cart2_t cart2_add(const cart2_t a, const cart2_t b) nogil
cart2_t cart2_substract(const cart2_t a, const cart2_t b) nogil
cart2_t cart2_scale(const double c, const cart2_t v) nogil
double cart2_dot(const cart2_t a, const cart2_t b) nogil
double cart2_normsq(const cart2_t a) nogil
double cart2norm(const cart2_t v) nogil
pol_t cart2pol(const cart2_t cart) nogil
sph_t pol2sph_equator(const pol_t pol) nogil
sph_t cart22sph(const cart2_t cart) nogil
sph_t cart12sph_zaxis(double z) nogil
cart3_t cart12cart3z(double z) nogil
cart3_t cart22cart3xy(const cart2_t a) nogil
cart2_t cart3xy2cart2(const cart3_t a) nogil
double cart3_dot(const cart3_t a, const cart3_t b) nogil
cart3_t cart3_vectorprod(const cart3_t a, const cart3_t b) nogil
double cart3_tripleprod(const cart3_t a, const cart3_t b, const cart3_t c) nogil
double cart3_normsq(const cart3_t a) nogil
double cart3norm(const cart3_t v) nogil
sph_t cart2sph(const cart3_t cart) nogil
cart3_t sph2cart(const sph_t sph) nogil
cart2_t pol2cart(const pol_t pol) nogil
cart3_t pol2cart3_equator(const pol_t pol) nogil
cart3_t cart3_add(const cart3_t a, const cart3_t b) nogil
cart3_t cart3_substract(const cart3_t a, const cart3_t b) nogil
cart3_t cart3_scale(const double c, const cart3_t v) nogil
cart3_t cart3_divscale( const cart3_t v, const double c) nogil
double cart3_dist(const cart3_t a, const cart3_t b) nogil
bint cart3_isclose(const cart3_t a, const cart3_t b, double rtol, double atol) nogil
ccart3_t ccart3_scale(const cdouble c, const ccart3_t v) nogil
ccart3_t ccart3_add(const ccart3_t a, const ccart3_t b) nogil
ccart3_t ccart3_substract(const ccart3_t a, const ccart3_t b) nogil
cdouble ccart3_dotnc(const ccart3_t a, const ccart3_t b) nogil
ccart3_t cart32ccart3(cart3_t c) nogil
csphvec_t csphvec_add(const csphvec_t a, const csphvec_t b) nogil
csphvec_t csphvec_substract(const csphvec_t a, const csphvec_t b) nogil
csphvec_t csphvec_scale(cdouble c, const csphvec_t v) nogil
cdouble csphvec_dotnc(const csphvec_t a, const csphvec_t b) nogil
sph_t sph_scale(double c, const sph_t s) nogil
csph_t sph_cscale(cdouble c, const sph_t s) nogil
ccart3_t csphvec2ccart(const csphvec_t sphvec, const sph_t at) nogil
ccart3_t csphvec2ccart_csph(const csphvec_t sphvec, const csph_t at) nogil
csphvec_t ccart2csphvec(const ccart3_t cartvec, const sph_t at) nogil
csph_t sph2csph(sph_t s) nogil
sph_t csph2sph(csph_t s) nogil
csph_t ccart2csph(const ccart3_t cart) nogil
csph_t cart2csph(const cart3_t cart) nogil
ccart3_t csph2ccart(const csph_t sph) nogil
void ccart3_kahaninit(ccart3_t *sum, ccart3_t *compensation) nogil
void csphvec_kahaninit(csphvec_t *sum, csphvec_t *compensation) nogil
void ccart3_kahanadd(ccart3_t *sum, ccart3_t *compensation, const ccart3_t input) nogil
void csphvec_kahanadd(csphvec_t *sum, csphvec_t *compensation, const csphvec_t input) nogil
double csphvec_norm(const csphvec_t a) nogil
double csphvec_reldiff_abstol(const csphvec_t a, const csphvec_t b, double tolerance) nogil
double csphvec_reldiff(const csphvec_t a, const csphvec_t b) nogil
sph_t anycoord2sph(anycoord_point_t p, qpms_coord_system_t t) nogil
cart3_t anycoord2cart3(anycoord_point_t p, qpms_coord_system_t t) nogil
double anycoord_norm(anycoord_point_t p, qpms_coord_system_t t) nogil
pol_t anycoord2pol(anycoord_point_t p, qpms_coord_system_t t) nogil
cart2_t anycoord2cart2(anycoord_point_t p, qpms_coord_system_t t) nogil
double anycoord2cart1(anycoord_point_t p, qpms_coord_system_t t) nogil
void qpms_array_coord_transform(void *dest, qpms_coord_system_t tdest,
const void *src, qpms_coord_system_t tsrc, size_t nmemb)
const void *src, qpms_coord_system_t tsrc, size_t nmemb) nogil
void anycoord_arr2something(void *dest, qpms_coord_system_t tdest,
const void *src, qpms_coord_system_t tsrc, size_t nmemb)
void cart3_to_double_array(double *a, cart3_t b)
cart3_t cart3_from_double_array(const double *a)
void cart2_to_double_array(double *a, cart2_t b)
cart2_t cart2_from_double_array(const double *a)
const void *src, qpms_coord_system_t tsrc, size_t nmemb) nogil
void cart3_to_double_array(double *a, cart3_t b) nogil
cart3_t cart3_from_double_array(const double *a) nogil
void cart2_to_double_array(double *a, cart2_t b) nogil
cart2_t cart2_from_double_array(const double *a) nogil
const cart2_t CART2_ZERO
cdef extern from "quaternions.h":
qpms_quat_t qpms_quat_2c_from_4d(qpms_quat4d_t q)
qpms_quat4d_t qpms_quat_4d_from_2c(qpms_quat_t q)
qpms_quat_t qpms_quat_mult(qpms_quat_t p, qpms_quat_t q)
qpms_quat_t qpms_quat_add(qpms_quat_t p, qpms_quat_t q)
qpms_quat_t qpms_quat_rscale(double s, qpms_quat_t q)
qpms_quat_t qpms_quat_conj(qpms_quat_t q)
double qpms_quat_norm(qpms_quat_t q)
double qpms_quat_imnorm(qpms_quat_t q)
qpms_quat_t qpms_quat_normalise(qpms_quat_t q)
qpms_quat_t qpms_quat_log(qpms_quat_t q)
qpms_quat_t qpms_quat_exp(qpms_quat_t q)
qpms_quat_t qpms_quat_pow(qpms_quat_t q, double exponent)
qpms_quat_t qpms_quat_2c_from_4d(qpms_quat4d_t q) nogil
qpms_quat4d_t qpms_quat_4d_from_2c(qpms_quat_t q) nogil
qpms_quat_t qpms_quat_mult(qpms_quat_t p, qpms_quat_t q) nogil
qpms_quat_t qpms_quat_add(qpms_quat_t p, qpms_quat_t q) nogil
qpms_quat_t qpms_quat_rscale(double s, qpms_quat_t q) nogil
qpms_quat_t qpms_quat_conj(qpms_quat_t q) nogil
double qpms_quat_norm(qpms_quat_t q) nogil
double qpms_quat_imnorm(qpms_quat_t q) nogil
qpms_quat_t qpms_quat_normalise(qpms_quat_t q) nogil
qpms_quat_t qpms_quat_log(qpms_quat_t q) nogil
qpms_quat_t qpms_quat_exp(qpms_quat_t q) nogil
qpms_quat_t qpms_quat_pow(qpms_quat_t q, double exponent) nogil
cdouble qpms_wignerD_elem(qpms_quat_t q, qpms_l_t l,
qpms_m_t mp, qpms_m_t m)
qpms_irot3_t qpms_irot3_mult(qpms_irot3_t p, qpms_irot3_t q)
qpms_irot3_t qpms_irot3_pow(qpms_irot3_t p, int n)
qpms_quat_t qpms_quat_from_rotvector(cart3_t v)
qpms_m_t mp, qpms_m_t m) nogil
qpms_irot3_t qpms_irot3_mult(qpms_irot3_t p, qpms_irot3_t q) nogil
qpms_irot3_t qpms_irot3_pow(qpms_irot3_t p, int n) nogil
qpms_quat_t qpms_quat_from_rotvector(cart3_t v) nogil
cdef extern from "groups.h":
struct qpms_finite_group_irrep_t:
@ -490,6 +493,7 @@ cdef extern from "tmatrices.h":
qpms_errno_t qpms_tmatrix_generator_interpolator(qpms_tmatrix_t *t, cdouble omega, const void *params)
qpms_errno_t qpms_tmatrix_generator_sphere(qpms_tmatrix_t *t, cdouble omega, const void *params)
qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t, cdouble omega, const void *params)
qpms_errno_t qpms_tmatrix_generator_zero(qpms_tmatrix_t *t, cdouble omega, const void *params)
struct qpms_tmatrix_generator_sphere_param_t:
qpms_epsmu_generator_t outside
qpms_epsmu_generator_t inside
@ -623,6 +627,7 @@ cdef extern from "scatsystem.h":
size_t *saecv_sizes
const qpms_finite_group_t *sym
qpms_scatsys_periodic_info_t per
qpms_trans_calculator *c
# per[] and other stuff not currently needed in cython
void qpms_scatsys_free(qpms_scatsys_t *s)
@ -638,6 +643,11 @@ cdef extern from "scatsystem.h":
cdouble omega, const qpms_tolerance_spec_t *tol)
qpms_scatsys_at_omega_t *qpms_scatsys_at_omega(const qpms_scatsys_t *ss, cdouble omega)
void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *ssw)
cdouble *qpms_scatsys_irrep_transform_matrix(cdouble *target, const qpms_scatsys_t *ss, qpms_iri_t iri)
cdouble *qpms_scatsys_irrep_pack_matrix_stupid(cdouble *target_packed,
const cdouble *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri)
cdouble *qpms_scatsys_irrep_unpack_matrix_stupid(cdouble *target_full,
const cdouble *orig_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, bint add)
cdouble *qpms_scatsys_irrep_pack_matrix(cdouble *target_packed,
const cdouble *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri)
cdouble *qpms_scatsys_irrep_unpack_matrix(cdouble *target_full,
@ -708,8 +718,20 @@ cdef extern from "scatsystem.h":
const cdouble *f_excitation_vector_full, cart3_t where) nogil
ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk, qpms_bessel_t btyp,
const cdouble *f_excitation_vector_full, cart3_t where) nogil
ccart3_t qpms_scatsyswk_scattered_E_e(const qpms_scatsys_at_omega_k_t *sswk, qpms_bessel_t btyp,
const cdouble *f_excitation_vector_full, cart3_t where, qpms_ewald_part parts) nogil
qpms_errno_t qpms_scatsys_scattered_field_basis(ccart3_t *target, const qpms_scatsys_t *ss,
qpms_bessel_t btyp, cdouble wavenumber, cart3_t where) nogil
qpms_errno_t qpms_scatsysw_scattered_field_basis(ccart3_t *target, const qpms_scatsys_at_omega_t *ssw,
qpms_bessel_t btyp, cart3_t where) nogil
qpms_errno_t qpms_scatsys_scattered_field_basis_pi(ccart3_t *target, const qpms_scatsys_t *ss,
qpms_ss_pi_t pi, qpms_bessel_t btyp, cdouble wavenumber, cart3_t where) nogil
qpms_errno_t qpms_scatsysw_scattered_field_basis_pi(ccart3_t *target, const qpms_scatsys_at_omega_t *ssw,
qpms_ss_pi_t pi, qpms_bessel_t btyp, cart3_t where) nogil
qpms_errno_t qpms_scatsyswk_scattered_field_basis(ccart3_t *target, const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t btyp, cart3_t where) nogil
qpms_errno_t qpms_scatsyswk_scattered_field_basis_e(ccart3_t *target, const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t btyp, cart3_t where, qpms_ewald_part parts) nogil
double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, cdouble wavenumber, const double *wavevector) nogil
@ -718,12 +740,6 @@ cdef extern from "ewald.h":
cdouble val
double err
ctypedef enum qpms_ewald_part:
QPMS_EWALD_LONG_RANGE
QPMS_EWALD_SHORT_RANGE
QPMS_EWALD_FULL
QPMS_EWALD_0TERM
struct qpms_ewald3_constants_t:
qpms_l_t lMax
qpms_y_t nelem_sc

View File

@ -5,6 +5,9 @@
*/
#ifndef QPMS_ERROR_H
#define QPMS_ERROR_H
#ifdef __cplusplus
extern "C" {
#endif
#include "optim.h"
@ -266,4 +269,7 @@ qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types);
}\
}
#ifdef __cplusplus
}
#endif
#endif

View File

@ -6,7 +6,7 @@ import numpy as np
from .qpms_c import *
ň = np.newaxis
import scipy
from .constants import ε_0, c, pi, π, e, , μ_0
from scipy.constants import epsilon_0 as ε_0, c, pi as π, e, hbar as , mu_0 as μ_0
eV = e
from scipy.special import lpmn, lpmv, spherical_jn, spherical_yn, poch, gammaln, factorial
import math

View File

@ -3,6 +3,10 @@
*/
#ifndef QPMS_SPECFUNC_H
#define QPMS_SPECFUNC_H
#ifdef __cplusplus
extern "C" {
#endif
#include "qpms_types.h"
#include <gsl/gsl_sf_legendre.h>
@ -11,25 +15,25 @@
******************************************************************************/
// TODO unify types
qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, complex double x, complex double *result_array);
qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, _Complex double x, _Complex double *result_array);
typedef struct {
qpms_l_t lMax;
double *akn; // coefficients as in DLMF 10.49.1
//complex double *bkn; // coefficients of the derivatives
//_Complex double *bkn; // coefficients of the derivatives
} qpms_sbessel_calculator_t;
qpms_sbessel_calculator_t *qpms_sbessel_calculator_init(void);
void qpms_sbessel_calculator_pfree(qpms_sbessel_calculator_t *c);
qpms_errno_t qpms_sbessel_calc_fill(qpms_sbessel_calculator_t *c, qpms_bessel_t typ, qpms_l_t lmax,
double x, complex double *result_array);
double x, _Complex double *result_array);
complex double qpms_sbessel_calc_h1(qpms_sbessel_calculator_t *c, qpms_l_t n, complex double x);
_Complex double qpms_sbessel_calc_h1(qpms_sbessel_calculator_t *c, qpms_l_t n, _Complex double x);
qpms_errno_t qpms_sbessel_calc_h1_fill(qpms_sbessel_calculator_t *c, qpms_l_t lmax,
complex double x, complex double *result_array);
_Complex double x, _Complex double *result_array);
/******************************************************************************
@ -107,4 +111,7 @@ double qpms_legendre0(int m, int n);
// Associated Legendre polynomial derivative at zero argument (DLMF 14.5.2)
double qpms_legendred0(int m, int n);
#ifdef __cplusplus
}
#endif
#endif // QPMS_SPECFUNC_H

View File

@ -3,10 +3,13 @@
*/
#ifndef QPMS_TYPES_H
#define QPMS_TYPES_H
#include <complex.h>
#include <stdbool.h>
#ifdef __cplusplus
extern "C" {
#endif
#include <stddef.h>
#include <stdint.h>
#include <stdbool.h>
#ifndef M_PI_2
#define M_PI_2 (1.570796326794896619231321691639751442098584699687552910487L)
@ -99,7 +102,9 @@ typedef enum {
* These bit flags are used by the functions declared in normalisation.h
* that return the appropriate convention-dependent factors.
*
* See @ref vswf_conventions for comparison of the various conventions used.
* \see @ref vswf_conventions for comparison of the various conventions used.
* \see @ref swf_conventions_qpms for description how the conventions are used internally and in the QPMS API.
*
*/
typedef enum {
QPMS_NORMALISATION_UNDEF = 0, ///< Convention undefined. This should not happen.
@ -189,7 +194,7 @@ typedef struct cart3_t {
/// 3D complex (actually 6D) coordinates. See also vectors.h.
typedef struct ccart3_t {
complex double x, y, z;
_Complex double x, y, z;
} ccart3_t;
/// 3D complex vector pair (represents the E, H fields).
@ -210,13 +215,13 @@ typedef struct sph_t {
/// Spherical coordinates with complex radial component. See also vectors.h.
typedef struct csph_t { // Do I really need this???
complex double r;
_Complex double r;
double theta, phi;
} csph_t;
/// 3D complex vector components in local spherical basis. See also vectors.h.
typedef struct csphvec_t {
complex double rc, thetac, phic;
_Complex double rc, thetac, phic;
} csphvec_t;
/// 2D polar coordinates. See also vectors.h.
@ -259,7 +264,7 @@ typedef enum {
* See quaternions.h for "methods".
*/
typedef struct qpms_quat_t {
complex double a, b;
_Complex double a, b;
} qpms_quat_t;
/// Quaternion type as four doubles.
@ -344,7 +349,7 @@ typedef struct qpms_tmatrix_t {
* there is also one with order \a m.
*/
const qpms_vswf_set_spec_t *spec;
complex double *m; ///< Matrix elements in row-major order.
_Complex double *m; ///< Matrix elements in row-major order.
bool owns_m; ///< Information wheter m shall be deallocated with qpms_tmatrix_free()
} qpms_tmatrix_t;
@ -420,11 +425,22 @@ typedef struct qpms_abstract_tmatrix_t {
/// A type holding electric permittivity and magnetic permeability of a material.
typedef struct qpms_epsmu_t {
complex double eps; ///< Relative permittivity.
complex double mu; ///< Relative permeability.
_Complex double eps; ///< Relative permittivity.
_Complex double mu; ///< Relative permeability.
} qpms_epsmu_t;
struct qpms_tolerance_spec_t; // See tolerances.h
typedef enum {
QPMS_EWALD_LONG_RANGE = 1,
QPMS_EWALD_SHORT_RANGE = 2,
QPMS_EWALD_0TERM = 4,
QPMS_EWALD_FULL = QPMS_EWALD_LONG_RANGE | QPMS_EWALD_SHORT_RANGE | QPMS_EWALD_0TERM,
} qpms_ewald_part;
#define lmcheck(l,m) assert((l) >= 1 && abs(m) <= (l))
#ifdef __cplusplus
}
#endif
#endif // QPMS_TYPES

View File

@ -14,6 +14,10 @@
*/
#ifndef QPMSBLAS_H
#define QPMSBLAS_H
#ifdef __cplusplus
extern "C" {
#endif
#define QPMS_BLAS_INDEX_T long long int
#ifndef CBLAS_H
@ -31,4 +35,7 @@ void qpms_zgemm(CBLAS_LAYOUT Order, CBLAS_TRANSPOSE TransA, CBLAS_TRANSPOSE Tran
const _Complex double *B, const QPMS_BLAS_INDEX_T ldb,
const _Complex double *beta, _Complex double *C, const QPMS_BLAS_INDEX_T ldc);
#ifdef __cplusplus
}
#endif
#endif //QPMSBLAS_H

View File

@ -3,6 +3,9 @@
*/
#ifndef QPMS_WIGNER_H
#define QPMS_WIGNER_H
#ifdef __cplusplus
extern "C" {
#endif
#include "qpms_types.h"
#include "vectors.h"
@ -93,7 +96,7 @@ static inline double qpms_quat_norm(const qpms_quat_t q) {
}
/// Test approximate equality of quaternions.
static inline bool qpms_quat_isclose(const qpms_quat_t p, const qpms_quat_t q, double atol) {
static inline _Bool qpms_quat_isclose(const qpms_quat_t p, const qpms_quat_t q, double atol) {
return qpms_quat_norm(qpms_quat_sub(p,q)) <= atol;
}
@ -119,7 +122,7 @@ static inline qpms_quat_t qpms_quat_standardise(qpms_quat_t p, double atol) {
}
/// Test approximate equality of "standardised" quaternions, so that \f$-q\f$ is considered equal to \f$q\f$.
static inline bool qpms_quat_isclose2(const qpms_quat_t p, const qpms_quat_t q, double atol) {
static inline _Bool qpms_quat_isclose2(const qpms_quat_t p, const qpms_quat_t q, double atol) {
return qpms_quat_norm(qpms_quat_sub(
qpms_quat_standardise(p, atol),
qpms_quat_standardise(q, atol))) <= atol;
@ -202,17 +205,17 @@ static inline qpms_quat_t qpms_quat_from_rotvector(cart3_t v) {
* The D matrix are calculated using formulae (3), (4), (6), (7) from
* http://moble.github.io/spherical_functions/WignerDMatrices.html
*/
complex double qpms_wignerD_elem(qpms_quat_t q, qpms_l_t l,
_Complex double qpms_wignerD_elem(qpms_quat_t q, qpms_l_t l,
qpms_m_t mp, qpms_m_t m);
/// A VSWF representation element of the O(3) group.
/**
* TODO more doc.
*/
complex double qpms_vswf_irot_elem_from_irot3(
_Complex double qpms_vswf_irot_elem_from_irot3(
const qpms_irot3_t q, ///< The O(3) element in the quaternion representation.
qpms_l_t l, qpms_m_t mp, qpms_m_t m,
bool pseudo ///< Determines the sign of improper rotations. True for magnetic waves, false otherwise.
_Bool pseudo ///< Determines the sign of improper rotations. True for magnetic waves, false otherwise.
);
@ -251,7 +254,7 @@ static inline qpms_irot3_t qpms_irot3_pow(const qpms_irot3_t p, int n) {
}
/// Test approximate equality of irot3.
static inline bool qpms_irot3_isclose(const qpms_irot3_t p, const qpms_irot3_t q, double atol) {
static inline _Bool qpms_irot3_isclose(const qpms_irot3_t p, const qpms_irot3_t q, double atol) {
return qpms_quat_isclose2(p.rot, q.rot, atol) && p.det == q.det;
}
@ -303,4 +306,7 @@ static inline qpms_irot3_t qpms_irot3_zrot_Nk(double N, double k) {
return qpms_irot3_zrot_angle(2 * M_PI * k / N);
}
#ifdef __cplusplus
}
#endif
#endif //QPMS_WIGNER_H

View File

@ -24,6 +24,7 @@
#include "tolerances.h"
#include "beyn.h"
#include "tiny_inlines.h"
#include "lattices.h"
#ifdef QPMS_SCATSYSTEM_USE_OWN_BLAS
#include "qpmsblas.h"
@ -2049,6 +2050,61 @@ ccart3_t qpms_scatsysw_scattered_E(const qpms_scatsys_at_omega_t *ssw,
cvf, where);
}
qpms_errno_t qpms_scatsys_scattered_field_basis_pi(
ccart3_t *target,
const qpms_scatsys_t *ss,
const qpms_ss_pi_t pi,
const qpms_bessel_t btyp,
const complex double k,
const cart3_t where
) {
qpms_ss_ensure_nonperiodic_a(ss, "qpms_scatsyswk_scattered_field_basis()");
QPMS_UNTESTED;
const qpms_vswf_set_spec_t *bspec = qpms_ss_bspec_pi(ss, pi);
csphvec_t *vswfs_sph; //Single particle contributions in spherical coordinates
QPMS_CRASHING_MALLOC(vswfs_sph, bspec->n * sizeof(*vswfs_sph));
const cart3_t particle_pos = ss->p[pi].pos;
const csph_t kr = sph_cscale(k, cart2sph(
cart3_substract(where, particle_pos)));
QPMS_ENSURE_SUCCESS(qpms_uvswf_fill(vswfs_sph, bspec, kr, btyp));
for(size_t i = 0; i < bspec->n; ++i)
target[i] = csphvec2ccart_csph(vswfs_sph[i], kr);
free(vswfs_sph);
return QPMS_SUCCESS;
}
qpms_errno_t qpms_scatsysw_scattered_field_basis_pi(
ccart3_t *target, const qpms_scatsys_at_omega_t *ssw, const qpms_ss_pi_t pi,
const qpms_bessel_t btyp, const cart3_t where) {
return qpms_scatsys_scattered_field_basis_pi(target, ssw->ss, pi, btyp,
ssw->wavenumber, where);
}
qpms_errno_t qpms_scatsys_scattered_field_basis(
ccart3_t *target,
const qpms_scatsys_t *ss,
const qpms_bessel_t btyp,
const complex double k,
const cart3_t where
) {
qpms_ss_ensure_nonperiodic_a(ss, "qpms_scatsyswk_scattered_field_basis()");
QPMS_UNTESTED;
for (qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi)
QPMS_ENSURE_SUCCESS(qpms_scatsys_scattered_field_basis_pi(
target + ss->fecv_pstarts[pi], ss, pi, btyp, k, where));
return QPMS_SUCCESS;
}
qpms_errno_t qpms_scatsysw_scattered_field_basis(
ccart3_t *target, const qpms_scatsys_at_omega_t *ssw,
const qpms_bessel_t btyp, const cart3_t where) {
return qpms_scatsys_scattered_field_basis(target, ssw->ss, btyp,
ssw->wavenumber, where);
}
const qpms_uvswfi_t ELDIPILIST[] = {6, 10, 14};
#define DIPSPECN 3 // We have three basis vectors
// Evaluates the regular electric dipole waves in the origin. The returned
@ -2058,11 +2114,15 @@ static inline const qpms_vswf_set_spec_t qpms_fill_regdipoles_0(
// bspec containing only electric dipoles
const qpms_vswf_set_spec_t dipspec = {
.n = DIPSPECN,
#if 0
.ilist = (qpms_uvswfi_t[]){
qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, -1, 1),
qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, 0, 1),
qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, +1, 1),
},
#else
.ilist = ELDIPILIST,
#endif
.lMax=1, .lMax_M=0, .lMax_N=1, .lMax_L=-1,
.capacity=0,
.norm = normalisation,
@ -2127,10 +2187,11 @@ ccart3_t qpms_scatsysw_scattered_E__alt(const qpms_scatsys_at_omega_t *ssw,
// For periodic lattices, we use directly the "alternative" implementation,
// using translation operator and regular dipole waves at zero
ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk,
ccart3_t qpms_scatsyswk_scattered_E_e(const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t btyp,
const complex double *cvf,
const cart3_t where
const cart3_t where,
const qpms_ewald_part parts
) {
QPMS_UNTESTED;
if (btyp != QPMS_HANKEL_PLUS)
@ -2162,13 +2223,13 @@ ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk,
const double maxR = sqrt(ss->per.unitcell_volume) * 64;
const double maxK = 2048 * 2 * M_PI / maxR;
QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32(
QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32_e(
ss->c, s, NULL,
&dipspec, 1, bspec, dipspec.n,
sswk->eta, sswk->ssw->wavenumber,
cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]),
cart2_from_double_array(sswk->k), cart3_substract(where, particle_pos) /*CHECKSIGN*/,
maxR, maxK));
maxR, maxK, parts));
for(size_t i = 0; i < bspec->n; ++i)
for(size_t j = 0; j < dipspec.n; ++j){
@ -2182,11 +2243,17 @@ ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk,
return res;
}
qpms_errno_t qpms_scatsyswk_scattered_field_basis(
ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t btyp, const complex double *cvf, const cart3_t where) {
return qpms_scatsyswk_scattered_E_e(sswk, btyp, cvf, where, QPMS_EWALD_FULL);
}
qpms_errno_t qpms_scatsyswk_scattered_field_basis_e(
ccart3_t *target,
const qpms_scatsys_at_omega_k_t *sswk,
const qpms_bessel_t btyp,
const cart3_t where
const cart3_t where,
const qpms_ewald_part parts
) {
QPMS_UNTESTED;
if (btyp != QPMS_HANKEL_PLUS)
@ -2220,13 +2287,13 @@ qpms_errno_t qpms_scatsyswk_scattered_field_basis(
const double maxR = sqrt(ss->per.unitcell_volume) * 64;
const double maxK = 2048 * 2 * M_PI / maxR;
QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32(
QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32_e(
ss->c, s, NULL,
&dipspec, 1, bspec, dipspec.n,
sswk->eta, sswk->ssw->wavenumber,
cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]),
cart2_from_double_array(sswk->k), cart3_substract(where, particle_pos) /*CHECKSIGN*/,
maxR, maxK));
maxR, maxK, parts));
for(size_t i = 0; i < bspec->n; ++i)
for(size_t j = 0; j < dipspec.n; ++j){
@ -2242,6 +2309,13 @@ qpms_errno_t qpms_scatsyswk_scattered_field_basis(
return QPMS_SUCCESS;
}
qpms_errno_t qpms_scatsyswk_scattered_field_basis(
ccart3_t *target, const qpms_scatsys_at_omega_k_t *sswk,
const qpms_bessel_t btyp, const cart3_t where) {
return qpms_scatsyswk_scattered_field_basis_e(
target, sswk, btyp, where, QPMS_EWALD_FULL);
}
#if 0
ccart3_t qpms_scatsys_scattered_E_irrep(const qpms_scatsys_t *ss,
qpms_iri_t iri, const complex double *cvr, cart3_t where) {

View File

@ -10,6 +10,10 @@
*/
#ifndef QPMS_SCATSYSTEM_H
#define QPMS_SCATSYSTEM_H
#ifdef __cplusplus
extern "C" {
#endif
#include "qpms_types.h"
#include "vswf.h"
#include "tmatrices.h"
@ -107,7 +111,7 @@ typedef struct qpms_ss_orbit_type_t {
*
* TODO doc.
*/
complex double *irbases;
_Complex double *irbases;
/// TODO doc.
size_t instance_count;
/// Cumulative sum of the preceding ot->siza * ot->instance_count;
@ -243,9 +247,9 @@ typedef struct qpms_scatsys_at_omega_t {
* i.e in the order corresponding to \a ss->tm.
*/
qpms_tmatrix_t **tm;
complex double omega; ///< Angular frequency
_Complex double omega; ///< Angular frequency
qpms_epsmu_t medium; ///< Background medium optical properties at the given frequency
complex double wavenumber; ///< Background medium wave number
_Complex double wavenumber; ///< Background medium wave number
} qpms_scatsys_at_omega_t;
@ -291,7 +295,7 @@ typedef struct qpms_scatsys_at_omega_t {
*
*/
qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const struct qpms_finite_group_t *sym,
complex double omega, const struct qpms_tolerance_spec_t *tol);
_Complex double omega, const struct qpms_tolerance_spec_t *tol);
/// Destroys the result of qpms_scatsys_apply_symmetry or qpms_scatsys_load.
void qpms_scatsys_free(qpms_scatsys_t *s);
@ -303,50 +307,50 @@ void qpms_scatsys_at_omega_free(qpms_scatsys_at_omega_t *ssw);
/// Evaluates scattering system T-matrices at a given frequency.
/** Free the result using qpms_scatsys_at_omega_free() when done. */
qpms_scatsys_at_omega_t *qpms_scatsys_at_omega(const qpms_scatsys_t *ss,
complex double omega);
_Complex double omega);
/// Creates a "full" transformation matrix U that takes a full vector and projects it onto an symmetry adapted basis.
/** Mostly as a reference and a debugging tool, as multiplicating these big matrices would be inefficient.
*
* TODO doc about shape etc.
*/
complex double *qpms_scatsys_irrep_transform_matrix(complex double *target_U,
_Complex double *qpms_scatsys_irrep_transform_matrix(_Complex double *target_U,
const qpms_scatsys_t *ss, qpms_iri_t iri);
/// Projects a "big" matrix onto an irrep (slow reference implementation).
/** TODO doc */
complex double *qpms_scatsys_irrep_pack_matrix_stupid(complex double *target_packed,
const complex double *orig_full, const qpms_scatsys_t *ss,
_Complex double *qpms_scatsys_irrep_pack_matrix_stupid(_Complex double *target_packed,
const _Complex double *orig_full, const qpms_scatsys_t *ss,
qpms_iri_t iri);
/// Transforms a big "packed" matrix into the full basis (slow reference implementation).
/** TODO doc */
complex double *qpms_scatsys_irrep_unpack_matrix_stupid(complex double *target_full,
const complex double *orig_packed, const qpms_scatsys_t *ss,
_Complex double *qpms_scatsys_irrep_unpack_matrix_stupid(_Complex double *target_full,
const _Complex double *orig_packed, const qpms_scatsys_t *ss,
qpms_iri_t iri, bool add);
/// Projects a "big" matrix onto an irrep.
/** TODO doc */
complex double *qpms_scatsys_irrep_pack_matrix(complex double *target_packed,
const complex double *orig_full, const qpms_scatsys_t *ss,
_Complex double *qpms_scatsys_irrep_pack_matrix(_Complex double *target_packed,
const _Complex double *orig_full, const qpms_scatsys_t *ss,
qpms_iri_t iri);
/// Transforms a big "packed" matrix into the full basis.
/** TODO doc */
complex double *qpms_scatsys_irrep_unpack_matrix(complex double *target_full,
const complex double *orig_packed, const qpms_scatsys_t *ss,
_Complex double *qpms_scatsys_irrep_unpack_matrix(_Complex double *target_full,
const _Complex double *orig_packed, const qpms_scatsys_t *ss,
qpms_iri_t iri, bool add);
/// Projects a "big" vector onto an irrep.
/** TODO doc */
complex double *qpms_scatsys_irrep_pack_vector(complex double *target_packed,
const complex double *orig_full, const qpms_scatsys_t *ss,
_Complex double *qpms_scatsys_irrep_pack_vector(_Complex double *target_packed,
const _Complex double *orig_full, const qpms_scatsys_t *ss,
qpms_iri_t iri);
/// Transforms a big "packed" vector into the full basis.
/** TODO doc */
complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full,
const complex double *orig_packed, const qpms_scatsys_t *ss,
_Complex double *qpms_scatsys_irrep_unpack_vector(_Complex double *target_full,
const _Complex double *orig_packed, const qpms_scatsys_t *ss,
qpms_iri_t iri, bool add);
/// Global translation matrix.
@ -354,31 +358,31 @@ complex double *qpms_scatsys_irrep_unpack_vector(complex double *target_full,
* The diagonal (particle self-) block are filled with zeros (even for regular Bessel waves).
* This may change in the future.
*/
complex double *qpms_scatsys_build_translation_matrix_full(
_Complex double *qpms_scatsys_build_translation_matrix_full(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target,
_Complex double *target,
const qpms_scatsys_t *ss,
complex double k ///< Wave number to use in the translation matrix.
_Complex double k ///< Wave number to use in the translation matrix.
);
/// Creates the full \f$ (I - WS) \f$ matrix of the periodic scattering system.
/**
* \returns \a target on success, NULL on error.
*/
complex double *qpms_scatsysw_build_modeproblem_matrix_full(
_Complex double *qpms_scatsysw_build_modeproblem_matrix_full(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target,
_Complex double *target,
const qpms_scatsys_at_omega_t *ssw
);
/// As qpms_scatsys_build_translation_full() but with choice of Bessel function type.
/** Might be useful for evaluation of cross sections and testing.
*/
complex double *qpms_scatsys_build_translation_matrix_e_full(
_Complex double *qpms_scatsys_build_translation_matrix_e_full(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target,
_Complex double *target,
const qpms_scatsys_t *ss,
complex double k, ///< Wave number to use in the translation matrix.
_Complex double k, ///< Wave number to use in the translation matrix.
qpms_bessel_t J
);
@ -387,12 +391,12 @@ complex double *qpms_scatsys_build_translation_matrix_e_full(
* The diagonal (particle self-) blocks are currently filled with zeros.
* This may change in the future.
*/
complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
_Complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target,
_Complex double *target,
const qpms_scatsys_t *ss,
qpms_iri_t iri,
complex double k, ///< Wave number to use in the translation matrix.
_Complex double k, ///< Wave number to use in the translation matrix.
qpms_bessel_t J
);
@ -400,9 +404,9 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed(
/**
* \returns \a target on success, NULL on error.
*/
complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed(
_Complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target,
_Complex double *target,
const qpms_scatsys_at_omega_t *ssw,
qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym
);
@ -410,9 +414,9 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed(
/**
* \returns \a target on success, NULL on error.
*/
complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR(
_Complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target,
_Complex double *target,
const qpms_scatsys_at_omega_t *ssw,
qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym
);
@ -420,9 +424,9 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR(
/**
* \returns \a target on success, NULL on error.
*/
complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial(
_Complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target,
_Complex double *target,
const qpms_scatsys_at_omega_t *ssw,
qpms_iri_t iri ///< Index of the irreducible representation in ssw->ss->sym
);
@ -435,7 +439,7 @@ typedef struct qpms_ss_LU {
bool full; ///< true if full matrix; false if irrep-packed.
qpms_iri_t iri; ///< Irrep index if `full == false`.
/// LU decomposition array.
complex double *a;
_Complex double *a;
/// Pivot index array, size at least max(1,min(m, n)).
int *ipiv;
} qpms_ss_LU;
@ -443,14 +447,14 @@ void qpms_ss_LU_free(qpms_ss_LU);
/// Builds an LU-factorised mode/scattering problem \f$ (I - TS) \f$ matrix from scratch. Nonperiodic systems only.
qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_full_LU(
complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
_Complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
const qpms_scatsys_at_omega_t *ssw
);
/// Builds an irrep-packed LU-factorised mode/scattering problem matrix from scratch.
qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU(
complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
_Complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
const qpms_scatsys_at_omega_t *ssw,
qpms_iri_t iri
@ -458,7 +462,7 @@ qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU(
/// Computes LU factorisation of a pre-calculated mode/scattering problem matrix, replacing its contents.
qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise(
complex double *modeproblem_matrix_full, ///< Pre-calculated mode problem matrix (I-TS). Mandatory.
_Complex double *modeproblem_matrix_full, ///< Pre-calculated mode problem matrix (I-TS). Mandatory.
int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
const qpms_scatsys_at_omega_t *ssw, ///< Must be filled for non-periodic systems.
const struct qpms_scatsys_at_omega_k_t *sswk ///< Must be filled for periodic systems, otherwise must be NULL.
@ -466,16 +470,16 @@ qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise(
/// Computes LU factorisation of a pre-calculated irrep-packed mode/scattering problem matrix, replacing its contents.
qpms_ss_LU qpms_scatsysw_modeproblem_matrix_irrep_packed_factorise(
complex double *modeproblem_matrix_irrep_packed, ///< Pre-calculated mode problem matrix (I-TS). Mandatory.
_Complex double *modeproblem_matrix_irrep_packed, ///< Pre-calculated mode problem matrix (I-TS). Mandatory.
int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
const qpms_scatsys_at_omega_t *ssw,
qpms_iri_t iri
);
/// Solves a (possibly partial, irrep-packed) scattering problem \f$ (I-TS)f = Ta_\mathrm{inc} \f$ using a pre-factorised \f$ (I-TS) \f$.
complex double *qpms_scatsys_scatter_solve(
complex double *target_f, ///< Target (full or irrep-packed, depending on `ludata.full`) array for \a f. If NULL, a new one is allocated.
const complex double *a_inc, ///< Incident field expansion coefficient vector \a a (full or irrep-packed, depending on `ludata.full`).
_Complex double *qpms_scatsys_scatter_solve(
_Complex double *target_f, ///< Target (full or irrep-packed, depending on `ludata.full`) array for \a f. If NULL, a new one is allocated.
const _Complex double *a_inc, ///< Incident field expansion coefficient vector \a a (full or irrep-packed, depending on `ludata.full`).
qpms_ss_LU ludata ///< Pre-factorised \f$ I - TS \f$ matrix data.
);
@ -495,33 +499,33 @@ typedef struct qpms_scatsys_at_omega_k_t {
/**
* \returns \a target on success, NULL on error.
*/
complex double *qpms_scatsyswk_build_modeproblem_matrix_full(
_Complex double *qpms_scatsyswk_build_modeproblem_matrix_full(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target,
_Complex double *target,
const qpms_scatsys_at_omega_k_t *sswk
);
/// Global translation matrix.
complex double *qpms_scatsys_periodic_build_translation_matrix_full(
_Complex double *qpms_scatsys_periodic_build_translation_matrix_full(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target,
_Complex double *target,
const qpms_scatsys_t *ss,
complex double wavenumber, ///< Wave number to use in the translation matrix.
_Complex double wavenumber, ///< Wave number to use in the translation matrix.
const cart3_t *wavevector, ///< Wavevector / pseudomomentum in cartesian coordinates.
double eta ///< Ewald parameter eta. Pass 0 or NaN to use the default value in \a ss.
);
/// Global translation matrix.
complex double *qpms_scatsyswk_build_translation_matrix_full(
_Complex double *qpms_scatsyswk_build_translation_matrix_full(
/// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
complex double *target,
_Complex double *target,
const qpms_scatsys_at_omega_k_t *sswk
);
/// Builds an LU-factorised mode/scattering problem \f$ (I - TS) \f$ matrix from scratch. Periodic systems only.
qpms_ss_LU qpms_scatsyswk_build_modeproblem_matrix_full_LU(
complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
_Complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated).
int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated).
const qpms_scatsys_at_omega_k_t *sswk
);
@ -539,7 +543,7 @@ struct beyn_result_t *qpms_scatsys_periodic_find_eigenmodes(
const qpms_scatsys_t *ss,
/// Wavevector in cartesian coordinates (must lie in the lattice plane).
const double k[3],
complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
_Complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
double omega_rr, ///< Real half-axis of the ellipse inside which the eigenfreqs are searched for.
double omega_ri, ///< Imaginary half-axis of the ellipse inside which the eigenfreqs are searched for.
size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number)
@ -561,12 +565,12 @@ struct qpms_finite_group_t;
/// Constructs a "full matrix action" of a point group element for an orbit type.
/** TODO detailed doc */
complex double *qpms_orbit_action_matrix(
_Complex double *qpms_orbit_action_matrix(
/// Target array. If NULL, a new one is allocated.
/** The size of the array is (orbit->size * bspec->n)**2
* (it makes sense to assume all the T-matrices share their spec).
*/
complex double *target,
_Complex double *target,
/// The orbit (type).
const qpms_ss_orbit_type_t *orbit,
/// Base spec of the t-matrices (we don't know it from orbit, as it has
@ -579,12 +583,12 @@ complex double *qpms_orbit_action_matrix(
/// Constructs a dense matrix representation of a irrep projector for an orbit type.
/** TODO detailed doc */
complex double *qpms_orbit_irrep_projector_matrix(
_Complex double *qpms_orbit_irrep_projector_matrix(
/// Target array. If NULL, a new one is allocated.
/** The size of the array is (orbit->size * bspec->n)**2
* (it makes sense to assume all the T-matrices share their spec).
*/
complex double *target,
_Complex double *target,
/// The orbit (type).
const qpms_ss_orbit_type_t *orbit,
/// Base spec of the t-matrices (we don't know it from orbit, as it has
@ -596,14 +600,14 @@ complex double *qpms_orbit_irrep_projector_matrix(
const qpms_iri_t iri);
/// TODO DOC!!!!!
complex double *qpms_orbit_irrep_basis(
_Complex double *qpms_orbit_irrep_basis(
/// Here theh size of theh basis shall be saved,
size_t *basis_size,
/// Target array. If NULL, a new one is allocated.
/** The size of the array is basis_size * (orbit->size * bspec->n)
* (it makes sense to assume all the T-matrices share their spec).
*/
complex double *target,
_Complex double *target,
/// The orbit (type).
const qpms_ss_orbit_type_t *orbit,
/// Base spec of the t-matrices (we don't know it from orbit, as it has
@ -618,10 +622,10 @@ complex double *qpms_orbit_irrep_basis(
/// Creates an incident field vector in the full basis, given a function that evaluates the field expansions at points.
/** TODO detailed doc!
* \returns target_full if target_full was not NULL, otherwise the newly allocated array. */
complex double *qpms_scatsys_incident_field_vector_full(
_Complex double *qpms_scatsys_incident_field_vector_full(
/// Target array. If NULL, a new one is allocated.
/** The length of the array is ss->fecv_size. */
complex double *target_full,
_Complex double *target_full,
const qpms_scatsys_t *ss,
qpms_incfield_t field_at_point,
const void *args, ///< Pointer passed as the last argument to (*field_at_point)()
@ -629,9 +633,9 @@ complex double *qpms_scatsys_incident_field_vector_full(
);
/// Applies T-matrices onto an incident field vector in the full basis.
complex double *qpms_scatsysw_apply_Tmatrices_full(
complex double *target_full, /// Target vector array. If NULL, a new one is allocated.
const complex double *inc_full, /// Incident field coefficient vector. Must not be NULL.
_Complex double *qpms_scatsysw_apply_Tmatrices_full(
_Complex double *target_full, /// Target vector array. If NULL, a new one is allocated.
const _Complex double *inc_full, /// Incident field coefficient vector. Must not be NULL.
const qpms_scatsys_at_omega_t *ssw
);
@ -650,7 +654,7 @@ struct beyn_result_t *qpms_scatsys_finite_find_eigenmodes(
const qpms_scatsys_t *ss,
/// A valid irrep index to search only in one irrep, or QPMS_NO_IRREP for solving the full system.
qpms_iri_t iri,
complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
_Complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
double omega_rr, ///< Real half-axis of the ellipse inside which the eigenfreqs are searched for.
double omega_ri, ///< Imaginary half-axis of the ellipse inside which the eigenfreqs are searched for.
size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number)
@ -673,7 +677,7 @@ struct beyn_result_t *qpms_scatsys_find_eigenmodes(
const qpms_scatsys_t *ss,
double eta, ///< Ewald sum parameter
const double *beta_, ///< k-vector of corresponding dimensionality, NULL/ignored for finite system.
complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
_Complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
double omega_rr, ///< Real half-axis of the ellipse inside which the eigenfreqs are searched for.
double omega_ri, ///< Imaginary half-axis of the ellipse inside which the eigenfreqs are searched for.
size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number)
@ -687,10 +691,10 @@ struct beyn_result_t *qpms_scatsys_find_eigenmodes(
#if 0
/// Creates a (partial) incident field vector in the symmetry-adapted basis, given a function that evaluates the field expansions at points.
/** TODO detailed doc! */
complex double *qpms_scatsys_incident_field_vector_irrep_packed(
_Complex double *qpms_scatsys_incident_field_vector_irrep_packed(
/// Target array. If NULL, a new one is allocated.
/** The length of the array is ss->fecv_size. */
complex double *target_full,
_Complex double *target_full,
const qpms_scatsys_t *ss,
const qpms_iri_t iri, ///< The index of given irreducible representation of ss->sym.
qpms_incfield_t field_at_point,
@ -706,7 +710,7 @@ complex double *qpms_scatsys_incident_field_vector_irrep_packed(
*
* \return Complex electric field at the point defined by \a evalpoint.
*
* \see qpms_scatsysw_scattered_E()
* \see qpms_scatsysw_scattered_E(), qpms_scatsys_scattered_field_basis()
*
* \see qpms_scatsyswk_scattered_E() for periodic systems.
*
@ -715,8 +719,8 @@ complex double *qpms_scatsys_incident_field_vector_irrep_packed(
ccart3_t qpms_scatsys_scattered_E(
const qpms_scatsys_t *ss,
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
complex double wavenumber, ///< Wavenumber of the background medium.
const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
_Complex double wavenumber, ///< Wavenumber of the background medium.
const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
@ -727,14 +731,96 @@ ccart3_t qpms_scatsys_scattered_E(
*
* \return Complex electric field at the point defined by \a evalpoint.
*
* \see qpms_scatsys_scattered_E()
* \see qpms_scatsys_scattered_E(), qpms_scatsysw_scattered_field_basis()
*
* \see qpms_scatsyswk_scattered_E() for periodic systems.
*/
ccart3_t qpms_scatsysw_scattered_E(
const qpms_scatsys_at_omega_t *ssw,
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
/// Evaluates a "basis" for electric field at a given point.
/**
* This function evaluates all the included VSWFs from the particles in the system, evaluated
* at a given point. Taking a linear combination of these with the coefficients \a scattcoeff_full[]
* would be equivalent to the result of qpms_scatsys_scattered_E().
*
* \see qpms_scatsysw_scattered_field_basis()
*
* \see qpms_scatsyswk_scattered_field_basis() for periodic systems.
*
*/
qpms_errno_t qpms_scatsys_scattered_field_basis(
ccart3_t *target, ///< Target array of length \a ss->fecv_size
const qpms_scatsys_t *ss,
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
_Complex double wavenumber, ///< Wavenumber of the background medium.
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
/// Evaluates a "basis" for electric field of a given particle at a given point.
/**
* This function evaluates all the included VSWFs from a particle in the system, evaluated
* at a given point.
*
* Finite systems only.
*
* \see qpms_scatsys_scattered_field_basis() for the whole-array equivalent.
*
* \see qpms_scatsysw_scattered_field_basis_pi()
*/
qpms_errno_t qpms_scatsys_scattered_field_basis_pi(
ccart3_t *target, ///< Target array of length at least `qpms_ss_bspec_pi(ssw->ss, pi)->n`
const qpms_scatsys_t *ss,
qpms_ss_pi_t pi, ///< Particle index
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
_Complex double wavenumber, ///< Wavenumber of the background medium
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
/// Evaluates a "basis" for electric field at a given point.
/**
* This function evaluates all the included VSWFs from the particles in the system, evaluated
* at a given point. Taking a linear combination of these with the coefficients \a scattcoeff_full[]
* would be equivalent to the result of qpms_scatsysw_scattered_E().
*
* Note that this might require a relatively large amount of memory. Depending
* on the application, qpms_scatsysw_scattered_field_basis_pi() might be a better
* alternative.
*
* \see qpms_scatsysw_scattered_field_basis_pi() for the single-particle equivalent.
*
* \see qpms_scatsys_scattered_field_basis()
*
* \see qpms_scatsyswk_scattered_field_basis() for periodic systems.
*
*/
qpms_errno_t qpms_scatsysw_scattered_field_basis(
ccart3_t *target, ///< Target array of length \a ss->fecv_size
const qpms_scatsys_at_omega_t *ssw,
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
/// Evaluates a "basis" for electric field of a given particle at a given point.
/**
* This function evaluates all the included VSWFs from a particle in the system, evaluated
* at a given point.
*
* Finite systems only.
*
* \see qpms_scatsysw_scattered_field_basis() for the whole-array equivalent.
*
* \see qpms_scatsys_scattered_field_basis_pi()
*/
qpms_errno_t qpms_scatsysw_scattered_field_basis_pi(
ccart3_t *target, ///< Target array of length at least `qpms_ss_bspec_pi(ssw->ss, pi)->n`
const qpms_scatsys_at_omega_t *ssw,
qpms_ss_pi_t pi, ///< Particle index
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
@ -750,8 +836,8 @@ ccart3_t qpms_scatsysw_scattered_E(
ccart3_t qpms_scatsys_scattered_E__alt(
const qpms_scatsys_t *ss,
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
complex double wavenumber, ///< Wavenumber of the background medium.
const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
_Complex double wavenumber, ///< Wavenumber of the background medium.
const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
@ -767,7 +853,7 @@ ccart3_t qpms_scatsys_scattered_E__alt(
ccart3_t qpms_scatsysw_scattered_E__alt(
const qpms_scatsys_at_omega_t *ssw,
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
@ -786,10 +872,18 @@ ccart3_t qpms_scatsysw_scattered_E__alt(
ccart3_t qpms_scatsyswk_scattered_E(
const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supported).
const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
ccart3_t qpms_scatsyswk_scattered_E_e(
const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supported).
const _Complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
cart3_t evalpoint, ///< A point \f$ \vect r \f$, at which the field is evaluated.
qpms_ewald_part parts
);
/// Evaluates "scattered" field basis functions in a periodic system.
/**
*
@ -802,10 +896,18 @@ qpms_errno_t qpms_scatsyswk_scattered_field_basis(
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the basis is evaluated.
);
qpms_errno_t qpms_scatsyswk_scattered_field_basis_e(
ccart3_t *target, ///< Target array of size sswk->ssw->ss->fecv_size
const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supponted).
cart3_t evalpoint, ///< A point \f$ \vect r \f$, at which the basis is evaluated.
qpms_ewald_part parts
);
/// Adjusted Ewadl parameter to avoid high-frequency breakdown.
// TODO DOC
double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, complex double wavenumber, const double wavevector[3]);
double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, _Complex double wavenumber, const double wavevector[3]);
#if 0
/** Evaluates partial scattered fields (corresponding to a given irrep-reduced excitation vector)
@ -815,10 +917,13 @@ double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, complex double wavenumber,
*/
ccart3_t qpms_scatsys_scattered_E_irrep(const qpms_scatsys_t *ss,
qpms_iri_t iri, ///< Irreducible representation
const complex double *coeff_vector, ///< A reduced excitation vector, corresponding to \a iri.
const _Complex double *coeff_vector, ///< A reduced excitation vector, corresponding to \a iri.
cart3_t where, ///< Evaluation point.
complex double k ///< Wave number.
_Complex double k ///< Wave number.
);
#endif
#ifdef __cplusplus
}
#endif
#endif //QPMS_SCATSYSTEM_H

48
qpms/scatsystem_dbg.c Normal file
View File

@ -0,0 +1,48 @@
#include "scatsystem_dbg.h"
#include "translations_dbg.h"
#include "indexing.h"
// analogue of qpms_scatsyswk_scattered_field_basis_e()
qpms_errno_t qpms_scatsyswk_test_sswf_basis_e(
complex double *target, ///< Target array of size (2 * sswk->ssw->ss->c->lMax + 1) * sswk->ssw->ss->p_count
const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t btyp, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supponted).
cart3_t where, ///< A point \f$ \vect r \f$, at which the basis is evaluated.
qpms_ewald_part parts
) {
QPMS_UNTESTED;
if (btyp != QPMS_HANKEL_PLUS)
QPMS_NOT_IMPLEMENTED("Only scattered field with first kind Hankel functions currently implemented.");
const qpms_scatsys_t *ss = sswk->ssw->ss;
if (ss->lattice_dimension != 2)
QPMS_NOT_IMPLEMENTED("Only 2D-periodic lattices implemented");
//ccart3_t res = {0,0,0};
//ccart3_t res_kc = {0,0,0}; // kahan sum compensation
const size_t particle_nelem = qpms_lMax2nelem_sc(2*ss->c->lMax + 1);
for (qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {
const cart3_t particle_pos = ss->p[pi].pos;
const cart3_t origin_cart = {0, 0, 0};
QPMS_ASSERT(sswk->k[2] == 0); // At least not implemented now
QPMS_ASSERT(ss->per.lattice_basis[0].z == 0);
QPMS_ASSERT(ss->per.lattice_basis[1].z == 0);
// Same choices as in qpms_ss_ppair_W32xy; TODO make it more dynamic
const double maxR = sqrt(ss->per.unitcell_volume) * 64;
const double maxK = 2048 * 2 * M_PI / maxR;
QPMS_ENSURE_SUCCESS(qpms_trans_calculator_test_sswf_e32(ss->c,
target + pi * particle_nelem, NULL,
sswk->eta, sswk->ssw->wavenumber,
cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]),
cart2_from_double_array(sswk->k), cart3_substract(where, particle_pos) /*CHECKSIGN*/,
maxR, maxK, parts));
}
return QPMS_SUCCESS;
}

10
qpms/scatsystem_dbg.h Normal file
View File

@ -0,0 +1,10 @@
#include "scatsystem.h"
qpms_errno_t qpms_scatsyswk_test_sswf_basis_e(
_Complex double *target, ///< Target array of size sswk->ssw->ss->fecv_size
const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supponted).
cart3_t evalpoint, ///< A point \f$ \vect r \f$, at which the basis is evaluated.
qpms_ewald_part parts
);

View File

@ -21,40 +21,44 @@
*/
#ifndef SYMMETRIES_H
#define SYMMETRIES_H
#ifdef __cplusplus
extern "C" {
#endif
#include "qpms_types.h"
#include <cblas.h>
/// Dense matrix representation of the z coordinate sign flip operation (xy-plane mirroring).
complex double *qpms_zflip_uvswi_dense(
complex double *target, ///< If NULL, a new array is allocated.
_Complex double *qpms_zflip_uvswi_dense(
_Complex double *target, ///< If NULL, a new array is allocated.
const qpms_vswf_set_spec_t *bspec);
/// Dense matrix representation of the y coordinate sign flip operation (xz-plane mirroring).
complex double *qpms_yflip_uvswi_dense(
complex double *target, ///< If NULL, a new array is allocated.
_Complex double *qpms_yflip_uvswi_dense(
_Complex double *target, ///< If NULL, a new array is allocated.
const qpms_vswf_set_spec_t *bspec);
/// Dense matrix representation of the x coordinate sign flip operation (yz-plane mirroring).
complex double *qpms_xflip_uvswi_dense(
complex double *target, ///< If NULL, a new array is allocated.
_Complex double *qpms_xflip_uvswi_dense(
_Complex double *target, ///< If NULL, a new array is allocated.
const qpms_vswf_set_spec_t *bspec);
/// Dense matrix representation of a rotation around the z-axis
complex double *qpms_zrot_uvswi_dense(
complex double *target, ///< If NULL, a new array is allocated.
_Complex double *qpms_zrot_uvswi_dense(
_Complex double *target, ///< If NULL, a new array is allocated.
const qpms_vswf_set_spec_t *bspec,
double phi ///< Rotation angle
);
/// Dense matrix representation of a "rational" rotation around the z-axis
/** Just for convenience. Corresponds to the angle \f$ \phi = 2\piw/N \f$.
*/
complex double *qpms_zrot_rational_uvswi_dense(
complex double *target, ///< If NULL, a new array is allocated.
_Complex double *qpms_zrot_rational_uvswi_dense(
_Complex double *target, ///< If NULL, a new array is allocated.
const qpms_vswf_set_spec_t *bspec,
int N,
int w
);
/// Dense matrix (uvswi-indexed) representation of any O(3) transformation.
complex double *qpms_irot3_uvswfi_dense(
complex double *target, ///< If NULL, a new array is allocated.
_Complex double *qpms_irot3_uvswfi_dense(
_Complex double *target, ///< If NULL, a new array is allocated.
const qpms_vswf_set_spec_t *bspec,
const qpms_irot3_t transf);
@ -74,5 +78,9 @@ size_t qpms_zero_roundoff_clean(double *arr, size_t nmemb, double atol);
* Works on real and imaginary parts separately.
* TODO doc.
*/
size_t qpms_czero_roundoff_clean(complex double *arr, size_t nmemb, double atol);
size_t qpms_czero_roundoff_clean(_Complex double *arr, size_t nmemb, double atol);
#ifdef __cplusplus
}
#endif
#endif // SYMMETRIES_H

View File

@ -1,6 +1,5 @@
import sympy
from sympy.combinatorics import Permutation, PermutationGroup
sympy.init_printing(perm_cyclic = True)
Permutation.print_cyclic = True
import cmath
from cmath import exp, pi
from math import sqrt

View File

@ -4,6 +4,7 @@
#ifndef TINY_INLINES_H
#define TINY_INLINES_H
#include <stdlib.h>
#include <complex.h>
static inline int min1pow(int pow) { return (pow % 2) ? -1 : 1; }
@ -19,8 +20,8 @@ static inline int min1pow_m_neg(int m) {
#if 0
#ifdef __GSL_SF_LEGENDRE_H__
static inline complex double
spharm_eval(gsl_sf_legendre_t P_normconv, int P_csphase, qpms_l_t l, qpms_m_t m, double P_n_abs_m, complex double exp_imf) {
static inline _Complex double
spharm_eval(gsl_sf_legendre_t P_normconv, int P_csphase, qpms_l_t l, qpms_m_t m, double P_n_abs_m, _Complex double exp_imf) {
return;
}
@ -28,9 +29,9 @@ spharm_eval(gsl_sf_legendre_t P_normconv, int P_csphase, qpms_l_t l, qpms_m_t m,
#endif
// this has shitty precision:
// static inline complex double ipow(int x) { return cpow(I, x); }
// static inline _Complex double ipow(int x) { return cpow(I, x); }
static inline complex double ipow(int x) {
static inline _Complex double ipow(int x) {
x = ((x % 4) + 4) % 4;
switch(x) {
case 0:

View File

@ -988,6 +988,12 @@ qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t,
return QPMS_SUCCESS;
}
qpms_errno_t qpms_tmatrix_generator_zero(qpms_tmatrix_t *t,
complex double omega_ignored, const void *params_ignored) {
memset(t->m, 0, SQ(t->spec->n)*sizeof(*(t->m)));
return QPMS_SUCCESS;
}
void qpms_tmatrix_operation_clear(qpms_tmatrix_operation_t *f) {
switch(f->typ) {
case QPMS_TMATRIX_OPERATION_NOOP:

View File

@ -3,6 +3,10 @@
*/
#ifndef TMATRICES_H
#define TMATRICES_H
#ifdef __cplusplus
extern "C" {
#endif
// #include "qpms_types.h" // included via materials.h
// #include <gsl/gsl_spline.h> // included via materials.h
#include "materials.h"
@ -14,7 +18,7 @@ struct qpms_finite_group_t;
typedef struct qpms_finite_group_t qpms_finite_group_t;
/// Returns a pointer to the beginning of the T-matrix row \a rowno.
static inline complex double *qpms_tmatrix_row(qpms_tmatrix_t *t, size_t rowno){
static inline _Complex double *qpms_tmatrix_row(qpms_tmatrix_t *t, size_t rowno){
return t->m + (t->spec->n * rowno);
}
@ -42,7 +46,7 @@ void qpms_tmatrix_free(qpms_tmatrix_t *t);
* This function actually checks for identical vswf specs.
* TODO define constants with "default" atol, rtol for this function.
*/
bool qpms_tmatrix_isclose(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2,
_Bool qpms_tmatrix_isclose(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2,
const double rtol, const double atol);
/// Creates a T-matrix from another matrix and a symmetry operation.
@ -53,7 +57,7 @@ bool qpms_tmatrix_isclose(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2,
*/
qpms_tmatrix_t *qpms_tmatrix_apply_symop(
const qpms_tmatrix_t *T, ///< the original T-matrix
const complex double *M ///< the symmetry op matrix in row-major format
const _Complex double *M ///< the symmetry op matrix in row-major format
);
/// Applies a symmetry operation onto a T-matrix, rewriting the original T-matrix data.
@ -64,7 +68,7 @@ qpms_tmatrix_t *qpms_tmatrix_apply_symop(
*/
qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace(
qpms_tmatrix_t *T, ///< the original T-matrix
const complex double *M ///< the symmetry op matrix in row-major format
const _Complex double *M ///< the symmetry op matrix in row-major format
);
/// Symmetrizes a T-matrix with an involution symmetry operation.
@ -75,7 +79,7 @@ qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace(
*/
qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution(
const qpms_tmatrix_t *T, ///< the original T-matrix
const complex double *M ///< the symmetry op matrix in row-major format
const _Complex double *M ///< the symmetry op matrix in row-major format
);
/// Creates a \f$ C_\infty \f$ -symmetrized version of a T-matrix.
@ -105,7 +109,7 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N(
*/
qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution_inplace(
qpms_tmatrix_t *T, ///< the original T-matrix
const complex double *M ///< the symmetry op matrix in row-major format
const _Complex double *M ///< the symmetry op matrix in row-major format
);
/// Creates a \f$ C_\infty \f$ -symmetrized version of a T-matrix, rewriting the original one.
@ -149,7 +153,7 @@ qpms_errno_t qpms_load_scuff_tmatrix(
double **freqs_su, ///< Frequencies in SCUFF units (optional).
/// The resulting T-matrices (optional).
qpms_tmatrix_t **tmatrices_array,
complex double **tmdata ///< The T-matrices raw contents
_Complex double **tmdata ///< The T-matrices raw contents
);
/// Tells whether qpms_load_scuff_tmatrix should crash if fopen() fails.
@ -161,7 +165,7 @@ qpms_errno_t qpms_load_scuff_tmatrix(
* This is desirable e.g. when used in Python (so that proper exception can
* be thrown).
*/
extern bool qpms_load_scuff_tmatrix_crash_on_failure;
extern _Bool qpms_load_scuff_tmatrix_crash_on_failure;
/// Loads a scuff-tmatrix generated file.
/** A simple wrapper over qpms_read_scuff_tmatrix() that needs a
@ -189,7 +193,7 @@ qpms_errno_t qpms_read_scuff_tmatrix(
* is accessed via
* (*tmdata)[bspec->n*bspec->n*fi + desti*bspec->n + srci].
*/
complex double ** tmdata
_Complex double ** tmdata
);
/// In-place application of point group elements on raw T-matrix data.
@ -200,7 +204,7 @@ qpms_errno_t qpms_read_scuff_tmatrix(
* TODO more doc.
*/
qpms_errno_t qpms_symmetrise_tmdata_irot3arr(
complex double *tmdata, const size_t tmcount,
_Complex double *tmdata, const size_t tmcount,
const qpms_vswf_set_spec_t *bspec,
size_t n_symops,
const qpms_irot3_t *symops
@ -213,7 +217,7 @@ qpms_errno_t qpms_symmetrise_tmdata_irot3arr(
* TODO more doc.
*/
qpms_errno_t qpms_symmetrise_tmdata_finite_group(
complex double *tmdata, const size_t tmcount,
_Complex double *tmdata, const size_t tmcount,
const qpms_vswf_set_spec_t *bspec,
const qpms_finite_group_t *pointgroup
);
@ -242,9 +246,9 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace(
);
/// Application of T-matrix on a vector of incident field coefficients, \f$ f = Ta \f$.
complex double *qpms_apply_tmatrix(
complex double *f_target, ///< Scattered field coefficient array of size T->spec->n; if NULL, a new one is allocated.
const complex double *a, ///< Incident field coefficient array of size T->spec->n.
_Complex double *qpms_apply_tmatrix(
_Complex double *f_target, ///< Scattered field coefficient array of size T->spec->n; if NULL, a new one is allocated.
const _Complex double *a, ///< Incident field coefficient array of size T->spec->n.
const qpms_tmatrix_t *T ///< T-matrix \a T to apply.
);
@ -254,10 +258,11 @@ complex double *qpms_apply_tmatrix(
* qpms_tmatrix_generator_interpolator()
* qpms_tmatrix_generator_sphere()
* qpms_tmatrix_generator_constant()
* qpms_tmatrix_generator_zero()
*/
typedef struct qpms_tmatrix_generator_t {
qpms_errno_t (*function) (qpms_tmatrix_t *t, ///< T-matrix to fill.
complex double omega, ///< Angular frequency.
_Complex double omega, ///< Angular frequency.
const void *params ///< Implementation dependent parameters.
);
const void *params; ///< Parameter pointer passed to the function.
@ -267,7 +272,7 @@ typedef struct qpms_tmatrix_generator_t {
qpms_tmatrix_t *qpms_tmatrix_init_from_generator(
const qpms_vswf_set_spec_t *bspec,
qpms_tmatrix_generator_t gen,
complex double omega);
_Complex double omega);
/// Implementation of qpms_matrix_generator_t that just copies a constant matrix.
@ -275,11 +280,20 @@ qpms_tmatrix_t *qpms_tmatrix_init_from_generator(
* the same base spec.
*/
qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t,
complex double omega,
_Complex double omega_ignored,
/// Source T-matrix, real type is (const qpms_tmatrix_t*).
const void *tmatrix_orig
);
/// Dummy implementation of qpms_tmatrix_generator_t, yielding a zero matrix.
/**
* This effectively represents a part of the background medium without
* any scatterer present.
*/
qpms_errno_t qpms_tmatrix_generator_zero(qpms_tmatrix_t *t,
_Complex double omega_ignored,
const void *params_ignored);
/* Fuck this, include the whole <gsl/gsl_spline.h>
typedef struct gsl_spline gsl_spline; // Forward declaration for the interpolator struct
typedef struct gsl_interp_type gsl_interp_type;
@ -295,7 +309,7 @@ extern const gsl_interp_type * gsl_interp_steffen;
// struct gsl_interp_accel; // use if lookup proves to be too slow
typedef struct qpms_tmatrix_interpolator_t {
const qpms_vswf_set_spec_t *bspec;
//bool owns_bspec;
//_Bool owns_bspec;
gsl_spline **splines_real; ///< There will be a spline object for each nonzero element
gsl_spline **splines_imag; ///< There will be a spline object for each nonzero element
// gsl_interp_accel **accel_real;
@ -317,7 +331,7 @@ qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t
qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Number of freqs and T-matrices provided.
const double *freqs, const qpms_tmatrix_t *tmatrices_array, ///< N.B. array of qpms_tmatrix_t, not pointers!
const gsl_interp_type *iptype
//, bool copy_bspec ///< if true, copies its own copy of basis spec from the first T-matrix.
//, _Bool copy_bspec ///< if true, copies its own copy of basis spec from the first T-matrix.
/*, ...? */);
@ -325,7 +339,7 @@ qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Num
/** As in qpms_tmatrix_interpolator_eval(), the imaginary part of frequency is discarded!
*/
qpms_errno_t qpms_tmatrix_generator_interpolator(qpms_tmatrix_t *t, ///< T-matrix to fill.
complex double omega, ///< Angular frequency.
_Complex double omega, ///< Angular frequency.
const void *interpolator ///< Parameter of type qpms_tmatrix_interpolator_t *.
);
@ -342,14 +356,14 @@ qpms_errno_t qpms_tmatrix_generator_interpolator(qpms_tmatrix_t *t, ///< T-matri
*
* TODO better doc.
*/
complex double *qpms_mie_coefficients_reflection(
complex double *target, ///< Target array of length bspec->n. If NULL, a new one will be allocated.
_Complex double *qpms_mie_coefficients_reflection(
_Complex double *target, ///< Target array of length bspec->n. If NULL, a new one will be allocated.
const qpms_vswf_set_spec_t *bspec, ///< Defines which of the coefficients are calculated.
double a, ///< Radius of the sphere.
complex double k_i, ///< Wave number of the internal material of the sphere.
complex double k_e, ///< Wave number of the surrounding medium.
complex double mu_i, ///< Relative permeability of the sphere material.
complex double mu_e, ///< Relative permeability of the surrounding medium.
_Complex double k_i, ///< Wave number of the internal material of the sphere.
_Complex double k_e, ///< Wave number of the surrounding medium.
_Complex double mu_i, ///< Relative permeability of the sphere material.
_Complex double mu_e, ///< Relative permeability of the surrounding medium.
qpms_bessel_t J_ext, ///< Kind of the "incoming" waves. Most likely QPMS_BESSEL_REGULAR.
qpms_bessel_t J_scat ///< Kind of the "scattered" waves. Most likely QPMS_HANKEL_PLUS.
);
@ -357,10 +371,10 @@ complex double *qpms_mie_coefficients_reflection(
/// Replaces the contents of an existing T-matrix with that of a spherical nanoparticle calculated using the Lorentz-mie theory.
qpms_errno_t qpms_tmatrix_spherical_fill(qpms_tmatrix_t *t, ///< T-matrix whose contents are to be replaced. Not NULL.
double a, ///< Radius of the sphere.
complex double k_i, ///< Wave number of the internal material of the sphere.
complex double k_e, ///< Wave number of the surrounding medium.
complex double mu_i, ///< Relative permeability of the sphere material.
complex double mu_e ///< Relative permeability of the surrounding medium.
_Complex double k_i, ///< Wave number of the internal material of the sphere.
_Complex double k_e, ///< Wave number of the surrounding medium.
_Complex double mu_i, ///< Relative permeability of the sphere material.
_Complex double mu_e ///< Relative permeability of the surrounding medium.
);
/// Parameter structure for qpms_tmatrix_generator_sphere().
@ -372,7 +386,7 @@ typedef struct qpms_tmatrix_generator_sphere_param_t {
/// T-matrix generator for spherical particles using Lorentz-Mie solution.
qpms_errno_t qpms_tmatrix_generator_sphere(qpms_tmatrix_t *t,
complex double omega,
_Complex double omega,
const void *params ///< Of type qpms_tmatrix_generator_sphere_param_t.
);
@ -380,10 +394,10 @@ qpms_errno_t qpms_tmatrix_generator_sphere(qpms_tmatrix_t *t,
static inline qpms_tmatrix_t *qpms_tmatrix_spherical(
const qpms_vswf_set_spec_t *bspec,
double a, ///< Radius of the sphere.
complex double k_i, ///< Wave number of the internal material of the sphere.
complex double k_e, ///< Wave number of the surrounding medium.
complex double mu_i, ///< Relative permeability of the sphere material.
complex double mu_e ///< Relative permeability of the surrounding medium.
_Complex double k_i, ///< Wave number of the internal material of the sphere.
_Complex double k_e, ///< Wave number of the surrounding medium.
_Complex double mu_i, ///< Relative permeability of the sphere material.
_Complex double mu_e ///< Relative permeability of the surrounding medium.
) {
qpms_tmatrix_t *t = qpms_tmatrix_init(bspec);
qpms_tmatrix_spherical_fill(t, a, k_i, k_e, mu_i, mu_e);
@ -395,8 +409,8 @@ qpms_errno_t qpms_tmatrix_spherical_mu0_fill(
qpms_tmatrix_t *t, ///< T-matrix whose contents are to be replaced. Not NULL.
double a, ///< Radius of the sphere.
double omega, ///< Angular frequency.
complex double epsilon_fg, ///< Relative permittivity of the sphere.
complex double epsilon_bg ///< Relative permittivity of the background medium.
_Complex double epsilon_fg, ///< Relative permittivity of the sphere.
_Complex double epsilon_bg ///< Relative permittivity of the background medium.
);
/// Convenience function to calculate T-matrix of a non-magnetic spherical particle using the permittivity values.
@ -404,8 +418,8 @@ static inline qpms_tmatrix_t *qpms_tmatrix_spherical_mu0(
const qpms_vswf_set_spec_t *bspec,
double a, ///< Radius of the sphere.
double omega, ///< Angular frequency.
complex double epsilon_fg, ///< Relative permittivity of the sphere.
complex double epsilon_bg ///< Relative permittivity of the background medium.
_Complex double epsilon_fg, ///< Relative permittivity of the sphere.
_Complex double epsilon_bg ///< Relative permittivity of the background medium.
) {
qpms_tmatrix_t *t = qpms_tmatrix_init(bspec);
qpms_tmatrix_spherical_mu0_fill(t, a, omega, epsilon_fg, epsilon_bg);
@ -460,7 +474,7 @@ qpms_arc_function_retval_t qpms_arc_sphere(double theta,
*/
qpms_errno_t qpms_tmatrix_axialsym_fill(
qpms_tmatrix_t *t, ///< T-matrix whose contents are to be replaced. Not NULL.
complex double omega, ///< Angular frequency.
_Complex double omega, ///< Angular frequency.
qpms_epsmu_t outside, ///< Optical properties of the outside medium.
qpms_epsmu_t inside, ///< Optical properties of the particle's material.
qpms_arc_function_t shape, ///< Particle surface parametrisation.
@ -474,7 +488,7 @@ qpms_errno_t qpms_tmatrix_axialsym_fill(
/// Creates a new T-matrix of a particle with \f$ C_\infty \f$ symmetry.
static inline qpms_tmatrix_t *qpms_tmatrix_axialsym(
const qpms_vswf_set_spec_t *bspec,
complex double omega, ///< Angular frequency.
_Complex double omega, ///< Angular frequency.
qpms_epsmu_t outside, ///< Optical properties of the outside medium.
qpms_epsmu_t inside, ///< Optical properties of the particle's material.
qpms_arc_function_t shape, ///< Particle surface parametrisation.
@ -501,22 +515,22 @@ typedef struct qpms_tmatrix_generator_axialsym_param_t {
/// qpms_tmatrix_axialsym for qpms_tmatrix_generator_t
qpms_errno_t qpms_tmatrix_generator_axialsym(qpms_tmatrix_t *t, ///< T-matrix to fill.
complex double omega, ///< Angular frequency.
_Complex double omega, ///< Angular frequency.
const void *params ///< Parameters of type qpms_tmatrix_generator_axialsym_param_t.
);
/// Computes the (reduced) transposed R or Q matrix for axially symmetric particle (useful for debugging).
qpms_errno_t qpms_tmatrix_generator_axialsym_RQ_transposed_fill(complex double *target,
complex double omega,
qpms_errno_t qpms_tmatrix_generator_axialsym_RQ_transposed_fill(_Complex double *target,
_Complex double omega,
const qpms_tmatrix_generator_axialsym_param_t *param,
qpms_normalisation_t norm,
qpms_bessel_t J
);
/// Computes the (reduced) transposed R or Q matrix for axially symmetric particle (useful mostly for debugging).
qpms_errno_t qpms_tmatrix_axialsym_RQ_transposed_fill(complex double *target,
complex double omega, qpms_epsmu_t outside, qpms_epsmu_t inside,
qpms_errno_t qpms_tmatrix_axialsym_RQ_transposed_fill(_Complex double *target,
_Complex double omega, qpms_epsmu_t outside, qpms_epsmu_t inside,
qpms_arc_function_t shape, qpms_l_t lMaxQR, qpms_normalisation_t norm,
qpms_bessel_t J ///< Use QPMS_BESSEL_REGULAR to calculate \f$ R^T\f$ or QPMS_HANKEL_PLUS to calculate \f$ Q^T\f$.
);
@ -537,7 +551,7 @@ typedef struct qpms_tmatrix_function_t {
/// Convenience function to create a new T-matrix from qpms_tmatrix_function_t.
// FIXME the name is not very intuitive.
static inline qpms_tmatrix_t *qpms_tmatrix_init_from_function(qpms_tmatrix_function_t f, complex double omega) {
static inline qpms_tmatrix_t *qpms_tmatrix_init_from_function(qpms_tmatrix_function_t f, _Complex double omega) {
return qpms_tmatrix_init_from_generator(f.spec, *f.gen, omega);
}
@ -558,9 +572,9 @@ typedef enum {
struct qpms_tmatrix_operation_lrmatrix {
/// Raw matrix data of \a M in row-major order.
/** The matrix must be taylored for the given bspec! */
complex double *m;
size_t m_size; ///< Total size of \a m matrix in terms of sizeof(complex double).
bool owns_m; ///< Whether \a m is owned by this;
_Complex double *m;
size_t m_size; ///< Total size of \a m matrix in terms of sizeof(_Complex double).
_Bool owns_m; ///< Whether \a m is owned by this;
};
struct qpms_tmatrix_operation_t; // Forward declaration for the composed operations.
@ -591,9 +605,9 @@ struct qpms_tmatrix_operation_compose_chain {
struct qpms_tmatrix_operation_scmulz {
/// Raw matrix data of \a M in row-major order.
/** The matrix must be taylored for the given bspec! */
complex double *m;
size_t m_size; ///< Total size of \a m matrix in terms of sizeof(complex double).
bool owns_m; ///< Whether \a m is owned by this.
_Complex double *m;
size_t m_size; ///< Total size of \a m matrix in terms of sizeof(_Complex double).
_Bool owns_m; ///< Whether \a m is owned by this.
};
/// Specifies a symmetrisation using a set of rotoreflections (with equal weights) for qpms_tmatrix_operation_t.
@ -602,7 +616,7 @@ struct qpms_tmatrix_operation_scmulz {
struct qpms_tmatrix_operation_irot3arr {
size_t n; ///< Number of rotoreflections;
qpms_irot3_t *ops; ///< Rotoreflection array of size \a n.
bool owns_ops; ///< Whether \a ops array is owned by this.
_Bool owns_ops; ///< Whether \a ops array is owned by this.
};
/// A generic T-matrix transformation operator.
@ -689,4 +703,7 @@ typedef struct qpms_abstract_particle_t {
typedef qpms_particle_tid_t qpms_abstract_particle_tid_t;
#endif // 0
#ifdef __cplusplus
}
#endif
#endif //TMATRICES_H

View File

@ -12,7 +12,8 @@ except ImportError:
use_moble_quaternion = False
import re
from .constants import hbar, eV, pi, c
from scipy import interpolate
from scipy.constants import hbar, e as eV, pi, c
from .cycommon import get_mn_y, get_nelem
from .cyquaternions import CQuat
ň = np.newaxis
@ -439,7 +440,6 @@ class TMatrix(TMatrixSpec):
TODO support for different/multiple interpolators
'''
def __init__(self, tmatrix_spec):
from scipy import interpolate
#self.specification = tmatrix_spec
self.lMax_override = tmatrix_spec.lMax_override
self.tmatrix_path = tmatrix_spec.tmatrix_path

View File

@ -1,6 +1,9 @@
/*! \file tolerances.h */
#ifndef QPMS_TOLERANCES_H
#define QPMS_TOLERANCES_H
#ifdef __cplusplus
extern "C" {
#endif
// TODO DOC
@ -12,4 +15,7 @@ typedef struct qpms_tolerance_spec_t {
/// A rather arbitrary default tolerance.
static const qpms_tolerance_spec_t QPMS_TOLERANCE_DEFAULT = {.atol = 1e-9, .rtol = 1e-8};
#ifdef __cplusplus
}
#endif
#endif // QPMS_TOLERANCES_H

View File

@ -1,5 +1,6 @@
#include <math.h>
#include "qpms_types.h"
#include <complex.h>
#include "qpms_specfunc.h"
#include "gaunt.h"
#include "translations.h"
@ -15,6 +16,10 @@
#include "normalisation.h"
#include "translations_inlines.h"
#if defined LATTICESUMS31 || defined LATTICESUMS32
#include "lattices.h"
#endif
/*
* Define macros with additional factors that "should not be there" according
* to the "original" formulae but are needed to work with my vswfs.
@ -825,7 +830,7 @@ int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_tr
)
{
const qpms_y_t nelem2_sc = qpms_lMax2nelem_sc(c->e3c->lMax);
const qpms_y_sc_t nelem2_sc = qpms_lMax2nelem_sc(c->e3c->lMax);
//const qpms_y_t nelem = qpms_lMax2nelem(c->lMax);
const bool doerr = Aerr || Berr;
const bool do_sigma0 = (particle_shift == 0)//DIFF21((particle_shift.x == 0) && (particle_shift.y == 0)); // FIXME ignoring the case where particle_shift equals to lattice vector
@ -846,15 +851,15 @@ int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_tr
QPMS_ENSURE_SUCCESS(ewald31z_sigma_short_points_and_shift(sigmas_short, serr_short, //DIFF21
c->e3c, eta, k, nRpoints, Rpoints, beta, particle_shift));
for(qpms_y_t y = 0; y < nelem2_sc; ++y)
for(qpms_y_sc_t y = 0; y < nelem2_sc; ++y)
sigmas_total[y] = sigmas_short[y] + sigmas_long[y];
if (doerr) for(qpms_y_t y = 0; y < nelem2_sc; ++y)
if (doerr) for(qpms_y_sc_t y = 0; y < nelem2_sc; ++y)
serr_total[y] = serr_short[y] + serr_long[y];
complex double sigma0 = 0; double sigma0_err = 0;
if (do_sigma0) {
QPMS_ENSURE_SUCCESS(ewald31z_sigma0(&sigma0, &sigma0_err, c->e3c, eta, k));
const qpms_l_t y = qpms_mn2y_sc(0,0);
const qpms_y_sc_t y = qpms_mn2y_sc(0,0);
sigmas_total[y] += sigma0;
if(doerr) serr_total[y] += sigma0_err;
}
@ -872,7 +877,7 @@ int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_tr
// TODO skip if ... (N.B. skip will be different for 31z and 32)
for(qpms_l_t q = 0; q <= qmax; ++q) {
const qpms_l_t p = n + nu - 2*q;
const qpms_y_t y_sc = qpms_mn2y_sc(mu_m, p);
const qpms_y_sc_t y_sc = qpms_mn2y_sc(mu_m, p);
const complex double multiplier = c->A_multipliers[i][q];
complex double sigma = sigmas_total[y_sc];
ckahanadd(&Asum, &Asumc, multiplier * sigma);
@ -887,7 +892,7 @@ int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_tr
double Bsumerr, Bsumerrc; if(Berr) kahaninit(&Bsumerr, &Bsumerrc);
for(qpms_l_t q = 0; q <= qmax; ++q) {
const qpms_l_t p_ = n + nu - 2*q + 1;
const qpms_y_t y_sc = qpms_mn2y_sc(mu_m, p_);
const qpms_y_sc_t y_sc = qpms_mn2y_sc(mu_m, p_);
const complex double multiplier = c->B_multipliers[i][q-BQ_OFFSET];
complex double sigma = sigmas_total[y_sc];
ckahanadd(&Bsum, &Bsumc, multiplier * sigma);
@ -937,7 +942,7 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
)
{
const qpms_y_t nelem2_sc = qpms_lMax2nelem_sc(c->e3c->lMax);
const qpms_y_sc_t nelem2_sc = qpms_lMax2nelem_sc(c->e3c->lMax);
//const qpms_y_t nelem = qpms_lMax2nelem(c->lMax);
const bool doerr = Aerr || Berr;
const bool do_sigma0 = ((particle_shift.x == 0) && (particle_shift.y == 0) && (particle_shift.z == 0)); // FIXME ignoring the case where particle_shift equals to lattice vector
@ -1006,17 +1011,17 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
PGen_destroy(&Rgen);
}
for(qpms_y_t y = 0; y < nelem2_sc; ++y)
for(qpms_y_sc_t y = 0; y < nelem2_sc; ++y)
sigmas_total[y] = ((parts & QPMS_EWALD_SHORT_RANGE) ? sigmas_short[y] : 0)
+ ((parts & QPMS_EWALD_LONG_RANGE) ? sigmas_long[y] : 0);
if (doerr) for(qpms_y_t y = 0; y < nelem2_sc; ++y)
if (doerr) for(qpms_y_sc_t y = 0; y < nelem2_sc; ++y)
serr_total[y] = ((parts & QPMS_EWALD_SHORT_RANGE) ? serr_short[y] : 0)
+ ((parts & QPMS_EWALD_LONG_RANGE) ? serr_long[y] : 0);
complex double sigma0 = 0; double sigma0_err = 0;
if (do_sigma0 && (parts & QPMS_EWALD_0TERM)) {
QPMS_ENSURE_SUCCESS(ewald3_sigma0(&sigma0, &sigma0_err, c->e3c, eta, k));
const qpms_l_t y = qpms_mn2y_sc(0,0);
const qpms_y_sc_t y = qpms_mn2y_sc(0,0);
sigmas_total[y] += sigma0;
if(doerr) serr_total[y] += sigma0_err;
}
@ -1034,7 +1039,7 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
// TODO skip if ...
for(qpms_l_t q = 0; q <= qmax; ++q) {
const qpms_l_t p = n + nu - 2*q;
const qpms_y_t y_sc = qpms_mn2y_sc(mu_m, p);
const qpms_y_sc_t y_sc = qpms_mn2y_sc(mu_m, p);
const complex double multiplier = c->A_multipliers[i][q];
complex double sigma = sigmas_total[y_sc];
ckahanadd(&Asum, &Asumc, multiplier * sigma);
@ -1049,7 +1054,7 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
double Bsumerr, Bsumerrc; if(Berr) kahaninit(&Bsumerr, &Bsumerrc);
for(qpms_l_t q = 0; q <= qmax; ++q) {
const qpms_l_t p_ = n + nu - 2*q + 1;
const qpms_y_t y_sc = qpms_mn2y_sc(mu_m, p_);
const qpms_y_sc_t y_sc = qpms_mn2y_sc(mu_m, p_);
const complex double multiplier = c->B_multipliers[i][q-BQ_OFFSET];
complex double sigma = sigmas_total[y_sc];
ckahanadd(&Bsum, &Bsumc, multiplier * sigma);

View File

@ -26,9 +26,12 @@
*/
#ifndef QPMS_TRANSLATIONS_H
#define QPMS_TRANSLATIONS_H
#ifdef __cplusplus
extern "C" {
#endif
#include "vectors.h"
#include "qpms_types.h"
#include <complex.h>
#include <stdbool.h>
#include <stddef.h>
@ -37,10 +40,10 @@
#endif
complex double qpms_trans_single_A(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
_Complex double qpms_trans_single_A(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
bool r_ge_d, qpms_bessel_t J);
complex double qpms_trans_single_B(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
_Complex double qpms_trans_single_B(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
bool r_ge_d, qpms_bessel_t J);
/// Structure holding the constant factors in normalisation operators.
@ -65,8 +68,8 @@ typedef struct qpms_trans_calculator {
qpms_normalisation_t normalisation;
qpms_l_t lMax;
qpms_y_t nelem;
complex double **A_multipliers;
complex double **B_multipliers;
_Complex double **A_multipliers;
_Complex double **B_multipliers;
#if 0
// Normalised values of the Legendre functions and derivatives
// for θ == π/2, i.e. for the 2D case.
@ -90,44 +93,44 @@ qpms_trans_calculator *qpms_trans_calculator_init(qpms_l_t lMax, ///< Truncation
/// Destructor for qpms_trans_calculator_t.
void qpms_trans_calculator_free(qpms_trans_calculator *);
complex double qpms_trans_calculator_get_A(const qpms_trans_calculator *c,
_Complex double qpms_trans_calculator_get_A(const qpms_trans_calculator *c,
qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
bool r_ge_d, qpms_bessel_t J);
complex double qpms_trans_calculator_get_B(const qpms_trans_calculator *c,
_Complex double qpms_trans_calculator_get_B(const qpms_trans_calculator *c,
qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
bool r_ge_d, qpms_bessel_t J);
int qpms_trans_calculator_get_AB_p(const qpms_trans_calculator *c,
complex double *Adest, complex double *Bdest,
_Complex double *Adest, _Complex double *Bdest,
qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, csph_t kdlj,
bool r_ge_d, qpms_bessel_t J);
int qpms_trans_calculator_get_AB_arrays(const qpms_trans_calculator *c,
complex double *Adest, complex double *Bdest,
_Complex double *Adest, _Complex double *Bdest,
size_t deststride, size_t srcstride,
csph_t kdlj, bool r_ge_d, qpms_bessel_t J);
// TODO update the types later
complex double qpms_trans_calculator_get_A_ext(const qpms_trans_calculator *c,
int m, int n, int mu, int nu, complex double kdlj_r,
_Complex double qpms_trans_calculator_get_A_ext(const qpms_trans_calculator *c,
int m, int n, int mu, int nu, _Complex double kdlj_r,
double kdlj_th, double kdlj_phi, int r_ge_d, int J);
complex double qpms_trans_calculator_get_B_ext(const qpms_trans_calculator *c,
int m, int n, int mu, int nu, complex double kdlj_r,
_Complex double qpms_trans_calculator_get_B_ext(const qpms_trans_calculator *c,
int m, int n, int mu, int nu, _Complex double kdlj_r,
double kdlj_th, double kdlj_phi, int r_ge_d, int J);
int qpms_trans_calculator_get_AB_p_ext(const qpms_trans_calculator *c,
complex double *Adest, complex double *Bdest,
int m, int n, int mu, int nu, complex double kdlj_r,
_Complex double *Adest, _Complex double *Bdest,
int m, int n, int mu, int nu, _Complex double kdlj_r,
double kdlj_th, double kdlj_phi, int r_ge_d, int J);
int qpms_trans_calculator_get_AB_arrays_ext(const qpms_trans_calculator *c,
complex double *Adest, complex double *Bdest,
_Complex double *Adest, _Complex double *Bdest,
size_t deststride, size_t srcstride,
complex double kdlj_r, double kdlj_theta, double kdlj_phi,
_Complex double kdlj_r, double kdlj_theta, double kdlj_phi,
int r_ge_d, int J);
// Convenience functions using VSWF base specs
qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator *c,
complex double *target,
_Complex double *target,
/// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
const qpms_vswf_set_spec_t *destspec, size_t deststride,
/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
@ -138,12 +141,12 @@ qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator *
/// and with automatic \a r_ge_d = `false`.
qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p(
const qpms_trans_calculator *c,
complex double *target,
_Complex double *target,
/// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
const qpms_vswf_set_spec_t *destspec, size_t deststride,
/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
const qpms_vswf_set_spec_t *srcspec, size_t srcstride,
complex double k, cart3_t destpos, cart3_t srcpos,
_Complex double k, cart3_t destpos, cart3_t srcpos,
qpms_bessel_t J
/// Workspace has to be at least 2 * c->neleme**2 long
);
@ -156,10 +159,10 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p(
int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
complex double *Adest, double *Aerr,
complex double *Bdest, double *Berr,
_Complex double *Adest, double *Aerr,
_Complex double *Bdest, double *Berr,
const ptrdiff_t deststride, const ptrdiff_t srcstride,
const double eta, const complex double wavenumber,
const double eta, const _Complex double wavenumber,
cart2_t b1, cart2_t b2,
const cart2_t beta,
const cart3_t particle_shift,
@ -167,10 +170,10 @@ int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
);
int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
complex double *Adest, double *Aerr,
complex double *Bdest, double *Berr,
_Complex double *Adest, double *Aerr,
_Complex double *Bdest, double *Berr,
const ptrdiff_t deststride, const ptrdiff_t srcstride,
const double eta, const complex double wavenumber,
const double eta, const _Complex double wavenumber,
cart2_t b1, cart2_t b2,
const cart2_t beta,
const cart3_t particle_shift,
@ -180,12 +183,12 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
// Convenience functions using VSWF base specs
qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculator *c,
complex double *target, double *err,
_Complex double *target, double *err,
/// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
const qpms_vswf_set_spec_t *destspec, size_t deststride,
/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
const qpms_vswf_set_spec_t *srcspec, size_t srcstride,
const double eta, const complex double wavenumber,
const double eta, const _Complex double wavenumber,
cart2_t b1, cart2_t b2,
const cart2_t beta,
const cart3_t particle_shift,
@ -193,12 +196,12 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculat
);
qpms_errno_t qpms_trans_calculator_get_trans_array_e32_e(const qpms_trans_calculator *c,
complex double *target, double *err,
_Complex double *target, double *err,
/// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
const qpms_vswf_set_spec_t *destspec, size_t deststride,
/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
const qpms_vswf_set_spec_t *srcspec, size_t srcstride,
const double eta, const complex double wavenumber,
const double eta, const _Complex double wavenumber,
cart2_t b1, cart2_t b2,
const cart2_t beta,
const cart3_t particle_shift,
@ -212,8 +215,8 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_e32_e(const qpms_trans_calcul
// e31z means that the particles are positioned along the z-axis;
// their positions and K-values are then denoted by a single z-coordinate
int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_trans_calculator *c,
complex double *Adest, double *Aerr,
complex double *Bdest, double *Berr,
_Complex double *Adest, double *Aerr,
_Complex double *Bdest, double *Berr,
const ptrdiff_t deststride, const ptrdiff_t srcstride,
const double eta, const double k,
const double unitcell_area, // just lattice period
@ -244,4 +247,7 @@ int qpms_cython_trans_calculator_get_AB_arrays_loop(
#endif //QPMS_COMPILE_PYTHON_EXTENSIONS
#ifdef __cplusplus
}
#endif
#endif // QPMS_TRANSLATIONS_H

149
qpms/translations_dbg.c Normal file
View File

@ -0,0 +1,149 @@
#include "qpms_types.h"
#include "qpms_specfunc.h"
#include "translations_dbg.h"
#include "indexing.h"
#include "tiny_inlines.h"
#include "qpms_error.h"
#include "normalisation.h"
#include "lattices.h"
int qpms_trans_calculator_test_sswf(const qpms_trans_calculator *c,
complex double *dest, const csph_t kdlj, const qpms_bessel_t J) {
if (0 == kdlj.r && J != QPMS_BESSEL_REGULAR) {
for (qpms_l_t l = 0; l <= 2*c->lMax+1; ++l)
for (qpms_m_t m = 0; m <= l; ++m) {
const qpms_y_t y_sc = qpms_mn2y_sc(m, l);
dest[y_sc] = NAN+I*NAN;
}
// TODO warn? different return value?
return 0;
}
double legendre_buf[gsl_sf_legendre_array_n(c->lMax+c->lMax+1)];
complex double bessel_buf[2*c->lMax+2]; // maximum order is 2n+1 for B coeffs, plus the zeroth.
double costheta = cos(kdlj.theta);
QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,2*c->lMax+1,
costheta,-1,legendre_buf));
QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(J, 2*c->lMax+1, kdlj.r, bessel_buf));
for (qpms_l_t l = 0; l <= 2*c->lMax+1; ++l)
for (qpms_m_t m = -l; m <= l; ++m) {
const qpms_y_t y_sc = qpms_mn2y_sc(m, l);
dest[y_sc] = bessel_buf[l]
* legendre_buf[gsl_sf_legendre_array_index(l, abs(m))]
* cexp(I * m * kdlj.phi);
}
return 0;
}
int qpms_trans_calculator_test_sswf_lc3p(const qpms_trans_calculator *c,
complex double *dest, const complex double k, const cart3_t r, const qpms_bessel_t J) {
csph_t kdlj = cart2csph(r);
kdlj.r *= k;
return qpms_trans_calculator_test_sswf(c, dest, kdlj, J);
}
#ifdef LATTICESUMS32
int qpms_trans_calculator_test_sswf_e32(const qpms_trans_calculator *c,
complex double * const sigmas_total, double * serr_total,
const double eta, const complex double k,
const cart2_t b1, const cart2_t b2,
const cart2_t beta,
const cart3_t particle_shift,
double maxR, double maxK,
const qpms_ewald_part parts)
{
const qpms_y_t nelem2_sc = qpms_lMax2nelem_sc(c->e3c->lMax);
//const qpms_y_t nelem = qpms_lMax2nelem(c->lMax);
const bool doerr = !!serr_total;
const bool do_sigma0 = ((particle_shift.x == 0) && (particle_shift.y == 0) && (particle_shift.z == 0)); // FIXME ignoring the case where particle_shift equals to lattice vector
complex double *sigmas_short = malloc(sizeof(complex double)*nelem2_sc);
complex double *sigmas_long = malloc(sizeof(complex double)*nelem2_sc);
//complex double *sigmas_total = malloc(sizeof(complex double)*nelem2_sc);
double *serr_short, *serr_long;//, *serr_total;
if(doerr) {
serr_short = malloc(sizeof(double)*nelem2_sc);
serr_long = malloc(sizeof(double)*nelem2_sc);
//serr_total = malloc(sizeof(double)*nelem2_sc);
} else serr_short = serr_long = serr_total = NULL;
const double unitcell_area = l2d_unitcell_area(b1, b2);
cart2_t rb1, rb2; // reciprocal basis
QPMS_ENSURE_SUCCESS(l2d_reciprocalBasis2pi(b1, b2, &rb1, &rb2));
if (parts & QPMS_EWALD_LONG_RANGE) {
PGen Kgen = PGen_xyWeb_new(rb1, rb2, BASIS_RTOL,
#ifdef GEN_KSHIFTEDPOINTS
beta,
#else
CART2_ZERO,
#endif
0, true, maxK, false);
QPMS_ENSURE_SUCCESS(ewald3_sigma_long(sigmas_long, serr_long, c->e3c, eta, k,
unitcell_area, LAT_2D_IN_3D_XYONLY, &Kgen,
#ifdef GEN_KSHIFTEDPOINTS
true,
#else
false,
#endif
cart22cart3xy(beta), particle_shift));
if(Kgen.stateData) // PGen not consumed entirely (converged earlier)
PGen_destroy(&Kgen);
}
if (parts & QPMS_EWALD_SHORT_RANGE) {
PGen Rgen = PGen_xyWeb_new(b1, b2, BASIS_RTOL,
#ifdef GEN_RSHIFTEDPOINTS
cart2_scale(-1 /*CHECKSIGN*/, cart3xy2cart2(particle_shift)),
#else
CART2_ZERO,
#endif
0, !do_sigma0, maxR, false);
#ifdef GEN_RSHIFTEDPOINTS // rather ugly hacks, LPTODO cleanup
if (particle_shift.z != 0) {
const cart3_t zshift = {0, 0, -particle_shift.z /*CHECKSIGN*/};
Rgen = Pgen_shifted_new(Rgen, zshift);
}
#endif
QPMS_ENSURE_SUCCESS(ewald3_sigma_short(sigmas_short, serr_short, c->e3c, eta, k,
particle_shift.z ? LAT_2D_IN_3D : LAT_2D_IN_3D_XYONLY, &Rgen,
#ifdef GEN_RSHIFTEDPOINTS
true,
#else
false,
#endif
cart22cart3xy(beta), particle_shift));
if(Rgen.stateData) // PGen not consumed entirely (converged earlier)
PGen_destroy(&Rgen);
}
for(qpms_y_t y = 0; y < nelem2_sc; ++y)
sigmas_total[y] = ((parts & QPMS_EWALD_SHORT_RANGE) ? sigmas_short[y] : 0)
+ ((parts & QPMS_EWALD_LONG_RANGE) ? sigmas_long[y] : 0);
if (doerr) for(qpms_y_t y = 0; y < nelem2_sc; ++y)
serr_total[y] = ((parts & QPMS_EWALD_SHORT_RANGE) ? serr_short[y] : 0)
+ ((parts & QPMS_EWALD_LONG_RANGE) ? serr_long[y] : 0);
complex double sigma0 = 0; double sigma0_err = 0;
if (do_sigma0 && (parts & QPMS_EWALD_0TERM)) {
QPMS_ENSURE_SUCCESS(ewald3_sigma0(&sigma0, &sigma0_err, c->e3c, eta, k));
const qpms_l_t y = qpms_mn2y_sc(0,0);
sigmas_total[y] += sigma0;
if(doerr) serr_total[y] += sigma0_err;
}
free(sigmas_short);
free(sigmas_long);
//free(sigmas_total);
if(doerr) {
free(serr_short);
free(serr_long);
//free(serr_total);
}
return 0;
}
#endif // LATTICESUMS32

28
qpms/translations_dbg.h Normal file
View File

@ -0,0 +1,28 @@
#ifndef TRANSLATIONS_DBG_H
// the following functions evaluate the SSWFs the same way as
// certain functions in translations.h, but they leave them
// as they are, without rearranging them into the translation
// operators.
#include "translations.h"
int qpms_trans_calculator_test_sswf(const qpms_trans_calculator *c,
_Complex double *dest, const csph_t kdlj, const qpms_bessel_t J);
int qpms_trans_calculator_test_sswf_lc3p(const qpms_trans_calculator *c,
_Complex double *dest, const _Complex double k, const cart3_t r, const qpms_bessel_t J);
#ifdef LATTICESUMS32
int qpms_trans_calculator_test_sswf_e32(const qpms_trans_calculator *c,
_Complex double * const sigmas_total, double * serr_total,
const double eta, const _Complex double k,
const cart2_t b1, const cart2_t b2,
const cart2_t beta,
const cart3_t particle_shift,
double maxR, double maxK,
const qpms_ewald_part parts);
#endif // LATTICESUMS32
#endif // TRANSLATIONS_DBG_H

View File

@ -6,12 +6,12 @@
/// Rearranges the default-ordered "A,B" array elements into "bspec"-defined matrix.
// TODO DOC
static inline void qpms_trans_array_from_AB(
complex double *t,
_Complex double *t,
const qpms_vswf_set_spec_t *const t_destspec,
const size_t t_deststride,
const qpms_vswf_set_spec_t *const t_srcspec,
const size_t t_srcstride,
const complex double *const A, const complex double *const B,
const _Complex double *const A, const _Complex double *const B,
/// A and B matrices' lMax.
/** This also determines their size and stride: they are assumed to
* be square matrices of size `nelem * nelem` where

View File

@ -1,4 +1,5 @@
#include "groups.h"
#include <complex.h>
/// Trivial group, with one (reduntant) generator.
/**

View File

@ -4,6 +4,7 @@
#ifndef VECTORS_H
#define VECTORS_H
#include <math.h>
#include <complex.h>
#ifndef M_PI_2
#define M_PI_2 (1.570796326794896619231321691639751442098584699687552910487)
#endif
@ -198,12 +199,12 @@ static inline double cart3_dist(const cart3_t a, const cart3_t b) {
return cart3norm(cart3_substract(a,b));
}
static inline bool cart3_isclose(const cart3_t a, const cart3_t b, double rtol, double atol) {
static inline _Bool cart3_isclose(const cart3_t a, const cart3_t b, double rtol, double atol) {
return cart3_dist(a,b) <= atol + rtol * (cart3norm(b) + cart3norm(a)) * .5;
}
/// Complex 3D vector scaling.
static inline ccart3_t ccart3_scale(const complex double c, const ccart3_t v) {
static inline ccart3_t ccart3_scale(const _Complex double c, const ccart3_t v) {
ccart3_t res = {c * v.x, c * v.y, c * v.z};
return res;
}
@ -221,7 +222,7 @@ static inline ccart3_t ccart3_substract(const ccart3_t a, const ccart3_t b) {
}
/// Complex 3D cartesian vector "dot product" without conjugation.
static inline complex double ccart3_dotnc(const ccart3_t a, const ccart3_t b) {
static inline _Complex double ccart3_dotnc(const ccart3_t a, const ccart3_t b) {
return a.x * b.x + a.y * b.y + a.z * b.z;
}
@ -244,13 +245,13 @@ static inline csphvec_t csphvec_substract(const csphvec_t a, const csphvec_t b)
}
/// Complex 3D vector (geographic coordinates) scaling.
static inline csphvec_t csphvec_scale(complex double c, const csphvec_t v) {
static inline csphvec_t csphvec_scale(_Complex double c, const csphvec_t v) {
csphvec_t res = {c * v.rc, c * v.thetac, c * v.phic};
return res;
}
/// Complex 3D vector (geographic coordinates) "dot product" without conjugation.
static inline complex double csphvec_dotnc(const csphvec_t a, const csphvec_t b) {
static inline _Complex double csphvec_dotnc(const csphvec_t a, const csphvec_t b) {
//N.B. no complex conjugation done here
return a.rc * b.rc + a.thetac * b.thetac + a.phic * b.phic;
}
@ -262,7 +263,7 @@ static inline sph_t sph_scale(double c, const sph_t s) {
}
/// "Complex spherical" coordinate system scaling.
static inline csph_t sph_cscale(complex double c, const sph_t s) {
static inline csph_t sph_cscale(_Complex double c, const sph_t s) {
csph_t res = {c * s.r, s.theta, s.phi};
return res;
}
@ -358,8 +359,9 @@ static inline csph_t ccart2csph(const ccart3_t cart) {
/// Real 3D cartesian to spherical (complex r) coordinates conversion. See @ref coord_conversions.
static inline csph_t cart2csph(const cart3_t cart) {
csph_t sph;
sph.r = cart3norm(cart);
sph.theta = sph.r ? acos(cart.z / sph.r) : M_PI_2;
const double r = cart3norm(cart);
sph.r = r;
sph.theta = sph.r ? acos(cart.z / r) : M_PI_2;
sph.phi = atan2(cart.y, cart.x);
return sph;
}
@ -578,7 +580,6 @@ static inline double anycoord2cart1(anycoord_point_t p, qpms_coord_system_t t) {
" to 1D not allowed.");
}
/// Coordinate conversion of point arrays (something to something).
/** The dest and src arrays must not overlap */
static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t tdest,
@ -589,31 +590,31 @@ static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t td
sph_t *d = (sph_t *) dest;
switch (tsrc & QPMS_COORDS_BITRANGE) {
case QPMS_COORDS_SPH: {
const sph_t *s = src;
const sph_t *s = (const sph_t*) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = s[i];
return;
} break;
case QPMS_COORDS_CART3: {
const cart3_t *s = src;
const cart3_t *s = (const cart3_t*) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = cart2sph(s[i]);
return;
} break;
case QPMS_COORDS_POL: {
const pol_t *s = src;
const pol_t *s = (const pol_t*) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = pol2sph_equator(s[i]);
return;
} break;
case QPMS_COORDS_CART2: {
const cart2_t *s = src;
const cart2_t *s = (const cart2_t*) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = cart22sph(s[i]);
return;
} break;
case QPMS_COORDS_CART1: {
const double *s = src;
const double *s = (const double *) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = cart12sph_zaxis(s[i]);
return;
@ -627,31 +628,31 @@ static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t td
cart3_t *d = (cart3_t *) dest;
switch (tsrc & QPMS_COORDS_BITRANGE) {
case QPMS_COORDS_SPH: {
const sph_t *s = src;
const sph_t *s = (const sph_t *) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = sph2cart(s[i]);
return;
} break;
case QPMS_COORDS_CART3: {
const cart3_t *s = src;
const cart3_t *s = (const cart3_t *) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = s[i];
return;
} break;
case QPMS_COORDS_POL: {
const pol_t *s = src;
const pol_t *s = (const pol_t *) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = pol2cart3_equator(s[i]);
return;
} break;
case QPMS_COORDS_CART2: {
const cart2_t *s = src;
const cart2_t *s = (const cart2_t *) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = cart22cart3xy(s[i]);
return;
} break;
case QPMS_COORDS_CART1: {
const double *s = src;
const double *s = (const double *)src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = cart12cart3z(s[i]);
return;
@ -669,13 +670,13 @@ static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t td
QPMS_PR_ERROR("Implicit conversion from 3D to 2D coordinates not allowed");
break;
case QPMS_COORDS_POL: {
const pol_t *s = src;
const pol_t *s = (const pol_t *) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = s[i];
return;
} break;
case QPMS_COORDS_CART2: {
const cart2_t *s = src;
const cart2_t *s = (const cart2_t *) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = cart2pol(s[i]);
return;
@ -696,13 +697,13 @@ static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t td
QPMS_PR_ERROR("Implicit conversion from 3D to 2D coordinates not allowed");
break;
case QPMS_COORDS_POL: {
const pol_t *s = src;
const pol_t *s = (const pol_t *) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = pol2cart(s[i]);
return;
} break;
case QPMS_COORDS_CART2: {
const cart2_t *s = src;
const cart2_t *s = (const cart2_t *) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = s[i];
return;
@ -727,7 +728,7 @@ static inline void qpms_array_coord_transform(void *dest, qpms_coord_system_t td
QPMS_PR_ERROR("Implicit conversion from 3D to 1D coordinates not allowed");
break;
case QPMS_COORDS_CART1: {
const double *s = src;
const double *s = (const double *) src;
for(size_t i = 0; i < nmemb; ++i)
d[i] = s[i];
return;
@ -899,6 +900,6 @@ static inline cart2_t cart2_from_double_array(const double a[]) {
typedef double matrix3d[3][3];
typedef double matrix2d[2][2];
typedef complex double cmatrix3d[3][3];
typedef complex double cmatrix2d[2][2];
typedef _Complex double cmatrix3d[3][3];
typedef _Complex double cmatrix2d[2][2];
#endif //VECTORS_H

View File

@ -1,2 +0,0 @@
#define QPMS_GIT_SHA1 "@GIT_SHA1@"
const char qpms_version_git_sha1[] = QPMS_GIT_SHA1;

View File

@ -1,4 +0,0 @@
#ifndef QPMS_VERSION_H
#define QPMS_VERSION_H
extern const char qpms_version_git_sha1[];
#endif //QPMS_VERSION_H

View File

@ -9,6 +9,7 @@
#include <string.h>
#include "qpms_error.h"
#include "normalisation.h"
#include <stdbool.h>
qpms_vswf_set_spec_t *qpms_vswf_set_spec_init() {
@ -52,7 +53,7 @@ qpms_errno_t qpms_vswf_set_spec_append(qpms_vswf_set_spec_t *s, const qpms_uvswf
return QPMS_SUCCESS;
}
bool qpms_vswf_set_spec_isidentical(const qpms_vswf_set_spec_t *a,
_Bool qpms_vswf_set_spec_isidentical(const qpms_vswf_set_spec_t *a,
const qpms_vswf_set_spec_t *b) {
if (a == b) return true;
if (a->norm != b->norm) return false;
@ -131,7 +132,7 @@ size_t *qpms_vswf_set_reindex(const qpms_vswf_set_spec_t *small, const qpms_vswf
while(big_pairs[bi].ui < small_pairs[si].ui)
++bi;
if(big_pairs[bi].ui == small_pairs[si].ui)
r[small_pairs[si].i_orig] = big_pairs[bi].i_orig;
r[small_pairs[si].i_orig] = big_pairs[si].i_orig;
else
r[small_pairs[si].i_orig] = ~(size_t)0;
}
@ -307,7 +308,7 @@ qpms_errno_t qpms_vswf_fill_alternative(csphvec_t *const longtarget, csphvec_t *
qpms_y_t nelem = qpms_lMax2nelem(lMax);
csphvec_t *a;
QPMS_CRASHING_MALLOC(a, 3*nelem*sizeof(csphvec_t))
csphvec_t * const a1 = a, * const a2 = a1 + nelem, * const a3 = a2 + 2 * nelem;
csphvec_t * const a1 = a, * const a2 = a1 + nelem, * const a3 = a1 + 2 * nelem;
QPMS_ENSURE_SUCCESS(qpms_vecspharm_fill(a1, a2, a3, lMax, kr, norm));
const csphvec_t *p1 = a1;
const csphvec_t *p2 = a2;
@ -478,7 +479,7 @@ qpms_errno_t qpms_planewave2vswf_fill_cart(cart3_t wavedir_cart /*allow complex
}
qpms_errno_t qpms_incfield_planewave(complex double *target, const qpms_vswf_set_spec_t *bspec,
const cart3_t evalpoint, const void *args, bool add) {
const cart3_t evalpoint, const void *args, _Bool add) {
QPMS_UNTESTED;
const qpms_incfield_planewave_params_t *p = args;

View File

@ -7,6 +7,10 @@
*/
#ifndef QPMS_VSWF_H
#define QPMS_VSWF_H
#ifdef __cplusplus
extern "C" {
#endif
#include <unistd.h> // ssize_t
#include "qpms_types.h"
#include <gsl/gsl_sf_legendre.h>
@ -16,7 +20,7 @@
/// Calculates the (regular VSWF) expansion coefficients of an external incident field.
typedef qpms_errno_t (*qpms_incfield_t)(
/// Target non-NULL array of the regular VSWF expansion coefficients of length bspec->n.
complex double *target,
_Complex double *target,
const qpms_vswf_set_spec_t *bspec,
const cart3_t evalpoint, ///< Point at which the VSWF expansion is made.
const void *args, ///< Pointer to additional function-specific arguments.
@ -82,7 +86,7 @@ qpms_errno_t qpms_uvswf_fill(
/** SVWF coefficients in \a coeffs must be ordered according to \a setspec->ilist.
*/
csphvec_t qpms_eval_uvswf(const qpms_vswf_set_spec_t *setspec,
const complex double *coeffs, ///< SVWF coefficient vector of size setspec->n.
const _Complex double *coeffs, ///< SVWF coefficient vector of size setspec->n.
csph_t kr, ///< Evaluation point.
qpms_bessel_t btyp);
@ -118,7 +122,7 @@ typedef struct qpms_incfield_planewave_params_t {
*/
qpms_errno_t qpms_incfield_planewave(
/// Target non-NULL array of the regular VSWF expansion coefficients of length bspec->n.
complex double *target,
_Complex double *target,
const qpms_vswf_set_spec_t *bspec,
const cart3_t evalpoint, ///< Point at which the VSWF expansion is made.
const void *args, ///< Pointer to additional function-specific arguments (converted to (const qpms_incfield_planewave_params_t *)).
@ -227,23 +231,26 @@ qpms_errno_t qpms_vecspharm_dual_fill(csphvec_t *const a1target, csphvec_t *cons
qpms_l_t lMax, sph_t dir, qpms_normalisation_t norm);
qpms_errno_t qpms_planewave2vswf_fill_cart(cart3_t wavedir, ccart3_t amplitude,
complex double *targt_longcoeff, complex double *target_mgcoeff, complex double *target_elcoeff,
_Complex double *targt_longcoeff, _Complex double *target_mgcoeff, _Complex double *target_elcoeff,
qpms_l_t lMax, qpms_normalisation_t norm);
qpms_errno_t qpms_planewave2vswf_fill_sph(sph_t wavedir, csphvec_t amplitude,
complex double *targt_longcoeff, complex double *target_mgcoeff, complex double *target_elcoeff,
_Complex double *targt_longcoeff, _Complex double *target_mgcoeff, _Complex double *target_elcoeff,
qpms_l_t lMax, qpms_normalisation_t norm);
csphvec_t qpms_eval_vswf(sph_t where,
complex double *longcoeffs, complex double *mgcoeffs, complex double *elcoeffs,
_Complex double *longcoeffs, _Complex double *mgcoeffs, _Complex double *elcoeffs,
qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm);
csphvec_t qpms_eval_vswf_csph(csph_t where,
complex double *longcoeffs, complex double *mgcoeffs, complex double *elcoeffs,
_Complex double *longcoeffs, _Complex double *mgcoeffs, _Complex double *elcoeffs,
qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm);
qpms_vswfset_sph_t *qpms_vswfset_make(qpms_l_t lMax, sph_t kdlj,
qpms_bessel_t btyp, qpms_normalisation_t norm);//NI
void qpms_vswfset_sph_pfree(qpms_vswfset_sph_t *);//NI
#ifdef __cplusplus
}
#endif
#endif // QPMS_VSWF_H

View File

@ -52,7 +52,7 @@ complex double qpms_wignerD_elem(const qpms_quat_t R,
}
complex double qpms_vswf_irot_elem_from_irot3(const qpms_irot3_t q,
const qpms_l_t l, const qpms_m_t mp, const qpms_m_t m, bool pseudo) {
const qpms_l_t l, const qpms_m_t mp, const qpms_m_t m, _Bool pseudo) {
#ifndef NDEBUG
qpms_irot3_checkdet(q);
#endif

View File

@ -111,8 +111,8 @@ cywaves = Extension('qpms.cywaves',
libraries = common_libs,
include_dirs = [numpy_includes],
)
cytranslations = Extension('qpms.cytranslations',
sources = ['qpms/cytranslations.pyx',
cytranslations = Extension('qpms.cytranslations_legacy',
sources = ['qpms/cytranslations_legacy.pyx',
'qpms/translations_python.c',
],
extra_compile_args=['-std=c99',
@ -150,7 +150,7 @@ qpms_c = Extension('qpms.qpms_c',
sources = ['qpms/qpms_c.pyx',],
libraries = common_libs,
include_dirs=['amos', 'qpms', numpy_includes],
extra_link_args=['-fopenmp'],
extra_link_args=['-fopenmp'], # needed?
extra_compile_args=['-fopenmp'],
)
@ -161,7 +161,7 @@ setup(name='qpms',
setup_requires=['cython>=0.28','numpy>=1.2','setuptools>=18.0'],
install_requires=['cython>=0.28',
#'quaternion','spherical_functions',
'sympy>=1.2', 'numpy>=1.2'],
'scipy>=0.18.0', 'sympy>=1.2', 'numpy>=1.2'],
#dependency_links=['https://github.com/moble/quaternion/archive/v2.0.tar.gz','https://github.com/moble/spherical_functions/archive/master.zip'],
ext_modules=cythonize([qpms_c, cywaves, cytranslations, cytmatrices, cycommon, cyquaternions, cybspec, cymaterials, cyewaldtest], include_path=['qpms', 'amos'], gdb_debug=True),
cmdclass = {'build_ext': build_ext},

View File

@ -0,0 +1,7 @@
find_package(Catch2 REQUIRED)
add_executable(catchtest all_includes.C t_vswf.C catch_aux.C)
# TODO ensure that the tests are linked against the "working tree" version of qpms,
# not the installed one
target_link_libraries(catchtest PRIVATE Catch2::Catch2WithMain qpms lapacke)

View File

@ -0,0 +1,30 @@
// Just to test that all headers are C++ compatible
#include <catch2/catch.hpp>
#include <vswf.h>
#include <beyn.h>
#include <gaunt.h>
#include <groups.h>
#include <kahansum.h>
#include <materials.h>
#include <optim.h>
#include <oshacks.h>
#include <parsing.h>
#include <qpms_error.h>
#include <qpms_specfunc.h>
#include <qpms_types.h>
#include <scatsystem.h>
#include <symmetries.h>
#include <tiny_inlines.h>
#include <tmatrices.h>
#include <tolerances.h>
#include <translations.h>
#include <vectors.h>
#include <pointgroups.h>
#include <lattices_types.h>
#include <indexing.h>
#include <ewald.h>
// C++ type conversion issues:
#include <quaternions.h>
//#include <lattices.h>
//#include <normalisation.h>
//#include <qpmsblas.h>

24
tests/catch/catch_aux.C Normal file
View File

@ -0,0 +1,24 @@
#include "catch_aux.h"
#include "complex.h"
namespace qpmstest{
/// Creates a new vector<double> from complex number array
std::vector<double> pointer2dvec(const _Complex double *arr, size_t siz){
std::vector<double> vec(2*siz);
for(size_t i = 0; i < siz; ++i) {
vec[2*i] = creal(arr[i]);
vec[2*i+1] = cimag(arr[i]);
}
return vec;
}
std::vector<double> pointer2dvec(const double *arr, size_t siz){
std::vector<double> vec(siz);
for(size_t i = 0; i < siz; ++i)
vec[i] = arr[i];
return vec;
}
std::vector<double> pointer2dvec(const csphvec_t *arr, size_t siz) {
return pointer2dvec(&arr->rc, 3*siz);
}
}

44
tests/catch/t_vswf.C Normal file
View File

@ -0,0 +1,44 @@
#include <catch2/catch.hpp>
#include <qpms_error.h>
#include <vswf.h>
#include "catch_aux.h"
using namespace qpmstest;
constexpr qpms_l_t lMax = 6;
static std::vector<double> vswfvals_dvector(double r, double theta, double phi,
qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm)
{
std::vector<double> vals(2 /*el and mg for now */ * (lMax * (lMax+2))
* 3 /* 3D vector components */ * 2 /*for real and imag*/);
qpms_vswf_fill(
nullptr, (csphvec_t *) &(vals[0]), (csphvec_t *) &(vals[(lMax*(lMax+2)) * 3 * 2]),
lMax, sph_t{r, theta, phi},
btyp, norm);
return vals;
}
static std::vector<double> vswfvals_dvector_alternative(double r, double theta, double phi,
qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm)
{
std::vector<double> vals(2 /*el and mg for now */ * (lMax * (lMax+2))
* 3 /* 3D vector components */ * 2 /*for real and imag*/);
qpms_vswf_fill_alternative(
nullptr, (csphvec_t *) &(vals[0]), (csphvec_t *) &(vals[(lMax*(lMax+2)) * 3 * 2]),
lMax, sph_t{r, theta, phi},
btyp, norm);
return vals;
}
TEST_CASE("VSWF alternative implementation") {
using Catch::Matchers::Approx;
auto norms = GENERATE(QPMS_NORMALISATION_DEFAULT, QPMS_NORMALISATION_NORM_POWER,
QPMS_NORMALISATION_NORM_SPHARM, QPMS_NORMALISATION_CONVENTION_SCUFF);
auto thetas = GENERATE(0., M_PI_2, M_PI_2/2, M_PI, 3*M_PI_2, 0.1, -0.2, 4., 12.666);
auto phis = GENERATE(0., M_PI_2, M_PI_2/2, M_PI, 3*M_PI_2, 0.1, -0.2, 4., 12.666);
auto rs_positive = GENERATE(0.0001, 0.001, 0.01, 0.1, 1., 10., 100.);
CHECK_THAT(
vswfvals_dvector(rs_positive, thetas, phis, lMax, QPMS_HANKEL_PLUS, norms),
Approx(vswfvals_dvector_alternative(rs_positive, thetas, phis, lMax, QPMS_HANKEL_PLUS, norms))
) ;
}

65
tests/test_ewald_consistency.py Executable file
View File

@ -0,0 +1,65 @@
#!/usr/bin/env python3
'''
Test 2D Ewald summation consistency for basis VSWFs, using varying values of eta.
'''
import sys
# Test parameters, rather arbitrary
npoints = 10
etafactors = [0.9, 1.1]
zpointmax = 0. if len(sys.argv) == 1 else float(sys.argv[1])
lattice_basis=[(580e-9,0),(0,580e-9)]
positions = [(0,0)] # particle positions
wavevector = [0.,0.,0.]
bgparam = 1.52
omega_eh = 2.1
lMax = 2
import numpy as np
from qpms import Particle, CTMatrix, lorentz_drude, EpsMuGenerator, TMatrixGenerator, EwaldPart, BaseSpec, FinitePointGroup, ScatteringSystem, TMatrixInterpolator, EpsMu, dbgmsg_enable, dbgmsg_disable, dbgmsg_active, BesselType,eV, hbar, c
from qpms.symmetries import point_group_info
eh = eV/hbar
np.random.seed(666)
omega = omega_eh * eh
bspec = BaseSpec(lMax=lMax)
medium = EpsMuGenerator(bgparam)
emg_dummy = EpsMuGenerator(10)
tmg_dummy = TMatrixGenerator.sphere(medium, emg_dummy, 50e-9)
particles = [Particle(pos, tmg_dummy, bspec) for pos in positions]
ss, ssw = ScatteringSystem.create(particles, medium, omega, latticebasis=lattice_basis)
wavevector = np.array(wavevector)
sswk = ssw._sswk(wavevector)
eta_orig = sswk.eta
print("default eta: %g" % (eta_orig,))
print("lattice basis:\n", ss.lattice_basis)
points = (np.random.rand(npoints) - .5)[:,None] * ss.lattice_basis[0] + (np.random.rand(npoints) - .5)[:,None] * ss.lattice_basis[1]
points += (np.random.rand(npoints) - .5)[:,None] * np.array([0.,0.,zpointmax])
fields_reference = sswk.scattered_field_basis(points)
fields = np.empty((len(etafactors),) + fields_reference.shape, dtype=complex)
fails = 0
for n, etafactor in enumerate(etafactors):
sswk.eta = eta_orig * etafactor
fields[n] = sswk.scattered_field_basis(points)
print(fields[n])
print(fields[n]-fields_reference)
if np.sum(1-np.isclose(fields[n], fields_reference)) > 0:
fails += 1
print(fails)
sys.exit(int(fails))