关于利用矩阵***术求解节点导纳矩阵的MATLAB编程

求大神帮忙看一下这个程序，解释一下！
%形成节点导纳矩阵
%考虑参考节点
%考虑互感支路
function [d_g,d_b,u_g,u_b]=calyy(itoi,noels_hgln,noelsnum_hgln,bs_noels,bs_cp,gb_hgln,gb_noels,gb_cp,lv_bs,lv_hgln,lv_noels,lv_cp,type_bs,iy,jy,remv_noels)
%lv_bs
d_g=zeros(1,lv_bs);
d_b=zeros(1,lv_bs);
u_g=zeros(1,length(jy));
u_b=zeros(1,length(jy));
type_noels=zeros(1,lv_noels);
for i=1:lv_hgln
    %%支路自阻抗
    lnnum=noelsnum_hgln(i);
    lnnum;
    gb_noels;
    for ii1=1:lnnum
        iln=noels_hgln(i,ii1);
        [d_g,d_b,u_g,u_b]=getlnselfgb(iln,bs_noels,type_bs,itoi,gb_noels,iy,jy,d_g,d_b,u_g,u_b,lv_bs);
        type_noels(iln)=1;
    end
    for ii1=1:lnnum
        iln=noels_hgln(i,ii1);
        for ii2=ii1+1:lnnum
            zln=noels_hgln(i,ii2);
            b=gb_hgln(i,ii1,ii2);
            [d_g,d_b,u_g,u_b]=getlnhugb(iln,zln,bs_noels,type_bs,itoi,0,b,iy,jy,d_g,d_b,u_g,u_b,lv_bs);
        end
    end
end
for i=1:lv_noels
    if remv_noels(i)==1%被除去线路
        continue;
    end
    ibs=bs_noels(i,1);
    zbs=bs_noels(i,2);
    if type_noels(i)==1
        continue;
    end
    type_noels(i)=1;
    iibs=itoi(ibs);
    zzbs=itoi(zbs);
    g=gb_noels(i,1);
    b=gb_noels(i,2);
    j=0;
    if iibs~=zzbs
        j=getlocal(iy,jy,iibs,zzbs,lv_bs,1);
    end
    bch=gb_noels(i,3);
    %分接头所在端的自导纳要+=（g+***）/t^2
    if bch<0
        bch=-bch;
        if gb_noels(i,4)==1%变比位于末端
            ibs=iibs;
            zbs=zzbs;
        else               %变比位于首段
            ibs=zzbs;
            zbs=iibs;
        end
        if ibs<=lv_bs
            d_b(ibs)=d_b(ibs)+b;
            d_g(ibs)=d_g(ibs)+g;
        end
        if zbs<=lv_bs
            d_b(zbs)=d_b(zbs)+b/bch/bch;
            d_g(zbs)=d_g(zbs)+g/bch/bch;
        end
        if j~=0
            u_g(j)=u_g(j)-g/bch;
            u_b(j)=u_b(j)-b/bch;
        end
        continue;
    end
    %disp('线路支路')；
    if iibs<=lv_bs
        d_b(iibs)=d_b(iibs)+b;
        d_b(iibs)=d_b(iibs)+bch;
        d_g(iibs)=d_g(iibs)+g;
    end
    if zzbs<=lv_bs
        d_b(zzbs)=d_b(zzbs)+b;
        d_b(zzbs)=d_b(zzbs)+bch;
        d_g(zzbs)=d_g(zzbs)+g;   
    end
    if j~=0
        u_g(j)=u_g(j)-g;
        u_b(j)=u_b(j)-b;  
    end
end
for i=1:lv_cp
    ibs=bs_cp(i,1);
    if ibs==0
        continue;
    end
    headbs=itoi(ibs);
    if headbs>lv_bs
        continue;
    end
    d_b(headbs)=d_b(headbs)+gb_cp(i,1);
end
k=iy(lv_bs);
u_g=u_g(1:k-1);
u_b=u_b(1:k-1);


function [d_g,d_b,u_g,u_b]=getlnselfgb(iln,bs_noels,type_bs,itoi,gb_noels,iy,jy,d_g,d_b,u_g,u_b,lv_bs)
g=gb_noels(iln,1);
b=gb_noels(iln,2);
bch=gb_noels(iln,3);
if bch<0
    disp('error!!!');
    return;
end
ibs=bs_noels(iln,1);
zbs=bs_noels(iln,2);
iibs=itoi(ibs);
zzbs=itoi(zbs);
if iibs~=lv_bs+1
    d_b(iibs)=d_b(iibs)+b+bch;
    d_g(iibs)=d_g(iibs)+g;
end
if zzbs~=lv_bs+1
    d_b(zzbs)=d_b(zzbs)+b+bch;
    d_g(zzbs)=d_g(zzbs)+g;
end

j=0;
if iibs~=zzbs
    j=getlocal(iy,jy,iibs,zzbs,lv_bs,1);
end

if j==0
    disp('找不到支路位置');
    %exit；
    return;
end
u_g(j)=u_g(j)-g;
u_b(j)=u_b(j)-b;


function [d_g,d_b,u_g,u_b]=getlnhugb(iln,zln,bs_noels,type_bs,itoi,g,b,iy,jy,d_g,d_b,u_g,u_b,lv_bs)
  ibsi=bs_noels(iln,1);%首
  zbsi=bs_noels(iln,2);%末
  ibsz=bs_noels(zln,1);%首
  zbsz=bs_noels(zln,2);%末
  first_head=itoi(ibsi);
  first_tail=itoi(zbsi);
  second_head=itoi(ibsz);
  second_tail=itoi(zbsz);
  if first_head==second_head&&first_head~=lv_bs+1
      disp('add');
      d_b(first_head)=d_b(first_head)+2*b;
      d_g(first_head)=d_g(first_head)+2*g;
  end
  if first_tail==second_tail&&first_tail~=lv_bs+1
      disp('add');
      d_b(first_tail)=d_b(first_tail)+2*b;
      d_g(first_tail)=d_g(first_tail)+2*g;
  end
  if first_head==second_tail&&first_head~=lv_bs+1
      disp('add');
      d_b(first_head)=d_b(first_head)-2*b;
      d_g(first_head)=d_g(first_head)-2*g;
  end
  if first_tail==second_head&&first_tail~=lv_bs+1
      disp('add');
      d_b(first_tail)=d_b(first_tail)-2*b;
      d_g(first_tail)=d_g(first_tail)-2*g;
  end

  [d_g,d_b,u_g,u_b]=addlnval(first_head,second_head,g,b,iy,jy,d_g,d_b,u_g,u_b,lv_bs);
  [d_g,d_b,u_g,u_b]=addlnval(first_head,second_tail,-g,-b,iy,jy,d_g,d_b,u_g,u_b,lv_bs);
  [d_g,d_b,u_g,u_b]=addlnval(second_head,first_tail,-g,-b,iy,jy,d_g,d_b,u_g,u_b,lv_bs);
  [d_g,d_b,u_g,u_b]=addlnval(second_tail,first_tail,g,b,iy,jy,d_g,d_b,u_g,u_b,lv_bs);

function [d_g,d_b,u_g,u_b]=addlnval(first_head,second_head,g,b,iy,jy,d_g,d_b,u_g,u_b,lv_bs)
    if first_head==second_head
        return;
    end
    j=getlocal(iy,jy,first_head,second_head,lv_bs,1);
    if j==0
        return;
    end
    u_b(j)=u_b(j)+b;
    u_g(j)=u_g(j)+g;