# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!<
!
! Recoginized types :: z2d,z3d,tau2d,tau3d
!
!
!
!
! MODULE CONTAINING CONVERSION FROM DEPTH
! TO WHAT EVER SPACE TOMOGRAPHY WILL BE PERFORMED IN
! TO ADD TYPE MAKE SURE TO INCLUDE IN CONV_TYPE
! AND ALL ACCESSING FUNCTIONS.
!
!
!>
module tomo_space_conv_mod
  use sep
  use velocity_mod
  use tau_z_conv_3d_mod
  use put_3d_mod
  implicit none
  real, pointer,private :: velocity_tomo(:,:,:),sigma_x(:,:,:),sigma_y&
    &(:,:,:),change(:)
  real, pointer, private :: mapz(:,:,:),mapt(:,:,:)
  character(len=10),save :: model_space
  integer,private,save :: beg(3),end(3),ntau
  real, private,save :: otau,dtau
  type(sep3d),private ,save:: vel,vel_tomo_struct,vel_z_tomo_struct
  logical, private, save :: tau
  interface add_change
    module procedure add_change_1d
    module procedure add_change_3d
  end interface
  contains
  logical function init_space_conv()
    integer :: ierr
    real, pointer :: tempr(:,:,:)
    real, pointer :: tempr2(:,:,:)
    real :: o(3),d(3)
    integer :: n(3)
    init_space_conv=.false.
    if (.not. vel_structure(vel)) then
      write(0,*) "vel space not initialized"
      return
    end if
    call set_n_o_d_3d(vel,n,o,d)
    allocate(tempr(n(1),n(2),n(3)))
    if (.not. vel_3d(tempr)) then
      write(0,*) "    trouble grabbing velocity function"
      return
    end if
    if (.not. grab_space_conv_pars()) then
      write(0,*) "    trouble grabbing space convert parameters"
      return
    end if
    if (.not. create_z_pad_s(vel,vel_z_tomo_struct)) then
      write(0,*) "    trouble creating padded velocity structure"
      return
    end if
    if (.not. pad_vel(tempr)) then
      write(0,*) "    trouble padding velocity velocity"
      return
    end if
    call init_sep3d(vel_z_tomo_struct,vel_tomo_struct,"SCRATCH","vel_tomo"&
      &)
    if (tau) then
      vel_tomo_struct%n(1)=ntau
      vel_tomo_struct%o(1)=otau
      vel_tomo_struct%d(1)=dtau
      if (.not. create_sigma_maps(tempr,sigma_x,sigma_y)) then
        write(0,*) "    trouble creating sigma maps"
        return
      end if
      allocate(velocity_tomo(ntau,size(tempr,2),size(tempr,3)))
      if (.not. convert_field_z_tau_3d(tempr,velocity_tomo)) then
        write(0,*) "    trouble converting velocity from z  to tau"
        return
      end if
    else
      allocate(velocity_tomo(size(tempr,1),size(tempr,2),size(tempr,3)&
        &))
      velocity_tomo=tempr
    end if
    allocate(change(size(velocity_tomo)))
    change=0.
    deallocate(tempr)
    init_space_conv=.true.
  end function
  logical function create_sigma_maps(velin,sigmax,sigmay)
    real :: velin(:,:,:)
    real,pointer :: sigmax(:,:,:)
    real,pointer :: sigmay(:,:,:)
    real, pointer :: temp1(:,:,:),temp2(:,:,:)
    integer :: nz(3),nt(3)
    real :: ot(3),oz(3)
    real :: dt(3) ,dz(3)
    create_sigma_maps=.false.
    call set_n_o_d_3d(vel_tomo_struct,nt,ot,dt)
    call set_n_o_d_3d(vel_z_tomo_struct,nz,oz,dz)
    if (.not. init_tau_z_conv_3d(nz,oz,dz,nt,ot,dt)) then
      write(0,*) "trouble initializing tau-z converter"
      return
    end if
    allocate(mapz(nz(1),nz(2),nz(3)))
    call add_mem(size(mapz))
    if (.not. create_tau_map_3d(velin,mapz)) then
      write(0,*) "trouble createing tau map"
      return
    end if
    if (.not. set_z_vel_map(velin,mapz)) then
      write(0,*) "trouble setting z vel map"
      return
    end if
    allocate(temp1(nz(1),nz(2),nz(3)))
    if (.not. calc_sigmax_3d(velin,temp1)) then
      write(0,*) "trouble calculating sigmax"
      return
    end if
    allocate(sigmax(nt(1),nt(2),nt(3)))
    call add_mem(size(sigmax))
    if (.not. convert_field_z_tau_3d(temp1,sigmax)) then
      write(0,*) "trouble converting sigmax to tau space"
      return
    end if
    if (nz(3)>1) then
      if (.not. calc_sigmay_3d(velin,temp1)) then
        write(0,*) "trouble calculating sigmax"
        return
      end if
      allocate(sigmay(nt(1),nt(2),nt(3)))
      call add_mem(size(sigmay))
      if (.not. convert_field_z_tau_3d(temp1,sigmay)) then
        write(0,*) "trouble converting sigmax to tau space"
        return
      end if
    end if
    create_sigma_maps=.true.
  end function
  logical function create_tau_struct()
    create_tau_struct=.true.
  end function
  logical function pad_vel(array)
    real, pointer :: array(:,:,:)
    real, pointer :: temp(:,:,:)
    integer :: ns(3),nb(3),i1,i2,i3
    real :: o(3),d(3)
    pad_vel=.false.
    if (any(beg.ne.0) .or. any(end.ne.0)) then
      call set_n_o_d_3d(vel,ns,o,d)
      call set_n_o_d_3d(vel_z_tomo_struct,nb,o,d)
      allocate(temp(ns(1),ns(2),ns(3)))
      temp=array
      deallocate(array)
      allocate(array(nb(1),nb(2),nb(3)))
      do i3=1,ns(3)
        do i2=1,ns(2)
          call pad_axis(beg(1),end(1),temp(:,i2,i3),array(:,i2+beg(2)&
            &,i3+beg(3)))
        end do
      end do
      deallocate(temp)
      do i3=1,size(array,3)
        do i1=1,size(array,1)
          call pad_axis2(beg(2),end(2),array(i1,:,i3))
        end do
      end do
      do i2=1,size(array,2)
        do i1=1,size(array,1)
          call pad_axis2(beg(3),end(3),array(i1,i2,:))
        end do
      end do
    end if
    pad_vel=.true.
  end function
  logical function create_z_pad_s(vel_s,velz_pad_s)
    type(sep3d) :: vel_s,velz_pad_s
    integer :: ndim
    create_z_pad_s=.false.
    ndim=sep3d_ndims(vel_s)
    call init_sep3d(vel_s,velz_pad_s,"SCRATCH","velz_pad_s")
    velz_pad_s%o(1:ndim)=vel_s%o(1:ndim)-beg(1:ndim)*vel_s%d(1:ndim)
    velz_pad_s%n(1:ndim)=beg(1:ndim)+end(1:ndim)+vel_s%n(1:ndim)
    call sep3d_set_sep3d(velz_pad_s)
    create_z_pad_s=.true.
  end function
  logical function convert_field_tomo_space(field)
    real, pointer :: field(:,:,:)
    integer :: nz(3),nt(3)
    real :: o(3),d(3)
