  % Title "Intersections Demo"
 '% Animated
 '% ScaleUp
 '% PlaceCentral
 '% AnchorCentral
  uses consoleG
  uses GLO2\GeoPlanar


  sub drawintersections()
  =======================


  type linex
    float x1,y1 'first point
    float x2,y2 'second point
    float le    'length
    int dl      'deletion flag
    int p1,p2   'point indexes
  end type

   % q=31

  string tab="    "
  point  p
  int    a
  float le

  seed=0x12379823

  static MoveableObject pn[100]
  static linex dn[10000]
  float r[16]
  indexbase 0


 int i,j,k
  '
  '
  'POLYGONS

  point qn[q]

  if opening
    qn[0]={0,0}
    CreatePolypoints qn[1],0,0,15,1,1,0
    CreatePolypoints qn[16],0.03,0.04,15,0.7,0.75,12.5

    MoveableObject *mp
    for i=0 to q-1
      @mp=@pn[i]
      mp.id=500+i*100
      mp.mode=0x111
      mp.x=10*qn[i].x
      mp.y=10*qn[i].y
      mp.s.x=1.0
      mp.s.y=1.0
      mp.s.z=1.0
      mp.sa.x=1.0
      mp.sa.y=1.0
      mp.sa.z=1.0
    next
  endif 'opening

  k=0
  float x,y
  int qh=(q-1)/2
  int mi,mj
  j=0
  'CLEAR HUB LINE REFS
  for i=0 to q-1
    pn[i].n=0
  next
  '
  'RADIALS
  for i=1 to qh
    mi=i+qh
    mj=0
    le=sqdist pn[mi], pn[mj]
    dn[k]={ pn[mi].x,pn[mi].y,pn[mj].x,pn[mj].y,le,0,mi,mj }
    k++
  next
  '
  'CIRCUMS
  '
  for i=1 to qh
    j=i+qh
    le=sqdist pn[i], pn[j]
    dn[k]={ pn[i].x,pn[i].y,pn[j].x,pn[j].y,le,0,i,j }
    k++
    j=i+qh-1 : if i=1 then j=qh+qh 'overlap
    le=sqdist pn[i], pn[j]
    dn[k]={ pn[i].x,pn[i].y,pn[j].x,pn[j].y,le,0,i,j }
    k++
    '
    j=i+1 : if i=qh then j=1 'overlap
    le=sqdist pn[i], pn[j]
    dn[k]={ pn[i].x,pn[i].y,pn[j].x,pn[j].y,le,0,i,j }
    k++
    mi=i+qh
    mj=j+qh
    le=sqdist pn[mi], pn[mj]
    dn[k]={ pn[mi].x,pn[mi].y,pn[mj].x,pn[mj].y,le,0,mi,mj }
    k++
  next
  '

  'shading
  flat 'default

  '
  color 0,1,1
  printl "lines " tab k
  '
  move 10,-10
  '
  color 0,1,1
  int e
  e=k-1
  linex dln[1000]
  int dc
  int p1,p2,n
  MoveableObject*ph
  for i=0 to e
    mi=i
    if not dn[mi].dl
      dln[dc]=dn[mi]
      p1=dln[dc].p1 'first point index
      p2=dln[dc].p2 'second point index
      @ph=@pn[p1] 'first hub
      n=ph.n
      ph.pp[n]=p2 'store other hub index
      ph.n++
      @ph=@pn[p2] 'second hub
      n=ph.n
      ph.pp[n]=p1 'store other hub index
      ph.n++
      color 0.7,0,0
      drawline @dln[dc] 'DELAUNAY LINES
      dc++ 'index next dln
    endif
  next

  '
  'SORT EACH HUB LINES IN ASCENDING ANGLE ORDER
  MoveableObject*pt
  float ang[32],an,xe,ye
  int m
  string s,t
  'HUBS
  for i=1 to q-1
    @pt=@pn[i]
    n=pt.n
    x=pt.x
    y=pt.y
    'COLLECT LINE ANGLES
    t=""
    for j=0 to n-1
      m=pt.pp[j]
      xe=pn[m].x
      ye=pn[m].y
      an=deg(atan(ye-y, xe-x))
      if an<0 then an+=360
      ang[j]=an
    next
    'SORT INTO ASCENDING ORDER
    for j=0 to n-2
      for k=j+1 to n-1
        if ang[j]>ang[k]
          swap ang[k],ang[j]
          swap pt.pp[k],pt.pp[j] 
        endif
      next
    next
    't="hub  " i "  " n " lines" cr
    for j=0 to n-1
      m=pt.pp[j]
      xe=pn[m].x
      ye=pn[m].y
      line l1={x,y,xe,ye}
      'drawline @l1
      't+=" " j "  " int(ang(j)) cr
    next
    s+=t
    'GET MIDPOINTS OF EACH ADJACENT LINE
    int m1,m2
    point qc1,qc2,qi
    point q1,q2,q3,q4,q5
    q1={x,y}
    pt.pp[n]=pt.pp[0] 'overlap
    ang[n]=ang[0]+360 'overlap
    for j=0 to n-1
      k=j+1
      m1=pt.pp[j]
      m2=pt.pp[k]
      q2={pn[m1].x, pn[m1].y}
      q3={pn[m2].x, pn[m2].y}
      'line l1={x,y,q2.x,q2.y}
      'drawline @l1
      midpoint q1, q2, qc1
      midpoint q1, q3, qc2
      pointsize 3
      drawpoint @qc1
      drawpoint @qc2
      '
      'EXCLUDE
      if ang[k]-ang[j]>=180
        continue for 'j
      end if
      'GET INTERSECTION POINT OF EACH LINE PAIR  PERPENDICULARS
      rot90  qc1,q1,q4
      rotm90 qc2,q1,q5
      line l1,l2,l3,l4
      l1={qc1.x,qc1.y,q4.x,q4.y}
      l2={qc2.x,qc2.y,q5.x,q5.y}
      a=intersect l1,l2,qi
      if a'=16
        if inrange -10,10,qi.x
          if inrange -10,10,qi.y
            float ac=deg(atan(qi.y-pt.y,qi.x-pt.x))
            if ac<0 then ac+=360
            if ac<ang[j] and ang[k]>=360 then ac+=360
            if not inrange ang[j],ang[k],ac
'mbox " " int(ang(j)) " " int(ang(k)) " " int(ac)
              l3={qc1.x,qc1.y,qc2.x,qc2.y}
              color 1,0,1
              thickness 2
              'drawline @l3 
              continue for
            endif
            'pointsize 4
            'color 1,0,0
            'drawpoint @qi
            'VORONOI LINES
            l3={qc1.x,qc1.y,qi.x,qi.y}
            color 0,1,1
            thickness 2
            drawline @l3 
            l4={qc2.x,qc2.y,qi.x,qi.y}
            drawline @l4
          endif
        endif  
      endif
    next 'j
  next 'i
'putfile "t.txt",s
 
          
 
  shading
  color 1,1,0
  for i=0 to q-1
    pushstate
      pn[i].act
      scale .5
      go sphere
    popstate
  next
  picklabel 0
  end sub


  sub main()
  ==========
  string s
  cls 0,0,0
  pushstate
  scale 2,2
  'typeface=1
  printl "Intersection Points"
  popstate
  '
  move 0,-2
  '
  PushState
  move 1,-1
  DrawIntersections
  PopState
  '
  
  end sub

  EndScript
