# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module tomo_oper_mod
  use process_types
  use model_types
  use ray_types
  implicit none
  implicit none
  integer, private :: tau_domain,approx,new
  integer, private :: ndim_mod,nv_mod,nx_mod,ny_mod
  real, private :: dv_mod,dx_mod,dy_mod,v0_mod,x0_mod,y0_mod
  integer, private :: ndim_ray,nt_ray,n_ray,nref_ray
  real, private :: t0_ray
  real, private :: dt_ray
  real, pointer, dimension (:,:,:,:), private :: traj_ray
  real, pointer, dimension (:,:,:,:), private :: rpar_ray
  real, pointer, dimension (:,:,:,:), private :: mod_ray
  real, pointer, dimension (:,:), private :: tref_ray
  real, pointer, dimension (:,:), private :: ttot_ray
  real, pointer, dimension (:,:,:), private :: sigma0_op,slow0_op
!integer,pointer,dimension(:,:), private :: zero_offset
  contains
  subroutine tomo_oper_init(param,model_param,ray,sigma0,slow0)
    type (tomo_param_type) param
    type (model_param_type) model_param
    type (ray_obj_type) ray
    real, pointer, dimension (:,:,:) :: sigma0,slow0
!integer, pointer, dimension (:,:) :: zero_o
    tau_domain=param%tau_domain
    approx=param%approx
    new=param%new
    ndim_mod=model_param%ndim_mod
    nv_mod=model_param%nv_mod
    nx_mod=model_param%nx_mod
    ny_mod=model_param%ny_mod
    v0_mod=model_param%v0_mod
!write(0,*) "VOOOOOOOOO",v0_mod
    x0_mod=model_param%x0_mod
    y0_mod=model_param%y0_mod
    dv_mod=model_param%dv_mod
    dx_mod=model_param%dx_mod
    dy_mod=model_param%dy_mod
    ndim_ray=ray%ndim_ray
    nt_ray=ray%nt_ray
    n_ray=ray%n_ray
    nref_ray=ray%nref_ray
    t0_ray=ray%t0_ray
    dt_ray=ray%dt_ray
    if (.not.associated (traj_ray, ray%traj_ray)) then
      traj_ray => ray%traj_ray
    end if
    if (.not.associated (rpar_ray, ray%rpar_ray)) then
      rpar_ray => ray%rpar_ray
    end if
    if (.not.associated (mod_ray, ray%mod_ray)) then
      mod_ray => ray%mod_ray
    end if
    if (.not.associated (ttot_ray, ray%ttot_ray)) then
      ttot_ray => ray%ttot_ray
    end if
    if (.not.associated (tref_ray, ray%tref_ray)) then
      tref_ray => ray%tref_ray
    end if
!if (.not.associated (zero_offset, zero_o)) zero_offset => zero_o
    if (.not.associated (slow0_op, slow0)) then
      slow0_op => slow0
! write(0,*) "good association",size(slow0),size(slow0_op)
    end if
    if (tau_domain .eq. 1) then
      if (.not.associated (sigma0_op, sigma0)) then
        sigma0_op => sigma0
      end if
    end if
  end subroutine
  function tau_tomo (adj, add, slow, time_ray) result (stat)
    logical, intent(in) :: adj,add
    integer :: stat
    real, dimension (:) :: slow
    real, dimension (:) :: time_ray
    stat=tau_tomo_op (adj, add, slow, time_ray)
  end function
  function tau_tomo_op (adj, add, slow, time_ray) result (stat)
    logical, intent(in) :: adj,add
    integer :: stat
    real, dimension (nv_mod,nx_mod,ny_mod) :: slow
    real, dimension (n_ray,nref_ray) :: time_ray
    real, dimension (nv_mod,nx_mod,ny_mod) :: dz,dtau
    real, dimension (nv_mod,nx_mod,ny_mod) :: dsdtau,dsigdtau
    double precision*16, dimension (nv_mod,nx_mod,ny_mod) :: dsigma&
      &,ddz
    integer i_ref,i_ray,it,i_dim,nt_ref,nt_inc,itemp
    integer iv,ix,iy,iz_m1,itau_m1,idz_ref,iv_up
    real dl,dt_ray_ref,vel_ray,sigma_ray,idt_ray,dt_ray_min,rtemp
    real pz_ref,pz_inc,delta_s,dz_ref,dz_ref_max,dz_ref_min
    real d_tau,z,zz,wz,z0,vp,vm,vel_m1
    real d_z0,tau,ttau,wtau,tau0,sp,sm,slow_m1
    real coeff1
    double precision*16 coeff2,coeff3,coeff4,delta_t
    double precision ddiff
    integer icoord_mod(3),icoord_ray(3)
    real id_mod(3),delta_mod(3)
    real id_ray(3),delta_ray(3)
    write(0,*) "alpha"
