動作環境

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)

関連: Jupyter | Matplotlib > 凹凸の激しい形状の可視化 > ADDA:IntField-Yファイルからの読込み v0.1 > A sphere | Chebyshev particle の表示

Jupyter | Matplotlib > 凹凸の激しい形状の可視化 > ADDA:IntField-Yファイルからの読込み v0.1 > A sphere | Chebyshev particle の表示
の続き。

This article is related to ADDA (light scattering simulator based on the discrete dipole approximation).

ADDAの出力ファイル（dipoleのx,y,z座標)を表示してみた。

失敗とその理由

v0.1において、下半球の値をとろうとした
- => 平坦な値になった

理由について気づいた。

ADDAの出力ファイルは球形の内部の座標も保持している。

// フィルタ条件

mz[iy][ix] = zar[idx]

上記の処理を繰り返した場合、iy,ixが同じ値を持つidxは複数あり、最後の値だけが残るようになる。

おそらく、ADDAの出力では、Z軸方向の値が大きいものが最後の値となるようになっており、それがたまたまv0.1のプログラムで上半球の座標取得に有効に働いた。

一方で、v0.1のプログラムをフィルタ条件で「zar[idx] > 0.0を排除」した場合、0.0付近の値がmzに入った。0.0付近の値がADDA出力の最後に来るようなファイル内配置のようだ。

v0.2 変更点

get_meshgrid_from_xyzArray()に引数 upperを追加
- upper = True: 上半球の処理
- upper = False: 下半球の処理
以下のフィルタを追加

if upper and zar[idx] < mz[iy][ix]:  # upper hemisphere
     continue
if not upper and zar[idx] > mz[iy][ix]:  # lower hemisphere
     continue

mzの初期値は以下とした。

DEFAULT_VAL = 1e-316
mz = np.empty([len(yuniq), len(xuniq)])
mz.fill(DEFAULT_VAL)

以上により、上半球、下半球、それぞれの「外型の」座標を取得できる。

入力ファイル

ADDAの出力ファイル（dipoleのx,y,z座標) を含むIntField-Y。

$ head IntField-Y 
x y z |E|^2 Ex.r Ex.i Ey.r Ey.i Ez.r Ez.i
-0.2166615618 -1.516630933 -5.416539045 0.5516786267 -0.02931061187 0.02652841941 0.1694307002 0.7044162375 -0.09245370677 -0.1290700275
0.2166615618 -1.516630933 -5.416539045 0.5516786268 0.0293106112 -0.02652841834 0.1694306999 0.7044162377 -0.09245370673 -0.1290700275
-1.083307809 -1.083307809 -5.416539045 0.5009708155 -0.1210777782 0.1210667177 0.3201499555 0.5959123603 -0.06902233265 -0.09634427508
-0.6499846854 -1.083307809 -5.416539045 0.5264383492 -0.11184864 0.05898299433 0.1712903828 0.6818062435 -0.06862119958 -0.1074254128
-0.2166615618 -1.083307809 -5.416539045 0.7409394521 -0.04313270512 0.006539303629 0.1179657327 0.8426239137 -0.06420657934 -0.1047988574
0.2166615618 -1.083307809 -5.416539045 0.7409394524 0.04313270474 -0.00653930278 0.1179657325 0.8426239139 -0.06420657931 -0.1047988574
0.6499846854 -1.083307809 -5.416539045 0.5264383496 0.1118486397 -0.0589829935 0.1712903821 0.6818062441 -0.0686211995 -0.1074254129
1.083307809 -1.083307809 -5.416539045 0.5009708155 0.1210777777 -0.121066717 0.3201499543 0.5959123613 -0.06902233252 -0.09634427512
-1.083307809 -0.6499846854 -5.416539045 0.6583579786 -0.0961687486 0.04927804902 0.3522909973 0.7189310502 -0.04061949379 -0.0637218901

code v0.2

showChebyshev_170909.ipynb

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from pylab import rcParams
import sys


