Program 4 : d1intr, intrpr

      program prog4 
      implicit none
      
      real bicubic,derx,dery,derxy,rms,ranf
      external bicubic,derx,dery,derxy,rms,ranf
      real f(10,10),x(10),y(10),hx(10),hy(10)
      real fx(10,10),fy(10,10),fxy(10,10)
      real s(16,16),sderx(10,10),sdery(10,10),sderxy(10,10)
      real fi(16,16),xi(16,16),yi(16,16)
      real wa(10),wb(10),wc(10),wd(10),we(10)
      real za(16),zb(16),zc(16),zd(16)
      real work(16,4)
      integer i,j

# define the coordinates of the rectangular grid
# (where f will be evaluated) and their differences
      do i=1,10
         x(i)=0.1*i**2-5
         y(i)=0.1*i**2-5
         if (i.gt.1) then
            hx(i-1)=x(i)-x(i-1)
            hy(i-1)=y(i)-y(i-1)
         endif
      enddo

      
# compute the function f and its partial derivatives      
      do i=1,10
         do j=1,10
            f(i,j)=bicubic(x(i),y(j))
            fx(i,j)=derx(x(i),y(j))
            fy(i,j)=dery(x(i),y(j))
            fxy(i,j)=derxy(x(i),y(j))
         enddo
      enddo

# define the coordinates of the rectangular grid
# where interpolation will be done

      do i=1,16
         do j=1,16
            xi(i,j)=(i-8)*0.3+ranf()*0.1
            yi(i,j)=(j-8)*0.3+ranf()*0.1
         enddo
      enddo

# compute the real value of f at interpolation points
      do i=1,16
         do j=1,16
            fi(i,j)=bicubic(xi(i,j),yi(i,j))
         enddo
      enddo

# interpolate at point(i,j)=(xi(i,j),yi(i,j))
      call intrpr(s,16,16,f,10,10,xi,yi,x,y,sderx,sdery,sderxy,hx,hy,
     * 1,wa,wb,wc,wd,we,work) 
      print *,'For a bicubic with 16 coefficients(up to x**3*y**3):'
      print *,'rms error of intrpr on:'
      print *,'  interpolation:',rms(s,fi,16,16)
      print *,'  x derivative:',rms(sderx,fx,10,10)
      print *,'  y derivative:',rms(sdery,fy,10,10)
      print *,'  x-y derivative:',rms(sderxy,fxy,10,10)

      
# here is an alternative to intrpr:
# if the derivatives are already known
# (for example, because intrp or intrpr was called)
# d1intr should be used 
# derivatives from intrpr
      call d1intr(s,16,16,f,10,10,sderx,sdery,sderxy,xi,
     * yi,x,y,hx,hy,za,zb,zc,zd)  
      print *,'rms error of d1intr on interpolation'
      print *,'  with derivative from intrpr:',rms(s,fi,16,16)


      stop
      end

      

# this polynomial has terms up to x**3*y**3      
      real function bicubic(x,y)
      real x,y,coef0,coef1,coef2,coef3
      coef3=((1.0*y+2.0)*y+3.0)*y+4.0
      coef2=((5.0*y+6.0)*y+7.0)*y+8.0
      coef1=((9.0*y+10.0)*y+11.0)*y+12.0
      coef0=((13.0*y+14.0)*y+15.0)*y+16.0
      bicubic=((coef3*x+coef2)*x+coef1)*x+coef0
      return
      end

# returns the partial derivative 
# computed analytically of the bicubic with respect to x 
      real function derx(x,y)
      real x,y,coef0,coef1,coef2
      coef2=((1.0*y+2)*y+3.0)*y+4.0
      coef1=((5.0*y+6.0)*y+7.0)*y+8.0
      coef0=((9.0*y+10.0)*y+11.0)*y+12.0
      derx=3*coef2*x**2+2*coef1*x+coef0
      return
      end

# this function returns the partial derivative 
# computed analytically of the bicubic with respect to y
      real function dery(x,y)
      real x,y,coef3,coef2,coef1
      coef3=3.0*y**2+4.0*y+3.0
      coef2=15.0*y**2+12.0*y+7.0
      coef1=27.0*y**2+20.0*y+11.0
      coef0=39.0*y**2+28.0*y+15.0
      dery=coef3*x**3+coef2*x**2+coef1*x+coef0
      return
      end

# this function returns the partial derivative 
# computed analytically of the bicubic with respect to x and y 
      real function derxy(x,y)
      real x,y,coef2,coef1,coef0
      coef2=3.0*y**2+4.0*y+3.0
      coef1=15.0*y**2+12.0*y+7.0
      coef0=27.0*y**2+20.0*y+11.0
      derxy=3.0*coef2*x**2+2.0*coef1*x+coef0
      return
      end

# computes the rms of the difference of arrays a and b
      real function rms(a,b,n,m)
      real a(n,m),b(n,m),sum
      integer i,j,n,m
      sum=0
      do i=1,n
         do j=1,m
            sum=sum+(a(i,j)-b(i,j))**2
         enddo
      enddo
      rms=sqrt(sum/n/m)
      return
      end

      real function ranf()      
      integer rand
      
      ranf=rand()/32767.0
      return
      end