
                                                                17-Feb-1993 14:31:35	VAX FORTRAN V5.3-50                 Page   1
                                                                 6-Jun-1990 17:43:53	POWD:[PUBLIC.SOURCE.GENERAL]HKLCAL.FOR;1    

0001    *   program hklmcal
0002    *
0003    C  PROGRAM TO COMPUTE 2-THETA'S AND HKLm'S FROM LATTICE PARAMETERS
0004    *  and modulation vector q*
0005    C
0006    C  PROGRAMMED BY: R. L. SNYDER
0007    C		  N.Y.S. COLLEGE OF CERAMICS
0008    C		  ALFRED UNIVERSITY
0009    C		  ALFRED, N.Y. 14802
0010    *  and modified to handle modulated structures by:
0011    *                 David P. Matheis
0012    *                 May 30, 1990
0013    C
0014          INTEGER H, HMAX, HMIN, inten(200)
0015          real jq1, jq2, jq3, DO(200), qINT(200)
0016          CHARACTER*75 NFILE, MFILE, title
0017          dimension d(2000), thet(2000), ind(2000,4)
0018          DATA WAVE / 1.5405981/, IREAD/5/, ITYP/6/
0019          HMIN = 0
0020          HMAX = 0
0021          DTR = 3.141592654 / 180.
0022          WRITE(ITYP,1)
0023     1     FORMAT(/' hklm, d* and 2-theta computation'//
0024         1' INPUT ONLY THE NEEDED LATTICE PARAMETERS FOR YOUR SYMMETRY'/
0025         1'$a, b, c, alpha, beta, gamma: ')
0026          READ(IREAD,50)A,B,C,AL,BE,GA
0027     50	FORMAT(6G10.5)
0028          write(ityp,2)
0029     2      format(' Enter the coefficients of the modulation vector q*'/
0030         1' The modulation vector is the form:'/
0031         1' q^* = r * a^* + s * b^* + t * c^*'/
0032         1' If there is no modulation then return.')
0033          read(iread,51)jq1, jq2, jq3
0034     51    format(3g10.5)
0035          WRITE(ITYP,3)
0036     3     FORMAT(/'$Input the wavelength (if other than Cu Ka1):')
0037          READ(IREAD,52)W
0038     52     format(g10.5)
0039          IF(W.NE.0.0)WAVE = W
0040          WRITE(ITYP,4)
0041     4     FORMAT(/'$Input lattice type, if known (P,I or F), (default=P):')
0042          IF(B.EQ.0.0)B = A
0043          IF(C.EQ.0.0)C = A
0044          IF(AL.EQ.0.0)AL = 90.
0045          IF(BE.EQ.0.0)BE = 90.
0046          IF(GA.EQ.0.0)GA = 90.
0047          READ(IREAD,53)LAT
0048     53    FORMAT(A)
0049          ILAT = 1
0050          IF(LAT.EQ.1HI.or.lat.eq.1hi)then
0051            ILAT = 2
0052          endif
0053          IF(LAT.EQ.1HF.or.lat.eq.1hf)then
0054            ILAT = 3
0055          endif
0056          WRITE(ITYP,5)
0057    5	FORMAT(/' Do you want to input a file of d values to',

HKLCAL$MAIN                                                     17-Feb-1993 14:31:35	VAX FORTRAN V5.3-50                 Page   2
                                                                 6-Jun-1990 17:43:53	POWD:[PUBLIC.SOURCE.GENERAL]HKLCAL.FOR;1    

0058         1' to compare'/' to the hkl''s calculated from the lattice',
0059         2' parameters?'/'$Enter <<0>> no,  <1> yes: ')
0060    	READ(IREAD,54)IFILE
0061    C
0062    C  OUTPUT TO A FILE OR TERMINAL?
0063    C
0064          WRITE(ITYP,6)
0065     6     FORMAT(/'$List output on <<0>> terminal, <1> a file :')
0066          IO30 = 6
0067          READ(IREAD,54)IO
0068     54    FORMAT(I)
0069          IF(IO .EQ. 1)IO30 = 30
0070          IF(IO .EQ. 0)GO TO 100
0071          IF(IFILE .EQ. 1) GO TO 100
0072    C
0073    C  GET NAME OF THE output FILE if no d/I input file
0074    C
0075          WRITE(ITYP,7)
0076     7     FORMAT(/' Enter the name of the output file: ')
0077          READ(IREAD,53)NFILE
0078     100  IF(IFILE .EQ. 0)GO TO 103
0079    C
0080    C  GET THE INPUT FROM THE MDI INPUT d-I FILE FORMAT
0081    C
0082    *        call mdiin
0083    *      subroutine mdiin
0084           write(ityp,8) 
0085     8      format(' The d/I input format is an mdi file.'/
0086         1' Enter the input filename. <.mdi>')
0087           read(ityp,53)nfile
0088           close(unit=10)
0089           open(unit=10, file=nfile, defaultfile='.mdi', 
0090         1 form='formatted', status='old', err=197)
0091           read(10,53,err=198)title
0092           read(10,55,err=198)nn
0093     55     format(i4)
0094           do i=1,nn
0095             read(10,56,err=198)do(i), qint(i)
0096             inten(i) = int(qint(i))
0097           enddo
0098     56    format(8x,f11.7,8x,f5.0)
0099           close(unit=10)
0100           goto 199
0101     197   write(6,9)
0102     9      format(' error opening d/I file.')
0103           go to 199
0104     198   write(6,10)
0105     10     format(' error reading d/I file.')
0106     199   IF(IO .EQ. 0)GO TO 105
0107           DO I=1,13
0108             IF(NFILE(I:I) .EQ. '.')GO TO 101
0109             if(nfile(i:i) .eq. ' ')go to 102
0110           enddo
0111     101   NFILE(I+1:I+3) = 'HKL'
0112           WRITE(ITYP,11)NFILE(1:i+3)
0113           l = i+3
0114           go to  104

