SHARED SUBROUTINES

The following smoothing subroutines are described in PVI and used in both PVI and BEI.  

subroutine boxconv( nb, nx, xx, yy)
# inputs:       nx,  xx(i), i=1,nx      the data
#               nb                      the box length
# output:       yy(i),i=1,nx+nb-1       smoothed data
integer nx, ny, nb, i
real xx(nx), yy(1)
temporary real bb(nx+nb)
if( nb < 1 || nb > nx)  call erexit('boxconv')  # "||" means .OR.
ny = nx+nb-1
do i= 1, ny
        bb(i) = 0.
bb(1) = xx(1)
do i= 2, nx
                bb(i) = bb(i-1) + xx(i)         # make B(Z) = X(Z)/(1-Z)
do i= nx+1, ny
                bb(i) = bb(i-1)
do i= 1, nb
                yy(i) = bb(i)
do i= nb+1, ny
                yy(i) = bb(i) - bb(i-nb)        # make Y(Z) = B(Z)*(1-Z**nb)
do i= 1, ny
                yy(i) = yy(i) / nb
return; end

 

# Convolve with triangle
#
subroutine triangle( nr, m1, n12, uu, vv)
# input:        nr      rectangle width (points) (Triangle base twice as wide.)
# input:        uu(m1,i2),i2=1,n12      is a vector of data.
# output:       vv(m1,i2),i2=1,n12      may be on top of uu
integer nr,m1,n12, i,np,nq
real uu( m1, n12),  vv( m1, n12)        
temporary real pp(n12+nr-1), qq(n12+nr+nr-2), tt(n12)
do i=1,n12 {    qq(i) = uu(1,i) }
if( n12 == 1 )
        do i= 1, n12
                tt(i) = qq(i)
else {
        call boxconv( nr, n12, qq, pp);         np = nr+n12-1
        call boxconv( nr, np , pp, qq);         nq = nr+np-1
        do i= 1, n12
                        tt(i) = qq(i+nr-1)
        do i= 1, nr-1                                   # fold back near end
                        tt(i) = tt(i) + qq(nr-i)
        do i= 1, nr-1                                   # fold back far end
                        tt(n12-i+1) = tt(n12-i+1) + qq(n12+(nr-1)+i)
        }
do i=1,n12 {    vv(1,i) = tt(i) }
return; end

 

# smooth by convolving with triangle in two dimensions.
#
subroutine triangle2( rect1, rect2, n1, n2, uu, vv)
integer i1,i2,        rect1, rect2, n1, n2
real uu(n1,n2), vv(n1,n2)
temporary real ss(n1,n2)
do i1= 1, n1
                call triangle( rect2, n1, n2, uu(i1,1), ss(i1,1))
do i2= 1, n2
                call triangle( rect1,  1, n1, ss(1,i2), vv(1,i2))
return; end