rm funtools for update

1 files changed, 0 insertions, 914 deletions

diff --git a/funtools/wcs/poly.c b/funtools/wcs/poly.c
deleted file mode 100644
index 7e88d95..0000000
--- a/funtools/wcs/poly.c
+++ /dev/null
@@ -1,914 +0,0 @@
- /*
- 				poly.c
-
-*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-*
-*	Part of:	A program using Polynomials
-*
-*	Author:		E.BERTIN (IAP)
-*
-*	Contents:	Polynomial fitting
-*
-*	Last modify:	08/03/2005
-*
-*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-*/
-
-#if HAVE_CONFIG_H
-#include	"conf.h"
-#endif
-
-#include	<math.h>
-#include	<stdio.h>
-#include	<stdlib.h>
-#include	<string.h>
-
-#include	"wcslib.h"
-
-
-#define	QCALLOC(ptr, typ, nel) \
-		{if (!(ptr = (typ *)calloc((size_t)(nel),sizeof(typ)))) \
-		  qerror("Not enough memory for ", \
-			#ptr " (" #nel " elements) !");;}
-
-#define	QMALLOC(ptr, typ, nel) \
-		{if (!(ptr = (typ *)malloc((size_t)(nel)*sizeof(typ)))) \
-		  qerror("Not enough memory for ", \
-			#ptr " (" #nel " elements) !");;}
-
-/********************************* qerror ************************************/
-/*
-I hope it will never be used!
-*/
-void	qerror(char *msg1, char *msg2)
-  {
-  fprintf(stderr, "\n> %s%s\n\n",msg1,msg2);
-  exit(-1);
-  }
-
-
-/****** poly_init ************************************************************
-PROTO   polystruct *poly_init(int *group, int ndim, int *degree, int ngroup)
-PURPOSE Allocate and initialize a polynom structure.
-INPUT   1D array containing the group for each parameter,
-        number of dimensions (parameters),
-        1D array with the polynomial degree for each group,
-        number of groups.
-OUTPUT  polystruct pointer.
-NOTES   -.
-AUTHOR  E. Bertin (IAP)
-VERSION 08/03/2003
- ***/
-polystruct	*poly_init(int *group, int ndim, int *degree, int ngroup)
-  {
-   void	qerror(char *msg1, char *msg2);
-   polystruct	*poly;
-   char		str[512];
-   int		nd[POLY_MAXDIM];
-   int		*groupt,
-		d,g,n,num,den;
-
-  QCALLOC(poly, polystruct, 1);
-  if ((poly->ndim=ndim) > POLY_MAXDIM)
-    {
-    sprintf(str, "The dimensionality of the polynom (%d) exceeds the maximum\n"
-		"allowed one (%d)", ndim, POLY_MAXDIM);
-    qerror("*Error*: ", str);
-    }
-
-  if (ndim)
-    QMALLOC(poly->group, int, poly->ndim);
-    for (groupt=poly->group, d=ndim; d--;)
-      *(groupt++) = *(group++)-1;
-
-  poly->ngroup = ngroup;
-  if (ngroup)
-    {
-    group = poly->group;	/* Forget the original *group */
-
-    QMALLOC(poly->degree, int, poly->ngroup);
-
-/*-- Compute the number of context parameters for each group */
-    memset(nd, 0, ngroup*sizeof(int));
-    for (d=0; d<ndim; d++)
-      {
-      if ((g=group[d])>=ngroup)
-        qerror("*Error*: polynomial GROUP out of range", "");
-      nd[g]++;
-      }
-    }
-
-/* Compute the total number of coefficients */
-  poly->ncoeff = 1;
-  for (g=0; g<ngroup; g++)
-    {
-    if ((d=poly->degree[g]=*(degree++))>POLY_MAXDEGREE)
-      {
-      sprintf(str, "The degree of the polynom (%d) exceeds the maximum\n"
-		"allowed one (%d)", poly->degree[g], POLY_MAXDEGREE);
-      qerror("*Error*: ", str);
-      }
-
-/*-- There are (n+d)!/(n!d!) coeffs per group, that is Prod_(i<=d) (n+i)/i */
-    for (num=den=1, n=nd[g]; d; num*=(n+d), den*=d--);
-    poly->ncoeff *= num/den;
-    }
-
-  QMALLOC(poly->basis, double, poly->ncoeff);
-  QCALLOC(poly->coeff, double, poly->ncoeff);
-
-  return poly;
-  }
-
-
-/****** poly_end *************************************************************
-PROTO   void poly_end(polystruct *poly)
-PURPOSE Free a polynom structure and everything it contains.
-INPUT   polystruct pointer.
-OUTPUT  -.
-NOTES   -.
-AUTHOR  E. Bertin (IAP, Leiden observatory & ESO)
-VERSION 09/04/2000
- ***/
-void	poly_end(polystruct *poly)
-  {
-  if (poly)
-    {
-    free(poly->coeff);
-    free(poly->basis);
-    free(poly->degree);
-    free(poly->group);
-    free(poly);
-    }
-  }
-
-
-/****** poly_func ************************************************************
-PROTO   double poly_func(polystruct *poly, double *pos)
-PURPOSE Evaluate a multidimensional polynom.
-INPUT   polystruct pointer,
-        pointer to the 1D array of input vector data.
-OUTPUT  Polynom value.
-NOTES   Values of the basis functions are updated in poly->basis.
-AUTHOR  E. Bertin (IAP)
-VERSION 03/03/2004
- ***/
-double	poly_func(polystruct *poly, double *pos)
-  {
-   double	xpol[POLY_MAXDIM+1];
-   double      	*post, *xpolt, *basis, *coeff, xval;
-   long double	val;
-   int		expo[POLY_MAXDIM+1], gexpo[POLY_MAXDIM+1];
-   int	       	*expot, *degree,*degreet, *group,*groupt, *gexpot,
-			d,g,t, ndim;
-
-/* Prepare the vectors and counters */
-  ndim = poly->ndim;
-  basis = poly->basis;
-  coeff = poly->coeff;
-  group = poly->group;
-  degree = poly->degree;
-  if (ndim)
-    {
-    for (xpolt=xpol, expot=expo, post=pos, d=ndim; --d;)
-      {
-      *(++xpolt) = 1.0;
-      *(++expot) = 0;
-      }
-    for (gexpot=gexpo, degreet=degree, g=poly->ngroup; g--;)
-      *(gexpot++) = *(degreet++);
-    if (gexpo[*group])
-      gexpo[*group]--;
-    }
-
-/* The constant term is handled separately */
-  val = *(coeff++);
-  *(basis++) = 1.0;
-  *expo = 1;
-  *xpol = *pos;
-
-/* Compute the rest of the polynom */
-  for (t=poly->ncoeff; --t; )
-    {
-/*-- xpol[0] contains the current product of the x^n's */
-    val += (*(basis++)=*xpol)**(coeff++);
-/*-- A complex recursion between terms of the polynom speeds up computations */
-/*-- Not too good for roundoff errors (prefer Horner's), but much easier for */
-/*-- multivariate polynomials: this is why we use a long double accumulator */
-    post = pos;
-    groupt = group;
-    expot = expo;
-    xpolt = xpol;
-    for (d=0; d<ndim; d++, groupt++)
-      if (gexpo[*groupt]--)
-        {
-        ++*(expot++);
-        xval = (*(xpolt--) *= *post);
-        while (d--)
-          *(xpolt--) = xval;
-        break;
-        }
-      else
-        {
-        gexpo[*groupt] = *expot;
-        *(expot++) = 0;
-        *(xpolt++) = 1.0;
-        post++;
-        }
-    }
-
-  return (double)val;
-  }
-
-
-/****** poly_fit *************************************************************
-PROTO   double poly_fit(polystruct *poly, double *x, double *y, double *w,
-        int ndata, double *extbasis)
-PURPOSE Least-Square fit of a multidimensional polynom to weighted data.
-INPUT   polystruct pointer,
-        pointer to the (pseudo)2D array of inputs to basis functions,
-        pointer to the 1D array of data values,
-        pointer to the 1D array of data weights,
-        number of data points,
-        pointer to a (pseudo)2D array of computed basis function values.
-OUTPUT  Chi2 of the fit.
-NOTES   If different from NULL, extbasis can be provided to store the
-        values of the basis functions. If x==NULL and extbasis!=NULL, the
-        precomputed basis functions stored in extbasis are used (which saves
-        CPU). If w is NULL, all points are given identical weight.
-AUTHOR  E. Bertin (IAP, Leiden observatory & ESO)
-VERSION 08/03/2005
- ***/
-void	poly_fit(polystruct *poly, double *x, double *y, double *w, int ndata,
-		double *extbasis)
-  {
-   void	qerror(char *msg1, char *msg2);
-   double	/*offset[POLY_MAXDIM],*/x2[POLY_MAXDIM],
-		*alpha,*alphat, *beta,*betat, *basis,*basis1,*basis2, *coeff,
-		*extbasist,*xt,
-		val,wval,yval;
-   int		ncoeff, ndim, matsize,
-		d,i,j,n;
-
-  if (!x && !extbasis)
-    qerror("*Internal Error*: One of x or extbasis should be "
-	"different from NULL\nin ", "poly_func()");
-  ncoeff = poly->ncoeff;
-  ndim = poly->ndim;
-  matsize = ncoeff*ncoeff;
-  basis = poly->basis;
-  extbasist = extbasis;
-  QCALLOC(alpha, double, matsize);
-  QCALLOC(beta, double, ncoeff);
-
-/* Subtract an average offset to maintain precision (droped for now ) */
-/*
-  if (x)
-    {
-    for (d=0; d<ndim; d++)
-      offset[d] = 0.0;
-    xt = x;
-    for (n=ndata; n--;)
-      for (d=0; d<ndim; d++)
-        offset[d] += *(xt++);
-    for (d=0; d<ndim; d++)
-      offset[d] /= (double)ndata;    
-    }
-*/ 
-/* Build the covariance matrix */
-  xt = x;
-  for (n=ndata; n--;)
-    {
-    if (x)
-      {
-/*---- If x!=NULL, compute the basis functions */
-      for (d=0; d<ndim; d++)
-        x2[d] = *(xt++)/* - offset[d]*/;     
-      poly_func(poly, x2);
-/*---- If, in addition, extbasis is provided, then fill it */
-      if (extbasis)
-        for (basis1=basis,j=ncoeff; j--;)
-          *(extbasist++) = *(basis1++);
-      }
-    else
-/*---- If x==NULL, then rely on pre-computed basis functions */
-      for (basis1=basis,j=ncoeff; j--;)
-        *(basis1++) = *(extbasist++);
-
-    basis1 = basis;
-    wval = w? *(w++) : 1.0;
-    yval = *(y++);
-    betat = beta;
-    alphat = alpha;
-    for (j=ncoeff; j--;)
-      {
-      val = *(basis1++)*wval;
-      *(betat++) += val*yval;
-      for (basis2=basis,i=ncoeff; i--;)
-        *(alphat++) += val**(basis2++);
-      }
-    }
-
-/* Solve the system */
-  poly_solve(alpha,beta,ncoeff);
-
-  free(alpha);
-
-/* Now fill the coeff array with the result of the fit */
-  betat = beta;
-  coeff = poly->coeff;
-  for (j=ncoeff; j--;)
-    *(coeff++) = *(betat++);
-/*
-  poly_addcste(poly, offset);
-*/
-  free(beta);
-
-  return;
-  }
-
-
-/****** poly_addcste *********************************************************
-PROTO   void poly_addcste(polystruct *poly, double *cste)
-PURPOSE Modify matrix coefficients to mimick the effect of adding a cst to
-	the input of a polynomial.
-INPUT   Pointer to the polynomial structure,
-        Pointer to the vector of cst.
-OUTPUT  -.
-NOTES   Requires quadruple-precision. **For the time beeing, this function
-	returns completely wrong results!!**
-AUTHOR  E. Bertin (IAP)
-VERSION 03/03/2004
- ***/
-void	poly_addcste(polystruct *poly, double *cste)
-  {
-   long double	*acoeff;
-   double	*coeff,*mcoeff,*mcoefft,
-		val;
-   int		*mpowers,*powers,*powerst,*powerst2,
-		i,j,n,p, denum, flag, maxdegree, ncoeff, ndim;
-
-  ncoeff = poly->ncoeff;
-  ndim = poly->ndim;
-  maxdegree = 0;
-  for (j=0; j<poly->ngroup; j++)
-    if (maxdegree < poly->degree[j])
-      maxdegree = poly->degree[j];
-  maxdegree++;		/* Actually we need maxdegree+1 terms */
-  QCALLOC(acoeff, long double, ncoeff);
-  QCALLOC(mcoeff, double, ndim*maxdegree);
-  QCALLOC(mpowers, int, ndim);
-  mcoefft = mcoeff;		/* To avoid gcc -Wall warnings */
-  powerst = powers = poly_powers(poly);
-  coeff = poly->coeff;
-  for (i=0; i<ncoeff; i++)
-    {
-    for (j=0; j<ndim; j++)
-      {
-      mpowers[j] = n = *(powerst++);
-      mcoefft = mcoeff+j*maxdegree+n;
-      denum = 1;
-      val = 1.0;
-      for (p=n+1; p--;)
-        {
-        *(mcoefft--) = val;
-        val *= (cste[j]*(n--))/(denum++);	/* This is C_n^p X^(n-p) */
-        }
-      }
-/*-- Update all valid coefficients */
-    powerst2 = powers;
-    for (p=0; p<ncoeff; p++)
-      {
-/*---- Check that this combination of powers is included in the series above */
-      flag = 0;
-      for (j=0; j<ndim; j++)
-        if (mpowers[j] < powerst2[j])
-	  {
-          flag = 1;
-          powerst2 += ndim;
-          break;
-          }
-      if (flag == 1)
-        continue;
-      val = 1.0;
-      mcoefft = mcoeff;
-      for (j=ndim; j--; mcoefft += maxdegree)
-        val *= mcoefft[*(powerst2++)];
-      acoeff[i] += val*coeff[p];
-/*
-printf("%g \n", val);
-*/
-      }
-    }
-
-/* Add the new coefficients to the previous ones */
-
-  for (i=0; i<ncoeff; i++)
-{
-/*
-printf("%g %g\n", coeff[i], (double)acoeff[i]);
-*/
-    coeff[i] = (double)acoeff[i];
-}
-
-  free(acoeff);
-  free(mcoeff);
-  free(mpowers);
-  free(powers);
-
-  return;
-  }
-
-/****** poly_solve ************************************************************
-PROTO   void poly_solve(double *a, double *b, int n)
-PURPOSE Solve a system of linear equations, using Cholesky decomposition or
-	SVD (if the former fails due to hidden correlation between variables).
-INPUT   Pointer to the (pseudo 2D) matrix of coefficients,
-        pointer to the 1D column vector,
-        matrix size.
-OUTPUT  -.
-NOTES   -.
-AUTHOR  E. Bertin (IAP, Leiden observatory & ESO)
-VERSION 21/09/2004
- ***/
-void	poly_solve(double *a, double *b, int n)
-  {
-   double	*vmat,*wmat;
-
-  if (cholsolve(a,b,n))
-    {
-    QMALLOC(vmat, double, n*n);
-    QMALLOC(wmat, double, n);
-    svdsolve(a, b, n,n, vmat,wmat);
-    free(vmat);
-    free(wmat);
-    }
-
-  return;
-  }
-
-/****** cholsolve *************************************************************
-PROTO   void cholsolve(double *a, double *b, int n)
-PURPOSE Solve a system of linear equations, using Cholesky decomposition.
-INPUT   Pointer to the (pseudo 2D) matrix of coefficients,
-        pointer to the 1D column vector,
-        matrix size.
-OUTPUT  -1 if the matrix is not positive-definite, 0 otherwise.
-NOTES   Based on Numerical Recipes, 2nd ed. (Chap 2.9). The matrix of
-        coefficients must be symmetric and positive definite.
-AUTHOR  E. Bertin (IAP, Leiden observatory & ESO)
-VERSION	28/10/2003
- ***/
-int	cholsolve(double *a, double *b, int n)
-  {
-   void	qerror(char *msg1, char *msg2);
-   double	*p, *x, sum;
-   int		i,j,k;
-
-/* Allocate memory to store the diagonal elements */
-  QMALLOC(p, double, n);
-
-/* Cholesky decomposition */
-  for (i=0; i<n; i++)
-    for (j=i; j<n; j++)
-      {
-      for (sum=a[i*n+j],k=i-1; k>=0; k--)
-        sum -= a[i*n+k]*a[j*n+k];
-      if (i==j)
-        {
-        if (sum <= 0.0)
-	  {
-          free(p);
-          return -1;
-          }
-        p[i] = sqrt(sum);
-        }
-      else
-        a[j*n+i] = sum/p[i];
-      }
-
-/* Solve the system */
-  x = b;		/* Just to save memory:  the solution replaces b */
-  for (i=0; i<n; i++)
-    {
-    for (sum=b[i],k=i-1; k>=0; k--)
-      sum -= a[i*n+k]*x[k];
-    x[i] = sum/p[i];
-    }
-
-  for (i=n-1; i>=0; i--)
-    {
-    for (sum=x[i],k=i+1; k<n; k++)
-      sum -= a[k*n+i]*x[k];
-    x[i] = sum/p[i];
-    }
-
-  free(p);
-
-  return 0;
-  }
-
-
-/****** svdsolve *************************************************************
-PROTO   void svdsolve(double *a, double *b, int m, int n, double *vmat,
-		double *wmat)
-PURPOSE General least-square fit A.x = b, based on Singular Value
-	Decomposition (SVD).
-	Loosely adapted from Numerical Recipes in C, 2nd Ed. (p. 671).
-INPUT   Pointer to the (pseudo 2D) matrix of coefficients,
-        pointer to the 1D column vector (replaced by solution in output),
-        number of matrix rows,
-	number of matrix columns,
-	pointer to the (pseudo 2D) SVD matrix,
-	pointer to the diagonal (1D) matrix of singular values.	
-OUTPUT  -.
-NOTES   Loosely adapted from Numerical Recipes in C, 2nd Ed. (p. 671). The a
-	and v matrices are transposed with respect to the N.R. convention.
-AUTHOR  E. Bertin (IAP)
-VERSION 26/12/2003
- ***/
-void svdsolve(double *a, double *b, int m, int n, double *vmat, double *wmat)
-  {
-#define MAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
-        (maxarg1) : (maxarg2))
-#define PYTHAG(a,b)     ((at=fabs(a)) > (bt=fabs(b)) ? \
-                                  (ct=bt/at,at*sqrt(1.0+ct*ct)) \
-                                : (bt ? (ct=at/bt,bt*sqrt(1.0+ct*ct)): 0.0))
-#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
-#define TOL             1.0e-11
-   void	qerror(char *msg1, char *msg2);
-
-   int                  flag,i,its,j,jj,k,l,mmi,nm, nml;
-   double               *w,*ap,*ap0,*ap1,*ap10,*rv1p,*vp,*vp0,*vp1,*vp10,
-                        *bp,*tmpp, *rv1,*tmp, *sol,
-			c,f,h,s,x,y,z,
-                        anorm, g, scale,
-                        at,bt,ct,maxarg1,maxarg2,
-                        thresh, wmax;
-
-  anorm = g = scale = 0.0;
-  if (m < n)
-    qerror("*Error*: Not enough rows for solving the system ", "in svdfit()");
-  
-  sol = b;	/* The solution overwrites the input column matrix */
-  QMALLOC(rv1, double, n);
-  QMALLOC(tmp, double, n);
-  l = nm = nml = 0;	/* To avoid gcc -Wall warnings */
-  for (i=0;i<n;i++)
-    {
-    l = i+1;
-    nml = n-l;
-    rv1[i] = scale*g;
-    g = s = scale = 0.0;
-    if ((mmi = m - i) > 0)
-      {
-      ap = ap0 = a+i*(m+1);
-      for (k=mmi;k--;)
-        scale += fabs(*(ap++));
-      if (scale)
-        {
-        for (ap=ap0,k=mmi; k--; ap++)
-          {
-          *ap /= scale;
-          s += *ap**ap;
-          }
-        f = *ap0;
-        g = -SIGN(sqrt(s),f);
-        h = f*g-s;
-        *ap0 = f-g;
-        ap10 = a+l*m+i;
-        for (j=nml;j--; ap10+=m)
-          {
-          for (s=0.0,ap=ap0,ap1=ap10,k=mmi; k--;)
-            s += *(ap1++)**(ap++);
-          f = s/h;
-          for (ap=ap0,ap1=ap10,k=mmi; k--;)
-            *(ap1++) += f**(ap++);
-          }
-        for (ap=ap0,k=mmi; k--;)
-          *(ap++) *= scale;
-        }
-      }
-    wmat[i] = scale*g;
-    g = s = scale = 0.0;
-    if (i < m && i+1 != n)
-      {
-      ap = ap0 = a+i+m*l;
-      for (k=nml;k--; ap+=m)
-        scale += fabs(*ap);
-      if (scale)
-        {
-        for (ap=ap0,k=nml;k--; ap+=m)
-          {
-          *ap /= scale;
-          s += *ap**ap;
-          }
-        f=*ap0;
-        g = -SIGN(sqrt(s),f);
-        h=f*g-s;
-        *ap0=f-g;
-        rv1p = rv1+l;
-        for (ap=ap0,k=nml;k--; ap+=m)
-          *(rv1p++) = *ap/h;
-        ap10 = a+l+m*l;
-        for (j=m-l; j--; ap10++)
-          {
-          for (s=0.0,ap=ap0,ap1=ap10,k=nml; k--; ap+=m,ap1+=m)
-            s += *ap1**ap;
-          rv1p = rv1+l;
-          for (ap1=ap10,k=nml;k--; ap1+=m)
-            *ap1 += s**(rv1p++);
-          }
-        for (ap=ap0,k=nml;k--; ap+=m)
-          *ap *= scale;
-        }
-      }
-    anorm=MAX(anorm,(fabs(wmat[i])+fabs(rv1[i])));
-    }
-
-  for (i=n-1;i>=0;i--)
-    {
-    if (i < n-1)
-      {
-      if (g)
-        {
-        ap0 = a+l*m+i;
-        vp0 = vmat+i*n+l;
-        vp10 = vmat+l*n+l;
-        g *= *ap0;
-        for (ap=ap0,vp=vp0,j=nml; j--; ap+=m)
-          *(vp++) = *ap/g;
-        for (j=nml; j--; vp10+=n)
-          {
-          for (s=0.0,ap=ap0,vp1=vp10,k=nml; k--; ap+=m)
-            s += *ap**(vp1++);
-          for (vp=vp0,vp1=vp10,k=nml; k--;)
-            *(vp1++) += s**(vp++);
-          }
-        }
-      vp = vmat+l*n+i;
-      vp1 = vmat+i*n+l;
-      for (j=nml; j--; vp+=n)
-        *vp = *(vp1++) = 0.0;
-      }
-    vmat[i*n+i]=1.0;
-    g=rv1[i];
-    l=i;
-    nml = n-l;
-    }
-
-  for (i=(m<n?m:n); --i>=0;)
-    {
-    l=i+1;
-    nml = n-l;
-    mmi=m-i;
-    g=wmat[i];
-    ap0 = a+i*m+i;
-    ap10 = ap0 + m;
-    for (ap=ap10,j=nml;j--;ap+=m)
-      *ap=0.0;
-    if (g)
-      {
-      g=1.0/g;
-      for (j=nml;j--; ap10+=m)
-        {
-        for (s=0.0,ap=ap0,ap1=ap10,k=mmi; --k;)
-              s += *(++ap)**(++ap1);
-        f = (s/(*ap0))*g;
-        for (ap=ap0,ap1=ap10,k=mmi;k--;)
-          *(ap1++) += f**(ap++);
-        }
-      for (ap=ap0,j=mmi;j--;)
-        *(ap++) *= g;
-      }
-    else
-      for (ap=ap0,j=mmi;j--;)
-        *(ap++)=0.0;
-    ++(*ap0);
-    }
-
-  for (k=n; --k>=0;)
-      {
-      for (its=0;its<100;its++)
-        {
-        flag=1;
-        for (l=k;l>=0;l--)
-          {
-          nm=l-1;
-          if (fabs(rv1[l])+anorm == anorm)
-            {
-            flag=0;
-            break;
-            }
-          if (fabs(wmat[nm])+anorm == anorm)
-            break;
-          }
-        if (flag)
-          {
-          c=0.0;
-          s=1.0;
-          ap0 = a+nm*m;
-          ap10 = a+l*m;
-          for (i=l; i<=k; i++,ap10+=m)
-            {
-            f=s*rv1[i];
-            if (fabs(f)+anorm == anorm)
-              break;
-            g=wmat[i];
-            h=PYTHAG(f,g);
-            wmat[i]=h;
-            h=1.0/h;
-            c=g*h;
-            s=(-f*h);
-            for (ap=ap0,ap1=ap10,j=m; j--;)
-              {
-              z = *ap1;
-              y = *ap;
-              *(ap++) = y*c+z*s;
-              *(ap1++) = z*c-y*s;
-              }
-            }
-          }
-        z=wmat[k];
-        if (l == k)
-          {
-          if (z < 0.0)
-            {
-            wmat[k] = -z;
-            vp = vmat+k*n;
-            for (j=n; j--; vp++)
-              *vp = (-*vp);
-            }
-          break;
-          }
-        if (its == 99)
-          qerror("*Error*: No convergence in 100 SVD iterations ",
-                "in svdfit()");
-        x=wmat[l];
-        nm=k-1;
-        y=wmat[nm];
-        g=rv1[nm];
-        h=rv1[k];
-        f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
-        g=PYTHAG(f,1.0);
-        f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
-        c=s=1.0;
-        ap10 = a+l*m;
-        vp10 = vmat+l*n;
-        for (j=l;j<=nm;j++,ap10+=m,vp10+=n)
-          {
-          i=j+1;
-          g=rv1[i];
-          y=wmat[i];
-          h=s*g;
-          g=c*g;
-          z=PYTHAG(f,h);
-          rv1[j]=z;
-          c=f/z;
-          s=h/z;
-          f=x*c+g*s;
-          g=g*c-x*s;
-          h=y*s;
-          y=y*c;
-          for (vp=(vp1=vp10)+n,jj=n; jj--;)
-            {
-            z = *vp;
-            x = *vp1;
-            *(vp1++) = x*c+z*s;
-            *(vp++) = z*c-x*s;
-            }
-          z=PYTHAG(f,h);
-          wmat[j]=z;
-          if (z)
-            {
-            z=1.0/z;
-            c=f*z;
-            s=h*z;
-            }
-          f=c*g+s*y;
-          x=c*y-s*g;
-          for (ap=(ap1=ap10)+m,jj=m; jj--;)
-            {
-            z = *ap;
-            y = *ap1;
-            *(ap1++) = y*c+z*s;
-            *(ap++) = z*c-y*s;
-            }
-          }
-        rv1[l]=0.0;
-        rv1[k]=f;
-        wmat[k]=x;
-        }
-      }
-
-  wmax=0.0;
-  w = wmat;
-  for (j=n;j--; w++)
-    if (*w > wmax)
-      wmax=*w;
-  thresh=TOL*wmax;
-  w = wmat;
-  for (j=n;j--; w++)
-    if (*w < thresh)
-      *w = 0.0;
-
-  w = wmat;
-  ap = a;
-  tmpp = tmp;
-  for (j=n; j--; w++)
-    {
-    s=0.0;
-    if (*w)
-      {
-      bp = b;
-      for (i=m; i--;)
-        s += *(ap++)**(bp++);
-      s /= *w;
-      }
-    else
-      ap += m;
-    *(tmpp++) = s;
-    }
-
-  vp0 = vmat;
-  for (j=0; j<n; j++,vp0++)
-    {
-    s=0.0;
-    tmpp = tmp;
-    for (vp=vp0,jj=n; jj--; vp+=n)
-      s += *vp**(tmpp++);
-    sol[j]=s;
-    }
-/* Free temporary arrays */
-  free(tmp);
-  free(rv1);
-
-  return;
-  }
-
-#undef SIGN
-#undef MAX
-#undef PYTHAG
-#undef TOL
-
-/****** poly_powers ***********************************************************
-PROTO   int *poly_powers(polystruct *poly)
-PURPOSE	Return an array of powers of polynom terms
-INPUT   polystruct pointer,
-OUTPUT  Pointer to an array of polynom powers (int *), (ncoeff*ndim numbers).
-NOTES   The returned pointer is mallocated.
-AUTHOR  E. Bertin (IAP)
-VERSION 23/10/2003
- ***/
-int	*poly_powers(polystruct *poly)
-  {
-   int		expo[POLY_MAXDIM+1], gexpo[POLY_MAXDIM+1];
-   int	       	*expot, *degree,*degreet, *group,*groupt, *gexpot,
-		*powers, *powerst,
-		d,g,t, ndim;
-
-/* Prepare the vectors and counters */
-  ndim = poly->ndim;
-  group = poly->group;
-  degree = poly->degree;
-  QMALLOC(powers, int, ndim*poly->ncoeff);
-  if (ndim)
-    {
-    for (expot=expo, d=ndim; --d;)
-      *(++expot) = 0;
-    for (gexpot=gexpo, degreet=degree, g=poly->ngroup; g--;)
-      *(gexpot++) = *(degreet++);
-    if (gexpo[*group])
-      gexpo[*group]--;
-    }
-
-/* The constant term is handled separately */
-  powerst = powers;
-  for (d=0; d<ndim; d++)
-    *(powerst++) = 0;
-  *expo = 1;
-
-/* Compute the rest of the polynom */
-  for (t=poly->ncoeff; --t; )
-    {
-    for (d=0; d<ndim; d++)
-      *(powerst++) = expo[d];
-/*-- A complex recursion between terms of the polynom speeds up computations */
-    groupt = group;
-    expot = expo;
-    for (d=0; d<ndim; d++, groupt++)
-      if (gexpo[*groupt]--)
-        {
-        ++*(expot++);
-        break;
-        }
-      else
-        {
-        gexpo[*groupt] = *expot;
-        *(expot++) = 0;
-        }
-    }
-
-  return powers;
-  }
-