==================
可視化（アニメ）
==================

.. toctree::
   :maxdepth: 1

惑星座標を使ってアニメーションを作ってみよう．
アニメーションであれば視覚的に動きを捉えることができる．

.. index:: ParaView

ParaViewの使用
=========================

ParaViewはオープンソース，マルチプラットフォームの可視化ソフトウェアである．
大規模なデータにも対応できるようである．

`>>ParaView <http://www.paraview.org/>`_ のDownloadからインストーラがダウンロードできる．
Windows,MacOS,Linuxで使える．

ParaView向けのデータ出力
===================================

ParaViewは様々な形式のデータを扱うことができるが，VTK XMLフォーマットを用いる．
中でもvtkUnstructuredGrid形式（拡張子は「.vtu」）を用いる．
以下のような形式である．
タグで囲んだ間に粒子の座標や速度などの情報を書き入れる．

::

    <VTKFile type=”UnstructuredGrid” ...>
    <UnstructuredGrid>
    <Piece NumberOfPoints=”#” NumberOfCells=”#”>
    <PointData>...</PointData>
    <CellData>...</CellData>
    <Points>...</Points>
    <Cells>...</Cells>
    </Piece>
    </UnstructuredGrid>
    </VTKFile>

.. note::
     〔 :title-reference:`粒子法入門`，越塚誠一ら（丸善出版，2014） 〕

ということで，プログラミングコードの中にvtkUnstructuredGrid形式を出力できるサブルーチンを加える．

.. code-block:: fortran
    :linenos:

    module output

    use globals

    implicit none
    
    contains
        subroutine write_vtk(step,c,v,m)
            integer,intent(in)::                step
            double precision,intent(in)::       c(ndim,nbody)
            double precision,intent(in)::       v(ndim,nbody)
            double precision,intent(in)::       m(nbody)

            integer::                   unit_vtu=11
            integer                     i,j
            double precision            v2
            double precision            ct(3)
            character(len=27)           file_vtu
            character(len=20)           step_str

            write(step_str,'(i20)')step
            file_vtu = step_str // '.vtu'
            file_vtu = 'pv_' // adjustl(file_vtu)
            file_vtu = trim(adjustl(file_vtu))
            open(unit=unit_vtu,file=file_vtu,status='replace')
            write(unit_vtu,'(a)')"<?xml version='1.0' encoding='UTF-8'?>"
            write(unit_vtu,'(a)')"<VTKFile xmlns='VTK' byte_order='LittleEndian' version='0.1' type='UnstructuredGrid'>"
            write(unit_vtu,'(a)')"<UnstructuredGrid>"
            write(unit_vtu,'(a,i0,a,i0,a)')"<Piece NumberOfCells='",nbody,"' NumberOfPoints='",nbody,"'>"
            write(unit_vtu,'(a)')" "

            write(unit_vtu,'(a)')"<Points>"
            write(unit_vtu,'(a,i0,a)')"<DataArray NumberOfComponents='",3,"' type='Float32' Name='Position' format='ascii'>"
            do i = 1,nbody
                ct(1:ndim) = c(1:ndim,i)
                do j = 1,3
                    if(ndim .eq. 2) then
                        ct(3) = 0.0d0
                    end if
                    write(unit_vtu,'(f8.3)',advance='no')ct(j)
                end do
                write(unit_vtu,*)
            end do
            write(unit_vtu,'(a)')"</DataArray>"

            write(unit_vtu,'(a)')"</Points>"
            write(unit_vtu,'(a)')"<PointData>"
            write(unit_vtu,'(a)')"<DataArray NumberOfComponents='1' type='Float32' Name='Velocity' format='ascii'>"
            do i = 1,nbody
                v2 = 0.0d0
                do j = 1,ndim
                    v2 = v2 + v(j,i)**2
                end do
                write(unit_vtu,'(f8.3)')sqrt(v2)
            end do
            write(unit_vtu,'(a)')"</DataArray>"

            write(unit_vtu,'(a)')"<DataArray NumberOfComponents='1' type='Float32' Name='Mass' format='ascii'>"
            do i = 1,nbody
                write(unit_vtu,'(f8.3)')log10(m(i)*10.0d8)
            end do
            write(unit_vtu,'(a)')"</DataArray>"
            write(unit_vtu,'(a)')"</PointData>"
            write(unit_vtu,'(a)')"<Cells>"
            write(unit_vtu,'(a)')"<DataArray type='Int32' Name='connectivity' format='ascii'>"
            do i = 1,nbody
                write(unit_vtu,'(i0)')i-1
            end do
            write(unit_vtu,'(a)')"</DataArray>"

            write(unit_vtu,'(a)')"<DataArray type='Int32' Name='offsets' format='ascii'>"
            do i = 1,nbody
                write(unit_vtu,'(i0)')i
            end do
            write(unit_vtu,'(a)')"</DataArray>"

            write(unit_vtu,'(a)')"<DataArray type='Int32' Name='types' format='ascii'>"
            do i = 1,nbody
                write(unit_vtu,'(i0)')1
            end do
            write(unit_vtu,'(a)')"</DataArray>"

            write(unit_vtu,'(a)')"</Cells>"
            write(unit_vtu,'(a)')"</Piece>"
            write(unit_vtu,'(a)')"</UnstructuredGrid>"
            write(unit_vtu,'(a)')"</VTKFile>"
            close(unit_vtu)
        end subroutine write_vtk
    end module output

