MATLAB Answers

ファブリペロー干渉計について

3 views (last 30 days)
yu
yu on 23 Dec 2019
Edited: yu on 27 Dec 2019
はじめてコードを書いた初心者です
ファブリぺロー干渉計の入力に対する出力を書きたいと思い書きました。
ですが、あまり上手くいきません。
間違えている場合や別の書き方があるならぜひ教えてください。
分かりにくい質問かもしれませんがよろしくお願いします。
式の内容はコメントに追記します
A=constant
R=constant
n=constant
c=constant
%
% This program simulates an Fabry-Perot
% Create an output structure similar to the input
OutputPort1 = InputPort1;
% Length variation
Length = Parameter0;
% calculate the optical signal
if(InputPort1.TypeSignal =='Optical' )
% verify how many sampled signals are in the structure
[ls, cs] = size(InputPort1.Sampled);
if( ls > 0 )
% caculate the at each signal
for counter1=1:cs
OutputPort1.Sampled(1, counter1).Signal = InputPort1.Sampled(1, counter1).Signal * ((1-A-R).^2/(1-R).^2+4*R*sin.^2*(2*pi*n*Length*f/c));
end
end
end

  4 Comments

Show 1 older comment
yu
yu on 24 Dec 2019
コメントありがとうございます。
現在、先端の反射膜と半透過ミラー部分の長さ変化で起こる反射特性(長さ-反射光強度)を考えています。
もしかしたら、自分で調べた式なので間違っているかもしれませんが、
現在T=((1-A-R).^2/((1-R).^2+4*R*sin.^2*(2*pi*n*Length*f/c))) という反射率の式に入力信号をかけることで反射光強度を求めようとしています。
A=損失 R=反射膜の反射率 n=屈折率 Length=反射膜と半透明膜間の距離(変数)
c=高速 f=周波数
この式から反射特性をプロットしたいと考えていました。
Akira Agata
Akira Agata on 26 Dec 2019
ファブリペロー干渉計の反射特性は、半透明膜の透過率と反射率にも依存するはずですが、この式にはそれらの項が含まれていないようです。念のため、いまいちど理論式をご確認ください。
yu
yu on 26 Dec 2019
コメント及びご助言ありがとうございます。
おっしゃられた通り、今までの式は恐らく反射膜と半透明膜が同一のものだった場合の式だったようです
P=(((Tx.^2+Rx.^2)*Ry-Rx).^2+4*Rx*Ry(Tx.^2+Rx.^2)*sin.^2*(φ/2)/(1-Rx*Ry) .^2+4*Rx*Ry*sin.^2(φ/2))*│Ei│.^2
φ=2Lω/c
Ei=レーザー光振幅 L=反射膜と半透明膜間の距離(変数)Rx,Tx.=半透明膜の反射率及び透過率 Ry,Ty=反射膜の反射率及び透過率 c=光速 ω=レーザーの角周波数
恐らくこちらの式が正しいようです。
この式からL(反射膜と半透明膜間の距離)で変化する反射特性をプロットしたいと考えます。

Sign in to comment.

Accepted Answer

Akira Agata
Akira Agata on 26 Dec 2019
Edited: Akira Agata on 27 Dec 2019
初期パラメータを以下のように想定して計算したところ、波長によっては透過率(=出力光電力/入力光電力)が 1 を超えるという結果になっています。記載頂いた式のどこかに誤りがあると思われますので、再度ご確認頂けないでしょうか。
% 設定パラメータ
Rx = 0.45; % 半透明膜の反射率
Tx = 0.45; % 半透明膜の透過率
Ry = 0.95; % 反射膜の反射率
L = 100.0e-6; % 反射膜と半透明膜間の距離 [m]
c = 3.0e8; % 光速 [1/m]
lambda =... % 波長範囲 [m]
linspace(1.5e-6,1.6e-6,1000);
% 設定パラメータからの換算値
f = c./lambda; % 周波数 [Hz]
omega = 2*pi*f; % 各周波数 [rad/s]
phi = 2*L*omega/c;
% ファブリペロー干渉計の透過率(Tfp)を計算 (ただしTfp = P_out/P_inと定義)
Tfp =...
((Tx^2+Rx^2)*Ry-Rx)^2 +...
(4*Rx*Ry*(Tx^2+Rx^2)*(sin(phi/2)).^2)./((1-Rx*Ry)^2) +...
4*Rx*Ry*(sin(phi/2)).^2;
% 計算結果をプロット
figure
plot(lambda*1e6,Tfp)
xlabel('Wavelength [\mum]','FontSize',12)
ylabel('Transmittance','FontSize',12)

  3 Comments

