      subroutine heapst (n,rb,is)                                       HEAP001
c ----------------------------------------------------------------------HEAP002
c     generic heap sort algorithm for sorting distances (after          HEAP003
c     NUMERICAL RECIPES                                                 HEAP004
c     XLIB Library (C) Physisoft Corp. 1988                             HEAP005
c ----------------------------------------------------------------------HEAP006
      parameter (m=8,ist=2500)                                          HEAP007
      dimension rb(ist,m),rrb(m)                                        HEAP008
                                                                        HEAP009
      l=n/2+1                                                           HEAP010
      ir=n                                                              HEAP011
10    continue                                                          HEAP012
        if(l.gt.1)then                                                  HEAP013
          l=l-1                                                         HEAP014
          rra=rb(l,is)                                                  HEAP015
          do 1001 ia=1,m                                                HEAP016
 1001        rrb(ia)=rb(l,ia)                                           HEAP017
        else                                                            HEAP018
          rra=rb(ir,is)                                                 HEAP019
          do 1002 ia=1,m                                                HEAP020
 1002        rrb(ia)=rb(ir,ia)                                          HEAP021
          rb(ir,is)=rb(1,is)                                            HEAP022
          do 1003 ia=1,m                                                HEAP023
 1003        rb(ir,ia)=rb(1,ia)                                         HEAP024
          ir=ir-1                                                       HEAP025
          if(ir.eq.1)then                                               HEAP026
            rb(1,is)=rra                                                HEAP027
            do 1004 ia=1,m                                              HEAP028
 1004          rb(1,ia)=rrb(ia)                                         HEAP029
            return                                                      HEAP030
          endif                                                         HEAP031
        endif                                                           HEAP032
        i=l                                                             HEAP033
        j=l+l                                                           HEAP034
20      if(j.le.ir)then                                                 HEAP035
          if(j.lt.ir)then                                               HEAP036
            if(rb(j,is).lt.rb(j+1,is))j=j+1                             HEAP037
          endif                                                         HEAP038
          if(rra.lt.rb(j,is))then                                       HEAP039
            rb(i,is)=rb(j,is)                                           HEAP040
            do 1005 ia=1,m                                              HEAP041
 1005          rb(i,ia)=rb(j,ia)                                        HEAP042
            i=j                                                         HEAP043
            j=j+j                                                       HEAP044
          else                                                          HEAP045
            j=ir+1                                                      HEAP046
          endif                                                         HEAP047
        go to 20                                                        HEAP048
        endif                                                           HEAP049
        rb(i,is)=rra                                                    HEAP050
        do 1006 ia=1,m                                                  HEAP051
 1006      rb(i,ia)=rrb(ia)                                             HEAP052
      go to 10                                                          HEAP053
      end                                                               HEAP054
