MODULE Shared
  USE, INTRINSIC :: ISO_Fortran_env, dp=>REAL64 !modern DOUBLE PRECISION
  IMPLICIT NONE
  INTEGER :: i, j, k
  INTEGER, PARAMETER :: i_max = 50000, BOUNDARY = 5, N = 10, path_dat = 10 !manage UNIT numbers as a resource
  REAL(dp) :: random, dy, y_new, dist, index1, index1_read, index2, index2_read, t_old, t_new
  REAL(dp), DIMENSION(:), ALLOCATABLE, ASYNCHRONOUS :: y, v
  REAL(dp), PARAMETER :: delta = 0.5, EIGHT = 8.0
  CHARACTER(*), PARAMETER :: my_format = '(i3, f5.1)' !dp not needed
END MODULE Shared
!1. dp?, D+00?
!2. i + 1?
!3. Stupid things with _plot?
PROGRAM fermat
  USE Shared
  IMPLICIT NONE
  CALL INIT_RANDOM_SEED()
  ALLOCATE(y(N + 1), v(N + 1))
  WRITE(6, '(a52)', advance='no') "Index of refraction in the left half of the system: "
  READ *, index1_read
  index1 = index1_read
  left: DO i = 1, BOUNDARY !speed of light in left half of system
    v(i) = 1.0_dp / index1
  END DO left
  WRITE(6, '(a53)', advance='no') "Index of refraction in the right half of the system: "
  READ *, index2_read
  index2 = index2_read
  right: DO i = BOUNDARY + 1, N !speed of light in right half of system
    v(i) = 1.0_dp / index2
  END DO right
  y(1) = 2.0_dp !fixed source
  y(N + 1) = EIGHT !fixed detector
  endpoints: DO i = 2, N !choose initial path between endpoints
    CALL RANDOM_NUMBER(random)
    y(i) = (y(N + 1) - y(1)) * random + y(1)
  END DO endpoints
  i = 1
  light: DO
    IF (i > i_max) EXIT
    CALL RANDOM_NUMBER(random)
    j = INT(FLOAT(N - 1) * random) + 2 !random integer between 2 and 11
    CALL RANDOM_NUMBER(random)
    y_new = y(j) + (2.0_dp * random - 1.0_dp) * delta
    dy = (y(j) - y(j - 1))**2
    dist = SQRT(1.0_dp + dy)
    t_old = dist / v(j - 1)
    dy = (y(j + 1) - y(j))**2
    dist = SQRT(1.0_dp + dy)
    t_old = t_old + dist / v(j)
    dy = (y_new - y(j + 1))**2
    dist = SQRT(1.0_dp + dy)
    t_new = dist / v(j)
    dy = (y_new - y(j - 1))**2
    dist = SQRT(1.0_dp + dy)
    t_new = t_new + dist / v(j - 1)
    IF (t_new < t_old) y(j) = y_new
    i = i + 1
  END DO light
  OPEN(path_dat, file='path.dat', access='sequential', status='unknown', action='write', asynchronous='yes')
  WRITE(path_dat, my_format, asynchronous='yes') (k, y(k), k=1, 11)
  CLOSE(path_dat)
  DEALLOCATE(y, v)
END PROGRAM fermat

SUBROUTINE INIT_RANDOM_SEED()
  USE ISO_Fortran_env, ONLY: INT64
  IMPLICIT NONE
  INTEGER, ALLOCATABLE :: seed(:)
  INTEGER :: i, n, un, istat, dt(8), pid
  INTEGER(INT64) :: t
  CALL RANDOM_SEED(size = n)
  ALLOCATE(seed(n))
  OPEN(newunit=un, file='/dev/urandom', access='stream', status='old', action='read', form='unformatted', iostat=istat)
  IF (istat == 0) THEN
    READ(un) seed
    CLOSE(un)
  ELSE
    CALL SYSTEM_CLOCK(t)
    IF (t == 0) THEN
      CALL DATE_AND_TIME(values = dt)
      t = (dt(1) - 1970) * 365_INT64 * 24 * 60 * 60 * 1000 + dt(2) * 31_INT64 * 24 * 60 * 60 * 1000 + dt(3) * 24_INT64 * 60 * 60   &
        * 1000 + dt(5) * 60 * 60 * 1000 + dt(6) * 60 * 1000 + dt(7) * 1000 + dt(8)
    END IF
    pid = GETPID()
    t = IEOR(t, INT(pid, KIND(t)))
    DO i = 1, n
      seed(i) = lcg(t)
    END DO
  END IF
  CALL RANDOM_SEED(put = seed)
  DEALLOCATE(seed)
CONTAINS
  FUNCTION lcg(s)
    INTEGER :: lcg
    INTEGER(INT64) :: s
    IF (s == 0) THEN
      s = 104729
    ELSE
      s = MOD(s, 4294967296_INT64)
    END IF
    s = MOD(s * 279470273_INT64, 4294967291_INT64)
    lcg = INT(MOD(s, INT(HUGE(0), INT64)), KIND(0))
  END FUNCTION lcg
END SUBROUTINE INIT_RANDOM_SEED