PROGRAM fermat
! Monte Carlo method for finding minimum optical path
DIM y(101),v(101)
CALL initial(y(),v(),N)
CALL change_path(y(),v(),N)
END

SUB initial(y(),v(),N)
    RANDOMIZE
    LET N = 10                    ! number of different media (even)
    ! speed of light in vacuum equal to unity
    LET v1 = 1.0
    LET index = 1.5               ! index of refraction of medium #2
    LET v2 = 1/index
    FOR i = 1 to N/2
        LET v(i) = v1             ! speed of light in left half
    NEXT i
    FOR i = N/2 + 1 to N
        LET v(i) = v2             ! speed of light in right half
    NEXT i
    LET y(1) = 2                  ! fixed source
    LET y(N) = 8                  ! fixed detector
    SET WINDOW 0.5,N+1,y(1)-0.5,y(N)+0.5
    BOX LINES 1,N,y(1),y(N)
    PLOT LINES: N/2,y(1);N/2,y(N)      ! boundary
    ! choose initial path at boundary between endpoints
    FOR i = 2 to N-1
        LET y(i) = (y(N) - y(1))*rnd + y(1)
    NEXT i
    SET COLOR "red"
    FOR i = 1 to N-1
        PLOT LINES: i,y(i);i+1,y(i+1)  ! initial path
    NEXT i
    SET COLOR "red/white"
END SUB

SUB change_path(y(),v(),N)
    LET delta = 0.5               ! maximum change in y
    DO
       ! choose random x not including x = 1 or N
       LET x = int((N-2)*rnd) + 2
       LET ytrial = y(x) + (2*rnd - 1)*delta     ! new y position
       LET dy2 = (y(x) - y(x+1))^2     ! vertical distance squared
       LET dist = sqr(1 + dy2)    ! horizontal distance is unity
       LET t_original = dist/v(x+1)
       LET dy2 = (y(x) - y(x-1))^2
       LET dist = sqr(1 + dy2)
       LET t_original = t_original + dist/v(x)
       LET dy2 = (ytrial - y(x+1))^2
       LET dist = sqr(1 + dy2)
       LET t_trial = dist/v(x+1)
       LET dy2 = (ytrial - y(x-1))^2
       LET dist = sqr(1 + dy2)
       LET t_trial = t_trial + dist/v(x)
       IF t_trial < t_original then     ! new position reduces time
          SET COLOR "white"
          PLOT LINES: x-1,y(x-1);x,y(x);x+1,y(x+1)
          SET COLOR "red"
          LET y(x) = ytrial
          PLOT LINES: x-1,y(x-1);x,y(x);x+1,y(x+1)
       END IF
    LOOP until key input
END SUB