function [x,y,z]=nvrtrf(u,v,w,pow,A,tol)
%
% function [x,y,z]=nvrtrf(u,v,w,pow,A,tol)
%
% A one-parameter trf u,v,w => x,y,z that maps a face of 
% the octahedron (u+v+w=1) to an octant of the sphere is:
%
%   bx=1-(1-u)^pow, by=1-(1-v)^pow, bz=1-(1-w)^pow,
%   xyz=sqrt(bx^2+by^2+bz^2), x=bx/xyz, y=by/xyz, z=bz/xyz
% 
% An approximate two-param inverse x,y,z => u,v,w is:
%
%    bu=1-(1-alf*x)^q          x^2+y^2+z^2 = 1
%    bv=1-(1-alf*y)^q  where   q=1/pow
%    bw=1-(1-alf*z)^q          alf = (x+y+z)-A*(xy+yz+zx)      
%
%    uvw=bu+bv+bw and u=bu/uvw, v=bv/uvw, w=bw/uvw
%
% What we do:   given (u,v,w) with u+v+w=1,
%               use Newton' method to find the (x,y,z) 
%               s.t. x,y,z => u,v,w (the approx inverse)
%
% if not specified, pow is 1.55         (try 0.33)
% if not specified, A is 1         (not bad for 1<pow<2)
% if not specified, tol is 1e-4   (max of |du|,|dv|,|dw|)
%

if nargin<3 help nvrtrf; return; end
if nargin<4 pow=1.55; end
if nargin<5 A=1; end
if nargin<6 tol=1e-4; end

B=5*(pow-1.55)/3; if B<=-1 B=-1; end;          % for first approximation 
s=w; w=abs(s);                                 % undo later

% now find x,y,z s.t approximate inverse => given u,v,w

ndx=find(s>eps&w+eps>u&u+eps>v)';    % do these, about 1/12 of the quadrant

% first approximation of x,y,z - see pnodes()

p=pow; q=1/p; r=1-q;
uu=u(ndx); vv=v(ndx); ww=w(ndx);                  % need them all 
bta=1-B*(uu.*vv+vv.*ww+ww.*uu);                   % better than bta=1
bx=1-(1-bta.*uu).^p;                              % x-bar
by=1-(1-bta.*vv).^p;                              % y-bar
bz=1-(1-bta.*ww).^p;                              % z-bar
alf=sqrt(bx.*bx+by.*by+bz.*bz);                   % the direct
xx=bx./alf; yy=by./alf; zz=bz./alf;               % transformation

% iterate to within tol

while 1    

   tta=xx+yy+zz; alf=tta-(tta.^2-1)*A/2;          % theta=x+y+z 
   alt=1-tta; tatpa=tta.*alt+alf;                 % tta*alf'(tta)+alf(tta)

   bu=1-(1-alf.*xx).^q; bv=1-(1-alf.*yy).^q;      % u-bar, ...
   bw=1-(1-alf.*zz).^q; bta=bu+bv+bw;             % done when bu=bta*u, ...

   if max(abs(ww-bw./bta))<tol break; end

   c11=xx.*alt-uu.*tatpa;                         % Newton - first row
   c12=(1-yy./zz).*c11; c11=alf+(1-xx./zz).*c11;

   c22=yy.*alt-vv.*tatpa;                         % Newton - second row
   c21=(1-xx./zz).*c22; c22=alf+(1-yy./zz).*c22;

   den=c11.*c22-c12.*c21;                         % far from zero

   s1=p*(1-alf.*xx).^r.*(bta.*uu-bu);             % first row
   s2=p*(1-alf.*yy).^r.*(bta.*vv-bv);             % second    

   dx=(c22.*s1-c12.*s2)./den; xx=xx+dx;           % by
   dy=(c11.*s2-c21.*s1)./den; yy=yy+dy;           % hand

   zz=sqrt(1-xx.*xx-yy.*yy);

end

% now fill in the rest, all rearrangements of what has been done

N=size(u); if N(1)~=N(2) return; end              % don't know what to do

NN=N(1)*N(2); kgt=reshape(1:NN,N);                % the key to it all
klr=fliplr(kgt);                                  % for permutations
ktr=kgt';                                         % ordinary transpose
kgt=reshape(NN:-1:1,N)';                          % goofy transpose

x=zeros(N); y=x; z=y;                             % 1/4 of the sphere 
x(ndx)=xx; y(ndx)=yy; z(ndx)=zz;                  % start to fill them in

k1=find(s>eps&w+eps>u&u>v+eps); k2=kgt(k1);       % reflect x <=> y 
z(k2)=z(k1); y(k2)=x(k1); x(k2)=y(k1);

k1=ktr(find(s<eps&u>w+eps&u+eps>v));                % this was hard  
k2=klr(sort(klr(find(s+eps>0&w>v+eps&w+eps>u))));   % really hard
x(k1)=z(k2); z(k1)=y(k2); y(k1)=x(k2);              % the permutation
 
k1=find(s+eps>0&u>v+eps&u>w+eps); k2=kgt(k1);     % reflect x <=> y 
z(k2)=z(k1); y(k2)=x(k1); x(k2)=y(k1);

k1=find(s>eps); k2=ktr(k1);                       % at last
x(k2)=x(k1); y(k2)=y(k1); z(k2)=-z(k1);           % the easy one