spline.src

来自「没有说明」· SRC 代码 · 共 1,362 行 · 第 1/3 页

                if isw == 0;
                    sinhm = (expx - 1 / expx) / (ax + ax) - 1;
                endif;
            else;
                xs = ax * ax;
                coshm = xs * ((((cp4 * xs + cp3) * xs + cp2) * xs + cp1) *
                    xs + cp0);
                if isw == 0;
                    sinhm = xs * (((sp13 * xs + sp12) * xs + sp11) * xs +
                        sp10);
                endif;
            endif;
        endif;
    else;
        if ax > 4.45;
            if ax > 7.65;
                if ax > 10.1;
                    sinhm = exp(ax) / (ax + ax) - 1;
                else;
                    xs = ax * ax;
                    sinhm = xs * (((sp43 * xs + sp42) * xs + sp41) * xs +
                        1) / ((sq42 * xs + sq41) * xs + sq40);
                endif;
            else;
                xs = ax*ax;
                sinhm = xs * (((sp33 * xs + sp32) * xs + sp31) * xs + 1) /
                    ((sq32 * xs + sq31) * xs + sq30);
            endif;
        else;
            xs = ax * ax;
            if ax > 2.3;
                sinhm = xs * ((((sp24 * xs + sp23) * xs + sp22) * xs +
                    sp21) * xs + sp20);
            else;
                sinhm = xs * (((sp13 * xs + sp12) * xs + sp11) * xs + sp10);
            endif;
        endif;
    endif;
    retp(sinhm,coshm);
endp;

proc _spl_intrvl(t,x);
/*
**    translated into GAUSS from Fortran by Ronald Schoenberg
**
**                    coded into Fortran by alan kaylor cline
**                           from fitpack -- january 26, 1987
**                        a curve and surface fitting package
**                      a product of pleasant valley software
**                  8603 altus cove, austin, texas 78759, usa
*/
    local n,il,ih,tt;
    n = rows(x);

    tt = t;
/*
** check for illegal i
*/
    if _spl_i >= n;
        _spl_i = n / 2;
    endif;

/*
** check old interval and extremes
*/
    if tt < x[_spl_i];
        if tt <= x[2];
            _spl_i = 1;
            retp(1);
        else;
            il = 2;
            ih = _spl_i;
        endif;
    elseif tt <= x[_spl_i+1];
        retp(_spl_i);
    elseif tt >= x[n-1];
        _spl_i = n - 1;
        retp(n - 1);
    else;
        il = _spl_i + 1;
        ih = n - 1;
    endif;

/*
** binary search loop
*/

A1:

    _spl_i = (il + ih) / 2;
    if tt < x[_spl_i];
        ih = _spl_i;
    elseif tt > x[_spl_i+1];
        il = _spl_i + 1;
    else;
        retp(_spl_i);
    endif;
    goto A1;
endp;

proc _spl_add(i,j,k,m,n);
    retp(m*n*(k-1)+m*(i-1)+j);
endp;

proc _spl_hermz(f1,f2,fp1,fp2,del1,del2,dels);
/* inline one dimensional cubic spline interpolation */
    retp((f2 * del1 + f1 * del2) / dels - del1 * del2 * (fp2 * (del1 +
        dels) + fp1 * (del2 + dels)) / (6.*dels) );
endp;

proc _spl_hermnz(f1,f2,fp1,fp2,sigmap,del1,del2,dels,sinhm1,sinhm2,sinhms);
/* inline one dimensional spline under tension interpolation */
    retp( (f2 * del1 + f1 * del2) / dels + (fp2 * del1 * (sinhm1 - sinhms)
        + fp1 * del2 * (sinhm2 - sinhms)) / (sigmap * sigmap * dels * (1 +
        sinhms)) );
endp;






proc _spl_curv2 (t,x,y,yp,sigma);

/*
**    translated into GAUSS from Fortran by Ronald Schoenberg
**
**                    coded into Fortran by alan kaylor cline
**                           from fitpack -- january 26, 1987
**                        a curve and surface fitting package
**                      a product of pleasant valley software
**                  8603 altus cove, austin, texas 78759, usa
**
*/

      local im1,i,n,sigmap,del1,del2,dels,sum,dummy,sigdel,ss,s1,s2;

      n = rows(x);
