moon_diffu.m

用fdtd方法 模拟集成光学离子交换中扩散离子浓度的分布 进而求得折射率的分布。更创新性加入了 可变扩散系数 条件（绝对原创）

function diffus=moon_diffu(halfwindow,t,high,halfwidth)
h=1;
half=halfwindow/h;m=high/h;n=halfwidth/h;

M0=zeros(m,n);
M0(1,[1:half])=1;
if t<=600
   for k=1:t
      D=DC_moondiffu(M0);
      DB=DC_moondiffu_b(M0);
      alfa=alfa_compute(D,DB);
      M(1,[1:half])=M0(1,[1:half]);
      M(1,[half+1:n-1])=M0(1,[half+1:n-1])+Side1_AA(M0,1,[half+1:n-1],1,D,DB,alfa,h);
      M(1,n)=M0(1,n)+corner1_AA(M0,1,n,1,D,DB,alfa,h);
      M([2:m-1],n)=M0([2:m-1],n)+Side2_AA(M0,[2:m-1],n,1,D,DB,alfa,h);
      M(m,n)=M0(m,n)+corner2_AA(M0,m,n,1,D,DB,alfa,h);
      M(m,[2:n-1])=M0(m,[2:n-1])+Side3_AA(M0,m,[2:n-1],1,D,DB,alfa,h);
      M(m,1)=M0(m,1)+corner3_AA(M0,m,1,1,D,DB,alfa,h);
      M([2:m-1],1)=M0([2:m-1],1)+Side4_AA(M0,[2:m-1],1,1,D,DB,alfa,h);
      M([2:m-1],[2:n-1])=M0([2:m-1],[2:n-1])+mid_AA(M0,[2:m-1],[2:n-1],1,D,DB,alfa,h);
      M0=M;
      diff_time=k
   end
   diffus=M0;
   ALFA=alfa
else if t<=43200
    for k=1:600
     D=DC_moondiffu(M0);
     DB=DC_moondiffu_b(M0);
     alfa=alfa_compute(D,DB);
     M(1,[1:half])=M0(1,[1:half]);
     M(1,[half+1:n-1])=M0(1,[half+1:n-1])+Side1_AA(M0,1,[half+1:n-1],1,D,DB,alfa,h);
     M(1,n)=M0(1,n)+corner1_AA(M0,1,n,1,D,DB,alfa,h);
     M([2:m-1],n)=M0([2:m-1],n)+Side2_AA(M0,[2:m-1],n,1,D,DB,alfa,h);
     M(m,n)=M0(m,n)+corner2_AA(M0,m,n,1,D,DB,alfa,h);
     M(m,[2:n-1])=M0(m,[2:n-1])+Side3_AA(M0,m,[2:n-1],1,D,DB,alfa,h);
     M(m,1)=M0(m,1)+corner3_AA(M0,m,1,1,D,DB,alfa,h);
     M([2:m-1],1)=M0([2:m-1],1)+Side4_AA(M0,[2:m-1],1,1,D,DB,alfa,h);
     M([2:m-1],[2:n-1])=M0([2:m-1],[2:n-1])+mid_AA(M0,[2:m-1],[2:n-1],1,D,DB,alfa,h);
     M0=M;
     diff_time=k
 end
  for k=600:5:t
       D=DC_moondiffu(M0);
     DB=DC_moondiffu_b(M0);
     alfa=alfa_compute(D,DB);
     M(1,[1:half])=M0(1,[1:half]);
     M(1,[half+1:n-1])=M0(1,[half+1:n-1])+Side1_AA(M0,5,[half+1:n-1],1,D,DB,alfa,h);
     M(1,n)=M0(1,n)+corner1_AA(M0,1,n,5,D,DB,alfa,h);
     M([2:m-1],n)=M0([2:m-1],n)+Side2_AA(M0,[2:m-1],n,5,D,DB,alfa,h);
     M(m,n)=M0(m,n)+corner2_AA(M0,m,n,5,D,DB,alfa,h);
     M(m,[2:n-1])=M0(m,[2:n-1])+Side3_AA(M0,m,[2:n-1],5,D,DB,alfa,h);
     M(m,1)=M0(m,1)+corner3_AA(M0,m,1,5,D,DB,alfa,h);
     M([2:m-1],1)=M0([2:m-1],1)+Side4_AA(M0,[2:m-1],1,5,D,DB,alfa,h);
     M([2:m-1],[2:n-1])=M0([2:m-1],[2:n-1])+mid_AA(M0,[2:m-1],[2:n-1],5,D,DB,alfa,h);
     M0=M;
     diff_time=k
 end
  diffus=M0;
end
end