PROGRAM projection
! Find coordinates of a 1-d quasicrystal using
! projection from 2-D lattice. The slope m should be irrational.
! Prints coordinate along line and perpendicular coodinate
! or prints out perpendicular coordinate distribution.
DIM perp(100)
PRINT "1 for list of coordinates or 2 for perpendicular space distribution"
INPUT output
IF output = 2 then
   INPUT prompt "# bins per unit length for perpendicular distribution = ": nbin
END IF
INPUT prompt "Slope of parallel space line between 0 and 1 = ? ": m
INPUT prompt "Intercept of parallel space line = ? ": b
INPUT prompt "Width of 2-d lattice = ? ": n
CALL coordinates(m,b,n,output,perp,nbin)
IF output = 2 then CALL distribution(nbin,perp)
END

SUB coordinates(m,b,n,output,perp(),nbin)
    LET db = 1 + m          ! vertical
    LET mp = -1/m           ! slope of line perpendicular to parallel space line
    FOR ix = 0 to n               ! x coordinate of 2-d lattice point
        LET iy = int(m*ix + b)    ! y coordinate of top lattice point in strip
        IF iy - 1 > m*ix + b-db then     ! Two values of iy for this ix
           CALL nextpoint(ix,iy-1,qold,m,mp,b,nbin,perp,originx,originy,output)
        END IF
        CALL  nextpoint(ix,iy,qold,m,mp,b,nbin,perp,originx,originy,output)
    NEXT ix
END SUB

SUB DISTRIBUTION(nbin,perp())
    ! print perpendicular space distribution
    PRINT
    PRINT "xperp number of points"
    FOR ibin = 1 to 2*nbin
        IF perp(ibin) > 0 then PRINT ibin/nbin,perp(ibin)
    NEXT ibin
END SUB

SUB nextpoint(ix,iy,qold,m,mp,b,nbin,perp(),originx,originy,output)
    ! 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)
    LET p = abs((iy-m*ix - b)/sqr(1+m*m))   ! perpendicular space coordinate
    LET bin = int(nbin*p + 1)
    LET perp(bin) = perp(bin) + 1   ! distribution of perpendicular space points
    IF originx = 0 and originy = 0 then
       LET originx = x
       LET originy = y
       IF output = 1 then
          PRINT
          PRINT "    q      p    q(n) - q(n-1)"
          PRINT
          PRINT using "---.###": 0,p
       END IF
    ELSE
       LET q = sqr((x-originx)^2 + (y-originy)^2)     ! quasicrystal coordinate
       IF output = 1 then PRINT using "---.###" : q,p,q-qold
    END IF
    LET qold = q
END SUB