!ASSUME THAT THE SIZE CHECKS HAVE ALREADY BEEN MAPE
    convert_field_tomo_space=.false.
    if (.not. pad_vel(field)) then
      write(0,*) "    ERROR: padding velocity velocity"
      return
    end if
    if (.not. field_z_to_tau(field)) then
      write(0,*) '    ERROR: converting field from z to tau '
      return
    end if
    convert_field_tomo_space=.true.
  end function
  logical function field_z_to_tau(field)
    real, pointer :: field(:,:,:)
    integer :: nz(3),nt(3)
    real :: o(3),d(3)
    real, pointer :: map(:,:,:),temp(:,:,:)
    field_z_to_tau=.false.
    call set_n_o_d_3d(vel_z_tomo_struct,nz,o,d)
    call set_n_o_d_3d(vel_tomo_struct,nt,o,d)
    allocate(temp(nz(1),nz(2),nz(3)))
    temp=field
    deallocate(field)
    allocate(field(nt(1),nt(2),nt(3)))
    if (.not. convert_field_z_tau_3d(temp,field)) then
      write(0,*) "    trouble converting velocity from z  to tau"
      return
    end if
    deallocate(temp)
    field_z_to_tau=.true.
  end function
  logical function convert_vel_tomo_space(field)
    real, pointer :: field(:,:,:)
    integer :: nz(3),nt(3)
    real :: o(3),d(3)
    real, pointer :: map(:,:,:),temp(:,:,:)
