# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!
!
!
! region_change
!
! 0 - complete variability
! 1 - v(z) change
! 2 - constant change
! 3 - no change
!
! region_init
! 0 - no modification;
! 1 - v(z)
! 2 - constant
!
! region_init_v0
! region_init_dz
! region_zero_at_layer
! region_init_minv
! region_init_maxv
!
! region_final_minv
! region_final_maxv
!
module rectangle_region_mod
  use adj_mod
  use param_check_mod
  use reflector_mod
  use tomo_space_conv_mod
  implicit none
  integer, private, save :: region_top(40)
  integer, private, save :: region_bot(40)
  integer, private, save,target :: region_init(40)
  real , private, save :: region_init_minv(40)
  real , private, save :: region_init_maxv(40)
  integer, private, save :: region_zero_layer(40)
  real , private, save :: region_init_dz(40)
  real , private, save :: region_init_v0(40)
  real , private, save :: region_final_minv(40)
  real , private, save :: region_final_maxv(40)
  integer, private, save,target :: region_change(40)
  integer, private, pointer,save :: region_type(:)
  integer, private, save :: nreg,n_mod
  integer, private, pointer, save :: itop(:,:,:)
  integer, private, pointer, save :: regions(:,:,:)
  integer, private, save :: ns(3),ns1,ns2,ns3
  real , private, save :: os(3),ds(3)
  integer, private, save :: region_ibeg(40),region_ndat(40)&
    &,region_nmod(40)
  logical, private, save :: write_regions
  contains
  logical function init_rectangle_regions(nmod)
    integer :: nmod
    real, pointer :: vel3d(:,:,:)
    type(sep3d) :: vel_s
    init_rectangle_regions=.false.
    if (.not. vel_tomo_structure(vel_s)) then
      write(0,*) "trouble grabbing velocity structure"
      return
    end if
    call set_n_o_d_3d(vel_s,ns,os,ds)
    if (.not. vel_tomo_3d_pt(vel3d)) then
      write(0,*) "trouble grabbing velocity"
      return
    end if
    if (.not. calc_regions(vel3d)) then
      write(0,*) "trouble calculating regions"
      return
    end if
    if (.not. region_calc_ndat()) then
      write(0,*) "trouble calculating the size of various regions"
      return
    end if
    if (.not. region_init_ds(vel3d)) then
      write(0,*) "trouble calculating initial ds"
      return
    end if
    if (.not. region_calc_npts(nmod,region_change)) then
      write(0,*) "trouble calculating number of pts"
      return
    end if
    if (.not. write_out_field("regions",vel_s,real(regions))) then
      write(0,*) "trouble writing out regions"
      return
    end if
    init_rectangle_regions=.true.
  end function
  logical function region_init_ds(vel3d)
    real :: vel3d(:,:,:)
    real :: temp(ns1,ns2,ns3)
    real :: acumx(40),acumy(40),acumxy(40),acumxx(40),acumyy(40)
    integer :: ireg
    real :: xy,yy,xx,x,y,b
    integer :: iloc(40),ns,i,i1,i2,i3,nsz
    real, pointer :: mod(:)
    temp=vel3d
    acumxx=0
    acumxy=0
    acumyy=0
    acumx=0
    acumy=0
    region_init_ds=.false.
!CALCULATE VZ ACCUMULATORS
    do i3=1,ns3
      do i2=1,ns2
        do i1=1,ns1
          ireg=regions(i1,i2,i3)
          if (temp(i1,i2,i3) > region_init_maxv(ireg)) then
            temp(i1,i2,i3)=region_init_maxv(ireg)
          end if
          if (temp(i1,i2,i3) < region_init_minv(ireg)) then
            temp(i1,i2,i3)=region_init_minv(ireg)
          end if
          acumx(ireg)= acumx(ireg)+(i1 -itop(i1,i2,i3))
          acumxx(ireg)= acumxx(ireg)+(i1 -itop(i1,i2,i3))**2
          acumy(ireg)= acumy(ireg)+vel3d(i1,i2,i3)
          acumxy(ireg)= acumxy(ireg)+(i1 -itop(i1,i2,i3))*vel3d(i1&
            &,i2,i3)
          acumyy(ireg)= acumyy(ireg)+vel3d(i1,i2,i3)**2
        end do
      end do
    end do
    acumx=acumx*ds(1)
    acumxy=acumx*ds(1)
    acumxx=acumx*ds(1)*ds(1)
