      PROGRAM fermat
*     Monte Carlo method for finding minimum optical path
      IMPLICIT NONE
      REAL v(101), y(101)
      INTEGER IR_, N
      CALL GWopen(IR_, 0)
      CALL initial(y,v,N)
      CALL change_path(y,v,N)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(y,v,N)
      IMPLICIT NONE
      REAL index, v1, v2, y(101), v(101)
      REAL dmy, rnd
      INTEGER IR_, i, N
      dmy = rnd(-1)
      N = 10                                ! number of different media (even)
*     speed of light in vacuum equal to unity
      v1 = 1.0D0
      index = 1.5D0                         ! index of refraction of medium #2
      v2 = 1/index
      DO i = 1, N/2
         v(i) = v1                          ! speed of light in left half
      END DO
      DO i = N/2 + 1, N
         v(i) = v2                          ! speed of light in right half
      END DO
      y(1) = 1                              ! fixed source
      y(N) = 8                              ! fixed detector
      CALL GWindow(IR_, -0.5, y(1)-0.5, N+1.0, y(N)+0.5)
      CALL GWrect(IR_, 1.0, y(1), REAL(N), y(N))
      CALL GWline(IR_, REAL(N/2), y(1), REAL(N/2), y(N))
                                            ! boundary
*     choose initial path at boundary between endpoints
      DO i = 2, N-1
         y(i) = (y(N) - y(1))*rnd(0) + y(1)
      END DO
      CALL GWsetpen(IR_, 13, -1, -1, -1)
      CALL GWanchor(IR_, 1)
      CALL GWsetogn(IR_, 0, 10)
      CALL GWplot1(IR_, -1, N, 1.0, REAL(N), 0.0,0.0,1.0,0.0, y)
      CALL GWsetogn(IR_, 0, -10)            ! initial path
      END
C
      SUBROUTINE change_path(y,v,N)
      IMPLICIT NONE
      REAL delta, dist, dy2, t_original, t_trial, ytrial, y(101),
     +       v(101)
      REAL rnd
      INTEGER IR_, x, N
      LOGICAL Ltmp1_, TBkeyinput
      delta = 0.5D0                         ! maximum change in y
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
*        choose random x not including x = 1 or N
         x = int((N-2)*rnd(0)) + 2
         ytrial = y(x) + (2*rnd(0) - 1)*delta   ! new y position
         dy2 = (y(x) - y(x+1))**2           ! vertical distance squared
         dist = sqrt(1 + dy2)               ! horizontal distance is unity
         t_original = dist/v(x+1)
         dy2 = (y(x) - y(x-1))**2
         dist = sqrt(1 + dy2)
         t_original = t_original + dist/v(x)
         dy2 = (ytrial - y(x+1))**2
         dist = sqrt(1 + dy2)
         t_trial = dist/v(x+1)
         dy2 = (ytrial - y(x-1))**2
         dist = sqrt(1 + dy2)
         t_trial = t_trial + dist/v(x)
         IF(t_trial .LT. t_original) THEN   ! new position reduces time
            y(x) = ytrial
            CALL GWplot1(IR_, -1, N, 1.0, REAL(N), 0.0,0.0,1.0,0.0, y)
            CALL GWflush(IR_, -10)
         END IF
         Ltmp1_ = .NOT.TBkeyinput()
      END DO
      END