!ASSUME THAT THE SIZE CHECKS HAVE ALREADY BEEN MAPE
    convert_vel_tomo_space=.false.
    if (.not. pad_vel(field)) then
      write(0,*) "    trouble padding velocity velocity"
      return
    end if
    allocate(map(size(field,1),size(field,2),size(field,3)))
    if (.not. create_tau_map_3d(field,map)) then
      write(0,*) "trouble createing tau map"
      return
    end if
    allocate(temp(size(field,1),size(field,2),size(field,3)))
    temp=field
    deallocate(field)
    allocate(field(size(velocity_tomo,1),size(velocity_tomo,2),size&
      &(velocity_tomo,3)))
    if (.not. convert_field_z_tau_3d(temp,field,map)) then
      write(0,*) "    trouble converting velocity from z  to tau"
      return
    end if
    if (.not. set_z_vel_map(velocity_tomo,mapz)) then
      write(0,*) "trouble setting z_vel_map"
      return
    end if
    deallocate(temp,map)
    convert_vel_tomo_space=.true.
  end function
  logical function vel_tomo_3d_pt(myvel,mymapz,mymapt,mysigmax,mysigmay&
    &)
    real,optional,pointer :: myvel(:,:,:),mysigmax(:,:,:),mymapz(:,:,:&
      &),mymapt(:,:,:),mysigmay(:,:,:)
    integer :: n(3)
    real :: o(3),d(3)
    vel_tomo_3d_pt=.false.
    if (.not. valid_structure(vel_tomo_struct)) then
      write(0,*) "invalid vel_tomo structure"
    else
      if (present(myvel)) then
        myvel=>velocity_tomo
      end if
      if (present(mysigmax) .or. present(mymapz) .or. present(mymapt))&
        & then
        if (.not. tau) then
          write(0,*) "requested tau component when only z built"
        else
          vel_tomo_3d_pt=.true.
          if (present(mymapt)) then
            mymapt=>mapt
          end if
          if (present(mymapz)) then
            mymapz=>mapz
          end if
          if (present(mysigmax)) then
            mysigmax=>sigma_x
          end if
          if (present(mysigmay)) then
            mysigmay=>sigma_y
          end if
        end if
      else
        vel_tomo_3d_pt=.true.
      end if
    end if
  end function
  logical function vel_tomo_3d(myvel,mymapz,mymapt,mysigmax,mysigmay)
    real,optional :: myvel(:,:,:),mysigmax(:,:,:),mymapz(:,:,:),mymapt&
      &(:,:,:),mysigmay(:,:,:)
    integer :: n(3)
    real :: o(3),d(3)
    vel_tomo_3d=.false.
    if (.not. valid_structure(vel_tomo_struct)) then
      write(0,*) "invalid vel_tomo structure"
    else
      if (present(myvel)) then
        myvel=velocity_tomo
      end if
      if (present(mysigmax) .or. present(mymapz) .or. present(mymapt))&
        & then
        if (.not. tau) then
          write(0,*) "requested tau component when only z built"
        else
          vel_tomo_3d=.true.
          if (present(mymapt)) then
            mymapt=mapt
          end if
          if (present(mymapz)) then
            mymapz=mapz
          end if
          if (present(mysigmax)) then
            mysigmax=sigma_x
          end if
          if (present(mysigmay)) then
            mysigmay=sigma_y
          end if
        end if
      else
        vel_tomo_3d=.true.
      end if
    end if
  end function
  logical function vel_tomo_2d(myvel,mymapz,mymapt,mysigmax)
    real,optional :: myvel(:,:),mysigmax(:,:),mymapz(:,:),mymapt(:,:)
    vel_tomo_2d=.false.
    if (.not. valid_structure(vel_tomo_struct)) then
      write(0,*) "invalid vel_tomo structure"
    else if (sep3d_ndims(vel_tomo_struct).ne.2) then
      write(0,*) "requested 2d vel when ndims=",sep3d_ndims(vel_tomo_struct&
        &)
    else
      if (present(myvel)) then
        myvel=velocity_tomo(:,:,1)
      end if
      if (present(mysigmax) .or. present(mymapz) .or. present(mymapt))&
        & then
        if (.not. tau) then
          write(0,*) "requested tau component when only z built"
        else
          vel_tomo_2d=.true.
          if (present(mymapt)) then
            mymapt=mapt(:,:,1)
          end if
          if (present(mymapz)) then
            mymapz=mapz(:,:,1)
          end if
          if (present(mysigmax)) then
            mysigmax=sigma_x(:,:,1)
          end if
        end if
      else
        vel_tomo_2d=.true.
      end if
    end if
  end function
  logical function update_vel_tomo(myvel)
    real :: myvel(:,:,:)
    update_vel_tomo=.false.
    if (size(myvel) .ne. size(velocity_tomo)) then
      write(0,*) "updated velocity size not the same as original size"
      return
    end if
    velocity_tomo=myvel
    update_vel_tomo=.true.
  end function
  logical function vel_tomo(myvel,mysigmax,mysigmay)
    real :: myvel(:)
    real,optional ::mysigmax(:),mysigmay(:)
    vel_tomo=.false.
    if (.not. valid_structure(vel_tomo_struct)) then
      write(0,*) "invalid vel_tomo structure"
    else
      myvel=reshape(velocity_tomo,(/size(velocity_tomo)/))
      if (present(mysigmax) .or. present(mysigmay)) then
        if ( tau) then
          write(0,*) "requested tau component when only z built"
        else
          vel_tomo=.true.
          if (present(mysigmax)) then
            mysigmax=reshape(sigma_x,(/size(sigma_x)/))
          end if
          if (present(mysigmay)) then
            mysigmay=reshape(sigma_y,(/size(sigma_y)/))
          end if
        end if
      else
        vel_tomo=.true.
      end if
    end if
  end function
  logical function vel_tomo_structure(myvel)
    type(sep3d) :: myvel
    if (.not. valid_structure(vel_tomo_struct)) then
      vel_tomo_structure=.false.
    else
      myvel=vel_tomo_struct
      vel_tomo_structure=.true.
    end if
  end function
  logical function vel_z_tomo_structure(myvel)
    type(sep3d) :: myvel
    if (.not. valid_structure(vel_z_tomo_struct)) then
      vel_z_tomo_structure=.false.
    else
      myvel=vel_z_tomo_struct
      vel_z_tomo_structure=.true.
    end if
  end function
  subroutine pad_axis2(beg,end,array)
    integer :: i1,beg,end
    real :: array(:)
    if (beg>0) then
      array(1:beg)=array(beg+1)
    end if
    if (end>0) then
      array(size(array)-end+1:)=array(size(array)-end)
    end if
  end subroutine
  subroutine pad_axis(beg,end,input,output)
    integer :: i1,beg,end
    real :: input(:),output(:)
    if (beg>0) then
      output(1:beg)=input(1)
    end if
    output(beg+1:beg+size(input))=input
    if (end>0) then
      output(beg+size(input)+1:)=input(size(input))
    end if
  end subroutine
  logical function return_padding(b,e)
    integer :: b(3),e(3)
    return_padding=.true.
    b=beg
    e=end
  end function
  logical function add_vel_change(tochange)
    real :: tochange(:,:,:)
    add_vel_change=.true.
    if (model_space(1:2).eq."ds") then
      change= change+reshape( 1./(tochange+velocity_tomo) -1./velocity_tomo&
        &,(/size(velocity_tomo)/)) !OLD VEL
    else if (model_space(1:2).eq."ds") then
      change= change+reshape( (tochange+velocity_tomo)**2 -velocity_tomo&
        &**2,(/size(velocity_tomo)/)) !OLD VEL
    else
      write(0,*) "unknown model space ",trim(model_space)
      return
    end if
    add_vel_change=.true.
  end function
  logical function add_change_3d(tochange,add_now)
    real :: tochange(:,:,:)
    logical :: add_now
    integer :: n(3)
    real :: o(3),d(3)
    add_change_3d=.false.
    call set_n_o_d_3d(vel_tomo_struct,n,o,d)
