タイトルとはまるで関係のない話、「放射性崩壊」から始めてみようと思います。にわかにして、巷にあふれるようになった放射能とか放射性崩壊といった言葉。このブログの読者の中にも、いろんな立場の方がいらっしゃることと思いますし、出来事の背景や是非をここで議論するつもりはありません。ただ、かつて高校で「放射性崩壊」について学んだ記憶が懐かしく思い出されたのは、まぎれもない事実でした。

#だって、平々凡々とした日常生活では、思い起こす機会がないですから、この知識。

放射性崩壊にはアルファ崩壊、ベータ崩壊、ガンマ崩壊の3種類があるというのが、僕の高校物理の記憶です。(その後、大学で聞いた話では、ベータ崩壊には、β+崩壊とβ-崩壊があるそうなんですが、残念ながら僕は詳しくないので、ここでは高校で習ったレベルで勘弁を。)

  • アルファ崩壊:原子核からアルファ粒子(=ヘリウム 4He2+ の原子核)が飛び出してくるもの。崩壊した原子核は、He の分だけ、つまり、質量数が4(中性子が2個、陽子が2個)減ります。
  • ベータ崩壊:中性子 n が陽子 p+に変化し、その時、電子 e- を放出します。中性子が陽子に変化するので、原子番号は一つ増加しますが、質量数は変化しません。
  • ガンマ崩壊:陽子や中性子の変化を伴いませんので質量数は変化しません。原子核内部で励起された電子がどこかのエネルギー順位に落ちるときに、それに等しいエネルギー(ガンマ線)が放出されます。

電離作用は、アルファ粒子>ベータ粒子>ガンマ線の順に強いので、曝露した際の影響はアルファ粒子が大きいのですが、一方で、透過力はアルファ粒子が一番弱くて、紙一枚透過できないほど(なんだそうです)。

もう一つよく耳にする言葉に「半減期」というのがあります。文字通り半分になるまでにかかる時間です。すっかり有名になった 131I(ヨウ素131)の半減期はおよそ8日。数時間とか数日だったら計測できる気がしますよね。

でも、238U(ウラン238)の半減期は44億年といわれています。到底測れるわけないです。これはどうやって導かれるのでしょう。

▼ 半減期

もちろん、44億年なんて計測できるわけありません。単純な系を考えれば、次のように計算できます。ポイントは、いずれの崩壊も原子核単独で起こっているということ。このことから、存在する原子核の数(N)とすると、単位時間に減る量(=崩壊の速度)(dN/dt)は N に比例すると考えることができますので、その比例定数(崩壊速度定数っていうのかな)を k (k>0) とすると、

-dN/dt = kN (k>0)

と原子核数 N の1次反応として表現することができます。1次反応の微分方程式を IGOR Pro で数値的に解くお話は、これまでにも紹介していますので、ここでは触れません。もうちょっと解析的に話を進めます。

この式は解析的に解くことができて、

N = N0* Exp(-kt)    (1)

が得られますよね。ここで、N0 は積分定数ですが、物理的には最初に何かぼわっと排出された原子核の総量みたいなものですね。半減期は N が N0/2 となるまでに要する時間 t1/2 のことですので、これを求めると、

t1/2 = ln(2)/k

となります。何度見てもこの式は素敵で、量や濃度には依存せず、核種(k)にだけ依存するわけです。ものすごくすっきりしています。2分の1 に限らず一般化すると、

t = ln(N0/N)/k     (2)

と表現できそうです。例えば、131I の場合、半減期 t1/2 は 8.02 日ですから、k = ln(2)/8.02 ≒ 0.0864 とわかりますので、100分の1になるまでの時間は、t1/100 = ln(100)/0.0864 ≒ 53.3 日と求めることができますね。


▼ 逆関数を求めずに値を求める

