spline.src

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

/*
** spline.src     Computes interpolated matrix
**
** (C) Copyright 1996  Aptech Systems, Inc.
** All Rights Reserved.
**
** This Software Product is PROPRIETARY SOURCE CODE OF APTECH
** SYSTEMS, INC.    This File Header must accompany all files using
** any portion, in whole or in part, of this Source Code.   In
** addition, the right to create such files is strictly limited by
** Section 2.A. of the GAUSS Applications License Agreement
** accompanying this Software Product.
**
** If you wish to distribute any portion of the proprietary Source
** Code, in whole or in part, you must first obtain written
** permission from Aptech Systems.
**
**
**  Format                   Purpose                                    Line
** --------------------------------------------------------------------------
**     { u,v,w } = spline2D(x,y,z,sigma,g)   2-D interpolatory spline     27
**       { u,v } = spline1D(x,y,d,s,Sigma,g)  1-D smoothing spline       123
*/


/*
**> spline2D
**
**  Format:  { u,v,w } = spline2D(x,y,z,sigma,g);
**
**  Input:   x       K x 1 vector, x-abscissae (x-axis values).
**
**           y       N x 1 vector, y-abscissae (y-axis values).
**
**           z       K x N matrix, ordinates (z-axis values).
**
**           sigma   scalar, tension factor.
**
**           G       scalar, grid size factor.
**
**
**  Output:  u      1 x K*G vector, x-abscissae, regularly spaced.
**
**           v      N*G x 1 vector, y-abscissae, regularly spaced.
**
**           w      K*G x N*G matrix, interpolated ordinates.
**
**
**  Remarks:
**
**    sigma contains the tension factor. This value indicates
**    the curviness desired. If sigma is nearly zero
**    (e. g. .001) the resulting surface is approximately the
**    tensor product of cubic splines.  If sigma is large
**    (e. g. 50.) the resulting surface is approximately
**    bi-linear. if sigma equals zero tensor products of
**    cubic splines result.  A standard value for sigma is
**    approximately 1.
**
**    G is the grid size factor.  It determines the fineness of
**    the output grid.  For G = 1, the output matrices are
**    identical to the input matrices.  For G = 2, the output
**    grid is twice as fine as the input grid, i.e., u will
**    have twice as many columns as x, and v will have twice as
**    many rows as y.
**
**-------------------**------------------**-------------------**-----------**
**-------------------**------------------**-------------------**-----------*/