!call srite("mtest.H",velocity_tomo,size(velocity_tomo)*4)
    if (add_now) then
      if (model_space(1:2).eq."ds") then
        velocity_tomo=1./velocity_tomo
        velocity_tomo= velocity_tomo+reshape(tochange,n)
        velocity_tomo=1./velocity_tomo
      else if (model_space(1:2).eq."v2") then
        write(0,*) "ADDING CHANGE",minval(tochange),maxval(tochange)
        velocity_tomo=velocity_tomo**2+reshape(tochange,n)
        where(velocity_tomo<0.)
          velocity_tomo=0. !AVOID NAN
end where
        velocity_tomo=sqrt(velocity_tomo)
      else
        write(0,*) "unknown model space ",trim(model_space)
        return
      end if
    else
      change= change+reshape(tochange,(/product(n)/))
    end if
!call srite("mtest.H",velocity_tomo,size(velocity_tomo)*4)
    add_change_3d=.true.
  end function
  logical function add_change_1d(tochange,add_now)
    real :: tochange(:)
    logical :: add_now
    integer :: n(3)
    real :: o(3),d(3)
    call set_n_o_d_3d(vel_tomo_struct,n,o,d)
    add_change_1d=add_change_3d(reshape(tochange,n),add_now)
  end function
  logical function return_change(mychange)
    real :: mychange(:)
    if (.not. associated(change)) then
      write(0,*) "change not associated"
      return_change=.false.
    else
      mychange=change
      return_change=.true.
    end if
  end function
  logical function transform_vel_orig_space(myv)
    real,pointer :: myv(:,:,:)
    transform_vel_orig_space=.false.
    allocate(myv(size(velocity_tomo,1),size(velocity_tomo,2),size&
      &(velocity_tomo,3)))
    myv=velocity_tomo
    if (tau) then
      if (.not. convert_vel_tau_z(myv)) then
        write(0,*) "trouble converting velocity to depth"
        return
      end if
    end if
