function [Tc,a,b]=sixim(g,op)


% g: shift angle between star windings (s1 and s2)
% op: vector containing open phase(s) number(s)


if nargin<2
    g=pi/3;
    op=[0 0 0 0];
end,
if length(op)>4
    op=op(1:4);
end,


al0=[1 cos(g) cos(2*pi/3) cos(2*pi/3+g) cos(4*pi/3) cos(4*pi/3+g)];
bt0=[0 sin(g) sin(2*pi/3) sin(2*pi/3+g) sin(4*pi/3) sin(4*pi/3+g)];
phi=[0  g  2*pi/3  2*pi/3+g  4*pi/3  4*pi/3+g];


kk1=op(1);kk2=op(2);kk3=op(3);kk4=op(4);
num=0;den=1e-12;
for i=1:6
    if ((i==kk1)|(i==kk2)|(i==kk3)|(i==kk4))
        num=num+sin(2*phi(i));
        den=den+cos(2*phi(i));
    end,
end,
Th0=-atan(num/den)/2;


al=[cos(0+Th0) cos(g+Th0) cos(2*pi/3+Th0) cos(2*pi/3+g+Th0) cos(4*pi/3+Th0) cos(4*pi/3+g+Th0)];
bt=[sin(0+Th0) sin(g+Th0) sin(2*pi/3+Th0) sin(2*pi/3+g+Th0) sin(4*pi/3+Th0) sin(4*pi/3+g+Th0)];


j=1;
for i=1:6
    if ((i~=kk1)&(i~=kk2)&(i~=kk3)&(i~=kk4))
            Th(j)=phi(i);
        t1(j)=al(i);
        t2(j)=bt(i);
        j=j+1;
    end,
end,
t1c=t1;
t2c=t2;


a=norm(t1c);
b=norm(t2c);
Tc=[t1c/a;t2c/b];


n=length(t1);
Tn_=[t1;t2];
if n>2
    Tz=null(Tn_)';
    Tn=[t1/norm(t1);t2/norm(t2);Tz];
else
    Tn=[t1/norm(t1);t2/norm(t2)];
end,