#$=head1 NAME
#$
#$lint2  -  2-D linear interpolation
#$
#$=head1 SYNOPSIS
#$
#$Initializer : C<call lint2_init(m1,m2,o1,d1,o2,d2,xy)>
#$
#$Operator    : C<ierr=binpull1_lop(adj,add,model,data)>
#$
#$=head1 PARAMETERS
#$
#$=over 4
#$
#$=item m1   - integer 
#$
#$      number of model points,axis 1
#$
#$=item m2   - integer 
#$
#$      number of model points,axis 2
#$
#$=item o1 -  real 
#$
#$      first sample of model space,axis 1
#$
#$=item d1 -  real 
#$
#$      sampling of model space,axis 1
#$
#$=item o1 -  real 
#$
#$      first sample of model space,axis 2
#$
#$=item d1 -  real 
#$
#$      sampling of model space,axis 2
#$
#$=item xy - C<real(:,2)> 
#$
#$      location of data
#$
#$=item adj,add,model,data - 
#$
#$      standard operator
#$
#$=back
#$
#$=head1 DESCRIPTION
#$
#$2-D linear  interpolation
#$
#$=head1 SEE ALSO
#$
#$L<binpull2>,L<binpull1>,L<lint1>
#$
#$=head1 LIBRARY
#$
#$B<geef90>
#$
#$=cut
#$
module lint2 {                                      # (Bi)Linear interpolation in 2-D
integer                        :: m1,m2
real                           :: o1,d1, o2,d2
real, dimension (:,:), pointer :: xy
#%  _init (     m1,m2, o1,d1, o2,d2, xy)
#%  _lop  ( mm (m1,m2), dd (:))
integer i, ix,iy, id
real    f, fx,gx, fy,gy
do id= 1, size(dd) {
   f = (xy(id,1)-o1)/d1; i=f; ix= 1+i; if( 1>ix .or. ix>=m1) cycle; fx=f-i; gx= 1.-fx
   f = (xy(id,2)-o2)/d2; i=f; iy= 1+i; if( 1>iy .or. iy>=m2) cycle; fy=f-i; gy= 1.-fy
                if( adj) {
                        mm(ix  ,iy  ) += gx * gy * dd(id)
                        mm(ix+1,iy  ) += fx * gy * dd(id)
                        mm(ix  ,iy+1) += gx * fy * dd(id)
                        mm(ix+1,iy+1) += fx * fy * dd(id)
			}
		else
	                dd(id) = dd(id) + gx * gy * mm(ix  ,iy  ) +
                                          fx * gy * mm(ix+1,iy  ) +
                                          gx * fy * mm(ix  ,iy+1) +
                                          fx * fy * mm(ix+1,iy+1)
		
   }
}