cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cc                                                                    cc
cc                 mnpola : polar test program                        cc
cc                                                                    cc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      real x(200),y(200),z(200),w(200),u(200),v(200),tu(30),tv(30),
     * c(300),exact(200),f(200),wrk1(15000),wrk2(5700)
      integer iopt(3),iwrk(500)
      integer i,is,ier,kwrk,l,lwrk1,lwrk2,l1,l2,m,m1,m2,nc,nu,nv,
     * nuest,nvest
      real s,fp,eps,sum,ermax,error,ai
      real abs,rad1,rad2,testpo,evapol
      external rad1,rad2
c  we fetch the number of data points.
      m1 = 200
      m2 =  90
c  we fetch and print the coordinates and function values of the data.
      write(6,900)
      write(6,905)
      l2 = 0
      do 10 i=1,50
         l1 = l2+1
         l2 = l2+4
         read(5,910) (x(l),y(l),z(l),l=l1,l2)
  10  continue
      write(6,915)(x(l),y(l),z(l),l=1,m1)
c  we calculate the exact function values and set up the weights w(i)=
c  (0.01)**(-1) (0.01 is an estimate for the standard deviation of the
c  error in z(i)). at the same time we calculate the mean and maximum
c  errors for the data values.
      sum = 0.
      ermax = 0.
      do 20 i=1,m1
         w(i) = 0.1e03
         exact(i) = testpo(x(i),y(i))
         error = abs(z(i)-exact(i))
         sum = sum+error
         if(error.gt.ermax) ermax = error
  20  continue
      ai = m1
      sum = sum/ai
      write(6,920) sum,ermax
c  we set up the dimension information
      nuest = 15
      nvest = 19
      lwrk1 = 15000
      lwrk2 = 5700
      kwrk = 500
c  we choose a value for eps
      eps = 0.1e-05
c  main loop for the different spline approximations
      do 400 is=1,5
        go to (110,120,130,140,160),is
c  we determine a number of smoothing spline approximations on the unit
c  disk x**2+y**2 <= 1.
c  all the data points are considered.
 110    m = m1
c  we set up the smoothing factor.
        s = 1500.
c  the approximations are not restricted at the boundaries of the disk
        iopt(3) = 0
c  we request c2-continuity at the origin.
        iopt(2) = 2
c  at the first call of polar iopt(1) must be zero.
        iopt(1) = 0
        go to 200
c  iopt(1) = 1 from the second call on
 120    iopt(1) = 1
        s = 200.
        go to 200
 130    s = 170.
        go to 200
c  we determine a smoothing spline approximation on the ellips
c  3*x**2+3*y**2-4*x*y<=1.
c  we only consider the data points inside this domain.
 140    m = m2
        ai = m
c  the given function has a constant value 0.4 at the boundary of the
c  ellips. we calculate new data values by substracting this constant
c  from the old ones.
        do 150 i=1,m
          z(i) = z(i)-0.4
 150    continue
c  given these data we will then determine approximations which are
c  identically zero at the boundary of the ellips.
        iopt(3) = 1
c  we still request c2-continuity at the origin.
        iopt(2) = 2
c  reinitialization for the knots.
        iopt(1) = 0
c  we set up the smoothing factor.
        s = 90.
        go to 250
c  at the last call we will determine the least-squares spline
c  approximation corresponding to the current set of knots
 160    iopt(1) = -1
        go to 250
c  determination of the spline approximation on the disk
 200    call polar(iopt,m,x,y,z,w,rad1,s,nuest,nvest,eps,nu,tu,
     *   nv,tv,u,v,c,fp,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ier)
        nc = (nu-4)*(nv-4)
c  we calculate the function values at the different points.
        do 220 i=1,m
            f(i) = evapol(tu,nu,tv,nv,c,rad1,x(i),y(i))
 220    continue
        write(6,925) s
        go to 300
c  determination of the spline approximation on the ellips.
 250    call polar(iopt,m,x,y,z,w,rad2,s,nuest,nvest,eps,nu,tu,
     *   nv,tv,u,v,c,fp,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ier)
c  we determine the b-spline coefficients for the spline approximations
c  of the given function.
        nc = (nu-4)*(nv-4)
        do 260 i=1,nc
            c(i) = c(i)+0.4
 260    continue
c  we calculate the function values at the different points.
        do 270 i=1,m
            f(i) = evapol(tu,nu,tv,nv,c,rad2,x(i),y(i))
 270    continue
        if(iopt(1).lt.0) go to 280
        write(6,930) s
        go to 300
 280    write(6,935)
 300    write(6,940) fp,ier
        write(6,945) nu
        write(6,950)
        write(6,955) (tu(i),i=1,nu)
        write(6,960) nv
        write(6,950)
        write(6,955) (tv(i),i=1,nv)
        write(6,965)
        write(6,970) (c(i),i=1,nc)
c  we determine mean and maximum errors.
        sum = 0.
        ermax = 0.
        do 350 i=1,m
          error = abs(f(i)-exact(i))
          sum = sum+error
          if(error.gt.ermax) ermax = error
 350    continue
        sum = sum/ai
        write(6,975)
        write(6,980)
        write(6,915)(x(l),y(l),f(l),l=2,m,3)
        write(6,920) sum,ermax
 400  continue
      stop
c  format statements
 900  format(15h1the input data)
 905  format(1h0,3(3x,1hx,6x,1hy,6x,1hz,5x))
 910  format(12f6.3)
 915  format(1h ,3(3f7.3,2x))
 920  format(14h0mean error = ,f7.4,5x,13hmax. error = ,f7.4)
 925  format(38h0smoothing spline on the disk with s =,f5.0)
 930  format(40h0smoothing spline on the ellips with s =,f5.0)
 935  format(35h0least-squares spline on the ellips)
 940  format(27h0sum of squared residuals =,e15.6,5x,12herror flag =,i5)
 945  format(1x,42htotal number of knots in the u-direction =,i3)
 950  format(1x,22hposition of the knots )
 955  format(5x,8f8.4)
 960  format(1x,42htotal number of knots in the v-direction =,i3)
 965  format(23h0b-spline coefficients )
 970  format(5x,8f9.4)
 975  format(33h0spline values at selected points)
 980  format(1h0,3(3x,1hx,6x,1hy,6x,1hf,5x))
      end
      real function rad1(v)
c  function program rad1 defines in polar coordinates, the boundary of
c  the approximation domain  x**2+y**2<=1.
      real v
      rad1 = 1.
      return
      end
      real function rad2(v)
c  function program rad2 defines in polar coordinates, the boundary of
c  the approximation domain  3*x**2+3*y**2-4*x*y<=1.
      real sqrt,sin,v
      rad2 = 1./sqrt(3.-2.*sin(v+v))
      return
      end
      real function testpo(x,y)
c  function program testpo evaluates the test function for the polar
c  package.
      real x,y
      testpo=(x**2+y**2)/((x+y)**2+0.5)
      return
      end