40行目「write(unit_vtu,'(f8.3)',advance='no')ct(j)」は書式指定付出力である．
「'(f8.3)'」はct(j)を文字幅8，小数点以下3桁の実数で印字せよの意味である．
つまり「_:_-1.263」「_:_:_0.850」のような数が印字される．
これまで数字を出力する際にこの書式のことは気にせず「*」としてきた．
この場合は適当な書式（桁数等）がデータの型に応じて自動的に選択される．
倍精度実数を使って計算をしていても図のプロットやアニメの作成にはそれほど多くの桁数は必要ない．
出力ファイルが大きくなるだけである．

.. note::
    advance='no'は印字後に改行するなという命令．

.. note::
    二次元の計算しかしていないのだが，3次元の座標データがないとParaViewがうまく使えなかった．34-40行目あたりでz座標0.0を無理矢理書き入れている．

このサブルーチンをメインループの中で時々呼び出す．

.. code-block:: fortran

        if(mod(i-1,paraviewout) .eq. 0) then
            call write_vtk(i,coord,veloc,mass)
        end if

paraviewoutは出力間隔である．
呼び出した瞬間の位置情報などがvtuファイルに書きこまれる．
出力数に応じた複数のvtuファイルが生成される．

ParaViewでアニメ作成
===================================

全く使いこなせていない．使用方法はwebで検索してなんとかしよう．
以下，メモ

+ Point Spliteで粒子を描きたければ，Tools→Manage Plugins→PointSplite_pluginをLoadしておく必要がある．Auto Loadにしておけば次回起動時はLoad不要．
+ File→Openで生成したvtuファイルを選択．複数ファイルを一括で選択できる．
+ Sources→Axesで軸を追加．
+ Tools→Lock View Size Customで描画領域を設定．Youtubeに投稿するために1280×720に設定したが，Save animationで893×543に調整されてしまう？

Paraviewの使用画面は以下のような感じ．

.. image:: ../img/340paraview.PNG
   :scale: 30 %
   :alt: Snapshot of ParaView

作成した動画をYoutubeに掲載した．
`>>Youtube1 <http://youtu.be/YpCpxvlTpo4>`_

粒子の軌跡を追加する場合

+ File→Openで生成したvtuファイルを選択
+ PropertiesのRepresentationはPointにする（Point Spliteは使わない）．
+ PointにしてからApplyを押すとFiltersが使用可能になる．
+ Filters→Alphabetical→Temporal Particles To Pathlinesを選択する
+ Pipeline Browserに追加されたPathlinesの目玉マークを押す．
+ PropertiesのMask Pointsを1に．1粒子ごとに（つまり全ての粒子の）軌跡を描くの意味．10だと10粒子ごとに1本の線になる．
+ Max Track Lengthは250位に，ながければ長い軌跡になる．
+ これで再生すると軌跡だけが書けるはず．Particlesはよくわからないので，描画しない（目玉マークを押さない）

きれいな粒子も描きたければ．．．

+ もう一度データを読み込む．(File→Openで生成したvtuファイルを選択)
+ 今度はPropertiesのRepresentationをPoint Spliteにする．

以下のようなパイプラインになった．

.. image:: ../img/341pipelines.PNG
   :scale: 50 %
   :alt: Snapshot of pipelies of ParaView

作成した動画をYoutubeに掲載した．
`>>Youtube2 <http://youtu.be/3J71tEJS8wc>`_

月が地球の軌道の外側と内側を交互に移動しながら周回するため2番目の軌跡の色が縞状になっている．