% This is a final document for black accretion 1D calculation
% === HISTORY ===
% ver1.0 : 2017.04.22 the first draft, S.Yamada

\documentclass{article} % 図をコメントアウトしない
% \documentclass[draft]{article} % 図と表をコメントアウトしてコンパイルしたい場合．

\usepackage{float}
\usepackage{color}
%%%%%%%% hyperref %%%%%%%%%%%%%%%%%%%
\usepackage[dvipdfm,
colorlinks=true,
bookmarks=true,
bookmarksnumbered=false,
bookmarkstype=toc]{hyperref}
%%%%%%%%%%%%%%%%%%%%%%%%%%%  

\usepackage{listings} % print out code 
\usepackage{setspace}

\pagestyle{empty}
\usepackage{array}
% \usepackage[dvips]{graphics} % 図をコメントアウトしない
% \usepackage[draft]{graphics} %  図をコメントアウトしてコンパイルしたい場合．
\usepackage{ascmac}

\setlength{\oddsidemargin}{0cm}
\setlength{\evensidemargin}{0cm}
\setlength{\textwidth}{16cm}
\setlength{\topmargin}{-1.5cm}
\setlength{\textheight}{24cm}

%\usepackage[dvips]{graphicx}
\usepackage[dvipdfmx]{graphicx}
%\usepackage{mediabb} % if pdf figures are included.
\usepackage{lscape}