/*
** determine interval
*/
      im1 = _spl_intrvl(t,x);
      i = im1 + 1;

/*
** denormalize tension factor
*/

      sigmap = abs(sigma) * (n - 1) / (x[n] - x[1]);

/*
** set up and perform interpolation
*/
      del1 = t - x[im1];
      del2 = x[i] - t;
      dels = x[i] - x[im1];
      sum = (y[i] * del1 + y[im1] * del2) / dels;
      if (sigmap == 0);
        retp(sum-del1*del2*(yp(i)*
            (del1+dels)+yp(im1)*(del2+dels))/(6*dels));
      endif;
      sigdel = sigmap * dels;
      { ss,dummy } = _spl_snhcsh(sigdel,-1);
      { s1,dummy } = _spl_snhcsh(sigmap*del1,-1);
      { s2,dummy } = _spl_snhcsh(sigmap*del2,-1);

      retp(sum+(yp[i]*del1*(s1-ss)+yp[im1]*del2*(s2-ss))/
                 (sigdel*sigmap*(1+ss)));
endp;




proc _spl_curv1(x,y,d,s,sigma);

/*
**    translated into GAUSS from Fortran by Ronald Schoenberg
**
**                    coded into Fortran by alan kaylor cline
**                           from fitpack -- january 26, 1987
**                        a curve and surface fitting package
**                      a product of pleasant valley software
**                  8603 altus cove, austin, texas 78759, usa
**
*/
      local n,ierr,eps,ys,ysp,td,tsd1,hd,hsd1,hsd2,rd,rsd1,rsd2,v;
      local f,wim1,wim2,h,p,g,step,sum,sl,su,delyi1,delyi,nm1,ibak,
         yspim2,rsd2i,alphap,betapp,hsd1p,hdim1,hdi,delxi,rsd1i,tui,wi,
         delxi1,dim1,di,rdim1,sigmap,nm3,i,alpha,beta,betap,tt;


      n = rows(x);


      if n < 2;
          retp(error(1));
      endif;

      s = abs(s);
      eps = sqrt(2/n);

      v = zeros(n,1);
      ysp = zeros(n,1);
      ierr = 0;
      p = 0;

      if n == 2;
          retp(ysp);
      endif;

      rsd1 = zeros(n,1);
      rd = zeros(n,1);
      rsd2 = zeros(n,1);
      ys = zeros(n,1);
      tsd1 = zeros(n,1);
      hsd1 = zeros(n,1);
      hsd2 = zeros(n,1);
      td = zeros(n,1);
      hd = zeros(n,1);

      rdim1 = 0;
      yspim2 = 0;

/*
** denormalize tension factor
*/

      sigmap = abs(sigma) * (n - 1) / (x[n] - x[1]);

/*
** form t matrix and second differences of y into ys
*/
      nm1 = n - 1;
      nm3 = n - 3;
      delxi1 = 1;
      delyi1 = 0;
      dim1 = 0;

      i = 1;
      do until i > nm1;
        delxi = x[i+1] - x[i];
        if delxi <= 0;
            retp(error(5));
        endif;
        delyi = (y[i+1]-y[i])/delxi;
        ys[i] = delyi-delyi1;
        { di,tt } = _spl_terms(sigmap,delxi);
        tsd1[i+1] = tt;
        td[i] = di + dim1;
        hd[i] = -(1/delxi+1/delxi1);
        hsd1[i+1] = 1/delxi;
        delxi1 = delxi;
        delyi1 = delyi;
        dim1 = di;
        i = i + 1;
     endo;

/*
** calculate lower and upper tolerances
*/

      sl = s * (1 - eps);
      su = s * (1 + eps);

      d = vec(d);
      if rows(d) == 1;
/*
** form h matrix - d constant
*/
          if d[1] <= 0;
              retp(error(5));
          endif;
          sl = d[1]*d[1]*sl;
          su = d[1]*d[1]*su;
          hsd1p = 0;
          hdim1 = 0;
          i = 2;
          do until i > nm1;
            hdi = hd[i];
            hd[i] = hsd1[i]*hsd1[i]+hdi*hdi+hsd1[i+1]*hsd1[i+1];
            hsd2[i] = hsd1[i]*hsd1p;
            hsd1p = hsd1[i];
            hsd1[i] = hsd1p*(hdi+hdim1);
            hdim1 = hdi;
            i = i + 1;
          endo;
     else;

