# include "data.h"

void sort (double *value, int N) {

  int l,j,ir,i;
  double temp;
  
  if (N < 2) return;
  l = N >> 1;
  ir = N - 1;
  for (;;) {
    if (l > 0) {
      temp = value[--l];
    }
    else {
      temp = value[ir];
      value[ir] = value[0];
      if (--ir == 0) {
	value[0] = temp;
	return;
      }
    }
    i = l;
    j = (l << 1) + 1;
    while (j <= ir) {
      if (j < ir && value[j] < value[j+1]) ++j;
      if (temp < value[j]) {
	value[i]=value[j];
	j += (i=j) + 1;
      }
      else j = ir + 1;
    }
    value[i] = temp;
  }
}

void fsort (float *value, int N) {

  int l,j,ir,i;
  float temp;
  
  if (N < 2) return;
  l = N >> 1;
  ir = N - 1;
  for (;;) {
    if (l > 0) {
      temp = value[--l];
    }
    else {
      temp = value[ir];
      value[ir] = value[0];
      if (--ir == 0) {
	value[0] = temp;
	return;
      }
    }
    i = l;
    j = (l << 1) + 1;
    while (j <= ir) {
      if (j < ir && value[j] < value[j+1]) ++j;
      if (temp < value[j]) {
	value[i]=value[j];
	j += (i=j) + 1;
      }
      else j = ir + 1;
    }
    value[i] = temp;
  }
}

void sortpair (double *value1, double *value2, int N) {

  int l,j,ir,i;
  double temp1, temp2;
  
  if (N < 2) return;
  l = N >> 1;
  ir = N - 1;
  for (;;) {
    if (l > 0) {
      l--;
      temp1 = value1[l];
      temp2 = value2[l];
    }
    else {
      temp1 = value1[ir];
      temp2 = value2[ir];
      value1[ir] = value1[0];
      value2[ir] = value2[0];
      if (--ir == 0) {
	value1[0] = temp1;
	value2[0] = temp2;
	return;
      }
    }
    i = l;
    j = (l << 1) + 1;
    while (j <= ir) {
      if (j < ir && value1[j] < value1[j+1]) ++j;
      if (temp1 < value1[j]) {
	value1[i]=value1[j];
	value2[i]=value2[j];
	j += (i=j) + 1;
      }
      else j = ir + 1;
    }
    value1[i] = temp1;
    value2[i] = temp2;
  }
}

void fsortpair (float *X, float *Y, int N) {

  int l,j,ir,i;
  float tX, tY;
  
  if (N < 2) return;
  l = N >> 1;
  ir = N - 1;
  for (;;) {
    if (l > 0) {
      l--;
      tX = X[l];
      tY = Y[l];
    }
    else {
      tX = X[ir];
      X[ir] = X[0];
      tY = Y[ir];
      Y[ir] = Y[0];
      if (--ir == 0) {
	X[0] = tX;
	Y[0] = tY;
	return;
      }
    }
    i = l;
    j = (l << 1) + 1;
    while (j <= ir) {
      if (j < ir && X[j] < X[j+1]) j++;
      if (tX < X[j]) {
	X[i] = X[j];
	Y[i] = Y[j];
	j += (i=j) + 1;
      }
      else j = ir + 1;
    }
    X[i] = tX;
    Y[i] = tY;
  }
}

void sortthree (float *X, float *Y, float *Z, int N) {

  int l,j,ir,i;
  float tX, tY, tZ;
  
  if (N < 2) return;
  l = N >> 1;
  ir = N - 1;
  for (;;) {
    if (l > 0) {
      l--;
      tX = X[l];
      tY = Y[l];
      tZ = Z[l];
    }
    else {
      tX = X[ir];
      X[ir] = X[0];
      tY = Y[ir];
      Y[ir] = Y[0];
      tZ = Z[ir];
      Z[ir] = Z[0];
      if (--ir == 0) {
	X[0] = tX;
	Y[0] = tY;
	Z[0] = tZ;
	return;
      }
    }
    i = l;
    j = (l << 1) + 1;
    while (j <= ir) {
      if (j < ir && X[j] < X[j+1]) j++;
      if (tX < X[j]) {
	X[i] = X[j];
	Y[i] = Y[j];
	Z[i] = Z[j];
	j += (i=j) + 1;
      }
      else j = ir + 1;
    }
    X[i] = tX;
    Y[i] = tY;
    Z[i] = tZ;
  }
}

void sort_lists (float *X, float *Y, int *S, int N) {

  int l,j,ir,i;
  double tX, tY, tS;
  
  if (N < 2) return;
  l = N >> 1;
  ir = N - 1;
  for (;;) {
    if (l > 0) {
      l--;
      tX = X[l];
      tY = Y[l];
      tS = S[l];
    }
    else {
      tX = X[ir];
      X[ir] = X[0];
      tY = Y[ir];
      Y[ir] = Y[0];
      tS = S[ir];
      S[ir] = S[0];
      if (--ir == 0) {
	X[0] = tX;
	Y[0] = tY;
	S[0] = tS;
	return;
      }
    }
    i = l;
    j = (l << 1) + 1;
    while (j <= ir) {
      if (j < ir && X[j] < X[j+1]) j++;
      if (tX < X[j]) {
	X[i] = X[j];
	Y[i] = Y[j];
	S[i] = S[j];
	j += (i=j) + 1;
      }
      else j = ir + 1;
    }
    X[i] = tX;
    Y[i] = tY;
    S[i] = tS;
  }
}

void dsort_lists (double *X, double *Y, int *S, int N) {

  int l,j,ir,i;
  double tX, tY, tS;
  
  if (N < 2) return;
  l = N >> 1;
  ir = N - 1;
  for (;;) {
    if (l > 0) {
      l--;
      tX = X[l];
      tY = Y[l];
      tS = S[l];
    }
    else {
      tX = X[ir];
      X[ir] = X[0];
      tY = Y[ir];
      Y[ir] = Y[0];
      tS = S[ir];
      S[ir] = S[0];
      if (--ir == 0) {
	X[0] = tX;
	Y[0] = tY;
	S[0] = tS;
	return;
      }
    }
    i = l;
    j = (l << 1) + 1;
    while (j <= ir) {
      if (j < ir && X[j] < X[j+1]) j++;
      if (tX < X[j]) {
	X[i] = X[j];
	Y[i] = Y[j];
	S[i] = S[j];
	j += (i=j) + 1;
      }
      else j = ir + 1;
    }
    X[i] = tX;
    Y[i] = tY;
    S[i] = tS;
  }
}
