断面積の取得

断面積を計算してみる

とりあえずデータがほしい人向けに。

やること

Janevの論文の式・データを使って、断面積(H原子の電子衝突電離)のエネルギー依存性の数値ファイルを得る。

参考文献

R.K. Janev, D. Reiter, U. Samm, "Collision Processes in Low-Temperature Hydrogen Plasmas", ATOMIC AND MOLECULAR PHYSICS S74, JUEL--4105(2003)

セットアップ

• 今回は電子衝突電離(低励起)とする
• 論文の該当ページを探す
  • 原理 11p
  • 係数データ Table 3, 130p
  • 計算したブラフ Figure 5, 164p
• 具体的に溶かしてみる

コード

• 反応閾値を threshold.f90 にメモ
• 各係数を reac001.dat 〜 reac005.dat (1s, 2s, 2p, n=2, n=3) にメモ
• エネルギー最小値 $10^{-3}$ eV から $10^4$ eV まで間隔 0.1 eVで計算
• 外部datファイル読み込み => E計算 => $\sigma$ 計算 => 出力

反応毎に分類したかったので、そういうコードになってます。

PROGRAM TEST
    IMPLICIT NONE
    
    INTEGER, PARAMETER :: NUM_REAC = 5      ! 考慮する反応の種類
    INTEGER, PARAMETER :: NUM_Esig = 100000 ! 断面積のエネルギーメッシュ数
    REAL(8) :: In(1:NUM_REAC)               ! スレッショルドエネルギー
    REAL(8) :: A(1:NUM_REAC,0:5)            ! 断面積を計算する係数
    REAL(8) :: E(0:NUM_Esig), dE=0.1d0      ! 考慮するエネルギー
    REAL(8) :: sig(1:NUM_REAC, 0:NUM_Esig)  ! 断面積 (cm^2)
    CHARACTER(len=128) :: reac_file
    
    INTEGER :: i, j, k
    E(0) = 1.d-3
    sig(:,:) = 0.d0
	
	
    ! Threshold Energy 
    OPEN( UNIT=100, FILE="threshold.dat", STATUS="OLD" )
        DO i = 1, NUM_REAC
            READ(100,*) In(i)
        END DO
    CLOSE(100)
    
    
    ! A coef
    DO i = 1, NUM_REAC     ! REAC NUMBER
        WRITE(reac_file, '("reac", i3.3, ".dat")'), i
        OPEN( UNIT=100, FILE=reac_file, STATUS="OLD" )
            DO j = 0, 5
                READ(100,*) A(i,j)
            END DO
        CLOSE(100)
    END DO
	
	
    ! Energy
    DO i = 1, NUM_Esig
        E(i) = E(i-1) + dE
    END DO
	
	
    ! Cross-section calculation
    DO i = 1, NUM_REAC
        DO j = 0, NUM_Esig
    	    DO k = 1,5         ! 係数の行数
    	        sig(i,j) = sig(i,j) + A(i,k)*(1.d0-In(i)/E(j))**k
    	    END DO
    	sig(i,j) = sig(i,j) + A(i,0)*log(E(j)/In(i))
    	sig(i,j) = sig(i,j) * 1d-13 / ( In(i) * E(j) )
        IF ( E(j) .lt. In(i) ) sig(i,j) = 0.d0
    	END DO
    END DO
	
	
    ! OUTPUT OF Cross-section
    DO i = 1, NUM_REAC
        DO j = 0, NUM_Esig
            WRITE(100+i, '(2ES16.3)') E(j), sig(i,j)
        END DO
    END DO
END PROGRAM TEST

コンパイルと実行

ifort test.f90
./a.out

結果

Janevらの論文のfigure 5 と同じような結果が得られた。