home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / listings / v_11_06 / 1106032a < prev    next >
Text File  |  1993-02-17  |  21KB  |  658 lines
  1. /******************************************************
  2.  * RPFT.C -- Curve fitting with optional extrapolation.
  3.  * Fits either polynomial or rational function. Based
  4.  * on algoritms defined in Press, et al., "Numerical
  5.  * Recipes", 2nd ed., Cambridge University Press, 1992
  6.  * Developed using the Borland C compiler.
  7.  *
  8.  *   Lowell Smith
  9.  *   3368 Crestwood Drive
  10.  *   Salt Lake City, UT  84109-3202
  11.  *****************************************************/
  12. #include <stdio.h>
  13. #include <math.h>
  14. #include <string.h>
  15. #include <dir.h>
  16. #include <ctype.h>
  17. #include <stddef.h>
  18. #include <stdlib.h>
  19.  
  20. #define NPFAC 8
  21. #define PI 3.141592653589793
  22. #define PIO2 (PI/2.0)
  23. #define SIGN(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a))
  24. void ratlsq(double xs[],double fs[],int npt,int mm,
  25.       int kk, double cof[], double *dev);
  26. double *dvector(long nl, long nh, const char *ner);
  27. double **dmatrix(long nrl,long nrh,long ncl,long nch,
  28.          const char *ner);
  29. void free_dvector(double *v, long ncl);
  30. void free_dmatrix(double **m, long nrl, long ncl);
  31.  
  32. /*****************************************************/
  33. void main(int argc, char **argv)
  34. { double a,b,*xs,*fs,cof[42],dev=0.,sumd,x,yy;
  35.   int i,j,nn,nd,ncof,dig,npt=0,xln,yln;
  36.   char nar[90];
  37.   char filenam[80];
  38.   void commd(int argc, char**argv, double *a,double *b,
  39.         int *nn, int *nd, int *dig, int *xln,
  40.         int *yln, char *filenam);
  41.   void inpt(int nn,int npt,double a, double b, int xln,
  42.         int yln, double *xs,double *fs,char filenam[]);
  43.   void pcshft(double a,double b,double d[],int n);
  44.   void chebpc(double c[],double d[],int n);
  45.  
  46.   commd(argc,argv,&a,&b,&nn,&nd,&dig,&xln,&yln,
  47.         filenam);
  48.   if (xln) { a=log(a);  b=log(b); }
  49.   if (nd)            /* Fit rational function */
  50.   { ncof=nn+nd+1;   npt=NPFAC*ncof;
  51.     xs=dvector(1L,(long)npt,"xs in main");
  52.     fs=dvector(1L,(long)npt,"fs in main");
  53.     inpt(0, npt, a, b, xln, yln, xs, fs, filenam);
  54.     ratlsq(xs,fs,npt,nn,nd,&cof[0],&dev);
  55.     fprintf(stderr,"Est. max error = %.7G\n",dev);
  56.     free_dvector(xs,1L);
  57.     free_dvector(fs,1L);
  58.   }
  59.   else               /* Fit polynomial */
  60.   { ncof=nn+5; xs=dvector(0L,(long)ncof,"xs in main");
  61.     fs=dvector(0L,(long)ncof,"fs in main");
  62.     inpt(ncof, npt, a, b, xln, yln, xs, fs, filenam);
  63.     chebpc(fs,cof,nn+1);
  64.     pcshft(a,b,cof,nn+1);
  65.     free_dvector(xs,0L); free_dvector(fs,0L);
  66.   }
  67.   for (i=0;i<=nn+nd;++i)
  68.   { sprintf(nar,"%.*G ",dig,cof[i]);
  69.     sscanf(nar,"%lf",&cof[i]);
  70.   }
  71.   if (nd) printf("\n      Rational function coefficien"
  72.      "ts\n   Numerator               Denominator\n\n");
  73.   else printf("\n  Polynomial\n  Coefficients\n\n");
  74.   for (i=0;i<max(nd,nn)+1;++i)
  75.   { if (i<=nn) printf("%-24.16G",cof[i]);
  76.     else printf("%*s",24," ");
  77.     if (i<nd) printf("%-24.16G\n",cof[i+nn+1]);
  78.     else printf("\n");
  79.   }
  80.   for(;;)            /* Calculate trial values */
  81.   { i=-1;
  82.     fprintf(stderr,
  83.           "Enter E to quit or a trial X value: ");
  84.     i=scanf("%lf",&x);
  85.     if (i<1) exit(1); if (xln) x=log(x);
  86.     yy=cof[nn]; for (j=nn-1;j>=0;j--) yy=yy*x+cof[j];
  87.     if (nd)
  88.     { for (sumd=0.,j=nn+nd;j>=nn+1;j--)
  89.                  sumd=(sumd+cof[j])*x;
  90.       yy=yy/(1.0+sumd);
  91.     }
  92.     if (yln) yy=exp(yy); if (xln) x=exp(x);
  93.     fprintf(stderr,"For X=%.8G  Y=%.8G\n",x,yy);
  94.   }
  95. }
  96. /******************************************************
  97.  * COMMD - Parses the command line
  98.  *****************************************************/
  99. void commd(int argc, char**argv, double *a,double *b,
  100.         int *nn, int *nd, int *dig, int *xln,
  101.         int *yln, char *filenam)
  102. { struct ffblk ffblk;
  103.   int i,j,k,l;
  104.   void help(char *msg);
  105.  
  106.   *a=1.11e300; *b=-1.11e300; *nn=-32768; *nd=-32768;
  107.   *dig=7, *xln=0, *yln=0;
  108.   if (argc<2) help(""); filenam[0]=0;
  109.   for (i=1,k=0;i<argc;k=0,++i)
  110.   { for (j=0; j<(l=strlen(argv[i])); ++j)
  111.     { argv[i][j]=toupper(argv[i][j]);
  112.       if (argv[i][j]=='=') ++k;
  113.     }
  114.     if (k)
  115.     { if (k>1 || (!isdigit(argv[i][l-1])
  116.                       && argv[i][l-1] !='.'))
  117.       { fprintf(stderr,"Invalid command line parameter"
  118.                  " %s\n",argv[i]);
  119.         help("");
  120.       }
  121.       if (!strncmp(argv[i],"LL=",3))
  122.       { k=sscanf(&argv[i][3],"%lf",a);
  123.         if (k<1) help("Invalid value for command line"
  124.                    " parameter LL\n");  continue;
  125.       }
  126.       if (!strncmp(argv[i],"UL=",3))
  127.       { k=sscanf(&argv[i][3],"%lf",b);
  128.         if (k<1) help("Invalid value for command line"
  129.                    " parameter UL\n");  continue;
  130.       }
  131.       if (!strncmp(argv[i],"ND=",3))
  132.       { k=sscanf(&argv[i][3],"%d",nd);
  133.         if (k<1) help("Invalid value for command line"
  134.                    " parameter ND\n");  continue;
  135.       }
  136.       if (!strncmp(argv[i],"NN=",3))
  137.       { k=sscanf(&argv[i][3],"%d",nn);
  138.         if (k<1) help("Invalid value for command line"
  139.                    " parameter NN\n");  continue;
  140.       }
  141.       if (!strncmp(&argv[i][0],"-DIG=",5))
  142.       { k=sscanf(&argv[i][6],"%d",dig);
  143.         if (k<1) help("Invalid value for optional"
  144.                    " command line parameter DIG\n");
  145.         continue;
  146.       }
  147.       fprintf(stderr,"Invalid command line parameter "
  148.                    "%s\n",argv[i]);
  149.       help("");
  150.     }
  151.     else
  152.     { if (!strncmp(&argv[i][0],"-XLN",4))
  153.       { *xln=1; continue; }
  154.       if (!strncmp(&argv[i][0],"-YLN",4))
  155.       { *yln=1; continue; }
  156.       if (!filenam[0] && argv[i][0] != '-')
  157.       { strcpy(filenam,argv[i]);
  158.         if (findfirst(filenam,&ffblk,0))
  159.         { fprintf(stderr,"Input file %s doesn't "
  160.             "exist\n",filenam); help("");
  161.         }
  162.       }
  163.       else
  164.       { fprintf(stderr,"Invalid command line parameter"
  165.                    " %s\n",argv[i]);   help("");
  166.       }
  167.     }
  168.   }
  169.   k=0;
  170.   if (*a>1.1e300)
  171.   { fprintf(stderr,"Parameter LL undefined\n"); ++k; }
  172.   if (*b<-1.1e300)
  173.   { fprintf(stderr,"Parameter UL undefined\n"); ++k; }
  174.   if (*nn==-32768)
  175.   { fprintf(stderr,"Parameter NN undefined\n"); ++k;
  176.     *nn=5;
  177.   }
  178.   if (*nd==-32768)
  179.   { fprintf(stderr,"Parameter ND undefined\n"); ++k;
  180.     *nd=5;
  181.   }
  182.   if (filenam[0]==0)
  183.     { fprintf(stderr,"Input file is not defined\n");
  184.       ++k;
  185.     }
  186.   if (*nn<1 || (*nd==0 && *nn>20) || (*nd && *nn>10))
  187.   { fprintf(stderr,"Invalid value %d for "
  188.             "command line parameter NN\n",*nn); ++k;
  189.   }
  190.   if (*nd<0 || *nd>10)
  191.   { fprintf(stderr,"Invalid value %d for "
  192.            "command line parameter ND\n",*nd); ++k;
  193.   }
  194.   if (*a>*b)
  195.   { fprintf(stderr,"Lower limit > upper limit\n");
  196.     ++k;
  197.   }
  198.   if (*dig<1 || *dig>20)
  199.   { fprintf(stderr,"Invalid value %d for "
  200.           "command line option DIG\n",*dig); ++k;
  201.   }
  202.   if (*xln && *a<=0.)
  203.   { fprintf(stderr,"Lower limit is <=0  with "
  204.         "logarithmic transform of X values\n"); ++k;
  205.   }
  206.   if (k) help("");
  207. }
  208. /******************************************************
  209.  * HELP - Prints the help screen
  210.  *****************************************************/
  211. void help(char *msg)
  212. { char *hlp[] =
  213.   { "USAGE: rpft LL=a UL=b NN=n ND=m [-DIG=n -XLN -YLN"
  214.     "] file","","   LL=a   a is the lower limit of the"
  215.     "region for fit","   UL=b   b is the upper limit o"
  216.     "f the region for fit","","   ND=m | m>0: rational"
  217.     " function, denominator degree m,","   NN=n |     "
  218.     " numerator degree n.  1<=m<=10  1<=n<=10","      "
  219.     "  | m=0: Polynomial degree n.  1<=n<=20",""," -DI"
  220.     "G=n   (optional) n = number of significant digits"
  221.     ,"              for coefficients on output. Defaul"
  222.     "t=7.","   -XLN   (optional) Transform X values to"
  223.     " ln(X)","   -YLN   (optional) Transform Y values "
  224.     "to ln(Y)","","   file   File with input X Y data "
  225.     "points, 1 per line","","All command line paramete"
  226.     "rs except DIG, XLN and YLN are ","required. They "
  227.     "may be in any order, upper or lower case."
  228.   };
  229.   int i,n=sizeof(hlp)/sizeof(char *)-1;
  230.  
  231.   fprintf(stderr,"%s\n",msg);
  232.   for (i=0; i<n; ++i) fprintf(stderr,"%s\n",hlp[i]);
  233.   fprintf(stderr,"%s",hlp[n]);  exit(1);
  234. }
  235. /******************************************************
  236.  * INPT - Read input data and calculate data for the
  237.  *        Chebyshev polynomial or rational polynomial
  238.  *        curve fit.
  239.  *****************************************************/
  240. void inpt(int nn, int npt, double a, double b, int xln,
  241.         int yln,double *xs, double *fs, char filenam[])
  242. { double bma,bpa,hth,yy,*xa,*ya,*c,*d;
  243.   int i,j,n=0;
  244.   char nar[82];
  245.   FILE *infile;
  246.   void ratint(double xa[], double ya[], double *c,
  247.         double *d, int n, double x, double *y);
  248.  
  249.   if ((infile=fopen(filenam,"r"))==NULL)
  250.     { fprintf(stderr,"Can't open file\n"); exit(1); }
  251.   while (!feof(infile) && fgets(nar,80,infile))
  252.   { i=sscanf(nar,"%lf%lf",&hth,&yy);
  253.   if (i<1)  break;
  254.     ++n;
  255.     if (i<2)
  256.     { fprintf(stderr,"Need 2 numbers on line %d\n",n);
  257.       exit(1);
  258.     }
  259.     if (xln && hth<=0.)
  260.     { fprintf(stderr,"Nonpositive X = %.8G on line %d "
  261.            "with ln(X) transformation\n",hth,n);
  262.       exit(1);
  263.     }
  264.     if (yln && yy<=0.)
  265.     { fprintf(stderr,"Nonpositive Y = %.8G on line %d "
  266.            "with ln(Y) transformation\n",yy,n);
  267.       exit(1);
  268.     }
  269.   }
  270.   if (n<3)
  271.   { fprintf(stderr,
  272.       "You provided only %d data points\n",n);
  273.     exit(1);
  274.   }
  275.   fseek(infile,0L,0);
  276.   xa=dvector(1L,(long)n,"xa in inpt");
  277.   ya=dvector(1L,(long)n,"ya in inpt");
  278.   c=dvector(1L,(long)n,"c in inpt");
  279.   d=dvector(1L,(long)n,"d in inpt");
  280.   for (i=1; i<=n; ++i)
  281.   { (void)fgets(nar,80,infile);
  282.     sscanf(nar,"%lf%lf",&xa[i],&ya[i]);
  283.     if (xln) xa[i]=log(xa[i]);
  284.     if (yln) ya[i]=log(ya[i]);
  285.   }
  286.   fclose(infile);
  287.   if (!nn)
  288.   /* Calculate arrays of abcissa and function values
  289.    * for rational function evaluation in ratlsq.
  290.    * Return in xs and fs
  291.    */
  292.   { for (i=1;i<=npt;i++)
  293.     { if (i < npt/2)
  294.       { hth=PIO2*(i-1)/(npt-1.0);
  295.         yy=sin(hth); xs[i]=a+(b-a)*yy*yy;
  296.       }
  297.       else
  298.       { hth=PIO2*(npt-i)/(npt-1.0);
  299.         yy=sin(hth); xs[i]=b-(b-a)*yy*yy;
  300.       }
  301.       ratint(xa, ya, c, d, n, xs[i], &yy);
  302.       fs[i]=yy;
  303.     }
  304.   }
  305.   else
  306.   /* Calculate Chebyshev coefficients and return in
  307.    * the array fs
  308.    */
  309.   { bma=0.5*(b-a); bpa=0.5*(b+a);
  310.     for (i=0;i<nn;i++)
  311.        {  ratint(xa,ya,c,d,n,
  312.               cos(PI*(i+0.5)/nn)*bma+bpa,&xs[i]); }
  313.     for (i=0;i<nn;i++)
  314.     { yy=0.; for (j=0;j<nn;j++)
  315.                   yy += xs[j]*cos(PI*i*(j+0.5)/nn);
  316.       fs[i]=2.0*yy/nn;
  317.     }
  318.   }
  319.   free_dvector(xa,1L); free_dvector(ya,1L);
  320.   free_dvector(c,1L); free_dvector(d,1L);  return;
  321. }
  322. /******************************************************
  323.  * RATLSQ - Returns in cof[0..mm+kk] the coefficients
  324.  * of a rational function approximation to the X,Y data
  325.  * points xs[1..npt],fs[1..npt]. The deviation of the
  326.  * approximation (insofar as is known) returned as dev.
  327.  *****************************************************/
  328. void ratlsq(double xs[],double fs[],int npt, int mm,
  329.         int kk, double cof[], double *dev)
  330. { void dsvbksb(double **u, double  w[], double **v,
  331.         int m, int n, double b[], double x[]);
  332.   void dsvdcmp(double **a, int m, int n, double w[],
  333.         double **v);
  334.   int i,it,j,ncof;
  335.   double devmax,e=0.0,power,sum,sumn,sumd,*bb,*coff,
  336.         *ee,**u,**v,*w,*wt;
  337.  
  338.   ncof=mm+kk+1;
  339.   bb=dvector(1L,(long)npt,"bb in ratlsq");
  340.   ee=dvector(1L,(long)npt,"ee in ratlsq");
  341.   coff=dvector(0L,(long)(ncof-1),"coff in ratlsq");
  342.   u=dmatrix(1L,(long)npt,1L,(long)ncof,"u in ratlsq");
  343.   v=dmatrix(1L,(long)ncof,1L,(long)ncof,"v in ratlsq");
  344.   w=dvector(1L,(long)ncof,"w in ratlsq");
  345.   wt=dvector(1L,(long)npt,"wt in ratlsq");
  346.   *dev=1.e50;
  347.   for (i=1;i<=npt;++i) { wt[i]=1.0; ee[i]=1.0; }
  348.   for (it=1;it<=5;it++)
  349.   { for (i=1;i<=npt;i++)
  350.     { power=wt[i]; bb[i]=power*(fs[i]+SIGN(e,ee[i]));
  351.       for (j=1;j<=mm+1;j++)
  352.              { u[i][j]=power; power *= xs[i]; }
  353.       power = -bb[i];
  354.       for (j=mm+2;j<=ncof;j++)
  355.              { power *= xs[i]; u[i][j]=power; }
  356.     }
  357.     dsvdcmp(u,npt,ncof,w,v);
  358.     dsvbksb(u,w,v,npt,ncof,bb,coff-1);
  359.     devmax=sum=0.0;
  360.     for (j=1;j<=npt;j++)
  361.     { e=xs[j]; for (sumn=coff[mm],i=mm-1;i>=0;i--)
  362.                        sumn=sumn*e+coff[i];
  363.       for (sumd=0.0,i=mm+kk;i>=mm+1;i--)
  364.                        sumd=(sumd+coff[i])*e;
  365.       ee[j]=sumn/(1.0+sumd)-fs[j];
  366.       wt[j]=fabs(ee[j]); sum += wt[j];
  367.       if (wt[j] > devmax) devmax=wt[j];
  368.     }
  369.     e=sum/npt;
  370.     if (devmax <= *dev) 
  371.     { for(j=0;j<ncof;j++)cof[j]=coff[j]; *dev=devmax; }
  372.     fprintf(stderr,"For ratlsq iteration %d max error="
  373.          " %.5G\n",it,devmax);
  374.   }
  375.   free_dvector(wt,1L); free_dvector(w,1L);
  376.   free_dmatrix(v,1L,1L); free_dmatrix(u,1L,1L);
  377.   free_dvector(ee,1L); free_dvector(coff,0L);
  378.   free_dvector(bb,1L);
  379. }
  380. /******************************************************
  381.  * DSVDCMP - Computes singular value decomposition of
  382.  * the matrix a[1..m][1..n].
  383.  *****************************************************/
  384. void dsvdcmp(double **a, int m, int n, double w[],
  385.         double **v)
  386. { double dpythag(double a, double b);
  387.   int flag,i,its,j,jj,k,l=0,nm=0;
  388.   double anorm,c,f,g,h,s,scale,x,y,z,*rv1;
  389.  
  390.   rv1=dvector(1L,(long)n,"rv1 in dsvdcmp");
  391.   g=scale=anorm=0.0;
  392.   for (i=1;i<=n;i++)
  393.   { l=i+1; rv1[i]=scale*g; g=s=scale=0.0;
  394.     if (i <= m)
  395.     { for (k=i;k<=m;k++) scale += fabs(a[k][i]);
  396.       if (scale)
  397.       { for (k=i;k<=m;k++)
  398.             { a[k][i] /= scale; s += a[k][i]*a[k][i]; }
  399.         f=a[i][i]; g = -SIGN(sqrt(s),f); h=f*g-s;
  400.         a[i][i]=f-g;
  401.         for (j=l;j<=n;j++)
  402.         { for(s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j];
  403.           f=s/h;
  404.           for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
  405.         }
  406.         for (k=i;k<=m;k++) a[k][i] *= scale;
  407.       }
  408.     }
  409.     w[i]=scale*g; g=s=scale=0.0;
  410.     if (i <= m && i != n)
  411.     { for (k=l;k<=n;k++) scale += fabs(a[i][k]);
  412.       if (scale)
  413.       { for (k=l;k<=n;k++)
  414.            { a[i][k] /= scale;  s += a[i][k]*a[i][k]; }
  415.         f=a[i][l];  g = -SIGN(sqrt(s),f); h=f*g-s;
  416.         a[i][l]=f-g;
  417.         for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
  418.         for (j=l;j<=m;j++)
  419.         { for(s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k];
  420.           for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
  421.         }
  422.         for (k=l;k<=n;k++) a[i][k] *= scale;
  423.       }
  424.     }
  425.     x=fabs(w[i])+fabs(rv1[i]); if (x>anorm) anorm=x;
  426.   }
  427.   for (i=n;i>=1;i--)
  428.   { if (i < n)
  429.     { if (g)
  430.       { for (j=l;j<=n;j++)
  431.               v[j][i]=(a[i][j]/a[i][l])/g;
  432.         for (j=l;j<=n;j++)
  433.         { for(s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j];
  434.           for (k=l;k<=n;k++) v[k][j] += s*v[k][i];
  435.         }
  436.       }
  437.       for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;
  438.     }
  439.     v[i][i]=1.0; g=rv1[i]; l=i;
  440.   }
  441.   for (i=min(m,n);i>=1;i--)
  442.   { l=i+1; g=w[i];
  443.     for (j=l;j<=n;j++) a[i][j]=0.0;
  444.     if (g)
  445.     { g=1.0/g;
  446.       for (j=l;j<=n;j++)
  447.       { for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j];
  448.         f=(s/a[i][i])*g;
  449.         for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
  450.       }
  451.       for (j=i;j<=m;j++) a[j][i] *= g;
  452.     }
  453.     else for (j=i;j<=m;j++) a[j][i]=0.0;
  454.     ++a[i][i];
  455.   }
  456.   for (k=n;k>=1;k--)
  457.   { for (its=1;its<=30;its++)
  458.     { flag=1;
  459.       for (l=k;l>=1;l--)
  460.       { nm=l-1;
  461.         if ((fabs(rv1[l])+anorm) == anorm)
  462.         { flag=0; break; }
  463.         if ((fabs(w[nm])+anorm) == anorm) break;
  464.       }
  465.       if (flag)
  466.       { c=0.0;  s=1.0;
  467.         for (i=l;i<=k;i++)
  468.         { f=s*rv1[i];  rv1[i]=c*rv1[i];
  469.           if ((fabs(f)+anorm) == anorm) break;
  470.           g=w[i]; h=dpythag(f,g); w[i]=h; h=1.0/h;
  471.           c=g*h; s = -f*h;
  472.           for (j=1;i<=m;j++)
  473.           { y=a[j][nm]; z=a[j][i]; a[j][nm]=y*c+z*s;
  474.             a[j][i]=z*c-y*s;
  475.           }
  476.         }
  477.       }
  478.       z=w[k];
  479.       if (l == k)
  480.       { if (z < 0.0)
  481.         { w[k]=-z;
  482.           for (j=1;j<=n;j++) v[j][k] = -v[j][k];
  483.         }
  484.           break;
  485.       }
  486.       if (its == 30)
  487.       { fprintf(stderr,"No convergence in 30 "
  488.                 "dsvdcmp iterations\n");
  489.         exit(1);
  490.       }
  491.       x=w[l]; nm=k-1; y=w[nm]; g=rv1[nm]; h=rv1[k];
  492.       f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
  493.       g=dpythag(f,1.0);
  494.       f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
  495.       c=s=1.0;
  496.       for  (j=l;j<=nm;j++)
  497.       { i=j+1; g=rv1[i]; y=w[i]; h=s*g; g=c*g;
  498.         z=dpythag(f,h); rv1[j]=z; c=f/z; s=h/z;
  499.         f=x*c+g*s; g=g*c-x*s; h=y*s; y *= c;
  500.         for (jj=1;jj<=n;jj++)
  501.         { x=v[jj][j]; z=v[jj][i]; v[jj][j]=x*c+z*s;
  502.           v[jj][i]=z*c-x*s;
  503.         }
  504.         z=dpythag(f,h); w[j]=z;
  505.         if (z) { z=1.0/z; c=f*z; s=h*z; }
  506.         f=c*g+s*y; x=c*y-s*g;
  507.         for (jj=1;jj<=m;jj++)
  508.         { y=a[jj][j]; z=a[jj][i]; a[jj][j]=y*c+z*s;
  509.           a[jj][i]=z*c-y*s;
  510.         }
  511.       }
  512.       rv1[l]=0.0; rv1[k]=f; w[k]=x;
  513.     }
  514.   }
  515.   free_dvector(rv1,1L);
  516. }
  517. /******************************************************
  518.  * DPYTHAG - Computes sqrt(a*a + b*b) without
  519.  * destructive underflow or overflow.
  520.  *****************************************************/
  521. double dpythag(double a, double b)
  522. { double aba,abb,xx,yy;
  523.   aba=fabs(a); abb=fabs(b); xx=abb/aba; yy=aba/abb;
  524.   if (aba > abb) return aba*sqrt(1.0+xx*xx);
  525.   else return (abb == 0.0 ? 0.0 : abb*sqrt(1.0+yy*yy));
  526. }
  527. /******************************************************
  528.  * DSVBKSB -  Solves A.X = B for a vector X, where A is
  529.  * specified by  the arrays u[1..m][1..n], w[1..n],
  530.  * v[1..n][1..n] as returned by dsvdcmp.
  531.  *****************************************************/
  532. void dsvbksb(double **u, double w[], double **v, int m,
  533.            int n, double b[], double x[])
  534. { int jj,j,i;
  535.   double s,*tmp;
  536.  
  537.   tmp=dvector(1L,(long)n,"tmp in dsvbksb");
  538.   for (j=1;j<=n;j++)
  539.   { s=0.0;
  540.     if (w[j])
  541.     { for (i=1;i<=m;i++) s += u[i][j]*b[i]; s /= w[j];}
  542.     tmp[j]=s;
  543.   }
  544.   for (j=1;j<=n;j++)
  545.   { s=0.0;
  546.     for (jj=1;jj<=n;jj++) s += v[j][jj]*tmp[jj];
  547.     x[j]=s;
  548.   }
  549.   free_dvector(tmp,1L);
  550. }
  551. /******************************************************
  552.  * CHEBPC - Transformation of Chebyshev coefficients
  553.  *    generated in INPT to the interval -1 to 1
  554.  *****************************************************/
  555. void chebpc(double c[],double d[],int n)
  556. { int k,j;
  557.   double sv,*dd;
  558.  
  559.   dd=dvector(0L,(long)n-1,"dd in chebpc");
  560.   for (j=0;j<n;j++) d[j]=dd[j]=0.0;
  561.   d[0]=c[n-1];
  562.   for (j=n-2;j>=1;j--)
  563.   { for (k=n-j;k>=1;k--)
  564.     { sv=d[k]; d[k]=2.0*d[k-1]-dd[k]; dd[k]=sv; }
  565.     sv=d[0]; d[0] = -dd[0]+c[j]; dd[0]=sv;
  566.   }
  567.   for (j=n-1;j>=1;j--) d[j]=d[j-1]-dd[j];
  568.   d[0] = -dd[0]+0.5*c[0];
  569.   free_dvector(dd,0L);
  570. }
  571. /******************************************************
  572.  * PCSHFT - Generate the polynomial coefficients
  573.  *   using the Chebyshev coefficients from CHEBPC
  574.  *****************************************************/
  575. void pcshft(double a,double b,double d[],int n)
  576. { int k,j;
  577.   double fac,cnst;
  578.  
  579.   cnst=2.0/(b-a);   fac=cnst;
  580.   for (j=1;j<n;j++) { d[j] *= fac; fac *= cnst; }
  581.   cnst=0.5*(a+b);
  582.   for (j=0;j<=n-2;j++)
  583.          for (k=n-2;k>=j;k--) d[k] -= cnst*d[k+1];
  584. }
  585. /******************************************************
  586.  * RATINT - Diagonal rational function interpolation in
  587.  * the arrays xa[1..n] and ya[1..n].
  588.  *****************************************************/
  589. void ratint(double xa[], double ya[], double *c,
  590.         double *d, int n, double x, double *y)
  591. { int m,i,ns=1;
  592.   double w,t,hh,h,dd;
  593.  
  594.   hh=fabs(x-xa[1]);
  595.   for (i=1;i<=n;i++)
  596.   { h=fabs(x-xa[i]);
  597.     if (h == 0.0)  { *y=ya[i]; return; }
  598.     else  if (h < hh) { ns=i; hh=h; }
  599.     c[i]=ya[i]; d[i]=ya[i]+1.e-50;
  600.   }
  601.   *y=ya[ns--];
  602.   for (m=1;m<n;m++)
  603.   { for (i=1;i<=n-m;i++)
  604.     { w=c[i+1]-d[i] ; h=xa[i+m]-x; t=(xa[i]-x)*d[i]/h;
  605.       dd=t-c[i+1];
  606.       if (dd == 0.0)
  607.       { fprintf(stderr,"Error in routine ratint\n");
  608.         exit(1); /* The function may have a pole  */
  609.       }          /* at the requested value of x.  */
  610.       dd=w/dd; d[i] =c [i+1]*dd; c[i]=t*dd;
  611.     }
  612.     *y += (2*ns < (n-m) ? c[ns+1] : d[ns--]);
  613.   }
  614.   return;
  615. }
  616. /******************************************************
  617.  *  Utility routines based on Numerical Recipes
  618.  *****************************************************/
  619. double *dvector(long nl, long nh, const char *ner)
  620. { double *v;
  621.  
  622.   v=(double *)malloc((size_t) ((nh-nl+1+1)*
  623.              sizeof(double)));
  624.   if (!v)
  625.   { fprintf(stderr,"Allocation failure for vector %s\n"
  626.           ,ner);    exit(1);
  627.   }
  628.   return  v-nl+1;
  629. }
  630. double **dmatrix(long nrl,long nrh,long ncl,long nch,
  631.                const char *ner)
  632. { long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
  633.   double **m;
  634.  
  635.   m=(double **) malloc((size_t)((nrow+1)*
  636.              sizeof(double*)));
  637.   if (m)
  638.   { m += 1;  m -=  nrl;
  639.     m[nrl]=(double *) malloc((size_t)((nrow*ncol+1)*
  640.              sizeof(double)));
  641.   }
  642.   if (m[nrl])
  643.   { m[nrl] += 1;  m[nrl] -= ncl;
  644.     for(i=nrl+1;i<=nrh;i++)   m[i]=m[i-1]+ncol;
  645.   }
  646.   if (!m || !m[nrl])
  647.   { fprintf(stderr,"Allocation failure for matrix %s\n"
  648.             ,ner);   exit(1);
  649.   }
  650.   return m;
  651. }
  652. void free_dvector(double *v, long nl)
  653. { free(v+nl-1); }
  654.  
  655. void free_dmatrix(double **m, long nrl, long ncl)
  656. { free(m[nrl]+ncl-1); free(m+nrl-1); }
  657. /************* End of file RPFT.C *******************/
  658.