Regridding

Because of the weighting $\bold W$,which is a function of the residual itself, the fitting problems (25) and (26) 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 subroutine lint2(). It is designed for data points at irregular locations, but we can use it for regular locations too. Subroutine 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.  

# Refine mesh. Input mm(m1,m2) is coarse.  Output dd(n1,n2) linear interpolated.
#
subroutine refine2( co1,cd1,co2,cd2, mm,m1,m2,    fo1,fd1,fo2,fd2, dd, n1,n2)
integer i1,i2, id,                      m1,m2,                         n1,n2
real                                 mm(m1,m2),                    dd( n1*n2)
real                co1,cd1,co2,cd2,              fo1,fd1, fo2,fd2
#     inputs= coarse o1  d1  o2  d2  mesh.      | Outputs= fine    mesh
real xmin,xmax, ymin,ymax, x,y
temporary       real                  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, 1) = x
        xy( id, 2) = y
        }}
call lintr2( 0,0, co1,cd1, co2,cd2, xy, mm,m1,m2, dd, n1*n2)
return; end

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

# Linear interpolation 2-D
#
subroutine lint2( adj, add, o1,d1,o2,d2, xy,       mm,m1,m2,  dd, nd)
integer i, ix,iy, adj, add,                   id,     m1,m2,      nd
real    f, fx,fy, gx,gy,    o1,d1,o2,d2, xy(2,nd), mm(m1,m2), dd( nd)
call adjnull(     adj, add,                        mm,m1*m2,  dd, nd)
do id=1,nd {
        f=(xy(1,id)-o1)/d1; i=f; ix=1+i; if( 1<=ix && ix<m1){  fx=f-i; gx=1.-fx
        f=(xy(2,id)-o2)/d2; i=f; iy=1+i; if( 1<=iy && iy<m2){  fy=f-i; gy=1.-fy
                if( adj == 0)                                   
                        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)
                else {                                  
                        mm(ix  ,iy  ) = mm(ix  ,iy  ) + gx * gy * dd(id)
                        mm(ix+1,iy  ) = mm(ix+1,iy  ) + fx * gy * dd(id)
                        mm(ix  ,iy+1) = mm(ix,  iy+1) + gx * fy * dd(id)
                        mm(ix+1,iy+1) = mm(ix+1,iy+1) + fx * fy * dd(id)
                        }
                }}
        }
return; end