! zero vectors if necessary
    if (.not. add) then
      if (adj) then
        slow=0.
      else
        time_ray=0.
      end if
    end if
    write(0,*) "zero"
    idt_ray=1./dt_ray
!write(0,*) "2ero"
    dt_ray_min=.01*dt_ray
!write(0,*) "3ero"
    id_mod(1)=1./dv_mod
!write(0,*) "4ero"
    id_mod(2)=1./dx_mod
!write(0,*) "5ero"
    id_mod(3)=1./dy_mod
!write(0,*) "6ero"
    id_ray(1)=1./dx_mod
    id_ray(2)=1./dv_mod
    id_ray(3)=1./dy_mod
    if (ndim_ray .eq. 2) then
      icoord_mod(3)=1
    end if
!i always use a 2-d ray field
    dsigma=0.0d0
    ddz=0.0d0
    if (.not. adj) then
! write(0,*) "in not adj",size(sigma0_op),size(slow0_op)
      do iy=1,ny_mod
        do ix=1,nx_mod
          do iv=1,nv_mod
! write(0,*) "looping",iv,ix,iy
            coeff4=( 2.*sigma0_op(iv,ix,iy)/slow0_op(iv,ix,iy))
            dsigma(iv,ix,iy)= coeff4 * slow(iv,ix,iy)
          end do
        end do
      end do
      do iy=1,ny_mod
        do ix=1,nx_mod
          do iv=2,nv_mod
! write(0,*) "itwo",iv,ix,iy
            ddz(iv,ix,iy)= ddz(iv-1,ix,iy) -.5*dv_mod*slow(iv-1,ix,iy)&
              &/ (slow0_op(iv-1,ix,iy)**2)
          end do
        end do
        do iv=1,nv_mod
          do ix=2,nx_mod-1
! write(0,*) "three",ix,iy
            coeff3=dble(-slow0_op(iv,ix,iy) * id_mod(2))
            dsigma(iv,ix,iy)= dsigma(iv,ix,iy)+ coeff3 * &
              &(ddz(iv,ix+1,iy)-ddz(iv,ix-1,iy))
          end do
          coeff3=dble(-2.*slow0_op(iv,1,iy) * id_mod(2))
          dsigma(iv,1,iy)= dsigma(iv,1,iy)+ coeff3 * (ddz(iv,2,iy)-ddz&
            &(iv,1,iy))
          coeff3=dble(-2.*slow0_op(iv,nx_mod,iy) * id_mod(2))
          dsigma(iv,nx_mod,iy)= dsigma(iv,nx_mod,iy)+coeff3 * (ddz(iv&
            &,nx_mod,iy)-ddz(iv,nx_mod-1,iy))
        end do
      end do
    end if
!#####
    write(0,*) "one"
    do iy=1,ny_mod
      do ix=1,nx_mod
        do iv=2,nv_mod-1
          dsdtau(iv,ix,iy)=.5*id_mod(1)*(slow0_op (iv+1,ix,iy)-slow0_op&
            & (iv-1,ix,iy))
          dsigdtau(iv,ix,iy)=.5*id_mod(1)* &
            &(sigma0_op(iv+1,ix,iy)-sigma0_op(iv-1,ix,iy))
        end do
        dsdtau(1,ix,iy)= dsdtau(2,ix,iy)
        dsdtau(nv_mod,ix,iy)= dsdtau(nv_mod-1,ix,iy)
        dsigdtau(1,ix,iy)= dsigdtau(2,ix,iy)
        dsigdtau(nv_mod,ix,iy)= dsigdtau(nv_mod-1,ix,iy)
      end do
    end do
    write(0,*) "two",size(dtau)
    dtau=0.
    if (.not. adj) then
      do iy=1,ny_mod
        do ix=1,nx_mod
          do iv=2,nv_mod
            dtau(iv,ix,iy)= dtau(iv-1,ix,iy) -.5*dv_mod*slow(iv-1,ix,iy&
              &)/(slow0_op(iv-1,ix,iy)**2)
          end do
          do iv=1,nv_mod
            dtau(iv,ix,iy)= - dtau(iv,ix,iy)*(2.*slow0_op(iv,ix,iy))
          end do
        end do
      end do
    end if
    write(0,*) "three",adj,add
