Help us understand the problem. What is going on with this article?

[obsolete] bisphere.f > tool > calc_qext_bisphere_180104.py > v0.1-v0.3 > calculate Qext from the result (R1, R2, CEXT)

More than 3 years have passed since last update.

(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
7of9
セブンオブナインです。Unimatrix 01の第三付属物 9の7という識別番号です。Star trek Voyagerの好きなキャラクターです。まとめ記事は後日タイトルから内容がわからなくなるため、title検索で見つかるよう個々の記事にしてます。いわゆるBorg集合体の有名なセリフから「お前たち(の知識)を吸収する。抵抗は無意味だ」。Thanks in advance.
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away