```/***
sort two array by the order of the first one
Qsort is implemented following the algorithm explained in Numrec C, pp 332-335
Under a certain size, a part of the array is sorted by the naive method (straight insertion)
***/

#include <stdlib.h>
#include <stdio.h>

#define NSTACK 1000                           /* the number of subarray that could be memorized */
#define SWITCH_METHOD 7                       /* under that size switch to a Insertion method */

#define SWAPf(a,b) {tmpf=(a); (a)=(b); (b)=tmpf;}
#define SWAPi(a,b) {tmpi=(a); (a)=(b); (b)=tmpi;}

void qsort2(float *array1, long *array2, long n){

float *Stack;       /* where remaining values of subarrays are temporarily stocked */
long nStack=0;     /* How many of these values are in the Stack. Is never odd */

long end=n-1,      /* the last value of the array to sort */
beg=0,        /* the first value ---  */
postbeg;      /* the second value ---  */

float val_postbeg;  /* the value of the postbeg position - the one which is used for partion-exchange */

long demi_len;     /* half of the distance between beg and end */

long i,            /* counter from postbeg to right */
j;            /* counter from end to left */

float val1_i,       /* for insertion stock temporarily a value */
val2_i;

float tmpf;          /* used for the SWAPf macro */
long tmpi;          /* used for the SWAPi macro */

Stack = (float *)malloc( (size_t) NSTACK*sizeof(float));
if(! Stack )fprintf(stderr , "qsort2: not enough memory for Stack"), exit(1) ;

while(1){

if( end-beg+1 > SWITCH_METHOD){

demi_len = (end-beg) >> 1 ;
postbeg = beg+1;

SWAPf( array1[beg+demi_len], array1[postbeg] );
SWAPi( array2[beg+demi_len], array2[postbeg] );

if(array1[beg] > array1[postbeg]){            /* rearrange to have  beg <= postbeg <= end */
SWAPf( array1[beg], array1[postbeg] );
SWAPi( array2[beg], array2[postbeg] );
}

if(array1[beg] > array1[end]){
SWAPf( array1[beg], array1[end] );
SWAPi( array2[beg], array2[end] );
}

if(array1[postbeg] > array1[end]){
SWAPf( array1[postbeg], array1[end] );
SWAPi( array2[postbeg], array2[end] );
}

i = postbeg;
j = end;

val_postbeg =  array1[postbeg];

while(1)                                   /* enter the partition exchange process */
{
do i++; while( array1[i] < val_postbeg );
do j--; while( array1[j] > val_postbeg );

if(j<i) break;

SWAPf( array1[i], array1[j] );
SWAPi( array2[i], array2[j] );
}

SWAPf( array1[postbeg] , array1[j] );   /* place the postbeg value into j */
SWAPi( array2[postbeg] , array2[j] );

if(nStack+2 > NSTACK)
fprintf(stderr, "qsort2: not enough Stack... sorry bye\n"),exit(1);

if(end-i+1 >= j-beg){
Stack[nStack++] = i;                /* stock in Stack the largest and go with the smallest */
Stack[nStack++] = end;

end = j-1;
}
else{
Stack[nStack++] = beg;
Stack[nStack++] = j-1;

beg = i;
}

}
else{                                      /* Under a certain size, switch to the straight insertion method */

for(i=beg+1; i <= end ; i++){

val1_i = array1[i];
val2_i = array2[i];

for(j=i-1;j>=beg;j--){

if(array1[j] < val1_i)break;

array1[j+1] = array1[j];
array2[j+1] = array2[j];
}
array1[j+1]=val1_i;
array2[j+1]=val2_i;
}

if(nStack==0)break;          /* this si the end  - exit the process */

end = Stack[--nStack];       /* go for the next round with the stacked parameters */
beg = Stack[--nStack];
}

}

free(Stack);

}
```