[ix,iy,n,kflag]=isort(ix,iy,n,kflag);

function [ix,iy,n,kflag]=isort(ix,iy,n,kflag);
%***BEGIN PROLOGUE  ISORT
%***PURPOSE  Sort an array and optionally make the same interchanges in
%            an auxiliary array.  The array may be sorted in increasing
%            or decreasing order.  A slightly modified QUICKSORT
%            algorithm is used.
%***LIBRARY   SLATEC
%***CATEGORY  N6A2A
%***TYPE      INTEGER (SSORT-S, DSORT-D, ISORT-I)
%***KEYWORDS  SINGLETON QUICKSORT, SORT, SORTING
%***AUTHOR  Jones, R. E., (SNLA)
%           Kahaner, D. K., (NBS)
%           Wisniewski, J. A., (SNLA)
%***DESCRIPTION
%
%   ISORT sorts array IX and optionally makes the same interchanges in
%   array IY.  The array IX may be sorted in increasing order or
%   decreasing order.  A slightly modified quicksort algorithm is used.
%
%   Description of Parameters
%      IX - integer array of values to be sorted
%      IY - integer array to be (optionally) carried along
%      N  - number of values in integer array IX to be sorted
%      KFLAG - control parameter
%            =  2  means sort IX in increasing order and carry IY along.
%            =  1  means sort IX in increasing order (ignoring IY)
%            = -1  means sort IX in decreasing order (ignoring IY)
%            = -2  means sort IX in decreasing order and carry IY along.
%
%***REFERENCES  R. C. Singleton, Algorithm 347, An efficient algorithm
%                 for sorting with minimal storage, Communications of
%                 the ACM, 12, 3 (1969), pp. 185-187.
%***ROUTINES CALLED  XERMSG
%***REVISION HISTORY  (YYMMDD)
%   761118  DATE WRITTEN
%   810801  Modified by David K. Kahaner.
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   891009  Removed unreferenced statement labels.  (WRB)
%   891009  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
%   901012  Declared all variables; changed X,Y to IX,IY. (M. McClain)
%   920501  Reformatted the REFERENCES section.  (DWL, WRB)
%   920519  Clarified error messages.  (DWL)
%   920801  Declarations section rebuilt and code restructured to use
%           IF-THEN-ELSE-ENDIF.  (RWC, WRB)
%***end PROLOGUE  ISORT
%     .. Scalar Arguments ..
persistent gt i ij il iu j k kk l m nn r t tt tty ty ; 

