Superposition T-matrix codes > bisphere.f > 2つの球の光散乱計算コード (analytical orientation averaging)

動作環境

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 Fortran (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

概要

sphere, bi-sphere, elllipsoidなどのいくつかの形状に対してはT-matrix法に基づく「analytical orientation averaging」の計算手法が確立している。

T-matrix + analytical orientation averagingによる手法は精度高く計算が可能である一方で、計算可能な形状に制限がある。計算が不可能な形状に対しては、T-matrix法とは別の手法を用いる場合がある(FDTD法やDDA法など)。それら手法に対してanalytical orientation averaging手法が利用できない場合にnumerical orientation averagingを使うことになる。

一方で、numerical orientation averagingの計算の精度を確認するために、上記のT-matrix + analytical orientation averagingの結果と比較する、ということを行うことがある。

ここでは、bi-sphereに対するT-matrix + analytical orientation averagingの計算コードの利用について確認をする。

bi-sphere code

code取得

bi-sphereの場合は以下のサイトからコードが利用可能である。
https://www.giss.nasa.gov/staff/mmishchenko/t_matrix.html

Superposition T-matrix codes
Double-precision superposition T-matrix code for randomly oriented two-sphere clusters with externally touching or separated components: bisphere.f

bisphere.fを適当なディレクトリに取得する。
bisphere.fにはヘッダー部分に計算結果例が掲載されている。このままでは実行できないため、下記の処理を行う。

$ head -n 89 bisphere.f > result.txt
$ cat bisphere.f | awk 'NR>90' > bisphere_work.f

bisphere_work.fを計算に使う。

実行

$ gfortran bisphere_work.f 
bisphere_work.f:1655:72: Warning: Deleted feature: PAUSE statement at (1)
bisphere_work.f:1670:72: Warning: Deleted feature: PAUSE statement at (1)

$ ./a.out 
FAC(165)=0.5424+296
 TIME =   34.85 MIN
Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL

bisphere.printというファイルが生成される。
result.txtの計算結果と比較を行い、同じ値であることを確認する。
(コンパイラの種類により、数値の表示が異なる場合はある(例: 0.94160と.94160)）

結果例

• 設定 (bisphere_work.fの160行あたりに設定されている)
  • N(1),K(1): 1.5 + 0.005i # 一つ目の球の複素屈折率
  • N(2),K(2): 1.5 + 0.005i # 二つ目の球の複素屈折率
  • R1: 5.0 # 一つ目の球の半径
  • R2: 5.0 # 二つ目の球の半径
  • R12: 20.0 # 二つの球の距離
  • LAM: 2pi # 波長 [um]

上記の設定に対する結果は下記となる。

bisphere.print

R1=0.5000D+01  R2=0.5000D+01  R12=0.2000D+02  LAM=0.6283D+01
X12=0.2000D+02  NODRT= 33
X1=0.500000D+01 N1=0.150000D+01 K1=0.500000D-02 NODR(1)= 14
X2=0.500000D+01 N2=0.150000D+01 K2=0.500000D-02 NODR(2)= 14
 CEXT=0.589017D+03  CSCA=0.567076D+03    W=0.962749D+00  <COS>=0.717082D+00
 TEST OF VAN DER MEE & HOVENIER IS SATISFIED
 
  S      ALPHA1      ALPHA2      ALPHA3      ALPHA4       BETA1       BETA2
  0     1.00000     0.00000     0.00000     0.94160     0.00000     0.00000
  1     2.15125     0.00000     0.00000     2.21374     0.00000     0.00000
  2     2.99437     4.00361     3.83857     2.89376    -0.13647     0.11745
  3     3.13877     3.55605     3.63214     3.20301    -0.11917     0.00216
  4     3.32190     3.77498     3.68288     3.26877    -0.15255     0.10368
  5     3.32936     3.46498     3.50235     3.40753    -0.19897    -0.04975
  6     3.18348     3.58119     3.45692     3.09234    -0.23832     0.00832
  7     2.95333     2.98997     3.11209     3.13450    -0.27378    -0.06001
  8     2.47715     2.91924     2.73441     2.37083    -0.43506    -0.33890
  9     1.94766     1.95517     1.99875     2.03174    -0.08178    -0.21840
 10     1.62692     1.85138     1.74819     1.56046    -0.15704    -0.30806
 11     1.20454     1.22355     1.22116     1.20713     0.04472    -0.06978
 12     1.27983     1.27994     1.27321     1.27673     0.01441    -0.00550
 13     1.39985     1.39649     1.39540     1.40073     0.01718     0.01608
 14     1.52960     1.52807     1.52819     1.53191     0.01704     0.02120
 15     1.66183     1.65843     1.65642     1.66209     0.02479     0.02294
 16     1.79873     1.79692     1.79303     1.79602     0.03754     0.01614
 17     1.88397     1.89037     1.88688     1.87990     0.04223     0.01045
 18     1.87931     1.89487     1.89122     1.87306     0.04009     0.00650
 19     1.77319     1.79585     1.79750     1.77103     0.02763     0.00827
 20     1.59409     1.62019     1.62268     1.59283     0.00972     0.00493
 21     1.37720     1.40242     1.40635     1.37792    -0.00433     0.00240
 22     1.15582     1.17892     1.18291     1.15756    -0.01639    -0.00189
 23     0.93894     0.95985     0.96444     0.94251    -0.02933    -0.00688
 24     0.73777     0.75635     0.75982     0.74154    -0.03810    -0.01711
 25     0.55240     0.56872     0.57050     0.55550    -0.03998    -0.02786
 26     0.38710     0.40096     0.40023     0.38810    -0.03397    -0.03530
 27     0.24605     0.25677     0.25475     0.24563    -0.02365    -0.03396
 28     0.14064     0.14780     0.14516     0.13914    -0.01241    -0.02673
 29     0.06992     0.07398     0.07220     0.06878    -0.00530    -0.01628
 30     0.03132     0.03330     0.03213     0.03048    -0.00145    -0.00876
 31     0.01211     0.01294     0.01241     0.01172    -0.00016    -0.00377
 32     0.00430     0.00460     0.00437     0.00411     0.00013    -0.00148
 33     0.00133     0.00143     0.00135     0.00127     0.00011    -0.00049
 34     0.00038     0.00041     0.00039     0.00036     0.00005    -0.00015
 35     0.00010     0.00011     0.00010     0.00009     0.00002    -0.00004
 36     0.00002     0.00003     0.00002     0.00002     0.00001    -0.00001
 37     0.00001     0.00001     0.00001     0.00000     0.00000    -0.00000
 
      <       F11       F22       F33       F44       F12       F34
   0.00  49.78455  49.77951  49.77951  49.77447   0.00000   0.00000
   5.00  35.98162  35.97666  35.97626  35.97159   0.10705   0.13306
  10.00  16.24586  16.24111  16.23788  16.23416   0.21533   0.23828
  15.00  10.67507  10.67056  10.65748  10.65496   0.39163   0.34560
  20.00   8.90278   8.89855   8.85575   8.85437   0.71344   0.48493
  25.00   4.67417   4.67021   4.59943   4.59892   0.73578   0.30617
  30.00   1.90103   1.89726   1.82180   1.82187   0.51177   0.01957
  35.00   1.11642   1.11257   1.03461   1.03518   0.33577  -0.19571
  40.00   1.10229   1.09805   1.03145   1.03264   0.15158  -0.32446
  45.00   1.13873   1.13393   1.09973   1.10167   0.00530  -0.26081
  50.00   1.00942   1.00424   0.99510   0.99766  -0.02457  -0.11589
  55.00   0.85571   0.85046   0.84877   0.85169   0.02002  -0.00392
  60.00   0.64929   0.64408   0.63356   0.63674   0.09647   0.04628
  65.00   0.43570   0.43039   0.40340   0.40693   0.14346   0.01536
  70.00   0.31993   0.31440   0.27708   0.28099   0.13643  -0.04728
  75.00   0.28821   0.28250   0.24356   0.24780   0.10457  -0.09298
  80.00   0.27305   0.26716   0.23409   0.23865   0.06754  -0.10614
  85.00   0.22690   0.22087   0.19891   0.20376   0.03623  -0.08551
  90.00   0.16608   0.15991   0.14797   0.15309   0.01953  -0.05383
  95.00   0.13383   0.12732   0.12010   0.12564   0.01720  -0.03445
 100.00   0.14155   0.13466   0.12760   0.13359   0.02517  -0.02894
 105.00   0.16635   0.15925   0.14929   0.15557   0.04271  -0.02521
 110.00   0.17408   0.16672   0.14829   0.15490   0.06852  -0.01610
 115.00   0.14870   0.14109   0.10365   0.11058   0.09019  -0.00940
 120.00   0.10957   0.10169   0.03337   0.04060   0.08931  -0.02033
 125.00   0.09606   0.08760  -0.01761  -0.00976   0.05673  -0.05878
 130.00   0.13177   0.12253  -0.01147  -0.00280  -0.00159  -0.11762
 135.00   0.20600   0.19550   0.05854   0.06853  -0.06144  -0.17075
 140.00   0.28155   0.26925   0.16422   0.17607  -0.08667  -0.18663
 145.00   0.31664   0.30306   0.25096   0.26414  -0.05212  -0.14696
 150.00   0.29863   0.28522   0.26592   0.27896   0.03527  -0.06289
 155.00   0.25450   0.24319   0.18699   0.19797   0.13353   0.03002
 160.00   0.23308   0.22382   0.03195   0.04092   0.19088   0.09261
 165.00   0.27032   0.26186  -0.15254  -0.14430   0.17895   0.10303
 170.00   0.35827   0.34853  -0.32017  -0.31056   0.10971   0.07009
 175.00   0.45586   0.43878  -0.43657  -0.41953   0.03254   0.02295
 180.00   0.50091   0.47807  -0.47807  -0.45522   0.00000   0.00000

備考

• 実行時にNote:floating-point exceptionsが表示される。
  • 計算結果は期待された値にはなっている。
• R12が20であり、二つの球が接触した計算ではない
  • R12=R(1)+R(2)とすることで接触した状態の計算ができる