module slo_
  use util_
  use outnull_
  use dottest_

  implicit none

  integer, pointer, private :: cc(:,:)
  real,    pointer, private :: ww(:,:),xx(:,:)
  integer,          private :: n
  type(myaxis),     private :: a1,a2,aa,al
  real,             private :: pi
  logical,          private :: verb

contains
  !----------------------------------------------------------------
  subroutine slo_init(a1in,a2in,aain,alin,verbin)
    type(myaxis), intent(in) :: a1in,a2in,aain,alin
    logical,      intent(in) :: verbin

    a1 = a1in
    a2 = a2in
    aa = aain
    al = alin

    verb=verbin

    allocate(cc(al%n,2),ww(al%n,2))
    pi   = 4*atan(1.)

    allocate(xx(2,2))
  end subroutine slo_init
  !----------------------------------------------------------------
  subroutine slo_reset(ia)
    integer,     intent(in) :: ia
    integer                 :: il
    real                    :: a,tana,sina,cosa
    real                    :: l,d1,d2

    a=aa%o+(ia-1)*aa%d
    a = a/180*pi
    tana = tan(a)
    sina = sin(a)
    cosa = cos(a)

    do il=1,al%n
       l=al%o+(il-1)*al%d

       d1=l*sina
       d2=l*cosa

       cc(il,1)=nint(d1/a1%d)
       cc(il,2)=nint(d2/a2%d)

       ww(il,1) = 1.-( d1/a1%d - cc(il,1))
       ww(il,2) = 1.-( d2/a2%d - cc(il,2))
    end do

  end subroutine slo_reset
  !----------------------------------------------------------------
  integer function slo_lop(adj,add,mm,dd) result(stat)
    logical, intent(in)    :: adj,add
    real,    intent(inout) :: mm(:),dd(:)

    call outnull(adj,add,mm,dd)
    stat = slo_lop2(adj,add,mm,dd)
  end function slo_lop
  !----------------------------------------------------------------
  integer function slo_lop2(adj,add,mm,dd)
    logical, intent(in)             :: adj,add
    real, dimension(a1%n,a2%n)      :: dd
    real, dimension(a1%n,a2%n,aa%n) :: mm
    integer                         :: ia,i1,i2

    if(adj) then
       do ia=1,aa%n    !! for all dip angles
          if(verb) write(0,*) "a",ia,aa%n
          call slo_reset(ia)
          write(0,*) "slo reset OK"
          do i1=1,a1%n !! for every location
             do i2=1,a2%n
                call slo_i2s(i1,i2,ia,mm,dd)
             end do
          end do
       end do
    else
       do ia=1,aa%n    !! for all dip angles
          if(verb) write(0,*) "a",ia,aa%n
          call slo_reset(ia)

          do i1=1,a1%n !! for every location
             do i2=1,a2%n
                call slo_s2i(i1,i2,ia,mm,dd)
             end do
          end do
       end do
    end if

  end function slo_lop2
  !----------------------------------------------------------------
  subroutine slo_i2s(i1,i2,ia,m,d)
    integer, intent(in)    :: i1,i2,ia
    real,    intent(inout) :: d(:,:), m(:,:,:)
    integer                :: j1,j2,il
    integer                :: k1,k2

!write(0,*) size(xx),size(ww)

    xx(1,1) =    ww(il,1)  *    ww(il,2)
    xx(2,1) = (1-ww(il,1)) *    ww(il,2)
    xx(1,2) =    ww(il,1)  * (1-ww(il,2))
    xx(2,2) = (1-ww(il,1)) * (1-ww(il,2))

!write(0,*) xx
    do il=1,al%n
       j1=i1+cc(il,1)
       j2=i2+cc(il,2)

       if(j1<1 .or. j1>a1%n-1) cycle
       if(j2<1 .or. j2>a2%n-1) cycle

       do k2=1,2
          do k1=1,2
             m(i1,i2,ia) = &
             m(i1,i2,ia) + &
             d(j1+k1-1,j2+k2-1) * xx(k1,k2)
          end do
       end do

    end do

  end subroutine slo_i2s
  !----------------------------------------------------------------
  subroutine slo_s2i(i1,i2,ia,m,d)
    integer, intent(in)    :: i1,i2,ia
    real,    intent(inout) :: d(:,:), m(:,:,:)
    integer                :: j1,j2,il
    integer                :: k1,k2

    xx(1,1) =    ww(il,1)  *    ww(il,2)
    xx(2,1) = (1-ww(il,1)) *    ww(il,2)
    xx(1,2) =    ww(il,1)  * (1-ww(il,2))
    xx(2,2) = (1-ww(il,1)) * (1-ww(il,2))

    do il=1,al%n
       j1=i1+cc(il,1)
       j2=i2+cc(il,2)

       if(j1<1 .or. j1>a1%n-1) cycle
       if(j2<1 .or. j2>a2%n-1) cycle

       do k2=1,2
          do k1=1,2
             d(j1+k1-1,j2+k2-1) = &
             d(j1+k1-1,j2+k2-1) + &
             m(i1,i2,ia) * xx(k1,k2)
          end do
       end do

    end do

  end subroutine slo_s2i
  !----------------------------------------------------------------
  subroutine slo_close()
    deallocate(cc,ww)
  end subroutine slo_close
end module slo_
