APPENDIX

subroutine recdfilt( conj, hilo, alpha, nt,nx, p, q)
#
#--------------------------------------
#	recursive dip filter
#--------------------------------------
#	conj  	: forward operator [ = 0 ] conjugate operator [ = 1 ]
#	hilo    : low-dip pass filter [ = 0 ] high-dip pass filter [ = 1 ]
#	alpha 	: cutoff dip
#   	nt,nx 	: dimension of array
#	p     	: input data
#	q     	: output data
#
#--------------------------------------
integer it,ix, nt,nx
real p(nt,nx), q(nt,nx)
real bl,br,endl,endr
real conj, hilo, alpha, beta
real p11,p12,p13,q11,q12,q13,p21,p22,p23,q21,q22,q23
temporary real r(nt,nx),px(nx),rhs(nx),sol(nx),e(nx),f(nx)

beta = 1./6.726

q11 =-1.-alpha 	    ;   q12 = -1./beta+2.+2.*alpha 	;q13 = -1.-alpha
p11 =-1.*(1.)	    ;   p12 = -1.*(1./beta-2.)		;p13 = -1.*(1.)
q21 = 1.-alpha	    ;   q22 = 1./beta-2.+2.*alpha 	;q23 = 1.-alpha
p21 = -1.*(-1.)	    ;   p22 = -1.*(-1./beta+2.)		;p23 = -1.*(-1.)
# boundary condition
bl = 0.;	br = 0.
endl = q23*bl + q22;	endr = q21*br + q22

if( conj == 0. ) {
	do it= 1, nt {	
	do ix= 1, nx {
		q(it,ix) = 0.; r(it,ix) = 0.;  rhs(ix)=0.
		}}
	call bicon(0.,bl,p11,p12,p13,p21,p22,p23,nt,nx,p,r) 
	do it= 1, nt {				
		if(it==1) { do ix= 1, nx { rhs(ix) = r(it,ix) } }
		else {  call tricon(bl,q11,q12,q13,nx,sol,rhs)
		 	do ix= 1, nx { rhs(ix) = r(it,ix) - rhs(ix)} 
			}
		call rtris(nx,endl,q21,q22,q23,endr,rhs,sol,e,f)
		do ix = 1,nx{ q(it,ix) = sol(ix) }
		}
	if( hilo==0.)
	do it= 1, nt {	do ix= 1, nx { q(it,ix) = q(it,ix) }}
	else
	do it= 1, nt {	do ix= 1, nx { q(it,ix) = p(it,ix) - q(it,ix) }}
	}
else	{
	do it= 1, nt {	
	do ix= 1, nx {
		p(it,ix) = 0.; r(it,ix) = 0.; rhs(ix)=0.
		}}
	do it= nt, 1, -1 {
		if(it==nt) { do ix= 1, nx { rhs(ix) = q(it,ix) } }
		else {  call tricon(bl,q11,q12,q13,nx,sol,rhs)
		 	do ix= 1, nx { rhs(ix) = q(it,ix) - rhs(ix)} 
			}
		call rtris(nx,endl,q21,q22,q23,endr,rhs,sol,e,f)
		do ix = 1,nx{ r(it,ix) = sol(ix) }
		}
	call bicon(1.,bl,p11,p12,p13,p21,p22,p23,nt,nx,p,r) 
	if( hilo==0.)
	do it= 1, nt {	do ix= 1, nx { p(it,ix) = p(it,ix) }}
	else
	do it= 1, nt {	do ix= 1, nx { p(it,ix) = q(it,ix) - p(it,ix) }}
	}
return; end

subroutine bicon(conj,bl,p11,p12,p13,p21,p22,p23,nt,nx,p,r) 
integer it,ix,nt,nx
real p(nt,nx),r(nt,nx)
real conj,bl,p11,p12,p13,p21,p22,p23
temporary real px(nx),dia(nx),offdia(nx)

do ix= 1, nx { dia(ix) = 0.; offdia(ix) = 0.  }
if( conj==0. )
	do it= 1, nt {				
		if( it==1 ) { 	do ix= 1, nx { px(ix) = p(it,ix) }
		    		call tricon(bl,p21,p22,p23,nx,px,dia)	 }
		else {  	do ix= 1, nx { px(ix) = p(it,ix) }
				call tricon(bl,p21,p22,p23,nx,px,dia) 
				do ix= 1, nx { px(ix) = p(it-1,ix) }
				call tricon(bl,p11,p12,p13,nx,px,offdia) }
		do ix= 1, nx { r(it,ix) = dia(ix) + offdia(ix) }
	}
else
	do it= nt, 1, -1{				
		if( it==nt ) { 	do ix= 1, nx { px(ix) = r(it,ix) }
		     		call tricon(bl,p21,p22,p23,nx,px,dia)  	 }
		else {  	do ix= 1, nx { px(ix) = r(it,ix) }
				call tricon(bl,p21,p22,p23,nx,px,dia) 
				do ix= 1, nx { px(ix) = r(it+1,ix) }
				call tricon(bl,p11,p12,p13,nx,px,offdia) }
		do ix= 1, nx { p(it,ix) = dia(ix) + offdia(ix) }
	}
return; end

subroutine tricon(bl,p1,p2,p3,nx,px,rhs)
integer ix,nx
real px(nx),rhs(nx)
real bl,p0,pnxp1,p1,p2,p3

p0 = bl*px(1);	pnxp1 = bl*px(nx)

rhs(1)  = p1*p0       + p2*px(1)  + p3*px(2)
do ix = 2, nx-1 {   rhs(ix) = p1*px(ix-1) + p2*px(ix) + p3*px(ix+1) }
		    rhs(nx) = p1*px(nx-1) + p2*px(nx) + p3*pnxp1
return;	 end

DOT-PRODUCT TEST RESULT

./Src/result.tex