1 / 52

授業の内容

授業の内容. 天文学は天体からの光を研究する学問です。 そこでこの授業では、「光」をどう扱うかの基礎を学びます。 授業計画は、 A.水素原子  B.エネルギー準位  C.熱平衡  D.線吸収  E.連続吸収   F.光のインテンシティ G.黒体輻射 H.等級  I.色等級図   J.光の伝達式 I  K.光の伝達式 II L.星のスペクトル という順で進めます。 最後まで行くと、星のスペクトルがどんな仕組みで決まっているかが判る、 というのが目標です。 AからEまでは光の吸収に関係する物理の話です。Fでは光の強さをきちん

gellert
Download Presentation

授業の内容

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript


  1. 授業の内容 天文学は天体からの光を研究する学問です。 そこでこの授業では、「光」をどう扱うかの基礎を学びます。 授業計画は、 A.水素原子  B.エネルギー準位  C.熱平衡  D.線吸収  E.連続吸収   F.光のインテンシティ G.黒体輻射 H.等級  I.色等級図   J.光の伝達式 I  K.光の伝達式 II L.星のスペクトル という順で進めます。 最後まで行くと、星のスペクトルがどんな仕組みで決まっているかが判る、 というのが目標です。 AからEまでは光の吸収に関係する物理の話です。Fでは光の強さをきちん と定義します。GからIは光の強さを天文学でどう使うかを示します。JからLは 光がガス中を伝わる様子を式に表わし、その式を解いて星のスペクトルを導き ます。それでは、始めましょう。

  2. D 線吸収 今回の内容 (D.1) 吸収断面積 σ   原子の吸収断面積 σ の意味。幾何学的な断面積との違いを学びます。   (D.2) 強制振動   摩擦により失われる振動エネルギーを強制外力による仕事が補って安定   振動になることを学びます。 (D.3) 双極子の振動   原子を電気双極子、光の電磁場を強制外力と考え、光から原子に注入され   るエネルギー率を求めます。それから原子の光吸収断面積 σ を求めます。  (D.4) 振動子強度 (= f - 値)   上で求めた吸収断面積は、吸収線の強度がみな同じになってしまいます。   吸収強度を表わす量としてf-値を導入します。 (D.5) 等値幅  吸収線の強さを観測的に表現する量です。 (D.6) Voight プロファイル   原子は運動しています。ドップラー効果を受けた原子吸収断面積がどうなる       かを調べます。様々な速度で運動する原子全体での平均吸収断面積の形を   Voight プロファイルと呼びます。

  3. 数学の準備 簡単な公式ですが、今回使う式をいくつか紹介します。 (1) を微分してみましょう。 (2) (3) (4)

  4. (4)の dy/dx, d2y/dx2が ia・y、-a2y になることを示して下さい 三角関数の公式も少し。

  5. dx σ σ S σ σ S σ σ σ D.1. 吸収断面積 σ と吸収係数 k σ:粒子断面積  n:粒子数密度 dx I(x)+dI I(x) σ ? dV=S・dx 正面(面積S)から見ると 総断面積 Σ Σ = σ n dV   = σ n S dx 被覆率 = C C = Σ/S = σ n dx

  6. エネルギー流 I の立場から見ると、 dx dI=-I・C =-I・κ・dx   =-I・σ・n・dx I I-IC エネルギー流の強度が弱っていく。 κ=吸収係数 σ =吸収断面積 κ=n・ σ  吸収係数 κ は単位体積中の総吸収断面積と考えればよい。 吸収粒子の立場から見ると、 σ =吸収断面積を持つ粒子は、強度 I の エネルギー流にさらされると、そこから I σ のエネルギーを吸収する。 I σ

  7. D.2.ばねの強制振動 ばねの運動を考えましょう。ばねの長さを z とし、ばね先端の点 Aの運動を 考えます。A点にはばねの長さ z に比例し、A点を中心方向に引き戻そうとする ばねの引力 -Kz が働いています。Aの質量をmとします。運動方程式は、 A -ω02z z 解 I は、 となります。 の解 II を とします。 この 解 II の実部が 解 I と一致することを確かめて下さい。

  8. 解に含まれる定数 A と B は、t=0の時の zと (dz/dt ) から決まります。 例えば、 の解で、t=0の時 z=5、dz・dt=0 だったら、 A=5, B=0 となり、解I では、 です。これはまた、解II の書き方では、 です。解 I と解II は異なる形で表現されていますが同じ解です。 この様に単純なばねのモデルでは、点Aが同じ周期、同じ 振幅でいつまでも振動します。 このような振動を単振動と言います。ω0は、固有角振動数と呼ばれます。

  9. バネのエネルギー で表わされるバネのエネルギー E を考えましょう。 =運動エネルギー E = K + U   =位置エネルギー  を使うと では、 m=1の時、今やった ですから、 単振動のエネルギー E は一定です。次ページのグラフで確認して下さい。

  10. バネのエネルギー E を一般の場合で扱うと、 単振動のエネルギー E は振幅Aの二乗に比例するのです。さっきの例では m=1、A=5, ωo=3 でしたから、 E=(1/2)・25・9=225/2 になって一致します。

  11. 摩擦のあるバネ 次に、振動子に摩擦が働いているとします。摩擦力は速度に比例し、速度と反対向きに働くと考え、ーD(dz/dt)  と表わします。 これを先ほどの単振動の式に加えると、 この式が摩擦入りの振動の式です。 この式を z = A exp[i (ω t+B)]  (A, B 実数)と仮定して解いてみましょう。 (1) z = A exp[i (ω t+B)]を上の運動方程式に代入し、ωを求めて下さい。    ただし、摩擦は小さく  γ < 2 ωoとします。 (2) 方程式の解は、 z = A exp[i (ω t+B)] の実数部分と考えます。     (1)で求めた ω を上の式に代入し、実数と虚数に分けて下さい。 (3) ω02 =9.25, γ=1 の解を求めて下さい。 (4) 初期値を t=0 で z=5, dz/dt = -2.5  としたとき、A と B はいくつですか (5) (4)を t = 0 – 8 でグラフにして下さい。

  12. 計算用I

  13. 計算用II

  14. 摩擦によるバネエネルギーの減衰 の解は、 この式は、cos(ω1t+B)という単振動の振幅が A・exp(-γt/2)で低下していくと読めます。  前に示したように、振幅が X のバネの総エネルギーEは (1/2)mX2ω2 です。 ですから、 Eoは t=0 の時振動子が持っていたエネルギーです。 上の式は          つまり、エネルギーEの振動子がエネルギーを 失う率は、 γE であることを示します。

  15. 外力による強制振動 問1(3)、(4)で見たように、単振動運動は摩擦力によって減衰していきます。 しかし、適当な振動外力(ω)が加わると、この外力が摩擦による減衰を丁度補うところで、安定な振動が可能になります。その振動を求めてみましょう。 摩擦力入りのばねの運動方程式は 外力を という形にしましょう。すると運動方程式は C/m = α とおくと、 安定な振動解を、z=A・ei ωt とおいて求めます。上式に代入して、問1 と同様に計算を進めると、

  16. 問1の(2)でやった時と同様に、この z の実部が求める解です。 強制振動の力のどれだけがばねに吸収されて、摩擦で消費されていくかはこのzを使って求めることができます。外力 F の仕事は F・dz ですから、仕事 率は F・(dz/dt) です。すると、強制振動がばねにする仕事率Wは、

  17. ですから、 sin(2ωt-δ) の項は平均すると消えてしまい、sin δ の項のみ残ります。 したがって、 これが前に求めた減衰率 γE と等しいことにも注意しておいて下さい。 通常問題にするのは、強制振動の各周波数 ω が固有角振動数ω0 に近い場合が多いので、 ω ≒ω0 の近似で上の式を少し変形します。 角振動数 ω よりも普通の振動数 ν の方がなじみやすいので、 ω = 2π ν を使って上の式を書き換えると、

  18. 強制振動の力が最も効率よくばねの振動に伝わるのは ω=ω0の時です。前の式を調べると、その時は δ =π/2 つまり、ばねの位置の位相が強制力と90度 ずれていることが判ります。しかし、その時には速度 dz/dt の位相が強制力と一致し、ばねが正方向に動くときには強制力が正方向に押してくれ、マイナス方向に動くときにはマイナスに押してくれるのです。 この様に強制力がばねの固有振動数と一致して働くことを「共鳴」と呼びます。 共鳴から離れると、強制振動の力がばねに伝わる効率は低下することも図から読み取れます。

  19. もう一度まとめると、 (1) 自然なままなら、exp (i ω0 t) で振動するばねがあるとします。  ω0のことを固有角振動数と呼びます。 (2) バネに-D(dz/dt)  という速度に比例する摩擦が働くと、ばねはエネ     ルギーを失い、次第に振幅が小さくなっていきます。 (3) しかしそこに外から C exp (i ω t) という力が働くと、ばねにエネルギーが     注入されます。 (4) ある振幅で外力からの注入エネルギーと摩擦により失われるエネ     ルギーが釣り合い、一定の振幅、一定の周波数で振動します。 (5) それが解 です。 この解を見ると幾つかの性質に気付きます。 (1) ばねは固有角振動数ω0ではなく、外力の角振動数ωで揺れます。 ばねが外力の振動に合わせて揺すられるので、「強制振動」と言います。 (2)  ω が ω0 に近づくと 強制振動の振幅が大きくなります。これを「共鳴」と     呼びます。 (3) ばねの振動は外力と同期していません。少しずれます。これを「位相シフ     ト」と呼びます。

  20. 「ぶらんこ」 共鳴のよい例はブランコでしょう。漕がない時のブランコは一定の周期で揺れています。これが「固有振動」です。 ブランコには摩擦が働くため、ブランコの揺れはだんだん小さくなっていきます。これが問1でやった「減衰振動」です。 ブランコを押してみましょう。これは「外力」にあたります。しかし、固有振動と 違うリズムで押されると、ブランコは仕方なく押される周期で行ったり来たりしますがあまり揺れません。これが「強制振動」です。 押すリズムが固有振動と合ってくるとブランコの揺れが大きくなります。揺れが最大になるのは、押す周期が固有振動と一致した時です。これが「共鳴」です。

  21. D.3.双極子としての原子 原子はプラスの電荷を持つ原子核とマイナス電荷の電子から成っています。 電子の分布の中心と原子核の位置が一致すれば外から見た原子は中性で す。 しかし、二つの位置がずれると原子は双極子として見えます。 双極子というのは、 左図のようにプラスの電荷とマイナスの電荷が距離 z 離れて並んでいるものの名前です。この時、二つの電荷の間に働く引力は距離に比例して大きくなります。 +q 原子核 z これはばねの力と、丁度同じ性質です。 ですから双極子の運動には C.2. でやった、ばねの運動の話がそのまま使えるのです。 電子 -q つまり、原子が双極子になるとそれは一定の固有角周波数 ω0 で振動します。 しかし、原子内部の摩擦により、放っておかれた原子双極子は次第に振動が減衰し、ついに z = 0 つまり中性原子に戻ってしまいます。 そこに角振動数ωの周期的な外力が働くと、原子双極子は強制振動状態となり、固有振動数ではなく、外力の振動数で揺れるようになります。

  22. 古典的電気双極子の吸収断面積 外力として、入射電磁波の振動電場 E = Eo eiωt を考えます。   z=原子から電子への位置ベクトル   m=電子質量(本当は電子と原子核の換算質量) ω=電磁波の角振動数   ω0=原子の角固有振動数     となるので、運動方程式は 光 原子 ここで、 とおくと、 ―e z これは強制振動の式です。 その解はもう分かっていて、 +e

  23. ここで、振動電場 E = Eo eiωtが一方向に進む平面電磁波によるものと考えます。この平面波の運ぶエネルギーフラックス Fは、cgs単位系で、 です。ここで、C-1の最後をもう一度見て下さい。電気双極子のエネルギー吸収を吸収断面積 σ で表わすと <W>=Fσ です。 これが古典電磁気学が与える電気双極子の吸収断面積です。

  24. 下の図はσの表式中の だけを示したものです。 ν0=100の場合につき、 γ=1とγ=3の例を図示しました。吸収断面積σが ν=ν0 のところでピークを持つことがわかります。これを共鳴吸収と呼びます。

  25. D.2.で摩擦によるバネの減衰が E=Eo・exp (-γt ) で表わされ、これはγが大きいほどエネルギーの散逸が大きいことを示していると学びました。 ところが、前ページの図を見ると、γが大きいと共鳴吸収のピークの高さは低下します。摩擦係数γが大きいとエネルギー散逸が大きくなり、電磁波のエネルギーを大量に吸収するような気がしますが、逆にピークが下がるのはなぜでしょう? 

  26. D.4.振動子強度 (Oscillator Strength、  f-value) 吸収強度の問題 D.3.で学んだ電気双極子の吸収断面積モデルでは、双極子は固有振動にごく近い振動数の光を吸収します。それがスペクトルの吸収線を産み出すのです。前節ではその振動数毎に吸収断面積がどう変わるかを学びました。 しかし、 (1) 一般に光は連続した振動数スペクトルを持っています。 (2) 原子の運動のドップラーシフトが吸収断面積をシフトさせます。 (3) 測定装置の波長分解能が吸収線の形をなまらせます。 などから、吸収線の測定は、 波長毎にどのくらい光が弱くなったというより、 連続光から吸収線全体ではどのくらいの光が奪われたかを調べる、 ことに主眼がおかれるケースが多いのです。

  27. σν Fν Fν X = ν ν ν (ドップラーシフト) X = = (分解能) X

  28. 線吸収の強度 Fν = Fν ν ν 上の図のように、ドップラーシフトや装置の分解能では、吸収線の面積は変わらないという性質があります。 D.1.でやったように、                  dF(ν) =N・σ(ν) ・F(ν) ですが、吸収線の巾程度では通常F(ν)の変化は小さいので、F(ν)=Fo=一定とします。すると、図の赤線部分は、 ΔF=∫N・F(ν)・σ(ν)dν=N・Fo・ ∫σ(ν)dν

  29. 下には をν0=100、 γ=1とγ=3の例を 図示しました。 青と赤の線の下の面積が等しいことが判りますか?

  30. 吸収断面積の積分からはγが消えます それを式で確認しましょう。 を代入して、 周波数積分では、 波長積分では、 数値を入れて実際の値を計算すると、

  31. σ(ν) a π ( )d = ( ) ∫σ ν ν σ ν f[2mc /(h /4 ) ] 3 2 γ π a 2 fc/ c 2 3 π α π λ =( q /mc)f 2 π o-2 /4 o- /4 o o+ /4 o+2 /4 ν γ π ν γ π ν ν γ π ν γ π 2 /4 γ π νo-2γ/4π νo-γ/4π νo νo+γ/4π νo+2γ/4π 吸収断面積の積分からγが消える理由は前に述べましたが、もう一度 下に図示します。 吸収断面積σ(ν) (e2/mc)(4π/γ) 積分値= (πq2/mc) はγに依らない。 γ/2π 結局、吸収が弱い近似で計算すると、 で、どの吸収線も強度は一定となります。しかし、実際には吸収線毎にその強度は様々な値を取る。古典的電気双極子モデルではこの違いを説明できませんでした。

  32. f=oscillator strength の導入 古典物理学では、電気双極子の吸収を計算すると吸収の強さが皆同じになり、原子の吸収線強度が様々であることと矛盾します。この問題は多くの科学者を悩ませましたが解決には量子力学の誕生を待たなければなりませんでした。 量子力学によって電気双極子の吸収を計算すると、古典電磁気学が与えた吸収断面積に f という係数をかければよいことが分かりました。 したがって、量子力学的双極子による吸収断面積は となります。 原子の吸収線毎に f は異なります。線の強さは f の大きさで決まります。そこで f=oscillator strength またはf-値( f-value)  と呼ばれます。 また、等値巾Wは

  33. 吸収断面積σpピーク値の概算 吸収線ピークの吸収断面積σpを概算する場合は、は線幅Dを使って、 σp=(πe2/mc) (λ2/c) f/D=2.654・10-2(cm2sec-1)f・(λ2/Dc) Hα: λ=0.65μ=0.6563・10-4cm D=0.0001μ=10-8cm c=2.998・1010cm/secf=0.6407   を代入すると、 Hβ: λ=0.4861μ=0.4861・10-4cm D=0.0001μ=10-8cm c=2.998・1010cm/secf=0.1193   を代入すると、

  34.  振動子強度の例 n=2 l=1 S=1/2 L=1 g=4 n=2 l=0 S=1/2 L=0 2P3/2 g=2 g=2 2P1/2 2S1/2 n=1 l=0 S=1/2 L=0 2S1/2 g=2 例1:Lα線 g (1s2S1/2) f(1s2S1/22p2P1/2)=0.2774, f(1s2S1/22p2P1/2) =0.1387 g (1s2S1/2) f(1s2S1/22p2P3/2)=0.5547, f(1s2S1/22p2P3/2) =0.2774 g (n=1) f(n=1n=2)=0.2774+0.5547=0.8321, f(n=1n=2) =0.4161 selection rules Δl=±1 ΔS=0、ΔL=0、±1、  ΔJ=0、±1   (J=0J=0、 L=0L=0を除く)

  35. g=6 g=4 g=4 3d2D5/2 g=2 3d2D3/2 3p2P3/2 g=2 3p2P1/2 3s2S1/2 2p2P3/2 2p2P1/2 2s2S1/2 g=4 g=2 g=2 例2:Hα レベル間遷移(ライン)のf-値 ターム間遷移(マルチプレット)のf-値 transition gLfLU gL fLU 2s2S1/23p2P1/2 0.5796 2 0.2898 2s2S1/23p2P3/2 1.1592 2 0.5796 2p2P1/23s2S1/2 0.05434 2 0.02717 2p2P3/23s2S1/2 0.10468 4 0.02717 2p2P1/23d2D3/2 2.782 2 1.391 2p2P3/23d2D3/2 0.5564 4 0.1392 2p2P3/23d2D5/2 5.008 4 1.252 transition gLfLU gL fLU 2s3p 0.8694 2 0.4347 2p3s 0.08151 6 0.01358 2p3d 4.6732 6 0.6955 Hα線のf-値 23 5.1241 8 0.6405

  36. D.5.等値幅 (Equivalent Width) 吸収線の強さをどう決めるか? の原子が数密度n、厚みLで存在するとします。 そこを通過した光 Fc は、nLσ<<1の時、吸収を受けて F(λ)=Fc(λ)[1-NnLσ(λ)] へと減衰します。 Fc(λ) F(λ) L λ λ 天体スペクトルの観測ではF(λ)から吸収線の強度を測定する必要があります。どんな量を使えばよいでしょうか?

  37. 吸収の深さを考える際に注意しなくてはならないのは、フラックスFの大きさは関係ないということです。 吸収の深さを考える際に注意しなくてはならないのは、フラックスFの大きさは関係ないということです。  Fλ A 上図で赤線のようなスペクトルに吸収線AとBがあった時に、へこみの大きさはAの方が大きく見えます。 図1 B λ Fλ 吸収線の強さは下の図、つまり、もとのフラックスFcを1として、そこから何割落ち込んだかを問題にするのです。 下の図では吸収線BがAよりずっと強いことは明らかです。 A B FC 図2 λ

  38. では、吸収線の強さを何割落ちたかで決めましょうか。では、吸収線の強さを何割落ちたかで決めましょうか。 吸収線の中心での波長をλ0とします。 吸収の深さを D0 =[Fc(λ0 )-F(λ0 )]/Fc(λ0 ) で与えるのはもっともな定義でしょう。 この D0 を吸収線強度として採用してはどうでしょう? このアイデアはかなりいいのですが難があります。 それは今問題にしているのは観測スペクトルであるということです。観測には分光器の性能に伴う波長分解能が付きまといます。 高分解分光器で 観測したスペクトル 低分解分光器で 観測したスペクトル オリジナルのスペクトル

  39. 等値巾W (Equivalent Width) そうすると、同じ星のスペクトルでも天文台毎に吸収線強度が異なってしまいます。これでは研究になりません。 そこで、吸収線の深さでなく、吸収線の面積を使うというアイデアが出てきました。その場合スペクトルは前にやった図2のようにFcで割って、吸収のない所が1になるように規格化しておきます。下の右図がそれです。 下図で吸収線を同じ面積で高さ1の長方形に直すと、その巾が吸収線面積を表わす事になります。それで、その巾を吸収線の等値幅 W と呼ぶのです。 Fc 1 F(λ) Wλ Fλ λ 0 λ λ

  40. 分光器によって吸収線の形が変わっても、この等値幅Wの値は同じという性質があります。なぜそうなるかの説明を以下に載せますが、飛ばしても構いません。要するに、分解能が悪くなると、吸収線の弱い光が周りに散らばり、逆に周りの強い光が吸収線部分に入ってきて、浅く広くなるが面積は同じということです。分光器によって吸収線の形が変わっても、この等値幅Wの値は同じという性質があります。なぜそうなるかの説明を以下に載せますが、飛ばしても構いません。要するに、分解能が悪くなると、吸収線の弱い光が周りに散らばり、逆に周りの強い光が吸収線部分に入ってきて、浅く広くなるが面積は同じということです。 星からのフラックス=G(λ)、観測で求められるスペクトル=F(λ)とします。 F(λ)には分光器の分解能 φ(λ) が次のような形で入ってきます。  F(λ)=∫G(t)φ(λ-t)dt        ここで、 ∫φ(t)dt=1 例 G(λ)=0 (-0.01<λ<0.01)   φ(λ)=1 (-0.5<λ<0.5) =5  その他          =0  その他       5 5 F(λ) G(λ) 4.9 φ(λ) × = -0.01 0.01 -0.5 0.5 -0.51 -0.49 0.49 0.51 W=0.02 等しい W=0.0002+0.0196+0.0002=0.02

  41. 一般に  F(λ)=∫G(t)φ(λ-t)dt        ここで、 ∫φ(t)dt=1 の時には、∫λF(λ)dλ= ∫λ[∫tG(t)φ(λ-t)dt]dλ= ∫t ∫λG(t)φ(λ-t)dλdt= ∫t G(t) dt したがって、 このように、等値幅Wは分光器の分解能によりません。 このため、特に分解能の低い観測ではWがよく用いられます。

  42. D.6. Voigt Profile 速度Vで動いている原子に、静止系で振動数νの光が当たる状況を考えましょう。原子はドップラーシフトした光の振動を νD と見ます。 v=V ν (νD-ν)/ν= V/cドップラーシフト νD=ν+(V/c)ν=ν+D 静止している原子の吸収断面積は、 速度分布 f(V) で動く原子の 平均吸収断面積σT(ν)は ? 1.速度Vの原子の吸収断面積 σV(ν)=σ(νD)   ここで、Vは光と同じ方向の速度成分であることに注意して下さい。

  43. 2.速度分布 f(V)、∫f(V)dV=1で規格化、 の原子の平均吸収断面積は σT(ν)=∫σV(ν)f(V)dV=∫σ(νD)f(V)dV   で与えられます。 D=(V/c)ν なので、 ただし、

  44. 3.σT(ν)=∫σ(νD)f(V)dV をDの積分で表示すると、3.σT(ν)=∫σ(νD)f(V)dV をDの積分で表示すると、 4.νDで規格化します。

  45. νD=熱運動に    よるドップラー巾 = Voigt function ∫V(u,a)du=1 =中心周波数との差 =吸収線自然巾 =ドップラーシフト 熱運動をする気体原子の平均吸収断面積σT(ν) σT(ν) ドップラー     核 ローレンツ    ウィング a=0 a=0.03 ν

  46. Voigt関数の性質 (1)a<<1の場合 (自然巾<<熱運動の巾、大抵の吸収線では成立) H or G (2) 1/(aπ) 2a x ‐u

  47. (3) H(u=0) << G(x=-u) 、 大体 u<≒1、 の領域では 原子の熱運動によるドップラーシフトが支配的でガウス型のプロファイルとなります。吸収線の中央部分なので、ドップラーコアとも呼ばれます。 (4) H(u=0) >> G(x=-u) 、大体 1<<u、の領域では 吸収線中心から離れるとドップラーシフトの影響が弱くなり、静止原子のローレンツ型プロファイルが再び出現します。

  48. 問(1)解答 (1) に (2) (3)

  49. 問(1)解答 (4) のグラフは以下のようとなります。 上の式は始めに見た単振動 z = 5 cos(3t) の振幅が5exp(-0.5t) で減衰していくことを表わしています。これは摩擦力のために振動 のエネルギーが失われて振幅が小さくなっていくからです。

More Related