C     SIERP.F
C     PROGRAM TO GENERATE THE SIEPINSKI TRIANGLE

C     ORIGINALLY WRITTEN 1994 WHILE A STUDENT AT UNIV OF WISCONSIN-
C     RIVER FALLS.

C     COPYRIGHT (C) 1994, 2007 JIM HALL, JHALL@FREEDOS.ORG
C Permission is hereby granted, free of charge, to any person obtaining
C a copy of this software and associated documentation files (the
C "Software"), to deal in the Software without restriction, including
C without limitation the rights to use, copy, modify, merge, publish,
C distribute, sublicense, and/or sell copies of the Software, and to
C permit persons to whom the Software is furnished to do so, subject to
C the following conditions:
C
C The above copyright notice and this permission notice shall be
C included in all copies or substantial portions of the Software.
C
C THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
C EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
C MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
C NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
C LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
C OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
C WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

      PROGRAM SIERP
      PARAMETER (ITERMAX = 100000)
      REAL X(3), Y(3)

C     INITIALIZE THE TRIANGLE

      X(1) = 0.
      Y(1) = 0.

      X(2) = 300.
      Y(2) = 600.

      X(3) = 600.
      Y(3) = 0.

      X0 = RAND(0) * 600.
      Y0 = RAND(0) * 600.

C     ITERATE THE TRIANGLE
C     PRINT OUTPUT IN GNUPLOT STYLE

      PRINT *, '## USE GNUPLOT TO DISPLAY THE OUTPUT'

      DO I = 1, ITERMAX, 1
         J = INT(RAND(0) * 3.) + 1
         X0 = (X0 + X(J)) / 2.
         Y0 = (Y0 + Y(J)) / 2.

         PRINT *, X0, Y0
      ENDDO
      END