# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module envelope_mod
  implicit none
  contains
  logical function calc_envelope(image,tc1,tc2)
    real :: image(:,:,:)
    integer :: tc1,tc2
    real, allocatable :: top(:,:),bot(:,:),data(:)
    complex,allocatable :: cdata(:),ctop(:)
    integer :: n1_2,n2_2,i2,i3,i1
    calc_envelope=.false.
    n1_2=1
    do while ( n1_2 < size(image,1))
      n1_2=n1_2*2
    end do
    n2_2=1
    do while ( n2_2 < size(image,2))
      n2_2=n2_2*2
    end do
    allocate(data(n1_2),top(n1_2,n2_2),bot(n1_2,n2_2),cdata(n1_2),ctop&
      &(n1_2))
    do i3=1,size(image,3)
      do i2= 1, size(image,2)
        cdata(1:size(image,1))=cmplx(image(:,i2,i3),0.)
        cdata(size(image,1)+1:)=cmplx(0,0.)
        call ftu( 1., n1_2, cdata)
        cdata(1) = cdata(1) / 2.
        cdata(1+n1_2/2) = cdata(1+n1_2/2) / 2.
        cdata(n1_2/2+2:)=cmplx(0,0.)
        do i1=1,n1_2
          ctop(i1) = cmplx(0., 2.*3.14159*(i1-1.)/n1_2) * cdata(i1)
        end do
        call ftu( -1., n1_2, cdata)
        call ftu( -1., n1_2, ctop)
        bot(:,i2)=conjg(cdata)*cdata
        top(:,i2)=conjg(cdata)*ctop/cmplx(0.,1.)
      end do
      call triangle2( tc1, tc2, n1_2, n2_2, bot, bot)
      call triangle2( tc1, tc2, n1_2, n2_2, top, top)
      do i2=1,n2_2
        do i1=1,n1_2
          if (bot(i1,i2)>0.) then
            top(i1,i2)=top(i1,i2)/bot(i1,i2)
          else
            top(i1,i2)=0.
          end if
        end do
      end do
      image(:,:,i3)=sqrt(bot(1:size(image),1:size(image,2)))
    end do
    deallocate(data,top,bot,cdata,ctop)
    calc_envelope=.true.
  end function
  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( :), vv( :)
    real pp(n12+nr-1), qq(n12+nr+nr-2), tt(n12)
    qq=uu(1:n12)
    tt=qq
    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)
    end do
    do i= 1, nr-1
    ! fold back near end
      tt(i) = tt(i) + qq(nr-i)
    end do
    do i= 1, nr-1
    ! fold back far end
      tt(n12-i+1) = tt(n12-i+1) + qq(n12+(nr-1)+i)
    end do
    vv=tt
  end subroutine
  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(:), yy(:)
    real :: bb(nx+nb)
    if ( nb < 1 .or. nb > nx) then
      call erexit('boxconv') ! "||" means .OR.
    end if
    ny = nx+nb-1
    do i= 1, ny
      bb(i) = 0.
    end do
    bb(1) = xx(1)
    do i= 2, nx
      bb(i) = bb(i-1) + xx(i) ! make B(Z) = X(Z)/(1-Z)
    end do
    do i= nx+1, ny
      bb(i) = bb(i-1)
    end do
    do i= 1, nb
      yy(i) = bb(i)
    end do
    do i= nb+1, ny
      yy(i) = bb(i) - bb(i-nb) ! make Y(Z) = B(Z)*(1-Z**nb)
    end do
    do i= 1, ny
      yy(i) = yy(i) / nb
    end do
  end subroutine
  subroutine triangle2( rect1, rect2, n1, n2, uu, vv)
    integer :: i1,i2, rect1, rect2, n1, n2
    real :: uu(:,:), vv(:,:)
    real, allocatable :: one(:),two(:)
    real :: ss(n1,n2)
    allocate(one(size(uu,2)))
    allocate(two(size(uu,2)))
    do i1= 1, n1
      one=uu(i1,:)
      call triangle( rect2, 1, n2, one,two)
      ss(i1,:)=two
    end do
    deallocate(one,two)
    allocate(one(size(uu,1)))
    allocate(two(size(uu,1)))
    do i2= 1, n2
      one=ss(:,i2)
      call triangle( rect1, 1, n1,one,two)
      vv(:,i1)=two
    end do
    deallocate(one,two)
  end subroutine
end module
