!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! THE CHAOTIC MOTION OF DYNAMICAL SYSTEMS: stadium billiard
!  billiard.f90
!  (reference system: origin in the middle point)
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

module bill
  implicit none
  public:: start,move,hit,output,reflect
  real,public::z,w,c,L,m,q,xr1,xr2,a,b,xrmax,xrmin
  !integer,public::i
  integer,public::n,ro
  real,dimension(:),public,allocatable::x,y,px,py
  real,dimension(4),public::xr

contains

  subroutine start()
    ! gives pos., mom., geometry, numbers of steps    
    print*," number of steps >"
    read(unit=5,fmt=*)n
    allocate (x(n))
    allocate (y(n))
    allocate (px(n))
    allocate (py(n))
    print*,"length of the straight side (i nunits of r) >"
    read(unit=5,fmt=*)L
    print*,"initial position (x0,y0) >"
    read(unit=5,fmt=*)x(1),y(1)
    print*,"initial momentum (px0,py0) >"
    read(unit=5,fmt=*)px(1),py(1)
    if((px(1)**2+py(1)**2)/=1.0) then        ! renormalize to 1:
       px(1)=px(1)/sqrt(px(1)**2+py(1)**2)
       py(1)=py(1)/sqrt(px(1)**2+py(1)**2)
    end if
    return 
  end subroutine start

  subroutine hit(i)
    !     given position and momentum of the particles, determine where
    !     it will reflect:
    !     1) top line (y=1,ro=1)
    !     2) bottom line (y=-1,ro=2)
    !     3) right semicircle with center in (L/2,0)(ro=3),
    !     4) left semicircle with center in (-L/2,0)(ro=4)

    integer,intent(in) :: i
    !       real,dimension(n+1)::x,y,px,py

    if(px(i)/=0.0) then
       !cccccccc case px(i)/=0
       m=py(i)/px(i)
       q =y(i)-m*x(i)  ! line through (x,y) and with 
       ! the direction of the momentum
       if(m>0.0) then
          !ccccccccccc subcase m>0
          if(py(i)>0.0) then
             if(((1.0-q)/m)>L/2) then
                ro=3
             else
                if(((1.0-q)/m)>(-L/2)) then
                   ro=1
                else
                   ro=4
                end if
             end if
          else
             if(((-1.0-q)/m)<(-L/2)) then
                ro=4
             else
                if(((-1.0-q)/m)<L/2) then
                   ro=2
                else
                   ro=3 
                end if
             end if
          end if
       end if

       if(m<0.0)then
          !ccccccccccc subcase m<0
          if(py(i)>0.0) then
             if(((1-q)/m)<(-L/2)) then
                ro=4
             else
                if(((1.0-q)/m)<L/2) then              
                   ro=1
                else
                   ro=3
                end if
             end if
          else
             if(((-1.0-q)/m)>L/2) then
                ro=3
             else
                if(((-1.0-q)/m)>(-L/2)) then
                   ro=2
                else
                   ro=4
                end if
             end if
          end if
       end if

       ! subcase m=0 (within the case px(i)=/0)
       if(m==0.0) then
          if(px(i)>0.0) then
             ro=3
          else
             ro=4
          end if
       end if

    else

       ! case px(i)=0
       if(py(i)>0.0) then
          if(abs(x(i))<L/2) then
             ro=1
          else
             if(x(i)>=L/2)ro=3
             if(x(i)<=(-L/2))ro=4
          end if
       else 
          if(abs(x(i))<L/2) then
             ro=2
          else
             if(x(i)>=L/2)ro=3
             if(x(i)<=(-L/2))ro=4
          end if
       end if

    end if

    return
  end subroutine hit

  subroutine reflect(i)
    ! calc. the abscissa of the reflection point, x(i+1), here called xr
    integer,intent(in) :: i

    if(px(i)==0.0) then

       xr(ro)=x(i)

    else

       if((ro==1.0).or.(ro==2.0)) then
          xr(ro)=((-1)**(ro+1) - q)/m
       end if
       if((ro==3.0).or.(ro==4.0)) then
          b=(1.0+m**2)*(q**2+(L**2)/4.0 -1.0)
          a=((m*q+((-1)**ro)*L/2)**2)-b
          xr1=(-(m*q+((-1)**ro)*L/2)+sqrt(a))/(1+m**2)
          xr2=(-(m*q+((-1)**ro)*L/2)-sqrt(a))/(1+m**2)
          xrmin=min(xr1,xr2)
          xrmax=max(xr1,xr2)
       end if
       if(ro==3.0) then
          if((m*py(i))<0.0) then
             xr(ro)=xrmin
          else
             xr(ro)=xrmax
          end if
       end if
       if(ro==4.0) then
          if((m*py(i))<=0.0) then
             xr(ro)=xrmin
          else
             xr(ro)=xrmax
          end if
       end if

    end if

    return
  end subroutine reflect

  subroutine move(i)
    !     determine the position of reflection and momenta after reflection
    integer,intent(in) :: i

    x(i+1)=xr(ro)

    if((ro==1.0).or.(ro==2.0)) then
       px(i+1)=px(i)
       py(i+1)= - py(i)
       y(i+1)= (-1)**(ro+1.0)
    end if

    if((ro==3.0).or.(ro==4.0)) then
       if(px(i)/=0.0) then
          y(i+1)= m*x(i+1)+q
       else
          y(i+1)=sign(1.,py(i))*sqrt(1-(x(i+1)+((-1)**ro)*L/2)**2) 
       end if
       a=(y(i+1)**2-(x(i+1)+((-1)**ro)*L/2)**2)*px(i)
       b=2*(x(i+1)+((-1)**ro)*L/2)*y(i+1)*py(i)
       px(i+1)=a-b          
       a=-2*(x(i+1)+((-1)**ro)*L/2)*y(i+1)*px(i)
       b=((x(i+1)+((-1)**ro)*L/2)**2-y(i+1)**2)*py(i)
       py(i+1)=a+b
    end if

    return 
  end subroutine move

  subroutine output(i)
    !     write positions and momenta on a file
    integer,intent(in) :: i
    write(unit=3,fmt=*)i,x(i),y(i),px(i),py(i)
    return
  end subroutine output

end module bill


program billiard
  use bill
  integer::j
  call start()
  open(unit=3,file="dati",status="replace",action="write")

  do j=1,n-1
     call hit(j)
     call reflect(j)
     call move(j)
     call output(j)
  end do

  deallocate (x)
  deallocate (y)
  deallocate (px)
  deallocate (py)

  close(3)
end program billiard