PROGRAM QuasiCrystal1
! Find coordinates of a 1-d quasicrystal using
! projection from 2-D lattice. The slope m should be irrational.
! Prints distance between neighboring points and plots
! "quasicrystal" lattice.
INPUT prompt "Enter slope m (irrational) between 0 and 1 = ": m
INPUT prompt "Enter intercept = ": b
INPUT prompt "Enter width of 2-d lattice = ": n
CLEAR
ASK PIXELS px,py
LET ar = px/py
LET ymax = int(m*(n+2) + b)
SET WINDOW -2*ar,(n+2)*ar,-2,n+2
FOR ix = 0 to n-1
    FOR iy = 0 to ymax
        BOX LINES ix,ix+1,iy,iy+1      ! draw lattice
    NEXT iy
NEXT ix
PLOT LINES:0,b;n+3,m*(n+3)+b      ! plot parallel space line
LET db = 1 + m                    ! vertical displacement of other side of strip
PLOT LINES :0,b-db;n+3,m*(n+3)+b-db    ! plot other side of strip
LET mp = -1/m
LET np = 0
FOR ix = 0 to n
    LET iy = int(m*ix + b)        ! y coordinate of lattice point
    IF iy-1 > m*ix + b-db then    ! Two values of iy for this ix
       LET np = np + 1
       CALL nextpoint(ix,iy-1,xold,yold,m,mp,b,np)
    END IF
    LET np = np + 1
    CALL nextpoint(ix,iy,xold,yold,m,mp,b,np)
NEXT ix
END

SUB nextpoint(ix,iy,xold,yold,m,mp,b,np)
    ! find location on parallel space coordinate
    LET x = (-mp*ix + iy - b)/(m - mp)    ! interception of parallel space line
    LET y = m*x+b                         ! and perpendicular line from (ix,iy)
    BOX CIRCLE x-.1,x+.1,y-.1,y+.1
    PLOT LINES: ix,iy;x,y
    ASK MAX CURSOR mr,mc
    IF ix > 0 and np < mr then
       SET CURSOR np,70
       LET d = sqr((x - xold)^2 + (y - yold)^2)
       PRINT d                    ! distance between parallel space points
    END IF
    LET xold = x
    LET yold = y
END SUB