! seed packing
!
SET COLOR MIX(0) .3,.3,.3

LET s = 400            ! set window size
SET WINDOW 0,3*s,0,3*s
!LET g = 1/(.5*(5^(.5)-1))   ! set spacing factor
!LET g = sqr(2)
LET g = 1/20
!LET g = .5*(exp(1)-1)
!Let g = .5*(pi-3)
LET n = 100             ! set number of seeds
LET a = 20.1          ! set size of (eqil triang) seed
LET i = 1              ! initialize
LET m = a/2
DIM h(200)              ! NOTE this must match n

DO while i < n+1
   LET x =  (i*g - int(i*g)) * s   
   LET gim = x - m     
      IF gim < 0 THEN  ! fix to be mod s
      LET gim = gim + s
      END IF
   LET gip = x + m
      IF gip > s THEN     ! fix to be mod s
      LET gip = gip - s
      END IF
   LET h(i) = m*sqr(3)   !h(i) will be APEX - set to mroot3
   LET j = 1
   DO while j < i
      LET xj = (j*g - int(j*g)) * s
      LET gjm = xj - m    
          IF gjm < 0 THEN
          LET gjm = gjm + s    ! fix to be mod s
          END IF
      LET gjp = xj + m      
          IF gjp > s THEN     ! fix to be mod s
          LET gjp = gjp - s
          END IF
      CALL hij( gim, gip, gjm, gjp)
      IF hnew > h(i) THEN
         LET h(i) = hnew
      END IF
      LET j = j + 1
   LOOP
   CALL tri ( x, h(i))
   LET i = i + 1
   Pause .1
LOOP

SUB hij(gim, gip, gjm, gjp)
   CALL dif (gip, gjm)
   IF d < 0 THEN
        LET hnew = m*sqr(3)
        EXIT SUB
   END IF
   CALL dif (gip, gjm + m)
   IF d < 0 THEN
        LET hnew = h(j) + sqr(3)*(gip-gjm) 
        EXIT SUB
   END IF
   CALL dif (gip, gjp)
   IF d < 0 THEN
        LET hnew = h(j) + m*sqr(3) 
        EXIT SUB
   END IF
   CALL DIF (gim, gjm + m)
   IF d < 0 THEN
        LET hnew = h(j) + m*sqr(3) 
        EXIT SUB
   END IF
   CALL DIF (gim,gjp)
   IF d < 0 THEN
        LET hnew = h(j) - sqr(3)*(gim-gjp) 
        EXIT SUB
   END IF
   LET hnew = m*sqr(3)
END SUB


SUB dif (a, b)
    LET d = a - b
    IF d < -s/2 THEN
       LET d = d + s
       EXIT SUB
    END IF
    IF d > s/2 THEN
       LET d = d - s
       EXIT SUB
    END IF
END SUB


SUB tri (p,q)
    PLOT s/10+p-m+a/10,s/10+q-m*sqr(3)+a/10;s/10+p,s/10+q-2*a/(10*sqr(3))
    PLOT s/10+p,s/10+q-2*a/(10*sqr(3));s/10+p+m-a/10,s/10+q-m*sqr(3)+a/10
    PLOT s/10+p-m+a/10,s/10+q-m*sqr(3)+a/10; s/10+p+m-a/10,s/10+q-m*sqr(3)+a/10
END SUB

GET key zz

END