/*
** form h matrix - d array
*/
          if d[1] <= 0 or d[2] <= 0;
              retp(error(5));
          endif;

          betapp = 0;
          betap = 0;
          alphap = 0;

          i = 2;
          do until i > nm1;
            alpha = hd[i]*d[i]*d[i];
            if d[i+1] <= 0;
                retp(error(5));
            endif;
            beta = hsd1[i+1]*d[i+1]*d[i+1];
            hd[i] = (hsd1[i] * d[i-1])^2 + alpha * hd[i] + beta * hsd1[i+1];
            hsd2[i] = hsd1[i]*betapp;
            hsd1[i] = hsd1[i]*(alpha+alphap);
            alphap = alpha;
            betapp = betap;
            betap = beta;
            i = i + 1;
          endo;

      endif;
/*
** top of iteration
** cholesky factorization of p*t+h into r
*/

A0:

      i = 2;
      do until i > nm1;
        rsd2i = hsd2[i];
        rsd1i = p*tsd1[i]+hsd1[i]-rsd2i*rsd1[i-1];
        rsd2[i] = rsd2i*rdim1;
        rdim1 = rd[i-1];
        rsd1[i] = rsd1i*rdim1;
        rd[i] = 1/(p*td[i]+hd[i]-rsd1i*rsd1[i]-rsd2i*rsd2[i]);
        ysp[i] = ys[i]-rsd1[i]*ysp[i-1]-rsd2[i]*yspim2;
        yspim2 = ysp[i-1];
        i = i + 1;
      endo;

/*
** back solve of r(transpose)* r * ysp = ys
*/

      ysp[nm1] = rd[nm1]*ysp[nm1];
      if n /= 3;
          ibak = 1;
          do until ibak > nm3;
              i = nm1 - ibak;
              ysp[i] = rd[i]*ysp[i]-rsd1[i+1]*ysp[i+1]-rsd2[i+2]*ysp[i+2];
              ibak = ibak + 1;
          endo;
      endif;
      sum = 0;
      delyi1 = 0;
      if rows(d) == 1;

/*
** calculation of residual norm
**  - d constant
*/
           i = 1;
           do until i > nm1;
                delyi = (ysp[i+1]-ysp[i])/(x[i+1]-x[i]);
                v[i] = delyi-delyi1;
                sum = sum+v[i]*(delyi-delyi1);
                delyi1 = delyi;
                i = i + 1;
          endo;
          v[n] = -delyi1;

      else;
/*
** calculation of residual norm
**  - d array
*/
          i = 1;
          do until i > nm1;
            delyi = (ysp[i+1]-ysp[i])/(x[i+1]-x[i]);
            v[i] = (delyi-delyi1)*d[i]*d[i];
            sum = sum+v[i]*(delyi-delyi1);
            delyi1 = delyi;
            i = i + 1;
          endo;
          v[n] = -delyi1*d[n]*d[n];
      endif;
      sum = sum-v[n]*delyi1;

/*
** test for convergence
*/

      if sum <= su;
          goto A1;
      endif;

/*
** calculation of newton correction
*/
      f = 0;
      g = 0;
      wim2 = 0;
      wim1 = 0;

      i = 2;
      do until i > nm1;
        tui = tsd1[i]*ysp[i-1]+td[i]*ysp[i]+tsd1[i+1]*ysp[i+1];
        wi = tui-rsd1[i]*wim1-rsd2[i]*wim2;
        f = f + tui * ysp[i];
        g = g + wi * wi * rd[i];
        wim2 = wim1;
        wim1 = wi;
        i = i + 1;
      endo;

      h = f - p * g;
      if h <= 0;
          goto A1;
      endif;

/*
** update p - newton step
*/
      step = (sum - sqrt(sum * sl)) / h;
      if sl /= 0;
           step = step * sqrt(sum / sl);
      endif;
      p = p + step;

      goto A0;

A1:

/*
** store smoothed y-values and second derivatives
*/
      retp(p*ysp);
endp;











proc(3) = spline(x,y,z,sigma,gridf);
    retp(spline2D(x,y,z,sigma,gridf));
endp;


proc(2) = curve(x,y,d,s,sigma,gridf);
    retp(spline1D(x,y,d,s,sigma,gridf));
endp;