!UNPA
    if (.not. unpad_vel(myv)) then
      write(0,*) "trouble unpadding velocity"
      return
    end if
    transform_vel_orig_space=.true.
  end function
  logical function convert_vel_tau_z(myv)
    real, pointer,dimension(:,:,:) :: myv,mapt,temp
    integer :: nz(3),nt(3)
    real :: oz(3),dz(3),ot(3),dt(3)
    convert_vel_tau_z=.false.
    call set_n_o_d_3d(vel_tomo_struct,nt,ot,dt)
    call set_n_o_d_3d(vel_z_tomo_struct,nz,oz,dz)
    allocate(mapt(nt(1),nt(2),nt(3)))
    allocate(temp(nt(1),nt(2),nt(3)))
    temp=myv
    deallocate(myv)
    allocate(myv(nz(1),nz(2),nz(3)))
    if (.not. create_z_map_3d(temp,mapt)) then
      write(0,*) "trouble createing tau map"
      return
    end if
    if (.not. convert_field_tau_z_3d(temp,myv,mapt)) then
      write(0,*) "trouble converting field tau to z"
      return
    end if
    deallocate(mapt,temp)
    convert_vel_tau_z=.true.
  end function
  logical function unpad_vel(array)
    integer :: n(3)
    real :: o(3),d(3)
    real, pointer :: array(:,:,:),temp(:,:,:)
    unpad_vel=.false.
    allocate(temp(size(array,1),size(array,2),size(array,3)))
    temp=array
    deallocate(array)
    call set_n_o_d_3d(vel,n,o,d)
    allocate(array(n(1),n(2),n(3)))
    array=temp(beg(1)+1:beg(1)+n(1),beg(2)+1:beg(2)+n(2),beg(3)+1:beg(3&
      &)+n(3))
    deallocate(temp)
    unpad_vel=.true.
  end function
  logical function update_vel()
    integer :: n(3)
    real :: o(3),d(3)
    update_vel=.false.
    call set_n_o_d_3d(vel_tomo_struct,n,o,d)
    if (model_space(1:2).eq."ds") then
      velocity_tomo=1./(1./velocity_tomo+reshape(change,(/n(1),n(2),n(3&
        &)/)))
    else if (model_space(1:2).eq."v2") then
      velocity_tomo=velocity_tomo**2+reshape(change,(/n(1),n(2),n(3)/)&
        &)
      where(velocity_tomo<0)
        velocity_tomo=0. !NAN bad
