program makecc
      implicit none
      character*80 fn,s80
      character*106 s106
      character*1 s1,s2
      integer*4 nat,NA1,NA2,it,ig,NN,ia,ix,iz,izf,nal,nt,ii
      real*8 bohr,x(3),r(3),u(3,3)
      logical lunit
      integer*4,allocatable::ll(:)
c
      bohr=0.529177d0
      lunit=.false.
      nat=0
      ig=0
      it=0
      write(6,6000)
6000  format(' Making big FILE.X movie file from CASTEP output',/)
      if(iargc().gt.0)then
       call getarg(1,fn)
       call getarg(2,s80)
       read(s80,*)NA1
       call getarg(3,s80)
       read(s80,*)NA2
       call getarg(4,s80)
       read(s80,*)NN
      else
       write(6,*)'Trajectory filename:'
       read(5,200)fn
200    format(a80)
       write(6,6001)
6001   format(
     1 'Interval of atoms, geom int (N1 N2 NN, "0 0 0" for all): ',$)
       read(5,*)NA1,NA2,NN
      endif

      if(NA1.lt.0)then
       write(6,*)'SOLU.LST requested by negative N1'
       open(9,file='SOLU.LST',status='old')
       read(9,*)nal
       allocate(ll(nal))
       read(9,*)ll
       close(9)
      endif


c
c     read separately the first geometry to find number of atoms:
      s1='X'
      open(4,file=fn)
11    read(4,900,err=99,end=99)s106
900   format(a106)
      s2=s106(106:106)
      if(s2.eq.'R')nat=nat+1
      if(s2.ne.'R'.and.s1.eq.'R')goto 99
      s1=s2
      if(s106(26:55).eq.'Total number of ions in cell =')then
       read(s106(56:60),*)nat
       goto 99
      endif
      goto 11
99    write(6,*)nat,' atoms'
      close(4)
      if(nat.eq.0)stop

      if(NN.eq.0)NN=1
      if(NA1.eq.0)NA1=1
      if(NA2.eq.0)NA2=nat
      if(NA1.lt.0)then
       nt=nal
      else
       nt=NA2-NA1+1
      endif

      open(4,file=fn)
      open(3,file='FILE.X')
      s1='X'
1     read(4,900,err=999,end=999)s106
c      
c     New format:      
      if(s106(39:47).eq.'Unit Cell')then
       lunit=.true.
       read(4,*)
       read(4,*)
       read(4,*)u(1,1),u(1,2),u(1,3)
       read(4,*)u(2,1),u(2,2),u(2,3)
       read(4,*)u(3,1),u(3,2),u(3,3)
      endif
      if(s106(39:69).eq.'Fractional coordinates of atoms')then
       if(.not.lunit)then
        write(6,*)'Fractions found, but unit cell unknown'
        stop
       endif       
       read(4,*)
       read(4,*)
       ig=ig+1
       if(mod(ig,NN).eq.0.or.ig.eq.1)then
        it=it+1
        write(6,*)ig,it
        write(3,203)ig,it
        write(3,*)nt
        do 3 ia=1,nat
        read(4,900)s106
        if(s106(16:16).eq.' ')then
         iz=izf(s106(17:18))
        else
         iz=izf(s106(16:17))
        endif
        read(s106(34:68),*)x
        do 4 ix=1,3
4       r(ix)=u(1,ix)*x(1)+u(2,ix)*x(2)+u(3,ix)*x(3)
        if(NA1.lt.0)then
         do 5 ii=1,nal
5        if(ll(ii).eq.ia)write(3,300)iz,(r(ix),ix=1,3)
300      format(i3,3f12.6)
        else
         if(ia.ge.NA1.and.ia.le.NA2)write(3,300)iz,(r(ix),ix=1,3)
        endif
3       continue
       endif
      endif

c
c     Old format:      
      s2=s106(106:106)
      if(s2.eq.'R'.and.s1.ne.'R')then
       ig=ig+1
       if((ig/NN)*NN.eq.ig)then
        it=it+1
        write(3,203)ig,it
203     format('CASTEP GEOMETRY ',2I12)
        write(3,*)nt
        backspace 4
        do 2 ia=1,nat
        read(4,900)s106
        read(s106(19:99),*)x
        iz=izf(s106(2:3))
        if(NA1.lt.0)then
         do 6 ii=1,nal
6        if(ll(ii).eq.ia)write(3,300)iz,(x(ix)*bohr,ix=1,3)
        else
         if(ia.ge.NA1.and.ia.le.NA2)write(3,300)iz,(x(ix)*bohr,ix=1,3)
        endif
2       continue
       else
        do 21 ia=1,nat-1
21      read(4,*)
       endif
      endif
      s1=s2

      goto 1
999   close(4)
      close(3)
      write(6,*)' FILE.X made.'
      write(6,600)ig,it
600   format(I8,' geometries found;',I8,' written.')
      stop
      end

      function izf(sy)
      implicit none
      character*2 sy
      integer*4 izf,i
      character*2 as(87)
      data as/
     1 'H ','He'
     1,'Li','Be','B ','C ','N ','O ','F ','Ne'
     1,'Na','Mg','Al','Si','P ','S ','Cl','Ar'
     1,'K ','Ca','Sc','Ti','V ','Cr','Mn','Fe','Co','Ni','Cu','Zn',
     1                                              'Ga','Ge','As'
     1,'Se','Br','Kr'
     1,'Rb','Sr','Y ','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd',
     1                                              'In','Sn','Sb'
     1,'Te','I ','Xe'
     1,'Cs','Ba'
     1,'La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er',
     1                                              'Tm','Yb','Lu'
     1,'Hf','Ta','W ','Re','Os','Ir','Pt','Au','Hg','Tl','Pb','Bi',
     1                                              'Po','At','Rn','L '/

      do 1 i=1,87
      if(sy.eq.as(i))then
       izf=i
       return
      endif
  1   continue
      write(6,*)'Unknown atom '//sy
      izf=0
      return
      end