MENU

溶けかけてるうさぎ HP GALLERY BLOG TOP RECENT ARTICLES POPULAR ARTICLES ABOUT THIS BLOG

CATEGORY

大学 (140) 仕事 (17) 航空宇宙 (104) 写真 (78) 旅行 (32) 飯・酒 (17) コンピュータ (119) その他 (44)

TAG

ARCHIVE

RECENT

【写真】撮影写真を Map 上に表示できるようにした 【カメラ】X100 シリーズが好きすぎる(主にリーフシャッタ) 【カメラ】X100V から X100VI に買い替えました 【自宅サーバー】Google Domains から Cloudflare にドメインを移管 【カメラ】FUJIFILM XF レンズのサイズ比較ができるようにしてみた

【波動光学】波動光学入門1 - システム応答の導出

事象発生日:2017-12-23

記事公開日:2017-12-23

アクセス数:57043

この記事は,東京大学航空宇宙工学科/専攻 Advent Calendar 2017向けの記事です.

 

今年の4~6月あたりに得た波動光学の基礎をまとめておく.

波動光学,和書があまりなく,毎回毎回洋書を参照するのもだるいので.

 

トップ画像の出典は

1.はじめに

モチベーション

」などの記事で,レンズの基本的な物理量についてみてきた.

F値や被写界深度は幾何光学,つまり光が進む経路(反射や屈折など)を考えるだけで十分だった.

 

しかしながら,MTFなどは幾何光学では説明不可能である.

のようにレンズメーカーのホームページにはMTF特性図が載っている.)

CanonレンズのMTF特性

そこで,ここでは光を波動として考える波動光学の基礎について,備忘録的にまとめたい.

(備忘録的,というのは,波動光学の考え方を卒業研究で多用し修士研究でも用いる予定なのだが,現在は卒業設計に時間をとられ,せっかく学んだ波動光学の基礎を忘れそう,という事情による.)

波動光学

波動光学とは,光を波動として考える光学であり,干渉・回折・偏光などが扱える.

 

中学,高校で学ぶレンズの公式などは幾何光学であり,ヤングの実験などが波動光学である.

さらに発展として,量子力学も考慮した量子光学という分野も存在する.

本記事の流れ

はじめに「」で波動光学の基礎式となる回折現象を定式化する.

次に「」で薄肉レンズでの像をそのインパルス応答を定式化することによって理解する.

さらに「」で一般の結像系について,そのシステム応答を考える.

そして続編の記事「」では,波面収差の取り扱いと実際の数値計算を行っている.

 

最初は式変形ばかりでわかりにくいが,続編で可視化したものを見れば雰囲気はわかるだろう.

2.波動光学の基本的な回折の定式化

幾何光学では,点光源をレンズを介して再び完全な点像へと集光できた.

しかし,実際にはレンズを通ったときの回折により一点へと集光させることは不可能である.

 

そこで,まずはじめに回折を簡単に記述するための波動光学の基礎近似式を導出する.

 

ここでは簡単のために完全な単色光を仮定する.

 

位置\(\boldsymbol{r}\)での正弦波複素振幅\(U(\boldsymbol{r}, t)\)は,波数ベクトル\(\boldsymbol{k}\)を用いて,

\begin{align*} U(\boldsymbol{r}, t) = A \exp \left[ i \left( \boldsymbol{k}\cdot \boldsymbol{r} - \omega t \right) \right] \tag{2-1} \end{align*}

となるのは流石にいいよね?(卒業設計後にこれを読み返すであろう自分を煽っておく.)

(まあ,ここではもっぱら時間平均しか考えないので,時間項は落ちるが.)

ホイヘンス - フレネル原理による回折

回折を考えるための座標系

上の座標系を用いて,ホイヘンス - フレネル原理に基づく回折式を導入する.

 

ホイヘンス - フレネル原理とは,高校の教科書ではホイヘンスの原理と呼ばれているあれで,

“前進波の波面の各点が二次波とよばれる新しい波の波源となり,全体としての前進波は全ての二次波を重ね合わせたものとなる.”

というもの.

ホイヘンス - フレネル原理

 

このとき,点\(P_0\)での複素振幅\(U(P_0)\)は,点\(P_0\),\(P_1\)間距離\(r_{01}\)を用いて