'''
v0.2 Sep. 09, 2017
   - fix bug > values of lower hemisphere is filled with the values near the equatorial planes
      + taking "the last" values of the [iy][ix]
v0.1 Sep. 09, 2017
   - read [IntField-Y] file

=== based on [toMeshgrid_170902.ipynb] ===

v0.1 Sep. 02, 2017
   - add get_meshgrid_from_xyzArray()
   - add func_z()
   - display 3D surface with lines
      + ref: https://stackoverflow.com/questions/30497737/applying-colormaps-to-custom-axis-in-matplotlib-3d-surface
'''

# coding rule: PEP8


rcParams['figure.figsize'] = 15, 10


def get_meshgrid_from_xyzArray(xar, yar, zar, upper):
    # mx, my, mz : in meshgrid
    #
    DEFAULT_VAL = 1e-316

    xuniq = np.unique(xar)
    yuniq = np.unique(yar)
    mz = np.empty([len(yuniq), len(xuniq)])
    mz.fill(DEFAULT_VAL)

    for ix in range(len(xuniq)):
        for iy in range(len(yuniq)):
            xx, yy = xuniq[ix], yuniq[iy]
            for idx in range(len(xar)):
                tx, ty = xar[idx], yar[idx]
                if abs(tx - xx) >= sys.float_info.epsilon:
                    continue
                if abs(ty - yy) >= sys.float_info.epsilon:
                    continue
                if upper and zar[idx] < mz[iy][ix]:  # upper hemisphere
                    continue
                if not upper and zar[idx] > mz[iy][ix]:  # lower hemisphere
                    continue

                mz[iy][ix] = zar[idx]
    mx, my = np.meshgrid(xuniq, yuniq)
    return mx, my, mz

INPFILE = 'IntField-Y'

dat = np.genfromtxt(INPFILE, delimiter=' ', skip_header=1)

xpar, ypar, zpar = [], [], []
for elem in dat:
    axp, ayp, azp = elem[:3]
    xpar += [axp]
    ypar += [ayp]
    zpar += [azp]

# 2. from x,y,z arrays
res1 = get_meshgrid_from_xyzArray(xpar, ypar, zpar, upper=True)
gx1, gy1, gz1 = res1
res2 = get_meshgrid_from_xyzArray(xpar, ypar, zpar, upper=False)
gx2, gy2, gz2 = res2

ax1 = plt.subplot2grid((1, 2), (0, 0), projection='3d')
surf1 = ax1.plot_surface(gx1, gy1, gz1, shade=False,
                         facecolors=plt.cm.Set2((gx1-gx1.min())/(gx1.max()-gx1.min()))
                         )

ax2 = plt.subplot2grid((1, 2), (0, 1), projection='3d')
surf2 = ax2.plot_surface(gx2, gy2, gz2, shade=False,
                         facecolors=plt.cm.Set2((gx2-gx2.min())/(gx2.max()-gx2.min()))
                         )

plt.draw()

#   draw lines on the surface
lines1 = np.array(surf1.get_edgecolor())
lines2 = np.array(surf2.get_edgecolor())
# make lines white, and keep alpha==1. It's an array of colors like this: [r,g,b,alpha]
surf1.set_edgecolor(lines2 * np.array([0, 0, 0, 0]) + 1)
surf2.set_edgecolor(lines2 * np.array([0, 0, 0, 0]) + 1)

plt.show()

結果 > spherical particles

2つのグラフに表示する形ではあるが、だいだいの形状は見て取れるようになった。

結果 > Chebyshev particles

結果 > Prism

所感

細かい凹凸のある粒子はMeshgridで描画するに適さないような気がする。

TODO

Cartesian coordinateをPolar coordinateに変換 > link
Polar coordinateにて外側の位置座標をMeshgridに拾う
plot_surfaceで表示

https://stackoverflow.com/questions/11140163/python-matplotlib-plotting-a-3d-cube-a-sphere-and-a-vector
https://stackoverflow.com/questions/31768031/plotting-points-on-the-surface-of-a-sphere-in-pythons-matplotlib