spline.src

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

            npmpj = npm + j;
            del2 = y[jp1] - y[j];

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

            deli = 1 / del2;
            i = 1;
            do until i > m;
                zp[_spl_add(i,jp1,1,m,n)] = deli * (z[i,jp1] - z[i,j]);
                i = i + 1;
            endo;
            zp[_spl_add(1,jp1,3,m,n)] = deli * (zp[_spl_add(1,jp1,2,m,n)]
                - zp[_spl_add(1,
                j,2,m,n)]);
            zp[_spl_add(m,jp1,3,m,n)] = deli * (temp[npmpj+1] - temp[npmpj]);

            { diag2, sdiag2 } = _spl_terms(sigmay,del2);

            diagin = 1./(diag1+diag2-sdiag1*temp[jm1]);

            i = 1;
            do until i > m;
                zp[_spl_add(i,j,1,m,n)] = diagin*(zp[_spl_add(i,jp1,1,m,n)]
                     -zp[_spl_add(i,
                    j,1,m,n)]-sdiag1*zp[_spl_add(i,jm1,1,m,n)]);
                i = i + 1;
            endo;

            zp[_spl_add(1,j,3,m,n)] = diagin * (zp[_spl_add(1,jp1,3,m,n)]
                      - zp[_spl_add(1,
                j,3,m,n)] - sdiag1*zp[_spl_add(1,jm1,3,m,n)]);
            zp[_spl_add(m,j,3,m,n)] = diagin * (zp[_spl_add(m,jp1,3,m,n)]
                   - zp[_spl_add(m,
                j,3,m,n)] - sdiag1*zp[_spl_add(m,jm1,3,m,n)]);
            temp[j] = diagin * sdiag2;
            diag1 = diag2;
            sdiag1 = sdiag2;
            j = j + 1;
        endo;

    endif;

    diagin = 1 / (diag1 - sdiag1 * temp[nm1]);
    i = 1;
    do until i > m;
        npi = n + i;
        zp[_spl_add(i,n,1,m,n)] = diagin * (temp[npi] - zp[_spl_add(i,n,1,m,n)]
                - sdiag1 * zp[_spl_add(i,nm1,1,m,n)]);
        i = i + 1;
    endo;

    zp[_spl_add(1,n,3,m,n)] = diagin * (zxy1ns - zp[_spl_add(1,n,3,m,n)]
         - sdiag1 * zp[_spl_add(1,nm1,3,m,n)]);
    temp[n] = diagin * (zxymns - zp[_spl_add(m,n,3,m,n)] - sdiag1
        * zp[_spl_add(m,nm1,3,m,n)]);

/*
** perform back substitution
*/

    j = 2;
    do until j > n;
        jbak = np1 - j;
        jbakp1 = jbak + 1;
        t = temp[jbak];
        i = 1;
        do until i > m;

            zp[_spl_add(i,jbak,1,m,n)] = zp[_spl_add(i,jbak,1,m,n)]
              - t * zp[_spl_add(i,jbakp1,1,m,n)];
            i = i + 1;
        endo;
        zp[_spl_add(1,jbak,3,m,n)] = zp[_spl_add(1,jbak,3,m,n)] - t
                 * zp[_spl_add(1,jbakp1,3,m,n)];
        temp[jbak] = zp[_spl_add(m,jbak,3,m,n)] - t * temp[jbakp1];
        j = j + 1;
    endo;

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

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

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

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

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

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

    temp[n+1] = diagi * sdiag1;

    if m /= 2;
        i = 2;
        do until i > mm1;
            im1 = i - 1;
            ip1 = i + 1;
            npi = n + i;
            del2 = x[ip1] - x[i];
            if del2 <= 0;
                ierr = 2;
                retp(error(0));
            endif;
            deli = 1 / del2;
            j = 1;
            do until j > n;
                zp[_spl_add(ip1,j,2,m,n)] = deli * (z[ip1,j] - z[i,j]);
                zp[_spl_add(ip1,j,3,m,n)] = deli * (zp[_spl_add(ip1,j,1,m,n)] -
                    zp[_spl_add(i,j,1,m,n)]);
                j = j + 1;
            endo;

            { diag2,sdiag2 } = _spl_terms(sigmax,del2);

            diagin = 1 / (diag1 + diag2 - sdiag1 * temp[npi-1]);

            j = 1;
            do until j > n;
                zp[_spl_add(i,j,2,m,n)] = diagin *
                     (zp[_spl_add(ip1,j,2,m,n)] -
                    zp[_spl_add(i,j,2,m,n)] -
                    sdiag1 * zp[_spl_add(im1,j,2,m,n)]);
                zp[_spl_add(i,j,3,m,n)] = diagin *
                   (zp[_spl_add(ip1,j,3,m,n)] -
                    zp[_spl_add(i,j,3,m,n)] - sdiag1 *
                    zp[_spl_add(im1,j,3,m,n)]);
                j = j + 1;
            endo;
            temp[npi] = diagin * sdiag2;
            diag1 = diag2;
            sdiag1 = sdiag2;

            i = i + 1;
        endo;
    endif;

    diagin = 1 / (diag1-sdiag1*temp[npm-1]);
    j = 1;
    do until j > n;
        npmpj = npm + j;
        zp[_spl_add(m,j,2,m,n)] = diagin * (temp[npmpj]
              - zp[_spl_add(m,j,2,m,n)] -
            sdiag1 * zp[_spl_add(mm1,j,2,m,n)]);
        zp[_spl_add(m,j,3,m,n)] = diagin * (temp[j]
         - zp[_spl_add(m,j,3,m,n)] -
            sdiag1 * zp[_spl_add(mm1,j,3,m,n)]);
        j = j + 1;
    endo;