!#####
    do i_ref=1,nref_ray
      do i_ray=1,n_ray
! write(0,*) i_ray,i_ref
        nt_inc=((tref_ray(i_ray,i_ref)-t0_ray)*idt_ray)+1
        nt_ref=((ttot_ray(i_ray,i_ref)-tref_ray(i_ray,i_ref))*idt_ray)&
          &+1
        do it=2,nt_inc + nt_ref +2
          do i_dim=1,ndim_ray
            delta_ray(i_dim)=(traj_ray(i_dim,it ,i_ray,i_ref)- &
              & traj_ray(i_dim,it-1,i_ray,i_ref))
            icoord_ray(i_dim) = (traj_ray(i_dim,it-1,i_ray,i_ref)+.5&
              &*delta_ray(i_dim)) *id_ray(i_dim) +1
! if(.not. adj && add) write(0,*) "coord",traj_ray(i_dim,it-1,i_ray,i_ref), delta_ray(i_dim) ,id_ray(i_dim)
          end do
! needed because rays are in (x,v) instead of (v,x)
          icoord_mod(1)=icoord_ray(2)
          icoord_mod(2)=icoord_ray(1)
          delta_mod(1)=delta_ray(2)
          delta_mod(2)=delta_ray(1)
          vel_ray=mod_ray(1,it-1,i_ray,i_ref)+.5*(mod_ray(1,it ,i_ray&
            &,i_ref) -mod_ray(1,it-1,i_ray,i_ref))
          if ((it /= nt_inc+1) .and. (it /= nt_inc+2)) then
            dt_ray_ref=dt_ray
          else
            if (it .eq. nt_inc+1) then
              dt_ray_ref=tref_ray(i_ray,i_ref)-(nt_inc-1)*dt_ray +&
                & t0_ray
            else
              dt_ray_ref=ttot_ray(i_ray,i_ref)-tref_ray(i_ray,i_ref) -&
                & (nt_ref-1)*dt_ray + t0_ray
            end if
          end if
          dl=dt_ray_ref*vel_ray
          if (dt_ray_ref < dt_ray_min) then
            dl=1.
            dt_ray_ref=1.
            delta_mod(1)=0.
            delta_mod(2)=0.
          end if
          coeff1=(delta_mod(2)*delta_mod(2)/dl)
          delta_t=0.0d0
!if(.not. adj && add) write(0,*) "b",i_ray,i_ref,size(time_ray)
          if (adj) then
            slow(icoord_mod(1),icoord_mod(2),icoord_mod(3)) = slow&
              &(icoord_mod(1),icoord_mod(2),icoord_mod(3))+coeff1 *&
              & time_ray(i_ray,i_ref)
          else
!if(.not. adj && add) write(0,*) "pos",icoord_mod
            time_ray(i_ray,i_ref) = time_ray(i_ray,i_ref) + &
              & delta_t +coeff1 * slow(icoord_mod(1),icoord_mod(2),icoord_mod(3))
          end if
        end do
      end do
    end do
    write(0,*) "four"
!#########
    if (adj) then
      ddz=0.0d0
      do iy=1,ny_mod
        do iv=1,nv_mod
          do ix=2,nx_mod-1
            if (dsigma(iv,ix,iy) .ne. 0.) then
              coeff3=dble(-slow0_op(iv,ix,iy) * id_mod(2))*dsigma(iv,ix&
                &,iy)
              ddz(iv,ix+1,iy) = ddz(iv,ix+1,iy) + coeff3
              ddz(iv,ix-1,iy) = ddz(iv,ix-1,iy) - coeff3
            end if
          end do
          if (dsigma(iv,1,iy) /= 0.) then
            coeff3=dble(-2.*slow0_op(iv,1,iy) * id_mod(2))*dsigma(iv,1&
              &,iy)
            dz(iv,2,iy) = dz(iv,2,iy) + coeff3
            dz(iv,1,iy) = dz(iv,1,iy) - coeff3
          end if
          if (dsigma(iv,nx_mod,iy) /= 0.) then
            coeff3=dble(-2.*slow0_op(iv,nx_mod,iy) * id_mod(2))*dsigma&
              &(iv,nx_mod,iy)
            dz(iv,nx_mod,iy) = dz(iv,nx_mod,iy) + coeff3
            dz(iv,nx_mod-1,iy) = dz(iv,nx_mod-1,iy) - coeff3
          end if
        end do
        do ix=1,nx_mod
          do iv=2,nv_mod
            if (ddz(iv,ix,iy) /= 0.) then
              do iv_up=iv-1,1,-1
                slow(iv_up,ix,iy) = slow(iv_up,ix,iy) - &
                  & ddz(iv,ix,iy)*.5*dv_mod/ (slow0_op(iv_up,ix,iy)**2)
              end do
            end if
          end do
        end do
        do iv=1,nv_mod
          do ix=1,nx_mod
            coeff4=( 2.*sigma0_op(iv,ix,iy)/slow0_op(iv,ix,iy))
            slow(iv,ix,iy) = slow(iv,ix,iy) + coeff4* dsigma(iv,ix,iy)
          end do
        end do
      end do
    end if