\newcommand{\mfw}{0.5\textwidth}
\newcommand{\sfw}{0.3\textwidth}
\newcommand{\red}[1]{\textcolor{red}{\bf #1}}

\def\figdir{./fig}

\begin{document}
%
% ヘッダ部
% ----------------
\begin{screen}
\rule{0cm}{1.5cm}
{\Huge\bf \hspace{0.5cm} ブラックホール降着流の計算方法}\\
\rule{0cm}{1.5cm}
%\hspace{8cm}
\begin{tabular}{l}
S.Yamada (Tokyo Metropolitan University) \\ 
-- History -- \\
ver 1.0, 2017/04/22, S.Yamada, 1st draft, 
\end{tabular}
\end{screen}
% ----------------

\tableofcontents
% \listoffigures % YT commented out
% \listoftables % YT commented out

\clearpage

%===== section 1 (start) =====

\section{はじめに}

2017年3月27日(月)に、金沢大学附属高校の渡會氏から、降着円盤の計算コードを伝授され、それをまとめたもの。
元々は、千葉大学の松元亮治先生から始まった計算コードで、それを色んな方々で引き継いで使われてきた。
主には山田の個人的な備忘録の為にまとめるのが目的であるが、誰かの為にもなれば幸いである。

\section{使い方}

\subsection{ファイルの種類}

\begin{itemize}
\item accdisk.f　：メインのプログラム
\item at.job　　：インプットのパラメータファイル
\item para.h　　：配列や必要なパラメータを設定しているファイル
\item disk　　　：コンパイルするためのファイル
\item c0.out　　：時間発展の初期条件のファイル
\end{itemize}

\subsection{コンパイル方法}

gcc の 3 系でコンパイルしないと、コンパイルが通っても、正常に動作しない。gcc4 系と gcc5 系では正常動作しない。
Ubuntu を使うと、gcc の 3 系が古くて素直には入らない可能性が高い。Linux であれば、CentOSの6.7, 6.8 であれば動作が確認できている。
gcc 3 系の導入は次のように行う。

まず、fortran のライブラリを検索する。
\begin{lstlisting}[basicstyle=\ttfamily\footnotesize, frame=single]
> yum search fortran
....
============================= N/S Matched: fortran =============================
compat-libgfortran-41.i686 : Compatibility Fortran 95 runtime library version
                           : 4.1.2
compat-libgfortran-41.x86_64 : Compatibility Fortran 95 runtime library version
                             : 4.1.2
.... 
compat-gcc-34-g77.x86_64 : Fortran 77 support for compatibility compiler
.... 
\end{lstlisting}

ここで、gcc-34 のライブラリも現れるはず。それを確認したら、

\begin{lstlisting}[basicstyle=\ttfamily\footnotesize, frame=single]
> sudo yum install compat-gcc-34-g77.x86_64
\end{lstlisting}

とすると、インストールされる。そのあとに、

\begin{lstlisting}[basicstyle=\ttfamily\footnotesize, frame=single]

> f77 -v
Reading specs from /usr/lib/gcc/x86_64-redhat-linux/3.4.6/specs
Configured with: ../configure --prefix=/usr --mandir=/usr/share/man 
--infodir=/usr/share/info --enable-shared --enable-threads=posix 
--disable-checking --with-system-zlib --enable-__cxa_atexit 
--disable-libunwind-exceptions --enable-languages=c,c++,f77 
--disable-libgcj --host=x86_64-redhat-linux
Thread model: posix
gcc version 3.4.6 20060404 (Red Hat 3.4.6-19.el6)
\end{lstlisting}
と表示されて、バージョンが gcc version 3.4.6 であれば動作する。

コンパイル方法は、
\begin{lstlisting}[basicstyle=\ttfamily\footnotesize, frame=single]
> f77 -o accdisk accdisk.f
\end{lstlisting}
のように、f77 で必ずコンパイルする必要がある。

もし、gfortran でコンパイルすると、
\begin{lstlisting}[basicstyle=\ttfamily\footnotesize, frame=single]
> gfortran -o accdisk accdisk.f
accdisk.f:58.72:

      pause                                                             
                                                                        1
Warning: Deleted feature: PAUSE statement at (1)
accdisk.f:69.72:

            pause                                                       
                                                                        1

\end{lstlisting}
のようなコンパイルのワーニングがでるだけでなく、なぜかsegmentation fault を起こして止まる。

\subsection{走らせ方}

\subsubsection{定常解} 
at.job の
\begin{itemize}
\item KI = 0  定常解を解く
\item KR = 0 標準降着円盤、KR = 1 ADAF
\end{itemize}
に設定し、走らせる。

quick に結果を見たい場合は、IMORE = 0 にして、
標準出力にすると、定常解の結果だけが出力されたファイルが生成。
gnuplot で u 2:5 で、半径 vs 温度の絵ができるので確認する。

\subsubsection{時間依存} 
KI = 1 にして、時間発展をとく。
ADAF の場合は、ROUT2 を 50や100 にして、その程度から解くようにする。


\subsection{実行結果}
\subsection{標準降着円盤}
\subsection{ADAF}

\section{コードの概要}
\subsection{コードの流れ}

\subsubsection{SUBROUTINE cnvrbs}

この中で、xkff を用いて、tauff を２回計算している。
\begin{lstlisting}[basicstyle=\ttfamily\footnotesize, frame=single]
grep -n tau ./accdisk.f   | grep write
1552:      write(6,*) "tau1, tauff, taues, xkff", tau, tauff, taues, xkff
1600:      write(6,*) "tau3, tauff, taues, xkff", tau, tauff, taues, xkff
\end{lstlisting}
のように、プリントアウトをさせると、
\begin{lstlisting}[basicstyle=\ttfamily\footnotesize, frame=single]
 LIN =   1.83693930435181E+00     R =   2.6912097E+00     IC =     0
 tau1, tauff, taues, xkff  23.9289483  8.89031566E-311  23.9289483 1.26320111E-312
 J, R, taup, te 1  2.7012097  197.164335  794624.689
 tau3, tauff, taues, xkff  129.613409  105.684461  23.9289483  1.50164214
 J, R, taup, te 1  2.7012097  197.164335  794624.689
 
 tau1, tauff, taues, xkff  136.4634  111.269821  25.1935788  1.50164214
 J, R, taup, te 2  2.7112097  199.16652  811724.998
 tau3, tauff, taues, xkff  131.951256  106.757677  25.1935788  1.44074848
 J, R, taup, te 2  2.7112097  199.16652  811724.998
\end{lstlisting}
のように、一番最初は、xkff $=  8.89031566E-311$ というものすごい小さい数字が入っている。
その直後に、圧力と密度から、温度を正しく計算する。
２回目以降は、一つ前の xkff を使って、最初も計算するので、変なことにはならない。
最初だけ、初期化されてない変数 xkff を使っているために、変なことが起こる。


\subsection{コーディングの詳細}

\subsection{配列}

変数をいれる配列は、ヘッダーファイル para.h の中に、長さが  parameter (jmax=500) で定義されており、

\begin{lstlisting}[basicstyle=\ttfamily\footnotesize, frame=single]
      parameter (jmax=500)
      common /BMESH/ RIN2,ROUT1,ROUT2,X0(jmax),X1(jmax),Y0(2,jmax),JX
      common /BPHY0/ W1(jmax,4),W0(jmax,4),U0(jmax,4)
      common /optic/ tauabso(jmax),tauo(jmax)
      common /shoki/ deini(jmax),vxini(jmax),vpini(jmax),teini(jmax)
      common /cauti/ krs(jmax),kes(jmax)
      common /coula/ nlmax,nl(jmax)
      common /energ/ es1(jmax)
      common /fluxr0/ drepsi(jmax),dgepsi(jmax),repsi(jmax),gepsi(jmax)
      common /flux1/ tepsi(jmax),dtepsi(jmax)
\end{lstlisting}
のように、ほとんどの変数が、jmax の長さ分だけ定義されている。
実際に確保される変数の長さは、at.job の中の MESH という変数で、これはデフォルトで 200 に
設定されている。

fortran で、配列の中身をプリントアウトしたい場合は、
\begin{lstlisting}[basicstyle=\ttfamily\footnotesize, frame=single]
 write(*,*) x0
\end{lstlisting}
とすると、全部出力される。

実際に、コンパイルして、x0 の中身を出力すると、
\begin{lstlisting}[basicstyle=\ttfamily\footnotesize, frame=single]
 end steady to test. x0 = 
  2000.  1935.0031  1872.11849  1811.27754  1752.41383  1695.46309  1640.36317
  1587.0539  1535.47711  1485.57648  1437.29755  1390.5876  1345.39566
  1301.67238  1259.37005  1218.44247  1178.84498  1140.53434  1103.46874
.....
  2.8012097  2.7912097  2.7812097  2.7712097  2.7612097  2.7512097  2.7412097
  2.7312097  2.7212097  2.7112097  2.7012097  2.6912097  0.  0.  0.  0.  0.  0.
  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
....
  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
\end{lstlisting}
のようになっており、$r$ の大きい数字から配列に保存されており、
確保された500個の配列要素のうち、200個(=MESH数)まで値が詰まっている。


\subsubsection{tauff と taup の違い}

\begin{lstlisting}[basicstyle=\ttfamily\footnotesize, frame=single]
      ratio = TAUFF/taup
      print *,"SADMN : TAUES=",TAUES," TAUFF=", TAUFF, " taup=", taup
      print *,"SADMN : ratio=", ratio
\end{lstlisting}
とコメントを入れて、表示させると、

\begin{lstlisting}[basicstyle=\ttfamily\footnotesize, frame=single]
.....
SADMN : TAUES=  5.09422324E-09 TAUFF=  5.97688026E-67 taup=  1.11504341E-66
SADMN : ratio=  0.536022202
SADMN : TAUES=  7.71452468E-06 TAUFF=  7.20877975E-48 taup=  1.34486589E-47
SADMN : ratio=  0.536022202
SADMN : TAUES=  0.00187276193 TAUFF=  1.47537411E-33 taup=  2.75244962E-33
SADMN : ratio=  0.536022202
.....
\end{lstlisting}
となるので、0.536022202 だけの違いなので、どこかで係数が約２倍ほど定義が違うだけだと思われる。



\subsection{注意点}

\subsubsection{大文字と小文字}

fortran は、大文字と小文字を区別しないので、
ソースコードを調べるときには、grep にオプション i をつけて、
\begin{lstlisting}[basicstyle=\ttfamily\footnotesize, frame=single]
grep -i xkff ./accdisk.f
\end{lstlisting}
としないと、出現している変数を調べきれない。

\section{基礎方程式}

\subsection{標準降着円盤モデル}

標準円盤モデルは、以下の近似および簡単化をしている。

\begin{enumerate}
\item 円筒座標系($r,\phi, z$) を用いる
\item 円盤は定常である : $\partial / \partial t \Longrightarrow 0 $
\item 円盤は軸対称である : $\partial / \partial \phi \Longrightarrow 0 $
\item 円盤は幾何学的に薄い。$H \ll r $ : ($r, z$) 2次元 $\Longrightarrow $ (1 + 1) 次元化
\item 円盤はケプラー回転 (遠心力 = 重力) する
\item 円盤は光学的に厚い : $\tau \gg 1$
\item 円盤は黒体放射する
\item 円盤表面単位面積からの単位時間当たりの放射エネルギー量は、放射拡散近似で、
\begin{equation}
F~[\rm{erg/cm}^{-2}\rm{s}^{-1}]= \int (ac/3\kappa \rho) (dT^4/dz) dx \approx 4 \sigma T^4 / 3 \tau 
\end{equation}
と書ける。ここで、
\begin{itemize}
\item 輻射密度定数 $a =  \pi^2 k^4 / 15 c^3 \hbar^3 = 7.57 \times 10^{-15}$ erg/ cm$^{-3}$ K$^{-4}$, 
\item シュテファン--ボルツマン定数 $\sigma =  ac/4 = 5.67 \times 10^{-5}$ erg cm$^{-2}$ K$^{-4}$ s$^{-1}$
\end{itemize}
である。
\item 円盤において、粘性加熱率と放射冷却率がバランスする。
\end{enumerate}

まとめると、
\begin{itemize}
\item 未知の物理量は 10 個 ($\rho, \Sigma, H, v_r, c_s, \Omega, P, T, \tau, \nu$) 
\item 基本方程式は 10 本
\item パラメータは、$M, \dot{M}, \alpha$ の３つ。
\end{itemize}
で式が閉じて、各物理量が $r$ の関数として決定できる。

\begingroup
\renewcommand{\arraystretch}{1.8}
\begin{table}[htb]
\caption{$\alpha$ 円盤モデルの基本方程式}
\begin{center}
\begin{tabular}{| l  l |}
\hline
1. $\rho = \Sigma / 2H$ & 連続の式 (鉛直方向) \\ 
2. $\dot{M} = - 2 \pi r v_r \Sigma $ & 連続の式 (動径方向) \\ 
3. $H = c_s / \Omega $ & 静水圧平衡($z$方向重力=圧力勾配力) \\ 
4. $\Omega = \sqrt{GM/r^3}$ & ケプラー回転 (動径方向の運動量の式) \\ 
5. $\nu \Sigma = \frac{\dot{M}}{3 \pi} \left( 1 - \sqrt{\frac{r_*}{r}} \right)$ & 角運動量保存 \\
6. $\frac{4 \sigma T^4}{ 3 \tau} (=F)= \frac{9}{8} \nu \Sigma \Omega^2$ & エネルギーバランス \\ 
7. $c_s^2 = P/\rho$ & 音速 \\
8. $P = P_{\rm{gas}} +  P_{\rm{rad}}  = \frac{\rho k T}{ \mu m_p} + \frac{4\sigma}{3c} T^4$ & 状態方程式 \\ 
9. $\tau = (\kappa_{abs} + \kappa_{es}) \rho H = (\kappa_{o} \rho T^{-3.5} + \kappa_{es} ) \rho H $  & 光学的厚み \\
10. $\nu = \frac{2}{3} \alpha c_s H$ & $\alpha$ 粘性 \\
\hline
\end{tabular}
\end{center}
\label{default}
\end{table}%
\endgroup


基本方程式は、一見すると複雑に見えるが、圧力と吸収係数を覗いて、掛け算と割り算のみで表現されている。
そこで、圧力は、ガス圧 $P_{\rm{gas}}$ か輻射圧 $P_{\rm{rad}}$、
free-free による吸収 $\kappa_{abs}$ か電子散乱 $\kappa_{es}$ か、それぞれどちらが支配的になるかで場合分けを
すると、領域ごとに、すべて掛け算と割り算で表現できる。そうすると、すべての式を対数をとると、
掛け算と割り算が、足し算と引き算で表現できて、代数方程式になる。最終的な結果は、次のようになる。

３つの領域に場合分けをする。
\begin{itemize}
\item 領域 a : $P \sim P_{rad}, \kappa \sim \kappa_{es}, r < r_{ab} $
\item 領域 b : $P \sim P_{gas}, \kappa \sim \kappa_{es}, r_{ab} < r < r_{bc} $
\item 領域 c : $P \sim P_{gas}, \kappa \sim \kappa_{abs}, r > r_{bc} $
\end{itemize}

定数を下記のように定義する。
\begin{eqnarray}
R           & \equiv & r/r_s \\
m           & \equiv & M/M_{solar} \\
\dot{m}   & \equiv & \dot{M}/ (L_E/c^2) \\
f              & \equiv & 1 - \sqrt{3 r_s/ r}\\
r_{ab} & \sim    & 18 \alpha^{2/21} m^{2/21} \dot{m}^{16/21} f^{16/21} r_s \\ 
r_{bc} & \sim    & 2500 \dot{m}^{2/3} f^{2/3} r_s 
\end{eqnarray}




\begingroup
\renewcommand{\arraystretch}{1.8}
\begin{table}[htb]
\caption{$\alpha$ 円盤モデルの内縁部の解 (領域a)}
\begin{center}
\begin{tabular}{| l  l |}
\hline
変数 & 単位 \\
\hline
$ H =  5.5 \times 10^4 m \dot{m} f                                                  $ & cm \\ 
$ \Sigma = 1.0 \times 10^2 \alpha^{-1} \dot{m}^{-1} R^{3/2} f^{-1} $ & g cm$^{-2}$ \\ 
$ | v_r | = 7.6 \times 10^8 \alpha \dot{m}^{2} R^{-5/2} f                  $ & cm s$^{-1}$ \\ 
$ T = 4.9 \times 10^7 \alpha^{-1/4} m^{-1/4} R^{-3/8}            $ & K \\ 
$ \tau_* = 8.4 \times 10^{-3} \alpha^{-17/16} m^{-1/16} \dot{m}^{-2} R^{93/32} f^{-2}             $ & -- \\ 
\hline
\end{tabular}
\end{center}
\label{default}
\end{table}%
\endgroup

\begingroup
\renewcommand{\arraystretch}{1.8}
\begin{table}[htb]
\caption{$\alpha$ 円盤モデルの中間部の解 (領域b)}
\begin{center}
\begin{tabular}{| l  l |}
\hline
変数 & 単位 \\
\hline
$ H =  2.7 \times 10^3 \alpha^{-1/10} m^{9/10} \dot{m}^{1/5} R^{21/20} f^{1/5}                                                  $ & cm \\ 
$ \Sigma = 4.3 \times 10^4 \alpha^{-4/5} m^{1/5}  \dot{m}^{3/5} R^{-3/5} f^{3/5}                                                $ & g cm$^{-2}$ \\ 
$ | v_r | = 1.7 \times 10^6 \alpha^{4/5} m^{-1/5} \dot{m}^{2/5} R^{-2/5} f^{-3/5}                                                  $ & cm s$^{-1}$ \\ 
$ T = 2.2 \times 10^8 \alpha^{-1/5} m^{-1/5} \dot{m}^{2/5} R^{-9/10} f^{2/5}                                                       $ & K \\ 
$ \tau_* = 2.4 \times 10^{1} \alpha^{-4/5} m^{1/5} \dot{m}^{1/10} R^{3/20} f^{1/10}             $ & -- \\ 
\hline
\end{tabular}
\end{center}
\label{default}
\end{table}%
\endgroup


\begingroup
\renewcommand{\arraystretch}{1.8}
\begin{table}[htb]
\caption{$\alpha$ 円盤モデルの外縁部の解 (領域c)}
\begin{center}
\begin{tabular}{| l  l |}
\hline
変数 & 単位 \\
\hline
$ H =  1.5 \times 10^3 \alpha^{-1/10} m^{9/10} \dot{m}^{3/20} R^{9/8} f^{3/20}                                                  $ & cm \\ 
$ \Sigma = 1.4 \times 10^5 \alpha^{-4/5} m^{1/5}  \dot{m}^{7/10} R^{-3/4} f^{7/10}                                                $ & g cm$^{-2}$ \\ 
$ | v_r | = 5.4 \times 10^5 \alpha^{4/5} m^{-1/5} \dot{m}^{3/10} R^{-1/4} f^{-7/10}                                                  $ & cm s$^{-1}$ \\ 
$ T = 6.9 \times 10^7 \alpha^{-1/5} m^{-1/5} \dot{m}^{3/10} R^{-3/4} f^{3/10}                                                       $ & K \\ 
$ \tau_* = 7.9 \times 10^{1} \alpha^{-4/5} m^{1/5} \dot{m}^{1/5} f^{1/5}             $ & -- \\ 
\hline
\end{tabular}
\end{center}
\label{default}
\end{table}%
\endgroup




\clearpage

\subsection{高温降着流 ADAF}

ADAFモデルは、以下の近似および簡単化をしている。

\begin{enumerate}
\item 円筒座標系($r,\phi, z$) を用いる
\item 円盤は定常である : $\partial / \partial t \Longrightarrow 0 $
\item 円盤は軸対称である : $\partial / \partial \phi \Longrightarrow 0 $
\item 円盤は幾何学的にやや厚めである。$H < r $ : ($r, z$) 2次元 $\Longrightarrow $ (1 + 1) 次元化
\item 円盤は非ケプラー回転 (遠心力 $ < $ 重力) する
\item 円盤は光学的に薄い : $\tau < 1$
\item 円盤は熱制動放射や逆コンプトン散乱, シンクロトロン放射で光る。
\item 円盤において、粘性加熱率と移流冷却率がバランスする。
\end{enumerate}


まとめると、
\begin{itemize}
\item 未知の物理量は 10 個 ($\rho, \Sigma, H, v_r, c_s, \Omega, P, T, \tau, \nu$) 
\item 基本方程式は 10 本
\item パラメータは、$M, \dot{M}, \alpha, l_{in}$ の４つ。(アウトフローがある場合はさらに $r_{disk}$ と $p$が加わる)
\end{itemize}

ここで、$l_{in}$は、ブラックホールに吸い込まれるガスのもつ、単位質量当たりの角運動量を表す。標準円盤との大きな違いは、
運動方程式(4)とエネルギー方程式(6)に現れる。ともに、微分を含んでいることが解析を困難にし、数値的に解く必要がある。
適当な境界条件を与えると、微分方程式が解けて、密度や温度などの物理量が $r$ の関数として決定できる。


\begingroup
\renewcommand{\arraystretch}{1.8}
\begin{table}[htb]
\caption{ADAFモデルの基本方程式}
\begin{center}
\begin{tabular}{| l  l |}
\hline
1. $\rho = \Sigma / 2H$ & 連続の式 (鉛直方向) \\ 
2. $\dot{M} = - 2 \pi r v_r \Sigma = \dot{M}_{input} (r/r_{disk})^p$ & 連続の式 (動径方向) \\ 
3. $H = c_s / \Omega $ & 静水圧平衡($z$方向重力=圧力勾配力) \\ 
4. $v_r \frac{d v_r}{dr}  = - \frac{1}{\rho}\frac{dP}{dr} + r \left( \Omega^2 -  GM/r^3 \right)$ & 動径方向運動量 \\ 
5. $\nu \Sigma = \frac{\dot{M}}{3 \pi} \left( 1 - \sqrt{\frac{l_{in}}{ \sqrt{GMr}  }} \right)$ & 角運動量保存 \\
6. $\frac{1}{2} v_r \Sigma T \frac{ds(\rho,T)}{dr} = \frac{9}{8} \nu \Sigma \Omega^2  - F(\rho,T) $ & エネルギー保存\\ 
7. $c_s^2 = P/\rho$ & 音速 \\
8. $P = P_{\rm{gas}} = \frac{\rho k T}{ \mu m_p} $ & 状態方程式 \\ 
9. $\tau = (\kappa_{abs} + \kappa_{es}) \rho H = (\kappa_{o} \rho T^{-3.5} + \kappa_{es} ) \rho H $  & 光学的厚み \\
10. $\nu = \frac{2}{3} \alpha c_s H$ & $\alpha$ 粘性 \\
\hline
\end{tabular}
\end{center}
\label{default}
\end{table}%
\endgroup

ここで、$F(\rho,T)$は、放射冷却率のことである。



\subsection{時間依存の方程式}

%%% YT Figure~1 の場所を移動しました．
%\begin{figure}[htbp]
%\begin{center}
%\includegraphics[angle=0,width=0.7\textwidth]{\figdirOverview/sxs_photo.pdf}
%\caption{Photo of the SXS.}
%\label{sxs-overview:fig:sxsphoto} 
%\end{center}
%\end{figure}


\newpage

\section{コードの詳細}

\subsection{定常解}
\begin{figure}[htbp]
\begin{center}
\includegraphics[angle=0,width=0.9\textwidth]{omnigraffle_bh/doc_calc_blackhole_accretion_20170520.pdf}
\caption{定常部分のコードのフローチャート}
\label{fig:block} 
\end{center}
\end{figure}




\newpage

\section{参考文献}

定常解に関しては、
\begin{itemize}
\item ブラックホール天文学 (新天文学ライブラリー), 嶺重 慎 
\end{itemize}
の記述に沿ってまとめた。

\end{document}