--- python-scattnlay-2.0.1.orig/README.md +++ python-scattnlay-2.0.1/README.md @@ -0,0 +1,128 @@ +![output example](/doc/OutputExample.png) +**Output example:** Field distribution inside layered Si\Ag\Si sphere +and Poynting vector distribution in Ag sphere with poweflow lines +calculated with Scattnlay. + +How to use scattnlay +==================== + +**Table of contents:** +- [Compile code](#compile-code) +- [Use](#use) +- [Papers](#papers) +- [Acknowledgment](#acknowledgment) +- [License](#license) + +Compile Code: +------------- +To compile the source you will need C++11 capable compiler. + +To compile the Python extension you also need the following packages: + - **python-numpy (>= 1.0)** + - **python-support (>= 0.6)** + - **python-all-dev (any version)** + - **python-numpy-dev (any version)** + +To compile the Debian package you also need the following packages: + - **debhelper (>=7.0.0)** + - **cdbs (>= 0.4.49)** + +Compilation options + - **make source** - Create source package for Python extension + - **make cython** - Convert Cython code to C++ + - **make python_ext** - Create Python extension using C++ code + - **make cython_ext** - Create Python extension using Cython code + - **make install** - Install Python extension on local system + - **make buildrpm** - Generate a rpm package for Python extension + - **make builddeb** - Generate a deb package for Python extension + - **make standalone** - Create standalone programs (scattnlay and fieldnlay) + - **make clean** - Delete temporal files + +Use: +---- + +1. Python library + * Use scattnlay directly + + ```python +from scattnlay import scattnlay, fieldnlay +... +x = ... +m = ... +coords = ... +terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(x, m) +terms, E, H = fieldnlay(x, m, coords) +... + ``` + + * Execute some of the test scripts (located in the folder 'tests/python') + Example: + + ```bash +./test01.py + ``` + +2. Standalone program + * Execute scattnlay directly: + + ```bash +scattnlay -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] [-t ti tf nt] [-c comment] + ``` + + * Execute fieldnlay directly: + + ```bash +fieldnlay -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] -p xi xf nx yi yf ny zi zf nz [-c comment] + ``` + + * Execute some of the test scripts (located in the folder 'tests/shell'): + + ```bash +./test01.sh > test01.csv + ``` + +3. C++ library + +```C++ + try { + MultiLayerMie multi_layer_mie; + multi_layer_mie.SetLayersSize(x); + multi_layer_mie.SetLayersIndex(m); + + multi_layer_mie.RunMieCalculation(); + + *Qsca = multi_layer_mie.GetQsca(); + *Qabs = multi_layer_mie.GetQabs(); + } catch(const std::invalid_argument& ia) { + // Will catch if multi_layer_mie fails or other errors. + std::cerr << "Invalid argument: " << ia.what() << std::endl; + throw std::invalid_argument(ia); + return -1; + } +``` + +Papers +------ + +1. "Scattering of electromagnetic radiation by a multilayered sphere" + O. Pena and U. Pal, Computer Physics Communications, vol. 180, + Nov. 2009, pp. 2348-2354. http://dx.doi.org/10.1016/j.cpc.2009.07.010 + +2. "Reduction of scattering using thin all-dielectric shells designed by stochastic optimizer" + Konstantin Ladutenko, Ovidio Peña-Rodríguez, Irina Melchakova, Ilya + Yagupov, and Pavel Belov J. Appl. Phys., vol. 116, pp. 184508, + 2014 http://dx.doi.org/10.1063/1.4900529 + +Acknowledgment +-------------- + +We expect that all publications describing work using this software, +or all commercial products using it, cite the following reference: +> O. Pena and U. Pal, "Scattering of electromagnetic radiation +> by a multilayered sphere," Computer Physics Communications, +> vol. 180, Nov. 2009, pp. 2348-2354. + +License +------- + +GPL v3+ --- python-scattnlay-2.0.1.orig/debian/source/options +++ python-scattnlay-2.0.1/debian/source/options @@ -0,0 +1,6 @@ +#https://www.debian.org/doc/manuals/maint-guide/dother.ru.html#sourceopt +# Don't store changes on autogenerated files +# extend-diff-ignore = "(^|/)(config\.sub|config\.guess|Makefile)$" + + +extend-diff-ignore = "\.png$" --- python-scattnlay-2.0.1.orig/scattnlay.pyx +++ python-scattnlay-2.0.1/scattnlay.pyx @@ -0,0 +1,150 @@ +# Copyright (C) 2009-2015 Ovidio Pena +# Copyright (C) 2013-2015 Konstantin Ladutenko +# +# This file is part of python-scattnlay +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# The only additional remark is that we expect that all publications +# describing work using this software, or all commercial products +# using it, cite the following reference: +# [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by +# a multilayered sphere," Computer Physics Communications, +# vol. 180, Nov. 2009, pp. 2348-2354. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +# distutils: language = c++ +# distutils: sources = nmie.cc + +from __future__ import division +import numpy as np +cimport numpy as np +from libcpp.vector cimport vector +from libcpp.vector cimport complex + +cdef inline double *npy2c(np.ndarray a): + assert a.dtype == np.float64 + + if not (a).flags["C_CONTIGUOUS"]: # Array is not contiguous, need to make contiguous copy + a = a.copy('C') + + # Return data pointer + return (a.data) + +cdef extern from "py_nmie.h": + cdef int ScattCoeffs(int L, int pl, vector[double] x, vector[complex] m, int nmax, double anr[], double ani[], double bnr[], double bni[]) + cdef int nMie(int L, int pl, vector[double] x, vector[complex] m, int nTheta, vector[double] Theta, int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, double S1r[], double S1i[], double S2r[], double S2i[]) + cdef int nField(int L, int pl, vector[double] x, vector[complex] m, int nmax, int nCoords, vector[double] Xp, vector[double] Yp, vector[double] Zp, double Erx[], double Ery[], double Erz[], double Eix[], double Eiy[], double Eiz[], double Hrx[], double Hry[], double Hrz[], double Hix[], double Hiy[], double Hiz[]) + +def scattcoeffs(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t, ndim = 2] m, np.int_t nmax, np.int_t pl = -1): + cdef Py_ssize_t i + + cdef np.ndarray[np.int_t, ndim = 1] terms = np.zeros(x.shape[0], dtype = np.int) + + cdef np.ndarray[np.complex128_t, ndim = 2] an = np.zeros((x.shape[0], nmax), dtype = np.complex128) + cdef np.ndarray[np.complex128_t, ndim = 2] bn = np.zeros((x.shape[0], nmax), dtype = np.complex128) + + cdef np.ndarray[np.float64_t, ndim = 1] anr + cdef np.ndarray[np.float64_t, ndim = 1] ani + cdef np.ndarray[np.float64_t, ndim = 1] bnr + cdef np.ndarray[np.float64_t, ndim = 1] bni + + for i in range(x.shape[0]): + anr = np.zeros(nmax, dtype = np.float64) + ani = np.zeros(nmax, dtype = np.float64) + bnr = np.zeros(nmax, dtype = np.float64) + bni = np.zeros(nmax, dtype = np.float64) + + terms[i] = ScattCoeffs(x.shape[1], pl, x[i].copy('C'), m[i].copy('C'), nmax, npy2c(anr), npy2c(ani), npy2c(bnr), npy2c(bni)) + + an[i] = anr.copy('C') + 1.0j*ani.copy('C') + bn[i] = bnr.copy('C') + 1.0j*bni.copy('C') + + return terms, an, bn + +def scattnlay(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t, ndim = 2] m, np.ndarray[np.float64_t, ndim = 1] theta = np.zeros(0, dtype = np.float64), np.int_t nmax = -1, np.int_t pl = -1): + cdef Py_ssize_t i + + cdef np.ndarray[np.int_t, ndim = 1] terms = np.zeros(x.shape[0], dtype = np.int) + + cdef np.ndarray[np.float64_t, ndim = 1] Qext = np.zeros(x.shape[0], dtype = np.float64) + cdef np.ndarray[np.float64_t, ndim = 1] Qabs = np.zeros(x.shape[0], dtype = np.float64) + cdef np.ndarray[np.float64_t, ndim = 1] Qsca = np.zeros(x.shape[0], dtype = np.float64) + cdef np.ndarray[np.float64_t, ndim = 1] Qbk = np.zeros(x.shape[0], dtype = np.float64) + cdef np.ndarray[np.float64_t, ndim = 1] Qpr = np.zeros(x.shape[0], dtype = np.float64) + cdef np.ndarray[np.float64_t, ndim = 1] g = np.zeros(x.shape[0], dtype = np.float64) + cdef np.ndarray[np.float64_t, ndim = 1] Albedo = np.zeros(x.shape[0], dtype = np.float64) + + cdef np.ndarray[np.complex128_t, ndim = 2] S1 = np.zeros((x.shape[0], theta.shape[0]), dtype = np.complex128) + cdef np.ndarray[np.complex128_t, ndim = 2] S2 = np.zeros((x.shape[0], theta.shape[0]), dtype = np.complex128) + + cdef np.ndarray[np.float64_t, ndim = 1] S1r + cdef np.ndarray[np.float64_t, ndim = 1] S1i + cdef np.ndarray[np.float64_t, ndim = 1] S2r + cdef np.ndarray[np.float64_t, ndim = 1] S2i + + for i in range(x.shape[0]): + S1r = np.zeros(theta.shape[0], dtype = np.float64) + S1i = np.zeros(theta.shape[0], dtype = np.float64) + S2r = np.zeros(theta.shape[0], dtype = np.float64) + S2i = np.zeros(theta.shape[0], dtype = np.float64) + + terms[i] = nMie(x.shape[1], pl, x[i].copy('C'), m[i].copy('C'), theta.shape[0], theta.copy('C'), nmax, &Qext[i], &Qsca[i], &Qabs[i], &Qbk[i], &Qpr[i], &g[i], &Albedo[i], npy2c(S1r), npy2c(S1i), npy2c(S2r), npy2c(S2i)) + + S1[i] = S1r.copy('C') + 1.0j*S1i.copy('C') + S2[i] = S2r.copy('C') + 1.0j*S2i.copy('C') + + return terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 + +def fieldnlay(np.ndarray[np.float64_t, ndim = 2] x, np.ndarray[np.complex128_t, ndim = 2] m, np.ndarray[np.float64_t, ndim = 2] coords, np.int_t nmax = -1, np.int_t pl = -1): + cdef Py_ssize_t i + + cdef np.ndarray[np.int_t, ndim = 1] terms = np.zeros(x.shape[0], dtype = np.int) + + cdef np.ndarray[np.complex128_t, ndim = 3] E = np.zeros((x.shape[0], coords.shape[0], 3), dtype = np.complex128) + cdef np.ndarray[np.complex128_t, ndim = 3] H = np.zeros((x.shape[0], coords.shape[0], 3), dtype = np.complex128) + + cdef np.ndarray[np.float64_t, ndim = 1] Erx + cdef np.ndarray[np.float64_t, ndim = 1] Ery + cdef np.ndarray[np.float64_t, ndim = 1] Erz + cdef np.ndarray[np.float64_t, ndim = 1] Eix + cdef np.ndarray[np.float64_t, ndim = 1] Eiy + cdef np.ndarray[np.float64_t, ndim = 1] Eiz + cdef np.ndarray[np.float64_t, ndim = 1] Hrx + cdef np.ndarray[np.float64_t, ndim = 1] Hry + cdef np.ndarray[np.float64_t, ndim = 1] Hrz + cdef np.ndarray[np.float64_t, ndim = 1] Hix + cdef np.ndarray[np.float64_t, ndim = 1] Hiy + cdef np.ndarray[np.float64_t, ndim = 1] Hiz + + for i in range(x.shape[0]): + Erx = np.zeros(coords.shape[0], dtype = np.float64) + Ery = np.zeros(coords.shape[0], dtype = np.float64) + Erz = np.zeros(coords.shape[0], dtype = np.float64) + Eix = np.zeros(coords.shape[0], dtype = np.float64) + Eiy = np.zeros(coords.shape[0], dtype = np.float64) + Eiz = np.zeros(coords.shape[0], dtype = np.float64) + Hrx = np.zeros(coords.shape[0], dtype = np.float64) + Hry = np.zeros(coords.shape[0], dtype = np.float64) + Hrz = np.zeros(coords.shape[0], dtype = np.float64) + Hix = np.zeros(coords.shape[0], dtype = np.float64) + Hiy = np.zeros(coords.shape[0], dtype = np.float64) + Hiz = np.zeros(coords.shape[0], dtype = np.float64) + + terms[i] = nField(x.shape[1], pl, x[i].copy('C'), m[i].copy('C'), nmax, coords.shape[0], coords[:, 0].copy('C'), coords[:, 1].copy('C'), coords[:, 2].copy('C'), npy2c(Erx), npy2c(Ery), npy2c(Erz), npy2c(Eix), npy2c(Eiy), npy2c(Eiz), npy2c(Hrx), npy2c(Hry), npy2c(Hrz), npy2c(Hix), npy2c(Hiy), npy2c(Hiz)) + + E[i] = np.vstack((Erx.copy('C') + 1.0j*Eix.copy('C'), Ery.copy('C') + 1.0j*Eiy.copy('C'), Erz.copy('C') + 1.0j*Eiz.copy('C'))).transpose() + H[i] = np.vstack((Hrx.copy('C') + 1.0j*Hix.copy('C'), Hry.copy('C') + 1.0j*Hiy.copy('C'), Hrz.copy('C') + 1.0j*Hiz.copy('C'))).transpose() + + return terms, E, H + --- python-scattnlay-2.0.1.orig/tests/c++/go-speed-test.sh +++ python-scattnlay-2.0.1/tests/c++/go-speed-test.sh @@ -0,0 +1,19 @@ +#!/bin/bash +path=$PWD +PROGRAM='scattnlay.bin' + +echo Compile with gcc +rm -f $PROGRAM + +#google profiler ######## FAST!!! +# file=speed-test.cc +# g++ -Ofast -std=c++11 $file ../../src/nmie.cc -lm -lrt -o $PROGRAM /usr/lib/libtcmalloc.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2 + +file=speed-test-applied.cc +g++ -Ofast -std=c++11 $file ../../src/nmie.cc ../../src/nmie-applied.cc -lm -lrt -o $PROGRAM /usr/lib/libtcmalloc.so.4 -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free -march=native -mtune=native -msse4.2 + +echo Should be: +echo test01, +1.41154e+00, +4.17695e-01, +9.93844e-01, +1.59427e-01, +1.25809e+00, +3.67376e-01, +2.95915e-01 +echo Result: +time ./$PROGRAM -l 5 0.4642 1.8000 1.7000 0.7114 0.8000 0.7000 0.7393 1.2000 0.0900 0.9168 2.8000 0.2000 1.0000 1.5000 0.4000 -t 0.0 90.0 5 -c test01 +# #result --- python-scattnlay-2.0.1.orig/tests/c++/speed-test.cc +++ python-scattnlay-2.0.1/tests/c++/speed-test.cc @@ -0,0 +1,237 @@ +//**********************************************************************************// +// Copyright (C) 2009-2015 Ovidio Pena // +// Copyright (C) 2013-2015 Konstantin Ladutenko // +// // +// This file is part of scattnlay // +// // +// This program is free software: you can redistribute it and/or modify // +// it under the terms of the GNU General Public License as published by // +// the Free Software Foundation, either version 3 of the License, or // +// (at your option) any later version. // +// // +// This program is distributed in the hope that it will be useful, // +// but WITHOUT ANY WARRANTY; without even the implied warranty of // +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // +// GNU General Public License for more details. // +// // +// The only additional remark is that we expect that all publications // +// describing work using this software, or all commercial products // +// using it, cite the following reference: // +// [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by // +// a multilayered sphere," Computer Physics Communications, // +// vol. 180, Nov. 2009, pp. 2348-2354. // +// // +// You should have received a copy of the GNU General Public License // +// along with this program. If not, see . // +//**********************************************************************************// + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +//sudo aptitude install libgoogle-perftools-dev +//#include +#include "../../src/nmie.h" + +timespec diff(timespec start, timespec end); +const double PI=3.14159265358979323846; +template inline T pow2(const T value) {return value*value;} + + +//***********************************************************************************// +// This is the main function of 'scattnlay', here we read the parameters as // +// arguments passed to the program which should be executed with the following // +// syntaxis: // +// ./scattnlay -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] [-t ti tf nt] [-c comment] // +// // +// When all the parameters were correctly passed we setup the integer L (the // +// number of layers) and the arrays x and m, containing the size parameters and // +// refractive indexes of the layers, respectively and call the function nMie. // +// If the calculation is successful the results are printed with the following // +// format: // +// // +// * If no comment was passed: // +// 'Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo' // +// // +// * If a comment was passed: // +// 'comment, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo' // +//***********************************************************************************// +int main(int argc, char *argv[]) { + try { + std::vector args; + args.assign(argv, argv + argc); + std::string error_msg(std::string("Insufficient parameters.\nUsage: ") + args[0] + + " -l Layers x1 m1.r m1.i [x2 m2.r m2.i ...] " + + "[-t ti tf nt] [-c comment]\n"); + enum mode_states {read_L, read_x, read_mr, read_mi, read_ti, read_tf, read_nt, read_comment}; + // for (auto arg : args) std::cout<< arg < x, Theta; + std::vector > m, S1, S2; + double Qext, Qabs, Qsca, Qbk, Qpr, g, Albedo; + + std::vector > mw, S1w, S2w; + double Qextw, Qabsw, Qscaw, Qbkw, Qprw, gw, Albedow; + + double ti = 0.0, tf = 90.0; + int nt = 0; + if (argc < 5) throw std::invalid_argument(error_msg); + + //strcpy(comment, ""); + // for (i = 1; i < argc; i++) { + int mode = -1; + double tmp_mr; + for (auto arg : args) { + // For each arg in args list we detect the change of the current + // read mode or read the arg. The reading args algorithm works + // as a finite-state machine. + + // Detecting new read mode (if it is a valid -key) + if (arg == "-l") { + mode = read_L; + continue; + } + if (arg == "-t") { + if ((mode != read_x) && (mode != read_comment)) + throw std::invalid_argument(std::string("Unfinished layer!\n") + +error_msg); + mode = read_ti; + continue; + } + if (arg == "-c") { + if ((mode != read_x) && (mode != read_nt)) + throw std::invalid_argument(std::string("Unfinished layer or theta!\n") + error_msg); + mode = read_comment; + continue; + } + // Reading data. For invalid date the exception will be thrown + // with the std:: and catched in the end. + if (mode == read_L) { + L = std::stoi(arg); + mode = read_x; + continue; + } + if (mode == read_x) { + x.push_back(std::stod(arg)); + mode = read_mr; + continue; + } + if (mode == read_mr) { + tmp_mr = std::stod(arg); + mode = read_mi; + continue; + } + if (mode == read_mi) { + m.push_back(std::complex( tmp_mr,std::stod(arg) )); + mode = read_x; + continue; + } + if (mode == read_ti) { + ti = std::stod(arg); + mode = read_tf; + continue; + } + if (mode == read_tf) { + tf = std::stod(arg); + mode = read_nt; + continue; + } + if (mode == read_nt) { + nt = std::stoi(arg); + Theta.resize(nt); + S1.resize(nt); + S2.resize(nt); + S1w.resize(nt); + S2w.resize(nt); + continue; + } + if (mode == read_comment) { + comment = arg; + has_comment = 1; + continue; + } + } + if ( (x.size() != m.size()) || (L != x.size()) ) + throw std::invalid_argument(std::string("Broken structure!\n") + +error_msg); + if ( (0 == m.size()) || ( 0 == x.size()) ) + throw std::invalid_argument(std::string("Empty structure!\n") + +error_msg); + + if (nt < 0) { + printf("Error reading Theta.\n"); + return -1; + } else if (nt == 1) { + Theta[0] = ti*PI/180.0; + } else { + for (i = 0; i < nt; i++) { + Theta[i] = (ti + (double)i*(tf - ti)/(nt - 1))*PI/180.0; + } + } + timespec time1, time2; + long cpptime_nsec, best_cpp; + long ctime_nsec, best_c; + long cpptime_sec, ctime_sec; + long repeats = 150; + //HeapProfilerStart("heapprof"); + do { + clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time1); + for (int i = 0; i 0) { + printf(" Theta, S1.r, S1.i, S2.r, S2.i\n"); + + for (i = 0; i < nt; i++) { + printf("%6.2f, %+.5e, %+.5e, %+.5e, %+.5e \n", Theta[i]*180.0/PI, S1w[i].real(), S1w[i].imag(), S2w[i].real(), S2w[i].imag()); + } + } + + } catch( const std::invalid_argument& ia ) { + // Will catch if multi_layer_mie fails or other errors. + std::cerr << "Invalid argument: " << ia.what() << std::endl; + return -1; + } + return 0; +} + + + +timespec diff(timespec start, timespec end) +{ + timespec temp; + if ((end.tv_nsec-start.tv_nsec)<0) { + temp.tv_sec = end.tv_sec-start.tv_sec-1; + temp.tv_nsec = 1000000000+end.tv_nsec-start.tv_nsec; + } else { + temp.tv_sec = end.tv_sec-start.tv_sec; + temp.tv_nsec = end.tv_nsec-start.tv_nsec; + } + return temp; +} --- python-scattnlay-2.0.1.orig/tests/c++/test-on-i5-4300U.log +++ python-scattnlay-2.0.1/tests/c++/test-on-i5-4300U.log @@ -0,0 +1,15 @@ +Compile with gcc +Should be: +test01, +1.41154e+00, +4.17695e-01, +9.93844e-01, +1.59427e-01, +1.25809e+00, +3.67376e-01, +2.95915e-01 +Result: +-- C++ time consumed 0.00700425 sec +-- C++ time consumed 0.0670615 sec +-- C++ time consumed 0.669772 sec + +test01, +1.41154e+00, +4.17695e-01, +9.93844e-01, +1.59427e-01, +1.25809e+00, +3.67376e-01, +2.95915e-01 + Theta, S1.r, S1.i, S2.r, S2.i + 0.00, +3.52885e-01, -4.27839e-01, +3.52885e-01, -4.27839e-01 + 22.50, +3.44554e-01, -4.14362e-01, +3.29828e-01, -3.94049e-01 + 45.00, +3.21195e-01, -3.76903e-01, +2.65943e-01, -3.03992e-01 + 67.50, +2.87219e-01, -3.23325e-01, +1.75076e-01, -1.85579e-01 + 90.00, +2.48609e-01, -2.63824e-01, +7.48392e-02, -6.95864e-02 --- python-scattnlay-2.0.1.orig/utils/bessel/bessel.cc +++ python-scattnlay-2.0.1/utils/bessel/bessel.cc @@ -0,0 +1,354 @@ +//**********************************************************************************// +// Copyright (C) 2009-2015 Ovidio Pena // +// Copyright (C) 2013-2015 Konstantin Ladutenko // +// // +// This file is part of scattnlay // +// // +// This program is free software: you can redistribute it and/or modify // +// it under the terms of the GNU General Public License as published by // +// the Free Software Foundation, either version 3 of the License, or // +// (at your option) any later version. // +// // +// This program is distributed in the hope that it will be useful, // +// but WITHOUT ANY WARRANTY; without even the implied warranty of // +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // +// GNU General Public License for more details. // +// // +// The only additional remark is that we expect that all publications // +// describing work using this software, or all commercial products // +// using it, cite the following reference: // +// [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by // +// a multilayered sphere," Computer Physics Communications, // +// vol. 180, Nov. 2009, pp. 2348-2354. // +// // +// You should have received a copy of the GNU General Public License // +// along with this program. If not, see . // +//**********************************************************************************// +#include "bessel.h" +#include +#include +#include +#include +#include +#include + +namespace nmie { + namespace bessel { + + void calcZeta(int n, std::complexz, std::vector< std::complex >& Zeta, + std::vector< std::complex >& dZeta) { + std::vector< std::complex > csj, cdj, csy, cdy; + int nm; + csphjy (n, z, nm, csj, cdj, csy, cdy ); + Zeta.resize(n+1); + dZeta.resize(n+1); + auto c_i = std::complex(0.0,1.0); + for (int i = 0; i < n+1; ++i) { + Zeta[i]=z*(csj[i] + c_i*csy[i]); + dZeta[i]=z*(cdj[i] + c_i*cdy[i])+csj[i] + c_i*csy[i]; + } + } // end of calcZeta() + + void calcPsi(int n, std::complexz, std::vector< std::complex >& Psi, + std::vector< std::complex >& dPsi) { + std::vector< std::complex > csj, cdj, csy, cdy; + int nm; + csphjy (n, z, nm, csj, cdj, csy, cdy ); + Psi.resize(n+1); + dPsi.resize(n+1); + for (int i = 0; i < n+1; ++i) { + Psi[i]=z*(csj[i]); + dPsi[i]=z*(cdj[i])+csj[i]; + } + } // end of calcPsi() + +// !*****************************************************************************80 +// +// C++ port of fortran code +// +// !! CSPHJY: spherical Bessel functions jn(z) and yn(z) for complex argument. +// ! +// ! Discussion: +// ! +// ! This procedure computes spherical Bessel functions jn(z) and yn(z) +// ! and their derivatives for a complex argument. +// ! +// ! Licensing: +// ! +// ! This routine is copyrighted by Shanjie Zhang and Jianming Jin. However, +// ! they give permission to incorporate this routine into a user program +// ! provided that the copyright is acknowledged. +// ! +// ! Modified: +// ! +// ! 01 August 2012 +// ! +// ! Author: +// ! +// ! Shanjie Zhang, Jianming Jin +// ! +// ! Reference: +// ! +// ! Shanjie Zhang, Jianming Jin, +// ! Computation of Special Functions, +// ! Wiley, 1996, +// ! ISBN: 0-471-11963-6, +// ! LC: QA351.C45. +// ! +// ! Parameters: +// ! +// ! Input, integer ( kind = 4 ) N, the order of jn(z) and yn(z). +// ! +// ! Input, complex ( kind = 8 ) Z, the argument. +// ! +// ! Output, integer ( kind = 4 ) NM, the highest order computed. +// ! +// ! Output, complex ( kind = 8 ) CSJ(0:N0, CDJ(0:N), CSY(0:N), CDY(0:N), +// ! the values of jn(z), jn'(z), yn(z), yn'(z). +// ! + void csphjy (int n, std::complexz, int& nm, + std::vector< std::complex >& csj, + std::vector< std::complex >& cdj, + std::vector< std::complex >& csy, + std::vector< std::complex >& cdy ) { + double a0; + csj.resize(n+1); + cdj.resize(n+1); + csy.resize(n+1); + cdy.resize(n+1); + std::complex cf, cf0, cf1, cs, csa, csb; + int m; + a0 = std::abs(z); + nm = n; + if (a0 < 1.0e-60) { + for (int k = 0; k < n+1; ++k) { + csj[k] = 0.0; + cdj[k] = 0.0; + csy[k] = -1.0e+300; + cdy[k] = 1.0e+300; + } + csj[0] = std::complex( 1.0, 0.0); + cdj[1] = std::complex( 0.3333333333333333, 0.0); + return; + } + csj[0] = std::sin ( z ) / z; + csj[1] = ( csj[0] - std::cos ( z ) ) / z; + + if ( 2 <= n ) { + csa = csj[0]; + csb = csj[1]; + int precision = 1; + m = msta1 ( a0, 200*precision); + if ( m < n ) nm = m; + else m = msta2 ( a0, n, 15*precision); + cf0 = 0.0; + cf1 = 1.0e-100; + for (int k = m; k>=0; --k) { + cf = ( 2.0 * k + 3.0 ) * cf1 / z - cf0; + if ( k <= nm ) csj[k] = cf; + cf0 = cf1; + cf1 = cf; + } + if ( std::abs ( csa ) <= std::abs ( csb ) ) cs = csb / cf0; + else cs = csa / cf; + for (int k = 0; k <= nm; ++k) { + csj[k] = cs * csj[k]; + } + } + cdj[0] = ( std::cos ( z ) - std::sin ( z ) / z ) / z; + for (int k = 1; k <=nm; ++k) { + cdj[k] = csj[k-1] - ( k + 1.0 ) * csj[k] / z; + } + csy[0] = - std::cos ( z ) / z; + csy[1] = ( csy[0] - std::sin ( z ) ) / z; + cdy[0] = ( std::sin ( z ) + std::cos ( z ) / z ) / z; + cdy[1] = ( 2.0 * cdy[0] - std::cos ( z ) ) / z; + for (int k = 2; k<=nm; ++k) { + if ( std::abs ( csj[k-2] ) < std::abs ( csj[k-1] ) ) { + csy[k] = ( csj[k] * csy[k-1] - 1.0 / ( z * z ) ) / csj[k-1]; + } else { + csy[k] = ( csj[k] * csy[k-2] - ( 2.0 * k - 1.0 ) / std::pow(z,3) ) + / csj[k-2]; + } + } + for (int k = 2; k<=nm; ++k) { + cdy[k] = csy[k-1] - ( k + 1.0 ) * csy[k] / z; + } + + return; + } + // function msta2 ( x, n, mp ) + + // !*****************************************************************************80 + // ! + // !! MSTA2 determines a backward recurrence starting point for Jn(x). + // ! + // ! Discussion: + // ! + // ! This procedure determines the starting point for a backward + // ! recurrence such that all Jn(x) has MP significant digits. + // ! + // ! Licensing: + // ! + // ! This routine is copyrighted by Shanjie Zhang and Jianming Jin. However, + // ! they give permission to incorporate this routine into a user program + // ! provided that the copyright is acknowledged. + // ! + // ! Modified: + // ! + // ! 08 July 2012 + // ! + // ! Author: + // ! + // ! Shanjie Zhang, Jianming Jin + // ! + // ! Reference: + // ! + // ! Shanjie Zhang, Jianming Jin, + // ! Computation of Special Functions, + // ! Wiley, 1996, + // ! ISBN: 0-471-11963-6, + // ! LC: QA351.C45. + // ! + // ! Parameters: + // ! + // ! Input, real ( kind = 8 ) X, the argument of Jn(x). + // ! + // ! Input, integer ( kind = 4 ) N, the order of Jn(x). + // ! + // ! Input, integer ( kind = 4 ) MP, the number of significant digits. + // ! + // ! Output, integer ( kind = 4 ) MSTA2, the starting point. + // ! + int msta2 ( double x, int n, int mp ) { + double a0, ejn, f, f0, f1, hmp; + int n0, n1, nn; + double obj; + a0 = std::abs ( x ); + hmp = 0.5 * mp; + ejn = envj ( n, a0 ); + if ( ejn <= hmp ) { + obj = mp; + n0 = static_cast ( 1.1 * a0 ); + } else { + obj = hmp + ejn; + n0 = n; + } + f0 = envj ( n0, a0 ) - obj; + n1 = n0 + 5; + f1 = envj ( n1, a0 ) - obj; + for (int it = 1; it < 21; ++it) { + nn = n1 - ( n1 - n0 ) / ( 1.0 - f0 / f1 ); + f = envj ( nn, a0 ) - obj; + if ( std::abs ( nn - n1 ) < 1 ) break; + n0 = n1; + f0 = f1; + n1 = nn; + f1 = f; + } + return nn + 10; + } + + + +// !*****************************************************************************80 +// ! +// !! MSTA1 determines a backward recurrence starting point for Jn(x). +// ! +// ! Discussion: +// ! +// ! This procedure determines the starting point for backward +// ! recurrence such that the magnitude of +// ! Jn(x) at that point is about 10^(-MP). +// ! +// ! Licensing: +// ! +// ! This routine is copyrighted by Shanjie Zhang and Jianming Jin. However, +// ! they give permission to incorporate this routine into a user program +// ! provided that the copyright is acknowledged. +// ! +// ! Modified: +// ! +// ! 08 July 2012 +// ! +// ! Author: +// ! +// ! Shanjie Zhang, Jianming Jin +// ! +// ! Reference: +// ! +// ! Shanjie Zhang, Jianming Jin, +// ! Computation of Special Functions, +// ! Wiley, 1996, +// ! ISBN: 0-471-11963-6, +// ! LC: QA351.C45. +// ! +// ! Parameters: +// ! +// ! Input, real ( kind = 8 ) X, the argument. +// ! +// ! Input, integer ( kind = 4 ) MP, the negative logarithm of the +// ! desired magnitude. +// ! +// ! Output, integer ( kind = 4 ) MSTA1, the starting point. +// ! + int msta1 ( double x, int mp ) { + double a0, f, f0, f1; + int n0, n1, nn; + a0 = std::abs ( x ); + n0 = static_cast (1.1 * a0 ) + 1; + f0 = envj ( n0, a0 ) - mp; + n1 = n0 + 5; + f1 = envj ( n1, a0 ) - mp; + for (int it = 1; it <= 20; ++it) { + nn = n1 - ( n1 - n0 ) / ( 1.0 - f0 / f1 ); + f = envj ( nn, a0 ) - mp; + if ( abs ( nn - n1 ) < 1 ) break; + n0 = n1; + f0 = f1; + n1 = nn; + f1 = f; + } + return nn; + } + // function envj ( n, x ) + + // !*****************************************************************************80 + // ! + // !! ENVJ is a utility function used by MSTA1 and MSTA2. + // ! + // ! Licensing: + // ! + // ! This routine is copyrighted by Shanjie Zhang and Jianming Jin. However, + // ! they give permission to incorporate this routine into a user program + // ! provided that the copyright is acknowledged. + // ! + // ! Modified: + // ! + // ! 14 March 2012 + // ! + // ! Author: + // ! + // ! Shanjie Zhang, Jianming Jin + // ! + // ! Reference: + // ! + // ! Shanjie Zhang, Jianming Jin, + // ! Computation of Special Functions, + // ! Wiley, 1996, + // ! ISBN: 0-471-11963-6, + // ! LC: QA351.C45. + // ! + // ! Parameters: + // ! + // ! Input, integer ( kind = 4 ) N, ? + // ! + // ! Input, real ( kind = 8 ) X, ? + // ! + // ! Output, real ( kind = 8 ) ENVJ, ? + // ! + double envj (int n, double x ) { + return 0.5 * std::log10(6.28 * n) - n * std::log10(1.36 * x / static_cast(n) ); + } + } // end of namespace bessel +} // end of namespace nmie --- python-scattnlay-2.0.1.orig/utils/bessel/bessel.h +++ python-scattnlay-2.0.1/utils/bessel/bessel.h @@ -0,0 +1,46 @@ +//**********************************************************************************// +// Copyright (C) 2009-2015 Ovidio Pena // +// Copyright (C) 2013-2015 Konstantin Ladutenko // +// // +// This file is part of scattnlay // +// // +// This program is free software: you can redistribute it and/or modify // +// it under the terms of the GNU General Public License as published by // +// the Free Software Foundation, either version 3 of the License, or // +// (at your option) any later version. // +// // +// This program is distributed in the hope that it will be useful, // +// but WITHOUT ANY WARRANTY; without even the implied warranty of // +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // +// GNU General Public License for more details. // +// // +// The only additional remark is that we expect that all publications // +// describing work using this software, or all commercial products // +// using it, cite the following reference: // +// [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by // +// a multilayered sphere," Computer Physics Communications, // +// vol. 180, Nov. 2009, pp. 2348-2354. // +// // +// You should have received a copy of the GNU General Public License // +// along with this program. If not, see . // +//**********************************************************************************// +#include +#include +namespace nmie { + namespace bessel { + void calcZeta(int n, std::complexz, + std::vector< std::complex >& Zeta, + std::vector< std::complex >& dZeta); + void calcPsi(int n, std::complexz, + std::vector< std::complex >& Psi, + std::vector< std::complex >& dPsi); + void csphjy (int n, std::complexz, int& nm, + std::vector< std::complex >& csj, + std::vector< std::complex >& cdj, + std::vector< std::complex >& csy, + std::vector< std::complex >& cdy ); + double envj (int n, double x ); + int msta2 ( double x, int n, int mp ); + int msta1 ( double x, int mp ); + } // end of namespace bessel +} // end of namespace nmie --- python-scattnlay-2.0.1.orig/utils/bessel/go.sh +++ python-scattnlay-2.0.1/utils/bessel/go.sh @@ -0,0 +1,3 @@ +#!/bin/bash +clang++ -g -O1 -fsanitize=address -fno-optimize-sibling-calls -fno-omit-frame-pointer -std=c++11 ../../bessel.cc ../../nmie.cc test_bessel.cc -lm -lrt -o test_bessel.bin +./test_bessel.bin --- python-scattnlay-2.0.1.orig/utils/bessel/test_bessel.cc +++ python-scattnlay-2.0.1/utils/bessel/test_bessel.cc @@ -0,0 +1,428 @@ +//**********************************************************************************// +// Copyright (C) 2009-2015 Ovidio Pena // +// Copyright (C) 2013-2015 Konstantin Ladutenko // +// // +// This file is part of scattnlay // +// // +// This program is free software: you can redistribute it and/or modify // +// it under the terms of the GNU General Public License as published by // +// the Free Software Foundation, either version 3 of the License, or // +// (at your option) any later version. // +// // +// This program is distributed in the hope that it will be useful, // +// but WITHOUT ANY WARRANTY; without even the implied warranty of // +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // +// GNU General Public License for more details. // +// // +// The only additional remark is that we expect that all publications // +// describing work using this software, or all commercial products // +// using it, cite the following reference: // +// [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by // +// a multilayered sphere," Computer Physics Communications, // +// vol. 180, Nov. 2009, pp. 2348-2354. // +// // +// You should have received a copy of the GNU General Public License // +// along with this program. If not, see . // +//**********************************************************************************// +#include +#include +#include +#include +#include + +#include "../../bessel.h" +#include "../../nmie.h" + +int nmax_ = -1; +//Temporary variables +std::vector > PsiZeta_; + +void calcD1D3(const std::complex z, std::vector >& D1, std::vector >& D3) { + D1[nmax_] = std::complex(0.0, 0.0); + const std::complex zinv = std::complex(1.0, 0.0)/z; + for (int n = nmax_; n > 0; n--) { + D1[n - 1] = static_cast(n)*zinv - 1.0/(D1[n] + static_cast(n)*zinv); + } + if (std::abs(D1[0]) > 100000.0) + throw std::invalid_argument("Unstable D1! Please, try to change input parameters!\n"); + + PsiZeta_.resize(nmax_ + 1); + + // Upward recurrence for PsiZeta and D3 - equations (18a) - (18d) + PsiZeta_[0] = 0.5*(1.0 - std::complex(std::cos(2.0*z.real()), std::sin(2.0*z.real())) + *std::exp(-2.0*z.imag())); + D3[0] = std::complex(0.0, 1.0); + for (int n = 1; n <= nmax_; n++) { + PsiZeta_[n] = PsiZeta_[n - 1]*(static_cast(n)*zinv - D1[n - 1]) + *(static_cast(n)*zinv - D3[n - 1]); + D3[n] = D1[n] + std::complex(0.0, 1.0)/PsiZeta_[n]; + } + } + + void calcPsiZeta(std::complex z, std::vector >& Psi, std::vector >& Zeta) { + std::complex c_i(0.0, 1.0); + std::vector > D1(nmax_ + 1), D3(nmax_ + 1); + calcD1D3(z, D1, D3); + Psi.resize(nmax_+1); + Zeta.resize(nmax_+1); + + Psi[0] = std::sin(z); + Zeta[0] = std::sin(z) - c_i*std::cos(z); + for (int n = 1; n <= nmax_; n++) { + Psi[n] = Psi[n - 1]*(static_cast(n)/z - D1[n - 1]); + Zeta[n] = Zeta[n - 1]*(static_cast(n)/z - D3[n - 1]); + } + } + + +int main(int argc, char *argv[]) { + + int n_small = 1, n_big = 15; + int n = n_big+2, n_print; + std::complex z_small(0.0012, 0.0034); + std::complex z_medium(1.0, 2.0); + std::complex z_big(18.0, 17.0); + std::complex z_big_real(81.0, 0.0); + std::complex z_big_imag(0.0, 13.0); + std::complex z; + + std::vector > D1, D3; + + int nm; + // csj and cdj - complex argument spherical bessel first kind and derivative + // csy and cdy - complex argument spherical bessel second kind and derivative + std::vector< std::complex > csj, cdj, csy, cdy, csh; + std::vector > Psi, Zeta, dPsi, dZeta; + + auto c_i = std::complex(0.0,1.0); + + printf("===== Small ========\n"); + z = z_small; + printf("z = %+10.5er%+10.5ei\n", z.real(), z.imag()); + nmie::bessel::csphjy (n, z, nm, csj, cdj, csy, cdy ); + csh.resize(csj.size()); + for (unsigned int i = 0; i < csj.size(); ++i) + csh[i] = csj[i]+c_i*csy[i]; + nmie::bessel::calcPsi(n, z, Psi, dPsi); + nmie::bessel::calcZeta(n, z, Zeta, dZeta); + + // n_print = n_small; + // printf("csj[%i] = %+10.5er%+10.5ei; csy[%i] = %+10.5er%+10.5ei\n", n_print, + // csj[n_print].real(), csj[n_print].imag(), n_print, + // csy[n_print].real(), csy[n_print].imag()); + // printf("WA csj = e r e i; WA csy = e r e i\n"); + // n_print = n_big; + // printf("csj[%i] = %+10.5er%+10.5ei; csy[%i] = %+10.5er%+10.5ei\n", n_big, + // csj[n_print].real(), csj[n_print].imag(), n_print, + // csy[n_print].real(), csy[n_print].imag()); + // printf(" WA csj = e r e i; WA csy = e r e i\n"); + + // n_print = n_small; + // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print, + // csh[n_print].real(), csh[n_print].imag()); + // printf("WA csj = e r e i;\n"); + // n_print = n_big; + // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print, + // csh[n_print].real(), csh[n_print].imag()); + // printf("WA csj = e r e i;\n"); + + n_print = n_small; + printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_print, + Psi[n_print].real(), Psi[n_print].imag(), n_print, + Zeta[n_print].real(), Zeta[n_print].imag()); + printf("WA Psi = e r e i; WA Zeta = e r e i\n"); + n_print = n_big; + printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_big, + Psi[n_print].real(), Psi[n_print].imag(), n_print, + Zeta[n_print].real(), Zeta[n_print].imag()); + printf(" WA Psi = e r e i; WA Zeta = e r e i\n"); + + n_print = n_small; + printf("dPsi[%i] = %+10.5er%+10.5ei; dZeta[%i] = %+10.5er%+10.5ei\n", n_print, + dPsi[n_print].real(), dPsi[n_print].imag(), n_print, + dZeta[n_print].real(), dZeta[n_print].imag()); + printf("WA dPsi = 0.800005e-03r+0.226667e-02i; WA dZeta = +4.8284 e+04r-5.98822e+04i\n"); + n_print = n_big; + printf("dPsi[%i] = %+10.5er%+10.5ei; dZeta[%i] = %+10.5er%+10.5ei\n", n_big, + dPsi[n_print].real(), dPsi[n_print].imag(), n_print, + dZeta[n_print].real(), dZeta[n_print].imag()); + printf(" WA dPsi = +1.75388e-53r-6.94429e-54i; WA dZeta = +8.58553e+55r+7.47397e+55i\n"); + + nmax_ = nm; + printf("----- Scattnlay (nmax_=%d)-----\n",nmax_); + calcPsiZeta(z, Psi, Zeta); + + n_print = n_small; + printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_print, + Psi[n_print].real(), Psi[n_print].imag(), n_print, + Zeta[n_print].real(), Zeta[n_print].imag()); + printf("WA Psi = e r e i; WA Zeta = e r e i\n"); + n_print = n_big; + printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_big, + Psi[n_print].real(), Psi[n_print].imag(), n_print, + Zeta[n_print].real(), Zeta[n_print].imag()); + printf(" WA Psi = e r e i; WA Zeta = e r e i\n"); + + + + + printf("===== Medium ========\n"); + z = z_medium; + printf("z = %+10.5er%+10.5ei\n", z.real(), z.imag()); + nmie::bessel::csphjy (n, z, nm, csj, cdj, csy, cdy ); + csh.resize(csj.size()); + for (unsigned int i = 0; i < csj.size(); ++i) + csh[i] = csj[i]+c_i*csy[i]; + nmie::bessel::calcPsi(n, z, Psi, dPsi); + nmie::bessel::calcZeta(n, z, Zeta, dZeta); + + // n_print = n_small; + // printf("csj[%i] = %+10.5er%+10.5ei; csy[%i] = %+10.5er%+10.5ei\n", n_print, + // csj[n_print].real(), csj[n_print].imag(), n_print, + // csy[n_print].real(), csy[n_print].imag()); + // printf("WA csj = e r e i; WA csy = e r e i\n"); + // n_print = n_big; + // printf("csj[%i] = %+10.5er%+10.5ei; csy[%i] = %+10.5er%+10.5ei\n", n_big, + // csj[n_print].real(), csj[n_print].imag(), n_print, + // csy[n_print].real(), csy[n_print].imag()); + // printf(" WA csj = e r e i; WA csy = e r e i\n"); + + // n_print = n_small; + // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print, + // csh[n_print].real(), csh[n_print].imag()); + // printf("WA csj = e r e i;\n"); + // n_print = n_big; + // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print, + // csh[n_print].real(), csh[n_print].imag()); + // printf("WA csj = e r e i;\n"); + + n_print = n_small; + printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_print, + Psi[n_print].real(), Psi[n_print].imag(), n_print, + Zeta[n_print].real(), Zeta[n_print].imag()); + printf("WA Psi = e r e i; WA Zeta = e r e i\n"); + n_print = n_big; + printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_big, + Psi[n_print].real(), Psi[n_print].imag(), n_print, + Zeta[n_print].real(), Zeta[n_print].imag()); + printf(" WA Psi = e r e i; WA Zeta = e r e i\n"); + + n_print = n_small; + printf("dPsi[%i] = %+10.5er%+10.5ei; dZeta[%i] = %+10.5er%+10.5ei\n", n_print, + dPsi[n_print].real(), dPsi[n_print].imag(), n_print, + dZeta[n_print].real(), dZeta[n_print].imag()); + printf("WA dPsi = 2.41792e r+1.27781 e i; WA dZeta = +1.99423e-01r-7.01483e-02i\n"); + n_print = n_big; + printf("dPsi[%i] = %+10.5er%+10.5ei; dZeta[%i] = %+10.5er%+10.5ei\n", n_big, + dPsi[n_print].real(), dPsi[n_print].imag(), n_print, + dZeta[n_print].real(), dZeta[n_print].imag()); + printf(" WA dPsi = -1.03313e-11r-1.13267e-11i; WA dZeta = -2.11402e+11r+8.34528e+10i\n"); + + nmax_ = nm; + printf("----- Scattnlay (nmax_=%d)-----\n",nmax_); + calcPsiZeta(z, Psi, Zeta); + + n_print = n_small; + printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_print, + Psi[n_print].real(), Psi[n_print].imag(), n_print, + Zeta[n_print].real(), Zeta[n_print].imag()); + printf("WA Psi = e r e i; WA Zeta = e r e i\n"); + n_print = n_big; + printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_big, + Psi[n_print].real(), Psi[n_print].imag(), n_print, + Zeta[n_print].real(), Zeta[n_print].imag()); + printf(" WA Psi = e r e i; WA Zeta = e r e i\n"); + + printf("===== Big ========\n"); + z = z_big; + printf("z = %+10.5er%+10.5ei\n", z.real(), z.imag()); + nmie::bessel::csphjy (n, z, nm, csj, cdj, csy, cdy ); + csh.resize(csj.size()); + for (unsigned int i = 0; i < csj.size(); ++i) + csh[i] = csj[i]+c_i*csy[i]; + nmie::bessel::calcPsi(n, z, Psi, dPsi); + nmie::bessel::calcZeta(n, z, Zeta, dZeta); + + // n_print = n_small; + // printf("csj[%i] = %+10.5er%+10.5ei; csy[%i] = %+10.5er%+10.5ei\n", n_print, + // csj[n_print].real(), csj[n_print].imag(), n_print, + // csy[n_print].real(), csy[n_print].imag()); + // printf("WA csj = e r e i; WA csy = e r e i\n"); + // n_print = n_big; + // printf("csj[%i] = %+10.5er%+10.5ei; csy[%i] = %+10.5er%+10.5ei\n", n_big, + // csj[n_print].real(), csj[n_print].imag(), n_print, + // csy[n_print].real(), csy[n_print].imag()); + // printf(" WA csj = e r e i; WA csy = e r e i\n"); + + // n_print = n_small; + // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print, + // csh[n_print].real(), csh[n_print].imag()); + // printf("WA csj = e r e i;\n"); + // n_print = n_big; + // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print, + // csh[n_print].real(), csh[n_print].imag()); + // printf("WA csj = e r e i;\n"); + + n_print = n_small; + printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_print, + Psi[n_print].real(), Psi[n_print].imag(), n_print, + Zeta[n_print].real(), Zeta[n_print].imag()); + printf("WA Psi = e r e i; WA Zeta = e r e i\n"); + n_print = n_big; + printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_big, + Psi[n_print].real(), Psi[n_print].imag(), n_print, + Zeta[n_print].real(), Zeta[n_print].imag()); + printf(" WA Psi = e r e i; WA Zeta = e r e i\n"); + + n_print = n_small; + printf("dPsi[%i] = %+10.5er%+10.5ei; dZeta[%i] = %+10.5er%+10.5ei\n", n_print, + dPsi[n_print].real(), dPsi[n_print].imag(), n_print, + dZeta[n_print].real(), dZeta[n_print].imag()); + printf("WA dPsi = e r e i; WA dZeta = -3.11025e-08r-2.90558e-08i\n"); + n_print = n_big; + printf("dPsi[%i] = %+10.5er%+10.5ei; dZeta[%i] = %+10.5er%+10.5ei\n", n_big, + dPsi[n_print].real(), dPsi[n_print].imag(), n_print, + dZeta[n_print].real(), dZeta[n_print].imag()); + printf(" WA dPsi = -2.47740e+05r+3.05899e+05i; WA dZeta = -6.13497e-07r-1.14967e-06i\n"); + + nmax_ = 100; + printf("----- Scattnlay (nmax_=%d)-----\n",nmax_); + calcPsiZeta(z, Psi, Zeta); + + n_print = n_small; + printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_print, + Psi[n_print].real(), Psi[n_print].imag(), n_print, + Zeta[n_print].real(), Zeta[n_print].imag()); + printf("WA Psi = e r e i; WA Zeta = e r e i\n"); + n_print = n_big; + printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_big, + Psi[n_print].real(), Psi[n_print].imag(), n_print, + Zeta[n_print].real(), Zeta[n_print].imag()); + printf(" WA Psi = e r e i; WA Zeta = e r e i\n"); + + + printf("===== Big real========\n"); + z = z_big_real; + printf("z = %+10.5er%+10.5ei\n", z.real(), z.imag()); + nmie::bessel::csphjy (n, z, nm, csj, cdj, csy, cdy ); + csh.resize(csj.size()); + for (unsigned int i = 0; i < csj.size(); ++i) + csh[i] = csj[i]+c_i*csy[i]; + nmie::bessel::calcPsi(n, z, Psi, dPsi); + nmie::bessel::calcZeta(n, z, Zeta, dZeta); + + // n_print = n_small; + // printf("csj[%i] = %+10.5er%+10.5ei; csy[%i] = %+10.5er%+10.5ei\n", n_print, + // csj[n_print].real(), csj[n_print].imag(), n_print, + // csy[n_print].real(), csy[n_print].imag()); + // printf("WA csj = e r e i; WA csy = e r e i\n"); + // n_print = n_big; + // printf("csj[%i] = %+10.5er%+10.5ei; csy[%i] = %+10.5er%+10.5ei\n", n_big, + // csj[n_print].real(), csj[n_print].imag(), n_print, + // csy[n_print].real(), csy[n_print].imag()); + // printf(" WA csj = e r e i; WA csy = e r e i\n"); + + // n_print = n_small; + // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print, + // csh[n_print].real(), csh[n_print].imag()); + // printf("WA csj = e r e i;\n"); + // n_print = n_big; + // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print, + // csh[n_print].real(), csh[n_print].imag()); + // printf("WA csj = e r e i;\n"); + + n_print = n_small; + printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_print, + Psi[n_print].real(), Psi[n_print].imag(), n_print, + Zeta[n_print].real(), Zeta[n_print].imag()); + // + //SphericalHankelH1(1,81)*81 + printf("WA Psi =-0.784462377 e e i; WA Zeta = -0.784462377r+0.620299278i\n"); + n_print = n_big; + printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_big, + Psi[n_print].real(), Psi[n_print].imag(), n_print, + Zeta[n_print].real(), Zeta[n_print].imag()); + printf(" WA Psi = 0.699946236744 e i; WA Zeta = 0.69994623 r+0.727239996i\n"); + + // n_print = n_small; + // printf("dPsi[%i] = %+10.5er%+10.5ei; dZeta[%i] = %+10.5er%+10.5ei\n", n_print, + // dPsi[n_print].real(), dPsi[n_print].imag(), n_print, + // dZeta[n_print].real(), dZeta[n_print].imag()); + // printf("WA dPsi = e r e i; WA dZeta = -6.20203e-01r-7.84344e-01i\n"); + // n_print = n_big; + // printf("dPsi[%i] = %+10.5er%+10.5ei; dZeta[%i] = %+10.5er%+10.5ei\n", n_big, + // dPsi[n_print].real(), dPsi[n_print].imag(), n_print, + // dZeta[n_print].real(), dZeta[n_print].imag()); + // printf(" WA dPsi = e r e i; WA dZeta = -7.13982e-01r+6.86858e-01i\n"); + + nmax_ = 100 ; + printf("----- Scattnlay (nmax_=%d)-----\n",nmax_); + calcPsiZeta(z, Psi, Zeta); + + n_print = n_small; + printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_print, + Psi[n_print].real(), Psi[n_print].imag(), n_print, + Zeta[n_print].real(), Zeta[n_print].imag()); + printf("WA Psi = e r e i; WA Zeta = e r e i\n"); + n_print = n_big; + printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_big, + Psi[n_print].real(), Psi[n_print].imag(), n_print, + Zeta[n_print].real(), Zeta[n_print].imag()); + printf(" WA Psi = e r e i; WA Zeta = e r e i\n"); + + + + printf("===== Big imag========\n"); + z = z_big_imag; + printf("z = %+10.5er%+10.5ei\n", z.real(), z.imag()); + nmie::bessel::csphjy (n, z, nm, csj, cdj, csy, cdy ); + csh.resize(csj.size()); + for (unsigned int i = 0; i < csj.size(); ++i) + csh[i] = csj[i]+c_i*csy[i]; + nmie::bessel::calcPsi(n, z, Psi, dPsi); + nmie::bessel::calcZeta(n, z, Zeta, dZeta); + + // n_print = n_small; + // printf("csj[%i] = %+10.5er%+10.5ei; csy[%i] = %+10.5er%+10.5ei\n", n_print, + // csj[n_print].real(), csj[n_print].imag(), n_print, + // csy[n_print].real(), csy[n_print].imag()); + // printf("WA csj = e r e i; WA csy = e r e i\n"); + // n_print = n_big; + // printf("csj[%i] = %+10.5er%+10.5ei; csy[%i] = %+10.5er%+10.5ei\n", n_big, + // csj[n_print].real(), csj[n_print].imag(), n_print, + // csy[n_print].real(), csy[n_print].imag()); + // printf(" WA csj = e r e i; WA csy = e r e i\n"); + + // n_print = n_small; + // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print, + // csh[n_print].real(), csh[n_print].imag()); + // printf("WA csj = e r e i;\n"); + // n_print = n_big; + // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print, + // csh[n_print].real(), csh[n_print].imag()); + // printf("WA csj = e r e i;\n"); + + // n_print = n_small; + // printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_print, + // Psi[n_print].real(), Psi[n_print].imag(), n_print, + // Zeta[n_print].real(), Zeta[n_print].imag()); + // printf("WA Psi = e r e i; WA Zeta = e r e i\n"); + // n_print = n_big; + // printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_big, + // Psi[n_print].real(), Psi[n_print].imag(), n_print, + // Zeta[n_print].real(), Zeta[n_print].imag()); + // printf(" WA Psi = e r e i; WA Zeta = e r e i\n"); + + n_print = n_small; + printf("dPsi[%i] = %+10.5er%+10.5ei; dZeta[%i] = %+10.5er%+10.5ei\n", n_print, + dPsi[n_print].real(), dPsi[n_print].imag(), n_print, + dZeta[n_print].real(), dZeta[n_print].imag()); + printf("WA dPsi = e r e i; WA dZeta = -2.47233e-22r-2.44758e-06i\n"); + n_print = n_big; + printf("dPsi[%i] = %+10.5er%+10.5ei; dZeta[%i] = %+10.5er%+10.5ei\n", n_big, + dPsi[n_print].real(), dPsi[n_print].imag(), n_print, + dZeta[n_print].real(), dZeta[n_print].imag()); + printf(" WA dPsi = e r e i; WA dZeta = 8.39449e-19r+1.28767e-02i\n"); + +} +