% Copyright 1996 Thomas Mc Avoy and Dong Dong

 function [lmd,res1,distnew,kkk]=pjord1(f,data,span,kkk);

% function to project data to f and reorder data according to the 
% projection.
% lmd is the projection index, distnew gives distance from data to 
% principal curse, and kkk keeps tract of data ordering
% data is ordered based on last calculation
% function [lmd,res1,distnew,kkk]=pjorder(f,data,span,kkk)
	[m,n]=size(data);
	span=round(span/2); % half the span points are checked as nearest points
% this speeds up the algorithm
	h=ones(size(f(1:span+1,:)));
    distnew=0.;
	lmd1=zeros(1,m);
% calculate distance down curve, lmd
        for i=2:m
           lmd1(i)=lmd1(i-1)+norm(f(i,:)-f(i-1,:));
        end
% now order points. try neighboring points within span/2 of each point
        for i=1:m
			if i<=round(span/2)
           	distf=sum(((f(1:span+1,:)-h*diag(data(i,:))).^2)')';
			[a2,b2]=min(distf);
		elseif i>=round(span/2)&i<=m-round(span/2)
			distf=sum(((f(i-span/2:i+span/2,:)-h*diag(data(i,:))).^2)')';
			[a2,b2]=min(distf);
		else
			distf=sum(((f(m-span:m,:)-h*diag(data(i,:))).^2)')';
	   		[a2,b2]=min(distf);
		end
% add min distance to performance index
	    distnew=distnew+a2;
% estimate lmd by linear fit
           aa=-f(b2,:)+data(i,:);
           if b2 < m
                bb=(-f(b2,:)+f(b2+1,:))/(lmd1(b2+1)-lmd1(b2));
           end
           if b2==m
                bb=(-f(b2-1,:)+f(b2,:))/(lmd1(b2)-lmd1(b2-1));
           end
	    if norm(bb)==0
	       d(i)=0;
	    else
           d(i)=(aa*bb')/norm(bb)^2;
	   	end
           lmd(i)=lmd1(b2)+d(i);
        end
        lmd=lmd';
        [lmd,record]=sort(lmd);
		res1=data(record,:);% res give data sorted by lmd
		kkk=kkk(record);% keep track of sorted data. need to unsort later
        distnew=distnew/(m*n);
% routine returns lmd, data sorted by lmd, sort kkk, and distnew