任意の濃度を与える時間を得るには、式1 を式2 のように変形する必要があります。この操作は逆関数を求める操作ですね。(微分方程式を解析的に解くこともそれなりに大変なのですが、)逆関数 x = g(y) を求めることが、簡単じゃないことだってありますよね。ここでは、IGOR Pro で逆関数を求めずに、任意の y を与える x を求める方法を紹介したいと思います。(やっと本題)

IGOR Pro で関数 f(x) を定義するには、引数にパラメータウェーブと x変数を持つ次のような関数を用意します。

Function myFunc(w,x)
    Wave w
    Variable x

    Return f(w,x)

End

このブログではよく出てくるプロシージャですね。例えば、放射性崩壊の例であれば、

Function myFunc(w,x)
    Wave w
    Variable x

    Return w[0]*Exp(-w[1]*x)

End

とします。初期の原子核の数 N0 を、例えば単位を気にせず 1000、k を先ほど得られた 0.0864 とするには、

Make/o/n=2 coef={1000, 0.0864}

という2ポイントのウェーブ coef を作っておきます。

なんと準備はこれで完了。あとは FindRoots 操作関数にお任せです。

FindRoots/Z=y-value 関数名, 係数ウェーブ

初期値が1000だとすると、100分の1は 10 ですので、

FindRoots/Z=10 myFunc, coef

を実行すると

*FindRoots/Z=10 myFunc, coef
  	Looking for f(x) = 10
  	Possible root found at 53.3006
  	Y value there is 10

と出力が得られます。先ほど手計算した 53.3 日が得られます。


▼ もう一歩進めて

関数によっては /Z で指定した値をいくつも持つ関数があり得ますよね。例えば、Sinc 関数。

fig1.png

/Z=0(←0の場合は省略可能) を考えると、当然、解がわんさか見つかります。放射性崩壊の例は単調減少する関数ですので、範囲の指定は不要ですが、一般に、どの範囲で解を求めるのかという指定も重要になります。例えば x = 0~5 で解を探すのであれば、こんな風にします。

まずは、関数をプロシージャウィンドウに定義。

Function mySinc(w, x)
	Wave w
	Variable x

	return w[0]+w[1]*sinc(w[2]*(x-w[3]))
End

次に、係数ウェーブやグラフの準備

Make/D/O SincCoefs={0, 1, 2, .5}
Make/D/O SincWave
SetScale/I x -10,10,SincWave
SincWave = mySinc(SincCoefs, x)
Display SincWave
ModifyGraph zero(left)=1
ModifyGraph minor(bottom)=1

最後に、/L フラグで左側、/H フラグで右側を指定して実行。

Findroots/Z=0/L=0/H=5 mySinc, SincCoefs
  	Looking for two roots...

	Results for first root:
  	Possible root found at 2.0708
  	Y value there is -2.99483e-11
  
	Results for second root:
  	Possible root found at 3.64159
  	Y value there is 3.4303e-11

x = 2.07 と 3.64 あたりに解があることが分かります。グラフを見ると、さもありなんという感じですね。

fig2.png

▼ まとめ

「放射性崩壊」なんて話から始めましたが、逆関数を求めずに欲しい値を得るというお話でした。Mathematica だったら、もっとエレガントにやってのけるんだろうなぁと思いましたが、敢えてIGOR Pro でお届けしました。

--
井上@技術部でした。

※この記事の内容は執筆者の個人的見解で、ヒューリンクスによる公式情報ではありません。[免責事項]

トラックバック

この記事へのトラックバックURL
http://blog.hulinks.co.jp/cgi/mt/mt-tb.cgi/470
内容に対しての関連性がみられないものは削除する場合があります

コメントの投稿

Emailアドレスは表示されません。は必須項目です。
ヒューリンクス取り扱い製品の内容や購入に関するお問い合わせはヒューリンクスサイト連絡先へお願いいたします。投稿前にその他の注意事項もご覧ください。

HULINKS サイトの新着情報