\begin{align*} U(P_0) = \frac{1}{i\lambda} \iint\limits_\Sigma U(P_1) \frac{e^{i k r_{01}}}{r_{01}} \cos\theta ds \tag{2-2} \end{align*}

となる.さらに

\begin{align*} \cos\theta = \frac{z}{r_{01}} \tag{2-3} \end{align*}

より,ホイヘンス - フレネル原理に基づく回折式は

\begin{align*} U(x,y) = \frac{z}{i\lambda} \iint\limits_\Sigma U(\xi,\eta) \frac{e^{i k r_{01}}}{r_{01}^2} d\xi d\eta \tag{2-4} \end{align*}

となる.

 

 

...と,まあいきなり(2-2)式を持ち出して変形したわけだが,一般的な光学入門書だと,(2-4)式で\(z \approx r_{01}\)より

\begin{align*} U(x,y) = \iint\limits_\Sigma U(\xi,\eta) \frac{e^{i k r}}{i \lambda r} d\xi d\eta \tag{2-5} \end{align*}

から議論がスタートする.スケール因子\(1/i\lambda\)の説明もなしで.

波動光学を初めて勉強したときはこの(2-5)式をいきなり出されて??となっていたが,定性的には,

\( U(\xi,\eta) \exp \left[ i k r \right]/r \)がホイヘンス - フレネル原理の二次波の球面波であり,それを回折開口で積分している,と理解できる.

(球面波とは,点波源より伝播する波面が球となる波.点波源(波面の中心)からの距離\(r\)を用いるとその振幅は\(\exp \left[ i k r \right] / r\)に比例.)

 

(2-2)式の導出はとても煩雑であるので,あたりを参照されたいが,概要は,波動方程式を開口によって決まる境界条件(キルヒホッフ境界条件)の下で解けば求まる.

フレネル回折

(2-4)式を近似していく.

\begin{align*} r_{01} &= \sqrt{z^2 + (x-\xi)^2 + (y-\eta)^2} \\ &\approx z \left[ 1 + \frac{1}{2}\left(\frac{x-\xi}{z}\right)^2 + \frac{1}{2}\left(\frac{y-\eta}{z}\right)^2 \right] \tag{2-6} \end{align*}

となるが,(2-4)式において,分母の\(r_{01}\)はこの二次の項を落とすことによる影響は小さいが,\(\exp\)内の\(r_{01}\)でこれを落とすのは影響が無視できない.

したがって,(2-4)式は

\begin{align*} U(x,y) = \frac{e^{ikz}}{i\lambda z} \iint\limits_{-\infty}^{\infty} U(\xi,\eta) e^{i\frac{k}{2z} \left[ (x-\xi)^2 + (y-\eta)^2 \right] } d\xi d\eta \tag{2-7} \end{align*}

となり,畳み込みの形で表される.これがフレネル回折積分と呼ばれるものである.

ここで積分範囲が無限だが,開口が有限範囲であるという条件は\(U(\xi,\eta)\)に含まれている.

 

さらに\(\xi, \eta\)がかかわらない項を外に出して,

\begin{align*} U(x,y) &= \frac{e^{ikz}}{i\lambda z} e^{ i\frac{k}{2z} \left(x^2 + y^2\right) } \\ & \qquad \iint\limits_{-\infty}^{\infty} \left\{ U(\xi,\eta) e^{ i\frac{k}{2z} \left( \xi^2 + \eta^2 \right) } \right\} e^{-j\frac{2\pi}{\lambda z}(x\xi+y\eta)} d\xi d\eta \\ &= \frac{e^{ikz}}{i\lambda z} e^{ i\frac{k}{2z} \left(x^2 + y^2\right) } \mathcal{F}\left\{ U(\xi,\eta) e^{ i\frac{k}{2z} \left( \xi^2 + \eta^2 \right) } \right\}_{f_{\xi}=\frac{x}{\lambda z},f_{\eta}=\frac{y}{\lambda z}} \tag{2-8} \end{align*}

とフーリエ変換\(\mathcal{F}\{\}\)で表すことができる.

フラウンホーファー回折

さらに厳しい近似(フラウンホーファー近似)

\begin{align*} z \gg \frac{k\left(\xi^2+\eta^2\right)_\max}{2} \end{align*}

を(2-8)式にほどこすと,