!write(0,*) "five"
    if (adj) then
      do iy=1,ny_mod
        do ix=1,nx_mod
          do iv=2,nv_mod
            if (dtau(iv,ix,iy) /= 0.) then
              dtau(iv,ix,iy)= - dtau(iv,ix,iy)*(2.*slow0_op(iv,ix,iy))
              do iv_up=iv-1,1,-1
                slow(iv_up,ix,iy) = slow(iv_up,ix,iy) - &
                  & dtau(iv,ix,iy) * .5*dv_mod/(slow0_op(iv_up,ix,iy)**2)
              end do
            end if
          end do
        end do
      end do
    end if
!write(0,*) "six"
    stat=0
  end function
  function tomo_oper (adj, add, slow, time_ray) result (stat)
    logical, intent(in) :: adj,add
    integer :: stat
    real, dimension (:) :: slow
    real, dimension (:) :: time_ray
    stat=tomo_oper_apply (adj, add, slow, time_ray)
  end function
  function tomo_oper_apply (adj, add, slow, time_ray) result (stat)
    logical, intent(in) :: adj,add
    integer :: stat
    real, dimension (nv_mod,nx_mod,ny_mod) :: slow,temp_a,temp_b
    real, dimension (n_ray,nref_ray) :: time_ray
    real, dimension (nv_mod,nx_mod,ny_mod) :: dz,dtau
    real, dimension (nv_mod,nx_mod,ny_mod) :: dsdtau,dsigdtau
    integer i_ref,i_ray,it,i_dim,nt_ref,nt_inc,itemp
    integer iv,ix,iy,iz_m1,itau_m1,idz_ref,iv_up
    real dl,dt_ray_ref,vel_ray,sigma_ray,idt_ray,dt_ray_min,rtemp
    real pz_ref,pz_inc,delta_s,dz_ref,dz_ref_max,dz_ref_min
    real d_tau,z,zz,wz,z0,vp,vm,vel_m1
    real d_z0,tau,ttau,wtau,tau0,sp,sm,slow_m1
    real coeff1
    integer icoord_mod(3),icoord_ray(3),one
    real id_mod(3),delta_mod(3)
    real id_ray(3),delta_ray(3),o_r(3),o_m(3)
! zero vectors if necessary
    if (.not. add) then
      if (adj) then
        slow=0.
      else
        time_ray=0.
      end if
    end if
    one=(0.-v0_mod)/dv_mod
    temp_a=slow
    slow(1:nv_mod-2*one,:,:)=temp_a(one+1:nv_mod-one,:,:)
!write(0,*) "in b"
    idt_ray=1./dt_ray
    dt_ray_min=.01*dt_ray
    id_mod(1)=1./dv_mod
    id_mod(2)=1./dx_mod
    id_mod(3)=1./dy_mod
!write(0,*) "in c"
    id_ray(1)=1./dx_mod
    id_ray(2)=1./dv_mod
    id_ray(3)=1./dy_mod
    if (ndim_ray .eq. 2) then
      icoord_mod(3)=1
    end if
    do i_ref=1,nref_ray
!loop over the reflectors
      do i_ray=1,n_ray
!loop over the rays
! write(0,*) tref_ray(i_ray,i_ref),ttot_ray(i_ray,i_ref),t0_ray,idt_ray
        nt_inc=((tref_ray(i_ray,i_ref)-t0_ray)*idt_ray)+1
        !index of reflection
        nt_ref=((ttot_ray(i_ray,i_ref)-tref_ray(i_ray,i_ref))*idt_ray)&
          &+1 !index of
! write(0,*) i_ray,i_ref,"loop",nt_ref,nt_inc
!total ray
        do it=2,nt_inc + nt_ref +2
!loop over the portion of the ray we actually care about
          do i_dim=1,ndim_ray
