Help me port this Matlab code......pls....(ported)[Solved]

[toshk]

Code:


#cas
gausss(M1,M3):=
BEGIN
//M1 for linedata: column1=bus From, column2=bus To, column3=Zpu, column4=Shunt Admittance [Q value] but can enter zero if not know , column5=Shunt Admittance [Bpu/2]
//M3 for Busdata: column1=buses, column2=Pg, column3=Qg, column4=PL, column5=QL, column6=Vpu, column7= Vpu Angle, column8=Bus Type
//Bus Type: 1.Slack Bus    2.PQ Bus    3.PV Bus
ss:=(col(M1,5)) .* ; // capacitor   
y:=(col(M1,3));
y:=approx(1 ./ y);
totalbuses:= MAX(MAX(col(M1,1)),MAX(col(M1,2)));
totalbranches:= length(col(M1,1));
M2:=makemat(totalbuses,totalbuses);

FOR b FROM 1 TO totalbranches DO
    M2(M1(b,1),M1(b,2)):=-y(b);
    M2((M1(b,2)),(M1(b,1))):=M2((M1(b,1)),(M1(b,2)));   
END

FOR c FROM 1 TO totalbuses DO
    FOR d FROM 1 TO totalbranches DO
             IF M1(d,1)==c OR M1(d,2) == c THEN
                M2(c,c) := M2(c,c) + y(d)+ss(d);
             END
    END
END
M2;//ybus matrix
busnumber:=col(M3,1);
pg:=col(M3,2);
qg:=col(M3,3);
pl:=col(M3,4);
ql:=col(M3,5);
L2:= col(M3,6);
vv:=L2;
ANG:= col(M3,7);
types:=col(M3,8);

p:= approx((pg - pl) ./100); //Sbase of 100MVA
L0:=approx((qg - ql) ./100);

toll:=1;
iter:=0;
kk:=0.01;//tolerance
alfa:=1.6;

WHILE toll > kk DO
    
    FOR j FROM 2 TO totalbuses DO
    
      YV := 0;
        
        FOR k FROM 1 TO totalbuses DO
            IF j≠k THEN
                YV:= YV + M2(j,k)* L2(k); 
            END;
        YV;
        END;

//Calculating Qi for PV bus
        IF M3(j,8)== 3 THEN
            L0(j) :=  -im(conj(L2(j))*(YV + M2(j,j)*L2(j)));
            M3(j,3):=L0(j);
        END; 
        v1:=L2(j);
        m2:=M2(j,j);
        pu:=p(j);
        qu:=L0(j);
        L2(j):= (1/m2)*((pu-*qu)/conj(v1) - YV); //Compute Bus Voltages.
        
//Calculating Corrected Voltage for PV bus       
       IF M3(j,8)== 3 THEN
           L4(j):=abs(vv(j))*(L2(j)/abs(L2(j)));
           M3(j,6):=L4(j);
           L2(j):=L4(j);
       END;

//Calculating Accelerated Voltage for PQ bus
       IF M3(j,8)== 2 THEN
          L5(j):= vv(j)+alfa*(L2(j)-vv(j));
          M3(j,6):=L5(j);
          L2(j):=L5(j);
       END;
  
     END;
    
iter := iter + 1;
L1:=vv;  
toll:= max(abs(abs(L2)- abs(L1)));   
vv:= L2;
END;
msgbox("check M3 col 6 for |V| and Angle; type M2 for Ybus ")
END;
#end