\begin{align*} U(x,y) &= \frac{e^{ikz} e^{ i\frac{k}{2z} \left(x^2 + y^2\right) }}{i\lambda z} \iint\limits_{-\infty}^{\infty} \left\{ U(\xi,\eta) \right\} e^{-j\frac{2\pi}{\lambda z}(x\xi+y\eta)} d\xi d\eta \\ &= \frac{e^{ikz} e^{ i\frac{k}{2z} \left(x^2 + y^2\right) }}{i\lambda z} \mathcal{F}\left\{ U(\xi,\eta) \right\}_{f_{\xi}=\frac{x}{\lambda z},f_{\eta}=\frac{y}{\lambda z}} \tag{2-10} \end{align*}

となり,これをフラウンホーファー回折積分といい,\(z\)が大きい,つまり無限遠で観測させる場合に有効である.

3.像形成

今までは,有限の開口を光が通過することによって生じる回折を定量的に記述した.

次は下のような座標系で像形成を考えよう.

なお,以下では定数係数は省略する.

像形成のための配置と座標

一般論として,線形時不変なシステム応答を考えた場合,出力は単位インパルス応答の畳み込みで表示される.

これを光学系へあてはめると,単位インパルス応答は単位振幅の点光源像が相当し,それを\(h\)とする.

するとある物体\(U_o\)が作る像\(U_i\)は

\begin{align*} U_i(u,v) &= h(u,v) \ast U_o(u,v) \\ &= \iint\limits_{-\infty}^{\infty} h(u-\xi, v-\eta) U_o(\xi,\eta) d\xi d\eta \tag{3-1} \end{align*}

となる.

後の式変形からわかるが,実際に求めると

\begin{align*} U_i(u,v) &= \iint\limits_{-\infty}^{\infty} h(u-M\xi, v-M\eta) U_o(\xi,\eta) d\xi d\eta \tag{3-2} \end{align*}

とスケール因子がかかってしまうので,ここでは\(h(u-\xi, v-\eta)\)を一般的に\(h(u,v;\xi,\eta)\)と書くことにする.

 

光学系に限った話ではなく,線形時不変システムに一般的な話だが,結像系の特性はインパルス応答\(h\)を特定することによって完全に記述できるので,この\(h\)を求める.

 

インパルス応答は点光源像となるため,物体を点\((\xi,\eta)\)においたデルタ関数とする.

この場合,レンズに入射するのは点\((\xi,\eta)\)を中心とする発散球面波であり,球面波は波源からの距離\(r\)を用いて

\begin{align*} \frac{e^{ i k r }}{r} \tag{3-3} \end{align*}

となるが,フレネル近似と同様に近軸近似をほどこし

\begin{align*} U_l(x,y) = \frac{1}{i\lambda z_1} e^{i\frac{k}{2z_1} \left[ (x-\xi)^2 + (y-\eta)^2 \right] } \tag{3-4} \end{align*}

となる((2-7)式の開口を無限小点にした場合と同値).

また,焦点距離\(f\)のレンズを通過した後の光波場分布は,導出は省略するが,

\begin{align*} U_l'(x,y) = U_l(x,y) P(x,y) e^{-i\frac{k}{2f} \left( x^2 + y^2 \right) } \tag{3-5} \end{align*}

となる(凸レンズは発散球面波を収束球面波へ変換すると考えれば自明).

なお,\(P\)とは瞳関数と呼ばれるものであり,レンズ開口内で\(1\)を,開口外で\(0\)をとる関数であり,レンズの形状を規定する.

最後に,距離\(z_2\)の伝播をフレネル回折の式(2-7)を用いると,

\begin{align*} h(u,v;\xi,\eta) = \frac{1}{i\lambda z_2} \iint\limits_{-\infty}^{\infty} U_l'(x,y) e^{i\frac{k}{2z_2} \left[ (u-x)^2 + (v-y)^2 \right] } dxdy \tag{3-6} \end{align*}

が得られる.

上記3式を組み合わせると,

