//generate a path between points on a sphere
//path given by indices in points(i,1:3) where points are on unit sphere
//measures distance between segments along the arc
//measures angle between segments in local rectangular frame
path=[10,126,25,90,8,76,78,22,89,88,125,32,116,114,26,110,91,111,1,73,51,72,15,57,52,17,53,69,93,28,92,97,9,99,109,108,20,70,71,59,7];//point indices in order of traversal

points_lin=zeros(length(path),2);//start at origin
p1=points(path(1),:);
p2=points(path(2),:);
//theta=acos(sum(p1.*p2)); //angle between unit vectors
//from math.stackexchange.com/questions/1143354/
theta=2*atan(norm(p1-p2),norm(p1+p2)); //better for small angles
dist_arc=theta; //distance along arc
dist_lin=sqrt(sum((p1-p2).^2)); //distance along line
tpl=0;
points_lin(2,:)=points_lin(1,:)+dist_arc*[sin(tpl) cos(tpl)];

for i=2:(length(path)-1)
    p0=points(path(i-1),:);//previous point
    p1=points(path(i),:);
    p2=points(path(i+1),:);
    //theta=acos(sum(p1.*p2)); //angle between unit vectors
    //from math.stackexchange.com/questions/1143354/
    theta=2*atan(norm(p1-p2),norm(p1+p2)); //better for small angles
    dist_arc=theta; //distance along arc
    dist_lin=sqrt(sum((p1-p2).^2)); //distance along line
    
    n0=cross(p0,p1); //previous plane normal of arc
    n1=cross(p1,p2); //next plane normal of arc
    n0=n0/norm(n0); //make unit vectors
    n1=n1/norm(n1);
    cc=sum(n0.*n1); //cosine of angle
    //formula from stackoverflow.com/questions/14066933
    //ss=n0(1)*n1(2)*p1(3)+n1(1)*p1(2)*n0(3)+p1(1)*n0(2)*n1(3)-p1(1)*n1(2)*n0(3)-n0(1)*p1(2)*n1(3)-n1(1)*n0(2)*p1(3); //sine of angle
    nd=cross(n0,n1);
    ss=sum(p1.*nd);
    theta=atan(ss,cc); //negative values means turning clockwise in plane, zero = straight line
    tpl=tpl-theta;
    points_lin(i+1,:)=points_lin(i,:)+dist_arc*[sin(tpl) cos(tpl)];
end

plot(points_lin(:,1),points_lin(:,2),'-o')
isoview
