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