# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module shoot_mod
  use double_mod
  implicit none
  integer, save,private ::verb
  integer, save,private :: nvel(2),nt,nangle,zero_offset,j,ntout
  real,save,private :: ovel(2),dvel(2),ot,dt,dtout
  real,save ,private :: delmin,delta
  real(kind=wp),save ,private:: maxvel(2)
  logical,save,private :: take_off_angles
  real,save,private :: oxsrc,dxsrc,oang,dang,ozsrc,dzsrc
  type axis
    integer n
    real o
    real d
  end type
  contains
  subroutine get_shot_params(vel_axis,verb_in,raypar_rite,raymod_rite&
    &,tau,j1)
    type(axis) :: vel_axis(2)
    integer :: nloc,verb_in,output_tau,j1
    logical :: raypar_rite,raymod_rite,tau
    integer,external :: auxout,getch
    character(len=128) :: temp_ch
    integer :: n1, n2, verb, output_tau, zero_offset, j1
    real :: o1,o2,d1,d2

    from history: integer n1:vel_axis(1)%n,n2:vel_axis(2)%n
    from history: real o1:vel_axis(1)%o,o2:vel_axis(2)%o
    from history: real d1:vel_axis(1)%d,d2:vel_axis(2)%d
    from param: integer verb:verb_in=0,output_tau=0,zero_offset=0,j1=1
    if (0.eq.getch("raypar","s",temp_ch)) then
      raypar_rite=.false.
    else
      raypar_rite=.true.
    end if
    if (0.eq.getch("raymod","s",temp_ch)) then
      raymod_rite=.false.
    else
      raymod_rite=.true.
    end if
    if (0.eq.output_tau) then
      tau=.false.
    else
      tau=.true.
    end if
    oxsrc=0
    dxsrc=1
    ozsrc=0
    dzsrc=0
    oang=0
    dang=1
  end subroutine
  subroutine init_shot(nv,ov,dv,n_t,d_t,nang)
    integer :: nv(2),n_t,nang
    real :: dv(2),d_t,ov(2)
    nvel=nv
    ovel=ov
    dvel=dv
    maxvel=ovel+dvel*(nvel-1)
    nt=n_t
    dt=d_t
    nangle=nang
  end subroutine
  subroutine init_shot_params(vel_axis,nloc,verb_in,sources,angles,j1)
    type(axis) :: vel_axis(2)
    integer :: verb_in,auxin,tempi,nloc,sreed,i1,getch,j1
    complex,pointer :: sources(:),angles(:)
    real :: tempr
    character(len=128) :: temp_ch
    integer :: nt
    real :: ot,dt

    nvel=(/vel_axis(1)%n,vel_axis(2)%n/)
    ovel=(/vel_axis(1)%o,vel_axis(2)%o/)
    dvel=(/vel_axis(1)%d,vel_axis(2)%d/)
    maxvel=ovel+dvel*(nvel-1)
    call from_param( ot,dt from param: integer nt
    verb=verb_in
    j=j1
    ntout=ceiling(1.0* nt/j)
    dtout=j*dt
!FIRST TRY FILES
    if (-1.ne.auxin("sources")) then
!WE HAVE SOURCES PROVIDED
  from aux: sources integer n1:nloc,esize:tempi
      if (tempi.ne.8) then
        call erexit("esize not equal to 8")
      end if
      allocate(sources(nloc),angles(nloc))
      if (nloc*8.ne. sreed("sources",sources,size(sources)*8)) then
        call erexit("trouble reading sources")
      end if
    else
      call from_param("nsrc",nloc)
      oxsrc=0.;dxsrc=0.;ozsrc=0.;dzsrc=0.
      allocate(sources(nloc),angles(nloc))
      do i1=1,nloc
        sources(i1)=cmplx(oxsrc+dxsrc*(i1-1),ozsrc+dzsrc*(i1-1))
      end do
    end if
    if (zero_offset.eq.0) then
      call from_param("nang",nangle)
      if (1.eq.getch("angles","s",temp_ch)) then
        if (-1.ne.auxin("angles")) then
          call from_aux ("angles","n1",tempi)
          if (tempi.ne.nloc) then
            call erexit("number of angles not equal number of sources"&
              &) call from_aux("angles","esize",tempi)
          end if
          if (tempi.ne.8) then
            call erexit("angles mus be complex")
          end if
          if (nloc*8.ne. sreed("angles",angles,size(angles)*8)) then
            call erexit("trouble reading ngles")
          end if
        end if
      else
        call from_param("oang",oang)
        call from_param("dang",dang)
        angles=cmplx(oang,oang+dang*(nangle-1))
      end if
    else
      nangle=1
      do i1=1,nloc-1
        tempr=atan((aimag(sources(i1+1))-aimag(sources(i1)))/(real(sources&
          &(i1+1))-real(sources(i1))))*180/pi
        if (zero_offset.eq.-1) then
          tempr=-tempr+180
        end if
        angles(i1)=cmplx(tempr,tempr*2)
      end do
      angles(nloc)=angles(nloc-1)
    end if
  end subroutine
  subroutine rite_shot_params(nrays,nloc,rite_raypar,rite_raymod)
    integer nloc,nrays
    logical :: rite_raypar,rite_raymod
    real :: xsou,dx_sou
    real :: osrc,dsrc
    if (dzsrc.ne.0. .and. dxsrc.ne.0.) then
      osrc=0
      dsrc=1.
    else if (dzsrc.ne.0.) then
      osrc=ozsrc
      dsrc=dzsrc
    else
      osrc=oxsrc
      dsrc=dxsrc
    end if
    to history: integer n1:ntout,esize:8
!CHANGE WHEN MATCHING
 to history: real o1:ot,d1:dtout,delmin:delmin,delta:delta,o2:oang,d2:dang,o3:osrc,d3:dsrc to history: integer n2:nangle,n3:nloc
    if (rite_raypar) then
      to aux: raypar integer n1:ntout,esize:8,n2:nangle,n3:nloc to aux:&
        & raypar real o1:ot,d1:dtout,o2:oang,d2:dang,o3:osrc,d3:dsrc
    end if
    if (rite_raymod) then
      to aux: raymod integer n1:ntout,esize:8,n2:nangle,n3:nloc to aux:&
        & raymod real o1:ot,d1:dtout,o2:oang,d2:dang,o3:osrc,d3:dsrc
    end if
  end subroutine
  subroutine read_vel(velocity)
    integer, external :: sreed
    real :: velocity(:,:)
    if (size(velocity)*4.ne.sreed("in",velocity,size(velocity)*4))&
      & then
      call erexit("trouble reading velocity")
    end if
  end subroutine
  subroutine ray_trace(vel,x_ray,z_ray,v_ray,i1,i2)
    real :: vel(:,:)
    real(kind=wp) :: p2_x,p2_z,p2_v,q_x,q_z,q_v,p1_x,p1_z,p1_v
    real(kind=wp) :: p3_x,p3_z,p3_v,a,b,c,r2,alpha,beta,ds1,ds2
    real(kind=wp) :: ll,x_ray(:,:),z_ray(:,:),v_ray(:,:)
    real(kind=wp) :: s2_x,s2_z,s1_x,s1_z,sina,cosa,tana
    integer :: i1,ss,sa,i2
    p2_x=x_ray(i1,i2)
    p2_z=z_ray(i1,i2)
    p2_v=v_ray(i1,i2)
    q_x=x_ray(i1-1,i2)
    q_z=z_ray(i1-1,i2)
    q_v=v_ray(i1-1,i2)
    LL=sqrt((P2_x-Q_x)*(P2_x-Q_x)+(P2_z-Q_z)*(P2_z-Q_z))
    if (P2_z.ne.Q_z) then
      tana=(P2_x-Q_x)/(P2_z-Q_z)
      if (tana>=0) then
        ss=1
      else
        ss=-1
      end if
      tana=abs(tana)
      sina=tana/sqrt(1+tana*tana)
      cosa=1/sqrt(1+tana*tana)
    else
      cosa=0
      sina=1
      ss=1
    end if
    P1_x=P2_x-LL*cosa
    P1_z=P2_z+ss*LL*sina
    P1_v=get_vel(vel,p1_x,p1_z)
    P3_x=P2_x+LL*cosa
    P3_z=P2_z-ss*LL*sina
    P3_v=get_vel(vel,P3_x,p3_z)
    a= P1_x-P3_x
    if (a>=0) then
      sa=1
    else
      sa=-1
    end if
    b= P1_z-P3_z
    c=(P1_v-P3_v)*dt
!
    R2=P2_v*dt
    ALPHA=c/(a*a+b*b)
    BETA=sa*sqrt(a*a+b*b-c*c)/(a*a+b*b)
    S1_x=P2_x-R2*(ALPHA*a-BETA*b)
    S1_z=P2_z-R2*(ALPHA*b+BETA*a)
    S2_x=P2_x-R2*(ALPHA*a+BETA*b)
    S2_z=P2_z-R2*(ALPHA*b-BETA*a)
    DS1=distance(S1_x,s1_z,Q_x,q_z)
    DS2=distance(S2_x,s2_z,Q_x,q_z)
    if (DS1>DS2) then
      x_ray(i1+1,i2)=S1_x
      z_ray(i1+1,i2)=s1_z
    else
      x_ray(i1+1,i2)=S2_x
      z_ray(i1+1,i2)=s2_z
    end if
    v_ray(i1+1,i2)=get_vel(vel,x_ray(i1+1,i2),z_ray(i1+1,i2))
  end subroutine
  subroutine shoot_rays_new(x_ray,z_ray,v_ray,vel,source,angle)
    real(kind=wp),dimension(:,:) :: x_ray,z_ray,v_ray
    real,dimension(:,:) :: vel
    type(axis) :: angle_axis
    integer :: i1,i2,ia
    real(kind=wp) :: p1_x,p2_x,p3_x,q_x
    real(kind=wp) :: p1_z,p2_z,p3_z,q_z
    real(kind=wp) ::dd
    complex :: source,angle
    logical :: zero
    call init_rays_new(x_ray,z_ray,v_ray,vel,source,angle)
    i1=2
    zero=.false.
    do while (i1<=nt-1 .and..not. zero)
      do i2=2,nangle-1
        p1_x=x_ray(i1,i2-1)
        p1_z=z_ray(i1,i2-1)
        p2_x=x_ray(i1,i2)
        p2_z=z_ray(i1,i2)
        p3_x=x_ray(i1,i2+1)
        p3_z=z_ray(i1,i2+1)
        q_x=x_ray(i1-1,i2)
        q_z=z_ray(i1-1,i2)
        dd=distance(P3_x,p3_z,P1_x,p1_z)
        if (iscusp(q_x,q_z,p1_x,p1_z,p2_x,p2_z,p3_x,p3_z)) then
          call ray_trace(vel,x_ray,z_ray,v_ray,i1,i2)
        else
          call wfttr(vel,x_ray,z_ray,v_ray,i1,i2)
        end if
      end do
!remove this check for rays where z can be negative
      if (all(z_ray(i1,:)<0.)) then
        zero=.true.
        do ia=i1+1,nt
          x_ray(ia,:)=x_ray(i1+1,:)
          z_ray(ia,:)=z_ray(i1+1,:)
          v_ray(ia,:)=v_ray(i1+1,:)
        end do
      end if
      i1= i1+1
    end do
  end subroutine
  logical function iscusp(q_x,q_z,p1_x,p1_z,p2_x,p2_z,p3_x,p3_z)
    real(kind=wp) :: q_x,q_z,p1_x,p1_z,p2_x,p2_z,p3_x,p3_z
    real :: j1,j3
    j1=(P2_x-P1_x)*(P2_z-Q_z)-(P2_x-Q_x)*(P2_z-P1_z)
    j3=(P2_x-P3_x)*(P2_z-Q_z)-(P2_x-Q_x)*(P2_z-P3_z)
    if (j1*j3 >=0) then
      iscusp=.true.
    else
      iscusp=.false.
    end if
  end function
  subroutine wfttr(vel,x_ray,z_ray,v_ray,i1,i2)
    real(kind=wp) :: x_ray(:,:),z_ray(:,:),v_ray(:,:)
    real(kind=wp) :: p2_x,p2_z,p2_v,q_x,q_z,q_v,p1_x,p1_z,p1_v
    real(kind=wp) :: p3_x,p3_z,p3_v,a,b,c,r2,alpha,beta,ds1,ds2
    real(kind=wp) :: s2_x,s2_z,s1_x,s1_z
    real :: vel(:,:)
    integer i1,i2,sa
    p1_x=x_ray(i1,i2-1)
    p1_z=z_ray(i1,i2-1)
    p1_v=v_ray(i1,i2-1)
    p2_x=x_ray(i1,i2)
    p2_z=z_ray(i1,i2)
    p2_v=v_ray(i1,i2)
    p3_x=x_ray(i1,i2+1)
    p3_z=z_ray(i1,i2+1)
    p3_v=v_ray(i1,i2+1)
    q_x=x_ray(i1-1,i2)
    q_z=z_ray(i1-1,i2)
    q_v=v_ray(i1-1,i2)
    a= P1_x-P3_x
    b= P1_z-P3_z
    c=(P1_v-P3_v)*dT
    if (a>=0) then
      sa=1
    else
      sa=-1
    end if
    r2=p2_v*dt
    ALPHA=c/(a*a+b*b)
    BETA=sa*sqrt(a*a+b*b-c*c)/(a*a+b*b)
    S1_x=P2_x-R2*(ALPHA*a-BETA*b)
    S1_z=P2_z-R2*(ALPHA*b+BETA*a)
    S2_x=P2_x-R2*(ALPHA*a+BETA*b)
    S2_z=P2_z-R2*(ALPHA*b-BETA*a)
    DS1=distance(S1_x,s1_z,Q_x,q_z)
    DS2=distance(S2_x,s2_z,Q_x,q_z)
    if (DS1>DS2) then
      x_ray(i1+1,i2)=s1_x
      z_ray(i1+1,i2)=s1_z
    else
      x_ray(i1+1,i2)=s2_x
      z_ray(i1+1,i2)=s2_z
    end if
! write(0,*) "BEFORE GET VEL 4"
    v_ray(i1+1,i2)=get_vel(vel,x_ray(i1+1,i2),z_ray(i1+1,i2))
  end subroutine
  subroutine init_rays_new(x_ray,z_ray,v_ray,vel,source,angle)
    real,dimension(:,:) ::vel
    real(kind=wp),dimension(:,:) :: x_ray,z_ray,v_ray
    real(kind=wp) :: ang,sx,sz,vv,d
    complex :: source,angle
    real :: oangle,dangle
    type(axis) :: angle_axis
    integer :: i1,n
!AAAA
    sx=real(source)
    sz=aimag(source)
    oangle=real(angle)
    if (nangle>1) then
      dangle=(aimag(angle)-oangle)/(nangle-1)
    else
      dangle=1
    end if
    x_ray=0
    v_ray=0
    z_ray=0
! write(0,*) "BEFORE GET VEL 5"
    vv=get_vel(vel,sx,sz)
    d=dt*vv
    n=nangle
    do i1=1,n
      x_ray(1,i1)=sx
      z_ray(1,i1)=sz
      v_ray(1,i1)=vv
      ang=(oangle+dangle*(i1-1))*pi/180
!write(0,*) "CHECK B1",i1,x_ray(1,i1),z_ray(1,i1),d,size(x_ray,1),size(x_ray,2),ang,sin(ang)
      x_ray(2,i1)=x_ray(1,i1)+d*sin(ang)
      z_ray(2,i1)=z_ray(1,i1)+d*cos(ang)
!write(0,*) "CHECK ALPHA",i1,x_ray(2,i1),z_ray(2,i1)
      v_ray(2,i1)=get_vel(vel,x_ray(2,i1),z_ray(2,i1))
    end do
    do i1=2,nt-1
      call ray_trace(vel,x_ray,z_ray,v_ray,i1,1)
    end do
    if (nangle>1) then
      do i1=2,nt-1
        call ray_trace(vel,x_ray,z_ray,v_ray,i1,n)
      end do
    end if
  end subroutine
  real(kind=wp) function get_vel(vel,sx,sz)
    real :: vel(:,:)
    real(kind=wp) :: sx,sz
    real(kind=wp) :: x,z,fractx,fractz
    integer :: ix,iz
    x=(sx-ovel(2))/dvel(2)+1.
    z=(sz-ovel(1))/dvel(1)+1.
!write(0,*) "CHECK",sx,sz,ovel,dvel
    if (x<1) then
      ix=1
      fractx=0.
    else if (x>nvel(2)) then
      ix=nvel(2)-1
      fractx=1.
    else
      ix=floor(x)
      fractx=x-ix
    end if
    if (z<1) then
      iz=1
      fractz=0.
    else if (z>nvel(1)) then
      iz=nvel(1)-1
      fractz=1.
    else
      iz=floor(z)
      fractz=z-iz
    end if
    if (ix<1 .or. ix >= size(vel,2) .or. iz <1 .or. iz>= size(vel,1))&
      & then
      write(0,*) "x and z",x,z,fractx,fractz
      write(0,*) "this is very bad",ix,iz,size(vel,1),size(vel,2)
    end if
! write(0,*) "BEFORE GET VEL 6"
    get_vel=(1.-fractx)*(1.-fractz)*vel(iz,ix)+(1.-fractx)*(fractz)*vel&
      &(iz+1,ix)+(fractx)*(1.-fractz)*vel(iz,ix+1)+(fractx)*(fractz)&
      &*vel(iz+1,ix+1)
  end function
  logical function find_corner(q_x,q_z,p1_x,p1_z,p2_x,p2_z,p3_x,p3_z)
    real(kind=wp) :: q_x,q_z
    real(kind=wp) :: p1_x,p1_z
    real(kind=wp) :: p2_x,p2_z
    real(kind=wp) :: p3_x,p3_z
    integer s1,s2
    s1=sign_vector_product(q_x,q_z,p1_x,p1_z,p2_x,p2_z)
    s2=sign_vector_product(q_x,q_z,p2_x,p2_z,p3_x,p3_z)
    if (s1.eq.s2) then
      find_corner=.false.
    else
      find_corner=.true.
    end if
  end function
  integer function sign_vector_product(a_x,a_z,b_x,b_z,c_x,c_z)
    real(kind=wp) :: a_x,a_z,b_x,b_z,c_x,c_z,product
    product=(b_x-a_x)*(c_z-a_z)-(c_x-a_x)*(b_z-a_z)
    if (product>=0.) then
      sign_vector_product=1
    else
      sign_vector_product=-1
    end if
  end function
  subroutine calc_raypar(x_ray,z_ray,v_ray,n,raypar,time_ref)
    real(kind=wp) :: x_ray(:,:),z_ray(:,:),v_ray(:,:)
    real, optional :: time_ref(:)
    complex :: raypar(:,:)
    integer :: i1,i2,n,bn
    real :: angle
    logical :: exit
    if (present(time_ref)) then
      exit=.true.
    else
      exit=.false.
    end if
    do i2=1,n
      angle=calc_angle(x_ray(1:2,i2),z_ray(1:2,i2))
      raypar(1,i2)=cmplx(sin(angle)/v_ray(1,i2),cos(angle)/v_ray(1,i2)&
        &)
      do i1=2,nt-1
        angle=(calc_angle(x_ray(i1-1:i1,i2),z_ray(i1-1:i1,i2))+calc_angle&
          &(x_ray(i1:i1+1,i2),z_ray(i1:i1+1,i2)))/2.
        raypar(i1,i2)=cmplx(sin(angle)/v_ray(i1,i2),cos(angle)/v_ray(i1&
          &,i2))
      end do
      angle=calc_angle(x_ray(nt-1:nt,i2),z_ray(nt-1:nt,i2))
      raypar(nt,i2)=cmplx(sin(angle)/v_ray(nt,i2),cos(angle)/v_ray(nt&
        &,i2))
    end do
    if (exit) then
      do i2=1,n
        bn=(time_ref(i2)-ot)/dt+2.
        angle=calc_angle(x_ray(bn-1:bn,i2),z_ray(bn-1:bn,i2))
        raypar(bn,i2)=cmplx(sin(angle)/v_ray(bn,i2),cos(angle)/v_ray(bn&
          &,i2))
        angle=calc_angle(x_ray(bn:bn+1,i2),z_ray(bn:bn+1,i2))
        raypar(bn+1,i2)=cmplx(sin(angle)/v_ray(bn+1,i2),cos(angle)/v_ray&
          &(bn+1,i2))
      end do
    end if
  end subroutine
  real(kind=wp) function calc_angle(x,z)
    real(kind=wp) :: x(2),z(2),dx,dz
    dx=x(2)-x(1)
    dz=z(2)-z(1)
    calc_angle=atan2(dx,dz)
  end function
  real(kind=wp) function distance(s2_x,s2_z,q_x,q_z)
    real(kind=wp) :: s2_x,s2_z,q_x,q_z
    distance=(s2_x-q_x)**2+(s2_z-q_z)**2
  end function
end module
