      subroutine stplrg (m,x,y,n,c,iflag,a)                             STPL001
c --------------------------------------------------------------        STPL002
c      fits a polynomial of degree n using direct lsq method.           STPL003
c      A was dimensioned parametric before. Needs MATINV                STPL004
c ---------------------------------------------------------------       STPL005
c      Generic version, last revised 23.2.89 ETH-LNS                    STPL006
c ---------------------------------------------------------------       STPL007
      dimension x(1), y(1), a(1), c(1)                                  STPL008
      n1 = n + 1                                                        STPL009
      n2 = n1 + 1                                                       STPL010
      nnpn = n*n1                                                       STPL011
      n1pn1 = n1 + n1                                                   STPL012
      n1n1 = nnpn + n1                                                  STPL013
                                                                        STPL014
c ------- clear 1st row and last two columns                            STPL015
                                                                        STPL016
      i1j = 1                                                           STPL017
      do 30 j = 1,n                                                     STPL018
         a(i1j) = 0.0                                                   STPL019
 30   i1j = i1j + n1                                                    STPL020
      l = nnpn + 1                                                      STPL021
      do 40 j = 1,n1pn1                                                 STPL022
         a(l) = 0.0                                                     STPL023
 40   l = l + 1                                                         STPL024
                                                                        STPL025
c -------- compute sums of powers                                       STPL026
                                                                        STPL027
      do 10 k = 1,m                                                     STPL028
         xpwr = x(k)                                                    STPL029
         i1j = n1 + 1                                                   STPL030
         do 20 j = 2,n1                                                 STPL031
            jnp2 = n1n1 + j                                             STPL032
            a(i1j) = a(i1j) + xpwr                                      STPL033
            a(jnp2) = a(jnp2) + xpwr*y(k)                               STPL034
         xpwr = xpwr*x(k)                                               STPL035
 20      i1j = i1j + n1                                                 STPL036
         do 50 i = 2,n1                                                 STPL037
            in1 = nnpn + i                                              STPL038
            a(in1) = a(in1) + xpwr                                      STPL039
 50      xpwr = xpwr*x(k)                                               STPL040
 10   a(n1n1 + 1) = a(n1n1 + 1) + y(k)                                  STPL041
                                                                        STPL042
c -------- fill in rest of matrix a                                     STPL043
                                                                        STPL044
      a(1) = m                                                          STPL045
      do 60 i = 2,n1                                                    STPL046
         ij = i                                                         STPL047
         do 60 j = 1,n                                                  STPL048
            im1jp1 = ij + n                                             STPL049
            a(ij) = a(im1jp1)                                           STPL050
 60   ij = ij + n1                                                      STPL051
                                                                        STPL052
c -------- solve equations for coefficients (call matrix inversion)     STPL053
                                                                        STPL054
      nc = -n2                                                          STPL055
      call matinv (iflag,iflag,n1,nc,a,n1,i,xpwr)                       STPL056
                                                                        STPL057
      k = n1n1 + n1                                                     STPL058
      do 70 i = 1,n1                                                    STPL059
         c(i) = a(k)                                                    STPL060
 70   k = k - 1                                                         STPL061
      return                                                            STPL062
      end                                                               STPL063
                                                                        STPL064
