% 计算弯损,单一弯曲半径
function [neff,Field2Int,delta]=NBendloss(modeOrder,lam,neffRef,a,dn,nCladd,rBend,sCore)
%Calculation of effective RI without bending
neff= real( fsolve(@(neff)NLP_m_mode(neff,lam,a,dn,nCladd,modeOrder),neffRef,optimset('TolFun',1e-25,'Display','off')));
[MFD, Aeff, Field2Int]=NMFDCalculation(neff,lam,a,dn,nCladd,modeOrder);
[CO,COout]=NCoefficientField(neff,lam,a,dn,nCladd,modeOrder);
%---------------------------------------------------------
% Elastooptical correction
n=nCladd+dn;
Ec=1.28; %弹光系数的设定
Rin=a(sCore); %取芯层半径
k0=2*pi/lam;
beta=k0*neff;
n_2=n(sCore+1);
gamma=k0*sqrt(neff^2-n_2^2);
%对后向倏逝场的积分
stepz=5e3;
Startz=5e5;
Nop=Startz/stepz; %积分的点数
rEff=rBend*Ec; %Reff,有效弯曲半径(乘以弹光系数)
%----------------------------------------------------------------
warning('off');
%对后向倏逝场的积分,使用函数trapz。
Hzner=zeros(1,Nop+1);
ze=zeros(1,Nop+1);
for nn=1:Nop+1
zeta=0+(nn-1)*stepz;
X2_Rin=(rEff/2/k0^2/n_2^2)^(2/3)*(beta^2+zeta^2-k0^2*n_2^2*(1+2*Rin/rEff));
X2_0=(rEff/2/k0^2/n_2^2)^(2/3)*(beta^2+zeta^2-k0^2*n_2^2*(1+2*0/rEff));
B1=airy(2,X2_Rin);%艾里函数Bi
A1=airy(0,X2_Rin);%艾里函数Ai
A0=airy(0,X2_0);%LP01 or LP02
G=NGetG_TM(neff,rEff,zeta,lam,a,dn,nCladd,sCore);
H2=pi*exp(-Rin*sqrt(gamma^2+zeta^2))/...
(G*B1+A1)/sqrt(gamma^2+zeta^2);
Hzner(nn) = H2*A0;
ze(nn)=zeta;
end
Int=2*trapz(ze,Hzner); %对各介模式的倏逝场的原式进行积分
alpha=-1/beta/Field2Int*imag(Int)*CO(2,sCore)^2;
%--------------------------------------------------------------
%Bending loss in the units of a dB/m
delta = 10/log(10)*(2*pi*rEff/Ec)*alpha ;% without bendeff more accurate comparing with FEM for BIF!!!!!!
end