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)


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]
v0.1 Jan. 04, 2018
'''

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)

# 1. find line
keyline = []
with open('result.txt') as fd:
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):
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

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

セブンオブナインです。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