HKLCAL$MAIN                                                     17-Feb-1993 14:31:35	VAX FORTRAN V5.3-50                 Page   3
                                                                 6-Jun-1990 17:43:53	POWD:[PUBLIC.SOURCE.GENERAL]HKLCAL.FOR;1    

0115     102   nfile(i:i+3) = '.hkl'
0116           l = i+3
0117     104   write(ityp,11)nfile(1:l)
0118     11     FORMAT(/'$Enter the name of the output file (default=',A,'):')
0119           READ(IREAD,53)MFILE
0120           IF(MFILE(1:4).NE. '    ')NFILE = MFILE
0121    C
0122    C  OPEN THE OutPUT FiLE
0123    C
0124     103   IF(IO .EQ. 0)GO TO 105
0125           OPEN(UNIT=IO30, NAME=NFILE, TYPE='NEW',FORM='FORMATTED')
0126    C
0127    C  GET LIMITS ON hkl'S
0128    C
0129     105  WRITE(ITYP,12)
0130     12    FORMAT(/' The program will generate all hkl''s for your',
0131         1' symmetry using integers'/' through 8.  This could produce',
0132         2' over 4000 hkl''s.'/' To reduce the output enter lower limits',
0133         3' on the hkl scan range.'//' Enter hmin,hmax,kmin,kmax,lmin,lmax',
0134         3',mmin,mmax',
0135         4' seperated by commas'/' (Defaults are integers ranging from 8 to'
0136         5,' 0 or -8 depending on your symmetry')
0137          READ(IREAD,57)HMIN,HMAX,KMIN,KMAX,LMIN,LMAX,mmin,mmax
0138     57    FORMAT(8I)
0139          write(ityp,30)
0140     30     format(' Enter the two-theta angle limit for'/
0141         1' calculated values. (default=75 degrees)')
0142          read(iread,60)alim
0143     60     format(g10.5)
0144          if(alim.eq.0)then
0145             alim = 75.0
0146          endif
0147          WRITE(IO30,13)A,B,C,AL,BE,GA
0148     13    FORMAT(//' Real Cell Dimensions'/
0149         1' a = ',F10.5,' b = ',F10.5,' c = ',F10.5/
0150         2' alpha= ',F10.3,' beta= ',F10.3,' gamma= ',F10.3)
0151          AL = AL * DTR
0152          BE = BE * DTR
0153          GA = GA * DTR
0154          S = AL+BE+GA
0155          S=S/2.0
0156          V=2.0*A*B*C*SQRT(SIN(S)*SIN(S-AL)*SIN(S-
0157         1BE)*SIN(S-GA))
0158          WRITE(IO30,14) V, LAT
0159     14    FORMAT(' Cell Volume = ',F10.5, 'cubic Angstroms',
0160         1'  Lattice Type = ',A1)
0161          AS = B*C*SIN(AL)/V
0162          BS = A*C*SIN(BE)/V
0163          CS = A * B* SIN(GA)/V
0164          V = 1./V
0165          AL = 180.- AL/DTR
0166          BE = 180.- BE/DTR
0167          GA = 180.- GA/DTR
0168          WRITE(IO30,15)AS,BS,CS,AL,BE,GA,V
0169     15    FORMAT(/' Reciprocal Cell Dimensions'/
0170         1' a* = ',F10.6,' b* = ',F10.6,' c* = ',F10.6/
0171         2' alpha* = ',F10.3,' beta* = ',F10.3,'gamma* = ',F10.6/

HKLCAL$MAIN                                                     17-Feb-1993 14:31:35	VAX FORTRAN V5.3-50                 Page   4
                                                                 6-Jun-1990 17:43:53	POWD:[PUBLIC.SOURCE.GENERAL]HKLCAL.FOR;1    

0172         3' Reciprocal Volume = ',G20.8,' A-1'//)
0173          AL = AL * DTR
0174          BE = BE * DTR
0175          GA = GA * DTR
0176          if(jq1.ne.0.or.jq2.ne.0.or.jq3.ne.0)then
0177            write(io30,31)jq1,jq2,jq3
0178     31      format(' The modulated wave vector is:'/
0179         1 3x'q* = ',f5.3,'a* + ',f5.3,'b* + ',f5.3,'c*')    
0180          endif
0181    *
0182    * calculate the mod vector values
0183    *
0184          if(jq1.ne.0.or.jq2.ne.0.or.jq3.ne.0)then
0185            qstsq = jq1**2 * as**2 + jq2**2 * bs**2 + jq3**2 * cs**2 +
0186         1   2*jq1*jq2*as*bs*cos(ga) + 2*jq1*jq3*as*cs*cos(be) +
0187         1   2*jq2*jq3*bs*cs*cos(al)
0188            adq = jq1*as**2 + jq2*as*bs*cos(ga) + jq3*as*cs*cos(be)
0189            bdq = jq1*as*bs*cos(ga) + jq2*bs**2 + jq3*bs*cs*cos(al)
0190            cdq = jq1*as*cs*cos(be) + jq2*bs*cs*cos(al) + jq3*cs**2
0191          else
0192            qstsq = 0.0
0193            adq = 0.0
0194            bdq = 0.0
0195            cdq = 0.0
0196          endif
0197    *
0198          if(jq1.ne.0.or.jq2.ne.0.or.jq3.ne.0)then
0199            if(ifile.eq.0)then
0200              WRITE(IO30,16)
0201            else
0202              write(io30,17)
0203            endif
0204          else
0205            if(ifile.eq.0)then
0206              write(io30,18)
0207            else
0208              write(io30,19)
0209            endif
0210          endif
0211     16   FORMAT(//'   h  k  l  m     d       2*Theta')
0212     18   format(//'   h  k  l    d        2*theta')
0213     17   FORMAT(//'   h  k  l  m   2*Theta      d       d-obs    I-obs')
0214     19   format(//'   h  k  l   2*Theta      d        d-obs   I-obs')
0215    C
0216    C  SET LIMITS FOR DIFFERENT CRYSTAL SYSTEMS
0217    C
0218    C  FIRST TRICLINIC
0219    C
0220    	IF(HMIN .NE. 0 .OR. HMAX .NE. 0)GO TO 890
0221    	HMIN = 0
0222    	HMAX = 8
0223    	KMIN = 0
0224    	KMAX = 8
0225    	IF(ABS(AL-90.) .GT. .00001)KMIN=-8
0226    C
0227    C  NOW MONOCLINIC
0228    C

HKLCAL$MAIN                                                     17-Feb-1993 14:31:35	VAX FORTRAN V5.3-50                 Page   5
                                                                 6-Jun-1990 17:43:53	POWD:[PUBLIC.SOURCE.GENERAL]HKLCAL.FOR;1    

0229    	LMIN = 0
0230    	LMAX =20
0231    	IF(ABS(BE-90.) .GT. .00001 .OR. ABS(GA-90.) .GT. 
0232         1 .00001 )LMIN = -8
0233    	N = 0
0234     890  line = 1
0235          if(jq1.eq.0.and.jq2.eq.0.and.jq3.eq.0)goto 896
0236          do 900 m=mmin,mmax
0237     896  DO 900 L=LMIN,LMAX
0238          DO 900 K=KMIN,KMAX
0239          DO 900 J=HMIN,HMAX
0240          IF(J .EQ.0.AND.K.EQ.0.AND.L.EQ.0)GO TO 900
0241    	GO TO(850,830,840),ILAT
0242    C
0243    C  TEST FOR BODY CENTERING
0244    C
0245    830	H = JMOD(J+K+L,2)
0246    	IF(H.EQ.0)GO TO 850
0247    	GO TO 900
0248    C
0249    C  TEST FOR F CENTERING CONDITION
0250    C
0251    840	H = JMOD(J,2) + JMOD(K,2) + JMOD(L,2)
0252    	IF(H.EQ.0.OR.H.EQ.3)GO TO 850
0253    	GO TO 900
0254    850      DSs = J*J*AS*AS + K*K*BS*BS + L*L*CS*CS + m*m*qstsq +
0255         1 + 2*J*K*AS*BS*COS(GA) + 2*j*m*adq
0256         2 + 2*J*L*AS*CS*COS(BE) + 2*k*m*bdq
0257         3 + 2*K*L*BS*CS*COS(AL) + 2*l*m*cdq
0258            DS = SQRT(DSs)
0259            D(line) = 1.0/DS
0260            ind(line,1) = j
0261            ind(line,2) = k
0262            ind(line,3) = l
0263            ind(line,4) = m
0264            T2 = WAVE / (2.0 * D(line))
0265            IF(T2.GT.1.0)GO TO 900
0266            THET(line) = ASIN(T2) * 2.0 / DTR
0267            line = line + 1
0268     900  CONTINUE
0269          line = line - 1
0270    *
0271    * sort and reorder the d's
0272    *
0273    910   do i=1,line-1
0274            max = i
0275            do j = i+1, line
0276              if(d(max).lt.d(j))max = j
0277            enddo
0278    * assign sorted d value
0279            td = d(i)
0280            d(i) = d(max)
0281            d(max) = td
0282    * assign sorted indicie values
0283            do k=1,4
0284              indt = ind(i,k)
0285              ind(i,k) = ind(max,k)

HKLCAL$MAIN                                                     17-Feb-1993 14:31:35	VAX FORTRAN V5.3-50                 Page   6
                                                                 6-Jun-1990 17:43:53	POWD:[PUBLIC.SOURCE.GENERAL]HKLCAL.FOR;1    

0286              ind(max,k) = indt
0287            enddo
0288    * assign sorted two-theta value
0289            tthet = thet(i)
0290            thet(i) = thet(max)
0291            thet(max) = tthet
0292          enddo
0293    * write out the sorted values to a file
0294    C
0295    C  TEST IF D FALLS WITHIN -/+ .02 OF CURRENT D
0296    C
0297           icount = 0
0298           if(jq1.ne.0.or.jq2.ne.0.or.jq3.ne.0)then
0299           do 250 i=1,line
0300             if(thet(i).gt.alim)goto 400
0301             do 200 j=1,nn
0302               IF(ABS(DO(j) - D(i)) .lT. .02)then
0303                WRITE(IO30,22)(ind(i,k),k=1,4),THET(i),D(i),DO(j),inten(j)
0304                icount = icount + 1
0305                go to 250
0306               endif
0307     200     continue
0308             WRITE(Io30,20)(ind(i,k),k=1,4),THET(i),D(i)
0309     250   continue
0310          else
0311           do 350 i=1,line
0312            if(thet(i).gt.alim)goto 400
0313            do 300 j=1,nn
0314              IF(ABS(DO(j) - D(i)) .lT. .02)then
0315                WRITE(IO30,23)(ind(i,k),k=1,3),THET(i),D(i),DO(j),inten(j)
0316                icount = icount + 1
0317                go to 350
0318              endif
0319     300    continue
0320            WRITE(Io30,21)(ind(i,k),k=1,3),THET(i),D(i)
0321     350   continue
0322          endif
0323     20   FORMAT( 1X,4I3,2F10.5)
0324     21   FORMAT( 1X,3I3,2F10.5)
0325     22   FORMAT( 1x,4i3,3f10.5,i8)
0326     23   FORMAT( 1x,3i3,3f10.5,i8)
0327     400  close(unit=io30)
0328          if(ifile.eq.1.and.io.eq.1)then
0329            write(ityp,24)icount,nn,nfile(1:l)
0330     24      format(i3,' out of ',i3,' matching d values have been'/
0331         1' found and written to ',a)
0332          endif
0333          stop 
0334          END

HKLCAL$MAIN                                                     17-Feb-1993 14:31:35	VAX FORTRAN V5.3-50                 Page   7
01                                                               6-Jun-1990 17:43:53	POWD:[PUBLIC.SOURCE.GENERAL]HKLCAL.FOR;1    

PROGRAM SECTIONS

    Name				 Bytes   Attributes

  0 $CODE                                 3763   PIC CON REL LCL   SHR   EXE   RD NOWRT LONG
  1 $PDATA                                1995   PIC CON REL LCL   SHR NOEXE   RD NOWRT LONG
  2 $LOCAL                               50944   PIC CON REL LCL NOSHR NOEXE   RD   WRT LONG

    Total Space Allocated                56702


ENTRY POINTS

    Address  Type  Name          

  0-00000000       HKLCAL$MAIN   


VARIABLES

    Address  Type  Name              Address  Type  Name              Address  Type  Name              Address  Type  Name          

  2-0000C5DC  R*4  A               2-0000C638  R*4  ADQ             2-0000C5E8  R*4  AL              2-0000C630  R*4  ALIM          
      **      R*4  AS              2-0000C5E0  R*4  B               2-0000C63C  R*4  BDQ             2-0000C5EC  R*4  BE            
      **      R*4  BS              2-0000C5E4  R*4  C               2-0000C640  R*4  CDQ                 **      R*4  CS            
      **      R*4  DS                  **      R*4  DSS                 **      R*4  DTR             2-0000C5F0  R*4  GA            
      **      I*4  H               2-0000C5C4  I*4  HMAX            2-0000C5C8  I*4  HMIN                **      I*4  I             
      **      I*4  ICOUNT          2-0000C600  I*4  IFILE           2-0000C5FC  I*4  ILAT                **      I*4  INDT          
  2-0000C608  I*4  IO              2-0000C604  I*4  IO30                **      I*4  IREAD               **      I*4  ITYP          
      **      I*4  J               2-0000C5CC  R*4  JQ1             2-0000C5D0  R*4  JQ2             2-0000C5D4  R*4  JQ3           
      **      I*4  K               2-0000C61C  I*4  KMAX            2-0000C618  I*4  KMIN            2-0000C614  I*4  L             
  2-0000C5F8  I*4  LAT                 **      I*4  LINE            2-0000C624  I*4  LMAX            2-0000C620  I*4  LMIN          
  2-0000C644  I*4  M                   **      I*4  MAX             2-0000C52B  CHAR MFILE           2-0000C62C  I*4  MMAX          
  2-0000C628  I*4  MMIN                **      I*4  N               2-0000C4E0  CHAR NFILE           2-0000C60C  I*4  NN            
  2-0000C634  R*4  QSTSQ               **      R*4  S                   **      R*4  T2                  **      R*4  TD            
  2-0000C576  CHAR TITLE               **      R*4  TTHET               **      R*4  V               2-0000C5F4  R*4  W             
  2-0000C5D8  R*4  WAVE          


ARRAYS

    Address  Type  Name                Bytes  Dimensions

  2-00000960  R*4  D                    8000  (2000)
  2-00000320  R*4  DO                    800  (200)
  2-000047E0  I*4  IND                 32000  (2000, 4)
  2-00000000  I*4  INTEN                 800  (200)
  2-00000640  R*4  QINT                  800  (200)
  2-000028A0  R*4  THET                 8000  (2000)

HKLCAL$MAIN                                                     17-Feb-1993 14:31:35	VAX FORTRAN V5.3-50                 Page   8
01                                                               6-Jun-1990 17:43:53	POWD:[PUBLIC.SOURCE.GENERAL]HKLCAL.FOR;1    

LABELS

    Address   Label       Address   Label       Address   Label       Address   Label       Address   Label       Address   Label   

  1-00000005  1'        1-00000090  2'        1-00000141  3'        1-00000176  4'        1-000001B2  5'        1-00000246  6'      
  1-00000278  7'        1-000002A0  8'        1-000002FA  9'        1-00000315  10'       1-00000330  11'       1-00000365  12'     
  1-0000051E  13'       1-0000057E  14'       1-000005B9  15'       1-0000068F  16'       1-000006D8  17'       1-000006B5  18'     
  1-00000710  19'       1-00000745  20'       1-00000751  21'       1-0000075D  22'       1-0000076B  23'       1-00000779  24'     
  1-000004C8  30'       1-00000649  31'       1-0000008A  50'       1-0000013B  51'       1-00000172  52'       1-000001B0  53'     
  1-00000276  54'       1-000002EC  55'       1-000002EF  56'       1-000004C4  57'       1-0000051A  60'       0-0000026A  100     
  0-000003D8  101       0-000003C4  102       0-0000047C  103       0-00000417  104       0-0000048A  105       0-00000384  197     
  0-0000036C  198       0-0000039A  199           **      200       0-00000D44  250           **      300       0-00000E50  350     
  0-00000E56  400       0-000009FF  830       0-00000A18  840       0-00000A48  850       0-00000989  890       0-000009B4  896     
  0-00000B6E  900           **      910     


FUNCTIONS AND SUBROUTINES REFERENCED

  Type  Name            Type  Name            Type  Name            Type  Name            Type  Name            Type  Name          

        FOR$CLOSE             FOR$OPEN         R*4  MTH$ASIN         R*4  MTH$COS          R*4  MTH$SIN          R*4  MTH$SQRT      


COMMAND QUALIFIERS

  FORT/NOLIST/LIS HKLCAL.FOR

  /CHECK=(NOBOUNDS,OVERFLOW,NOUNDERFLOW)
  /DEBUG=(NOSYMBOLS,TRACEBACK)
  /SHOW=(NODICTIONARY,NOINCLUDE,MAP,NOPREPROCESSOR,SINGLE)
  /STANDARD=(NOSEMANTIC,NOSOURCE_FORM,NOSYNTAX)
  /WARNINGS=(NODECLARATIONS,GENERAL,NOULTRIX,NOVAXELN)
  /CONTINUATIONS=19  /NOCROSS_REFERENCE  /NOD_LINES  /NOEXTEND_SOURCE
  /F77  /NOG_FLOATING  /I4  /NOMACHINE_CODE  /OPTIMIZE  /NOPARALLEL
  /NOANALYSIS_DATA
  /NODIAGNOSTICS
  /LIST=POWD:[PUBLIC.SOURCE.GENERAL]HKLCAL.LIS;1                                                            
  /OBJECT=POWD:[PUBLIC.SOURCE.GENERAL]HKLCAL.OBJ;1                                                            


COMPILATION STATISTICS

  Run Time:           4.65 seconds
  Elapsed Time:       6.68 seconds
  Page Faults:        607
  Dynamic Memory:     820 pages
