Program 8 : bilin, coupe, defvec

 
      program prog8
      implicit none

      real bilinear,bilin,rms
      external bilinear,rms

      real z(12,16)
      real x(14),y(14),deltax,deltay,zi(14),li(14)
      real x1,y1,x2,y2
      real result1,result2,left,right,result3

      logical ok
      integer i,j

# initialize the field z
      do i=1,12
         do j=1,16
            z(i,j)=bilinear(float(i),float(j))
         enddo
      enddo   

# compute the value of z at (7.36,2.57) with three methods
      ok=.true.
      result1=bilin(z,12,16,7.36,2.57,ok)
      
      left=(1-0.57)*z(7,2)+0.57*z(7,3)
      right=(1-0.57)*z(8,2)+0.57*z(8,3)
      result2=(1-0.36)*left+0.36*right

      result3=bilinear(7.36,2.57)

      if (ok) then
         print *,'Interpolation by bilin successful.'
         print *,'The value of z at indexes (7.36,2.57):',result1,
     *    ' ( error:',abs(result3-result1)+abs(result3-result2),' )'
      else
         print *,'Interpolation failed: indexes out of limits.'
      endif

# coupe test
# define vector of coordinates
      x1=1.0
      y1=1.2
      x2=5.5
      y2=7.2
      deltax=(x2-x1)/13.0
      deltay=(y2-y1)/13.0
      call defvec(x,14,deltax,x1)
      call defvec(y,14,deltay,y1)
      
# compute the value of z at these coordinates
      do i=1,14
         zi(i)=bilinear(x(i),y(i))
      enddo
      
      call coupe(li,14,z,12,16,x1,y1,x2,y2,ok)

      if (ok) then
         print *,'Cross-section by coupe successful.'
         print *,'The rms error of coupe on interpolation is:',
     *    rms(li,zi,14)

      else
         print *,'Cross-section failed: indexes out of bound'
      endif


      stop
      end

      real function bilinear(x,y)
      real x,y
      bilinear=5.0*x*y+6.0*x+7.0*y+8.0
      return
      end

      real function rms(a,b,n)
      real a(n),b(n),sum
      integer i,n
      sum=0
      do i=1,n
         sum=sum+(a(i)-b(i))**2
      enddo
      rms=sqrt(sum/n)
      return
      end