void order(double *p, double *q, int *i, int *j)
{
  int temp1;
  double temp2;
  
  if (*p > *q) {
     temp1=*i; *i=*j; *j=temp1;
     temp2=*p; *p=*q; *q=temp2;
  }
}

void bubble(int n, double *a, double *c, int *b)
/* 
   sort a in ascending order using bubblesort algorithm, 
   put the resulting vector in c
   and put the permutation of indexes to b
   n is the dimension of a,b and c
*/
{ 
  int i, j;
  for (i=1; i<=n; i++)
     c[i-1]=a[i-1];

  for (i=0; i<n-1; ++i)
     for (j=n-1; i<j; --j)
        order(&c[j-1], &c[j], &b[j-1],&b[j]);
}