function [out,out1,out2]=basis(varargin)

%BASIS  coordinate transformation utilities
%
%	Syntax: see below
%
%       p = basis(x,y) P is the orthonormal 3 by 3 matrix whose first column
%          is in the direction of X, second column is along Y (but orthogonal
%          to X), and third column forms a direct orthonormal basis
%
%	   If x and y are collinear y is selected along the smallest component
%	   of x. A message is passed unless a third argument exists (call of 
%          the form basis(x,y,1))
%
%	[nodeGlobal,bas] = basis(nodeLocal,bas) transforms nodeLocal to a 
%          standard nodeGlobal node definition matrix. BAS is transformed to
%          eliminate all recursive definitions.
%
%	nodeGlobal = basis('gnode',bas,nodeLocal) given an single coordinate
%          system definition BAS, associated nodes nodeLocal are transform to
%          the global coordinate system.
%
%	cGL = basis('trans [ ,l][ ,t]',bas,node,DOF) given a set of coordinate
%          system definitions BAS, the displacement coordinate defintions
%          in NODE(:,3) are used to create a local to global coordinate
%          DOF transformation. The modifier L is used to specify that nodes
%          are given in local rather than global coordinates.
%
%	[p,nodeL]  = basis(node) With 2 outputs, NODE gives coordinates of a 
%          QUAD4 element, P is the local element basis, and NODEL the node 
%          positions in the element coordinate system
%
%	See also : doc basis, node, elt


%	Etienne Balmes
%       Copyright (c) 2001-2003 by INRIA and SDTools, All Rights Reserved.
%       Use under OpenFEM trademark.html license and LGPL.txt library license.

if comstr(varargin{1},'cvs')
 out='$Revision: 1.20 $  $Date: 2006/05/15 20:00:36 $'; return;
end

x=varargin{1}; carg=2;

