Regridding

Because of the weighting $\bold W$,which is a function of the residual itself, the fitting problems (22) and (23) are nonlinear. Thus a nonlinear solver is required. Unlike linear solvers, nonlinear solvers need a good starting approximation so they do not land in a false minimum. (Linear solvers benefit too by converging more rapidly when started from a good approximation.) I chose the starting solution $\bold h_0$beginning from median binning on a coarse mesh. Then I refined the mesh with linear interpolation.

The regridding chore reoccurs on many occasions so I present reusable code. When a continuum is being mapped to a mesh, it is best to allocate to each mesh point an equal area on the continuum. Thus we take an equal interval between each point, and a half an interval beyond the end points. Given n points, there are n-1 intervals between them, so we have

min = o - d/2
max = o + d/2 + (n-1)*d

which may be back solved to

d = (max-min)/n
o = (min*(n-.5) + max/2)/n

which is a memorable result for d and a less memorable one for o. With these not-quite-trivial results, we can invoke the linear interpolation operator lint2. It is designed for data points at irregular locations, but we can use it for regular locations too. Operator refine2 defines pseudoirregular coordinates for the bin centers on the fine mesh and then invokes lint2 to carry data values from the coarse regular mesh to the pseudoirregular finer mesh. Upon exiting from refine2, the data space (normally irregular) is a model space (always regular) on the finer mesh.  

#$=head1 NAME
#$
#$refine2 -refine mesh
#$
#$=head1 SYNOPSIS
#$
#$C<call refine2_init(co1,cd1,co2,cd2,m1,m2,fo1,fd1,fo2,fd2,n1,n2)>
#$  
#$C<ierr= refine2_lop(adj,model,data)>  
#$  
#$C<call refine2_close()>
#$  
#$=head1 PARAMETERS
#$  
#$=over 4
#$
#$=item co1,co2 - real
#$  
#$      Coarse grain origin
#$
#$=item cd1,cd2 - real
#$  
#$      Coarse grain sampling
#$  
#$=item ro1,ro2 - real
#$
#$      Refined grain origin
#$  
#$=item rd1,rd2 - real
#$
#$      Refined grain sampling
#$  
#$=item m1,m2   - integer
#$  
#$      Size of coarse grain input
#$  
#$=item n1,n2   - integer
#$  
#$      Size of fine grain output
#$
#$=item adj,add,model,data -
#$     
#$      Standard operator interface
#$     
#$=back
#$
#$=head1 DESCRIPTION
#$
#$ Refine a mesh by linear interpolation
#$  
#$=head1 SEE ALSO
#$
#$L<lint2>
#$
#$  
#$=head1 LIBRARY
#$   
#$B<geef90>
#$  
#$=cut
#$

module refine2 { # Refine mesh.
# Input mm(m1,m2) is coarse.  Output dd(n1,n2) linear interpolated.
#
use lint2
real, dimension( :, :), pointer, private :: xy
#% _init( co1,cd1,co2,cd2, m1,m2, fo1,fd1,fo2,fd2, n1,n2)
  integer, intent( in)  :: m1,m2, n1,n2
  real,    intent( in)  :: co1,cd1,co2,cd2                 # coarse grid
  real,    intent( out) :: fo1,fd1,fo2,fd2                 # fine grid
  integer               :: i1,i2, id
  real                  :: xmin,xmax, ymin,ymax, x,y
  allocate (xy( n1*n2, 2))
  xmin = co1-cd1/2;   xmax = co1 +cd1*(m1-1) +cd1/2        # Great formula!
  ymin = co2-cd2/2;   ymax = co2 +cd2*(m2-1) +cd2/2
  fd1= (xmax-xmin)/n1;  fo1= (xmin*(n1-.5) + xmax/2)/n1    # Great formula!
  fd2= (ymax-ymin)/n2;  fo2= (ymin*(n2-.5) + ymax/2)/n2
  do i2=1,n2 {  y = fo2 + fd2*(i2-1)
  do i1=1,n1 {  x = fo1 + fd1*(i1-1)
        id = i1+n1*(i2-1)
        xy( id, :) = (/ x, y /)
        }}
  call lint2_init( m1,m2, co1,cd1, co2,cd2, xy)
#% _lop (mm, dd) 
  integer stat1
  stat1 = lint2_lop( adj, .true., mm, dd)
#% _close
deallocate (xy)
}

Finally, here is the 2-D linear interpolation operator lint2, which is a trivial extension of the 1-D version lint1 .  

#$=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)
   }
}