#$
#$=head1 NAME
#$
#$nhelicon - non stationary convolution
#$
#$=head1 SYNOPSIS
#$
#$Initializer - C<call nhelicon_init(aa)>
#$
#$Operator    - C<ierr=nhelicon_lop(adj,add,xx,yy)>
#$
#$=head1 PARAMETERS
#$
#$=over 4
#$
#$=item aa - type(nfilter)  
#$
#$      nhelix filter to perform convolution with
#$
#$=item adj,add,xx,yy - 
#$
#$      standard operators parameters
#$
#$=back
#$
#$=head1 DESCRIPTION
#$
#$ Nonstationary convolution, inverse to deconvolution.
#$ Requires the filter be causal with an implicit "1." at the onset.
#$
#$=head1 SEE ALSO
#$
#$L<nhelix>,L<nhconest>,L<npolydiv>,L<helicon>
#$
#$=head1 LIBRARY
#$
#$B<geef90>
#$
#$=cut

module nhelicon {  # Nonstationary convolution, inverse to deconvolution.
#             Requires the filter be causal with an implicit "1." at the onset.
use nhelix
type( nfilter) :: aa
#%  _init( aa)
#%  _lop ( xx,  yy)
integer                         :: iy, ix, ia, ip
integer, dimension( :), pointer :: lag
real,    dimension( :), pointer :: flt
if( adj)              # zero lag
	xx += yy
else
	yy += xx
do iy = 1, size( yy) {    	
	if( associated( aa%mis)) { if( aa%mis( iy)) cycle}
        ip = aa%pch( iy); lag => aa%hlx( ip)%lag; flt => aa%hlx( ip)%flt
	do ia = 1, size( lag) {
            ix = iy - lag( ia); if( ix < 1) cycle;
            if( adj)
                        xx(ix) += yy(iy) * flt( ia)
	    else	
                        yy(iy) += xx(ix) * flt( ia)
	    }
      }
}