\begin{align*} h(u,v;\xi,\eta) &= \frac{1}{\lambda^2 z_1 z_2} e^{i\frac{k}{2z_2} \left(u^2 + v^2\right)} e^{i\frac{k}{2z_1} \left(\xi^2 + \eta^2\right)} \\ & \qquad \times \iint\limits_{-\infty}^{\infty} P(x,y) e^{i\frac{k}{2} \left(\frac{1}{z_1}+\frac{1}{z_2}-\frac{1}{f}\right) \left(x^2 + y^2\right)} \\ & \qquad \times e^{-jk \left[ \left( \frac{\xi}{z_1}+\frac{u}{z_2} \right)x + \left( \frac{\eta}{z_1}+\frac{v}{z_2} \right)y \right]} dxdy \tag{3-7} \end{align*}

となるが,とても煩雑である.

詳細な議論はなどに任せるが,結像面をレンズの公式をみたす距離に配置すれば

\begin{align*} \frac{1}{z_1}+\frac{1}{z_2}-\frac{1}{f} = 0 \tag{3-8} \end{align*}

とでき,

像面の強度分布のみに注目する場合で,雑煮関連する位相分布は重要ではない.
像の光波場分布に注目するが,像の観察は,光軸が薄肉レンズを貫く点を中心とする半径\(z_2\)の球面上で行う.

のどちらかをみたす場合,

\begin{align*} e^{i\frac{k}{2z_2} \left(u^2 + v^2\right)} \approx 1 \tag{3-9} \end{align*}

の近似が成立し,さらに

光軸が薄肉レンズを貫く点を中心とする半径\(z_1\)の球面上に物体がある.
光軸がレンズを貫く点に向かって収束する球面波で物体が照射される.
この2次位相因子の位相が,特定の像点\((u,v)\)の光波場に大きく寄与する物体部分の内側ではほとんど変化しない.

のいずれかをみたす場合,

\begin{align*} e^{i\frac{k}{2z_1} \left(\xi^2 + \eta^2\right)} \approx 1 \tag{3-10} \end{align*}

の近似が成立する.

さらにシステムの倍率

\begin{align*} M = - \frac{z_2}{z_1} \tag{3-11} \end{align*}

を定義すると,(3-7)式は

\begin{align*} & h(u,v;\xi,\eta) \\ & \qquad \approx \frac{1}{\lambda^2 z_1 z_2} \iint\limits_{-\infty}^{\infty} P(x,y) e^{-i\frac{2\pi}{\lambda z_2} \left[ (u-M\xi)x+(v-M\eta)y \right] } dxdy \tag{3-12} \end{align*}

となる.

したがって,レンズの公式がみたされる場合,インパルス応答は(スケール因子\(1/\lambda z_1\)も含め),像座標\((u=M\xi, v=M\eta)\)を中心に位置するレンズ開口のフラウンホーファー回折像と等しくなっている.

これは少し考えてみればそれほど意外な結果ではない.

つまりレンズの公式をみたすように\(z_2\)を設定したので,像平面はレンズ通過後の球面波が収束する点であり,この収束点近傍の光波場分布は,球面波の拡がりを制限するレンズ開口のまさにフラウンホーファー回折像である.

4.光学系のシステム応答解析

さらに結像系を一般化し,結像系を線形時不変システムとして描写していく.

これ以降の議論は和書ではほとんど触れられておらず,仕方がないのでなどで学んだ.

 

ここではまず光学系の一般化モデルを定義し,次に回折限界(幾何光学的には点光源(発散球面波)が点像(収束球面波)へと移されるシステム.波動光学的には回折の影響により,完全な点像とはならない.)における結像システムを定式化する.

さらにその後にその周波数特性についてみていく.

光学系の一般化モデル

再びCanonのレンズを例に取るが,下図の用に光学系は多数のレンズ,鏡群によって構成されている.

Canonレンズのレンズ構成

しかし物体上の点からは球面波が放射され,そして結像系を経て(球面波となり)最終的に像として結像すると仮定すれば,結像系を一つのブラックボックスにまとめてしまっても,そのシステム特性を記述する上では問題ない.

ここでは,下図のように物体,入射瞳,射出瞳,像を定義し,また座標系も下図のように設定する.

光学系の一般化モデル

コヒーレント結像系

空間的にコヒーレントな光学系の場合,像の複素振幅\(U_i\)は物体の複素振幅\(U_o\) と点像振幅分布関数 (Amplitude Spread Function / ASF) \(h\)の畳み込み

\begin{align*} U_i(u,v) &= \iint\limits_{-\infty}^{\infty} h(u-\xi, v-\eta) U_o(\xi,\eta) d\xi d\eta \\ &= h(u,v) \ast U_o(u,v) \tag{4-1} \end{align*}

