(Update Jan. 06, 2018)
Following code mistakenly read [result.txt] not [bisphere.print]. Use v0.4 instead of the code written below
GeForce GTX 1070 (8GB)
ASRock Z170M Pro4S [Intel Z170chipset]
Ubuntu 16.04 LTS desktop amd64
TensorFlow v1.2.1
cuDNN v5.1 for Linux
CUDA v8.0
Python 3.5.2
IPython 6.0.0 -- An enhanced Interactive Python.
gcc (Ubuntu 5.4.0-6ubuntu1~16.04.4) 5.4.0 20160609
GNU bash, version 4.3.48(1)-release (x86_64-pc-linux-gnu)
scipy v0.19.1
geopandas v0.3.0
MATLAB R2017b (Home Edition)
ADDA v.1.3b6
Related: Superposition T-matrix codes > bisphere.f > 2つの球の光散乱計算コード (analytical orientation averaging) in Japanese
In order to compare results from analytical orientation averaging (T-matrix method) and those from numerical orinetation averaging (ADDA + pySpherepts), both results need to have the same unit (i.e. Qext: Extinction efficiency).
Cext and Qext for bisphere
- Cext: Extinction cross section
- Qext: Extinction efficiency
The above two are related to:
Q_{ext} = \frac{C_{ext}}{\pi {r_{eff}}^2} ... (1)
where $r_{eff}$ is the equivalent volume sphere radius (i.e. the radius of the sphere having the same volume with the bisphere).
For a bisphere with composing sphere radii (R1 and R2), $r_{eff}$ can be obtained based on:
\frac{4}{3}\pi r_1^3 + \frac{4}{3}\pi r_2^3 = \frac{4}{3}\pi r_{eff}^3 ... (2-1)
r_{eff}^3 = r_1^3 + r_2^3 ... (2-2)
r_{eff} = \sqrt[3]{r_1^3 + r_2^3} ... (2-3)
code v0.1-v0.3
bisphere.f outputs files as the following (named as result.txt)
This test result was computed on an IBM RISC model 397 workstation
for the current code setting.
R1= .5000D+01 R2= .5000D+01 R12= .2000D+02 LAM= .6283D+01
X12= .2000D+02 NODRT= 33
X1= .500000D+01 N1= .150000D+01 K1= .500000D-02 NODR(1)= 14
X2= .500000D+01 N2= .150000D+01 K2= .500000D-02 NODR(2)= 14
CEXT= .589017D+03 CSCA= .567076D+03 W= .962749D+00 <COS>= .717082D+00
TEST OF VAN DER MEE & HOVENIER IS SATISFIED
S ALPHA1 ALPHA2 ALPHA3 ALPHA4 BETA1 BETA2
0 1.00000 .00000 .00000 .94160 .00000 .00000
1 2.15125 .00000 .00000 2.21375 .00000 .00000
2 2.99437 4.00361 3.83857 2.89376 -.13647 .11745
...
Code to calculate the Qext is as follows:
import numpy as np
import sys
'''
v0.3 Jan. 04, 2018
- fix bug: calc_qext() uses reff not reff^2
v0.2 Jan. 04, 2018
- rename to [calc_qext_bisphere_180104.py]
- was [calc_qext_180104.py]
- read_param_fromfile() uses convert_to_value()
- add convert_to_value()
- add calc_qext()
- add calc_reff()
v0.1 Jan. 04, 2018
- add read_param_fromfile()
'''
def convert_to_value(valstr):
# to avoid:
# TypeError: unsupported operand type(s) for
# ** or pow(): 'str' and 'int'
work = valstr.replace("D", "E")
return np.float64(work)
def read_param_fromfile(akey):
# 1. find line
keyline = []
with open('result.txt') as fd:
lines = fd.readlines()
for aline in lines:
if not akey + "=" in aline:
continue
keyline += [aline]
if len(keyline) == 0:
return None
# 2. extract value next to the keyword
found = False
for elem in keyline[0].split(' '):
if found:
return convert_to_value(elem)
if akey in elem:
found = True
return None # just in case
def calc_reff(r1, r2):
# volume equivalent sphere radius
return (r1**3 + r2**3)**(1./3.)
def calc_qext(r1, r2, cext):
reff = calc_reff(r1, r2)
denom = np.pi * reff * reff
if abs(denom < sys.float_info.epsilon):
return 0.0
return cext / denom
r1 = read_param_fromfile("R1")
r2 = read_param_fromfile("R2")
cext = read_param_fromfile("CEXT")
reff = calc_reff(r1, r2)
qext = calc_qext(r1, r2, cext)
print("r1, r2, reff, cext")
print("%.3e %.3e %.5e %.3e" % (r1, r2, reff, cext))
print("qext")
print("%.7e" % qext)
Example run
$ python3 calc_qext_bisphere_180104.py
r1, r2, reff, cext
5.000e+00 5.000e+00 6.29961e+00 5.890e+02
qext
4.7244503e+00