% basis construction (classical function) - - - - - - - - - - - 
if prod(size(x))==3  

 y = varargin{2}; 

 x=x(:);y=y(:);
 x=x/sqrt(x'*x);
 y=y-(x'*y)*x;
 if norm(y)<eps*norm(x)
   if nargin==2; warning('collinear vectors'); end
   [y,i1]=min(abs(x)); y=zeros(3,1);y(i1(1))=1;y=y-(x'*y)*x;
 end
 y=y/sqrt(y'*y);
 z=-[y(2)*x(3)-y(3)*x(2);y(3)*x(1)-y(1)*x(3);y(1)*x(2)-y(2)*x(1)];
 out=[x y z];

% element basis given faces ----------------------------
elseif nargin==1 %& size(x,1)==4  

 if size(x,2)==7;     x=x(:,5:7);
 elseif size(x,2)~=3; error('Not a valid node matrix');end

 if size(x,1)>3&size(x,1)~=6 % quad (4 or b)
   if size(x,1)==8 x=x(1:4,:);end
   [out,out1,out2]=sp_util('basis',x');

if 1==2
   out2 = mean(x,1);node = x(1:4,:)-out2(ones(4,1),:);
   x = node(3,:)-node(1,:); x=x/sqrt(x*x');
   y = node(4,:)-node(2,:); y=y/sqrt(y*y');
   x = x-y; x=x'/sqrt(x*x');
   y=y';y=y-(x'*y)*x; y=y/sqrt(y'*y);
   z=-[y(2)*x(3)-y(3)*x(2);y(3)*x(1)-y(1)*x(3);y(1)*x(2)-y(2)*x(1)];
   p=[x y z]; out=p; out1=node*p; x=node;
end   

   if nargout>2 %return surface
      %FEw=[0 0 0 4];
      %xi = [-1 -1 0;1 -1 0;1 1 0;-1 1 0];
      %FEnar=[xi(1:4,ones(1,size(FEw,1))*1)'.*(1+FEw(:,2)*xi(:,2)')/4 ];
      %FEnas=[xi(1:4,ones(1,size(FEw,1))*2)'.*(1+FEw(:,1)*xi(:,1)')/4 ];
      FEnar=[-.25 .25 .25 -.25]; FEnas=[-.25 -.25 .25 .25];
      xr = FEnar*x(:,1); xs = FEnas*x(:,1);
      yr = FEnar*x(:,2); ys = FEnas*x(:,2);
      zr = FEnar*x(:,3); zs = FEnas*x(:,3);
      jdet = sqrt((xr.*ys-xs.*yr)^2+(yr.*zs-ys.*zr)^2+(zr.*xs-zs.*xr)^2);
      out2=jdet;

   end

 elseif size(x,1)==3|size(x,1)==6 % triangle
   out2 = mean(x,1); 
   node = x-out2(ones(size(x,1),1),:);
   x = node(2,:)-node(1,:); x=x/sqrt(x*x');
   y = node(3,:)-node(2,:); r1=norm(y);y=y/r1;
   x = (x-y)'; r1=norm(x); 
   if r1==0; [r1,i1]=min(abs(y)); x(i1)=1; else x=x/r1;end
   y=y';y=y-(x'*y)*x; 
   r1=norm(y); if r1==0; [r1,i1]=min(abs(x)); y(i1)=1; else y=y/r1;end
   z=-[y(2)*x(3)-y(3)*x(2);y(3)*x(1)-y(1)*x(3);y(1)*x(2)-y(2)*x(1)];
   p=[x y z]; 
   out=p; out1=node*p;
   if nargout>2 %return surface
      FEw=[1/3 1/3 0 1/2];  % ASSUMES 3 NODE  or straight 6 node TRIANGLE
      xi = [0 0 0;1 0 0;0 1 0];
      FEnar=ones(size(FEw,1),1)*[-1 1 0];
      FEnas=ones(size(FEw,1),1)*[-1 0 1];
      xr = FEnar*node(1:3,1); xs = FEnas*node(1:3,1);
      yr = FEnar*node(1:3,2); ys = FEnas*node(1:3,2);
      zr = FEnar*node(1:3,3); zs = FEnas*node(1:3,3);
      jdet = sqrt((xr.*ys-xs.*yr)^2+(yr.*zs-ys.*zr)^2+(zr.*xs-zs.*xr)^2);
      out2=jdet;
   end
 end


% --------------------------------------------------------------
% transform nodes from local to global -------------------------
elseif (nargin==2&size(x,2)==7)|(ischar(x)&comstr(x,'nodebas'))

if ischar(x)
 model=varargin{carg};carg=carg+1;
 if ~isempty(stack_get(model,'SE')); [FEnode,bas]=fesuper('fnode',model);
 else;FEnode=model.Node; bas=model.bas;
 end
else;
 FEnode=x;  bas=varargin{carg};carg=carg+1;
end
% types : CORD2R 1, CORD2C 2, CORD2S 3, CORD1R 4, CORD1C 5 , CORD1S 6 
% CORD Final is [CID RID(0) Type Origin ex ey ez]

if isa(bas,'cell') 
  r1=[];for j2=1:length(bas);  r1(j2,1:length(bas{j2}))=bas{j2}; end
elseif isa(bas,'double'); r1=bas;
else error('Not a proper basis definition'); end
if isempty(bas); out1=[]; out=FEnode;return;end

if size(r1,2)<15; r1(1,15)=0; end

% work the coordinate changes to global coordinates

%NNode(FEnode(:,1))=1:size(FEnode,1);
NNode=sparse(FEnode(:,1),1,1:size(FEnode,1));

if ~isempty(r1) 

 i4=[1:size(r1,1)]'; i5=find(i4&r1(:,3)==0); % systems that can be transformed
 i1=length(i5); 
 while(any(i4)) % not all have been transformed
 r0=r1;nind=sparse(FEnode(:,2)+2,1:size(FEnode,1),FEnode(:,2), ...
   max(max(r1(:,1)),max(i5))+2,size(FEnode,1));
 for j1 = i5(:)' % loop on systems defined from global

   i2=r1(j1,:);
   % go to everything given form if needed
   if any(i2(7:12))&~any(i2(13:15)) % Ai Bi Ci given
      i3 = basis(i2(7:9)-i2(4:6),i2(10:12)-i2(4:6)); i3 = i3(:,[2 3 1]);
      r1(j1,7:15)=i3(:)';
      i4(j1)=0;
   elseif ~any(i2(7:15)) % node numbers are given
     i2=FEnode(NNode(r1(j1,4:6)),2);
     i2=[find(r1(:,1)==i2(1)) find(r1(:,1)==i2(2)) find(r1(:,1)==i2(3))];
     if isempty(i2)| ...
        (~any(r1(i2,3))&all(i4(i2)==0)) % refering to transformed basis
      x = FEnode(NNode(r1(j1,4:6)),5:7);
      r3 = basis(x(2,:)-x(1,:),x(3,:)-x(1,:)); r3 = r3(:,[2 3 1]);
      r1(j1,1:15) = [r1(j1,1:2)  0 x(1,:) r3(:)'];
      i4(j1)=0;
     end
   else i4(j1)=0; % everything is given
   end

   % transform nodes that depend on this system
   if i4(j1)==0
     ind=find(nind(r1(j1,1)+2,:)); %ind=find(FEnode(:,2)==r1(j1,1));
     %ind=[];if r1(j1,1)<=size(indb,1) ind=find(indb(r1(j1,1),:));end
     if ~isempty(ind)
        FEnode(ind,:)=basis('gnode',r1(j1,:),FEnode(ind,:));
     end
   end

   % transform systems that depend on this one
   for j3= find(r1(:,3)==r1(j1,1))'
    i2=r1(j3,:);
    if any(i2(7:12))&~any(i2(13:15)) % Ai Bi Ci given
     r2=basis('gnode',r1(j1,:),reshape(r1(j3,4:12),3,3)');
     i2=[0 0 0 r2(1,5:7)  r2(2,5:7)  r2(3,5:7)];
     i3 = basis(i2(7:9)-i2(4:6),i2(10:12)-i2(4:6)); i3 = i3(:,[2 3 1]);
     r1(j3,[3:15])=[0 r2(1,5:7) i3(:)'];
    else error('Not a supported basis transformation');
    end
   end

 end % j1 : systems defined from global
 i5=find(i4&r1(:,3)==0);if i1==length(i5)&isequal(r0,r1);
  warning(['Reference Coord Systems ' sprintf(' %i', ...
   unique(r1(find(r1(:,3)),3))) ...
   ' are not defined']);
  break
 else i1=length(i5); end
end;end

out1=r1; out=FEnode;

% --------------------------------------------------------------
% transform nodes given coordinate definition ------------------
elseif ischar(x) & comstr(x,'gnode')

% base definition
% [Id Type RefId OriginXYZ Xvector Yvector Zvector]
r1  = varargin{carg};carg=carg+1; node=varargin{carg};carg=carg+1;
try;out=sp_util('gnode',r1,node);return;end

if size(node,2)==3; node=[zeros(size(node,1),4) node];end

if size(r1,2)<16|r1(16)==0; scale=1; else scale=r1(16);end

if r1(2)==1 % rectangular
  node(:,5:7)=r1(ones(size(node,1),1),4:6)+ ...
           node(:,5:7)*reshape(r1(7:15)*scale,3,3)';
elseif r1(2)==2 % CORD2C
  node(:,5:7)=r1(ones(size(node,1),1),4:6)+ ...
           [node(:,5).*cos(node(:,6)/180*pi) ...
            node(:,5).*sin(node(:,6)/180*pi) node(:,7)] ...
            *reshape(r1(7:15)*scale,3,3)';
elseif r1(2)==3 % CORD2S
  node(:,6:7)=node(:,6:7)*pi/180;
  node(:,5:7)=r1(ones(size(node,1),1),4:6)+ node(:,[5 5 5]).* ...
           [sin(node(:,6)).*cos(node(:,7)) sin(node(:,6)).*sin(node(:,7)) ...
            cos(node(:,6))] ...
            *reshape(r1(7:15)*scale,3,3)';
end
out=node;

% --------------------------------------------------------------
% build the local to global coordinate transformation matrix
elseif ischar(x) & comstr(x,'trans')

[CAM,Cam]=comstr(x,6); if isempty(Cam); Cam=' ';end

bas=varargin{carg};carg=carg+1;
FEnode=varargin{carg};carg=carg+1;
if carg<=nargin; odof=varargin{carg};carg=carg+1; else odof=[]; end

if any(Cam=='l') % nodes are given in local coordinates
  [FEnode,bas]=basis(FEnode,bas);
end

mdof = FEnode(:,1);
if any(Cam=='t')
   mdof=[mdof+.01 mdof+.02 mdof+.03]'; ND=3;
else
   mdof=[mdof+.01 mdof+.02 mdof+.03 mdof+.04 mdof+.05 mdof+.06]';ND=6;
end
mdof=mdof(:);

T1=[ones(length(mdof)*3,2) zeros(length(mdof)*3,1)]; i3=0;
if ~isempty(bas); if size(bas,2)<15; bas(1,15)=0;end;bas=bas(:,1:15);end
bas=[0 1 0  0 0 0  1 0 0  0 1 0  0 0 1;bas];

i1=find(sparse(FEnode(:,3)+1,1,FEnode(:,3)+1))-1;

for j1=i1(:)'  % loop on coordinates
  r1=find(bas(:,1)==j1); if ~isempty(r1); r1 = bas(r1,:);end
  if isempty(r1)
   sdtw('Coordinate sytem %i is undefined',j1);
  else
   if ~any(r1(13:15));  
    error('Basis definitions must be transformed use ''TransL''');
   end
   tr= reshape(r1(7:15),3,3);
   i2 = find(FEnode(:,3)==j1);
   switch r1(2)
   case 1 % rectangular
     ind=ND*(i2-1)+1;in1=i3+[1:length(ind)];
     T1(in1,1:2)=[ind   ind  ];T1(in1,3)=tr(1);in1=in1+length(ind);
     T1(in1,1:2)=[ind+1 ind  ];T1(in1,3)=tr(2);in1=in1+length(ind);
     T1(in1,1:2)=[ind+2 ind  ];T1(in1,3)=tr(3);in1=in1+length(ind);
     T1(in1,1:2)=[ind   ind+1];T1(in1,3)=tr(4);in1=in1+length(ind);
     T1(in1,1:2)=[ind+1 ind+1];T1(in1,3)=tr(5);in1=in1+length(ind);
     T1(in1,1:2)=[ind+2 ind+1];T1(in1,3)=tr(6);in1=in1+length(ind);
     T1(in1,1:2)=[ind   ind+2];T1(in1,3)=tr(7);in1=in1+length(ind);
     T1(in1,1:2)=[ind+1 ind+2];T1(in1,3)=tr(8);in1=in1+length(ind);
     T1(in1,1:2)=[ind+2 ind+2];T1(in1,3)=tr(9);
     if ND==6
      in1=in1+length(ind);ind=ind+3;
      T1(in1,1:2)=[ind   ind  ];T1(in1,3)=tr(1);in1=in1+length(ind);
      T1(in1,1:2)=[ind+1 ind  ];T1(in1,3)=tr(2);in1=in1+length(ind);
      T1(in1,1:2)=[ind+2 ind  ];T1(in1,3)=tr(3);in1=in1+length(ind);
      T1(in1,1:2)=[ind   ind+1];T1(in1,3)=tr(4);in1=in1+length(ind);
      T1(in1,1:2)=[ind+1 ind+1];T1(in1,3)=tr(5);in1=in1+length(ind);
      T1(in1,1:2)=[ind+2 ind+1];T1(in1,3)=tr(6);in1=in1+length(ind);
      T1(in1,1:2)=[ind   ind+2];T1(in1,3)=tr(7);in1=in1+length(ind);
      T1(in1,1:2)=[ind+1 ind+2];T1(in1,3)=tr(8);in1=in1+length(ind);
      T1(in1,1:2)=[ind+2 ind+2];T1(in1,3)=tr(9);
     end
     i3=in1(end);
   case 2 % cylindrical
     for j2=i2(:)'
      r2=(FEnode(j2,5:7)-r1(4:6))*tr;
      r3=atan2(r2(2),r2(1)); 
      r3=tr*[cos(r3) -sin(r3) 0;sin(r3) cos(r3) 0;0 0 1];
      T1(i3+[1:9],1:3) = [ND*(j2-1)+1+[0;1;2;0;1;2;0;1;2] ...
                          ND*(j2-1)+1+[0;0;0;1;1;1;2;2;2] r3(:)];
      i3=i3+9;
      if ND==6 
        T1(i3+[1:9],1:3) = [ND*(j2-1)+4+[0;1;2;0;1;2;0;1;2] ...
                            ND*(j2-1)+4+[0;0;0;1;1;1;2;2;2] r3(:)];
      i3=i3+9;
      end
     end
   case 3 % spherical
     error(1)
     for j2=i2(:)'
      r2=(FEnode(j2,5:7)-r1(4:6))*tr;
      r3=atan2(r2(2),r2(1)); 
      r3=tr*[cos(r3) -sin(r3) 0;sin(r3) cos(r3) 0;0 0 1];
      ind = (j2-1)*6+[1:3]; T1(ind,ind)=r3; ind=ind+3;T1(ind,ind)=r3;
     end
  otherwise error('Unsupported coordinate transformation');
  end
 end % the basis was found
end

T1=sparse(T1(:,1),T1(:,2),T1(:,3));
if      isempty(odof); out=T1;
elseif size(odof)==size(mdof) & norm(odof-mdof)<.01; out=T1;
else   out=fe_c(odof,mdof,T1); end
out1=mdof;

% --------------------------------------------------------------

end

if 1==2 % INFO ON LMS EULER ANGLE COORDINATE TRANSFORM


%      [psi teta phi]  also named teta_xy \theta_{xz} \theta_{yz}
% in order rotation around global Z, rotation around intermediate y,
% rotation around intermediate x

% nast [    phi  teta] 

a=[0 0 15]/180*pi;
cA=cos(a(1)); cB=cos(a(2)); cC=cos(a(3));
sA=sin(a(1)); sB=sin(a(2)); sC=sin(a(3));

T=[cA*cB     -sA*cC-cA*sB*sC      sA*sC-cA*sB*cC;
   sA*cB      cA*cC-sA*sB*sC     -cA*sC-sA*sB*cC;
   sB         cB*sC               cB*cC;   
  ];

end