end where
      velocity_tomo=sqrt(velocity_tomo)
    else
      write(0,*) "unknown model space method",trim(model_space)
      return
    end if
    update_vel=.true.
  end function
  logical function vel_change(dvel)
    real :: dvel(:,:,:)
    vel_change=.false.
    if (model_space(1:2).eq."ds") then
      if (.not. add_change( (1./(velocity_tomo+dvel) - 1./velocity_tomo&
        &),.true.)) then
        write(0,*) "trouble adding delta slow change"
        return
      end if
    else if (model_space(1:2).eq."v2") then
      if (.not. add_change( (velocity_tomo+dvel)**2-velocity_tomo**2&
        &,.true.)) then
        write(0,*) "trouble adding vel^2 change"
        return
      end if
    else
      write(0,*) "unknown model_space",trim(model_space)
      return
    end if
    vel_change=.true.
  end function
  subroutine init_tomo_coord()
    integer :: n(3)
    real :: o(3),d(3)
    call set_n_o_d_3d(vel_tomo_struct,n,o,d)
    call init_coord(n,o,d)
  end subroutine
  logical function return_pad(b,e)
    real :: b(3),e(3)
    e=end
    b=beg
    return_pad=.true.
  end function
  logical function grab_space_conv_pars()
    grab_space_conv_pars=.false.
    call from_param("beg1_pad",beg(1),0)
    call from_param("beg2_pad",beg(2),0)
    call from_param("end3_pad",end(3),0)
    call from_param("beg3_pad",beg(3),0)
    call from_param("end1_pad",end(1),0)
    call from_param("end2_pad",end(2),0)
    call from_param("tau",tau,.false.)
    call from_param("model_space",model_space,"ds")
    if (model_space(1:2) .ne. "ds" .and. model_space(1:2) .ne. "v2")&
      & then
      write(0,*) "      Only understand model space as ds or v2 saw="&
        &,trim(model_space)
      return
    end if
    if (tau) then
      call from_param("otau",otau,0.)
      call from_param("ntau",ntau,vel%n(1))
      call from_param("dtau",dtau,4.5/ntau)
    end if
    if (sep3d_ndims(vel).eq.2 .and. (end(3) >0 .or. beg(3)>0)) then
      write(0,*) "can't pad axis 3 when we are doing 2-D tomography"
      return
    end if
    grab_space_conv_pars=.true.
  end function
end module
