Planned maintenance on Wednesday, 2021-01-20, 17:00-18:00. Expect some interruptions during that time

Commit ee2af9f5 authored by Martin Glatzle's avatar Martin Glatzle

Merge branch 'mie' into dev

parents d026bc19 7f010ad4
......@@ -8,3 +8,6 @@
This diff is collapsed.
The contents of this file are based to a large extent on the source code of
Cloudy 17.01, which has been published with the below license.
Martin Glatzle
Cloudy -- simulations of non-equilibrium plasmas and their spectra
Copyright (C) 1978-2017 Gary J. Ferland and others.
This software is provided 'as-is', without any express or implied
warranty. In no event will the author be held liable for any damages
arising from the use of this software.
Permission is granted to anyone to use this software for any purpose,
including commercial applications, and to alter it and redistribute it
freely, subject to the following restrictions:
1. The origin of this software must not be misrepresented; you must not
claim that you wrote the original software. If you use this software
in a product, an acknowledgment in the product documentation would be
appreciated but is not required.
2. Altered source versions must be plainly marked as such, and must not be
misrepresented as being the original software.
3. This notice may not be removed or altered from any source distribution.
Gary J. Ferland
This is an Open Source License as approved by the Open Source Initiative,
see for details.
#ifndef CMIE_HPP
#define CMIE_HPP
Conversion of complex dielectric function to complex refractive
index. Utility function.
#ifdef __cplusplus
extern "C" {
void dftori(/*@out@*/ double *nre,
/*@out@*/ double *nim,
double eps1,
double eps2);
#ifdef __cplusplus
Conversion of complex refractive index to complex dielectric
function. Utility function.
#ifdef __cplusplus
extern "C" {
void ritodf(double nre,
double nim,
/*@out@*/ double *eps1,
/*@out@*/ double *eps2);
#ifdef __cplusplus
Error values for mie_cs.
#define ANOSUCCESS 1
#define MIESUCCESS 0
#define TOOSMALL -1
#define TOOLARGE -2
#define FAILURE -4
Calculate absorption and scattering cross sections as well as scattering
asymmetry parameter g (aka cosb) for a given wavelength and a single grain
defined by its size and complex refractive index .
At first, a Mie computation is attempted, if that fails, anomalous
diffraction theory is employed (this results in cosb being invalid and a
corresponding error value being returned). If this also fails, a failure
error value is returned.
nr1 is Re(n)-1 with the complex refractive index n. It is only used in
anomalous diffraction theory computations and is a separate parameter (in
addtition to nre) to allow for increased numerical precision. For situations
that do not require the anomalous diffraction approximation, it can safely be
set to nre-1.
The minimum/maximum grain sizes accepted are 0.0001/10 micron. The lower
limit is physically motivated, since such a small grain is not really a solid
any more. The upper limit is motivated by convergence problems in the series
expansion used in Mie theory. For grains sizes outside this range,
corresponding error values are returned.
Checks on the computed results are performed if a calculation is successful
and if unphysical values are detected, a corresponding error value is
#ifdef __cplusplus
extern "C" {
void mie_cs(double wavlen, /* micron */
double grain_size, /* micron */
double nre,
double nim,
double nr1,
/*@out@*/ double *cs_abs, /* cm^2 */
/*@out@*/ double *cs_sct, /* cm^2 */
/*@out@*/ double *cosb,
/*@out@*/ int *error);
#ifdef __cplusplus
#include <Python.h>
#include <math.h>
#include <numpy/ndarraytypes.h>
#include <numpy/ufuncobject.h>
#include "cmie.hpp"
/* Module name is cmie, ufunc name is cs. */
/* function definitions (possibly for various dtypes) */
static void cs_double(char** args, npy_intp* dimensions,
npy_intp* steps, void* data)
npy_intp n = dimensions[0]; /* number of wavelengths*number of grain sizes */
char* wave = args[0];
char* size = args[1];
char* n_re = args[2];
char* n_im = args[3];
char* nr_1 = args[4];
char* cs_a = args[5];
char* cs_s = args[6];
char* cosb = args[7];
char* exit = args[8];
npy_intp wave_step = steps[0];
npy_intp size_step = steps[1];
npy_intp n_re_step = steps[2];
npy_intp n_im_step = steps[3];
npy_intp nr_1_step = steps[4];
npy_intp cs_a_step = steps[5];
npy_intp cs_s_step = steps[6];
npy_intp cosb_step = steps[7];
npy_intp exit_step = steps[8];
npy_intp i;
for (i = 0; i < n; i++) {
mie_cs(*(double*)wave, *(double*)size, *(double*)n_re, *(double*)n_im,
(double*)cs_a, (double*)cs_s, (double*)cosb, (int*)exit);
wave += wave_step;
size += size_step;
n_re += n_re_step;
n_im += n_im_step;
nr_1 += nr_1_step;
cs_a += cs_a_step;
cs_s += cs_s_step;
cosb += cosb_step;
exit += exit_step;
/* array of pointers to all cs functions */
PyUFuncGenericFunction cs_funcs[1] = {&cs_double};
/* input and return dtypes for all cs function */
static char cs_types[9] = {NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE,
static void* cs_data[1] = {NULL};
module initialization
static PyMethodDef CmieMethods[] = {
static struct PyModuleDef moduledef = {
PyMODINIT_FUNC PyInit_cmie(void)
PyObject* m, * cs, * d;
m = PyModule_Create(&moduledef);
if (!m) {
return NULL;
/* create the ufunc */
cs = PyUFunc_FromFuncAndData(cs_funcs,
1, /* number of cs functions */
5, /* number of input args for each */
4, /* number of output args for each */
d = PyModule_GetDict(m);
PyDict_SetItemString(d, "cs", cs);
return m;
end module initialization
# small makefile for quick test compilations
# NOTE: in order to run the produced executable, need to set LD_LIBRARY_PATH
# since libcmie is dynamic.
CXXFLAGS=-Wall -Wextra -Wconversion -fmax-errors=1 -pedantic -cpp
all: $(EXE)
$(EXE): makefile $(MAIN) $(LIB)
@echo "# building exe"
$(CXX) $(CXXFLAGS) -L./ -lcmie -o $(EXE) test.cpp
$(LIB): lib
@echo "# linking lib"
$(CXX) $(CXXFLAGS) -shared -Wl,-soname,$(LIB) -o $(LIB) $(OBJECTS)
$(OBJECTS): makefile
%.o: %.cpp
@echo "# compiling $*.cpp"
$(CXX) -fPIC $(CXXFLAGS) -o $@ -c $<
.PHONY: clean
$(RM) *.o *~ $(EXE) $(LIB) 2> /dev/null; true
#include "cmie.hpp"
#include <iostream>
int main(){
double wave, a, nre, nim, nr1;
double ca, cs, g;
int err;
wave = 0.1;
a = 0.1;
nre = 1;
nim = 0.3;
nr1 = 0;
wave, a, nre, nim, nr1,
&ca, &cs, &g, &err
std::cout << ca << " " << cs << " " << " " << g << " " << err << "\n";
return 0;
......@@ -6,10 +6,12 @@ with them easier.
useful things and the API _will_ change.
## Installation
- `numpy` and `setuptools` are required for installation, get them from `pip`:
`pip install numpy setuptools`
- Clone the repo: `git clone`
- `cd cosmic_dustbox`
- `python install`
- `pip install .`
## Examples
Note: `cosmic_dustbox` uses `astropy.units` throughout.
......@@ -19,8 +19,8 @@ class Crefin(object):
micron. It shall return a 2-D array of the complex refractive index
evaluated at the given sizes and wavelengths. The first axis of this
array shall correspond to grain size and the second axis to
wavelength. Note that this axis ordering is reversed wrt to the output
of `np.meshgrid`.
wavelength. Note that this axis ordering is reversed wrt the output of
def __init__(self, f):
......@@ -119,9 +119,9 @@ class Crefin(object):
2-D interpolation of complex function with input-sorted output.
Returns interpolator for R**2 -> C function. When the interpolator is
called, it sorts the output array according to how the input arrays for
`x` and `y` are sorted.
Returns interpolator for :math:`\mathcal{R}**2` -> :math:`\mathcal{C}`
function. When the interpolator is called, it sorts the output array
according to how the input arrays for `x` and `y` are sorted.
......@@ -188,8 +188,7 @@ class SGPAHCrefin(Crefin):
for j, path in enumerate(paths):
size, data = cls.parseCrefinFile(path)
if size in a:
raise ValueError("Grain size "
+ str(size) + "provided twice.")
raise ValueError("Grain size " + str(size) + "provided twice.")
if j == 0:
......@@ -197,9 +196,9 @@ class SGPAHCrefin(Crefin):
if not _np.array_equal(lam, data[:, 0]):
raise ValueError(
"Wavelengths in file " + path
+ "differ from those in the other files visited "
"so far.")
"Wavelengths in file " + path + "differ from those "
"in the other files visited so far."
nt = data[:, 3] + 1 + 1j*data[:, 4]
if len(a) > 1:
from setuptools import setup
from setuptools import setup, Extension
import metadata as md
import numpy as np
cmiemod = Extension(
["Mie/cloudy/cmiemodule.c", "Mie/cloudy/cmie.cpp"],
......@@ -22,4 +29,5 @@ setup(
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment