Skip to content

Commit

Permalink
updated function names and return dictionaries for better usability, …
Browse files Browse the repository at this point in the history 
…added dsharp composition
  • Loading branch information
@mlietzow
mlietzow committed Sep 5, 2024
1 parent 8142147 commit 6e92a4e
Show file tree
Hide file tree
Showing 17 changed files with 936 additions and 1,117 deletions.
106 changes: 32 additions & 74 deletions README.md
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
# MIEX for Python
# MIEX-Python

is a Mie scattering code for large grains written in Python and based on [MIEX](https://ui.adsabs.harvard.edu/abs/2018ascl.soft10019W) by [Wolf & Voshchinnikov (2004)](https://ui.adsabs.harvard.edu/abs/2004CoPhC.162..113W).

Expand Down Expand Up @@ -29,7 +29,7 @@ For the original source code of MIEX written in `FORTRAN90`, we refer to

## Requirements

To run MIEX for Python, at least **Python 3.6** and the following packages are required:
To run MIEX-Python, at least **Python 3.6** and the following packages are required:
- numba
- numpy

Expand All @@ -42,55 +42,15 @@ pip3 install -r requirements.txt

## Installation

The package can be installed via `pip`:
The package can be installed via `pip`, change into the directory `MIEX-Python` and run:
```
pip install MIEX-Python/
pip install .
```

## Run MIEX

[run_miex.py](run_miex.py) executes MIEX with the parameters of the input file:
## Use MIEX-Python

```bash
python3 run_miex.py example1.input
python3 run_miex.py example2.input
```

The results are stored in the `results` directory.

The input file is organized as follows:

| Parameter | Type |
| ------------------------------------------------------------------------------ | ----- |
| Real refractive index of the surrounding medium | float |
| Number of wavelengths | int |
| Number of components | int |
| Name of the dust data files (lambda/n/k data) | |
|   1. component | str |
|   ... | str |
|   n. component | str |
| Relative abundances of the different components [%] <sup>1</sup> | |
| &ensp; 1. component | float |
| &ensp; ... | float |
| &ensp; n. component | float |
| Single grain size (1) or grain size distribution (2) | int |
| Grain size [micron] / Minimum grain size [micron] | float |
| Maximum grain size [micron] <sup>2</sup> | float |
| Size distribution exponent <sup>2</sup> | float |
| Number of size bins <sup>2</sup> | int |
| Calculate scattering matrix elements (0: no / 1: yes) | int |
| Number of scattering angles in the interval [0, 180]; odd number! <sup>3</sup> | int |
| Project name for the output files | str |
| Save results in separate files (0: no / 1: yes) | int |

<sup>1</sup> only for multiple compositions, omit for single composition;

<sup>2</sup> only for a grain size distribution, omit for single grain size;

<sup>3</sup> only if scattering matrix elements are calculated, omit if not.


### Run MIEX via streamlit
### Run MIEX-Python via streamlit

[miex_app.py](miex_app.py) can be used as a simple Streamlit web app: https://miex-python.streamlit.app/.

Expand All @@ -107,14 +67,13 @@ streamlit run miex_app.py
```


## Import MIEX

[miex.py](src/miex.py) can be imported and used in any python script (see [jupyter notebook 1](miex_notebook_1.ipynb) or [jupyter notebook 2](miex_notebook_2.ipynb)).
### Import MIEX-Python
After installation, MIEX-Python can be imported via `import miex` and used in any python script (see [jupyter notebook 1](miex_notebook_1.ipynb)).

To calculate the efficiency factors and scattering amplitude functions (optionally), use e.g.,

```python
result = miex.shexqnn2(x=1.0, m=complex(1.5, 0.0), nang=91, doSA=True)
result = miex.get_mie_coefficients(x=1.0, m=complex(1.5, 0.0), nang=91, doSA=True)
```

| Variable | Input Parameter | Type |
Expand All @@ -127,31 +86,34 @@ result = miex.shexqnn2(x=1.0, m=complex(1.5, 0.0), nang=91, doSA=True)
| `eps=1.0e-20` | Accuracy to be achieved | float, optional |
| `xmin=1.0e-6` | Minimum size parameter | float, optional |

| Variable | Output Parameter | Type | Index |
| ------------- | ------------------------------| ------------------- | ----- |
| $Q_{\rm ext}$ | Extinction efficiency | float | 0 |
| $Q_{\rm abs}$ | Absorption efficiency | float | 1 |
| $Q_{\rm sca}$ | Scattering efficiency | float | 2 |
| $Q_{\rm bk}$ | Backscattering efficiency | float | 3 |
| $Q_{\rm pr}$ | Radiation pressure efficiency | float | 4 |
| $A$ | Single scattering albedo | float | 5 |
| $g_{\rm sca}$ | Scattering asymmetry factor | float | 6 |
| $S_{1}$ | Scattering amplitude function | complex, array-like | 7 |
| $S_{2}$ | Scattering amplitude function | complex, array-like | 8 |
The function returns a dictionary with the following entries:

| Entry | Output Parameter | Type |
| ------------- | ------------------------------| ------------------- |
| `Q_ext` | Extinction efficiency | float |
| `Q_abs` | Absorption efficiency | float |
| `Q_sca` | Scattering efficiency | float |
| `Q_bk` | Backscattering efficiency | float |
| `Q_pr` | Radiation pressure efficiency | float |
| `Albedo` | Single scattering albedo | float |
| `g_sca` | Scattering asymmetry factor | float |
| `SA_1` | Scattering amplitude function | ndarray, complex |
| `SA_2` | Scattering amplitude function | ndarray, complex |
| `theta` | Scattering angles [rad] | ndarray |

The scattering amplitude functions are an array with size `2*nang-1`.

To calculate the scattering matrix elements, use

```python
sca_mat = miex.scattering_matrix_elements(S1, S2)
scat_mat = miex.get_scattering_matrix_elements(S1, S2)
```

where $S_1$ and $S_2$​ are the scattering amplitude functions.
The function returns the scattering matrix elements $S_{11}$​, $S_{12}$​, $S_{33}$​, $S_{34}$​ (in this order).
The function returns a dictionary with the scattering matrix elements $S_{11}$​, $S_{12}$​, $S_{33}$​, $S_{34}$ as entries for instance `scat_mat["S_11"]`.


## Test MIEX
### Test MIEX-Python

[test_miex.py](test_miex.py) includes some test routines. The results are compared with results by [Bohren & Huffman (1998)](https://doi.org/10.1002/9783527618156) and by [Wiscombe (1979)](https://doi.org/10.5065/D6ZP4414):

Expand All @@ -163,17 +125,13 @@ python3 test_miex.py
## Project structure

.
├── ri-data # Input data used by MIEX
├── ri-data # Input data used by MIEX-Python
├── README.me
├── example1.input # Exemplary input file
├── example2.input # Exemplary input file
├── src
│ └── miex.py # Source code of MIEX
├── miex_app.py # Python script to run MIEX via Streamlit
├── miex_notebook_1.ipynb # Jupyter notebook on how to use MIEX (basic)
├── miex_notebook_2.ipynb # Jupyter notebook on how to use MIEX (advanced)
├── requirements.txt # Required python packages for MIEX
├── run_miex.py # Python script to run MIEX with input files
├── miex
│ └── miex.py # Source code of MIEX-Python
├── miex_app.py # Python script to run MIEX-Python via Streamlit
├── miex_notebook_1.ipynb # Jupyter notebook on how to use MIEX-Python
├── requirements.txt # Required python packages for MIEX-Python
└── test_miex.py # Python script for test purposes


Expand Down
18 changes: 0 additions & 18 deletions example1.input
View file

This file was deleted.

14 changes: 0 additions & 14 deletions example2.input
View file

This file was deleted.

1 change: 1 addition & 0 deletions miex/__init__.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
from .miex import get_mie_coefficients, get_scattering_matrix_elements
128 changes: 109 additions & 19 deletions src/miex/miex.py → miex/miex.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,32 +103,32 @@ def shexqnn2(x, ri, nang=2, doSA=False, nterm=2e7, eps=1.0e-20, xmin=1.0e-06):
Returns
-------
0: Q_ext : float
Q_ext : float
extinction efficiency
1: Q_abs : float
Q_abs : float
absorption efficiency
2: Q_sca : float
Q_sca : float
scattering efficiency
3: Q_bk : float
Q_bk : float
backscattering efficiency
4: Q_pr : float
Q_pr : float
radiation pressure efficiency
5: albedo : float
albedo : float
single scattering albedo
6: g_sca : float
g_sca : float
scattering asymmetry factor
7: SA_1 : ndarray
scattering amplitude function. The length of the array is 2*nang-1
8: SA_2 : ndarray
SA_1, SA_2 : ndarrays
scattering amplitude function. The length of the array is 2*nang-1
theta : ndarray
scattering angles [rad]
Raises
------
Expand Down Expand Up @@ -218,10 +218,10 @@ def shexqnn2(x, ri, nang=2, doSA=False, nterm=2e7, eps=1.0e-20, xmin=1.0e-06):

# scattering amplitude functions
nang2 = 2 * nang - 1
mu = np.cos(np.linspace(0, np.pi, nang2))
SA_1 = np.zeros(nang2, dtype=np.complex128)
SA_2 = np.zeros(nang2, dtype=np.complex128)
if doSA:
mu = np.cos(np.linspace(0, np.pi, nang2))
fpi0 = np.zeros(nang2)
fpi1 = np.ones(nang2)

Expand Down Expand Up @@ -329,20 +329,103 @@ def shexqnn2(x, ri, nang=2, doSA=False, nterm=2e7, eps=1.0e-20, xmin=1.0e-06):
return Q_ext, Q_abs, Q_sca, Q_bk, Q_pr, albedo, g_sca, SA_1, SA_2


def scattering_matrix_elements(SA_1, SA_2):
"""Calculations of the scattering matrix elements
def get_mie_coefficients(x, ri, nang=2, doSA=False, nterm=2e7, eps=1.0e-20, xmin=1.0e-6):
"""Calculations of the efficiencies and scattering amplitude functions
Parameters
----------
SA_1 : array_like
scattering amplitude function
x : float
size parameter = 2 * pi * radius / wavelength
ri : complex float
complex refractive index
nang : int, optional, default = 2
half number of scattering angles theta in the intervall 0...pi/2 (equidistantly distributed)
doSA : bool, optional, default = False
calculation of the scattering amplitudes
nterm : int, optional, default = 2e7
Maximum number of terms to be considered
eps : float, optional, default = 1.0e-20
Accuracy to be achieved
xmin : float, optional, default = 1.0e-06
Minimum size parameter
Returns
-------
Dictionary with the following entries:
Q_ext : float
extinction efficiency
SA_2 : array_like
Q_abs : float
absorption efficiency
Q_sca : float
scattering efficiency
Q_bk : float
backscattering efficiency
Q_pr : float
radiation pressure efficiency
albedo : float
single scattering albedo
g_sca : float
scattering asymmetry factor
SA_1, SA_2 : ndarrays
scattering amplitude function. The length of the array is 2*nang-1
theta : ndarray
scattering angles [rad]
"""

Q_ext, Q_abs, Q_sca, Q_bk, Q_pr, albedo, g_sca, SA_1, SA_2 = shexqnn2(
x,
ri,
nang=nang,
doSA=doSA,
nterm=nterm,
eps=eps,
xmin=xmin
)

results = {
"Q_ext": Q_ext,
"Q_abs": Q_abs,
"Q_sca": Q_sca,
"Q_bk": Q_bk,
"Q_pr": Q_pr,
"Albedo": albedo,
"g_sca": g_sca,
"SA_1": SA_1,
"SA_2": SA_2,
"theta": np.linspace(0, np.pi, 2 * nang - 1),
}

return results


def get_scattering_matrix_elements(SA_1, SA_2):
"""Calculations of the scattering matrix elements
Parameters
----------
SA_1, SA_2 : ndarrays
scattering amplitude function
Returns
-------
S_11, S_12, S_33, S_34 : array_like
Dictionary with the following entries:
S_11, S_12, S_33, S_34 : ndarrays
scattering matrix elements
"""

Expand All @@ -361,4 +444,11 @@ def scattering_matrix_elements(SA_1, SA_2):
S_12[-1] = 0.0
S_34[-1] = 0.0

return S_11, S_12, S_33, S_34
results = {
"S_11": S_11,
"S_12": S_12,
"S_33": S_33,
"S_34": S_34,
}

return results
Loading

0 comments on commit 6e92a4e

Please sign in to comment.