で表される.ASFは,「」での議論により,理想像(幾何光学的に求まる無回折時の像)点付近の光振幅は射出瞳のフラウンホーファー回折像であるため,

\begin{align*} h(u,v;\xi,\eta) = \frac{A}{\lambda z_i} \iint\limits_{-\infty}^{\infty} P(x,y) e^{-i\frac{2\pi}{\lambda z_i} \left[ (u-M\xi)x+(v-M\eta)y \right] } dxdy \tag{4-2} \end{align*}

となった.

ただし\(z_i\)は射出瞳―像平面間距離,\(P\)は瞳関数(開口内側で\(1\),外側で\(0\))である.

ここで

\begin{align*} \xi = M\xi, \qquad \eta = M\eta \tag{4-3} \end{align*}

と座標を取り直すとASFは,

\begin{align*} h(u-\xi,v-\eta) = \frac{A}{\lambda z_i} \iint\limits_{-\infty}^{\infty} P(x,y) e^{-i\frac{2\pi}{\lambda z_i} \left[ (u-\xi)x+(v-\eta)y \right] } dxdy \tag{4-4} \end{align*}

すなわち

\begin{align*} h(u,v) &= \frac{A}{\lambda z_i} \iint\limits_{-\infty}^{\infty} P(x,y) e^{-i\frac{2\pi}{\lambda z_i} \left(ux+vy \right) } dxdy \\ &= \frac{A}{\lambda z_i} \mathcal{F} \left\{P(x,y)\right\}_{f_x=\frac{u}{\lambda z_i},f_{y}=\frac{v}{\lambda z_i}} \tag{4-5} \end{align*}

と瞳関数のフーリエ変換で表される.

 

以上より,一般的な回折限界の光学系の場合の像は,幾何光学的像と,射出瞳のフラウンホーファー回折像であるインパルス応答との畳み込みであるとみなせる.

つまり,コヒーレント結像系では,複素振幅に対して線形である.

インコヒーレント結像系

空間的にインコヒーレントな光学系の場合,観測されるのは像の複素振幅\(U_i\)ではなく強度\(I_i\)である.

像の強度\(I_i\)は物体の強度\(I_o\)と点拡がり関数 (Point Spread Function / PSF) \(s\)の畳み込み

\begin{align*} I_i(u,v) &= \iint\limits_{-\infty}^{\infty} s(u-\xi, v-\eta) I_o(\xi,\eta) d\xi d\eta \\ &= s(u,v) \ast I_o(u,v) \tag{4-6} \end{align*}

で表される.

\(I\)は時間平均演算子\(\langle\rangle\)を用いて

\begin{align*} I(u,v) = \langle\left| U(u,v,t) \right|^2\rangle = \langle U(u,v,t) U^{\ast}(u,v,t) \rangle \tag{4-7} \end{align*}

となることより,PSFはASFの絶対値の二乗

\begin{align*} s(u,v) = \left| h(u,v) \right|^2 = \left| \frac{A}{\lambda z_i} \mathcal{F} \left\{P(x,y)\right\}_{f_x=\frac{u}{\lambda z_i},f_{y}=\frac{v}{\lambda z_i}} \right|^2 \tag{4-8} \end{align*}

となる.

 

以上より,インコヒーレント結像系は強度に関して線形であり,そのような結像系のインパルス応答は,複素インパルス応答の大きさの二乗で表される.

コヒーレント結像系での周波数応答

(4-1)式のフーリエ変換をとると,畳み込みの公式より,

\begin{align*} \mathcal{F}\left\{U_i(u,v)\right\} &= \mathcal{F}\left\{h(u,v)\right\} \mathcal{F}\left\{U_o(u,v)\right\} \tag{4-9} \end{align*}

となり,

\begin{align*} H(f_x,f_y) = \mathcal{F}\left\{h(u,v)\right\} \tag{4-10} \end{align*}

は振幅伝達関数 (amplitude transfer function / ATF) と呼ばれ,システムの周波数応答を表す.

ATFは

\begin{align*} H(f_x,f_y) &= \mathcal{F}\left\{h(u,v)\right\} \\ &= \mathcal{F}\left\{ \frac{A}{\lambda z_i} \mathcal{F} \left\{P(x,y)\right\}_{f_x=\frac{u}{\lambda z_i},f_{y}=\frac{v}{\lambda z_i}} \right\} \\ &= A\lambda z_i P(-\lambda z_i f_x, -\lambda z_i f_y) \tag{4-11} \end{align*}

