home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / listings / v_08_01 / 8n01054a < prev    next >
Text File  |  1990-02-19  |  2KB  |  77 lines
  1.  
  2. *****Coeff_Cubic (Listing 1)*****
  3.  
  4.  
  5. void coeff_cubic (a,p,q,x,y,lambda,k)
  6. /*
  7.  *  Set up the equations for the Hermite cubic approximation.
  8.  */
  9. double *a,*x,*y,*lambda;
  10. int p,q,k;
  11. {
  12. double d,alpha,beta,d3,alpha3;
  13. int i,j,col;
  14. for (i=0; i<p; i++)
  15.   {
  16.     j = interval (lambda,x[i],k);
  17.     col = SUB(i,2*(j-1),q);
  18.     d = lambda[j] - lambda[j-1];
  19.     alpha = x[i]-lambda[j-1];
  20.     beta = d-alpha;
  21.     d3 = d*d*d;
  22.     alpha3 = alpha*alpha*alpha;
  23.     *(a+col) = (d3-3.0*d*alpha*alpha+2.0*alpha3)/d3;
  24.     *(a+col+1) = d*alpha*beta*beta/d3;
  25.     *(a+col+2) = (3.0*d-2.0*alpha)*alpha*alpha/d3;
  26.     *(a+col+3) = -d*alpha*alpha*beta/d3;
  27.   }
  28. }
  29.  
  30. int interval (x,v,n)
  31. /*
  32.  *  Given a value v find the interval j such that v is in the interval
  33.  *  x[j-1] to x[j], where x is an increasing set of n values.
  34.  */
  35.  double x[],v;
  36.  int n;
  37.   {
  38.   int j = 0, found = 0;
  39.       if (v == x[0]) return(1);
  40.   while (!found && ++j<n)
  41.               found =( v<=x[j] && v> x[j-1]) ? 1:0;
  42.     return(j);
  43.    } 
  44.  
  45.  
  46.  
  47.  
  48. double app_cubic (x,j,lambda,res)
  49. /*
  50.  *  Given the result res[] from the routine Hermite() find the value
  51.  *  of y for the given x value.
  52.  */
  53. double x,*lambda,*res;
  54. int j;
  55. {
  56. double d,alpha,beta,d3,alpha3,sum,val[4];
  57. int i,col;
  58.     col = 2*(j-1);
  59.     d = lambda[j] - lambda[j-1];
  60.     alpha = x-lambda[j-1];
  61.     beta = d-alpha;
  62.     d3 = d*d*d;
  63.     alpha3 = alpha*alpha*alpha;
  64.     val[0] = (d3-3.0*d*alpha*alpha+2.0*alpha3)/d3;
  65.     val[1] = d*alpha*beta*beta/d3;
  66.     val[2] = (3.0*d-2.0*alpha)*alpha*alpha/d3;
  67.     val[3] = -d*alpha*alpha*beta/d3;
  68.     for (sum=0.0,i=0; i<4; i++) sum += val[i]*res[col+i];
  69. return (sum);
  70. }
  71.  
  72.   
  73.  
  74.         
  75.  
  76.  
  77.