SUBROUTINE VPNFN2 (
*
*  inputs
*
     :                  X, NX, X0, XS, IPOWER,
*
*  outputs
*
     :                  AFUNC)
*
*  Module number:
*
*  Module name:
*
*  Keyphrase:
*  ----------
*  construct basis functions for the (up to two dimensional) polynomial
*  function
*
*  Description:
*  ------------
*  basis functions for polynomial in (X-X0)/XS, where X is a vector of
*  length NDIM
*
*  Maximum power in dimension IDIM is IPOWER(IDIM), e.g. if
*  IPOWER(1)=2, IPOWER(2)=1 then afunc is
*  1, X1, X1**2, X2, X2*X1, X2*X1**2
*
*
*  FORTRAN name: VPNFN2.FOR
*
*  Keywords of accessed files and tables:
*  --------------------------------------
*  Name                         I/O     Description / Comments
*
*  Subroutines Called:
*  -------------------
*  CDBS:
*       None
*  SDAS:
*       None
*  Others:
*       None
*
*  History:
*  --------
*  Version      Date        Author          Description
*     1       12-21-88     J.-C. HSU        coding
*-------------------------------------------------------------------------------
*
*== input:
*                                       --input vector of length NX
        DOUBLE PRECISION        X(1),
*                                       --zero point offset of X
     :                          X0(1),
*                                       --scaling factor of X
     :                          XS(1)
*                                       --dimensions of X
        INTEGER                 NX,
*                                       --polynomial power in each dimension
*                                       --of X
     :                          IPOWER(1)
*
*== output:
*                                       --basis functions
        DOUBLE PRECISION        AFUNC(1)
*
*== local:
*
        DOUBLE PRECISION        DX
*                                       --loop indices
        INTEGER                 I, J, IP, IDIM, IFIRST, INUM
*------------------------------------------------------------------------------
*
*  set up the basis functions for the first dimension
*
        AFUNC(1) = 1.0D0
*
        DO 10 I = 2, IPOWER(1)+1
            AFUNC(I) = (X(1)-X0(1)) / XS(1) * AFUNC(I-1)
   10   CONTINUE
*
*  set up the basis functions for the second dimension, if any
*
        I = IPOWER(1) + 1
*
        DO 40 IDIM = 2, NX
            IFIRST = 1
            INUM = I
            DX = (X(IDIM) - X0(IDIM)) / XS(IDIM)
*
            DO 30 IP = 2, IPOWER(IDIM)+1
                DO 20 J = IFIRST, IFIRST+INUM-1
                    I = I + 1
                    AFUNC(I) = AFUNC(J) * DX
   20           CONTINUE
                IFIRST = IFIRST + INUM
   30       CONTINUE
   40   CONTINUE
*
        RETURN
        END