      subroutine expand(ia2,iatom,ipt,st,icnt,cel,cosi,                 EXPA001
     &                  itot,dmax,did,iortho,t)                         EXPA002
c ----------------------------------------------------------------------EXPA003
c     generates translation vectors for cell expansion                  EXPA004
c     and expands until all distances < dmax are found                  EXPA005
c ----------------------------------------------------------------------EXPA006
c     This program module is part of the PhysiSoft X-ray Library        EXPA007
c     and may not be sold or reproduced for commercial purposes.        EXPA008
c     Although it has been tested carefully, there is no warranty       EXPA009
c     regarding accuracy or freedom from side-effects.                  EXPA010
c                   (C) 1986-1990 Bernhard Rupp                         EXPA011
c ----------------------------------------------------------------------EXPA012
      parameter (istor=2500,iat=20)                                     EXPA013
                                                                        EXPA014
      dimension st(iat,192,4),ipt(iat),p(3),did(istor,8),cel(6),cosi(3) EXPA015
      dimension xo(3),xa(3),t(3,3)                                      EXPA016
c --- expand unit cell around atom --                                   EXPA017
 1001 itot=0                                                            EXPA018
      iend=0                                                            EXPA019
c --- define the central atom --                                        EXPA020
      do 2001 i=1,3                                                     EXPA021
 2001    xo(i)=st(ia2,icnt,i)*cel(i)                                    EXPA022
                                                                        EXPA023
 1002 continue                                                          EXPA024
      do 1010 ia=1,iatom                                                EXPA025
c ---    generate the quadrant translation vectors -                    EXPA026
         do 1011 ik=1,ipt(ia)                                           EXPA027
            is=iend-1                                                   EXPA028
            icyc=0                                                      EXPA029
            do 1012 it1= -iend,iend                                     EXPA030
               do 1013 it2= -iend,iend                                  EXPA031
                  do 1014 it3= -iend,iend                               EXPA032
c ---                skip if position already generated                 EXPA033
                     if (abs(it1).le.is.and.abs(it2).le.is.and.         EXPA034
     &                  abs(it3).le.is) goto 2014                       EXPA035
c ---                add translation vector to each equipoint --        EXPA036
                     p(1)=it1+st(ia,ik,1)                               EXPA037
                     p(2)=it2+st(ia,ik,2)                               EXPA038
                     p(3)=it3+st(ia,ik,3)                               EXPA039
c ---                calculate the distance to newly generated point -  EXPA040
                     do 2003 k=1,3                                      EXPA041
                        xa(k)=p(k)*cel(k)-xo(k)                         EXPA042
 2003                continue                                           EXPA043
                     dd=distan (iortho,xa,cosi)                         EXPA044
c ---                skip if distance larger than d(max) else add --    EXPA045
                     if(dd.gt.dmax) goto 2014                           EXPA046
                     itot=itot+1                                        EXPA047
                     icyc=icyc+1                                        EXPA048
                     do 2029 ij=1,3                                     EXPA049
                        did(itot,ij)=p(ij)                              EXPA050
 2029                continue                                           EXPA051
                     did(itot,4)=dd                                     EXPA052
                     did(itot,5)=ia                                     EXPA053
c ---                calculate cartesian coordinates of atom --         EXPA054
                     call tracor (p,t)                                  EXPA055
                     do 2030 k=1,3                                      EXPA056
                        iix= k+5                                        EXPA057
                        did(itot,iix)=p(k)                              EXPA058
 2030                continue                                           EXPA059
c ---                restart with smaller dmax if storage exceeded --   EXPA060
                     if (itot.ge.istor) then                            EXPA061
                        print*,' Required storage exceeds memory limits'EXPA062
                        print*,' reducing largest bond distance d(max)' EXPA063
                        dmax=dmax/2.0                                   EXPA064
                        print*,' New dmax = ',dmax                      EXPA065
                        goto 1001                                       EXPA066
                     end if                                             EXPA067
 2014                continue                                           EXPA068
 1014             continue                                              EXPA069
 1013          continue                                                 EXPA070
 1012       continue                                                    EXPA071
 1011    continue                                                       EXPA072
 1010 continue                                                          EXPA073
      if(icyc.eq.0.and.iend.gt.3) goto 1003                             EXPA074
      iend=iend+1                                                       EXPA075
      goto 1002                                                         EXPA076
 1003 continue                                                          EXPA077
      return                                                            EXPA078
      end                                                               EXPA079
