Seasat optimization

To prepare to minimize the fitting residuals (2), we begin by restating the residual as the output of the operator that takes a model perturbation $\Delta \bold h$and produces the weighted residual perturbation $\Delta \bold r$. 

$$\left[ \begin{array} {c} \ \\ \Delta \bold r \\ \ ... ...on \bold I \end{array} \right] \ \left[ \Delta \bold h \right]$$ (4)

where we have also split the track operator $\bold T$ into two parts, $\bold N$ making the northwest tracks and $\bold S$ making the southwest tracks.

To find the best-fitting topography $\bold h = h(x,y)$,we prepare the solver subroutine seatopo(). Begin by loading the residual vector with the negative of the data $\bold r=-\bold d$and load the solution vector with zeros $\bold h=\bold 0$.Then proceed as we did with earlier weighted residuals like pef1() or pef2() . First we revise the data to weighted data by a single application of subroutine trakwit() . Then begin the usual conjugate-direction iteration containing first the adjoint of operator (4) and then the forward operator. The only thing new here is that the filter weighting cannot have the output and input share the same computer memory, as is allowed by the more usual scalar weighting routines, ident() and diag() , so an extra array dr2 must be declared. (In a computer language better than Fortran there would be no need to list separate programs for each application of weighted least squares--one routine should meet all needs.)  

# Build topography from satellite tracks and PEF of prior topo.
#
subroutine seatopo( dkno,  dd,n1,n2, niter, eps, aa,a1,a2,  hh )
integer n12, r123, iter,      n1,n2, niter,         a1,a2
real        dkno(n1,n2,2), dd(n1,n2,2)    , eps, aa(a1,a2), hh(n1,n2)
temporary real             rr( n1*n2*3)
temporary real dh( n1*n2), dr( n1*n2*3), dr2(n1*n2*3)
temporary real sh( n1*n2), sr( n1*n2*3)
temporary real             wr( n1*n2*3)
               n12=n1*n2;              r123=n1*n2*3

call null( hh, n12);    call null( rr, r123)
call         negcopy( n12 * 2, dd, rr          )
call         trakwit( 0, aa,a1,a2,          rr, n1,n2,  wr )
do iter= 0, niter {
        call trakwit( 1, aa,a1,a2,          dr, n1,n2,  wr )
        call   track( 1,0, eps,  dh,n1,n2,  dr,              dkno )
        call   track( 0,0, eps,  dh,n1,n2,  dr,              dkno )
        call trakwit( 0, aa,a1,a2,          dr, n1,n2,  dr2 )
        call cgplus( iter,  n12,  hh, dh,  sh, _
                            r123, wr, dr2, sr )
        }
return; end

Finally, subroutine seastat() builds a topographic map from satellite data tracks; the subroutine uses either of two different statistical assumptions. The first assumption is that the 2-D prediction-error filter (PEF) for the topography is Laplace's operator. The alternate assumption is that a previous map is available which can be used by pef2() to compute the required PEF. The previous map is presumed to be poor where the tracks are missing. To put this presumption into practice, we set the map to zero where there are no data tracks, so pef2() will handle such regions as missing data. Then we define the template of free parameters with subroutine setfree() ; and we call subroutine pef2() to get the PEF. With the inputs now all prepared, we invoke subroutine seatopo().  

# Find the topo map from data and either a previous topo map or a PEF.
#
subroutine seastat( dd,n1,n2, niter, a1,a2, hh, laplace)
integer i1,i2,i3,      n1,n2, niter, a1,a2,     laplace,      lag1,lag2, iter
real                dd(n1,n2,2),            hh(n1,n2), eps
temporary real    dkno(n1,n2,2), afre(a1,a2)
temporary real     res(n1,n2  ),   aa(a1,a2)
eps= 1.;        lag1=1+a1/2;    lag2=1
do i3=1, 2{
do i2=1,n2{
do i1=1,n1  {   if( dd(i1,i2,i3) == 0.)   dkno(i1,i2,i3) = 0.0
                else                      dkno(i1,i2,i3) = 1.0 
                }}};           call null( aa, a1*a2)
if( laplace > 0 ) {                     aa(lag1-1,2) =  .25
                aa(lag1+0,1) = .25;     aa(lag1  ,2) = -1.0;    aa(lag1,3)= .25
                                        aa(lag1+1,2) =  .25                    }
else {  do i1=1,n1{
        do i2=1,n2   { 
	  if( dkno(i1,i2,1) == 0. || dkno(i1,i2,2) == 0.  )   hh(i1,i2) = 0.  }}
        call setfree(                      lag1,lag2, 1, afre,    a1,a2,1, iter)
        call pef2( hh, res, n1,n2,   iter, lag1,lag2,    afre, aa,a1,a2 )
        }
call seatopo( dkno, dd,n1,n2, niter, eps, aa,a1,a2,  hh )
return; end