Filter löschen
Filter löschen

Determining Coupling coefficent value for mode coupling in a dual core fiber solving characterstics equation

8 Ansichten (letzte 30 Tage)
The graph obtained is not as expected. Can someone have a look on my code and figure out mistakes in it?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lambda=(1:0.01:2); %unit: um
%constants and variables
pi=3.14159; %Greek letter pi
c=2.99792458e14; % speed of light (unit: um/s)
a=1.95; % radius of core (unit: um)
d=2.52*4; % distance of separation between cores (unit: um)
d1=1.1; % air hole diameter (unit:um)
P=2.52; %pitch or hole-hole spacing
ko=2*pi./lambda;
dbar=d/a;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Calculate refractive index of core material (n1) using Sellmeier's equation
G1=0.6961663; % Sellmeier's coefficients for pure silica
G2=0.4079426;
G3=0.8974794;
lambda1=.0684043; %wavelength for pure silica
lambda2=0.1162414;
lambda3=9.896161;
npower2oflambda=1+(G1*lambda.^2./(lambda.^2-power(lambda1,2)))+(G2*lambda.^2./(lambda.^2-power(lambda2,2)))+...
(G3*lambda.^2./(lambda.^2-power(lambda3,2)));
noflambda=sqrt(npower2oflambda);
n1=noflambda;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%calculate n2 in terms of n1
NA=0.20;
n2=sqrt(n1.^2-NA^2);
% calculate V-number
V=ko.*a.*(n1.^2-n2.^2 ).^(1/2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SOLVING BESSEL FUNCTION EIGEN VALUE EQUATION
syms U W
k1=(ko.*n1); %maxm value of beta
k2=(ko.*n2); %minimum value of beta
% beta=linspace(k2,k1,101);
beta=real(k2):real(k1-k2)/1e2:real(k1);
% n_eff=real(n2):real(n1-n2)/1e2:real(n1);
% beta=ko.*n_eff;
U=a.*sqrt(ko.^2.*n1.^2-beta.^2); %transverse mode paramter (phase constant)
W=a.*sqrt(beta.^2-ko.^2.*n2.^2); %transverse mode parameter (attenuation coefficient)
% h=U.*besselj(1,U)./besselj(0,U);
% m=W.*besselk(1,W)./besselk(0,W);
%solution=solve(U.*besselj(1,U)./besselj(0,U)==W.*besselk(1,W)./besselk(0,W), W);
%solution=double(solution)
figure(1)
plot(beta,h,'g','linewidth',2)
hold on
plot(beta,m,'r','linewidth',2)
legend('J_1(U)/U*J_0(U)','K_1(W)/W*K_0(W)')
xlabel('\delta\beta(um)')
grid on
hold off
betaroot=9.034; % first root from the plot 1
U1=a*sqrt(ko.^2.*n1.^2-betaroot.^2); %transverse mode parameter (phase constant)
W1=a*sqrt(betaroot.^2-ko.^2.*n2.^2);
%Uroot=imag(U1);
%Wroot=abs(W1);
K_0=besselk(0,W1.*dbar); %Zeroth order modified bessel function of 2nd kind
K_1=besselk(1,W1); %first order modified bessel function of 2nd kind
% %calculate coupling coefficient
R=sqrt(2*del)/a.*((U1.*U1)./(V.*V.*V)); % 1st part of coupling coefficient equation
Q=K_01./(K_1.*K_1); % 2nd part of coupling coefficient equation
q=abs(Q);
CC=R.*q;%Q;
figure(2)
plot(lambda,CC)
xlabel('wavelength (um)'); ylabel('coupling coeff. (/m)')
title('Coupling coeff. Solving characterstics equation');

Antworten (0)

Kategorien

Mehr zu Bessel functions finden Sie in Help Center und File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by