# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module max_pick_mod
  use sep3d_struct_mod
  use sep
  use envelope_mod
  use quantile_mod
  implicit none
  type(sep3d),private :: d
  integer,pointer, dimension(:,:),private :: lbound,ubound
  real,pointer, private :: data(:,:,:),dip(:,:)
  real, private, save :: stack_rat,pow
  real, pointer,private :: l_bound(:,:),u_bound(:,:),stack(:,:,:),env(:&
    &,:,:),edge(:)
  contains
  logical function init_max_params(data_s,dat,st,mp,wt)
    type(sep3d) :: data_s,bound
    real,pointer :: mp(:,:),wt(:,:),dat(:,:,:),st(:,:,:),temp(:)
    real, allocatable :: t(:,:)
    real :: u_bound_top,u_bound_bot,myt
    real :: l_bound_top,l_bound_bot,f,myp,pow_dip,beg0,beg1,end0,end1
    integer, external :: getch,sreed
    integer :: i,tc1,tc2
    init_max_params=.false.
    d=data_s
    allocate(data(size(dat,1),size(dat,2),size(dat,3)))
    data=dat
    allocate(dip(data_s%n(1),data_s%n(2)))
    allocate(wt(data_s%n(1),data_s%n(2)))
    allocate(env(size(data,1),size(data,2),1))
    i= sreed("image",env,size(env)*4)
    call from_param("env_smooth1",tc1,8)
    call from_param("pow",pow,2.)
    call from_param("env_smooth2",tc2,3)
    if (.not. calc_envelope(env,tc1,tc2)) then
      call seperr("trouble calculating evelope")
    end if
    allocate(edge(data_s%n(1)))
    call from_param("beg0",beg0,data_s%o(1)-data_s%d(1))
    call from_param("beg1",beg1,data_s%o(1))
    call from_param("end0",end0,data_s%o(1)+data_s%n(1) *data_s%d(1))
    call from_param("end0",end1,data_s%o(1)+data_s%n(1) *data_s%d(1)&
      &-data_s%d(1))
    edge=1
    do i=1,size(edge)
      myt=data_s%o(1)+data_s%d(1)*(i-1)
      if (myt < beg0) then
        edge(i)=0.
      else if (myt > end0) then
        edge(i)=0.
      end if
    end do
    dip=1.
    if (1.eq.getch("dip","d",i)) then
      i= sreed("dip",dip,size(dip)*4)
      dip=atan(dip) !dip is now angle in radians
      call from_param("pow_dip",pow_dip,2.)
      dip=cos(dip)**pow_dip
    end if
    allocate(temp(size(env)))
    temp= reshape(env,(/size(env)/))
!myp=quantile(nint(size(env)*.98),temp)
    myp=maxval(env)*.94
    deallocate(temp)
    env=env/myp
    allocate(mp(data_s%n(1),data_s%n(2)))
    allocate(l_bound(data_s%n(1),d%n(2)),u_bound(data_s%n(1),d%n(2)))
    allocate(lbound(data_s%n(1),d%n(2)),ubound(data_s%n(1),d%n(2)))
    if (1.eq.getch("bounds","d",i)) then
      call init_sep3d("bounds",bound,"INPUT")
      if (bound%n(1).ne.d%n(1) .and. bound%n(2).ne.2) then
        write(0,*) "expecting bounds to be nt,2"
        return
      end if
      allocate(t(bound%n(1),2))
      if (bound%n(1)*bound%n(2)*4.ne.sreed("bounds",t,size(t)*4)) then
        write(0,*) "trouble reading bounds"
        return
      end if
      l_bound(:,1)=t(:,1)
      u_bound(:,1)=t(:,2)
      deallocate(t)
    else
      call from_param("l_bound_top",l_bound_top,data_s%o(4))
      call from_param("l_bound_bot",l_bound_bot,data_s%o(4))
      call from_param("u_bound_top",u_bound_top,data_s%o(4)+data_s%d(4&
        &)*(data_s%n(4)-1))
      call from_param("u_bound_bot",u_bound_bot,data_s%o(4)+data_s%d(4&
        &)*(data_s%n(4)-1))
      do i=1,data_s%n(1)
        f=real(i-1)/real(data_s%n(1))
        l_bound(i,1)=(1.-f)*l_bound_top+f*l_bound_bot
        u_bound(i,1)=(1.-f)*u_bound_top+f*u_bound_bot
      end do
      do i=1,data_s%n(2)
        u_bound(:,i)=u_bound(:,1)
        l_bound(:,i)=l_bound(:,1)
      end do
    end if
    call from_param("env_rat",stack_rat,.5)
    init_max_params=.true.
  end function
  logical function max_pick(mp,wt)
    real,pointer :: mp(:,:),wt(:,:)
    integer :: i1,i2,imax1(1),imax2(1)
    real :: my_max
    max_pick=.false.
    do i2=1,d%n(2)
      do i1=1,d%n(1)
        lbound(i1,i2)=max(1,min(d%n(4),nint((l_bound(i1,i2)-d%o(4))/d%d&
          &(4)+1.)))
        ubound(i1,i2)=max(1,min(d%n(4),nint((u_bound(i1,i2)-d%o(4))/d%d&
          &(4)+1.)))
      end do
    end do
    do i2=1,size(data,2)
      do i1=1,size(data,1)
        imax1=maxloc(data(i1,i2,:))
        imax2=maxloc(data(i1,i2,lbound(i1,i2):ubound(i1,i2)))+lbound(i1&
          &,i2)-1
        wt(i1,i2)=data(i1,i2,imax2(1))*env(i1,i2,1)*edge(i1)
        if (imax1(1).ne.imax2(1)) then
          wt(i1,i2)=wt(i1,i2)/2.
        end if
        mp(i1,i2)=(imax2(1)-1)*d%d(4)+d%o(4)
      end do
    end do
    wt=wt**pow*dip
    max_pick=.true.
  end function
  logical function reset_bounds(base,delta)
    real :: base(:),delta
    integer :: i1,i2,i
    reset_bounds=.false.
    i=0
    do i2=1,size(lbound,2)
      do i1=1,size(lbound,1)
        i= i+1
        l_bound(i1,i2)=max(l_bound(i1,i2),base(i)-delta)
        u_bound(i1,i2)=min(u_bound(i1,i2),base(i)+delta)
      end do
    end do
    reset_bounds=.true.
  end function
end module