yu
yu on 27 Dec 2019
自分の知識不足及び説明不足で手間取らせてしまい大変申し訳ございません。
│Ei│.^2の前の(1-Rx*Ry) .^2+4*Rx*Ry*sin.^2(φ/2))が全て分子でした。
そして前回Pの記述がなくとEiに誤りがありました。
P=反射光強度であり、Ei=レーザー光の電場振幅でした。
% 設定パラメータ
Rx = 0.45; % 半透明膜の反射率
Tx = 0.45; % 半透明膜の透過率
Ry = 0.95; % 反射膜の反射率
L = 100.0e-6; % 反射膜と半透明膜間の距離 [m]
c = 3.0e8; % 光速 [1/m]
lambda =... % 波長範囲 [m]
linspace(1.5e-6,1.6e-6,1000);
% 設定パラメータからの換算値
f = c./lambda; % 周波数 [Hz]
omega = 2*pi*f; % 各周波数 [rad/s]
phi = 2*L*omega/c;
% ファブリペロー干渉計の透過率(Tfp)を計算 (ただしTfp = P_out/P_inと定義)
Tfp =...
((Tx^2+Rx^2)*Ry-Rx)^2 +...
(4*Rx*Ry*(Tx^2+Rx^2)*(sin(phi/2)).^2)./(((1-Rx*Ry)^2) +...
4*Rx*Ry*(sin(phi/2)).^2);
% 計算結果をプロット
figure
plot(lambda*1e6,Tfp)
xlabel('Wavelength [\mum]','FontSize',12)
ylabel('Transmittance','FontSize',12)
と直してみたところ1を超えない範囲になりました。
横軸をLの変化で見る場合は波長範囲を固定し、Lを変化させるようにすればよろしいのでしょうか
Akira Agata
Akira Agata on 27 Dec 2019
ご確認頂き、ありがとうございます。
ここまで来れば、あとはご指摘のとおり、波長範囲を固定してLを変化させれば問題ありません。
forループを使って計算しても良いのですが、せっかくですので行列演算に強いMATLABの特徴を活用してみてください。たとえば以下のようにすると、forループを使わずに計算することができます(lambdaとLを、それぞれ行ベクトルと列ベクトルにするところがポイントです)。
% 設定パラメータ
Rx = 0.45; % 半透明膜の反射率
Tx = 0.45; % 半透明膜の透過率
Ry = 0.95; % 反射膜の反射率
L = linspace(100.0e-6,105e-6)'; % 反射膜と半透明膜間の距離 [m]
c = 3.0e8; % 光速 [1/m]
lambda =... % 波長範囲 [m]
linspace(1.5e-6,1.55e-6,1000);
% 設定パラメータからの換算値
f = c./lambda; % 周波数 [Hz]
omega = 2*pi*f; % 各周波数 [rad/s]
phi = 2*L*omega./c;
% ファブリペロー干渉計の透過率(Tfp)を計算 (ただしR = P_out/P_inと定義)
Tfp =...
((Tx^2+Rx^2)*Ry-Rx)^2 +...
(4*Rx*Ry*(Tx^2+Rx^2)*(sin(phi/2)).^2)./(((1-Rx*Ry)^2) +...
4*Rx*Ry*(sin(phi/2)).^2);
% 計算結果をプロット
figure
surf(lambda*1e6,L*1e6,Tfp,...
'EdgeColor','none')
xlabel('Wavelength [\mum]',...
'FontSize',12,...
'Rotation',-20)
ylabel('Cavity length [\mum]',...
'FontSize',12,...
'Rotation',45)
zlabel('Transmittance','FontSize',12)
view(30,60)
yu
yu on 27 Dec 2019
本当に大変お世話になりました
forループを使わない計算なども大変参考になりました。
今後もコードの勉強に励んでいきたいと思います。

Sign in to comment.

More Answers (1)

Hiroumi Mita
Hiroumi Mita on 23 Dec 2019
構造体を使っているようですが、プログラムの基本は
部分から全体、簡単なものから複雑なものへ変化させることなので、基本から進めましょう。
#1. まず、この式の正誤を確認しましょう。
OutputPort1.Sampled(1, counter1).Signal = InputPort1.Sampled(1, counter1).Signal * ((1-A-R).^2/(1-R).^2+4*R*sin.^2*(2*pi*n*Length*f/c));
#2. あっているようであれば
適当な数値を各パラメータにあてはめて期待値がでるか確認しましょう。
3, 関数にします。
4.入出力を構造体にします。
このように、確実に階段を昇っていくことがプログラムのコツです。

  1 Comment

yu
yu on 24 Dec 2019
ご回答ありがとうございます。
まだおぼつかないですが、言われた通り1つずつこなしていきます。

Sign in to comment.

Sign in to answer this question.

Products


Release

R2019a