*   program hklmcal
*
C  PROGRAM TO COMPUTE 2-THETA'S AND HKLm'S FROM LATTICE PARAMETERS
*  and modulation vector q*
C
C  PROGRAMMED BY: R. L. SNYDER
C		  N.Y.S. COLLEGE OF CERAMICS
C		  ALFRED UNIVERSITY
C		  ALFRED, N.Y. 14802
*  and modified to handle modulated structures by:
*                 David P. Matheis
*                 May 30, 1990
C
      INTEGER H, HMAX, HMIN, inten(200)
      real jq1, jq2, jq3, DO(200), qINT(200)
      CHARACTER*75 NFILE, MFILE, title
      dimension d(2000), thet(2000), ind(2000,4)
      DATA WAVE / 1.5405981/, IREAD/5/, ITYP/6/
      HMIN = 0
      HMAX = 0
      DTR = 3.141592654 / 180.
      WRITE(ITYP,1)
 1     FORMAT(/' hklm, d* and 2-theta computation'//
     1' INPUT ONLY THE NEEDED LATTICE PARAMETERS FOR YOUR SYMMETRY'/
     1'$a, b, c, alpha, beta, gamma: ')
      READ(IREAD,50)A,B,C,AL,BE,GA
 50	FORMAT(6G10.5)
      write(ityp,2)
 2      format(' Enter the coefficients of the modulation vector q*'/
     1' The modulation vector is the form:'/
     1' q^* = r * a^* + s * b^* + t * c^*'/
     1' If there is no modulation then return.')
      read(iread,51)jq1, jq2, jq3
 51    format(3g10.5)
      WRITE(ITYP,3)
 3     FORMAT(/'$Input the wavelength (if other than Cu Ka1):')
      READ(IREAD,52)W
 52     format(g10.5)
      IF(W.NE.0.0)WAVE = W
      WRITE(ITYP,4)
 4     FORMAT(/'$Input lattice type, if known (P,I or F), (default=P):')
      IF(B.EQ.0.0)B = A
      IF(C.EQ.0.0)C = A
      IF(AL.EQ.0.0)AL = 90.
      IF(BE.EQ.0.0)BE = 90.
      IF(GA.EQ.0.0)GA = 90.
      READ(IREAD,53)LAT
 53    FORMAT(A)
      ILAT = 1
      IF(LAT.EQ.1HI.or.lat.eq.1hi)then
        ILAT = 2
      endif
      IF(LAT.EQ.1HF.or.lat.eq.1hf)then
        ILAT = 3
      endif
      WRITE(ITYP,5)
5	FORMAT(/' Do you want to input a file of d values to',
     1' to compare'/' to the hkl''s calculated from the lattice',
     2' parameters?'/'$Enter <<0>> no,  <1> yes: ')
	READ(IREAD,54)IFILE
C
C  OUTPUT TO A FILE OR TERMINAL?
C
      WRITE(ITYP,6)
 6     FORMAT(/'$List output on <<0>> terminal, <1> a file :')
      IO30 = 6
      READ(IREAD,54)IO
 54    FORMAT(I)
      IF(IO .EQ. 1)IO30 = 30
      IF(IO .EQ. 0)GO TO 100
      IF(IFILE .EQ. 1) GO TO 100
C
C  GET NAME OF THE output FILE if no d/I input file
C
      WRITE(ITYP,7)
 7     FORMAT(/' Enter the name of the output file: ')
      READ(IREAD,53)NFILE
 100  IF(IFILE .EQ. 0)GO TO 103
C
C  GET THE INPUT FROM THE MDI INPUT d-I FILE FORMAT
C
*        call mdiin
*      subroutine mdiin
       write(ityp,8) 
 8      format(' The d/I input format is an mdi file.'/
     1' Enter the input filename. <.mdi>')
       read(ityp,53)nfile
       close(unit=10)
       open(unit=10, file=nfile, defaultfile='.mdi', 
     1 form='formatted', status='old', err=197)
       read(10,53,err=198)title
       read(10,55,err=198)nn
 55     format(i4)
       do i=1,nn
         read(10,56,err=198)do(i), qint(i)
         inten(i) = int(qint(i))
       enddo
 56    format(8x,f11.7,8x,f5.0)
       close(unit=10)
       goto 199
 197   write(6,9)
 9      format(' error opening d/I file.')
       go to 199
 198   write(6,10)
 10     format(' error reading d/I file.')
 199   IF(IO .EQ. 0)GO TO 105
       DO I=1,13
         IF(NFILE(I:I) .EQ. '.')GO TO 101
         if(nfile(i:i) .eq. ' ')go to 102
       enddo
 101   NFILE(I+1:I+3) = 'HKL'
       WRITE(ITYP,11)NFILE(1:i+3)
       l = i+3
       go to  104
 102   nfile(i:i+3) = '.hkl'
       l = i+3
 104   write(ityp,11)nfile(1:l)
 11     FORMAT(/'$Enter the name of the output file (default=',A,'):')
       READ(IREAD,53)MFILE
       IF(MFILE(1:4).NE. '    ')NFILE = MFILE
