home *** CD-ROM | disk | FTP | other *** search
/ Nebula / nebula.bin / SourceCode / Classes / ContourView / splin2.c < prev    next >
C/C++ Source or Header  |  1992-12-03  |  3KB  |  158 lines
  1. /* 2-D bicubic spline */
  2.  
  3. #include <stdio.h>
  4. #include <stdlib.h>
  5. #include "splin2.h"
  6.  
  7. #define APPNAME        "ContourView"
  8.  
  9. /* Given x1a, x2a, ya (Data), and y2a (computed by ssplie2),
  10.  * this routine computes interpolated value y at point x1, x2 by
  11.  * bicubic spline.  From Press et al., Num Rec for C.
  12.  */
  13. void splin2(float *x1a, float *x2a, float **ya, float **y2a,
  14.         int m, int n, float x1, float x2, float *y)
  15. {
  16. int j;
  17. float *ytmp,*yytmp;
  18.  
  19.     ytmp=fvector(1,n);
  20.     yytmp=fvector(1,n);
  21.     for (j=1;j<=m;j++)
  22.         splint(x2a,ya[j],y2a[j],n,x2,&yytmp[j]);
  23.     spline(x1a,yytmp,m,1.0e30,1.0e30,ytmp);
  24.     splint(x1a,yytmp,ytmp,m,x1,y);
  25.     free_fvector(yytmp,1,n);
  26.     free_fvector(ytmp,1,n);
  27. }
  28.  
  29.  
  30. /* From Press, et al.  Num Rec for C.
  31.  *    Returns second derivatives in array y2a[1..m][1..n]
  32.  */
  33. void splie2(float *x1a, float *x2a, float **ya, int m, int n, float **y2a)
  34. {
  35. int j;
  36.     for (j=1;j<=m;j++)
  37.         spline(x2a,ya[j],n,1.0e30,1.0e30,y2a[j]);
  38. }
  39.  
  40.  
  41.  
  42.  
  43. void splint(float *xa, float *ya, float *y2a, int n, float x, float *y)
  44. {
  45. int klo,khi,k;
  46. float h,b,a;
  47.  
  48.     klo=1;
  49.     khi=n;
  50.     while (khi-klo > 1) {
  51.         k=(khi+klo) >> 1;
  52.         if (xa[k] > x) khi=k;
  53.         else klo=k;
  54.     }
  55.     h=xa[khi]-xa[klo];
  56.     if (h == 0.0)
  57.         fprintf(stderr,"%s: Bad XA input to routine splint()\n", APPNAME);
  58.     a=(xa[khi]-x)/h;
  59.     b=(x-xa[klo])/h;
  60.     *y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
  61. }
  62.  
  63.  
  64.  
  65. void spline(float *x, float *y, int n, float yp1, float ypn, float *y2)
  66. {
  67. int i,k;
  68. float p,qn,sig,un,*u;
  69.  
  70.     u=fvector(1,n-1);
  71.     if (yp1 > 0.99e30)
  72.         y2[1]=u[1]=0.0;
  73.     else {
  74.         y2[1] = -0.5;
  75.         u[1]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);
  76.     }
  77.     for (i=2;i<=n-1;i++) {
  78.         sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
  79.         p=sig*y2[i-1]+2.0;
  80.         y2[i]=(sig-1.0)/p;
  81.         u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
  82.         u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
  83.     }
  84.     if (ypn > 0.99e30)
  85.         qn=un=0.0;
  86.     else {
  87.         qn=0.5;
  88.         un=(3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]));
  89.     }
  90.     y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
  91.     for (k=n-1;k>=1;k--)
  92.         y2[k]=y2[k]*y2[k+1]+u[k];
  93.     free_fvector(u,1,n-1);
  94. }
  95.  
  96.  
  97. /* memory allocation routines */
  98.  
  99. /* ###### 1-d float */
  100.  
  101. float *fvector(int nl, int nh)
  102. {
  103.     float *v;
  104.     v=(float *)malloc((unsigned)(nh-nl+1)*sizeof(float));
  105.     if(!v) fprintf(stderr, "allocation failure in fvector()\n");
  106.     return(v-nl);
  107. }
  108.  
  109.  
  110. void free_fvector(float *v, int nl, int nh)
  111. {
  112.     free((void *)(v+nl));
  113. }
  114.  
  115.  
  116. /* ##### 2-d float */
  117.  
  118. float **fmatrix(int nrl, int nrh, int ncl, int nch)
  119. {
  120.    int i=2,j;
  121.    int error=0;
  122.    float **m;
  123.     m=(float **) malloc((unsigned) (nrh-nrl+1)*sizeof(float *));
  124.     if (!m) fprintf(stderr, "allocation failure 1 in fmatrix()\n");
  125.     else{
  126.        m -= nrl;
  127.        i=nrl;
  128.        do{
  129.         m[i]=(float *) malloc((unsigned) (nch-ncl+1)*sizeof(float));
  130.         if (!m[i]) {
  131.             fprintf(stderr, "allocation failure 2 in fmatrix()\n");
  132.             error=1;
  133.         }else
  134.             m[i] -= ncl;
  135.         i++;
  136.            }while( (i<=nrh) && (!error));
  137.     }
  138.     if(error){
  139.         fprintf(stderr, "Freeing allocated submatrices...\n");
  140.         for(j=i-2;j>=nrl;j--)           /* free those allocated */
  141.             free((void *)(m[j]+ncl));
  142.         free((void *)(m+nrl));
  143.         m=NULL;
  144.      }
  145.     return m;
  146. }
  147.  
  148.  
  149. void free_fmatrix(float **m, int nrl, int nrh, int ncl, int nch)
  150. {
  151.     int i;
  152.     for(i=nrh;i>=nrl;i--)free((void *)(m[i]+ncl));
  153.     free((void *) (m+nrl));
  154.  
  155. }
  156.  
  157.  
  158.