void transformgfs()
{ /* do stereographic projection transformation on gene frequencies to
    get variables that come closer to independent Brownian motions */
  long i, j, k, l, m, n, maxalleles;
  double f, sum;
  double *sumprod, *sqrtp, *pbar;
  phenotype3 *c;

  sumprod = (double *)Malloc(loci*sizeof(double));
  sqrtp = (double *)Malloc(totalleles*sizeof(double));
  pbar = (double *)Malloc(totalleles*sizeof(double));
  for (i = 0; i < totalleles; i++) {  /* get mean gene frequencies */
    pbar[i] = 0.0;
    for (j = 0; j < spp; j++)
      pbar[i] += x[j][i];
    pbar[i] /= spp;
    if (pbar[i] == 0.0)
      sqrtp[i] = 0.0;
    else
      sqrtp[i] = sqrt(pbar[i]);
  }
  for (i = 0; i < spp; i++) {
    for (j = 0; j < loci; j++)    /* for each locus, sum of root(p*x) */
      sumprod[j] = 0.0;
    for (j = 0; j < totalleles; j++)
      if ((pbar[j]*x[i][j]) >= 0.0)
        sumprod[locus[j]-1] += sqrtp[j]*sqrt(x[i][j]);
    for (j = 0; j < totalleles; j++) { /* the projection to tangent plane */
      f = (1.0 + sumprod[locus[j]-1])/2.0;
      if (x[i][j] == 0.0)
        x[i][j] = (2.0/f - 1.0)*sqrtp[j];
      else
        x[i][j] = (1.0/f)*sqrt(x[i][j]) + (1.0/f - 1.0)*sqrtp[j];
    }
  }
  maxalleles = 0;
  for (i = 0; i < loci; i++)
    if (alleles[i] > maxalleles)
      maxalleles = alleles[i];
  c = (phenotype3 *)Malloc(maxalleles*sizeof(phenotype3));
  for (i = 0; i < maxalleles; i++)  /* enough room for any locus's contrasts */
    c[i] = (double *)Malloc(maxalleles*sizeof(double));
  m = 0;        
  for (j = 0; j < loci; j++) {            /* do this for each locus */
    for (k = 0; k < alleles[j]-1; k++) {  /* one fewer than # of alleles */
      c[k][0] = 1.0;
      for (l = 0; l < k; l++) {          /* for contrasts 1 to k make it ... */
        sum = 0.0;
        for (n = 0; n <= l; n++)
          sum += c[k][n]*c[l][n];
        if (fabs(c[l][l+1]) > 0.000000001) /* ... orthogonal to those ones */
          c[k][l+1] = -sum / c[l][l+1];    /* set coeff to make orthogonal */
        else
          c[k][l+1] = 1.0;
      }
      sum = 0.0;
      for (l = 0; l <= k; l++) /* make it orthogonal to vector of sqrtp's */
        sum += c[k][l]*sqrtp[m+l];
      if (sqrtp[m+k+1] > 0.0000000001)
        c[k][k+1] = - sum / sqrtp[m+k+1];  /* ... setting last coeff */
      else {
        for (l = 0; l <= k; l++)
          c[k][l] = 0.0;
        c[k][k+1] = 1.0;
      }
      sum = 0.0;
      for (l = 0; l <= k+1; l++)
        sum += c[k][l]*c[k][l];
      sum = sqrt(sum);
      for (l = 0; l <= k+1; l++)
        if (sum > 0.0000000001)
          c[k][l] /= sum;
    }
    for (i = 0; i < spp; i++) {     /* the orthonormal axes in the plane */
      for (l = 0; l < alleles[j]-1; l++) { /* compute the l-th one */
        c[maxalleles-1][l] = 0.0;  /* temporarily store it ... */
        for (n = 0; n <= l+1; n++) 
          c[maxalleles-1][l] += c[l][n]*x[i][m+n];
      }
      for (l = 0; l < alleles[j]-1; l++)
        x[i][m+l] = c[maxalleles-1][l];   /* replace the gene freqs by it */
    }
    m += alleles[j];
  }
  for (i = 0; i < maxalleles; i++)
    free(c[i]);
  free(c);
  free(sumprod);
  free(sqrtp);
  free(pbar);
} /* transformgfs */