C
C  OPEN THE OutPUT FiLE
C
 103   IF(IO .EQ. 0)GO TO 105
       OPEN(UNIT=IO30, NAME=NFILE, TYPE='NEW',FORM='FORMATTED')
C
C  GET LIMITS ON hkl'S
C
 105  WRITE(ITYP,12)
 12    FORMAT(/' The program will generate all hkl''s for your',
     1' symmetry using integers'/' through 8.  This could produce',
     2' over 4000 hkl''s.'/' To reduce the output enter lower limits',
     3' on the hkl scan range.'//' Enter hmin,hmax,kmin,kmax,lmin,lmax',
     3',mmin,mmax',
     4' seperated by commas'/' (Defaults are integers ranging from 8 to'
     5,' 0 or -8 depending on your symmetry')
      READ(IREAD,57)HMIN,HMAX,KMIN,KMAX,LMIN,LMAX,mmin,mmax
 57    FORMAT(8I)
      write(ityp,30)
 30     format(' Enter the two-theta angle limit for'/
     1' calculated values. (default=75 degrees)')
      read(iread,60)alim
 60     format(g10.5)
      if(alim.eq.0)then
         alim = 75.0
      endif
      WRITE(IO30,13)A,B,C,AL,BE,GA
 13    FORMAT(//' Real Cell Dimensions'/
     1' a = ',F10.5,' b = ',F10.5,' c = ',F10.5/
     2' alpha= ',F10.3,' beta= ',F10.3,' gamma= ',F10.3)
      AL = AL * DTR
      BE = BE * DTR
      GA = GA * DTR
      S = AL+BE+GA
      S=S/2.0
      V=2.0*A*B*C*SQRT(SIN(S)*SIN(S-AL)*SIN(S-
     1BE)*SIN(S-GA))
      WRITE(IO30,14) V, LAT
 14    FORMAT(' Cell Volume = ',F10.5, 'cubic Angstroms',
     1'  Lattice Type = ',A1)
      AS = B*C*SIN(AL)/V
      BS = A*C*SIN(BE)/V
      CS = A * B* SIN(GA)/V
      V = 1./V
      AL = 180.- AL/DTR
      BE = 180.- BE/DTR
      GA = 180.- GA/DTR
      WRITE(IO30,15)AS,BS,CS,AL,BE,GA,V
 15    FORMAT(/' Reciprocal Cell Dimensions'/
     1' a* = ',F10.6,' b* = ',F10.6,' c* = ',F10.6/
     2' alpha* = ',F10.3,' beta* = ',F10.3,'gamma* = ',F10.6/
     3' Reciprocal Volume = ',G20.8,' A-1'//)
      AL = AL * DTR
      BE = BE * DTR
      GA = GA * DTR
      if(jq1.ne.0.or.jq2.ne.0.or.jq3.ne.0)then
        write(io30,31)jq1,jq2,jq3
 31      format(' The modulated wave vector is:'/
     1 3x'q* = ',f5.3,'a* + ',f5.3,'b* + ',f5.3,'c*')    
      endif
*
* calculate the mod vector values
*
      if(jq1.ne.0.or.jq2.ne.0.or.jq3.ne.0)then
        qstsq = jq1**2 * as**2 + jq2**2 * bs**2 + jq3**2 * cs**2 +
     1   2*jq1*jq2*as*bs*cos(ga) + 2*jq1*jq3*as*cs*cos(be) +
     1   2*jq2*jq3*bs*cs*cos(al)
        adq = jq1*as**2 + jq2*as*bs*cos(ga) + jq3*as*cs*cos(be)
        bdq = jq1*as*bs*cos(ga) + jq2*bs**2 + jq3*bs*cs*cos(al)
        cdq = jq1*as*cs*cos(be) + jq2*bs*cs*cos(al) + jq3*cs**2
      else
        qstsq = 0.0
        adq = 0.0
        bdq = 0.0
        cdq = 0.0
      endif
*
      if(jq1.ne.0.or.jq2.ne.0.or.jq3.ne.0)then
        if(ifile.eq.0)then
          WRITE(IO30,16)
        else
          write(io30,17)
        endif
      else
        if(ifile.eq.0)then
          write(io30,18)
        else
          write(io30,19)
        endif
      endif
 16   FORMAT(//'   h  k  l  m     d       2*Theta')
 18   format(//'   h  k  l    d        2*theta')
 17   FORMAT(//'   h  k  l  m   2*Theta      d       d-obs    I-obs')
 19   format(//'   h  k  l   2*Theta      d        d-obs   I-obs')
C
C  SET LIMITS FOR DIFFERENT CRYSTAL SYSTEMS
C
C  FIRST TRICLINIC
C
	IF(HMIN .NE. 0 .OR. HMAX .NE. 0)GO TO 890
	HMIN = 0
	HMAX = 8
	KMIN = 0
	KMAX = 8
	IF(ABS(AL-90.) .GT. .00001)KMIN=-8
C
C  NOW MONOCLINIC
C
	LMIN = 0
	LMAX =20
	IF(ABS(BE-90.) .GT. .00001 .OR. ABS(GA-90.) .GT. 
     1 .00001 )LMIN = -8
	N = 0
 890  line = 1
      if(jq1.eq.0.and.jq2.eq.0.and.jq3.eq.0)goto 896
      do 900 m=mmin,mmax
 896  DO 900 L=LMIN,LMAX
      DO 900 K=KMIN,KMAX
      DO 900 J=HMIN,HMAX
      IF(J .EQ.0.AND.K.EQ.0.AND.L.EQ.0)GO TO 900
	GO TO(850,830,840),ILAT
C
C  TEST FOR BODY CENTERING
C
830	H = JMOD(J+K+L,2)
	IF(H.EQ.0)GO TO 850
	GO TO 900
C
C  TEST FOR F CENTERING CONDITION
C
840	H = JMOD(J,2) + JMOD(K,2) + JMOD(L,2)
	IF(H.EQ.0.OR.H.EQ.3)GO TO 850
	GO TO 900
850      DSs = J*J*AS*AS + K*K*BS*BS + L*L*CS*CS + m*m*qstsq +
     1 + 2*J*K*AS*BS*COS(GA) + 2*j*m*adq
     2 + 2*J*L*AS*CS*COS(BE) + 2*k*m*bdq
     3 + 2*K*L*BS*CS*COS(AL) + 2*l*m*cdq
        DS = SQRT(DSs)
        D(line) = 1.0/DS
        ind(line,1) = j
        ind(line,2) = k
        ind(line,3) = l
        ind(line,4) = m
        T2 = WAVE / (2.0 * D(line))
        IF(T2.GT.1.0)GO TO 900
        THET(line) = ASIN(T2) * 2.0 / DTR
        line = line + 1
 900  CONTINUE
      line = line - 1
*
* sort and reorder the d's
*
910   do i=1,line-1
        max = i
        do j = i+1, line
          if(d(max).lt.d(j))max = j
        enddo
* assign sorted d value
        td = d(i)
        d(i) = d(max)
        d(max) = td
* assign sorted indicie values
        do k=1,4
          indt = ind(i,k)
          ind(i,k) = ind(max,k)
          ind(max,k) = indt
        enddo
* assign sorted two-theta value
        tthet = thet(i)
        thet(i) = thet(max)
        thet(max) = tthet
      enddo
* write out the sorted values to a file
C
C  TEST IF D FALLS WITHIN -/+ .02 OF CURRENT D
C
       icount = 0
       if(jq1.ne.0.or.jq2.ne.0.or.jq3.ne.0)then
       do 250 i=1,line
         if(thet(i).gt.alim)goto 400
         do 200 j=1,nn
           IF(ABS(DO(j) - D(i)) .lT. .02)then
            WRITE(IO30,22)(ind(i,k),k=1,4),THET(i),D(i),DO(j),inten(j)
            icount = icount + 1
            go to 250
           endif
 200     continue
         WRITE(Io30,20)(ind(i,k),k=1,4),THET(i),D(i)
 250   continue
      else
       do 350 i=1,line
        if(thet(i).gt.alim)goto 400
        do 300 j=1,nn
          IF(ABS(DO(j) - D(i)) .lT. .02)then
            WRITE(IO30,23)(ind(i,k),k=1,3),THET(i),D(i),DO(j),inten(j)
            icount = icount + 1
            go to 350
          endif
 300    continue
        WRITE(Io30,21)(ind(i,k),k=1,3),THET(i),D(i)
 350   continue
      endif
 20   FORMAT( 1X,4I3,2F10.5)
 21   FORMAT( 1X,3I3,2F10.5)
 22   FORMAT( 1x,4i3,3f10.5,i8)
 23   FORMAT( 1x,3i3,3f10.5,i8)
 400  close(unit=io30)
      if(ifile.eq.1.and.io.eq.1)then
        write(ityp,24)icount,nn,nfile(1:l)
 24      format(i3,' out of ',i3,' matching d values have been'/
     1' found and written to ',a)
      endif
      stop 
      END
