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

はじめてコードを書いた初心者です

ファブリぺロー干渉計の入力に対する出力を書きたいと思い書きました。

ですが、あまり上手くいきません。

間違えている場合や別の書き方があるならぜひ教えてください。

分かりにくい質問かもしれませんがよろしくお願いします。

式の内容はコメントに追記します

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

[Akira Agata]

図からの想像ですが、光ファイバの端面に反射膜と半透明膜を付加したファブリペロー共振器の「入力に対する出力」というのは、反射特性（反射率 vs 波長、又は反射光強度 vs 波長）のことでしょうか？　より具体的には、①反射膜の反射率、半透明膜の②透過率と③反射率、④反射膜と半透明膜間の距離、の４つの値から理論式による反射特性をプロットしたい、ということでしょうか？

[yu]

コメントありがとうございます。

現在、先端の反射膜と半透過ミラー部分の長さ変化で起こる反射特性（長さ-反射光強度）を考えています。

もしかしたら、自分で調べた式なので間違っているかもしれませんが、

現在T＝((1-A-R).^2/((1-R).^2+4*R*sin.^2*(2*pi*n*Length*f/c))) という反射率の式に入力信号をかけることで反射光強度を求めようとしています。

A=損失　R=反射膜の反射率　ｎ＝屈折率　Length＝反射膜と半透明膜間の距離（変数）

ｃ＝高速　f=周波数

この式から反射特性をプロットしたいと考えていました。

[Akira Agata]

ファブリペロー干渉計の反射特性は、半透明膜の透過率と反射率にも依存するはずですが、この式にはそれらの項が含まれていないようです。念のため、いまいちど理論式をご確認ください。

[yu]

コメント及びご助言ありがとうございます。

おっしゃられた通り、今までの式は恐らく反射膜と半透明膜が同一のものだった場合の式だったようです

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ω/ｃ

Ei＝レーザー光振幅　L=反射膜と半透明膜間の距離（変数）Rx,Tx.=半透明膜の反射率及び透過率　Ry,Ty=反射膜の反射率及び透過率　ｃ＝光速　ω＝レーザーの角周波数

恐らくこちらの式が正しいようです。

この式からL（反射膜と半透明膜間の距離）で変化する反射特性をプロットしたいと考えます。

Answer 1

[Akira Agata]

初期パラメータを以下のように想定して計算したところ、波長によっては透過率（＝出力光電力／入力光電力）が 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)

[yu]

自分の知識不足及び説明不足で手間取らせてしまい大変申し訳ございません。

│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]

ご確認頂き、ありがとうございます。

ここまで来れば、あとはご指摘のとおり、波長範囲を固定して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]

本当に大変お世話になりました

forループを使わない計算なども大変参考になりました。

今後もコードの勉強に励んでいきたいと思います。

Answer 2

[Hiroumi Mita]

構造体を使っているようですが、プログラムの基本は

部分から全体、簡単なものから複雑なものへ変化させることなので、基本から進めましょう。

#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．入出力を構造体にします。

このように、確実に階段を昇っていくことがプログラムのコツです。

[yu]

ご回答ありがとうございます。

まだおぼつかないですが、言われた通り1つずつこなしていきます。