%%%%%  晶體為Z切。
%%%%%  坐標原點設在波導左端面靠下的一條邊界的中點上,厚度方向為z軸，寬度方向為x軸，光的傳播方向為y軸。

clear
clc
format long;

Numx=15;   %%%% 把波導寬度分為Numx等分
Numz=30;   %%%% 把波導厚度分為Numz等分
%syms x z real
lambda=1.053;
dz=5;  %%%% 擴散深度   由擴散時間和擴散溫度決定
dx=8;  %%% 擴散寬度 由擴散時間和擴散溫度決定
Tao=0.1; %%%% 鈦條初始厚度
w=7;  %%%% 鈦條初始寬度
o_B0=6.53*10^-2-3.15*10^-2*lambda+7.09*10^-3*lambda^2;
o_B1=0.478+0.464*lambda-0.348*lambda^2;
e_B0=0.385-0.430*lambda+0.171*lambda^2;
e_B1=9.13+3.85*lambda-2.49*lambda^2;
ao=0.53;
ae=0.83;
o_deltans=(o_B0+o_B1*Tao/dz)*(Tao/dz)^ao;
e_deltans=(e_B0+e_B1*Tao/dz)*(Tao/dz)^ae;

x=linspace(-dx/2,dx/2,Numx);
z=linspace(0,dz,Numz);
for i=1:size(x,2)
    for j=1:size(z,2)
        u=z(j)/dz;
        s=2*x(i)/w;
        fu=exp(u^2);
        gs=0.5*(erf(w/(2*dx)*(1+s))+erf(w/(2*dx)*(1-s)));
        o_deltanus=o_deltans*fu*gs;
        e_deltanus=e_deltans*fu*gs;
        nop=(4.9048-0.11768/(0.0475-lambda^2-0.027169*lambda^2))^0.5;
        nep=(4.528-0.099169/(0.044432-lambda^2-0.02195*lambda^2))^0.5;
        no(i,j)=nop+o_deltanus*exp(-x(i)^2/dx^2);  %%%% o光
        ne(i,j)=nep+e_deltanus*exp(-x(i)^2/dx^2);  %%%% e光
    end
end
no-ne