!CALCULATE V0 AND DZ
    do ireg=1,nreg
      b=(region_ndat(ireg)*acumxy(ireg)-acumx(ireg)*acumy(ireg))/ &
        & region_ndat(ireg)*acumxx(ireg)-acumx(ireg)*acumx(ireg)
      if (region_init(ireg).eq.2) then
        b=0
      end if
      if (region_init_v0(ireg).eq. default_real) then
        region_init_v0=(acumy(ireg)-b*acumx(ireg))/region_ndat(ireg)
        region_init_dz=b
      end if
    end do
    region_init_dz=region_init_dz*ds(1)
    if (.not. region_calc_npts(nsz,region_init)) then
      write(0,*) "trouble calculating number of pts"
      return
    end if
    i=0
    allocate(mod(nsz))
    iloc=region_ibeg
!MAP TO THE MODEL SPACE
    do i3=1,ns3
      do i2=1,ns2
        do i1=1,ns1
          ireg=regions(i1,i2,i3)
          i= i+1
          if (region_init(ireg) .eq.0) then
            iloc(ireg)= iloc(ireg)+1
            mod(iloc(ireg))=temp(i1,i2,i3)
          else if (region_init(ireg).eq.1 .and. iloc(ireg) .eq.region_ibeg&
            &(ireg)) then
            mod(iloc(ireg)+1)=region_init_v0(ireg)
            mod(iloc(ireg)+2)=region_init_dz(ireg)
            iloc(ireg)= iloc(ireg)+2
          else if (region_init(ireg).eq.2 .and. iloc(ireg) .eq.region_ibeg&
            &(ireg)) then
            mod(iloc(ireg)+1)=region_init_v0(ireg)
            iloc(ireg)= iloc(ireg)+1
          end if
        end do
      end do
    end do
    n_mod=nsz