となる.

 

したがって,コヒーレント結像系ではその周波数特性は瞳関数と一致し,その通過周波数では振幅や位相を歪めることなく通過させる.

インコヒーレント結像系での周波数応答

インコヒーレント結像系でもコヒーレント結像系の場合と同様にシステムの周波数応答\(\mathcal{H}\)を(0周波数の値で正規化して)

\begin{align*} \mathcal{H}(f_x,f_y) &= \frac{\mathcal{F}\left\{s(u,v)\right\}}{\iint\limits_{-\infty}^{\infty}s(u,v)dudv} \\ &= \frac{\mathcal{F}\left\{\left|h(u,v)\right|^2\right\}}{\iint\limits_{-\infty}^{\infty}\left|h(u,v)\right|^2dudv} \tag{4-12} \end{align*}

と求めることができる.

\(\mathcal{H}\)は光伝達関数 (optical transfer function / OTF) といい,その絶対値\(|\mathcal{H}|\)を変調伝達関数 (modulation transfer function / MTF) という.

 

なお,OTFはパーセバルの定理を用いて

\begin{align*} \mathcal{H}(f_x,f_y) &= \frac{\mathcal{F}\left\{\left|h(u,v)\right|^2\right\}}{\iint\limits_{-\infty}^{\infty}\left|h(u,v)\right|^2dudv} \\ &= \frac{\iint\limits_{-\infty}^{\infty} H^{\ast}(p-f_x,q-f_y)H(p,q)dpdq}{\iint\limits_{-\infty}^{\infty} \left|H(p,q)\right|^2dpdq} \\ &= \frac{\iint\limits_{-\infty}^{\infty} P^{\ast}(x-\lambda z_i f_x,y-\lambda z_i f_y)P(x,y)dxdy}{\iint\limits_{-\infty}^{\infty} \left|P(x,y)\right|^2dxdq} \tag{4-13} \end{align*}

となり,瞳関数の自己相関関数であることがわかる.

5.つづき

」で,この記事のつづきとして,

波面収差が存在する場合
具体的な数値計算結果

を掲載する.

6.出典

Canon. EF24-70mm F4L IS USM 仕様. Retrieved December 24, 2017, from http://cweb.canon.jp/ef/lineup/standard-zoom/ef24-70-f4l/spec.html
Jason D. Schmidt. Numerical Simulation of Optical Wave Propagation with example in MATLAB®. SPIE PRESS, 3 edition, 2013.
Joseph W. Goodman. Introduction to Fourier Optics. W. H. Freeman, 3 edition, 2004.
Joseph W. Goodman. Trans. by 尾崎 義治, 朝倉 利光. フーリエ光学. 森北出版株式会社, 第3版, 2012.
James B. Breckinridge. Basic Optics for the ASTRONOMICAL SCIENCES. SPIE PRESS, 2012.

関連記事

コメント

おにぎり   [2019/12/11 14:26]

分かりやすい説明ですごくためになりました。

 

自分は4f系の顕微鏡(焦点距離の違う2枚レンズで観察対象をフーリエ変換して逆フーリエ変換する)をシミュレーションをしたい思っているのですが、まず最初にコヒーレントな場を考えようとしています。そこで質問なのですが、レンズでフーリエ変換をしたいとき(z1=z2=f)インパルス応答などを考える必要があるんでしょうか?

(3-8)式が成り立たないのでこの場合不適切なのでしょうか?

 

溶けかけてるうさぎ(管理者)   [2019/12/14 00:58]

おにぎりさん.

光学の専門家ではないので,はっきりしたことは言えませんが,レンズが何枚だろうと基本的には1つの光学系(システム)として取り扱えるはずなので,結像させるのであれば(3-8)式が成立しないことはないと思います.

そして,(線形時不変を仮定できる場合は)インパルス応答が求まっていれば,どんな入力に対しても(任意の形状の入力に対しても)その出力を求められます.

 

おにぎり   [2019/12/14 01:27]

返信ありがとうございます。

勉強になります!

コメントを投稿

名前

Email (※公開されることはありません)

コメント