Using a PEF on the Galilee residual

For simplicity, we begin with a simple gradient for the PEF of the map. We have

$$\begin{array} {lll} 0 &\approx & \bold A ( \bold B \bold h -\bold d ) \\ 0 &\approx & \epsilon \nabla \bold h \end{array}$$ (6)

This is our first fitting system that involves all the raw data. Previous ones have involved the data only after binning. Dealing with all the raw data, we can expect even more difficulty with impulsive and erratic noises. The way to handle such noise is via the weighting functions of Chapter . Including such a weighting function gives us the map-fitting goals,  

$$\begin{array} {lll} 0 &\approx & \bold W \bold A ( \bold B ... ...\bold d ) \\ 0 &\approx & \epsilon \nabla \bold h \end{array}$$ (7)

The weighting function depends on the residual itself, via such equations as () or (). The first time I solved the problem, I chose the PEF $\bold A$ to be the derivative (1,-1) along the data track. To avoid an excessive influence of noise bursts, I applied row-normalized PEF estimation using subroutine rnpef1() . I tried various scalings between the first and second parts of (7) and various values of the average residual $\bar r$ required in the weights () and (). Some sample noise PEFs I found were

Noise PEFs in the problem:   min ||WA(Bm-d)|| + ||Gm||
 sum     PEF of residual.
 0.10    1.00 -0.34 -0.45 -0.06  0.01 -0.03 -0.01 -0.03 -0.00  0.00
 0.09    1.00 -0.28 -0.48 -0.09  0.01 -0.03 -0.01 -0.03 -0.00 -0.00
 0.05    1.00 -0.44 -0.46 -0.05  0.03 -0.03 -0.00 -0.02  0.02  0.02

Involving as it does so many aspects of this book, the solver potato() is quite a clutter. Sorry about that. Maybe it will look cleaner after I learn more about all the features in Fortran 90.  

subroutine potato( xy,       dd,nd, o1,d1,o2,d2, mm,m1,m2, niter)
integer iter, i,      nd,                           m1,m2, niter
real               xy(nd,2), dd(nd),o1,d1,o2,d2, mm(m1,m2), rbar, aa(100)
integer nr, rd,rm, na
temporary real               rr( nd + 2*m1*m2), ad(nd), wr(nd), ww(nd)
temporary real  dm(m1*m2),   dr( nd + 2*m1*m2)
temporary real  sm(m1*m2),   sr( nd + 2*m1*m2)
                             nr= nd + 2*m1*m2;   rd=1;  rm=1+nd
na=4;  aa(1)=1.;  aa(2)= -.44;  aa(3)= -.46;  aa(4)=-.05
call scaleit( 2.,  na,  aa)
do iter= 0, niter {
        call negcopy(     nd,          dd,          ad       ) #        -d
        call dpbin2( 0, 1, o1,d1,o2,d2,xy,mm,m1,m2, ad    ,nd) # r_d=   -d+Bm
        rbar=.5; rbar=2.5; rbar=1.0
        do i=1,nd
                ww(i) = 1./sqrt( 1. + (ad(i)/rbar)**2 )
        call icai1(  0, 0, 1, aa,na,        ad ,nd, rr(rd)   ) # r_d= A(-d+Bm)
        call diag(   0, 0, ww, nd,          rr(rd), rr(rd)   ) # r_d=WA(-d+Bm)
        call igrad2( 0, 0,   m1,m2,       mm,       rr(rm)   ) # r_m =    G m
        call diag(   1, 0, ww, nd,        wr,       rr(rd)   ) #          Wr_d
        call icai1(  1, 0, 1, aa,na,      ad ,nd,   wr       ) # dm =   A'Wr_d
        call dpbin2( 1, 0, o1,d1,o2,d2,xy,dm,m1,m2, ad    ,nd) # dm = B'A'r_d
        call igrad2( 1, 1,   m1,m2,       dm,       rr(rm)   ) # dm= dm + G'r_m
        call dpbin2( 0, 0, o1,d1,o2,d2,xy,dm,m1,m2, ad    ,nd) # dr_d =  B dm
        call icai1(  0, 0, 1, aa,na,      ad,nd,    wr       ) # dr_d = AB dm
        call diag(   0, 0, ww, nd,        wr,       dr(rd)   ) # dr_d =WAB dm
        call igrad2( 0, 0,   m1,m2,       dm,       dr(rm)   ) # dr_m =  G dm
        call cgplus( iter, m1*m2, mm, dm, sm, _
                              nr, rr, dr, sr )
        }
return; end

Results are in Figure 7.

 

potato
Figure 7 Gradient of the Galilee map from (7).

It is pleasing to see the ship's tracks gone at last.