(Update Jan. 06, 2018)
Following code mistakenly read [result.txt] not [bisphere.print]. Use v0.4 instead of the code written below

bisphere.f > tool > calc_qext_bisphere_180104.py > v0.4 > calculate Qext from the result (R1, R2, CEXT) [bug fixed]

Environment

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)

result.txt(excerpt)

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:

calc_qext_bisphere_180104.py

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