if isempty(gt), gt=zeros(1,3); end;
%     .. Array Arguments ..
ix_shape=size(ix);ix=reshape(ix,1,[]);
iy_shape=size(iy);iy=reshape(iy,1,[]);
%     .. Local Scalars ..
if isempty(r), r=0; end;
if isempty(i), i=0; end;
if isempty(ij), ij=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(kk), kk=0; end;
if isempty(l), l=0; end;
if isempty(m), m=0; end;
if isempty(nn), nn=0; end;
if isempty(t), t=0; end;
if isempty(tt), tt=0; end;
if isempty(tty), tty=0; end;
if isempty(ty), ty=0; end;
%     .. Local Arrays ..
if isempty(il), il=zeros(1,21); end;
if isempty(iu), iu=zeros(1,21); end;
%     .. External Subroutines ..
%     .. Intrinsic Functions ..
%***FIRST EXECUTABLE STATEMENT  ISORT
nn = fix(n);
if( nn<1 )
xermsg('SLATEC','ISORT','The number of values to be sorted is not positive.',1,1);
ix_shape=zeros(ix_shape);ix_shape(:)=ix(1:numel(ix_shape));ix=ix_shape;
iy_shape=zeros(iy_shape);iy_shape(:)=iy(1:numel(iy_shape));iy=iy_shape;
return;
else;
%
kk = fix(abs(kflag));
if( kk~=1 && kk~=2 )
xermsg('SLATEC','ISORT','The sort control parameter, K, is not 2, 1, -1, or -2.',2,1);
ix_shape=zeros(ix_shape);ix_shape(:)=ix(1:numel(ix_shape));ix=ix_shape;
iy_shape=zeros(iy_shape);iy_shape(:)=iy(1:numel(iy_shape));iy=iy_shape;
return;
else;
%
%     Alter array IX to get decreasing order if needed
%
if( kflag<=-1 )
for i = 1 : nn;
ix(i) = fix(-ix(i));
end; i = fix(nn+1);
end;
%
if( kk==2 )
%
%     Sort IX and carry IY along
%
m = 1;
i = 1;
j = fix(nn);
r = 0.375e0;
%
gt(:)=0;
while( true );
if(gt(3)==0)
if(gt(2)==0)
if(gt(1)==0)
if( i==j )
gt(3)=1;
continue;
end;
if( r<=0.5898437e0 )
r = r + 3.90625e-2;
else;
r = r - 0.21875e0;
end;
%
end;
gt(1)=0;
k = fix(i);
%
%     Select a central element of the array and savemlv it in location T
%
ij = fix(i + fix((j-i).*r));
t = fix(ix(ij));
ty = fix(iy(ij));
%
%     If first element of array is greater than T, interchange with T
%
if( ix(i)>t )
ix(ij) = fix(ix(i));
ix(i) = fix(t);
t = fix(ix(ij));
iy(ij) = fix(iy(i));
iy(i) = fix(ty);
ty = fix(iy(ij));
end;
l = fix(j);
%
%     If last element of array is less than T, interchange with T
%
if( ix(j)<t )
ix(ij) = fix(ix(j));
ix(j) = fix(t);
t = fix(ix(ij));
iy(ij) = fix(iy(j));
iy(j) = fix(ty);
ty = fix(iy(ij));
%
%        If first element of array is greater than T, interchange with T
%
if( ix(i)>t )
ix(ij) = fix(ix(i));
ix(i) = fix(t);
t = fix(ix(ij));
iy(ij) = fix(iy(i));
iy(i) = fix(ty);
ty = fix(iy(ij));
end;
end;
%
%     Find an element in the second half of the array which is smaller
%     than T
%
while( true );
l = fix(l - 1);
if( ix(l)<=t )
%
%     Find an element in the first half of the array which is greater
%     than T
%
while( true );
k = fix(k + 1);
if( ix(k)>=t )
break;
end;
end;
%
%     Interchange these elements
%
if( k>l )
break;
end;
tt = fix(ix(l));
ix(l) = fix(ix(k));
ix(k) = fix(tt);
tty = fix(iy(l));
iy(l) = fix(iy(k));
iy(k) = fix(tty);
end;
end;
%
%     Save upper and lower subscripts of the array yet to be sorted
%
if( l-i>j-k )
il(m) = fix(i);
iu(m) = fix(l);
i = fix(k);
m = fix(m + 1);
else;
il(m) = fix(k);
iu(m) = fix(j);
j = fix(l);
m = fix(m + 1);
end;
%
end;
gt(2)=0;
if( j-i>=1 )
gt(1)=1;
continue;
end;
if( i==1 )
continue;
end;
i = fix(i - 1);
%
while( true );
i = fix(i + 1);
if( i==j )
break;
end;
t = fix(ix(i+1));
ty = fix(iy(i+1));
if( ix(i)>t )
k = fix(i);
%
while( true );
ix(k+1) = fix(ix(k));
iy(k+1) = fix(iy(k));
k = fix(k - 1);
if( t>=ix(k) )
break;
end;
end;
ix(k+1) = fix(t);
iy(k+1) = fix(ty);
end;
end;
%
%     Begin again on another portion of the unsorted array
%
end;
gt(3)=0;
m = fix(m - 1);
if( m==0 )
break;
end;
i = fix(il(m));
j = fix(iu(m));
gt(2)=1;
end;
if( kflag<=-1 )
for i = 1 : nn;
ix(i) = fix(-ix(i));
end; i = fix(nn+1);
end;
ix_shape=zeros(ix_shape);ix_shape(:)=ix(1:numel(ix_shape));ix=ix_shape;
iy_shape=zeros(iy_shape);iy_shape(:)=iy(1:numel(iy_shape));iy=iy_shape;
return;
else;
%
%     Sort IX only
%
m = 1;
i = 1;
j = fix(nn);
r = 0.375e0;
end;
end;
end;
%
gt(:)=0;
while( true );
if(gt(3)==0)
if(gt(2)==0)
if(gt(1)==0)
if( i==j )
gt(3)=1;
continue;
end;
if( r<=0.5898437e0 )
r = r + 3.90625e-2;
else;
r = r - 0.21875e0;
end;
%
end;
gt(1)=0;
k = fix(i);
%
%     Select a central element of the array and savemlv it in location T
%
ij = fix(i + fix((j-i).*r));
t = fix(ix(ij));
%
%     If first element of array is greater than T, interchange with T
%
if( ix(i)>t )
ix(ij) = fix(ix(i));
ix(i) = fix(t);
t = fix(ix(ij));
end;
l = fix(j);
%
%     If last element of array is less than than T, interchange with T
%
if( ix(j)<t )
ix(ij) = fix(ix(j));
ix(j) = fix(t);
t = fix(ix(ij));
%
%        If first element of array is greater than T, interchange with T
%
if( ix(i)>t )
ix(ij) = fix(ix(i));
ix(i) = fix(t);
t = fix(ix(ij));
end;
end;
%
%     Find an element in the second half of the array which is smaller
%     than T
%
while( true );
l = fix(l - 1);
if( ix(l)<=t )
%
%     Find an element in the first half of the array which is greater
%     than T
%
while( true );
k = fix(k + 1);
if( ix(k)>=t )
break;
end;
end;
%
%     Interchange these elements
%
if( k>l )
break;
end;
tt = fix(ix(l));
ix(l) = fix(ix(k));
ix(k) = fix(tt);
end;
end;
%
%     Save upper and lower subscripts of the array yet to be sorted
%
if( l-i>j-k )
il(m) = fix(i);
iu(m) = fix(l);
i = fix(k);
m = fix(m + 1);
else;
il(m) = fix(k);
iu(m) = fix(j);
j = fix(l);
m = fix(m + 1);
end;
%
end;
gt(2)=0;
if( j-i>=1 )
gt(1)=1;
continue;
end;
if( i==1 )
continue;
end;
i = fix(i - 1);
%
while( true );
i = fix(i + 1);
if( i==j )
break;
end;
t = fix(ix(i+1));
if( ix(i)>t )
k = fix(i);
%
while( true );
ix(k+1) = fix(ix(k));
k = fix(k - 1);
if( t>=ix(k) )
break;
end;
end;
ix(k+1) = fix(t);
end;
end;
%
%     Begin again on another portion of the unsorted array
%
end;
gt(3)=0;
m = fix(m - 1);
if( m==0 )
break;
end;
i = fix(il(m));
j = fix(iu(m));
gt(2)=1;
end;
%
%     Clean up
%
if( kflag<=-1 )
for i = 1 : nn;
ix(i) = fix(-ix(i));
end; i = fix(nn+1);
end;
ix_shape=zeros(ix_shape);ix_shape(:)=ix(1:numel(ix_shape));ix=ix_shape;
iy_shape=zeros(iy_shape);iy_shape(:)=iy(1:numel(iy_shape));iy=iy_shape;
end %subroutine isort
%DECK ISSBCG