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

LET s = 400         ! set window size
SET WINDOW 0,s,0,s
!LET g = .5*(5^(.5)-1)   ! set spacing factor
!LET g = sqr(2)
LET g = exp(1)
!LET g = .5*(exp(1)-1)
LET n = 100             ! set number of seeds
LET a = 30             ! set size of (eqil triang) seed
LET i = 1              ! initialize
LET m = a/2
DIM h(150)              ! 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 (a,b)
    PLOT a-m,b-m*sqr(3);a,b
    PLOT a,b;a+m,b-m*sqr(3)
    PLOT a-m,b-m*sqr(3); a+m,b-m*sqr(3)
END SUB

GET key zz

END