declare matrix _spl_i = 1;


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

    local zp0,xx,yy,zz,i,j,vv,xxc,yyr;

    if rows(x) /= rows(z) or rows(y) /= cols(z);
        if not trapchk(1);
            errorlog "arguments not conformable";
            end;
        endif;
        retp(error(0),error(0),error(0));
    endif;

    if rows(x) > rows(y);
        vv = x;
        x = y;
        y = vv;
        z = z';
        vv = 1;
    else;
        vv = 0;
    endif;

    xxc = floor(gridf*rows(x));
    yyr = floor(gridf*rows(y));

    zp0 = _spl_surf1(x,y,z,sigma);

    xx = seqa(x[1],(x[rows(x)]-x[1])./(xxc-1),xxc);
    yy = seqa(y[1],(y[rows(y)]-y[1])./(yyr-1),yyr);

    zz = zeros(xxc,yyr);

    i = 1;
    do until i > xxc;
        j = 1;
        do until j > yyr;
            zz[i,j] = _spl_surf2(xx[i],yy[j],x,y,z,zp0,sigma);
            j = j + 1;
        endo;
        i = i + 1;
    endo;

    if vv;
        retp(yy,xx,zz');
    else;
        retp(xx,yy,zz);
    endif;
endp;







/*
**> spline1D
**
**  Format:  { u,v } = spline1D(x,y,d,s,sigma,g);
**
**  Input:  x      K x 1 vector, x-abscissae (x-axis values).
**
**          y      K x 1 vector, y-ordinates (y-axis values).
**
**          d      K x 1 vector or scalar, observation weights.
**
**          s      scalar, smoothing parameter, if s == 0, curve
**                 performs an interpolation.  If d contains
**                 standard deviation estimates, a reasonable
**                 value for s is K.
**
**          sigma  scalar, tension factor.
**
**          G      scalar, grid size factor.
**
**
**  Output: u      K*G x 1 vector, x-abscissae, regularly spaced.
**
**          v      K*G x 1 vector, interpolated ordinates, regularly spaced.
**
**
**  Remarks:
**
**   sigma contains the tension factor. This value indicates
**   the curviness desired. If sigma is nearly zero
**   (e. g. .001) the resulting curve is approximately the
**   tensor product of cubic splines.  If sigma is large
**   (e. g. 50.) the resulting curve is approximately
**   bi-linear. if sigma equals zero tensor products of
**   cubic splines result.  A standard value for sigma is
**   approximately 1.
**
**   G is the grid size factor.  It determines the fineness of
**   the output grid.  For G = 1, the output matrices are
**   identical to the input matrices.  For G = 2, the output
**   grid is twice as fine as the input grid, i.e., u will
**   have twice as many columns as x, and v will have twice as
**   many rows as y.
*/


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

    local yp0,xx,yy,i,j,xxc,ysp0;

    if rows(x) /= rows(y);
        if not trapchk(1);
            errorlog "arguments not conformable";
            end;
        endif;
        retp(error(0),error(0));
    endif;
    xxc = floor(gridf*rows(x));

    ysp0 = _spl_curv1(x,y,d,s,sigma);

    if scalerr(ysp0);
        if not trapchk(4);
            errorlog "CURVE: spline calculation failed";
            end;
        else;
            retp(error(0),error(0));
        endif;
    endif;

    xx = seqa(x[1],(x[rows(x)]-x[1])./(xxc-1),xxc);
    yy = zeros(xxc,1);

    i = 1;
    do until i > xxc;
        yy[i] = _spl_curv2(xx[i],x,y,ysp0,sigma);
        i = i + 1;
    endo;

    retp(xx,yy);
endp;





proc _spl_surf1 (x,y,z,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 m,n,islpsw,zp,temp,ierr,mm1,m1,mp1,nm1,np1,npm,sigmay,sigmax,
        dely1,dely2, i,j,c1,c2,c3,delyn,delynm,delx1,delx2,delxm,zxy1ns,
        zxymns,jbak,jbakp1,delxmm, del1,deli,diag1,sdiag1,diagi,jm1,jp1,
        npmpj,del2,diag2,diagin,sdiag2,npi,im1,t, ip1,ibak,ibakp1,npibak;

    m = rows(x);
    n = rows(y);
    zp = zeros(3*m*n,1);    /* return array */
    temp = zeros(m+2*n,1);

    mm1 = m - 1;
    mp1 = m + 1;
    nm1 = n - 1;
    np1 = n + 1;
    npm = n + m;
    ierr = 0;

    if n <= 1 or m <= 1;
        ierr = 1;
        retp(error(0));
    endif;

    if y[n] <= y[1];
        ierr = 2;
        retp(error(0));
    endif;

/*
** denormalize tension factor in y-direction
*/

    sigmay = abs(sigma)*(n-1)/(y[n]-y[1]);

/*
** obtain y-partial derivatives along y = y(1)
*/

    dely1 = y[2] - y[1];
    dely2 = dely1 + dely1;

    if n > 2;
        dely2 = y[3] - y[1];
    endif;

    if dely1 <= 0 or dely2 <= dely1;
        ierr = 2;
        retp(error(0));
    endif;

    { c1, c2, c3 } = _spl_ceez(dely1,dely2,sigmay,n);

    i = 1;
    do until i > m;
        zp[_spl_add(i,1,1,m,n)] = c1 * z[i,1] + c2 * z[i,2];
        i = i + 1;
    endo;

    if n /= 2;
        i = 1;
        do until i > m;
            zp[_spl_add(i,1,1,m,n)] = zp[_spl_add(i,1,1,m,n)] + c3 * z[i,3];
            i = i + 1;
        endo;
    endif;

/*
** obtain y-partial derivatives along y = y(n)
*/

    delyn = y[n] - y[nm1];
    delynm = delyn + delyn;

    if n > 2;
        delynm = y[n] - y[n-2];
    endif;

    if delyn <= 0 or delynm <= delyn;
        ierr = 2;
        retp(error(0));
    endif;

    { c1, c2, c3 } = _spl_ceez(-delyn,-delynm,sigmay,n);

    i = 1;
    do until i > m;
        temp[n+i] = c1 * z[i,n] + c2*z[i,nm1];
        i = i + 1;
    endo;

    if n /= 2;
        i = 1;
        do until i > m;
            temp[n+i] = temp[n+i] + c3*z[i,n-2];
            i = i + 1;
        endo;
    endif;

    if x[m] <= x[1];
        ierr = 2;
        retp(error(0));
    endif;

/*
** denormalize tension factor in x-direction
*/

    sigmax = abs(sigma) * (m-1) / (x[m]-x[1]);

/*
** obtain x-partial derivatives along x = x(1)
*/

    delx1 = x[2] - x[1];
    delx2 = delx1 + delx1;

    if m > 2;
        delx2 = x[3] - x[1];
    endif;
    if delx1 <= 0 or delx2 <= delx1;
        ierr = 2;
        retp(error(0));
    endif;

    { c1, c2, c3 } = _spl_ceez(delx1,delx2,sigmax,n);

/*
** obtain x-y-partial derivative at (x(1),y(1))
*/

    zp[2*m*n+1] = c1 * zp[1] + c2 * zp[m+1];

    if m > 2;
        zp[2*m*n+1] = zp[2*m*n+1] + c3 * zp[2*m+1];
    endif;

/*
** obtain x-y-partial derivative at (x(1),y(n))
*/

    zxy1ns = c1 * temp[n+1] + c2 * temp[n+2];

    if m > 2;
        zxy1ns = zxy1ns + c3 * temp[n+3];
    endif;

/*
** obtain x-partial derivative along x = x(m)
*/

    delxm = x[m] - x[mm1];
    delxmm = delxm + delxm;

    if m > 2;
        delxmm = x[m] - x[m-2];
    endif;

    if delxm <= 0 or delxmm < delxm;
        ierr = 2;
        retp(error(0));
    endif;

    { c1, c2, c3 } = _spl_ceez(-delxm,-delxmm,sigmax,n);

/*
** obtain x-y-partial derivative at (x(m),y(1))
*/

    zp[_spl_add(m,1,3,m,n)] = c1 * zp[_spl_add(m,1,1,m,n)]
             + c2 * zp[_spl_add(mm1,1,1,m,n)];
    if m > 2;
        zp[_spl_add(m,1,3,m,n)] = zp[_spl_add(m,1,3,m,n)]
              + c3 * zp[_spl_add(m-2,1,1,m,n)];
    endif;

/*
** obtain x-y-partial derivative at (x(m),y(n))
*/

    zxymns = c1 * temp[npm] + c2 * temp[npm-1];
    if m > 2;
        zxymns = zxymns + c3 * temp[npm-2];
    endif;

/*
** set up right hand sides and tridiagonal system for y-grid
** perform forward elimination
*/

    del1 = y[2] - y[1];

    if del1 <= 0;
        ierr = 2;
        retp(error(0));
    endif;

    deli = 1 / del1;
    i = 1;
    do until i > m;
        zp[_spl_add(i,2,1,m,n)] = deli * (z[i,2] - z[i,1]);
        i = i + 1;
    endo;

    zp[_spl_add(1,2,3,m,n)] = deli * (zp[_spl_add(1,2,2,m,n)]
          - zp[_spl_add(1,1,2,m,n)]);
    zp[_spl_add(m,2,3,m,n)] = deli * (temp[npm+2] - temp[npm+1]);

    { diag1, sdiag1 } = _spl_terms(sigmay,del1);

    diagi = 1 / diag1;
    i = 1;
    do until i > m;
        zp[_spl_add(i,1,1,m,n)] = diagi * (zp[_spl_add(i,2,1,m,n)]
                  - zp[_spl_add(i,1,1,m,n)]);
        i = i + 1;
    endo;

    zp[_spl_add(1,1,3,m,n)] = diagi * (zp[_spl_add(1,2,3,m,n)]
             - zp[_spl_add(1,1,3,m,n)]);
    zp[_spl_add(m,1,3,m,n)] = diagi * (zp[_spl_add(m,2,3,m,n)]
       - zp[_spl_add(m,1,3,m,n)]);
    temp[1] = diagi * sdiag1;

    if n /= 2;
        j = 2;
        do until j > nm1;
            jm1 = j - 1;
            jp1 = j + 1;