/*
** perform back substitution
*/

    i = 2;
    do until i > m;

        ibak = mp1 - i;
        ibakp1 = ibak + 1;
        npibak = n + ibak;
        t = temp[npibak];

        j = 1;
        do until j > n;

            zp[_spl_add(ibak,j,2,m,n)] =
                 zp[_spl_add(ibak,j,2,m,n)] -
                  t * zp[_spl_add(ibakp1,j,2,m,n)];
            zp[_spl_add(ibak,j,3,m,n)] =
                 zp[_spl_add(ibak,j,3,m,n)] -
                  t * zp[_spl_add(ibakp1,j,3,m,n)];
            j = j + 1;
        endo;
        i = i + 1;
    endo;
    retp(zp);
endp;

proc _spl_surf2(xx,yy,x,y,z,zp,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,sigmax,sigmay,jm1,j,i,del1,del2,dels,im1,zim1,zi,zxxim1,zxxi,
        dummy,sinhm1,sinhm2,sinhms;

    m = rows(x);
    n = rows(y);

/*
** denormalize tension factor in x and y direction
*/
    sigmax = abs(sigma) * (m - 1) / (x[m] - x[1]);
    sigmay = abs(sigma) * (n - 1) / (y[n] - y[1]);

/*
** determine y interval
*/

    jm1 = _spl_intrvl(yy,y);

    j = jm1 + 1;

/*
** determine x interval
*/

    im1 = _spl_intrvl(xx,x);
    i = im1 + 1;
    del1 = yy - y[jm1];
    del2 = y[j] - yy;
    dels = y[j]-y[jm1];

/*
** perform four interpolations in y-direction
*/
    if sigmay == 0;
        zim1 = _spl_hermz(z[i-1,j-1],z[i-1,j],zp[_spl_add(i-1,j-1,1,m,n)],
            zp[_spl_add(i-1,j,1,m,n)],del1,del2,dels);
        zi = _spl_hermz(z[i,j-1],z[i,j],zp[_spl_add(i,j-1,1,m,n)],
                  zp[_spl_add(i,j,1,m,n)],del1,del2,dels);
        zxxim1 = _spl_hermz(zp[_spl_add(i-1,j-1,2,m,n)],
                 zp[_spl_add(i-1,j,2,m,n)],zp[_spl_add(i-1,j-1,3,m,n)],
                 zp[_spl_add(i-1,j,3,m,n)],del1,del2,dels);
        zxxi = _spl_hermz(zp[_spl_add(i,j-1,2,m,n)],zp[_spl_add(i,j,2,m,n)],
              zp[_spl_add(i,j-1,3,m,n)],zp[_spl_add(i,j,3,m,n)],del1,del2,dels);
    else;

        { sinhm1,dummy } = _spl_snhcsh(sigmay*del1,-1);
        { sinhm2,dummy } = _spl_snhcsh(sigmay*del2,-1);
        { sinhms,dummy } = _spl_snhcsh(sigmay*dels,-1);

        zim1 = _spl_hermnz(z[i-1,j-1],z[i-1,j],zp[_spl_add(i-1,j-1,1,m,n)],
            zp[_spl_add(i-1,j,1,m,n)],sigmay,del1,del2,dels,sinhm1,sinhm2,
            sinhms);
        zi = _spl_hermnz(z[i,j-1],z[i,j],zp[_spl_add(i,j-1,1,m,n)],
             zp[_spl_add(i,j,1,m,n)],sigmay,del1,del2,dels,sinhm1,
             sinhm2,sinhms);
        zxxim1 = _spl_hermnz(zp[_spl_add(i-1,j-1,2,m,n)],
             zp[_spl_add(i-1,j,2,m,n)],zp[_spl_add(i-1,j-1,3,m,n)],
             zp[_spl_add(i-1,j,3,m,n)],sigmay, del1,del2,
            dels,sinhm1,sinhm2,sinhms);
        zxxi = _spl_hermnz(zp[_spl_add(i,j-1,2,m,n)],zp[_spl_add(i,j,2,m,n)],
            zp[_spl_add(i,j-1,3,m,n)],zp[_spl_add(i,j,3,m,n)],sigmay, del1,
             del2,dels,sinhm1,sinhm2,sinhms);
    endif;

/*
** perform final interpolation in x-direction
*/
    del1 = xx - x[im1];
    del2 = x[i] - xx;
    dels = x[i] - x[im1];

    if sigmax == 0;
        retp(_spl_hermz(zim1,zi,zxxim1,zxxi,del1,del2,dels));
    endif;
    { sinhm1,dummy } = _spl_snhcsh(sigmax*del1,-1);
    { sinhm2,dummy } = _spl_snhcsh(sigmax*del2,-1);
    { sinhms,dummy } = _spl_snhcsh(sigmax*dels,-1);

    retp(_spl_hermnz(zim1,zi,zxxim1,zxxi,sigmax,del1,del2,dels,
             sinhm1,sinhm2,sinhms));
endp;

proc(3) = _spl_ceez (del1,del2,sigma,n);
/*
**    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 del,delp,delm,denom,dummy,coshm1,coshm2,sinhmp,sinhmm;

    if n == 2;
/*
** two coefficients
*/
        retp(-1/del1,1/del1,0);
    endif;

    if sigma == 0;
/*
** tension .eq. 0.
*/
        del = del2 - del1;
        retp(-(del1+del2)/(del1*del2),del2/(del1*del),-del1/(del2*del));
    else;
/*
** tension .ne. 0.
*/
        { dummy,coshm1 } = _spl_snhcsh(sigma*del1,1);
        { dummy,coshm2 } = _spl_snhcsh(sigma*del2,1);
        delp = sigma * (del2 + del1) / 2;
        delm = sigma * (del2 - del1) / 2;
        { sinhmp,dummy } = _spl_snhcsh(delp,-1);
        { sinhmm,dummy } = _spl_snhcsh(delm,-1);

        denom = coshm1 * (del2 - del1) - 2 * del1 * delp * delm * (1 +
            sinhmp) * (1 + sinhmm);
        retp(2*delp*delm*(1+sinhmp)*(1+sinhmm)/denom,-coshm2/denom,
            coshm1/denom);
    endif;
endp;

proc(2) = _spl_terms(sigma,del);
/*
**                                 coded 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
**
** this subroutine computes the diagonal and superdiagonal
** terms of the tridiagonal linear system associated with
** spline under tension interpolation.
**
** on input--
**
**   sigma contains the tension factor.
**
** and
**
**   del contains the step size.
**
** on output--
**
**                sigma*del*cosh(sigma*del) - sinh(sigma*del)
**   diag = del*--------------------------------------------.
**                     (sigma*del)**2 * sinh(sigma*del)
**
**                   sinh(sigma*del) - sigma*del
**   sdiag = del*----------------------------------.
**                (sigma*del)**2 * sinh(sigma*del)
**
** and
**
**   sigma and del are unaltered.
**
** this subroutine references package module _spl_snhcsh.
*/
    local sigdel,denom,sinhm,coshm;

    if sigma == 0;
        retp(del/3,del/6);
    else;
        sigdel = sigma * del;
        { sinhm,coshm } = _spl_snhcsh(sigdel,0);
        denom = sigma * sigdel * (1 + sinhm);
        retp((coshm-sinhm)/denom,sinhm/denom);
    endif;

endp;

proc(2) = _spl_snhcsh(x,isw);
/*
**    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 data,sinhm,coshm,ax,xs,expx;

    local sp13,sp12,sp11,sp10,sp24,sp23,sp22,sp21,sp20,sp33,sp32,sp31,sq32,
        sq31,sq30,sp43,sp42,sp41,sq42,sq41,sq40,cp4,cp3,cp2,cp1,cp0;

    clear sinhm, coshm;
    sp13 = .3029390e-5;
    sp12 = .1975135e-3;
    sp11 = .8334261e-2;
    sp10 = .1666665e0;
    sp24 = .3693467e-7;
    sp23 = .2459974e-5;
    sp22 = .2018107e-3;
    sp21 = .8315072e-2;
    sp20 = .1667035e0;
    sp33 = .6666558e-5;
    sp32 = .6646307e-3;
    sp31 = .4001477e-1;
    sq32 = .2037930e-3;
    sq31 = -.6372739e-1;
    sq30 = .6017497e1;
    sp43 = .2311816e-4;
    sp42 = .2729702e-3;
    sp41 = .9868757e-1;
    sq42 = .1776637e-3;
    sq41 = -.7549779e-1;
    sq40 = .9110034e1;
    cp4 = .2982628e-6;
    cp3 = .2472673e-4;
    cp2 = .1388967e-2;
    cp1 = .4166665e-1;
    cp0 = .5000000e0;

    ax = abs(x);
    if isw >= 0;
        if isw > 2;
            xs = ax * ax;
            if ax > 2.3;
                expx = exp(ax);
                coshm = ( (expx + 1 / expx - xs ) / 2 - 1) / xs;
                if isw == 3;
                    sinhm = (expx - 1 / expx) / (ax + ax) - 1;
                endif;
            else;
                coshm = xs*(((cp4*xs+cp3)*xs+cp2)*xs+cp1);
                if isw == 3;
                    sinhm = xs * (((sp13 * xs + sp12) * xs + sp11) * xs +
                        sp10);
                endif;
            endif;
        else;
            if ax > 2.3;
                expx = exp(ax);
                coshm = (expx + 1 / expx) / 2 - 1;