!CALL THE FORWARD OPERATOR
    region_type=>region_init
    i=region_des_3d_op(.false.,mod,temp)
    if (.not. add_vel_change(temp-vel3d)) then
      write(0,*) "troublad adding velocity change"
      return
    end if
    deallocate(mod)
    region_type=>region_change
    region_init_ds=.true.
  end function
  integer function rectangle_regions_inv_op(adj,add,model,data)
    logical, intent(in) :: adj,add
    real :: model(:),data(:)
    integer :: ireg,i1,i2,i3,i
    integer :: iloc(40)
    call adjnull(adj,add,model,data)
    i=0
    iloc=region_ibeg !
    if (.not. adj) then
      do i3=1,ns3
        do i2=1,ns2
          do i1=1,ns1
            i= i+1
            ireg=regions(i1,i2,i3)
            if (region_type(ireg).eq.0) then
              iloc(ireg)= iloc(ireg)+1
              data(iloc(ireg))= data(iloc(ireg))+model(i)
            end if
          end do
        end do
      end do
    else
      do i3=1,ns3
        do i2=1,ns2
          do i1=1,ns1
            i= i+1
            ireg=regions(i1,i2,i3)
            if (region_type(ireg).eq.0) then
              iloc(ireg)= iloc(ireg)+1
              model(i)= model(i)+data(iloc(ireg))
            end if
          end do
        end do
      end do
    end if
    rectangle_regions_inv_op=0
  end function
  integer function region_des_3d_op(adj,model,data)
    real :: model(n_mod),data(ns1,ns2,ns3)
    logical, intent(in) :: adj
    integer :: i1,i2,i3
    integer :: iloc(40),ireg
    region_des_3d_op=0
    iloc=region_ibeg
    if (.not. adj) then
      do i3=1,ns3
        do i2=1,ns2
          do i1=1,ns1
            ireg=regions(i1,i2,i3)
            if (ireg.eq.0) then
              iloc(ireg)= iloc(ireg)+1
              data(i1,i2,i3)= data(i1,i2,i3)+model(iloc(ireg))
            else if (ireg.eq.1) then
              data(i1,i2,i3)= data(i1,i2,i3)+model(iloc(ireg)&
                &+1)+(i1-itop(i2,i3,ireg))*model(iloc(ireg+2))
            else if (ireg.eq.2) then
              data(i1,i2,i3)= data(i1,i2,i3)+model(iloc(ireg)&
                &+1)
            end if
          end do
        end do
      end do
    else
      do i3=1,ns3
        do i2=1,ns2
          do i1=1,ns1
            ireg=regions(i1,i2,i3)
            if (ireg.eq.0) then
              iloc(ireg)= iloc(ireg)+1
              model(iloc(ireg))= model(iloc(ireg))+data(i1,i2&
                &,i3)
            else if (ireg.eq.1) then
              model(iloc(ireg)+1)= model(iloc(ireg)+1)+data(i1&
                &,i2,i3)
              model(iloc(ireg+2))= model(iloc(ireg+2))+data(i1&
                &,i2,i3)*(i1-itop(i2,i3,ireg))
            else if (ireg.eq.2) then
              model(iloc(ireg)+1)= model(iloc(ireg)+1)+data(i1&
                &,i2,i3)
            end if
          end do
        end do
      end do
    end if
  end function
  logical function rect_region_clean(vel)
    real :: vel(:,:,:)
    integer :: i1,i2,i3,ireg
    rect_region_clean=.false.
    do i3=1,ns3
      do i2=1,ns2
        do i1=1,ns1
          ireg=regions(i1,i2,i3)
          vel(i1,i2,i3)=max(region_final_maxv(ireg),min(region_final_minv&
            &(ireg),vel(i1,i2,i3)))
        end do
      end do
    end do
    rect_region_clean=.true.
  end function
  logical function region_calc_ndat()
    integer :: i1,i2,i3
    region_calc_ndat=.false.
    if (any(regions <1) .or. any(regions > nreg)) then
      write(0,*) "region range doesn't make sense should be",1,nreg
      write(0,*) "region range is",minval(regions),maxval(regions)
      return
    end if
    region_ndat=0
    do i3=1,ns3
      do i2=1,ns2
        do i1=1,ns1
          region_ndat(regions(i1,i2,i3))= region_ndat(regions(i1&
            &,i2,i3))+1
        end do
      end do
    end do
    region_calc_ndat=.true.
  end function
  logical function region_calc_npts(nmod,itype)
    integer :: i,itype(:),nmod
    region_calc_npts=.false.
    nmod=0
    do i=1,nreg
      if (itype(i) .eq. 0) then
        region_nmod(i)=region_ndat(i)
      else if (itype(i).eq.1) then
        region_nmod(i)=2
      else if (itype(i).eq.2) then
        region_nmod(i)=1
      else if (itype(i).eq.3) then
        region_nmod(i)=0
      end if
    end do
    region_ibeg(1)=0
    do i=2,nreg
      region_ibeg(i)=region_ibeg(i-1)+region_nmod(i)
    end do
    n_mod=nmod
    region_calc_npts=.true.
  end function
  logical function calc_regions(vel3d)
    type(sep3d) :: ref_s,vel_s
    complex, pointer :: ref2d(:,:)
    real, pointer :: ref3d(:,:,:)
    real, pointer :: vel3d(:,:,:)
    integer :: ireg,i
    calc_regions=.false.
    if (.not. return_ref(ref_s)) then
      write(0,*) "trouble grabbing reflector structure"
      return
    end if
    if (ns3.eq.1) then
      allocate(ref2d(ref_s%n(1),ref_s%n(2)))
      if (.not. return_2d_conv_refs(ref2d,.true.)) then
        write(0,*) "trouble grabbing 2-d converted refs"
        return
      end if
      nullify(ref3d)
    else
      allocate(ref3d(ref_s%n(1),ref_s%n(2),ref_s%n(3)))
      if (.not. return_3d_conv_refs(ref3d,.true.)) then
        write(0,*) "trouble grabbing 3-d converted refs"
        return
      end if
      if (.not. extend_3d_ref(ref_s,ref3d)) then
        write(0,*) "trouble extending reflector"
        return
      end if
      nullify(ref2d)
    end if
    allocate(regions(ns1,ns2,ns3))
    call add_mem(size(regions))
    allocate(itop(ns1,ns2,nreg))
    itop=-1
    call add_mem(size(itop))
    regions=1
    do ireg=2,nreg
      if (.not. add_region(ireg,itop(:,:,ireg),ref2d,ref3d)) then
        write(0,*) "trouble adding region",ireg
      end if
    end do
    if (associated(ref2d)) then
      deallocate(ref2d)
    end if
    if (associated(ref3d)) then
      deallocate(ref3d)
    end if
    calc_regions=.true.
  end function
  logical function add_region(ireg,itop,ref2d,ref3d)
    integer :: ireg
    integer :: itop(:,:)
    complex, pointer :: ref2d(:,:)
    real, pointer :: ref3d(:,:,:)
    integer :: map(ns1,ns2,ns3)
    integer :: i2,i3,i1,i,ib,ia
    add_region=.false.
    if (region_top(ireg).eq.0) then
      map(1,:,:)=0
    else if (associated(ref2d)) then
      if (.not. map_2d_ref(ref2d(:,region_top(ireg)),map,ireg)) then
        write(0,*) "trouble mapping 2-D reflector"
        return
      end if
    else
      if (.not. map_3d_ref(ref3d(:,:,region_top(ireg)),map,ireg)) then
        write(0,*) "trouble mapping 3-D reflector"
        return
      end if
    end if
    if (region_bot(ireg).eq.0) then
      map(ns1,:,:)=0
    else if (associated(ref2d)) then
      if (.not. map_2d_ref(ref2d(:,region_bot(ireg)),map,-ireg)) then
        write(0,*) "trouble mapping 2-D reflector"
        return
      end if
    else
      if (.not. map_3d_ref(ref3d(:,:,region_bot(ireg)),map,-ireg))&
        & then
        write(0,*) "trouble mapping 3-D reflector"
        return
      end if
    end if
    do i3=1,size(map,3)
      do i2=1,size(map,2)
        ia=-1
        do i1=1,size(map,1)
          if (map(i1,i2,i3).eq.ireg) then
            ia=i1
            itop(i2,i3)=ia
          end if
          if (map(i1,i2,i3).eq.-ireg) then
            if (ia.ne.-1) then
              regions(ia:i1,i2,i3)=ireg
            end if
          end if
        end do
      end do
    end do
    if (region_zero_layer(ireg).eq.0) then
      itop=0
    end if
    add_region=.true.
  end function
  logical function extend_3d_ref(ref_s,ref3d)
    type(sep3d) :: ref_s
    real :: o(3),d(3)
    real, pointer :: ref3d(:,:,:),temp(:,:,:)
    integer :: n(3)
    integer :: i2,i3
    real :: f2,f3
    real :: x,y
    real :: fx,fy
    integer :: ix,iy
    extend_3d_ref=.false.
    call set_n_o_d_3d(ref_s,n,o,d)
    allocate(temp(n(1),n(2),n(3)))
    temp=ref3d
    deallocate(ref3d)
    allocate(ref3d(ns1,ns2,ns3))
    ref3d=os(1)-ds(1)
    do i3=1,ns3
      y=os(3)+ds(3)*(i3-1)
      fy=(y-o(2))/d(2)+1.
      iy=floor(fy)
      if (iy>0 .and. iy < n(2)) then
        do i2=1,ns2
          x=os(2)+ds(2)*(i2-1)
          fx=(x-o(1))/d(1)+1.
          ix=floor(fx)
          if (ix>0 .and. ix < n(2)) then
            ref3d(i2,i3,:)= (1.-fx)*(1.-fy)*temp(ix ,iy ,:) +( &
              & fx)*(1.-fy)*temp(ix+1 ,iy ,:) +(1.-fx)*( fy)*temp&
              &(ix ,iy+1 ,:) +( fx)*( fy)*temp(ix+1 ,iy+1 ,:)
          end if
        end do
      end if
    end do
    deallocate(temp)
    extend_3d_ref=.true.
  end function
  logical function map_3d_ref(ref3d,map,ival)
    real :: ref3d(:,:)
    integer :: ival
    integer :: map(:,:,:)
    real :: fx,fy,z
    integer :: i,ixb,iz,ixe,ilastx,iyb,iye,ilasty,i1,i2,i3
    logical :: found
    map_3d_ref=.false.
    ref3d=(ref3d-os(1))/ds(1)+1.
    do i3=1,ns3
      do i2=1,ns2
        iz=nint(ref3d(i2,i3))
        if (iz >0 .and. iz <= ns1) then
          map(iz,i2,i3)=ival
        end if
      end do
    end do
    map_3d_ref=.true.
  end function
  logical function map_2d_ref(ref2d,map,ival)
    complex :: ref2d(:)
    integer :: ival
    integer :: map(:,:,:)
    real :: f,z
    integer :: i,ixb,iz,ixe,ilast
    map_2d_ref=.false.
    map_2d_ref=.false.
    do i=1,size(ref2d)
      ref2d(i)=cmplx((real(ref2d(i))-os(2))/ds(2)+1.,(aimag(ref2d(i)-os&
        &(1))/ds(1)+1.))
    end do
