# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!
! Inversion functions
!
! Notation:
!
! W_d = Data covariance
! T = Back projection operator
! R = Regularization operator
! P = Preconditioning operator
! K = Model mask operator
! M = Model description operator
! E = Extra operator
! e_r = Epsilon scaling for model styling
! e_e = Epsilon scaling for extra operator
!
module inv_basic_func_mod
  use tomo_chain_mod
  use tomo_space_conv_mod
  use basic_util_mod
  use punish_change_mod
  implicit none
  contains
  logical function tomo_grab_dat0(dat0)
    real, pointer :: dat0(:),sc(:)
    real, pointer :: temp(:,:,:),temp1(:),temp2(:,:,:)
    character(len=128) :: temp_ch
    type(sep3d) :: velm_s,vel_s,vz,vt
    integer :: n(3)
    real :: o(3),d(3)
    integer,external :: auxin
    integer :: ierr,nr,nmod,ndat
    tomo_grab_dat0=.false.
    if (.not. load_tomo_oper(temp1)) then
!d_0
      write(0,*) "        ERROR: loading operator"
      return
    end if
    if (.not. init_data_covar_chain(size(temp1),data_covar_nrange(size&
      &(temp1)))) then
      write(0,*) "trouble initializing data chain size"
      return
    end if
    allocate(dat0(data_covar_nrange(size(temp1))))
    dat0=0
    ndat=size(dat0)
    ierr=data_covar_op(.false.,.false.,temp1,dat0) !
    deallocate(temp1)
!NOW DEAL WITH A POTENTIAL vel_mig
    call from_param("vel_mig",temp_ch,"none")
    if (.not. vel_structure(vel_s)) then
      write(0,*) "trouble grabbing velocity structure"
      return
    end if
    call set_n_o_d_3d(vel_s,n,o,d)
    if (temp_ch(1:4) .ne. "none") then
      allocate(temp(n(1),n(2),n(3)))
      if (-1.eq.auxin("vel_mig")) then
        write(0,*) "specified vel_mig  but not a seplib file"
        return
      end if
      call auxclose("vel_mig")
      call init_sep3d("vel_mig",velm_s,"INPUT")
      if (.not. sep3d_conform(vel_s,velm_s)) then
        write(0,*) "vel_mig=",trim(temp_ch)," not the same  size as &
          & input velocity"
        return
      end if
      if (.not. sreed_field("vel_mig",temp)) then
        write(0,*) "trouble reading in vel_mig"
        return
      end if
      if (.not. convert_vel_tomo_space(temp)) then
        write(0,*) "trouble converting field to tau"
        return
      end if
      if (.not. vel_tomo_3d_pt(myvel=temp2)) then
        write(0,*) "trouble grabbing tomograhy velocity pointer"
        return
      end if
      temp1= temp1+reshape(1./temp2-1./temp,(/size(temp)/))
      deallocate(temp)
    end if
    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,n,o,d)
    allocate(temp1(product(n)))
    if (.not. return_change(temp1)) then
      write(0,*) "trouble grabbing current model space delta"
      return
    end if
    dat0=-dat0
    ierr=data_cov_tomo_op(.false.,.true.,temp1,dat0)
    ! -r_d = W_d d_0 - W_d T (s_0 - s_mig)
    dat0=-dat0
    temp1=0.
    if (.not. add_change(temp1,.true.)) then
      write(0,*) "trouble resetting m0 to zero"
      return
    end if
    deallocate(temp1)
    tomo_grab_dat0=.true.
  end function
end module