!loop over x,z, and possibly y
            delta_ray(i_dim)=(traj_ray(i_dim,it ,i_ray,i_ref)- &
              & traj_ray(i_dim,it-1,i_ray,i_ref)) !calculate ray change
!in position
            icoord_ray(i_dim) = (traj_ray(i_dim,it-1,i_ray,i_ref)+.5&
              &*delta_ray(i_dim)) *id_ray(i_dim) +1 !index location inside
!slowness model for this ray
!element
          end do
! needed because rays are in (x,v) instead of (v,x)
          icoord_mod(1)=icoord_ray(2)
          icoord_mod(2)=icoord_ray(1)
          delta_mod(1)=delta_ray(2)
          delta_mod(2)=delta_ray(1)
          vel_ray=mod_ray(1,it-1,i_ray,i_ref)+.5*(mod_ray(1,it ,i_ray&
            &,i_ref) -mod_ray(1,it-1,i_ray,i_ref))
!what is mod_ray again?
          if ((it /= nt_inc+1) .and. (it /= nt_inc+2)) then
!if we aren't at the reflection
            dt_ray_ref=dt_ray
          else
!this calculates the time of the ray segment near the reflection
            if (it .eq. nt_inc+1) then
              dt_ray_ref=tref_ray(i_ray,i_ref)-(nt_inc-1)*dt_ray +&
                & t0_ray
            else
              dt_ray_ref=ttot_ray(i_ray,i_ref)-tref_ray(i_ray,i_ref) -&
                & (nt_ref-1)*dt_ray + t0_ray
            end if
          end if
          dl=dt_ray_ref*vel_ray !now the length of the ray segment
! icoord_mod=(/2,2,1/)
! write(0,*) "icoord",icoord_mod
          if (adj) then
            slow(icoord_mod(1),icoord_mod(2),icoord_mod(3))= slow&
              &(icoord_mod(1),icoord_mod(2),icoord_mod(3))+ time_ray(i_ray,i_ref)*dl
          else
!the change in the time equals the change in slowness time the length
!of the ray going through the grid cell
            time_ray(i_ray,i_ref) = time_ray(i_ray,i_ref) + dl&
              &*slow(icoord_mod(1),icoord_mod(2),icoord_mod(3))
          end if
        end do
!calculate the location of the reflection
        do i_dim=1,ndim_ray
          icoord_ray(i_dim)=(traj_ray(i_dim,nt_inc+1,i_ray,i_ref))*id_ray&
            &(i_dim)+1
        end do
        icoord_mod(1)=icoord_ray(2)
        icoord_mod(2)=icoord_ray(1)
        pz_inc= rpar_ray(2,nt_inc+1,i_ray,i_ref)
        pz_ref= rpar_ray(2,nt_inc+2,i_ray,i_ref)
        ! ATTENTION approximation
        if (adj) then
          dz(icoord_mod(1),icoord_mod(2),icoord_mod(3)) = dz(icoord_mod&
            &(1),icoord_mod(2),icoord_mod(3)) + time_ray(i_ray&
            &,i_ref) * (pz_inc - pz_ref)
          d_tau=.5*dz(iv,icoord_mod(2),icoord_mod(3))/slow0_op(iv&
            &,icoord_mod(2),icoord_mod(3))
          do iv_up=iv-1,1,-1
            slow(iv_up,icoord_mod(2),icoord_mod(2))= slow(iv_up&
              &,icoord_mod(2),icoord_mod(3)) + 2.*d_tau*dv_mod
          end do
        else
          do iv=2,nv_mod
            dz(iv,icoord_mod(2),icoord_mod(3))=dz(iv-1,icoord_mod(2)&
              &,icoord_mod(3))+2.* slow(iv-1,icoord_mod(2),icoord_mod(3))*dv_mod
          end do
          dz_ref= -dz(icoord_mod(1),icoord_mod(3),icoord_mod(3))*(.5&
            &/slow0_op(iv,icoord_mod(3), icoord_mod(3)))
          idz_ref=dz_ref*id_mod(1)+1
          time_ray(i_ray,i_ref) = time_ray(i_ray,i_ref) + dz_ref &
            &* (pz_inc - pz_ref)
        end if
      end do
    end do
    temp_b=slow
    slow=temp_a
!write(0,*) "one",one
    slow(one+1:nv_mod-one,:,:)=slow(one+1:nv_mod-one,:,:)+temp_b(1:nv_mod&
      &-2*one,:,:)
    stat=0
  end function
end module