!ref2d=cmplx((real(ref2d)-os(2))/ds(2)+1.,(aimag(ref2d)-os(1))/ds(1)+1.)
    ixb=nint(real(ref2d(1)))
    ixe=nint(real(ref2d(size(ref2d))))
    ilast=1
    do i=ixb,ixe
      do while (ilast < ixe .and. i > real(ref2d(ilast)))
        ilast= ilast+1
      end do
      ilast= ilast-1
      if (ilast.eq.0) then
        ilast=1
        f =0.
        z=aimag(ref2d(1))
      else
        f=(real(i)-real(ref2d(ilast)))/(real(ref2d(ilast+1)-real(ref2d&
          &(ilast))))
        z=(1.-f)*aimag(ref2d(ilast))+f*aimag(ref2d(ilast+1))
      end if
      iz=nint(z)
      if (iz >0 .and. iz <=ns1) then
        map(iz,i,:)=ival
      end if
    end do
    map_2d_ref=.true.
  end function
  integer function rectangle_regions_mask_op(adj,add,model,data)
    logical,intent(in) :: adj,add
    real :: model(:),data(:)
    call adjnull(adj,add,model,data)
    if (.not. adj) then
      data= data+model
    else
      model= model+data
    end if
    rectangle_regions_mask_op=0
  end function
  integer function rectangle_regions_des_op(adj,add,model,data)
    logical,intent(in) :: adj,add
    real :: model(:),data(:)
    call adjnull(adj,add,model,data)
    rectangle_regions_des_op=region_des_3d_op(adj,model,data)
    rectangle_regions_des_op=0
  end function
  logical function return_no_smooth_rectangle_regions(region,nosmooth)
    integer :: region(:)
    logical :: nosmooth(:)
    integer :: i1,i2,i3,i
    return_no_smooth_rectangle_regions=.false.
    if (size(nosmooth) .ne. size(region)) then
      write(0,*) "region and nosmooth not the same size"
      return
    end if
    nosmooth=reshape(regions,(/size(regions)/))
    i=0
    do i3=1,ns3
      do i2=1,ns2
        do i1=1,ns1
          i= i+1
          if (region_type(regions(i1,i2,i3).eq.0)) then
            nosmooth(i)=.false.
          else
            nosmooth(i)=.true.
          end if
        end do
      end do
    end do
    return_no_smooth_rectangle_regions=.true.
  end function
  logical function grab_pars_rectangle_regions()
    type(sep3d) :: ref_s
    integer :: i,nref,ireg
    grab_pars_rectangle_regions=.false.
    if (.not. return_ref(ref_s)) then
      write(0,*) "trouble grabbing reflector structure"
      return
    end if
    nref=ref_s%n(sep3d_ndims(ref_s))
    region_top=default_int
    region_bot=default_int
    region_zero_layer=0
    region_change=0
    region_init=0
    region_init_minv=.1
    region_init_maxv=999999.
    region_init_dz=default_real
    region_init_v0=default_real
    region_final_minv=.1
    region_final_maxv=99999.
    call from_param("region_top",region_top,region_top)
    call from_param("region_bot",region_bot,region_bot)
    call from_param("region_zero_layer",region_zero_layer,region_zero_layer&
      &)
    call from_param("region_change",region_change,region_change)
    call from_param("region_init_minv",region_init_minv,region_init_minv&
      &)
    call from_param("region_init_maxv",region_init_maxv,region_init_maxv&
      &)
    call from_param("region_init_dz",region_init_dz,region_init_dz)
    call from_param("region_init_v0",region_init_v0,region_init_v0)
    call from_param("region_final_minv",region_final_minv,region_final_minv&
      &)
    call from_param("region_final_maxv",region_final_maxv,region_final_maxv&
      &)
    nreg=1
    do i=2,40
      if (region_top(i) .ne. default_int) then
        if (region_bot(i) .ne. default_int) then
          if (region_top(i) > -1 .and. region_top(i) <= nref) then
            write(0,*) "region_top",i," must be between 0 and",nref,"&
              & not ",region_top(i)
            return
          end if
          if (region_bot(i) > -1 .and. region_bot(i) <= nref) then
            write(0,*) "region_bot",i," must be between 0 and",nref,"&
              & not ",region_bot(i)
            return
          end if
          nreg= nreg+1
        else
          write(0,*) "region_top",i," specified region_bot not"
          return
        end if
      else
        write(0,*) "region_top",i," specified region_bot not"
        return
      end if
    end do
    if (region_change(1) .ne. 0 .or. region_init(1).ne.0) then
      write(0,*) "the first region needs to be variable"
      write(0,*) "if you want vz,etc specify top and bottom as region&
        & 2"
      return
    end if
    do ireg=2,nreg
      if ((region_init_v0(ireg) .eq. default_real .and. region_init_dz&
        &(ireg) .ne. default_real ) .or. (region_init_v0(ireg) .eq.&
        & default_real .and. region_init_dz(ireg) .ne. default_real )) then
        write(0,*) "must specify both or nether region_init_v0 and&
          & region_init_dz region=",ireg
        return
      end if
    end do
    if (any(region_init(2:nreg) < 0) .or. any(region_init(2:nreg) >2)&
      &) then
      write(0,*) "region_init must be between 0 and 2",minval(region_init&
        &(2:nreg)),maxval(region_init(2:nreg))
      return
    end if
    if (any(region_change(2:nreg) < 0) .or. any(region_change(2:nreg) &
      & >3)) then
      write(0,*) "region_change must be between 0 and 3",minval(region_change&
        &(2:nreg)),maxval(region_change(2:nreg))
      return
    end if
    write_regions=check_write("regions")
    if (write_regions) then
      call setup_reg_file("regions",ns,os,ds,1)
    end if
    grab_pars_rectangle_regions=.true.
  end function
end module
