design.m

来自「计算光纤模场分布和传播常数」· M 代码 · 共 35 行

function [index] = design(beta,E,r,k0,dr)

format long

num=length(E);
num1=num-2;

temp=beta.^2*k0.^-2;
A=zeros(num1,num1);

for i=2:num1-1
    s=i+1;
    A(i-1,i)=(2*E(s)-E(s+1)+E(s-1))/(2*dr^2);
    A(i,i)=(E(s+1)+E(s-1)-4*E(s))/dr^2+beta^2*E(s);
    A(i,i+1)=(2*E(s)+E(s+1)-E(s-1))/(2*dr^2);
    Bb(i)=temp*((E(s+1)-E(s-1))/(2*dr*r(s))-E(s));
end

A(1,1)=(2*E(2)-E(3)+E(1))/(2*dr^2)+(E(31)+E(1)-4*E(2))/dr^2+beta^2*E(2);
A(1,2)=(2*E(s)+E(3)-E(1))/(2*dr^2);
Bb(1)=temp*((E(3)-E(1))/(2*dr*r(2))-E(2));

A(num1,num1)=(E(num1+2)+E(num1)-4*E(num1+1))/dr^2+beta^2*E(num1+1)+(2*E(num1+1)+E(num1+2)-E(num1))/(2*dr^2);
A(num1-1,num1)=(2*E(num1+1)-E(num1+2)+E(num1))/(2*dr^2);
Bb(num1)=temp*((E(num1+2)-E(num1))/(2*dr*r(num1+1))-E(num1+1));


[L,U]=lu(A);

X=U\(L\Bb');

index=abs(sqrt(X)+sqrt(temp));

plot(index);