source: trunk/fsource/DIFFaXsubs/DIFFaXsubs.for @ 2205

Last change on this file since 2205 was 2205, checked in by vondreele, 6 years ago

implement sequential stacking fault calculations with plotting
fix PlotXY

File size: 193.4 KB
Line 
1************************************************************************
2*                                                                      *
3*     Copyright 1987-2010 Michael M. J. Treacy and Michael W. Deem     *
4*                                                                      *
5************************************************************************
6************************************************************************
7*****************      Source file DIFFaXsubs.for       ****************
8************************************************************************
9************************************************************************
10********************* version 1.813, 19th May, 2010 ********************
11******* Trimmed for GSAS-II use by R. Von Dreele March, 2016 ***********
12************************************************************************
13************************************************************************
14*
15* ______________________________________________________________________
16* Title: block data
17* Author: MMJT
18* Date: 23 Oct 1988
19* Description: Sets up some global constants which help make the
20* code easier to read.
21*
22*      COMMON VARIABLES:
23*
24*        modifies:  CENTRO, ELECTN, NEUTRN, NONE, GAUSS,
25*                   LORENZ, PS_VGT, PV_GSS, PV_LRN, X_RAY
26* ______________________________________________________________________
27*
28      block data
29      implicit none
30*
31      integer*4 NONE, CENTRO, GAUSS, LORENZ, PS_VGT, PV_GSS, PV_LRN,
32     |          X_RAY, NEUTRN, ELECTN
33*
34* The common block 'consts' also occurs in the file 'DIFFaX.inc'
35      common /consts/ NONE, CENTRO, GAUSS, LORENZ, PS_VGT, PV_GSS,
36     |                PV_LRN, X_RAY, NEUTRN, ELECTN
37*
38      data NONE, CENTRO /0, 1/
39      data GAUSS, LORENZ, PS_VGT, PV_GSS, PV_LRN /1, 2, 3, 4, 5/
40      data X_RAY, NEUTRN, ELECTN /0, 1, 2/
41*
42      end
43*
44* ______________________________________________________________________
45* Title: AGLQ16
46* Author: MWD
47* Date: 18 Aug 1988
48* Description:  This routine does adaptive 16-point Gauss-Legendre
49* quadrature in an interval (h,k,a) to (h,k,b) of reciprocal space.
50* ok is returned as .TRUE. if GLQ16 hasn't blown its stack. The
51* integrated result is returned in AGLQ16.
52*
53*      ARGUMENTS:
54*            h   -  reciprocal lattice vector h-component. (input).
55*            k   -  reciprocal lattice vector k-component. (input).
56*            a   -  l-value of the lower bound of reciprocal
57*                   lattice integration region. (input).
58*            b   -  l-value of the upper bound of reciprocal
59*                   lattice integration region. (input).
60*            ok  -  logical flag indicating all went well. (output).
61*
62*      AGLQ16 returns the adaptively integrated value.
63* ______________________________________________________________________
64*
65      real*8 function AGLQ16(h, k, a, b, ok)
66      include 'DIFFaX.par'
67*     save
68*
69      integer*4 h, k
70      real*8 a, b
71      logical ok
72*
73      integer*4 maxstk, stp, n, n2
74      parameter(maxstk = 200)
75      real*8 sum, sum1, sum2, sum3, epsilon, epsilon2, GLQ16,
76     |         stk(maxstk), d1, d2, d3, x
77      parameter(epsilon = FIVE * eps4)
78*
79* external function
80      external GLQ16
81*
82      AGLQ16 = ZERO
83* initalize stack; top is at highest index
84      stp = maxstk
85* get first integration
86      sum3 = GLQ16(h, k, a, b, ok)
87      if(.not.ok) goto 999
88      n2 = 2
89      n = 0
90      d1 = a
91      d3 = b
92      sum = ZERO
93   20 d2 = HALF * (d1 + d3)
94        n = n + 2
95        sum1 = GLQ16(h, k, d1, d2, ok)
96        if(.not.ok) goto 999
97        sum2 = GLQ16(h, k, d2, d3, ok)
98        if(.not.ok) goto 999
99        x = sum1 + sum2
100* determine figure of merit
101        epsilon2 = max(epsilon, abs(x) * epsilon)
102        if(abs(x - sum3).gt.epsilon2) then
103* the area of these two panels is not accurately known
104* check for stack overflow
105          if(stp.lt.3) goto 980
106* push right one onto stack
107          stk(stp) = sum2
108          stk(stp-1) = d2
109          stk(stp-2) = d3
110          stp = stp - 3
111          d3 = d2
112          sum3 = sum1
113        else
114* this panel has been accurately integrated; save its area
115          sum = sum + x
116* get next panel
117* check for stack underflow--happens when no panels left to pop off
118          if(stp.eq.maxstk) goto 30
119          d3 = stk(stp+1)
120          d1 = stk(stp+2)
121          sum3 = stk(stp+3)
122          stp = stp + 3
123        endif
124        if(n.eq.n2) then
125          n2 = n2 * 2
126        endif
127      goto 20
128   30 AGLQ16 = sum
129      return
130  980 write(op,402) 'Stack overflow in AGLQ16.'
131      return
132  999 write(op,402) 'GLQ16 returned an error to AGLQ16.'
133      write(op,403) h, k, d1, d2
134      return
135  402 format(1x, a)
136  403 format(3x,'at: h = ',i3,' k = ',i3,' l1 = ',g12.5,' l2 = ',g12.5)
137      end
138*
139* ______________________________________________________________________
140* Title: APPR_F
141* Author: MMJT
142* Date: 11 Feb 1989
143* Description: This subroutine returns polynomial approximations of f
144*  at 16 points, for use in GLQ16. The f's are returned in a MAX_L x 16
145*  array. The order of the polynomial is n-1, where n is input. The 16
146*  points in reciprocal space for which f's are needed are given by
147*  (h,k,ag_l(i)), where h, k and ag_l are input.
148*  The ll are the n sampling points, whose f's must be calculated
149*  by direct calls to GET_F, and from which the interpolations are
150*  calculated. list contains the indices of those n points. There is no
151*  need to interpolate for these values since we know them already!
152*
153*      ARGUMENTS:
154*            f      -  Array of layer form factors. (output).
155*            h      -  reciprocal lattice vector h-component. (input).
156*            k      -  reciprocal lattice vector k-component. (input).
157*            ll     -  an array of n l-values to be used in the
158*                      interpolation. (input).
159*            ag_l   -  an array of 16 l_values at which interpolated
160*                      form factors are required. (input).
161*            n      -  the order of the polynomial approximation.
162*                                                              (input).
163*            list   -  A list of the indices of the n ll points
164*                      entered. Interpolation is not needed at these
165*                      input values. (input).
166*            ok     -  logical flag indicating all went well. (output).
167*
168*      COMMON VARIABLES:
169*            uses:     a0, b0, c0, d0, n_layers
170* ______________________________________________________________________
171*
172      subroutine APPR_F(f,h,k,ll,ag_l,n,list,ok)
173      include 'DIFFaX.par'
174      include 'DIFFaX.inc'
175*
176      integer*4 n, list(n)
177      real*8 ll(n), ag_l(16)
178      integer*4 h, k
179      complex*16 f(MAX_L,16)
180      logical ok
181*
182      logical know_f
183      integer*4 i, j, m, p, max_poly
184      parameter (max_poly = 10)
185      real*8 Q2, l
186      complex*16 ff(MAX_L,max_poly), fa(MAX_L), f_ans, f_error
187*
188* external subroutines (Some compilers need them declared external)
189*      external POLINT, GET_F
190*
191* statement function
192* Q2 is the value of 1/d**2 at hkl
193      Q2(h,k,l) = h*h*a0 + k*k*b0 + l*l*c0 + h*k*d0
194*
195* sample GET_F n times
196      do 10 i = 1, n
197        call GET_F(ff(1,i), Q2(h,k,ll(i)), ll(i))
198   10 continue
199*
200* store the n sampled f values for each layer i in fa.
201* call POLINT for each layer in turn
202* do this 16 times.
203      do 20 m = 1, 16
204* check to see that we haven't just calculated f at l(m)
205        know_f = .false.
206        do 30 i = 1, n
207          if(m.eq.list(i)) then
208            p = i
209            know_f = .true.
210          endif
211   30   continue
212* if we have, then simply copy it
213        if(know_f) then
214          do 40 i = 1, n_layers
215            f(i,m) = ff(i,p)
216   40     continue
217        else
218* else, use polynomial interpolation.
219          do 50 i = 1, n_layers
220            do 60 j = 1, n
221              fa(j) = ff(i,j)
222   60       continue
223            call POLINT(ll,fa,n,ag_l(m),f_ans,f_error,ok)
224            if(.not.ok) goto 999
225            f(i,m) = f_ans
226   50     continue
227        endif
228   20 continue
229*
230      return
231  999 write(op,100) 'POLINT returned an error to APPR_F.'
232      return
233  100 format(1x, a)
234      end
235* ______________________________________________________________________
236* Title: BINPOW
237* Author: MMJT
238* Date: 18 Mar 1990
239* Description:  This function breaks down the number 'n' into its
240* binary representation. The result is stored in the global array 'pow'.
241* This is used for efficiently multiplying a square matrix by itself
242* n times if the RECURSIVE option was chosen for a finite number of
243* layers.
244*
245* n must be such that n <= RCSV_MAX+1 <= 2**(MAX_BIN+1) - 1
246*
247*      ARGUMENTS:
248*            n   -  number to binarize. (input).
249*
250*      COMMON VARIABLES:
251*            uses:  max_pow
252*
253*        modifies:  max_pow, pow
254*
255* BINPOW returns logical TRUE if no problems were encountered.
256* ______________________________________________________________________
257*
258      logical function BINPOW(n)
259      include 'DIFFaX.par'
260      include 'DIFFaX.inc'
261*
262      integer*4 n
263*
264      integer*4 i, j, itmp
265*
266      BINPOW = .false.
267*
268      itmp = n
269      max_pow = 0
270      i = MAX_BIN + 1
271*
272   10 i = i - 1
273        j = itmp / 2**(i-1)
274* j should be either 1 or zero if n was within bounds
275        if(j.eq.1) then
276          pow(i) = 1
277          itmp = itmp - 2**(i-1)
278          if(max_pow.eq.0) max_pow = i
279        else if(j.eq.0) then
280          pow(i) = 0
281        else
282          goto 999
283        endif
284      if(i.gt.1) goto 10
285*
286      BINPOW = .true.
287      return
288  999 write(op,100) 'ERROR in BINPOW: Invalid exponent ', n
289      write(op,100) ' Maximum No. of layers allowed = ', RCSV_MAX
290      write(op,100) ' Maximum supported by DIFFaX is', 2**(MAX_BIN+1)-1
291      return
292  100 format(1x, a, i8)
293      end
294*
295* ______________________________________________________________________
296* Title: BOUNDS
297* Authors: MMJT
298* Date: 24 Feb 1990
299* Description: This function translates the value x so that it lies
300* within the range 0 to 1.
301*
302*      ARGUMENTS:
303*                  x  -  real number to convert. (input).
304*
305*      COMMON VARIABLES:
306*
307*                  No COMMON variables are used
308*
309*      BOUNDS returns the translated value of x. x is not modified.
310* ______________________________________________________________________
311*
312      real*8 function BOUNDS(x)
313      include 'DIFFaX.par'
314*
315      real*8 x
316*
317      real*8 y
318*
319      y = x - int(x) + ONE
320      y = y - int(y)
321*
322* allow for rounding error
323      if((ONE - y).lt.eps5) y = ZERO
324*
325      BOUNDS = y
326*
327      return
328      end
329*
330* ______________________________________________________________________
331* Title: CHK_SYM
332* Author: MMJT
333* Date: 15 Aug 1989; 21 Jan 1995
334* Checks the user's assertions in the data file about the symmetry of
335* the diffraction pattern. The symmetry is used to define the smallest
336* possible volume in reciprocal space over which integration can
337* occur and be representative of the whole pattern. CHK_SYM gives the
338* user a crude 'goodness of fit', or tolerance, within which a random
339* sampling of intensities fit the asserted symmetry. CHK_SYM does
340* not override the user's judgment if the fit is poor, unless the cell
341* dimensions or angle are inconsistent with the requested symmetry. The
342* intention of CHK_SYM is to alert the user that the data does not
343* conform to the expected symmetry. Useful for debugging datafiles.
344*
345*      ARGUMENTS:
346*            ok  -  logical flag indicating all went well. (output).
347*
348*      COMMON VARIABLES:
349*            uses:  cell_a, cell_b, cell_gamma, max_angle, SymGrpNo,
350*                   pnt_grp, tolerance, PI, PI2, PS_VGT, RAD2DEG
351*
352*        modifies:  max_var, h_mirror, k_mirror, check_sym, SymGrpNo
353* ______________________________________________________________________
354*
355      subroutine CHK_SYM(ok)
356      include 'DIFFaX.par'
357      include 'DIFFaX.inc'
358*
359      logical ok
360*
361      logical diad, triad, tetrad
362      logical TST_ROT, TST_MIR
363      logical cell90, cell120, eq_sides
364      integer*4 GET_SYM, LENGTH, idum
365      real*8 tmp
366*
367* external functions
368      external TST_ROT, TST_MIR, GET_SYM, LENGTH
369* external subroutine (Some compilers need them declared external)
370*
371* reinitialize random numbers in RAN3
372      idum = -1
373*
374      diad   = .false.
375      triad  = .false.
376      tetrad = .false.
377      cell90 =  abs(cell_gamma - HALF*PI) .lt. HALF*PI*eps6
378      cell120 = abs(cell_gamma - PI2/THREE) .lt. PI2*eps6/THREE
379*
380* sample reciprocal space to get an idea of the sort of intensities
381* that are out there.
382* 360 degrees, no symmetry (-1)
383* there is nothing to check. The center of symmetry is a given.
384      if(SymGrpNo.eq.1) then
385        max_var = ZERO
386        goto 900
387      endif
388*
389* 180 degrees, rotation only (2/M, 1st setting)
390      if(SymGrpNo.eq.2) then
391        diad = TST_ROT(2, idum, ok)
392        if(.not.ok) goto 999
393        goto 900
394      endif
395*
396* 180 degrees, vertical mirror (2/M, 2nd setting)
397      if(SymGrpNo.eq.3) then
398        h_mirror = TST_MIR(1, idum, ok)
399        if(.not.ok) goto 999
400        tmp = max_var
401        max_var = ZERO
402        k_mirror = TST_MIR(2, idum, ok)
403        if(.not.ok) goto 999
404        max_var = max(tmp, max_var)
405        goto 900
406      endif
407*
408* 90 degrees, vertical mirrors (MMM)
409      if(SymGrpNo.eq.4) then
410        if(.not.cell90) goto 910
411        diad = TST_ROT(2, idum, ok)
412        if(.not.ok) goto 999
413        tmp = max_var
414        max_var = ZERO
415        h_mirror = TST_MIR(1, idum, ok)
416        if(.not.ok) goto 999
417        tmp = max(tmp, max_var)
418        max_var = ZERO
419        k_mirror = TST_MIR(2, idum, ok)
420        if(.not.ok) goto 999
421        max_var = max(tmp, max_var)
422        goto 900
423      endif
424*
425* the following point groups require equi-sided cells
426      eq_sides = abs(cell_a - cell_b) .le. HALF*eps6*(cell_a + cell_b)
427      if(.not.eq_sides) goto 920
428*
429*
430* 120 degrees, rotation only (-3)
431      if(SymGrpNo.eq.5) then
432        if(.not.cell120) goto 910
433        triad = TST_ROT(3, idum, ok)
434        if(.not.ok) goto 999
435        tmp = max_var
436        max_var = ZERO
437        h_mirror = TST_MIR(1, idum, ok)
438        if(.not.ok) goto 999
439        tmp = max(tmp, max_var)
440        max_var = ZERO
441        hk_mirror = TST_MIR(3, idum, ok)
442        if(.not.ok) goto 999
443        max_var = max(tmp, max_var)
444        goto 900
445      endif
446*
447* 60 degrees, vertical mirrors (-3M)
448      if(SymGrpNo.eq.6) then
449        if(.not.cell120) goto 910
450        triad = TST_ROT(3, idum, ok)
451        if(.not.ok) goto 999
452        tmp = max_var
453        max_var = ZERO
454        h_mirror = TST_MIR(1, idum, ok)
455        if(.not.ok) goto 999
456        tmp = max(tmp, max_var)
457        max_var = ZERO
458        hk_mirror = TST_MIR(3, idum, ok)
459        if(.not.ok) goto 999
460        max_var = max(tmp, max_var)
461        goto 900
462      endif
463*
464* 90 degrees, rotation (4/M)
465      if(SymGrpNo.eq.7) then
466        if(.not.cell90) goto 910
467        tetrad = TST_ROT(4, idum, ok)
468        if(.not.ok) goto 999
469        goto 900
470      endif
471*
472* 45 degrees, vertical mirrors (4/MMM)
473      if(SymGrpNo.eq.8) then
474        if(.not.cell90) goto 910
475        tetrad = TST_ROT(4, idum, ok)
476        if(.not.ok) goto 999
477        tmp = max_var
478        max_var = ZERO
479        h_mirror = TST_MIR(1, idum, ok)
480        if(.not.ok) goto 999
481          if(.not.h_mirror) then
482          tmp = max(tmp, max_var)
483          max_var = ZERO
484          hk_mirror = TST_MIR(3, idum, ok)
485          if(.not.ok) goto 999
486          endif
487        max_var = max(tmp, max_var)
488        goto 900
489      endif
490*
491* 60 degrees, rotation (6/M)
492      if(SymGrpNo.eq.9) then
493        if(.not.cell120) goto 910
494        diad = TST_ROT(2, idum, ok)
495        if(.not.ok) goto 999
496        tmp = max_var
497        max_var = ZERO
498        triad = TST_ROT(3, idum, ok)
499        if(.not.ok) goto 999
500        max_var = max(tmp, max_var)
501        goto 900
502      endif
503*
504* 30 degrees, vertical mirrors (6/MMM)
505      if(SymGrpNo.eq.10) then
506        if(.not.cell120) goto 910
507        diad = TST_ROT(2, idum, ok)
508        if(.not.ok) goto 999
509        tmp = max_var
510        max_var = ZERO
511        triad = TST_ROT(3, idum, ok)
512        if(.not.ok) goto 999
513        tmp = max(tmp, max_var)
514        max_var = ZERO
515        h_mirror = TST_MIR(1, idum, ok)
516        if(.not.ok) goto 999
517        if(.not.h_mirror) then
518          tmp = max(tmp, max_var)
519          max_var = ZERO
520          hk_mirror = TST_MIR(3, idum, ok)
521          if(.not.ok) goto 999
522        endif
523        max_var = max(tmp, max_var)
524      endif
525*
526      if(SymGrpNo.gt.10) goto 500
527*
528  900 write(op,202)
529     | 'The diffraction data fits the point group symmetry ''',
530     | pnt_grp(1:LENGTH(pnt_grp)),''''
531      if(max_var.gt.eps6 .and. max_var.le.eps1) then
532        write(op,203) '  with a tolerance of one part in ',
533     |                  nint(ONE / max_var)
534      else if(max_var.gt.eps1) then
535        write(op,204) '  with a tolerance of one part in ',
536     |                  ONE / max_var
537      else
538        write(op,100)
539     |      '  with a tolerance better than one part in a million.'
540      endif
541*
542  500 return
543*
544* The user's guess is inconsistent with cell_gamma.
545* Override the user.
546  910 write(op,200) 'The cell angle of',cell_gamma * RAD2DEG,
547     |            ' degrees,'
548      write(op,202) ' is inconsistent with point group symmetry ''',
549     |            pnt_grp(1:LENGTH(pnt_grp)),''''
550      write(op,300)
551      check_sym = .false.
552      SymGrpNo = GET_SYM(ok)
553      if(.not.ok) goto 999
554      write(op,205) pnt_grp(1:LENGTH(pnt_grp))
555      if(tolerance.gt.eps6 .and. tolerance.le.eps1) then
556        write(op,203) '  with a tolerance of one part in ',
557     |                  nint(ONE / tolerance)
558      else if(tolerance.gt.eps1) then
559        write(op,204) '  with a tolerance of one part in ',
560     |                  ONE / tolerance
561      else
562        write(op,100)
563     |      '  with a tolerance better than one part in a million.'
564      endif
565      return
566*
567* The user's guess is inconsistent with cell dimensions.
568* Override the user.
569  920 write(op,201) 'The cell a and b dimensions, ',
570     |            cell_a,' Angstroms by ',cell_b,' Angstroms,'
571      write(op,202) '   are inconsistent with point group symmetry ''',
572     |            pnt_grp(1:LENGTH(pnt_grp)),''''
573      write(op,300)
574* reset check_sym flag, since we are now evaluating from scratch
575      check_sym = .false.
576      max_var = ZERO
577      SymGrpNo = GET_SYM(ok)
578      if(.not.ok) goto 999
579      write(op,205) pnt_grp(1:LENGTH(pnt_grp))
580*
581      if(tolerance.gt.eps6 .and. tolerance.le.eps1) then
582        write(op,203) '  with a tolerance of one part in ',
583     |                  nint(ONE / tolerance)
584      else if(tolerance.gt.eps1) then
585        write(op,204) '  with a tolerance of one part in ',
586     |                  ONE / tolerance
587      else
588        write(op,100)
589     |      '  with a tolerance better than one part in a million.'
590      endif
591*
592      return
593*
594  999 write(op,100) 'ERROR in CHK_SYM'
595      return
596  100 format(1x, a)
597  200 format(1x, a, f7.3, a)
598  201 format(1x, 2(a, f7.3), a)
599  202 format(1x, 3a)
600  203 format(1x, a, i6)
601  204 format(1x, a, f3.1)
602  205 format(1x, 'Diffraction point symmetry is found to be ''',a,'''')
603  300 format(1x, 'Re-evaluating diffraction symmetry')
604      end
605*
606* ______________________________________________________________________
607* Title: CHWDTH
608* Author: MMJT
609* Date: 6 Mar 1995; 31st Oct 1996
610* Description:  This routine adds shape broadening caused by the finite
611* lateral width of the layers. This routine does not add the broadening
612* caused by finite width in the stacking direction. That is handled in
613* a different manner by the routine INTEN2 and associated routines.
614* The broadening handled here is added phenomenologically by applying a
615* Lorentzian profile to the computed intensity at each peak in the h-k
616* plane. If we are on the 00l axis the shape broadening is not
617* symmetrical. For a Lorentzian shape broadening intensity profile in
618* the h-k plane, the curvature of the Ewald sphere ensures a sharp
619* onset of intensity with a slowly decaying tail at higher l values. For
620* a symmetrical disk, the integrated intensity decays logarithmically.
621* In the event that the crystal width is different in the h, and h-
622* perpendicular directions, the disk of confusion becomes elongated
623* into a streak perpendicular to l, and the tail becomes more
624* Lorentzian-like. This is modelled in a phenomenological fashion, by
625* mixing Lorentzian and logarithmic terms in a manner that depends
626* on the ratio Wa/Wb. The costly logarithm function is avoided by
627* using its derivative.
628* When off the 00l axis, the broadening is modelled as a symmetric
629* Lorentzian whose half-width depends on the angle that the Ewald
630* sphere intercepts the disk of confusion (controlled by l). If
631* the lateral dimensions are not equal, then the half width is also
632* dependent on h and k. The Lorentzian is pre-computed in OPTIMZ to gain
633* computational speed.
634* This routine is called by GETSPC.
635*
636*      ARGUMENTS:
637*            h       -  Reciprocal space index. (input).
638*            k       -  Reciprocal space index. (input).
639*            l0      -  Lower bound of the l reciprocal space index
640*                       that is being integrated over. (input).
641*            l1      -  Lower bound of the l reciprocal space index
642*                       that is being integrated over. (input).
643*            x       -  The integrated intensity value along the
644*                       line defined by h,k,l0,l1. (input).
645*            m       -  The current index of the array 'spec'
646*                       corresponding to h,k,l0,l1. (input).
647*          max_indx  -  Maximum array value in spec that will be
648*                       accessed. (input).
649*
650*      COMMON VARIABLES:
651*            uses:      spec, brd_spec, FFACT_SIZE, formfactor, d_theta
652*                       ffact_scale, ffhkcnst, ffwdth
653*
654*        modifies:      spec, brd_spec
655* ______________________________________________________________________
656*
657      subroutine CHWDTH(h,k,l0,l1,x,m,max_indx)
658      include 'DIFFaX.par'
659      include 'DIFFaX.inc'
660*
661      integer*4 h, k, m, max_indx
662*
663      real*8 l0, l1, x
664*
665      integer*4 n, p, i, indx
666*
667      real*8 S, h_wdth, n_hw, d_hk, norm, l, scale, avg, xx, dx, tmp
668*
669* indx indexes into the arrays spec and brd_spec
670* n indexes into the array formfactor
671* p is the index of the centroid of the array formfactor
672* h_wdth contains the effective half-width of the size broadening
673* ffhkcnst depends only on h and k, and was computed in GETSPC
674* scale contains the calibration term for the formfactor array
675* d_hk is the average radius in reciprocal Angstroms that the
676* Ewald sphere of radius theta+d_theta intercepts the 00l plane
677* brd_spc is used for temporary storage.
678* We are only accessing half of the symmetrical formfactor array
679*
680* statement functions
681* S is the value of 1/d**2 at hkl
682      S(h,k,l) = h*h*a0 + k*k*b0 + l*l*c0 + h*k*d0
683*
684************************************************************************
685* special case if we are on the l axis
686* Don't bother if the intensity is negligible in the first place
687      if(h.eq.0 .and. k.eq.0 .and. x.gt.TEN*tiny_inty) then
688        l = HALF*(l0+l1)
689        d_hk = TWO*l*d_theta / (lambda*cell_c*ffwdth*ffwdth)
690        norm = ZERO
691        indx = 0
692        xx = Wa / Wb
693        if(xx.gt.ONE) xx = ONE / xx
694   10   indx = indx + 1
695          tmp = ONE / (ONE + indx*d_hk)
696* balance streak, versus disk of confusion (phenomenological treatment)
697* xx=1 means disk, xx=0 means streak. Intensity falls off more
698* rapidly for a streak
699          dx = ((ONE-xx)*sqrt(dble(indx))*tmp + xx)*tmp
700          if(m+indx-1.le.max_indx) brd_spc(m+indx-1) = dx
701          norm = norm + dx
702* eps5 is reasonable. However, it may be worth experimenting more.   
703          if(dx.lt.eps5) goto 20
704          goto 10
705   20   continue
706*
707      norm = x / norm
708      do 30 i = 0, indx - 1
709        if(m+i.le.max_indx) spec(m+i) = spec(m+i) + norm*brd_spc(m+i)
710   30 continue
711*
712* We were on the 00l axis, we can exit now
713        return
714      endif
715*
716************************************************************************
717* We are not on the l-axis. Broadening is handled differently.
718* scale relates the formfactor array indices to the spec array indices.
719* h_wdth and ffact_scale should never be zero. In case they are,
720* make scale large enough so that n.ge.p-1 and the loop below is
721* exited early
722      h_wdth = ffhkcnst / sqrt(S(h,k,HALF*(l0+l1)))
723      if(h_wdth.gt.ZERO .and. ffact_scale.gt.ZERO) then
724        scale = d_theta / (ffact_scale * h_wdth)
725      else
726        scale = FFACT_SIZE
727      endif
728*
729      p = FFACT_SIZE/2 + 1
730      norm = ONE
731      brd_spc(m) = ONE
732      indx = 0
733   40 indx = indx + 1
734        n_hw = indx*scale
735        n = int(n_hw)
736        if(n.ge.p-1) goto 50
737* linear interpolation of the pre-computed pseudo-Lorentzian
738        xx = n_hw - n
739        avg = (ONE-xx)*formfactor(p+n) + xx*formfactor(p+n+1)
740        if(m+indx.le.max_indx) brd_spc(m+indx) = avg
741        if(m-indx.gt.0)        brd_spc(m-indx) = avg
742* intensity x is being redistributed. We will need to normalize later
743        norm = norm + TWO*avg
744        goto 40
745   50 continue
746*
747      norm = x / norm
748      spec(m) = spec(m) + norm*brd_spc(m)
749      do 60 i = 1, indx - 1
750        if(m+i.le.max_indx) spec(m+i) = spec(m+i) + norm*brd_spc(m+i)
751        if(m-i.gt.0)        spec(m-i) = spec(m-i) + norm*brd_spc(m-i)
752   60 continue
753*
754      return
755      end
756*
757************************************************************************
758****************************LINPACK ROUTINES****************************
759************************************************************************
760* 
761* The following are the standard Linpack routines for solving complex
762* simultaneous equations. They were found to reduce DIFFaX run time by
763* significant amount (30% in one case) compared with the Numerical
764* Recipes routines LUDCMP and LUBKSB. The only changes are
765*                     
766*                         complex -> complex*16
767*                         real    -> real*8
768*                         real()  -> dble()
769*                         aimag   -> dimag
770*
771************************************************************************
772* ______________________________________________________________________
773* Title: CGESL (LINPACK ROUTINE)
774* Author: cleve moler, university of new mexico, argonne national lab.
775* Date: linpack. this version dated 08/14/78
776* Description:
777*     CGESL solves the complex system
778*     a * x = b  or  ctrans(a) * x = b
779*     using the factors computed by cgeco or CGEFA.
780*
781*     on entry
782*
783*        a       complex(lda, n)
784*                the output from cgeco or CGEFA.
785*
786*        lda     integer
787*                the leading dimension of the array  a .
788*
789*        n       integer
790*                the order of the matrix  a .
791*
792*        ipvt    integer(n)
793*                the pivot vector from cgeco or CGEFA.
794*
795*        b       complex(n)
796*                the right hand side vector.
797*
798*        job     integer
799*                = 0         to solve  a*x = b ,
800*                = nonzero   to solve  ctrans(a)*x = b  where
801*                            ctrans(a)  is the conjugate transpose.
802*
803*     on return
804*
805*        b       the solution vector  x .
806*
807*     error condition
808*
809*        a division by zero will occur if the input factor contains a
810*        zero on the diagonal.  technically this indicates singularity
811*        but it is often caused by improper arguments or improper
812*        setting of lda .  it will not occur if the subroutines are
813*        called correctly and if cgeco has set rcond .gt. 0.0
814*        or CGEFA has set info .eq. 0 .
815*
816*     to compute  inverse(a) * c  where  is a matrix
817*     with  p  columns
818*           call cgeco(a,lda,n,ipvt,rcond,z)
819*           if (rcond is too small) go to ...
820*           do 10 j = 1, p
821*              call CGESL(a,lda,n,ipvt,c(1,j),0)
822*        10 continue
823*
824*     subroutines and functions
825*
826*     blas CAXPY,CDOTC
827*     fortran conjg
828*
829* ______________________________________________________________________
830*
831      subroutine CGESL(a,lda,n,ipvt,b,job)
832      implicit none
833*
834      integer lda,n,ipvt(1),job
835      complex*16 a(lda,1),b(1)
836*
837      complex*16 CDOTC,t
838      integer k,kb,l,nm1
839* MMJT: external subroutine
840      external CAXPY
841* external function
842      external CDOTC
843*
844      nm1 = n - 1
845      if (job .ne. 0) go to 50
846*
847*        job = 0 , solve  a * x = b
848*        first solve  l*y = b
849*
850         if (nm1 .lt. 1) go to 30
851         do 20 k = 1, nm1
852            l = ipvt(k)
853            t = b(l)
854            if (l .eq. k) go to 10
855               b(l) = b(k)
856               b(k) = t
857   10       continue
858            call CAXPY(n-k,t,a(k+1,k),1,b(k+1),1)
859   20    continue
860   30    continue
861*
862*        now solve  u*x = y
863*
864         do 40 kb = 1, n
865            k = n + 1 - kb
866            b(k) = b(k)/a(k,k)
867            t = -b(k)
868            call CAXPY(k-1,t,a(1,k),1,b(1),1)
869   40    continue
870      go to 100
871   50 continue
872*
873*        job = nonzero, solve  ctrans(a) * x = b
874*        first solve  ctrans(u)*y = b
875*
876         do 60 k = 1, n
877            t = CDOTC(k-1,a(1,k),1,b(1),1)
878            b(k) = (b(k) - t)/conjg(a(k,k))
879   60    continue
880*
881*        now solve ctrans(l)*x = y
882*
883         if (nm1 .lt. 1) go to 90
884         do 80 kb = 1, nm1
885            k = n - kb
886            b(k) = b(k) + CDOTC(n-k,a(k+1,k),1,b(k+1),1)
887            l = ipvt(k)
888            if (l .eq. k) go to 70
889               t = b(l)
890               b(l) = b(k)
891               b(k) = t
892   70       continue
893   80    continue
894   90    continue
895  100 continue
896      return
897      end
898* ______________________________________________________________________
899* Title: CAXPY (LINPACK ROUTINE)
900* Author: jack dongarra
901* Date: linpack, 3/11/78
902* Description: constant times a vector plus a vector.
903* ______________________________________________________________________
904*
905      subroutine CAXPY(n,ca,cx,incx,cy,incy)
906      implicit none
907*
908      complex*16 cx(1),cy(1),ca
909      integer n,incx,incy
910*
911      integer i,ix,iy
912*
913      if(n.le.0)return
914      if (abs(dble(ca)) + abs(dimag(ca)) .eq. 0.0 ) return
915      if(incx.eq.1.and.incy.eq.1)go to 20
916*
917*        code for unequal increments or equal increments
918*          not equal to 1
919*
920      ix = 1
921      iy = 1
922      if(incx.lt.0)ix = (-n+1)*incx + 1
923      if(incy.lt.0)iy = (-n+1)*incy + 1
924      do 10 i = 1,n
925        cy(iy) = cy(iy) + ca*cx(ix)
926        ix = ix + incx
927        iy = iy + incy
928   10 continue
929      return
930*
931*        code for both increments equal to 1
932*
933   20 do 30 i = 1,n
934        cy(i) = cy(i) + ca*cx(i)
935   30 continue
936      return
937      end
938* ______________________________________________________________________
939* Title: CDOTC (LINPACK ROUTINE)
940* Author: jack dongarra
941* Date: linpack,  3/11/78.
942* Description:
943*     forms the dot product of two vectors, conjugating the first
944*     vector.
945* ______________________________________________________________________
946*
947      complex*16 function CDOTC(n,cx,incx,cy,incy)
948      implicit none
949*
950      complex*16 cx(1),cy(1)
951      integer incx,incy,n
952*
953      complex*16 ctemp
954      integer i,ix,iy
955*
956      ctemp = (0.0,0.0)
957      CDOTC = (0.0,0.0)
958      if(n.le.0)return
959      if(incx.eq.1.and.incy.eq.1)go to 20
960*
961*        code for unequal increments or equal increments
962*          not equal to 1
963*
964      ix = 1
965      iy = 1
966      if(incx.lt.0)ix = (-n+1)*incx + 1
967      if(incy.lt.0)iy = (-n+1)*incy + 1
968      do 10 i = 1,n
969        ctemp = ctemp + conjg(cx(ix))*cy(iy)
970        ix = ix + incx
971        iy = iy + incy
972   10 continue
973      CDOTC = ctemp
974      return
975*
976*        code for both increments equal to 1
977*
978   20 do 30 i = 1,n
979        ctemp = ctemp + conjg(cx(i))*cy(i)
980   30 continue
981      CDOTC = ctemp
982      return
983      end
984* ______________________________________________________________________
985* Title: CGEFA (LINPACK ROUTINE)
986* Author: cleve moler, university of new mexico, argonne national lab.
987* Date: linpack. this version dated 08/14/78
988* Description:
989*
990*     CGEFA factors a complex matrix by gaussian elimination.
991*
992*     CGEFA is usually called by cgeco, but it can be called
993*     directly with a saving in time if  rcond  is not needed.
994*     (time for cgeco) = (1 + 9/n)*(time for CGEFA) .
995*
996*     on entry
997*
998*        a       complex(lda, n)
999*                the matrix to be factored.
1000*
1001*        lda     integer
1002*                the leading dimension of the array  a .
1003*
1004*        n       integer
1005*                the order of the matrix  a .
1006*
1007*     on return
1008*
1009*        a       an upper triangular matrix and the multipliers
1010*                which were used to obtain it.
1011*                the factorization can be written  a = l*where
1012*                l  is a product of permutation and unit lower
1013*                triangular matrices and  is upper triangular.
1014*
1015*        ipvt    integer(n)
1016*                an integer vector of pivot indices.
1017*
1018*        info    integer
1019*                = 0  normal value.
1020*                = k  if  u(k,k) .eq. 0.0 .  this is not an error
1021*                     condition for this subroutine, but it does
1022*                     indicate that CGESL or cgedi will divide by zero
1023*                     if called.  use  rcond  in cgeco for a reliable
1024*                     indication of singularity.
1025*
1026*     subroutines and functions
1027*
1028*     blas CAXPY,CSCAL,ICAMAX
1029*     fortran abs,aimag,real
1030* ______________________________________________________________________
1031*
1032      subroutine CGEFA(a,lda,n,ipvt,info)
1033      implicit none
1034*
1035      integer lda,n,ipvt(1),info
1036      complex*16 a(lda,1)
1037*
1038      complex*16 t
1039      integer ICAMAX,j,k,kp1,l,nm1
1040*
1041      complex*16 zdum
1042      real*8 cabs1
1043*
1044* MMJT: external subroutine
1045      external CSCAL, CAXPY
1046*
1047* external function
1048      external ICAMAX
1049* statement function
1050      cabs1(zdum) = abs(dble(zdum)) + abs(dimag(zdum))
1051*
1052*     gaussian elimination with partial pivoting
1053*
1054      info = 0
1055      nm1 = n - 1
1056      if (nm1 .lt. 1) go to 70
1057      do 60 k = 1, nm1
1058         kp1 = k + 1
1059*
1060*        find l = pivot index
1061*
1062         l = ICAMAX(n-k+1,a(k,k),1) + k - 1
1063         ipvt(k) = l
1064*
1065*        zero pivot implies this column already triangularized
1066*
1067         if (cabs1(a(l,k)) .eq. 0.0e0) go to 40
1068*
1069*           interchange if necessary
1070*
1071            if (l .eq. k) go to 10
1072               t = a(l,k)
1073               a(l,k) = a(k,k)
1074               a(k,k) = t
1075   10       continue
1076*
1077*           compute multipliers
1078*
1079            t = -(1.0e0,0.0e0)/a(k,k)
1080            call CSCAL(n-k,t,a(k+1,k),1)
1081*
1082*           row elimination with column indexing
1083*
1084            do 30 j = kp1, n
1085               t = a(l,j)
1086               if (l .eq. k) go to 20
1087                  a(l,j) = a(k,j)
1088                  a(k,j) = t
1089   20          continue
1090               call CAXPY(n-k,t,a(k+1,k),1,a(k+1,j),1)
1091   30       continue
1092         go to 50
1093   40    continue
1094            info = k
1095   50    continue
1096   60 continue
1097   70 continue
1098      ipvt(n) = n
1099      if (cabs1(a(n,n)) .eq. 0.0e0) info = n
1100      return
1101      end
1102* ______________________________________________________________________
1103* Title: CSCAL (LINPACK ROUTINE)
1104* Author: jack dongarra
1105* Date: linpack,  3/11/78.
1106* Description: scales a vector by a constant.
1107* ______________________________________________________________________
1108*
1109      subroutine  CSCAL(n,ca,cx,incx)
1110      implicit none
1111*
1112      complex*16 ca,cx(1)
1113      integer incx,n
1114*
1115      integer i,nincx
1116*
1117      if(n.le.0)return
1118      if(incx.eq.1)go to 20
1119*
1120*        code for increment not equal to 1
1121*
1122      nincx = n*incx
1123      do 10 i = 1,nincx,incx
1124        cx(i) = ca*cx(i)
1125   10 continue
1126      return
1127*
1128*        code for increment equal to 1
1129*
1130   20 do 30 i = 1,n
1131        cx(i) = ca*cx(i)
1132   30 continue
1133      return
1134      end
1135* ______________________________________________________________________
1136* Title: ICAMAX (LINPACK ROUTINE)
1137* Author: jack dongarra
1138* Date: linpack, 3/11/78
1139* Description:
1140*     finds the index of element having max. absolute value.
1141* ______________________________________________________________________
1142*
1143      integer function ICAMAX(n,cx,incx)
1144      implicit none
1145*
1146      complex*16 cx(1)
1147      integer incx,n
1148*
1149      real*8 smax
1150      integer i,ix
1151      complex*16 zdum
1152*
1153* statement function
1154      real*8 cabs1
1155      cabs1(zdum) = abs(dble(zdum)) + abs(dimag(zdum))
1156*
1157      ICAMAX = 0
1158      if( n .lt. 1 ) return
1159      ICAMAX = 1
1160      if(n.eq.1)return
1161      if(incx.eq.1)go to 20
1162*
1163*        code for increment not equal to 1
1164*
1165      ix = 1
1166      smax = cabs1(cx(1))
1167      ix = ix + incx
1168      do 10 i = 2,n
1169         if(cabs1(cx(ix)).le.smax) go to 5
1170         ICAMAX = i
1171         smax = cabs1(cx(ix))
1172    5    ix = ix + incx
1173   10 continue
1174      return
1175*
1176*        code for increment equal to 1
1177*
1178   20 smax = cabs1(cx(1))
1179      do 30 i = 2,n
1180         if(cabs1(cx(i)).le.smax) go to 30
1181         ICAMAX = i
1182         smax = cabs1(cx(i))
1183   30 continue
1184      return
1185      end
1186*
1187* ______________________________________________________________________
1188* Title: DETUN
1189* Author: MMJT
1190* Date: 27 Jan 1989
1191* Description: This subroutine detunes the sharp peaks so that
1192* they can be integrated. In effect, it modifies the stacking
1193* probability matrix such that the determinant will never be
1194* exactly zero.
1195*
1196*      ARGUMENTS:
1197*            No input arguments.
1198*
1199*      COMMON VARIABLES:
1200*            uses:    n_layers
1201*
1202*        modifies:    detune
1203* ______________________________________________________________________
1204*
1205      subroutine DETUN()
1206      include 'DIFFaX.par'
1207      include 'DIFFaX.inc'
1208*
1209      integer*4 i, j
1210      real*8 delta
1211*
1212      delta = eps3
1213*
1214* A value of delta = 0.001 is found to be optimum.
1215* If preferred, user can specify 'delta' interactively
1216* using the following three lines.
1217C   30 write(op,400) 'Enter detune parameter'
1218C      read(cntrl,*,err=30,end=999) delta
1219C      if(CFile) write(op,401) delta
1220*
1221      do 10 i = 1, n_layers
1222        do 20 j = 1, n_layers
1223          detune(j,i) = ONE - abs(delta)
1224   20   continue
1225   10 continue
1226*
1227      return
1228999 stop 'ERROR: Bad delta value. DIFFaX aborted.'
1229400 format(1x, a)
1230401 format(1x, g12.5)
1231      end
1232*
1233* ______________________________________________________________________
1234* Title: EQUALB
1235* Author: MMJT
1236* Date: 13 Mar 1990; 21 July 1997; 3 July 2005
1237* Description:  This routine determines if all of the stacking
1238* uncertainty parameters are identical. There are six arrays to be
1239* tested, namely r_B11, r_B22, r_B33, r_B12, r_B23 and r_B31. These are
1240* passed one at a time from OPTIMZ as r_B. The average value of the
1241* r_B parameters is returned in a_B.
1242* The test fails if the user changes the sign, but not the amplitude, of
1243* some of the B11, B22 or B33. EQUALB returns false, and the calculation
1244* is then done the hard way.
1245*
1246*      ARGUMENTS:
1247*            r_B  -  Array of stacking uncertainty parameters. (input).
1248*            av_B  -  Average of r_B. (output).
1249*
1250*      COMMON VARIABLES:
1251*            uses the array 'there'. n_layers
1252* ______________________________________________________________________
1253*
1254      logical function EQUALB(r_B, av_B)
1255      include 'DIFFaX.par'
1256      include 'DIFFaX.inc'
1257*
1258      real*8 r_B(MAX_L,MAX_L)
1259      real*8 av_B
1260*
1261      integer*4 i, j, m
1262      real*8 error
1263*
1264      av_B = ZERO
1265      m = 0
1266      do 10 i = 1, n_layers
1267        do 20 j = 1, n_layers
1268* Examine only those transitions that actually occur
1269          if(there(j,i)) then
1270            m = m + 1
1271            av_B = av_B + r_B(j,i)
1272          endif
1273   20   continue
1274   10 continue
1275*
1276* Take average
1277      if(m.ne.0) av_B = av_B / dble(m)
1278*
1279      error = ZERO
1280* find absolute deviation of coefficients
1281      do 30 i = 1, n_layers
1282        do 40 j = 1, n_layers
1283          if(there(j,i)) error = error + abs(r_B(j,i) - av_B)
1284   40   continue
1285   30 continue
1286*
1287      EQUALB = abs(error).le.abs(eps3*av_B)
1288*
1289      return
1290      end
1291*
1292* ______________________________________________________________________
1293* Title: GETLAY
1294* Author: MMJT
1295* Date: 4 Oct 1989
1296* Description: This function generates a random sequence of layers
1297* weighted by the stacking probabilities. Needed only when the
1298* 'EXPLICIT' and 'RANDOM' options were specified in the 'STACKING'
1299* description.
1300*
1301*      ARGUMENTS:
1302*           No arguments are used. All data is in 'COMMON'.
1303*
1304*      COMMON VARIABLES:
1305*            uses:      l_cnt, l_g, n_layers, l_seq, l_alpha
1306*
1307*        modifies:      l_seq
1308*
1309*      GETLAY returns logical .true. if all went well.
1310* ______________________________________________________________________
1311*
1312      logical function GETLAY()
1313      include 'DIFFaX.par'
1314      include 'DIFFaX.inc'
1315*
1316      logical okay
1317      character*80 messge
1318      integer*4 i, j, idum
1319      real*8 RAN3, x, sum
1320* external function
1321      external RAN3
1322*
1323      GETLAY = .false.
1324      okay = .true.
1325*
1326* initialize random numbers in RAN3
1327      idum = -1
1328*
1329      write(op,100)
1330     |      'Generating a random sequence of ', l_cnt, ' layers.'
1331*
1332* Get the first layer. Even though some stacking transition
1333* probabilities to and from the layer are non-zero, the layer itself
1334* may not actually exist anywhere in the crystal! Must use the
1335* existence probabilities, l_g.
1336      x = RAN3(idum)
1337      if(x.eq.ONE) x = x - eps7
1338      sum = ZERO
1339      i = 1
1340   10 sum = sum + l_g(i)
1341      if(x.gt.sum) then
1342        i = i + 1
1343        if(i.le.n_layers) goto 10
1344* if all goes well, we should not reach here.
1345        messge = 'GETLAY could not generate the first layer.$'
1346        goto 999
1347      endif
1348      l_seq(1) = i
1349*
1350* Now generate the remaining l_cnt-1 layers
1351      j = 2
1352   20 x = RAN3(idum)
1353      if(x.eq.ONE) x = x - eps7
1354      sum = ZERO
1355      i = 1
1356   30 sum = sum + l_alpha(i,l_seq(j-1))
1357      if(x.gt.sum) then
1358        i = i + 1
1359        if(i.le.n_layers) goto 30
1360* if all goes well, we should not reach here.
1361        write(messge,101) 'GETLAY could not generate layer ', j, '.$'
1362        goto 999
1363      endif
1364      l_seq(j) = i
1365      j = j + 1
1366      if(j.le.l_cnt) goto 20
1367*
1368      GETLAY = .true.
1369      return
1370  999 write(op,102) messge(1:index(messge,'$')-1)
1371      return
1372  100 format(1x, a, i4, a)
1373  101 format(a, i4, a)
1374  102 format(1x, 'ERROR: ', a)
1375      end
1376*
1377* ______________________________________________________________________
1378* Title: GETSAD
1379* Author: MMJT
1380* Date: 29 Oct 1989; 16th April 1999
1381* Description: This subroutine generates the selected area diffraction
1382* pattern (SADP). GETSAD returns .true. if all went well. The pattern
1383* is stored in the linear array 'spec'. 'spec' is read by WRTSAD
1384* to write the diffraction pattern image to disk as a binary file.
1385*
1386*      ARGUMENTS:
1387*            FN      -  Function name passed by reference. The
1388*                       choice is between GLQ16 (non-adaptive
1389*                       Gauss-Legendre integration), and AGLQ16
1390*                       (adaptive Gauss-Legendre integration). (input).
1391*            view    -  Choice of beam direction (input).
1392*                              1  =  normal to the plane k = 0
1393*                              2  =  normal to the plane h = 0
1394*                              3  =  normal to the plane h = k
1395*                              4  =  normal to the plane h = -k
1396*            l_upper -  Upper limit of l. (input).
1397*            hk_lim  -  Upper limit of h (or k). (output).
1398*            infile  -  The name of the input data file. (input).
1399*            ok      -  logical flag indicating all went well.
1400*                                                      (output).
1401*
1402*      COMMON VARIABLES:
1403*            uses:      a0, b0, c0, d0, lambda, has_l_mirror, sadblock,
1404*                       spec, loglin, brightness, X_RAY, rad_type
1405*
1406*        modifies:      scaleint
1407* ______________________________________________________________________
1408*
1409      subroutine GETSAD(FN, view, l_upper, hk_lim, infile, ok)
1410      include 'DIFFaX.par'
1411      include 'DIFFaX.inc'
1412*
1413      real*8 FN, l_upper
1414      integer*4 view, hk_lim
1415      character*(*) infile
1416      logical ok
1417*
1418      integer*4 h, k, i, j, n, info_step, info, cnt, LENGTH, origin
1419      real*8 x, S, S_value, ANGLE, W4, PNTINT, theta, Q2, l
1420      real*8 l_lower, dl, high1, high2, intervals
1421      parameter (intervals = TWENTY)
1422*
1423* external functions (FN is either GLQ16 or AGLQ16)
1424      external FN, LENGTH, PNTINT
1425*
1426* external subroutines (Some compilers need them declared external)
1427*      external XYPHSE, PRE_MAT
1428*
1429* statement functions
1430* S is the value of 1/d**2 at hkl
1431      S(h,k,l) = h*h*a0 + k*k*b0 + l*l*c0 + h*k*d0
1432* ANGLE is the Bragg angle (in radians) of the h,k,l plane
1433      ANGLE(h,k,l) = asin(HALF * lambda * sqrt(S(h,k,l)))
1434* W4 is the X-ray polarization factor
1435      W4(theta) = HALF * (ONE + (cos(TWO*theta))**2)
1436*
1437      Q2 = FOUR / (lambda**2)
1438*
1439* check angles are meaningful
1440      S_value = S(0,0,l_upper)
1441      if(S_value.le.ZERO) then
1442        write(op,101)
1443     |     'ERROR: Illegal value in GETSAD: 1/d**2 = ',S_value
1444        ok = .false.
1445        goto 990
1446      endif
1447      if(S_value.gt.Q2) then
1448        write(op,100)
1449        l_upper = TWO / (lambda*sqrt(c0)) - eps10
1450        write(op,101) 'Upper bound reduced to ', l_upper
1451        S_value = S(0,0,l_upper)
1452      endif
1453*
1454* Set h and k limit
1455      if(view.eq.1) then
1456        hk_lim =  int(l_upper * sqrt(c0 / a0))
1457      else if(view.eq.2) then
1458        hk_lim =  int(l_upper * sqrt(c0 / b0))
1459      else if(view.eq.3) then
1460        hk_lim =  int(l_upper * sqrt(c0 / (a0 + b0 + d0)))
1461      else if(view.eq.4) then
1462        hk_lim =  int(l_upper * sqrt(c0 / (a0 + b0 - d0)))
1463      endif
1464*
1465* Get l increment
1466      dl = l_upper / dble(SADSIZE/2)
1467* reset l_upper so that our integration window, of width dl,
1468* straddles l = 0.
1469      l_upper = l_upper - HALF*dl
1470*
1471* Get diffraction pattern.
1472* First check scan limits, and make sure we are not going to overflow
1473* the array 'spec' which will be used to hold the scan information
1474* prior to writing to file in WRTSAD.
1475      if(has_l_mirror) then
1476        l_lower = ZERO - HALF*dl
1477        sadblock = SADSIZE/2
1478        if((hk_lim + 1)*(SADSIZE/2) .gt. MAX_SP) then
1479          hk_lim = MAX_SP/(SADSIZE/2) - 1
1480        endif
1481      else
1482        l_lower = -l_upper
1483        sadblock = SADSIZE
1484        if((hk_lim + 1)*SADSIZE .gt. MAX_SP) then
1485          hk_lim = MAX_SP/SADSIZE - 1
1486        endif
1487      endif
1488      info_step = nint( l_upper / (dl * intervals) )
1489      if(info_step.le.0) info_step = 1
1490*
1491      cnt = 0
1492      do 10 i = 0, hk_lim
1493* set h and k
1494        if(view.eq.1) then
1495          h =  i
1496          k =  0
1497        else if(view.eq.2) then
1498          h =  0
1499          k =  i
1500        else if(view.eq.3) then
1501          h =  i
1502          k =  i
1503        else if(view.eq.4) then
1504          h =  i
1505          k = -i
1506        else
1507          write(op,105) 'ERROR: Illegal view value in GETSAD', view
1508          goto 990
1509        endif
1510*
1511* write progress to screen
1512        if (debug) write(op,102) h, k, infile(1:LENGTH(infile))
1513*
1514        call XYPHSE(h, k)
1515        call PRE_MAT(h, k)
1516*
1517* The following piece of monkey business is here because if there is
1518* no mirror then the l-loop only calculates SADSIZE-1 values. If there
1519* is a mirror then the l-loop makes SADSIZE/2 calculations. We wish
1520* the final binary data file to be in a square array SADSIZE x SADSIZE,
1521* so we pad out the first value with a zero.
1522*
1523        if(.not.has_l_mirror) then
1524          cnt = cnt + 1
1525          spec(cnt) = ZERO
1526        endif
1527*
1528        info = 0
1529        do 10 l = l_lower, l_upper-eps10, dl
1530* If we are at the origin, don't integrate, we'll probably overflow.
1531          if(i.eq.0 .and. abs(l+dl).le.dl+eps10) then
1532            x = ZERO
1533            origin = cnt + 1
1534          else if(S(h,k,l+dl).gt.Q2) then
1535            x = ZERO
1536          else
1537            x = FN(h,k,l,l+dl,ok)
1538            if(.not.ok) goto 999
1539            if(rad_type.eq.X_RAY) x = x*W4(ANGLE(h,k,l+HALF*dl))
1540          endif
1541          cnt = cnt + 1
1542* make sure we do not overflow
1543          if(cnt.gt.MAX_SP) goto 998
1544          spec(cnt) = x
1545          if(mod(info,info_step).eq.0 .and. debug) then
1546            if(loglin.eq.0) then
1547              if(ONE+x.gt.ZERO) then
1548* write out log(1+x), since x may get to be quite small
1549                write(op,103) 'l = ',l,' log(intensity) = ',log(ONE+x)
1550              else
1551                write(op,103) 'l = ', l, ' log(intensity) = ', ZERO
1552              endif
1553            else
1554              write(op,103) 'l = ', l, ' intensity = ', x
1555            endif
1556          endif
1557          info = info + 1
1558   10 continue
1559* check cnt
1560      if(cnt.lt.2) goto 980
1561*
1562* patch the intensity at the origin to put a well-defined peak there
1563      if(has_l_mirror) then
1564* origin = 1, make origin slightly bigger than proceeding value
1565        spec(origin) = (ONE+eps4) * spec(origin+1)
1566      else
1567* make it slightly larger than the biggest adjacent value
1568        spec(origin) = (ONE+eps4) * max(spec(origin-1),spec(origin+1))
1569      endif
1570*
1571* we need to find the second highest peak intensity so as to scale the
1572* data. The highest intensity should be at the origin, which we discard.
1573      high1 = ZERO
1574      high2 = ZERO
1575      do 30 i = 0, hk_lim
1576        do 30 j = 1, sadblock-1
1577          n = i*sadblock + j
1578* check scaling type. If logarithmic, make sure values are always +ve
1579          if(loglin.eq.0) then
1580            if(ONE+spec(n).gt.ZERO) then
1581              spec(n) = log(ONE + spec(n))
1582            else
1583              spec(n) = ZERO
1584            endif
1585          endif
1586* check if origin is the first value. If so, preset high value.
1587          if(n.eq.1 .and. origin.eq.1) then
1588            high1 = spec(origin)
1589            goto 30
1590          endif
1591          x = spec(n)
1592          if(j.eq.1) then
1593            if(x.gt.spec(n+1)) then
1594              if(x.gt.high1) then
1595                high2 = high1
1596                high1 = x
1597              else if(x.gt.high2) then
1598                high2 = x
1599              endif
1600            endif
1601          else
1602            if(x.gt.spec(n-1) .and. x.gt.spec(n+1)) then
1603              if(x.gt.high1) then
1604                high2 = high1
1605                high1 = x
1606              else if(x.gt.high2) then
1607                high2 = x
1608              endif
1609            endif
1610          endif
1611   30 continue
1612*
1613      if(loglin.ne.0 .and. high2.le.ZERO) then
1614        write(op,101)
1615     |  'ERROR in intensity scaling in GETSAD. ''scale factor'' = ',
1616     |                                                      high2
1617        ok = .false.
1618        goto 990
1619      endif
1620*
1621      if(loglin.eq.0 .and. high1.le.ZERO) then
1622        write(op,101)
1623     |  'ERROR in intensity scaling in GETSAD. ''scale factor'' = ',
1624     |                                                      high1
1625        ok = .false.
1626        goto 990
1627      endif
1628*
1629* If logarithmic, scale to the brightest peak
1630* If linear, scale to the 2nd brightest peak
1631* Note: Intensity scaling can be modified by the user-defined
1632* brightness value
1633      if(loglin.eq.0) then
1634        scaleint = brightness * (maxsad - ONE) / high1
1635      else
1636        scaleint = brightness * (maxsad - ONE) / high2
1637      endif
1638*
1639  990 return
1640  980 write(op,105)
1641     |     'Error in GETSAD: loop counter is too small. cnt = ', cnt
1642      ok = .false.
1643      return
1644  998 write(op,104) 'ERROR in GETSAD: spectrum array overflow at h = ',
1645     |                                    h,', k = ',k,', l = ',l
1646      ok = .false.
1647      return
1648  999 write(op,104) 'ERROR in GETSAD at h = ',h,', k = ',k,', l = ',l
1649      return
1650  100 format(1x, 'Upper bound exceeds 180 degrees!')
1651  101 format(1x, a, g12.5)
1652  102 format(1x, 'h = ', i3, ' k = ', i3, 10x, '''', a, '''')
1653  103 format(1x, a, f10.5, a, g12.5)
1654  104 format(1x, 2(a, i3), a, f10.5)
1655  105 format(1x, a, i3)
1656      end
1657*
1658* ______________________________________________________________________
1659* Title: GETSPC
1660* Authors: MWD and MMJT
1661* Date: 17 Mar 1989; Last tweaked on 7 Mar 1995
1662* Description: This subroutine calculates the spectrum.
1663*
1664*      ARGUMENTS:
1665*            FN      -  Function name passed by reference. The
1666*                       choice is between GLQ16 (non-adaptive
1667*                       Gauss-Legendre), and AGLQ16
1668*                       (adaptive Gauss-Legendre). (input).
1669*            infile  -  The name of the input data file. (input).
1670*
1671*      COMMON VARIABLES:
1672*            uses:      CFile, ELECTN, NEUTRN, SymGrpNo, X_RAY, a0
1673*                       any_sharp, b0, bnds_wt, c0, cntrl, d0, d_theta
1674*                       full_brd, lambda, mltplcty, rad_type, rot_only,
1675*                       spec, th2_max, th2_min, theta1, theta2
1676*
1677*        modifies:      full_shrp
1678*
1679*      GETSPC returns logical .true. if all went well.
1680* ______________________________________________________________________
1681*
1682      logical function GETSPC(FN,infile)
1683      include 'DIFFaX.par'
1684      include 'DIFFaX.inc'
1685*
1686      real*8 FN
1687      character*(*) infile
1688*
1689      logical ok, SHARP, on_bndry, l_axis, shrp
1690      integer*4 h, k, h_lower, h_upper, k_lower, k_upper
1691      integer*4 m, i, max_indx
1692      integer*4 LENGTH
1693      real*8 S, Q, theta, tmp, tmp2, tmp3, fact, h_val, k_val
1694      real*8 HKANGL, LL, ANGLE, AGLQ16
1695      real*8 l, hk_th, x, GLQ16, l_max, min_th, max_th
1696      real*8 W1, l1, l0, d_l, INTENS, L_STEP, W2, W3, l00
1697      complex*16 f(MAX_L)
1698*
1699* external functions
1700      external FN, GLQ16, AGLQ16, INTENS, SHARP, L_STEP, LENGTH
1701* external subroutines (Some compilers need them declared external)
1702*      external XYPHSE, PRE_MAT, GET_F, CHWDTH
1703*
1704* statement functions
1705* S is the value of 1/d**2 at hkl
1706      S(h,k,l) = h*h*a0 + k*k*b0 + l*l*c0 + h*k*d0
1707* LL is the maximum allowable l value for a given h,k and theta
1708      LL(theta,h,k) = sqrt((fact * (sin(theta))**2
1709     |                    - h*h*a0 - k*k*b0 - h*k*d0)/ c0)
1710* ANGLE is the Bragg angle (in radians) of the h,k,l plane
1711      ANGLE(h,k,l) = asin(HALF * lambda * sqrt(S(h,k,l)))
1712* HKANGL is the angle between the vector (h_val,k_val,0) and (1,0,0)
1713      HKANGL(k_val,h_val) = atan2(k_val*sqrt(a0*b0 - d0*d0*QUARTER),
1714     |                              (h_val*a0 + k_val*d0*HALF))
1715* These factors are for powder diffraction patterns
1716* W1 is the polarization factor for weighting x-ray intensities
1717* it also incorporates the Lorentz factor
1718      W1(theta) = (ONE + cos(TWO*theta) * cos(TWO*theta)) /
1719     |                  (sin(theta) * sin(TWO*theta))
1720* W2 is the neutron weighting factor--the Lorentz factor
1721      W2(theta) = ONE / (sin(theta) * sin(TWO*theta))
1722* W3 is the electron weighting factor--the Lorentz factor
1723      W3(theta) = ONE / (sin(theta) * sin(TWO*theta))
1724*
1725      GETSPC = .FALSE.
1726      ok = .true.
1727*
1728* Make sure we are within bounds. If not, adjust.
1729      min_th = HALF * th2_min
1730      max_th = HALF * th2_max
1731      max_indx = int((max_th - min_th) / d_theta + 1)
1732*      if(max_indx.gt.MAX_SP) then
1733*        d_theta = (max_th - min_th) / (MAX_SP - 1)
1734*        max_indx = int((max_th - min_th) / d_theta + 1)
1735*        write(op,300) ''
1736*        write(op,250) 'd_theta is too small and has been adjusted to ',
1737*     |                   TWO*d_theta*RAD2DEG
1738*      endif
1739*
1740* 1234 write(op,300)
1741*     | 'Enter 1 for adaptive quadrature over all l values'
1742*      write(op,300) 'on rows with "sharp" spots'
1743*      read(cntrl,*,err=1234) full_shrp
1744*      if(CFile) write(op,400) full_shrp
1745* zero out spectra
1746      do 10 i = 1, MAX_SP
1747        spec(i) = ZERO
1748   10 continue
1749* See if there is a chance of any sharp peaks.
1750* If so, an appropriate step along l is found, and any_sharp is .true.
1751      d_l = L_STEP(ok)
1752      if(.not.ok) goto 999
1753* If no sharp peaks were unambiguously detected, override user.
1754      if(d_l.eq.ZERO) full_shrp = 1
1755* determine extreme values of h
1756      Q = TWO * sin(max_th) / lambda
1757      fact = TWO / lambda
1758      fact = fact * fact
1759* h_upper is the largest value that the index h can have,
1760* consistent with the value of Q. (When cell_gamma is not 90 degrees,
1761* k is not necessarily zero at this extreme).
1762* In case the incredible happens, immune system to rescue.
1763      tmp3 = FOUR * a0 * b0 - d0 * d0
1764      if(tmp3.le.ZERO) goto 990
1765      tmp3 = TWO * Q * sqrt(ONE / tmp3)
1766      h_upper = int(tmp3 * sqrt(b0))
1767      h_lower = -h_upper
1768      k_upper = int(tmp3 * sqrt(a0))
1769      k_lower = -k_upper
1770* scan along h-axis from h_lower to h_upper
1771      do 20 h = h_lower, h_upper
1772* determine limits along k for a given h
1773        do 30 k = k_lower, k_upper
1774* if out of bounds, cycle
1775          if(S(h,k,ZERO).gt.Q*Q) goto 30
1776          l_axis = h.eq.0 .and. k.eq.0
1777          hk_th = theta1
1778          if(.not.l_axis) hk_th = HKANGL(dble(k), dble(h))
1779* see if in wedge to be scanned
1780          if((theta1-hk_th)*(theta2-hk_th).le.eps3 .or.
1781     |            SymGrpNo.eq.1) then
1782* if rotational symmetry only, do not take points on upper wedge plane
1783            if(rot_only .and. (theta2-hk_th).le.eps3 .and.
1784     |            SymGrpNo.ne.1) goto 30
1785            if(SymGrpNo.eq.11 .and. .not.l_axis) goto 30
1786*            write(op,200) 'Integrating along l at ',h,k,
1787*     |            '''',infile(1:LENGTH(infile)),''''
1788            on_bndry = abs(hk_th-theta1).le.eps3 .or.
1789     |                 abs(hk_th-theta2).le.eps3
1790* set up the phases in the structure factors and stacking vectors
1791* which depend on h and k only
1792            call XYPHSE(h, k)
1793            call PRE_MAT(h, k)
1794* assign a corrected shape-broadening half-width
1795            if(finite_width) then
1796              tmp2 = (h + k*cos(PI-cell_gamma))/(Wa*sin(PI-cell_gamma))
1797              tmp3 = k / Wb
1798              ffhkcnst = QUARTER*lambda*sqrt(a0*tmp2*tmp2+b0*tmp3*tmp3)
1799            endif
1800* get starting value of l
1801            if(l_axis) then
1802              tmp = min(d_theta, max_th)
1803              if(tmp.lt.min_th) tmp = min_th
1804              l1 = LL(tmp, h, k)
1805              shrp = any_sharp
1806            else
1807              tmp = ANGLE(h, k, ZERO)
1808              if(tmp.lt.min_th) then
1809                l1 = LL(min_th,h,k)
1810                tmp = ANGLE(h, k, l1)
1811              else
1812                l1 = ZERO
1813              endif
1814              if(any_sharp .and. full_shrp.ne.1) then
1815                shrp = SHARP(h, k, d_l)
1816                if(.not.ok) goto 999
1817              else
1818                shrp = any_sharp
1819              endif
1820            endif
1821* m indexes into the array spec
1822            m = int((tmp - min_th) / d_theta) + 1
1823            if(.not.shrp .or. full_shrp.eq.1) then
1824* broad streak or full adaptive integration over sharp spots
1825*              if(full_shrp.eq.1 .or. full_brd.eq.1)
1826*     |                write(op,300) 'Full adaptive integration'
1827* integrate each d_theta's range of reciprocal space
1828              do 40 theta = tmp, max_th-eps10, d_theta
1829                l0 = l1
1830                tmp2 = min(d_theta, max_th-theta)
1831                l1 = LL(theta+tmp2, h, k)
1832* sharp spots; do not use knowledge of where they are
1833                if(shrp) then
1834                  x = AGLQ16(h, k, l0, l1, ok)
1835                else
1836* broad streaks
1837                  x = FN(h, k, l0, l1, ok)
1838                endif
1839                if(.not.ok) goto 110
1840*
1841* include weighting factors for radiation type
1842                if(rad_type.eq.X_RAY) then
1843                  x = TWO * x * W1(theta + HALF * tmp2)
1844                else if(rad_type.eq.NEUTRN) then
1845                  x = TWO * x * W2(theta + HALF * tmp2)
1846                else if(rad_type.eq.ELECTN) then
1847                  x = TWO * x * W3(theta + HALF * tmp2)
1848                else
1849                  ok = .false.
1850                  goto 130
1851                endif
1852*
1853* see if not on l-axis
1854                if(.not.l_axis) then
1855* apply multiplicity factor
1856                  x = x * mltplcty
1857* if on boundary, apply appropriate weighting (mirror vs rotation only)
1858                  if(on_bndry) x = x * bnds_wt
1859                endif
1860                if(finite_width) then
1861                  call CHWDTH(h,k,l0,l1,x,m,max_indx)
1862                else
1863                  spec(m) = spec(m) + x
1864                endif
1865                m = m + 1
1866   40         continue
1867            else
1868* line of sharp spots--detuned delta functions
1869* use knowledge of where spots are
1870* make sure we do all l values a multiple of d_l
1871* starting with spot on hk-plane
1872              write(op,300) 'which is a line of sharp spots'
1873              l00 = ZERO
1874              if(l_axis) then
1875                l00 = d_l
1876   50           if(l00.ge.th2_min) goto 60
1877                  l00 = l00 + d_l
1878                  goto 50
1879   60           continue
1880              endif
1881              l_max = LL(max_th, h, k)
1882* avoid trouble by ignoring l = l_max
1883              do 70 l = l00, l_max, d_l
1884                if(l.eq.l_max) goto 70
1885                theta = ANGLE(h,k,l)
1886                call GET_F(f, S(h,k,l), l)
1887                tmp = INTENS(f, h, k, l, ok) * eps8
1888* find width of peak
1889                x = eps10
1890   80             if(.not.ok) goto 120
1891                  x = TWO * x
1892                  call GET_F(f, S(h,k,l+x), l+x)
1893                  if(INTENS(f, h, k, l+x, ok).gt.tmp .and.
1894     |                  x.le.eps2*d_l) goto 80
1895                if(.not.ok) goto 120
1896                l0 = max(l - x, ZERO)
1897                l1 = min(l + x, l_max)
1898                x = AGLQ16(h, k, l0, l1, ok)
1899                if(.not.ok) goto 110
1900*
1901* include weighting factors for radiation type
1902                if(rad_type.eq.X_RAY) then
1903                  x = TWO * x * W1(theta)
1904                else if(rad_type.eq.NEUTRN) then
1905                  x = TWO * x * W2(theta)
1906                else if(rad_type.eq.ELECTN) then
1907                  x = TWO * x * W3(theta)
1908                else
1909                  ok = .false.
1910                  goto 130
1911                endif
1912*
1913* see if not on l-axis
1914                if(.not.l_axis) then
1915* apply multiplicity factor
1916                  x = x * mltplcty
1917* if on boundary, apply appropriate weighting (mirror vs rotation only)
1918                  if(on_bndry) x = x * bnds_wt
1919                endif
1920                m = int(theta / d_theta) + 1
1921                if(finite_width) then
1922                  call CHWDTH(h,k,l0,l1,x,m,max_indx)
1923                else
1924                  spec(m) = spec(m) + x
1925                endif
1926   70         continue
1927            endif
1928          endif
1929   30   continue
1930   20 continue
1931      GETSPC = .true.
1932      return
1933  110 write(op,300) 'GLQ16 returned error in GETSPC.'
1934      return
1935  120 write(op,300) 'INTENS returned error in GETSPC.'
1936      return
1937  130 write(op,300) 'ERROR: Radiation type is undefined in GETSPC'
1938  999 return
1939  990 write(op,300) 'Illegal cell parameters in GETSPC.'
1940      write(op,250) '4*a0*b0-d0*d0 = ', FOUR * a0 * b0 - d0 * d0
1941      return
1942  200 format(1x, a, 2i4, 6x, 3a)
1943  250 format(1x, a, g12.5)
1944  300 format(1x, a)
1945  400 format(1x, i3)
1946      end
1947*
1948* ______________________________________________________________________
1949* Title: GET_BDS
1950* Author: MMJT
1951* Date: 1 July 1989
1952* Description:  This routine assigns reciprocal space vectors in the
1953* h,k plane within which integration can be confined. Weightings are
1954* assigned to off-axis spot intensities to allow for their
1955* multiplicity relative to spots on the 00l axis. Spots that occur on
1956* one of the boundaries are assigned special weighting depending on
1957* whether or not the boundary is also a mirror plane.
1958*
1959*      ARGUMENTS:
1960*           No arguments are used. All data is in 'COMMON'.
1961*
1962*      COMMON VARIABLES:
1963*            uses:      SymGrpNo, h_mirror, k_mirror
1964*
1965*        modifies:      h_start, k_start, h_end, k_end, mltplcty,
1966*                       bnds_wt, rot_only, pnt_grp
1967* ______________________________________________________________________
1968*
1969      subroutine GET_BDS()
1970      include 'DIFFaX.par'
1971      include 'DIFFaX.inc'
1972*
1973* 360 degrees, no symmetry (-1)
1974* note, the scan vectors are not used in this instance since there is
1975* an ambiguity between 0 and 360 degrees. We assign values anyway, so
1976* that the variables are at least initialized in a controlled way.
1977      if(SymGrpNo.eq.1) then
1978        h_start  =  ONE
1979        k_start  =  ZERO
1980        h_end    =  ONE
1981        k_end    =  ZERO
1982        mltplcty =  ONE
1983        bnds_wt  =  ONE
1984        rot_only = .true.
1985      endif
1986*
1987* 180 degrees, rotation only (2/M, 1st setting)
1988      if(SymGrpNo.eq.2) then
1989        h_start  =  ONE
1990        k_start  =  ZERO
1991        h_end    = -ONE
1992        k_end    =  ZERO
1993        mltplcty =  TWO
1994        bnds_wt  =  ONE
1995        rot_only = .true.
1996      endif
1997*
1998* 180 degrees, vertical mirror (2/M, 2nd setting)
1999* we need to know which mirror plane.
2000      if(SymGrpNo.eq.3) then
2001        if(h_mirror .and. .not.k_mirror) then
2002          h_start  =  ONE
2003          k_start  =  ZERO
2004          h_end    = -ONE
2005          k_end    =  ZERO
2006          mltplcty =  TWO
2007          bnds_wt  =  HALF
2008          rot_only = .false.
2009        else if(k_mirror .and. .not.h_mirror) then
2010          h_start  =  ZERO
2011          k_start  =  ONE
2012          h_end    =  ZERO
2013          k_end    = -ONE
2014          mltplcty =  TWO
2015          bnds_wt  =  HALF
2016          rot_only = .false.
2017        else
2018* In theory, the following should never be needed, but just in case,
2019* let's bolster DIFFaX's immune system.
2020          write(op,400) 'DIFFaX is confused about vertical mirrors.'
2021          write(op,400) 'To be safe, symmetry is being set to -1'
2022          SymGrpNo = 1
2023          pnt_grp = '-1'
2024          h_start  =  ONE
2025          k_start  =  ZERO
2026          h_end    =  ONE
2027          k_end    =  ZERO
2028          mltplcty =  ONE
2029          bnds_wt  =  ONE
2030          rot_only = .true.
2031        endif
2032      endif
2033*
2034* 90 degrees, vertical mirrors (MMM)
2035      if(SymGrpNo.eq.4) then
2036        h_start  =  ONE
2037        k_start  =  ZERO
2038        h_end    =  ZERO
2039        k_end    =  ONE
2040        mltplcty =  FOUR
2041        bnds_wt  =  HALF
2042        rot_only = .false.
2043      endif
2044*
2045* 120 degrees, rotation only (-3)
2046      if(SymGrpNo.eq.5) then
2047        h_start  =  ONE
2048        k_start  =  ZERO
2049        h_end    = -ONE
2050        k_end    =  ONE
2051        mltplcty =  THREE
2052        bnds_wt  =  ONE
2053        rot_only = .true.
2054      endif
2055*
2056* 60 degrees, vertical mirrors (-3M)
2057      if(SymGrpNo.eq.6) then
2058        h_start  =  ONE
2059        k_start  =  ZERO
2060        h_end    =  ZERO
2061        k_end    =  ONE
2062        mltplcty =  SIX
2063        bnds_wt  =  HALF
2064        rot_only = .false.
2065      endif
2066*
2067* 90 degrees, rotation (4/M)
2068      if(SymGrpNo.eq.7) then
2069        h_start  =  ONE
2070        k_start  =  ZERO
2071        h_end    =  ZERO
2072        k_end    =  ONE
2073        mltplcty =  FOUR
2074        bnds_wt  =  ONE
2075        rot_only = .true.
2076      endif
2077*
2078* 45 degrees, vertical mirrors (4/MMM)
2079      if(SymGrpNo.eq.8) then
2080        h_start  =  ONE
2081        k_start  =  ZERO
2082        h_end    =  ONE
2083        k_end    =  ONE
2084        mltplcty =  EIGHT
2085        bnds_wt  =  HALF
2086        rot_only = .false.
2087      endif
2088*
2089* 60 degrees, rotation (6/M)
2090      if(SymGrpNo.eq.9) then
2091        h_start  =  ONE
2092        k_start  =  ZERO
2093        h_end    =  ZERO
2094        k_end    =  ONE
2095        mltplcty =  SIX
2096        bnds_wt  =  ONE
2097        rot_only = .true.
2098      endif
2099*
2100* 30 degrees, vertical mirrors (6/MMM)
2101      if(SymGrpNo.eq.10) then
2102        h_start  =  ONE
2103        k_start  =  ZERO
2104        h_end    =  ONE
2105        k_end    =  ONE
2106        mltplcty =  TWELVE
2107        bnds_wt  =  HALF
2108        rot_only = .false.
2109      endif
2110*
2111* integrate along 0 0 l only (axial)
2112* the following are somewhat arbitrary in this case. Assign values
2113* anyway just to make sure they are initialized
2114      if(SymGrpNo.eq.11) then
2115        h_start  =  ONE
2116        k_start  =  ZERO
2117        h_end    =  ONE
2118        k_end    =  ZERO
2119        mltplcty =  ONE
2120        bnds_wt  =  ONE
2121        rot_only = .true.
2122      endif
2123*
2124      return
2125  400 format(1x, a)
2126      end
2127*
2128* ______________________________________________________________________
2129* Title: GET_F
2130* Author: MWD and MMJT
2131* Date: 22 Mar 1989
2132* Description: This routine calculates the form factors for each layer.
2133* Since this routine is the main holdup for complex structures, it
2134* attempts to make use of shortcuts detected in OPTIMZ. The Debye-
2135* Waller diffuse background is assumed to be small and is not
2136* calculated in this version.
2137*
2138*      ARGUMENTS:
2139*            f   -  Array of layer form factors. (output).
2140*            Q2  -  Value of 1/d**2 at h,k,l. (input).
2141*            l   -  reciprocal lattice vector l-component. (input).
2142*
2143*      COMMON VARIABLES:
2144*            uses:  x_sf, n_sf, e_sf, rad_type, n_atoms, l_symmetry,
2145*                   one_B, l_n_atoms, a_type, hx_ky, a_pos, a_occup,
2146*                   a_B, l_actual, CENTRO, ELECTN, NEUTRN, X_RAY
2147*                   n_actual, n_layers
2148*
2149*        modifies:  no COMMON variables are modified
2150* ______________________________________________________________________
2151*
2152      subroutine GET_F(f, S2, l)
2153      include 'DIFFaX.par'
2154      include 'DIFFaX.inc'
2155*
2156      real*8 S2, l
2157      complex*16 f(MAX_L)
2158*
2159      integer*4 i, j, m, n, type
2160      real*8 fact(MAX_TA), tmp(MAX_TA), tmp_sum, dot, e_factor, Q2
2161      parameter(e_factor = 0.023934D0)
2162      complex*16 ctmp(MAX_TA), f_uniq(MAX_L), ctmp_sum
2163*
2164      Q2 = QUARTER * S2
2165* Q2 = sin(theta)**2 / lamba**2
2166* determine scattering factors for each atom type
2167      if(rad_type.eq.X_RAY .or. rad_type.eq.ELECTN) then
2168        do 10 i = 1, n_atoms
2169* This empirical formula comes from p. 71 of
2170* "International Tables for X-ray Crystallography, Vol. IV"
2171* (The Kynoch Press: Birmingham, England), 1974.
2172          fact(i) = x_sf(1,i) * exp(-x_sf(2,i) * Q2) +
2173     |              x_sf(3,i) * exp(-x_sf(4,i) * Q2) +
2174     |              x_sf(5,i) * exp(-x_sf(6,i) * Q2) +
2175     |              x_sf(7,i) * exp(-x_sf(8,i) * Q2) +
2176     |              x_sf(9,i)
2177   10 continue
2178      else if(rad_type.eq.NEUTRN) then
2179        do 20 i = 1, n_atoms
2180          fact(i) = n_sf(i)
2181   20   continue
2182      endif
2183*
2184* get electron scattering factor from x-ray scattering factor
2185* s = 4 pi sin(theta) / lambda
2186* f_electron(s) = (8 pi**2 m e**2 / h**2) {Z - fx(s)} / s**2
2187*
2188*       = 0.023934 lambda**2 {Z - fx(s)} / sin(theta)**2
2189      if(rad_type.eq.ELECTN) then
2190        do 30 i = 1, n_atoms
2191          fact(i) = e_factor * ( dble(e_sf(i)) - fact(i) ) / Q2
2192   30   continue
2193      endif
2194*
2195      do 80 m = 1, n_actual
2196        tmp_sum  = ZERO
2197        ctmp_sum = C_ZERO
2198        do 35 n = 1, n_atoms
2199          tmp(n)  = ZERO
2200          ctmp(n) = C_ZERO
2201   35   continue
2202*
2203* First calculate the scattering factors of the unique layers.
2204* Check to see if f_uniq(m) will all be real and if Debye-Waller is
2205* invariant
2206* Note: hx_ky(j,m) contains h*a_pos(1,j,m) + k*a_pos(2,j,m)
2207*
2208        if(l_symmetry(m).eq.CENTRO .and. one_B(m)) then
2209          do 40 j = 1, l_n_atoms(m)
2210            type = a_type(j,m)
2211            dot = hx_ky(j,m) + l*a_pos(3,j,m)
2212            tmp(type) = tmp(type) + a_occup(j,m) * cos(dot)
2213   40     continue
2214          do 45 j = 1, n_atoms
2215            tmp_sum = tmp_sum + tmp(j) * fact(j)
2216   45     continue
2217* NOTE: twice since centrosymmetric
2218          f_uniq(m) = TWO * exp(-a_B(1,m) * Q2) * tmp_sum
2219*
2220* Debye-Waller is not invariant
2221        else if(l_symmetry(m).eq.CENTRO) then
2222          do 50 j = 1, l_n_atoms(m)
2223            type = a_type(j,m)
2224            dot = hx_ky(j,m) + l*a_pos(3,j,m)
2225            tmp(type) = tmp(type) + a_occup(j,m) *
2226     |              exp(-a_B(j,m) * Q2) * cos(dot)
2227   50     continue
2228          do 55 j = 1, n_atoms
2229            tmp_sum = tmp_sum + tmp(j) * fact(j)
2230   55     continue
2231* NOTE: twice since centrosymmetric
2232          f_uniq(m) = TWO * tmp_sum
2233*
2234* check if Debye-Waller is the only invariant
2235* f(i) will be complex
2236        else if(one_B(m)) then
2237          do 60 j = 1, l_n_atoms(m)
2238            type = a_type(j,m)
2239            dot = hx_ky(j,m) + l*a_pos(3,j,m)
2240            ctmp(type) = ctmp(type) + a_occup(j,m) *
2241     |               dcmplx(cos(dot), sin(dot))
2242   60     continue
2243          do 65 j = 1, n_atoms
2244            ctmp_sum = ctmp_sum + ctmp(j) * fact(j)
2245   65     continue
2246          f_uniq(m) = exp(-a_B(1,m) * Q2) * ctmp_sum
2247*
2248* Nothing is invariant
2249        else
2250          do 70 j = 1, l_n_atoms(m)
2251            type = a_type(j,m)
2252            dot = hx_ky(j,m) + l*a_pos(3,j,m)
2253            ctmp(type) = ctmp(type) + a_occup(j,m) *
2254     |             exp(-a_B(j,m) * Q2) * dcmplx(cos(dot), sin(dot))
2255   70     continue
2256          do 75 j = 1, n_atoms
2257            ctmp_sum = ctmp_sum + ctmp(j) * fact(j)
2258   75     continue
2259          f_uniq(m) = ctmp_sum
2260        endif
2261   80 continue
2262*
2263* Now assign scattering factors to all the layers
2264      do 90 i = 1, n_layers
2265        f(i) = f_uniq(l_actual(i))
2266   90 continue
2267*
2268      return
2269      end
2270*
2271* ______________________________________________________________________
2272* Title: GET_G
2273* Author: MWD and MMJT
2274* Date: 18 Aug 1988; 15 Mar 1995
2275* Description: This function determines g_i, the a-priori probability
2276* that a layer of type i, will actually occur somewhere within the
2277* crystal.
2278* 'cnt' counts the l_alpha(i,i) = 1.0 terms. Then, l_g(i) = 1.0/cnt.
2279*
2280*      ARGUMENTS:
2281*            No input arguments.
2282*
2283*      COMMON VARIABLES:
2284*            uses:   n_layers, l_alpha
2285*
2286*        modifies:   l_g
2287*
2288*      GET_G returns logical .true. if all went well.
2289* ______________________________________________________________________
2290*
2291      logical function GET_G()
2292      include 'DIFFaX.par'
2293      include 'DIFFaX.inc'
2294*
2295      logical singular, LUDCMP
2296      integer*4 i, j, cnt, index(MAX_L)
2297      real*8 sum, g_mat(MAX_L,MAX_L), Det
2298*
2299* external function
2300      external LUDCMP
2301* external subroutine (Some compilers need them declared external)
2302*      external LUBKSB
2303*
2304      GET_G = .false.
2305* set up matrix that represents the probabilities
2306* only n-1 equations are independent
2307      do 10 i = 1, n_layers - 1
2308        l_g(i) = ZERO
2309        sum = ZERO
2310        do 20 j = 1, n_layers
2311          sum = sum + l_alpha(j,i)
2312   20   continue
2313        sum = ONE / sum
2314* sum should actually be ONE
2315        do 30 j = 1, n_layers
2316          g_mat(i,j) = sum * l_alpha(i,j)
2317   30   continue
2318        g_mat(i,i) = g_mat(i,i) - ONE
2319   10 continue
2320      l_g(n_layers) = ONE
2321*
2322* the sum of the g's must be 1
2323      do 40 i = 1, n_layers
2324        g_mat(n_layers, i) = ONE
2325   40 continue
2326*
2327* before we invert the matrix, let's catch the pathological values
2328      cnt = 0
2329      do 50 i = 1, n_layers
2330        if(l_alpha(i,i).eq.ONE) cnt = cnt + 1
2331   50 continue
2332*
2333      if(cnt.ne.0) then
2334        singular = .true.
2335      else
2336        singular = .false.
2337      endif
2338*
2339      if(singular) then
2340        do 60 i = 1, n_layers
2341          if(l_alpha(i,i).eq.ONE) then
2342* arbitrarily assume such layers occur with equal probability
2343            l_g(i) = ONE / dble(cnt)
2344          else
2345            l_g(i) = ZERO
2346          endif
2347   60   continue
2348      else
2349* solve the matrix
2350        if(.not. LUDCMP(g_mat,index,n_layers,MAX_L,Det)) goto 100
2351        call LUBKSB(g_mat,l_g,index,n_layers,MAX_L)
2352      endif
2353*
2354* If we are here then all went well
2355      GET_G = .true.
2356*
2357* Are some of the layers non-existent?
2358      do 70 i = 1, n_layers
2359        if(l_g(i).lt.eps6) then
2360            write(op,410) 'WARNING: Layer ', i,
2361     |      ' does not occur in any significant quantity.'
2362          endif
2363   70 continue
2364*
2365      return
2366  100 write(op,400)
2367     | 'ERROR: Stacking probabilities give a singular matrix in GET_G'
2368      return
2369  400 format(1x, a)
2370  410 format(1x, a, i2, a)
2371      end
2372*
2373* ______________________________________________________________________
2374* Title: GET_MAT
2375* Author: MMJT
2376* Date: 21 Mar 1990; 14th July 1995; 21st July 1997
2377* Description:  This subroutine calculates the elements of 'mat', the
2378* stacking transition matrix ie. {alpha * exp(2*pi*u.R)}ij. The h and k
2379* components were pre-calculated in PRE_MAT. GET_MAT calculates the
2380* l-component to complete the matrix. However, the mat(i,i) terms
2381* have yet to have 1 subtracted from them. This is done in GET_S.
2382*
2383*      ARGUMENTS:
2384*            h   -  reciprocal vector h-component. (input).
2385*            k   -  reciprocal vector k-component. (input).
2386*            l   -  reciprocal lattice vector l-component. (input).
2387*
2388*      COMMON VARIABLES:
2389*            uses:  same_rz, n_layers, same_Bs, Bs_zero, mat1, c0, bc0,
2390*                   ca0, r_B33, r_B23, r_B31, l_r, PI2, there,
2391*                   fatsWalla_hk
2392*
2393*        modifies:  mat
2394* ______________________________________________________________________
2395*
2396      subroutine GET_MAT(h, k, l)
2397      include 'DIFFaX.par'
2398      include 'DIFFaX.inc'
2399*
2400      integer*4 h, k
2401      real*8 l
2402*
2403      real*8 dot, twopi_l, fatsWaller
2404      integer*4 i, j
2405*
2406* set up matrix that represents the sequences
2407* Note: mat is in 'i,j' order.
2408      twopi_l = PI2 * l
2409      if(same_Bs) then
2410        if(all_Bs_zero) then
2411* all the fatsWalla_hk terms equal 1
2412          do 10 i = 1, n_layers
2413            do 20 j = 1, n_layers
2414              if(there(j,i)) then
2415                dot = twopi_l * l_r(3,j,i)
2416                mat(i,j) = mat1(i,j) * dcmplx( cos(dot), sin(dot) )
2417              else
2418                mat(i,j) = C_ZERO
2419              endif
2420   20       continue
2421   10     continue
2422        else
2423* all the fatsWalla terms are identical, but are less than 1
2424* fatsWalla_hk already contains the h-k component computed in PRE_MAT
2425          fatsWaller = fatsWalla_hk * exp(-(l*(QUARTER*a_B33*c0*l
2426     |         + HALF*(a_B23*bc0*k + a_B31*ca0*h))))
2427          do 30 i = 1, n_layers
2428            do 40 j = 1, n_layers
2429              if(there(j,i)) then
2430                dot = twopi_l * l_r(3,j,i)
2431                mat(i,j) = mat1(i,j) * fatsWaller
2432     |                   * dcmplx( cos(dot), sin(dot) )
2433              else
2434                mat(i,j) = C_ZERO
2435              endif
2436   40       continue
2437   30     continue
2438        endif
2439      else
2440* the fatsWalla terms differ. mat1 already contains the h-k component
2441        do 50 i = 1, n_layers
2442          do 60 j = 1, n_layers
2443            if(there(j,i)) then
2444              dot = twopi_l * l_r(3,j,i)
2445              if(Bs_zero(j,i)) then
2446                mat(i,j) = mat1(i,j) * dcmplx( cos(dot), sin(dot) )
2447              else
2448                mat(i,j) = mat1(i,j) * dcmplx( cos(dot), sin(dot) )
2449     |            * exp( -(l*(QUARTER*r_B33(j,i)*c0*l
2450     |            + HALF*(r_B23(j,i)*bc0*k + r_B31(j,i)*ca0*h))) )
2451              endif
2452            else
2453              mat(i,j) = C_ZERO
2454            endif
2455   60     continue
2456   50   continue
2457      endif
2458*
2459      return
2460      end
2461*
2462* ______________________________________________________________________
2463* Title: GET_S
2464* Author: MWD and MMJT
2465* Date: 5th Aug 1991
2466* Description:  This function determines the S's--the average scattered
2467* wave functions of each layer at h, k, l for crystals with an
2468* infinite number of layers. GET_MAT should be called prior to GET_S.
2469* Note, since GET_S is called billions of times from deep within the
2470* inner loops of DIFFaX's bowels, part of the matrix mat1 has been
2471* precalculated in PRE_MAT in order to speed things up a bit.
2472*
2473*      ARGUMENTS:
2474*            f   -  Array of layer form factors. (input).
2475*            s   -  Array of average layer wavefunctions. (output).
2476*            h   -  reciprocal vector h-component. (input).
2477*            k   -  reciprocal vector k-component. (input).
2478*            l   -  reciprocal lattice vector l-component. (input).
2479*
2480*      COMMON VARIABLES:
2481*            uses:  n_layers
2482*
2483*        modifies:  mat
2484*
2485*      GET_S returns logical .true. if all went well.
2486* ______________________________________________________________________
2487*
2488      logical function GET_S(f, s, h, k, l)
2489      include 'DIFFaX.par'
2490      include 'DIFFaX.inc'
2491*
2492      integer*4 h, k
2493      real*8 l
2494      complex*16 f(MAX_L), s(MAX_L)
2495*
2496* i_ok is used by Linpack routines
2497      integer i_ok, index(MAX_L)
2498      integer*4 i
2499      complex*16 Det, s_tmp(2)
2500* external subroutines (Some compilers need them declared external)
2501* CGEFA and CGESL are Linpack routines
2502      external CGEFA, CGESL
2503*
2504      GET_S = .false.
2505*
2506* subtract identity matrix (need do this for diagonal terms only).
2507      do 10 i = 1, n_layers
2508        mat(i,i) = mat(i,i) - C_ONE
2509   10 continue
2510*
2511* Now solve the system of equations.
2512      if(n_layers.gt.2) then
2513* now call LINPACK routines
2514        call CGEFA(mat, MAX_L, n_layers, index, i_ok)
2515        if(i_ok.ne.0) goto 999
2516        call CGESL(mat, MAX_L, n_layers, index, s, 0)
2517      else if(n_layers.eq.2) then
2518* its a simple 2 x 2, so solve it directly
2519        Det = mat(1,1) * mat(2,2) - mat(1,2) * mat(2,1)
2520        if(Det.eq.C_ZERO) goto 999
2521* copy s (remember, if infinitely thick, s = -f)
2522        s_tmp(1) = -s(1)
2523        s_tmp(2) = -s(2)
2524        s(1) = ( mat(1,2) * s_tmp(2) - mat(2,2) * s_tmp(1) ) / Det
2525        s(2) = ( mat(2,1) * s_tmp(1) - mat(1,1) * s_tmp(2) ) / Det
2526      else if(n_layers.eq.1) then
2527* only one layer, so solve it immediately
2528        s(1) = -f(1) / mat(1,1)
2529      endif
2530*
2531      GET_S = .true.
2532      return
2533  999 write(op,400) 'Solving for sequence produces a singular matrix.'
2534      write(op,401) h, k, l
2535      do 50 i = 1, n_layers
2536        write(op,402) i, f(i)
2537   50 continue
2538      return
2539  400 format(1x, 'GET_S:', a)
2540  401 format(1x, 'h = ', i3, ' k = ', i3, ' l = ', g12.5)
2541  402 format(5x, 'f(', i2, ') = (', g12.5, ',', g12.5, ')')
2542      end
2543*
2544* ______________________________________________________________________
2545* Title: GET_S2
2546* Author: MMJT
2547* Date: 5 Feb 1990
2548* Description:  This function determines the S's--the average scattered
2549* wave functions of each layer at h, k, l for crystals with only a
2550* finite number of layers, l_cnt. The equation being solved is
2551*
2552*   inv(Ident-T) * ((N+1)*Ident - inv(Ident-T)*(Ident-T**(N+1)) * F / N
2553*
2554*  where N = l_cnt, and T is the stacking probability matrix
2555*
2556*      ARGUMENTS:
2557*            f   -  Array of layer form factors. (input).
2558*            s   -  Array of average layer wavefunctions. (output).
2559*            h   -  reciprocal vector h-component. (input).
2560*            k   -  reciprocal vector k-component. (input).
2561*            l   -  reciprocal lattice vector l-component. (input).
2562*
2563*      COMMON VARIABLES:
2564*            uses:  n_layers, l_cnt
2565*
2566*        modifies:  mat
2567*
2568*      GET_S2 returns logical .true. if all went well.
2569* ______________________________________________________________________
2570*
2571      logical function GET_S2(f, s, h, k, l)
2572      include 'DIFFaX.par'
2573      include 'DIFFaX.inc'
2574*
2575      integer*4 h, k
2576      real*8 l
2577      complex*16 f(MAX_L), s(MAX_L)
2578*
2579      logical ok, GET_S, MAT2N
2580      integer*4 i, j
2581      complex*16 ctmp, mat_n(MAX_L,MAX_L), tmp_mat(MAX_L,MAX_L)
2582* external functions
2583      external GET_S, MAT2N
2584*
2585      GET_S2 = .false.
2586*
2587* get matrix mat multiplied by itself l_cnt+1 times
2588      ok = MAT2N(mat_n)
2589      if(.not.ok) goto 990
2590*
2591* subtract identity matrix, and make a copy of mat.
2592      do 10 i = 1, n_layers
2593        do 20 j = 1, n_layers
2594          tmp_mat(j,i) = mat(j,i)
2595   20   continue
2596        mat_n(i,i) = mat_n(i,i) - C_ONE
2597   10 continue
2598*
2599* Multiply out -(Ident - T**(l_cnt+1))F.
2600      do 30 i = 1, n_layers
2601        ctmp = C_ZERO
2602        do 40 j = 1, n_layers
2603          ctmp = ctmp + mat_n(i,j) * f(j)
2604   40   continue
2605        s(i) = ctmp
2606   30 continue
2607* Next, solve. ie. inv(Ident - T) * (Ident - T**(l_cnt+1))*F
2608* where now s = (Ident - T**(l_cnt+1))*F
2609      ok = GET_S(f, s, h, k, l)
2610      if(.not.ok) goto 999
2611*
2612* Use result to build a new vector, and solve again.
2613* First, reconstruct mat.
2614      do 50 i = 1, n_layers
2615        s(i) = (s(i) - f(i) * dble(l_cnt + 1)) / dble(l_cnt)
2616        do 60 j = 1, n_layers
2617          mat(j,i) = tmp_mat(j,i)
2618   60   continue
2619   50 continue
2620* Solve with new RHS vector
2621      ok = GET_S(f, s, h, k, l)
2622      if(.not.ok) goto 999
2623*
2624      GET_S2 = .true.
2625      return
2626  990 write(op,400) 'MAT2N returned an error in GET_S2.'
2627      write(op,401) h, k, l
2628      return
2629  999 write(op,400) 'Solving for sequence produces a singular matrix.'
2630      write(op,401) h, k, l
2631      do 70 i = 1, n_layers
2632        write(op,402) i, f(i)
2633   70 continue
2634      return
2635  400 format(1x, 'GET_S2:', a)
2636  401 format(1x, 'h = ', i3, ' k = ', i3, ' l = ', g12.5)
2637  402 format(5x, 'f(', i2, ') = (', g12.5, ',', g12.5, ')')
2638      end
2639*
2640* ______________________________________________________________________
2641* Title: GET_SYM
2642* Author: MMJT
2643* Date: 19 June 1990
2644* Determines the symmetry of the diffraction pattern, thereby
2645* defining the smallest volume of reciprocal space which needs to
2646* be integrated over. There are only 10 kinematical diffraction point
2647* groups in the presence of streaking. Friedel's law ensures there
2648* is always a center of symmetry, and the possibility of streaking
2649* prohibits the cubic point groups. The point group symmetry number
2650* is returned as GET_SYM.
2651* The 10 groups are:
2652*
2653*              GET_SYM          point group
2654*           -------------------------------
2655*                 1:       |        -1
2656*                 2:       |        2/M(1) (diad along streaks)
2657*                 3:       |        2/M(2) (mirror contains streaks)
2658*                 4:       |        MMM
2659*                 5:       |        -3
2660*                 6:       |        -3M
2661*                 7:       |        4/M
2662*                 8:       |        4/MMM
2663*                 9:       |        6/M
2664*                10:       |        6/MMM
2665*
2666* The point group symbol is returned in the global character string
2667* 'pnt_grp'.
2668*
2669*      ARGUMENTS:
2670*            ok  -  logical flag indicating all went well. (output).
2671*
2672*      COMMON VARIABLES:
2673*            uses:  cell_gamma, DoSymDump, cell_a, cell_b, PI, PI2
2674*                   RAD2DEG
2675*
2676*        modifies:  max_var, pnt_grp, h_mirror, k_mirror
2677*
2678*      GET_SYM returns one of the ten symmetry flags listed above.
2679* ______________________________________________________________________
2680*
2681      integer*4 function GET_SYM(ok)
2682      include 'DIFFaX.par'
2683      include 'DIFFaX.inc'
2684*
2685      logical ok
2686*
2687      logical diad, triad, tetrad, hexad
2688      logical TST_ROT, TST_MIR
2689      logical cell90, cell120, eq_sides
2690      integer*4 idum, rot_sym
2691      real*8 tmp_var
2692*
2693* external functions
2694      external TST_ROT, TST_MIR
2695* external subroutine (Some compilers need them declared external)
2696*
2697* initialize random numbers in RAN3
2698      idum = -1
2699*
2700* initialize function
2701      GET_SYM = 0
2702*
2703      max_var = ZERO
2704      diad   = .false.
2705      triad  = .false.
2706      tetrad = .false.
2707      hexad  = .false.
2708      cell90 =  abs(cell_gamma - HALF*PI) .lt. HALF*PI*eps6
2709      cell120 = abs(cell_gamma - PI2/THREE) .lt. PI2*eps6/THREE
2710*
2711* sample reciprocal space to get an idea of the sort of intensities
2712* that are out there.
2713* test for vertical mirror (equivalent to a 2-fold, Friedel's law)
2714      tmp_var = max_var
2715      rot_sym = 2
2716      diad = TST_ROT(rot_sym, idum, ok)
2717      if(.not.ok) goto 997
2718      if(.not.diad) max_var = tmp_var
2719*
2720* if the cell angle is neither 90 nor 120 degrees, then the symmetry
2721* has to be either -1 or 2/M.
2722      if( .not.cell90 .and. .not.cell120 ) then
2723        if(DoSymDump) then
2724          write(sy,'(a)') ' '
2725          write(sy,220) cell_gamma * RAD2DEG
2726          write(sy,221)
2727        endif
2728*
2729        if(diad) then
2730          GET_SYM = 2
2731          pnt_grp = '2/M(1)'
2732        else
2733          GET_SYM = 1
2734          pnt_grp = '-1'
2735        endif
2736        goto 900
2737      endif
2738      eq_sides = abs(cell_a - cell_b) .le. HALF*eps6*(cell_a + cell_b)
2739      if(.not.eq_sides .and. DoSymDump) then
2740        write(sy,'(a)') ' '
2741        write(sy,225) cell_a, cell_b
2742        write(sy,226)
2743        write(sy,221)
2744      endif
2745*
2746* cell_gamma is either 90 or 120 degrees.
2747* if cell_a = cell_b, higher rotational symmetry is possible
2748      if(eq_sides) then
2749* Note, it is quite possible for an oblique cell (whose cell_gamma is
2750* not equal to 120 degrees) to have 3-fold symmetry. We do not test
2751* for this.
2752        tmp_var = max_var
2753        if(cell120) then
2754          rot_sym = 3
2755          triad = TST_ROT(rot_sym, idum, ok)
2756          if(.not.ok) goto 997
2757        else
2758          triad = .false.
2759        endif
2760        if(.not.triad) max_var = tmp_var
2761        hexad = diad.and.triad
2762        if(hexad.and.DoSymDump) write(sy,200)
2763        if(diad.and.cell90 .and. (.not.triad)) then
2764          tmp_var = max_var
2765          rot_sym = 4
2766          tetrad = TST_ROT(rot_sym, idum, ok)
2767          if(.not.ok) goto 997
2768          if(.not.tetrad) max_var = tmp_var
2769        else
2770          tetrad = .false.
2771        endif
2772      endif
2773*
2774* Now test for mirrors.
2775      tmp_var = max_var
2776      h_mirror = TST_MIR(1, idum, ok)
2777      if(.not.ok) goto 998
2778      if(.not.h_mirror) max_var = tmp_var
2779      tmp_var = max_var
2780      k_mirror = TST_MIR(2, idum, ok)
2781      if(.not.ok) goto 999
2782      if(.not.k_mirror) max_var = tmp_var
2783      tmp_var = max_var
2784      hk_mirror = TST_MIR(3, idum, ok)
2785      if(.not.ok) goto 999
2786      if(.not.hk_mirror) max_var = tmp_var
2787*
2788* If h_mirror does not equal k_mirror, then there cannot be a higher
2789* rotation symmetry than a diad. If, by some bizarre freak, this is
2790* inconsistent with the result of TST_ROT, choose the lower symmetry.
2791      if(h_mirror.neqv.k_mirror) then
2792        if(diad) then
2793          GET_SYM = 2
2794          pnt_grp = '2/M(1)'
2795        else
2796          GET_SYM = 3
2797          pnt_grp = '2/M(2)'
2798        endif
2799        goto 900
2800      endif
2801* Now check for combinations of mirrors and rotation axes.
2802*
2803* 6-fold
2804      if(hexad) then
2805        if(h_mirror.or.hk_mirror) then
2806          GET_SYM = 10
2807          pnt_grp = '6/MMM'
2808        else
2809          GET_SYM = 9
2810          pnt_grp = '6/M'
2811        endif
2812        goto 900
2813      endif
2814* 4-fold
2815      if(tetrad) then
2816        if(h_mirror.or.hk_mirror) then
2817          GET_SYM = 8
2818          pnt_grp = '4/MMM'
2819        else
2820          GET_SYM = 7
2821          pnt_grp = '4/M'
2822        endif
2823        goto 900
2824      endif
2825* 3-fold
2826      if(triad.and.(.not.diad)) then
2827        if(h_mirror.or.hk_mirror) then
2828          GET_SYM = 6
2829          pnt_grp = '-3M'
2830        else
2831          GET_SYM = 5
2832          pnt_grp = '-3'
2833        endif
2834        goto 900
2835      endif
2836* 2-fold
2837      if(diad) then
2838* handle special case of non-orthogonal mesh which has a diad,
2839* no triad, and one mirror. Diad prevails, vertical mirrors ignored.
2840        if((h_mirror.or.hk_mirror).and.cell90) then
2841          GET_SYM = 4
2842          pnt_grp = 'MMM'
2843        else
2844          GET_SYM = 2
2845          pnt_grp = '2/M(1)'
2846        endif
2847        goto 900
2848      endif
2849* if no symmetry has been detected opt for lowest symmetry
2850      GET_SYM = 1
2851      pnt_grp = '-1'
2852*
2853  900 return
2854  997 write(op,228) 'ERROR in GET_SYM: error returned by TST_ROT'
2855      write(op,229) '      while testing for ', rot_sym, '-fold axis.'
2856        return
2857  998 write(op,228) 'ERROR in GET_SYM: error returned by TST_MIR'
2858      write(op,228) '   while testing for mirror about the a - c plane'
2859        return
2860  999 write(op,228) 'ERROR in GET_SYM: error returned by TST_MIR'
2861      write(op,228) '   while testing for mirror about the b - c plane'
2862      return
2863  200 format(1x,'THE 2-FOLD AND 3-FOLD IMPLY 6-FOLD ROTATION SYMMETRY')
2864  220 format(1x, 'Cell_gamma = ', g12.5, ' degrees')
2865  221 format(1x, 'Rotational symmetry higher than 2-fold is unlikely.')
2866  225 format(1x, 'cell-a = ', g12.5, ' cell-b = ', g12.5)
2867  226 format(1x, 'Cell sides are not equal.')
2868  228 format(1x, a)
2869  229 format(1x, a, i2, a)
2870      end
2871*
2872* ______________________________________________________________________
2873* Title: GAUSSN
2874* Author: MMJT
2875* Date: 17 Feb 1990; 7 Mar 1995
2876* Description: This subroutine simulates Gaussian instrumental
2877* broadening. std_dev is in degrees. The algorithm used does not
2878* conserve intensity when std_dev is comparable to d_theta. Intensities
2879* at the extreme ends of the spectrum are corrupted slightly.
2880*
2881*      ARGUMENTS:
2882*            th2_low  -  lowest 2theta angle to consider. (input).
2883*
2884*      COMMON VARIABLES:
2885*            uses:  NONE, PI2, RAD2DEG, blurring, d_theta, FWHM
2886*                   th2_max
2887*
2888*        modifies:  brd_spc, spec
2889*
2890* ______________________________________________________________________
2891*
2892      subroutine GAUSSN(th2_low)
2893      include 'DIFFaX.par'
2894      include 'DIFFaX.inc'
2895*
2896      real*8 th2_low
2897*
2898      integer*4 i, j, n_low, n_high, m
2899      real*8 k1, k2, k3, const, gss, std_dev, tmp, tmp1, tmp2
2900*
2901      if(FWHM.le.ZERO) goto 999
2902      std_dev = FWHM / sqrt(EIGHT * log(TWO))
2903*
2904* check that cut-off is reasonable
2905      if(th2_low.lt.ZERO .or. th2_low.ge.th2_max) then
2906        write(op,101) 'GAUSSN: Cut-off angle ', th2_low,
2907     |        ' is out of bounds. Angle reset to zero.'
2908        th2_low = ZERO
2909      endif
2910*
2911* th2_low is the angle relative to th2_min
2912* 2*d_theta is the angular step size
2913      n_low  = int(HALF*th2_low/d_theta) + 1
2914      n_high = int(HALF*(th2_max-th2_min)/d_theta) + 1
2915*
2916      const = TWO * RAD2DEG * d_theta
2917      k1 = const / (sqrt(PI2) * std_dev)
2918      k2 = HALF * (const / std_dev)**2
2919*
2920      do 10 i = 1, n_high
2921        brd_spc(i) = ZERO
2922   10 continue
2923*
2924* go out to 40 standard deviations, or to the end of the spectrum
2925      m = nint(TWO * TWENTY * std_dev / const)
2926      if(m.gt.n_high) m = n_high
2927      do 20 i = 0, m
2928        k3 = k2*dble(i*i)
2929        gss = k1*exp(-k3)
2930        do 30 j = n_low+1, n_high
2931          tmp1 = ZERO
2932          tmp2 = ZERO
2933          if((j-i).gt.n_low)  tmp1 = spec(j-i)
2934          if((j+i).le.n_high) tmp2 = spec(j+i)
2935          tmp = tmp1 + tmp2
2936          if(i.eq.0) tmp = HALF * tmp
2937          brd_spc(j) = brd_spc(j) + gss * tmp
2938   30   continue
2939   20 continue
2940      return
2941  999 write(op,101) 'Illegal FWHM ', FWHM, ' in GAUSSN.'
2942      write(op,100)'Gaussian instrumental broadening not added'
2943* kill blurring option
2944      blurring = NONE
2945      return
2946  100 format(1x, a)
2947  101 format(1x, a, g12.5, a)
2948      end
2949*
2950* ______________________________________________________________________
2951* Title: GLQ16
2952* Authors: MWD and MMJT
2953* Date: 6 April 1989; 13th July 95
2954*  This routine performs 16-point Gauss-Legendre quadrature on an
2955*  interval in reciprocal space. The interval is (h,k,a) to (h,k,b).
2956*  The routine calls INTENS at each of the 16 points, and if all goes
2957*  well, returns .true..
2958*  In the interests of speed, this routine has the option of calling
2959*  APPR_F, which returns interpolated f values for each of the 16
2960*  points. This modification slows down the procedure for structures
2961*  with few atoms per layer, but speeds it up if there are many atoms
2962*  per layer.
2963*
2964*      ARGUMENTS:
2965*            h   -  reciprocal lattice vector h-component. (input).
2966*            k   -  reciprocal lattice vector k-component. (input).
2967*            a   -  l-value of the lower bound of reciprocal
2968*                   lattice integration region. (input).
2969*            b   -  l-value of the upper bound of reciprocal
2970*                   lattice integration region. (input).
2971*            ok  -  logical flag indicating all went well. (output).
2972*
2973*      COMMON VARIABLES:
2974*            uses:  a0, b0, c0, d0, recrsv
2975*
2976*        modifies:  intp_F
2977*
2978*      GLQ16 returns the integrated value.
2979* ______________________________________________________________________
2980*
2981      real*8 function GLQ16(h, k, a, b, ok)
2982      include 'DIFFaX.par'
2983      include 'DIFFaX.inc'
2984*
2985      logical ok
2986      integer*4 h, k
2987      real*8 a, b
2988*
2989      logical o, too_close
2990      real*8 INTENS, INTEN2
2991      real*8 c1, c2, x1, x2, x3, x4, x5, x6, x7, x8
2992      real*8         w1, w2, w3, w4, w5, w6, w7, w8
2993*
2994      parameter (x1 = 0.095012509837637440185D0)
2995      parameter (x2 = 0.281603550779258913230D0)
2996      parameter (x3 = 0.458016777657227386342D0)
2997      parameter (x4 = 0.617876244402643748447D0)
2998      parameter (x5 = 0.755404408355003033895D0)
2999      parameter (x6 = 0.865631202387831743880D0)
3000      parameter (x7 = 0.944575023073232576078D0)
3001      parameter (x8 = 0.989400934991649932596D0)
3002*
3003      parameter (w1 = 0.189450610455068496285D0)
3004      parameter (w2 = 0.182603415044923588867D0)
3005      parameter (w3 = 0.169156519395002538189D0)
3006      parameter (w4 = 0.149595988816576732081D0)
3007      parameter (w5 = 0.124628971255533872052D0)
3008      parameter (w6 = 0.095158511682492784810D0)
3009      parameter (w7 = 0.062253523938647892863D0)
3010      parameter (w8 = 0.027152459411754094852D0)
3011*
3012      integer*4 i, j
3013* f is approximated by a polynomial of order (n-1)
3014      integer*4 n
3015      parameter (n = 3)
3016      integer*4 list(n)
3017*
3018      real*8 Q2, l
3019      real*8 ag_l(16), samp_l(n)
3020      complex*16 f(MAX_L,16)
3021*
3022* external functions
3023      external INTENS, INTEN2
3024* external subroutines (Some compilers need them declared external)
3025*      external APPR_F, GET_F
3026*
3027* statement function
3028* Q2 is the value of 1/d**2 at hkl
3029      Q2(h,k,l) = h*h*a0 + k*k*b0 + l*l*c0 + h*k*d0
3030*
3031* initialize, to suppress compiler warnings
3032      GLQ16 = ZERO
3033*
3034* check that the integration range is legitimate
3035      if(b.lt.a) goto 999
3036* catch pathological values of a and b - i.e. zero integration range!
3037      if(b.eq.a) then
3038        ok = .true.
3039        goto 900
3040      endif
3041*
3042* let's interpolate ('hard-wired' in this version)
3043      intp_F = .true.
3044*
3045      c1 = HALF*(b - a)
3046      c2 = c1 + a
3047*
3048* set up the 16 l-values to sample
3049      ag_l(1) = -c1*x8 + c2
3050      ag_l(2) = -c1*x7 + c2
3051      ag_l(3) = -c1*x6 + c2
3052      ag_l(4) = -c1*x5 + c2
3053      ag_l(5) = -c1*x4 + c2
3054      ag_l(6) = -c1*x3 + c2
3055      ag_l(7) = -c1*x2 + c2
3056      ag_l(8) = -c1*x1 + c2
3057*
3058      ag_l( 9) = c1*x1 + c2
3059      ag_l(10) = c1*x2 + c2
3060      ag_l(11) = c1*x3 + c2
3061      ag_l(12) = c1*x4 + c2
3062      ag_l(13) = c1*x5 + c2
3063      ag_l(14) = c1*x6 + c2
3064      ag_l(15) = c1*x7 + c2
3065      ag_l(16) = c1*x8 + c2
3066*
3067      if(intp_F) then
3068*
3069* choose special values to sample (3 point interpolation)
3070        list(1) = 1
3071        list(2) = 8
3072        list(3) = 16
3073        samp_l(1) = ag_l(list(1))
3074        samp_l(2) = ag_l(list(2))
3075        samp_l(3) = ag_l(list(3))
3076*
3077* Deal with very rare cases when the spread in l values is too small
3078        too_close = (samp_l(1).eq.samp_l(3)) .or.
3079     |              (samp_l(1).eq.samp_l(2)) .or.
3080     |              (samp_l(2).eq.samp_l(3))
3081*
3082        if(.not.too_close) then
3083          call APPR_F(f, h, k, samp_l, ag_l, n, list, ok)
3084          if(.not.ok) goto 990
3085        else
3086* sample f values once, and set all 16 values over l-range to be equal
3087          call GET_F(f(1,1), Q2(h,k,ag_l(1)), ag_l(1))
3088          do 10 j = 1, n_layers
3089            do 20 i = 2, 16
3090              f(j,i) = f(j,1)
3091   20       continue
3092   10     continue
3093        endif
3094*
3095      else
3096* do not interpolate
3097        do 30 i = 1, 16
3098          call GET_F(f(1,i), Q2(h,k,ag_l(i)), ag_l(i))
3099   30   continue
3100*
3101      endif
3102*
3103* sum intensities
3104*
3105      o = .true.
3106      if(recrsv) then
3107*
3108      GLQ16 = c1 * (
3109     |w8*(INTENS(f(1,1),h,k,ag_l(1),o)+INTENS(f(1,16),h,k,ag_l(16),o))+
3110     |w7*(INTENS(f(1,2),h,k,ag_l(2),o)+INTENS(f(1,15),h,k,ag_l(15),o))+
3111     |w6*(INTENS(f(1,3),h,k,ag_l(3),o)+INTENS(f(1,14),h,k,ag_l(14),o))+
3112     |w5*(INTENS(f(1,4),h,k,ag_l(4),o)+INTENS(f(1,13),h,k,ag_l(13),o))+
3113     |w4*(INTENS(f(1,5),h,k,ag_l(5),o)+INTENS(f(1,12),h,k,ag_l(12),o))+
3114     |w3*(INTENS(f(1,6),h,k,ag_l(6),o)+INTENS(f(1,11),h,k,ag_l(11),o))+
3115     |w2*(INTENS(f(1,7),h,k,ag_l(7),o)+INTENS(f(1,10),h,k,ag_l(10),o))+
3116     |w1*(INTENS(f(1,8),h,k,ag_l(8),o)+INTENS(f(1, 9),h,k,ag_l( 9),o)))
3117*
3118      else
3119*
3120      GLQ16 = c1 * (
3121     |w8*(INTEN2(f(1,1),h,k,ag_l(1),o)+INTEN2(f(1,16),h,k,ag_l(16),o))+
3122     |w7*(INTEN2(f(1,2),h,k,ag_l(2),o)+INTEN2(f(1,15),h,k,ag_l(15),o))+
3123     |w6*(INTEN2(f(1,3),h,k,ag_l(3),o)+INTEN2(f(1,14),h,k,ag_l(14),o))+
3124     |w5*(INTEN2(f(1,4),h,k,ag_l(4),o)+INTEN2(f(1,13),h,k,ag_l(13),o))+
3125     |w4*(INTEN2(f(1,5),h,k,ag_l(5),o)+INTEN2(f(1,12),h,k,ag_l(12),o))+
3126     |w3*(INTEN2(f(1,6),h,k,ag_l(6),o)+INTEN2(f(1,11),h,k,ag_l(11),o))+
3127     |w2*(INTEN2(f(1,7),h,k,ag_l(7),o)+INTEN2(f(1,10),h,k,ag_l(10),o))+
3128     |w1*(INTEN2(f(1,8),h,k,ag_l(8),o)+INTEN2(f(1, 9),h,k,ag_l( 9),o)))
3129*
3130      endif
3131      ok = o
3132  900 return
3133  999 ok = .false.
3134      write(op,100) 'GLQ16: Illegal integration interval!'
3135      write(op,101) h, k, a, b
3136      return
3137  990 ok = .false.
3138      write(op,100) 'GLQ16: ERROR returned from APPR_F'
3139      write(op,101) h, k, a, b
3140      return
3141  100 format(1x, a)
3142  101 format(1x, 'h = ',i4,', k = ',i4,', l0 = ',g12.5,', l1 = ',g12.5)
3143      end
3144*
3145* ______________________________________________________________________
3146* Title: HKL_LIM
3147* Author: MMJT
3148* Date: 2 August 1989
3149* Obtains upper limits of h, k, and l given the global variable
3150* 'max_angle' (in radians).
3151* The limits are returned in the global variables h_bnd, k_bnd and
3152* l_bnd. HKL_LIM may need to decrease the value of lambda if h_bnd
3153* and k_bnd are too small to allow adequate off-axis symmetry testing.
3154* lambda is restored in OPTIMZ after symmetry testing.
3155*
3156*      ARGUMENTS:
3157*           No arguments are used. All data is in 'COMMON'.
3158*
3159*      COMMON VARIABLES:
3160*            uses:  a0, b0, c0, d0, lambda
3161*
3162*        modifies:  h_bnd, k_bnd, l_bnd, lambda
3163* ______________________________________________________________________
3164*
3165      subroutine HKL_LIM()
3166      include 'DIFFaX.par'
3167      include 'DIFFaX.inc'
3168*
3169      integer*4 HKLUP
3170      real*8 y
3171*
3172* HKLUP returns the maximum value of h, k or l given 'max_angle'
3173      HKLUP(y) =  int(TWO * sin(HALF*max_angle) / (lambda*sqrt(y)))
3174*
3175* define upper h, k, l values consistent with max_angle
3176    1 h_bnd = HKLUP(a0)
3177      k_bnd = HKLUP(b0)
3178      l_bnd = HKLUP(c0)
3179*
3180* make sure bounds are not too small. This could occur for a small
3181* unit cell, a small value of max_angle, or a long wavelength.
3182      if(h_bnd.lt.2 .or. k_bnd.lt.2) then
3183        lambda = HALF * lambda
3184        goto 1
3185      endif
3186*
3187* Make sure bounds are not too large either
3188      if(h_bnd.gt.10) h_bnd = 10
3189      if(k_bnd.gt.10) k_bnd = 10
3190      if(l_bnd.gt.10) l_bnd = 10
3191*
3192      return
3193      end
3194*
3195* ______________________________________________________________________
3196* Title: INTEGR
3197* Author: MMJT and MWD
3198* Date: 15 Feb 1990; 7 Mar 1995; 28th May 1996
3199* Description:  This routine integrates intensity from
3200*               h,k,l0 to h,k,l1.
3201*
3202*      ARGUMENTS:
3203*            FN   -  Function name passed by reference. The
3204*                    choice is between GLQ16 (non-adaptive
3205*                    Gauss-Legendre integration), and AGLQ16
3206*                    (adaptive Gauss-Legendre integration). (input).
3207*            ok   -  logical flag indicating all went well. (output).
3208*
3209*      COMMON VARIABLES:
3210*            uses:   a0, b0, c0, d0, lambda, cntrl, CFile, xplcit,
3211*                    X_RAY, rad_type, th2_max
3212*
3213*        modifies:   no COMMON variables are modified
3214* ______________________________________________________________________
3215*
3216      subroutine INTEGR(FN, ok)
3217      include 'DIFFaX.par'
3218      include 'DIFFaX.inc'
3219*
3220      real*8 FN
3221      logical ok
3222*
3223      logical divided
3224      integer*4 h, k
3225      real*8 l0, l1, max_th, x, W4, ANGLE, theta, l, S, Q2
3226      real*8 t1, sum, tmp, LL, fact, d_th, l_tmp
3227*
3228* external function, passed by reference
3229      external FN
3230* external subroutines (Some compilers need them declared external)
3231*      external XYPHSE, PRE_MAT, GET_F
3232*
3233* statement functions
3234* S is the value of 1/d**2 at hkl
3235      S(h,k,l) = h*h*a0 + k*k*b0 + l*l*c0 + h*k*d0
3236* ANGLE is the Bragg angle (in radians) of the h,k,l plane
3237      ANGLE(h,k,l) = asin(HALF * lambda * sqrt(S(h,k,l)))
3238* LL is the maximum l value for a given h,k
3239      LL(theta,h,k) = sqrt((Q2 * sin(theta) * sin(theta)
3240     |                    - h*h*a0 - k*k*b0 - h*k*d0)/ c0)
3241* W4 is the X-ray polarization factor
3242      W4(theta) = HALF * (ONE + (cos(TWO*theta))**2)
3243*
3244      max_th = HALF * th2_max
3245      fact = TWO / lambda
3246      Q2 = fact * fact
3247* set increment to a safe value, in case l1-l0 is too large
3248      d_th = eps3 * DEG2RAD
3249*
3250   10 write(op,400) 'Enter h, k, l0, l1'
3251      read(cntrl,*,err=10)  h, k, l0, l1
3252      if(CFile) write(op,401) h, k, l0, l1
3253* check values
3254      if(l1.eq.l0) then
3255        write(op,400) 'Illegal input: l1 equals l0'
3256        if(CFile) then
3257          l1 = l0 + d_th
3258          write(op,403) 'l1 is set to ', l1
3259        else
3260          goto 10
3261        endif
3262      endif
3263* make sure we are not going to blow up at the origin
3264      if(h.eq.0 .and. k.eq.0 .and. rad_type.eq.ELECTN) then
3265        if(l0*l1.le.ZERO) then
3266          write(op,400)
3267     |   'Cannot integrate across the origin for electron radiation'
3268          write(op,400) 'Re-enter. . .'
3269          goto 10
3270        endif
3271      endif
3272* Finally, check angles are meaningful
3273      if(S(h,k,l0).gt.Q2 .or. S(h,k,l1).gt.Q2) then
3274        if(S(h,k,l0).gt.Q2) write(op,402) h, k, l0,
3275     |            ' exceeds 180 degree scattering angle!'
3276        if(S(h,k,l1).gt.Q2) write(op,402) h, k, l1,
3277     |            ' exceeds 180 degree scattering angle!'
3278        goto 10
3279      endif
3280* get angles corresponding to start and stop
3281* check if we need to break the integration region into two parts
3282* because h,k,l, and h,k,-l may subtend the same angle
3283      divided  = .false.
3284      if(l0.le.ZERO .and. l1.le.ZERO) then
3285* use Friedel's law, and keep l +ve. Swap l0 and l1
3286        h = -h
3287        k = -k
3288        tmp = -l0
3289        l0 = -l1
3290        l1 = tmp
3291      else if(l0.lt.ZERO .and. l1.gt.ZERO) then
3292        h = -h
3293        k = -k
3294        l_tmp = l1
3295        l1 = -l0
3296        l0 = ZERO
3297        divided = .true.
3298      endif
3299* swap if in reverse order
3300      if(l0.gt.l1) then
3301        tmp = l0
3302        l0 = l1
3303        l1 = tmp
3304      endif
3305*
3306      sum = ZERO
3307   30 max_th = ANGLE(h,k,l1)
3308      t1     = ANGLE(h,k,l0)
3309      l1 = l0
3310      call XYPHSE(h, k)
3311      call PRE_MAT(h, k)
3312* integrate each d_th's range of reciprocal space
3313      do 40 theta = t1, max_th-eps14, d_th
3314        l0 = l1
3315        tmp = min(d_th, max_th-theta)
3316        l1 = LL(theta+tmp,h,k)
3317        x = FN(h,k,l0,l1,ok)
3318        if(.not.ok) goto 999
3319        if(rad_type.eq.X_RAY) x = x * W4(theta + HALF*tmp)
3320        sum = sum + x
3321   40 continue
3322* do we need to integrate the other part?
3323      if(divided) then
3324* goto -h,-k,-l and continue
3325        h = -h
3326        k = -k
3327        l0 = ZERO
3328        l1 = l_tmp
3329        divided = .false.
3330        goto 30
3331      endif
3332      write(op,403) 'Integrated intensity = ', sum
3333  999 return
3334  400 format(1x, a)
3335  401 format(1x, 2i3, 2g12.5)
3336  402 format(1x, 2i3, g12.5, a)
3337  403 format(1x, a, g15.8)
3338      end
3339*
3340* ______________________________________________________________________
3341* Title: INTEN2
3342* Author: MMJT
3343* Date: 8 Apr 1992
3344* Description: This routine determines the intensity at (h,k,l)
3345* in reciprocal space, for an explicitly defined stacking sequence.
3346*
3347*      ARGUMENTS:
3348*            f   -  Layer form factors. (input).
3349*            h   -  reciprocal lattice vector h-component. (input).
3350*            k   -  reciprocal lattice vector k-component. (input).
3351*            l   -  reciprocal lattice vector l-component. (input).
3352*            ok  -  logical flag indicating all went well. (output).
3353*
3354*      COMMON VARIABLES:
3355*            uses:  same_layer, l_seq, l_r, n_layers, l_phi, l_cnt
3356*                   PI2
3357*
3358*        modifies:  wavefn
3359*
3360*      INTEN2 returns the intensity at h, k, l
3361* ______________________________________________________________________
3362*
3363      real*8 function INTEN2(f, h, k, l, ok)
3364      include 'DIFFaX.par'
3365      include 'DIFFaX.inc'
3366*
3367      logical ok
3368      integer*4 h, k
3369      real*8 l
3370      complex*16 f(MAX_L)
3371*
3372      integer*4 i, j, m
3373      real*8 twopi_l, dot, tmp
3374      complex*16 phi(MAX_L, MAX_L), z, z_to_n
3375*
3376      twopi_l = PI2 * l
3377*
3378* 'ok' is not used, but is included for compatibility with INTENS
3379      ok = .true.
3380*
3381* Is there only one layer? If so, let's get it over with.
3382      if(l_cnt.eq.1) then
3383        wavefn = f(l_seq(1))
3384        goto 900
3385      endif
3386*
3387* Check for an obvious optimization when all layers are identical.
3388      if(same_layer) then
3389        i = l_seq(1)
3390        dot = PI2*(h*l_r(1,i,i) + k*l_r(2,i,i) + l*l_r(3,i,i))
3391        z = dcmplx( cos(dot), sin(dot) )
3392        tmp = dot / PI2
3393* check we are not about to execute 0.0 / 0.0
3394        if(abs(tmp - nint(tmp)) .le. eps5) then
3395          wavefn = f(i) * l_cnt
3396        else
3397* sum the series
3398          dot = dot * l_cnt
3399          z_to_n = dcmplx( cos(dot), sin(dot) )
3400          wavefn = f(i) * (C_ONE - z_to_n) / (C_ONE - z)
3401        endif
3402        goto 900
3403      endif
3404*
3405* Else do it the long way
3406* Get phases
3407      do 10 i = 1, n_layers
3408        do 20 j = 1, n_layers
3409          dot = twopi_l * l_r(3,j,i)
3410          phi(j,i) = l_phi(j,i) * dcmplx( cos(dot), sin(dot) )
3411   20   continue
3412   10 continue
3413* Count down to the first layer (we know l_cnt is greater than 1)
3414* Initialize wavefunction to the scattering factor of the last layer
3415      wavefn = f(l_seq(l_cnt))
3416      do 30 m = l_cnt - 1, 1, -1
3417        i = l_seq(m)
3418        j = l_seq(m+1)
3419        wavefn = f(i) + wavefn * phi(j,i)
3420   30 continue
3421*
3422* Normalize to the number of layers
3423  900 INTEN2 = wavefn * conjg(wavefn) / l_cnt
3424*
3425      return
3426      end
3427*
3428* ______________________________________________________________________
3429* Title: INTENS
3430* Author: MWD and MMJT
3431* Date: 10 Feb 1989
3432* Description: This routine determines the intensity at (h,k,l)
3433* in reciprocal space, for recursive stacking. For this function
3434* to be called, 'rcrsv' must be TRUE.
3435* Note: The diffuse background is handled automatically by the
3436* recursion algorithm.
3437*
3438*      ARGUMENTS:
3439*            f   -  Layer form factors. (input).
3440*            h   -  reciprocal lattice vector h-component. (input).
3441*            k   -  reciprocal lattice vector k-component. (input).
3442*            l   -  reciprocal lattice vector l-component. (input).
3443*            ok  -  logical flag indicating all went well. (output).
3444*
3445*      COMMON VARIABLES:
3446*            uses:  only_real, l_g, n_layers, inf_thick
3447*
3448*        modifies:  no COMMON variables are modified
3449*
3450*      INTENS returns the intensity at h, k, l
3451* ______________________________________________________________________
3452*
3453      real*8 function INTENS(f, h, k, l, ok)
3454      include 'DIFFaX.par'
3455      include 'DIFFaX.inc'
3456*
3457      logical ok
3458      integer*4 h, k
3459      real*8 l
3460      complex*16 f(MAX_L)
3461*
3462      logical GET_S, GET_S2
3463      integer*4 i
3464      real*8 sum, x
3465      complex*16 s(MAX_L)
3466*
3467* external function
3468      external GET_S, GET_S2
3469* external subroutines (Some compilers need them declared external)
3470*      external GET_MAT
3471*
3472      sum = ZERO
3473      call GET_MAT(h, k, l)
3474      if(inf_thick) then
3475* initialize s to -f, since mat is -(1 - T)
3476        do 10 i = 1, n_layers
3477          s(i) = - f(i)
3478   10   continue
3479        ok = GET_S(f, s, h, k, l)
3480      else
3481* s is initialized inside GET_S2, from where GET_S is called
3482        ok = GET_S2(f, s, h, k, l)
3483      endif
3484      if(ok) then
3485* only use real part of f(i) if that's all that's there
3486        if(only_real) then
3487          do 20 i = 1, n_layers
3488            sum = sum + l_g(i) * dble(f(i)) * dble(s(i))
3489   20     continue
3490          sum = TWO * sum
3491          do 30 i = 1, n_layers
3492            x = dble(f(i))
3493            sum = sum - l_g(i) * x*x
3494   30     continue
3495        else
3496* must use complex part of f(i)
3497          do 40 i = 1, n_layers
3498            sum = sum + l_g(i) * dble(conjg(f(i)) * s(i))
3499   40     continue
3500          sum = TWO * sum
3501          do 50 i = 1, n_layers
3502            x = abs(f(i))
3503            sum = sum - l_g(i) * x*x
3504   50     continue
3505        endif
3506      endif
3507*
3508      INTENS = sum
3509*
3510      return
3511      end
3512*
3513* ______________________________________________________________________
3514* Title: LENGTH
3515* Author: MMJT
3516* Date: 20 Nov 1989
3517* Description:  This function returns the length of the first group
3518* of contiguous, non-space characters in the passed string. If the
3519* string has no blanks, the declared length of the string is returned.
3520*
3521*      ARGUMENTS:
3522*            string  -  Character string whose length is needed.
3523*                                                            (input).
3524*      LENGTH returns the string length.
3525* ______________________________________________________________________
3526*
3527      integer*4 function LENGTH(string)
3528      implicit none
3529*
3530      character string*(*)
3531*
3532      integer*4 i
3533*
3534      i = index(string,' ')
3535      if(i.eq.0) then
3536        LENGTH = len(string)
3537      else
3538        LENGTH = i-1
3539      endif
3540
3541      return
3542      end
3543*
3544* ______________________________________________________________________
3545* Title: LORNZN
3546* Author: MMJT
3547* Date: 17 Feb 1990; 7 Mar 1995
3548* Description: This subroutine performs the Lorentzian
3549* instrumental broadening. FWHM is in degrees. Does not conserve
3550* intensity well when FWHM is comparable to d_theta. Data at the
3551* extreme ends of the spectrum are corrupted slightly.
3552*
3553*      ARGUMENTS:
3554*            th2_low  -  lowest 2theta angle to consider. (input).
3555*
3556*      COMMON VARIABLES:
3557*            uses:  th2_max, d_theta, FWHM, spec, NONE, PI, RAD2DEG
3558*                   blurring
3559*
3560*        modifies:  brd_spc
3561* ______________________________________________________________________
3562*
3563      subroutine LORNZN(th2_low)
3564      include 'DIFFaX.par'
3565      include 'DIFFaX.inc'
3566*
3567      real*8 th2_low
3568*
3569      integer*4 i, j, n_low, n_high
3570      real*8 k1, k2, k3, const, lrnz, tmp, tmp1, tmp2
3571*
3572      if(FWHM.le.ZERO) goto 999
3573* check that cut-off is reasonable
3574      if(th2_low.lt.ZERO .or. th2_low.ge.th2_max) then
3575        write(op,101) 'LORNZN: Cut-off angle ', th2_low,
3576     |        ' is out of bounds. Angle reset to zero.'
3577        th2_low = ZERO
3578      endif
3579*
3580* th2_low is the angle relative to th2_min
3581* 2*d_theta is the angular step size
3582      n_low  = int(HALF*th2_low/d_theta) + 1
3583      n_high = int(HALF*(th2_max-th2_min)/d_theta) + 1
3584*
3585      const = TWO * RAD2DEG * d_theta
3586      k1 = const * TWO / (PI * FWHM)
3587      k2 = (const * TWO / FWHM)**2
3588*
3589      do 10 i = 1, n_high
3590        brd_spc(i) = ZERO
3591   10 continue
3592*
3593      do 20 i = 0, n_high
3594        k3 = ONE + k2*i*i
3595        lrnz = k1 / k3
3596        do 30 j = n_low+1, n_high
3597          tmp1 = ZERO
3598          tmp2 = ZERO
3599          if((j-i).gt.n_low)  tmp1 = spec(j-i)
3600          if((j+i).le.n_high) tmp2 = spec(j+i)
3601          tmp = tmp1 + tmp2
3602          if(i.eq.0) tmp = HALF * tmp
3603          brd_spc(j) = brd_spc(j) + lrnz * tmp
3604   30   continue
3605   20 continue
3606      return
3607  999 write(op,101) 'Illegal FWHM ', FWHM, ' in LORNZN()'
3608      write(op,100) 'Lorentzian instrumental broadening not added'
3609* kill blurring option
3610      blurring = NONE
3611      return
3612  100 format(1x, a)
3613  101 format(1x, a, g12.5, a)
3614      end
3615*
3616* ______________________________________________________________________
3617* Title: LUBKSB
3618* Author: MWD, adapted from
3619* "Numerical Recipes: The Art of Scientific Computing."
3620* Date: 18 Aug 1988
3621*  Description: Solves the set of linear equations ax = b,
3622*  where a, x and b contain real variables.
3623*  Here, a is input as the LU-decomposition of a, determined by the
3624*  routine LUDCMP. index is input as the permutation vector returned
3625*  by LUDCMP. b is input as the right hand side vector, and returns
3626*  with the solution x. a, n, MAX_N and index are not modified by this
3627*  routine. In DIFFaX, LUDCMP and LUBKSB are used to solve for l_g(i),
3628*  the a-priori probability that layer i exists within the crystal.
3629*
3630*      ARGUMENTS:
3631*            a      -  LU-decomposed square matrix of real numbers.
3632*                                                           (input).
3633*            b      -  vector of real numbers, the right hand side of
3634*                      ax = b, is input. The solution x is output.
3635*            index  -  vector containing the record of the row
3636*                      permutation effected by the partial pivoting
3637*                      carried out in CLUDCM (input).
3638*            n      -  size of the square matrix. (input).
3639*            MAX_N  -  physical dimension of a (MAX_N x MAX_N).(input).
3640* ______________________________________________________________________
3641*
3642      subroutine LUBKSB(a,b,index,n,MAX_N)
3643      include 'DIFFaX.par'
3644*     save
3645*
3646      integer*4 n, MAX_N, index(MAX_N)
3647      real*8 a(MAX_N,MAX_N), b(MAX_N)
3648*
3649      integer*4 i, i2, j, row
3650      real*8 sum
3651*
3652      i2 = 0
3653      do 20 i = 1, n
3654        row = index(i)
3655        sum = b(row)
3656        b(row) = b(i)
3657        if(i2.ne.0) then
3658          do 10 j = i2, i-1
3659            sum = sum - a(i,j) * b(j)
3660   10     continue
3661        else if(abs(sum).ne.ZERO) then
3662          i2 = i
3663        endif
3664        b(i) = sum
3665   20 continue
3666      do 40 i = n, 1, -1
3667        sum = b(i)
3668        do 30 j = i+1, n
3669          sum = sum - a(i,j) * b(j)
3670   30   continue
3671        b(i) = sum / a(i,i)
3672   40 continue
3673      return
3674      end
3675*
3676* ______________________________________________________________________
3677* Title: LUDCMP
3678* Author: MWD, adapted from
3679* "Numerical Recipes: The Art of Scientific Computing."
3680* Date: 18 Aug 1988
3681*  Description: This is an LU decomposition routine, and accepts
3682*  real*8 variables.
3683*  Given an n x n matrix a, with physical dimension MAX_N, this
3684*  routine replaces it by the LU decomposition of a rowwise permutation
3685*  of itself. a and n are input. a is the LU decomposed output; index
3686*  is an output vector which records the row permutation affected by
3687*  the partial pivoting; Det is the determinant of a. This routine is
3688*  used in combination with LUBKSB to solve linear equations. In DIFFaX,
3689*  these routines are used to solve for l_g(i), the a-priori
3690*  probability that layer i exists within the crystal.
3691*  LUDCMP returns .false. if the matrix turns out to be singular.
3692*
3693*      ARGUMENTS:
3694*            a      -  Square matrix of real numbers to LU-decompose is
3695*                      input. a is then replaced by the
3696*                      LU-decomposed result.
3697*            index  -  output vector which records the row permutation
3698*                      affected by the partial pivoting. (output).
3699*            n      -  size of the square matrix. (input).
3700*            MAX_N  -  physical dimension of a (MAX_N x MAX_N).(input).
3701*            Det    -  determinant of the matrix. (output).
3702*
3703*      LUDCMP returns logical .true. if all proceeded happily.
3704* ______________________________________________________________________
3705*
3706      logical function LUDCMP(a,index,n,MAX_N,Det)
3707      include 'DIFFaX.par'
3708*     save
3709*
3710      integer*4 n, MAX_N, index(MAX_N)
3711      real*8 a(MAX_N,MAX_N), Det
3712*
3713      integer*4 L_MAX, i, j, m, row
3714      parameter (L_MAX = 100)
3715      real*8 tiny, tmp(L_MAX), sum, max, tmp2
3716      parameter (tiny = 1.0D-20)
3717*
3718      LUDCMP = .false.
3719      Det = ONE
3720      if(n.gt.L_MAX) then
3721        write(op,400) 'Matrix too large for LUDCMP'
3722        return
3723      endif
3724      do 20 i = 1, n
3725          max = ZERO
3726        do 10 j = 1, n
3727          if(abs(a(i,j)).gt.max) max = abs(a(i,j))
3728   10   continue
3729        if(max.eq.ZERO) goto 200
3730        tmp(i) = ONE / max
3731   20 continue
3732      do 90 j = 1, n
3733        do 40 i = 1, j-1
3734          sum = a(i,j)
3735          do 30 m = 1, i-1
3736            sum = sum - a(i,m) * a(m,j)
3737   30     continue
3738          a(i,j) = sum
3739   40   continue
3740        max = ZERO
3741        do 60 i = j, n
3742          sum = a(i,j)
3743          do 50 m = 1, j-1
3744            sum = sum - a(i,m) * a(m,j)
3745   50     continue
3746          a(i,j) = sum
3747          tmp2 = tmp(i)*abs(sum)
3748          if(abs(tmp2).ge.max) then
3749            row = i
3750            max = tmp2
3751          endif
3752   60   continue
3753        if(j.ne.row) then
3754          do 70 m = 1, n
3755            tmp2 = a(row,m)
3756            a(row,m) = a(j,m)
3757            a(j,m) = tmp2
3758   70     continue
3759          Det = -Det
3760          tmp(row) = tmp(j)
3761        endif
3762        index(j) = row
3763        if(abs(a(j,j)).eq.ZERO) a(j,j) = tiny
3764        tmp2 = ONE / a(j,j)
3765        do 80 i = j+1, n
3766          a(i,j) = a(i,j) * tmp2
3767   80   continue
3768        Det = Det * a(j,j)
3769   90 continue
3770      LUDCMP = .true.
3771      return
3772  200 continue
3773      return
3774  400 format(1x, a)
3775      end
3776*
3777* ______________________________________________________________________
3778* Title: L_STEP
3779* Author: MMJT
3780* Date: 13 Aug 1989
3781* Description:  L_STEP attempts to determine whether or not there
3782* are any sharp peaks, and if so, their spacing. The algorithm inspects
3783* the sum of certain permutations of the stacking vector z_components.
3784* The components are chosen so that their values should be independent
3785* of the layer origins chosen by the user. ie. Rz(i,i) and
3786* Rz(i,j) + Rz(j,i) and Rz(i,j) + Rz(j,k) + Rz(k,i), etc...
3787* In the interests of clarity, the permutations of the indices are
3788* written out explicity, instead of being generated by a subroutine
3789* (which would considerably shorten the code). We stop searching after
3790* combinations of 4 stacking vectors, since then the structure is too
3791* complicated to be attempting to trick DIFFaX into believing that we
3792* have found sharp peaks. Under these circumstances, work out the
3793* diffraction pattern the long, but reliable, way.
3794*
3795*      ARGUMENTS:
3796*            ok  -  logical flag indicating all went well. (output).
3797*
3798*      COMMON VARIABLES:
3799*            uses:  n_layers, there, l_r
3800*
3801*        modifies:  any_sharp
3802*
3803*      L_STEP returns the l-increment that sharp peaks are likely to
3804*      be found at.
3805* ______________________________________________________________________
3806*
3807      real*8 function L_STEP(ok)
3808      include 'DIFFaX.par'
3809      include 'DIFFaX.inc'
3810*
3811      logical ok
3812*
3813      real*8 tmp, z_step
3814      logical YRDSTK, resonant, decided
3815      integer*4 i1, i2, i3, i4
3816*
3817* external function
3818      external YRDSTK
3819*
3820* initialize return value
3821      L_STEP = ZERO
3822*
3823      resonant = .true.
3824      decided  = .false.
3825      z_step   = ZERO
3826*
3827* Check z-components of Rii stacking vectors
3828* if any of the transitions do not exist, set resonant to false.
3829      do 10 i1 = 1, n_layers
3830        if(resonant) then
3831          if(there(i1,i1)) then
3832            decided = .true.
3833            tmp = l_r(3,i1,i1)
3834            resonant = resonant .and. YRDSTK(z_step, tmp, ok)
3835            if(.not.ok) goto 990
3836          endif
3837        endif
3838   10 continue
3839*
3840* Rii terms do not occur (ie. no layer i will stack to another layer i)
3841* We must therefore check z-components of Rij + Rji sequences (i.ne.j).
3842      if((n_layers.gt.1) .and. .not.decided) then
3843        do 20 i1 = 1, n_layers
3844          do 30 i2 = i1 + 1, n_layers
3845            if(resonant) then
3846              if(there(i2,i1).and.there(i1,i2)) then
3847                decided = .true.
3848                tmp = l_r(3,i2,i1) + l_r(3,i1,i2)
3849                resonant = resonant.and.YRDSTK(z_step, tmp, ok)
3850                if(.not.ok) goto 991
3851              endif
3852            endif
3853   30   continue
3854   20 continue
3855      endif
3856*
3857* No Rij + Rji sequences occur.
3858* Check z-components of Rij + Rjk + Rki sequences (where i.ne.j.ne.k).
3859      if((n_layers.gt.2) .and. .not.decided) then
3860        do 40 i1 = 1, n_layers
3861          do 50 i2 = i1 + 1, n_layers
3862            do 60 i3 = i2 + 1, n_layers
3863              if(resonant) then
3864* There are 2 permutations
3865                if(there(i2,i1).and.
3866     |             there(i3,i2).and.
3867     |             there(i1,i3)) then
3868                   decided = .true.
3869                   tmp = l_r(3,i2,i1) + l_r(3,i3,i2) + l_r(3,i1,i3)
3870                   resonant = resonant .and. YRDSTK(z_step, tmp, ok)
3871                   if(.not.ok) goto 992
3872                endif
3873                if(there(i3,i1).and.
3874     |             there(i2,i3).and.
3875     |             there(i1,i2).and.resonant) then
3876                   decided = .true.
3877                   tmp = l_r(3,i3,i1) + l_r(3,i2,i3) + l_r(3,i1,i2)
3878                   resonant = resonant .and. YRDSTK(z_step, tmp, ok)
3879                   if(.not.ok) goto 993
3880                endif
3881              endif
3882   60       continue
3883   50     continue
3884   40   continue
3885      endif
3886*
3887* No Rij + Rjk + Rki sequences occur.
3888* Check z-components of Rij + Rjk + Rkl + Rli sequences
3889* (where i.ne.j.ne.k.ne.l).
3890      if((n_layers.gt.3) .and. .not.decided) then
3891        do 70 i1 = 1, n_layers
3892          do 80 i2 = i1 + 1, n_layers
3893            do 90 i3 = i2 + 1, n_layers
3894              do 100 i4 = i3 + 1, n_layers
3895                if(resonant) then
3896* There are 6 permutations
3897                  if(there(i2,i1).and.
3898     |               there(i3,i2).and.
3899     |               there(i4,i3).and.
3900     |               there(i1,i4)) then
3901                       decided = .true.
3902                       tmp = l_r(3,i2,i1) + l_r(3,i3,i2) +
3903     |                       l_r(3,i4,i3) + l_r(3,i1,i4)
3904                       resonant = resonant.and.YRDSTK(z_step,tmp,ok)
3905                       if(.not.ok) goto 994
3906                  endif
3907                  if(there(i2,i1).and.
3908     |               there(i4,i2).and.
3909     |               there(i3,i4).and.
3910     |               there(i1,i3).and.resonant) then
3911                       decided = .true.
3912                       tmp = l_r(3,i2,i1) + l_r(3,i4,i2) +
3913     |                        l_r(3,i3,i4) + l_r(3,i1,i3)
3914                       resonant = resonant.and.YRDSTK(z_step,tmp,ok)
3915                       if(.not.ok) goto 995
3916                  endif
3917                  if(there(i3,i1).and.
3918     |               there(i2,i3).and.
3919     |               there(i4,i2).and.
3920     |               there(i1,i4).and.resonant) then
3921                       decided = .true.
3922                       tmp = l_r(3,i3,i1) + l_r(3,i2,i3) +
3923     |                       l_r(3,i4,i2) + l_r(3,i1,i4)
3924                       resonant = resonant.and.YRDSTK(z_step,tmp,ok)
3925                       if(.not.ok) goto 996
3926                  endif
3927                  if(there(i3,i1).and.
3928     |               there(i4,i3).and.
3929     |               there(i2,i4).and.
3930     |               there(i1,i2).and.resonant) then
3931                       decided = .true.
3932                       tmp = l_r(3,i3,i1) + l_r(3,i4,i3) +
3933     |                       l_r(3,i2,i4) + l_r(3,i1,i2)
3934                       resonant = resonant.and.YRDSTK(z_step,tmp,ok)
3935                       if(.not.ok) goto 997
3936                  endif
3937                  if(there(i4,i1).and.
3938     |               there(i2,i4).and.
3939     |               there(i3,i2).and.
3940     |               there(i1,i3).and.resonant) then
3941                       decided = .true.
3942                       tmp = l_r(3,i4,i1) + l_r(3,i2,i4) +
3943     |                       l_r(3,i3,i2) + l_r(3,i1,i3)
3944                       resonant = resonant.and.YRDSTK(z_step,tmp,ok)
3945                       if(.not.ok) goto 998
3946                  endif
3947                  if(there(i4,i1).and.
3948     |               there(i3,i4).and.
3949     |               there(i2,i3).and.
3950     |               there(i1,i2).and.resonant) then
3951                       decided = .true.
3952                       tmp = l_r(3,i4,i1) + l_r(3,i3,i4) +
3953     |                       l_r(3,i2,i3) + l_r(3,i1,i2)
3954                       resonant = resonant.and.YRDSTK(z_step,tmp,ok)
3955                       if(.not.ok) goto 999
3956                  endif
3957                endif
3958  100         continue
3959   90       continue
3960   80     continue
3961   70   continue
3962      endif
3963*
3964* If there is no stacking sequence that can bring us back to a layer
3965* similar to that at the origin after 4 layers, then this
3966* structure is sufficiently complicated that we may be better
3967* off doing adaptive integration anyway. (d_l still equals 0.0.)
3968*
3969      if(decided.and.resonant .and. (tmp.ne.ZERO)) then
3970        L_STEP = ONE / tmp
3971        any_sharp = .true.
3972      else
3973        L_STEP = ZERO
3974        any_sharp = .false.
3975      endif
3976*
3977      return
3978  990 write(op,200)
3979      write(op,201) i1, i1
3980      return
3981  991 write(op,202)
3982      write(op,203) i1, i2, i2, i1
3983      return
3984  992 write(op,202)
3985      write(op,204) i1, i2, i2, i3, i3, i1
3986      return
3987  993 write(op,202)
3988      write(op,204) i1, i3, i3, i2, i2, i1
3989      return
3990  994 write(op,202)
3991      write(op,205) i1, i2, i2, i3, i3, i4, i4, i1
3992      return
3993  995 write(op,202)
3994      write(op,205) i1, i2, i2, i4, i4, i3, i3, i1
3995      return
3996  996 write(op,202)
3997      write(op,205) i1, i3, i3, i2, i2, i4, i4, i1
3998      return
3999  997 write(op,202)
4000      write(op,205) i1, i3, i3, i4, i4, i2, i2, i1
4001      return
4002  998 write(op,202)
4003      write(op,205) i1, i4, i4, i2, i2, i3, i3, i1
4004      return
4005  999 write(op,202)
4006      write(op,205) i1, i4, i4, i3, i3, i2, i2, i1
4007      return
4008  200 format(1x,'L_STEP: Non-physical z-component of stacking vector.')
4009  201 format(1x,'Rz(',i2,',',i2,') = 0.0')
4010  202 format(1x,'L_STEP:Non-physical z-components of stacking vectors')
4011  203 format(1x,'Rz(',i2,',',i2,') + Rz(',i2,',',i2,') = 0.0')
4012  204 format(1x,'Rz(',i2,',',i2,')',2(' + Rz(',i2,',',i2,')'),' = 0.0')
4013  205 format(1x,'Rz(',i2,',',i2,')',3(' + Rz(',i2,',',i2,')'),' = 0.0')
4014      end
4015*
4016* ______________________________________________________________________
4017* Title: MATMUL
4018* Author: MMJT
4019* Date: 5 Feb 1990
4020* Description:  This subroutine multiplies the complex matrices
4021* 'a' and 'b', of logical size n x n. Result is returned in 'a'.
4022*
4023*      ARGUMENTS:
4024*            a   -  Complex*16 array to store result. (input and output).
4025*            b   -  Complex*16 array. (input).
4026*            n   -  Logical size of matrices. (input).
4027*
4028* ______________________________________________________________________
4029*
4030      subroutine MATMUL(a, b, n)
4031      include 'DIFFaX.par'
4032*     save
4033*
4034      integer*4 n
4035      complex*16 a(MAX_L,MAX_L), b(MAX_L,MAX_L)
4036*
4037      integer*4 i, j, m
4038      complex*16 c(MAX_L,MAX_L), ctmp
4039*
4040* first copy a into c
4041      do 10 j = 1, n
4042        do 20 i = 1, n
4043          c(i,j) = a(i,j)
4044   20   continue
4045   10 continue
4046*
4047      do 30 j = 1, n
4048        do 40 i = 1, n
4049          ctmp = C_ZERO
4050          do 50 m = 1, n
4051            ctmp = ctmp + c(i,m) * b(m,j)
4052   50     continue
4053          a(i,j) = ctmp
4054   40   continue
4055   30 continue
4056*
4057      return
4058      end
4059*
4060* ______________________________________________________________________
4061* Title: MAT2N
4062* Author: MMJT
4063* Date: 5 Feb 1990
4064* Description:  This function multiplies the matrix 'mat' by itself
4065* n = l_cnt+1 times. In order to speed up the process, n has been broken
4066* down into its binary representation (by the function BINPOW, which is
4067* called once in OPTIMZ), and the result is given by
4068*
4069*         mat**n = mat**n0 * mat**n1 * mat**n2 * mat**n3 * etc
4070*
4071* where ni = 2**i
4072* and n = n0 + n1 + n2 + n3 + . . .
4073*
4074* mat**ni is given by (mat**n(i-1)) * (mat**n(i-1))
4075* ie. mat**8 = (mat**4) * (mat**4)
4076* and similarly mat**4 = (mat**2) * (mat**2)
4077*
4078* n must be such that n <= RCSV_MAX+1 <= 2**(MAX_BIN+1) - 1
4079*
4080*      ARGUMENTS:
4081*            a   -  Complex*16 array to store result. (output).
4082*
4083*      COMMON VARIABLES:
4084*            uses:  n_layers, pow, max_pow
4085*
4086*        modifies:  No COMMON variable are modified.
4087*
4088*      MAT2N returns TRUE if all went well.
4089* ______________________________________________________________________
4090*
4091      logical function MAT2N(a)
4092      include 'DIFFaX.par'
4093      include 'DIFFaX.inc'
4094*
4095      complex*16 a(MAX_L,MAX_L)
4096*
4097      integer*4 i, j
4098      complex*16 tmp_mat(MAX_L,MAX_L,MAX_BIN)
4099*
4100* external subroutine (Some compilers need them declared external)
4101      external MATSQR, MATMUL
4102*
4103      MAT2N = .false.
4104*
4105* copy mat into the first 2-dimensional tmp_mat array. Initialize a
4106* to be the identity matrix.
4107      do 20 j = 1, n_layers
4108        do 30 i = 1, n_layers
4109          a(i,j) = C_ZERO
4110          tmp_mat(i,j,1) = mat(i,j)
4111   30   continue
4112        a(j,j) = C_ONE
4113   20 continue
4114*
4115      do 40 i = 1, max_pow-1
4116        if(pow(i).eq.1) then
4117          call MATMUL(a, tmp_mat(1,1,i), n_layers)
4118        endif
4119        call MATSQR(tmp_mat(1,1,i+1), tmp_mat(1,1,i), n_layers)
4120   40 continue
4121      if(pow(max_pow).eq.1)
4122     |         call MATMUL(a, tmp_mat(1,1,max_pow), n_layers)
4123*
4124      MAT2N = .true.
4125*
4126      return
4127      end
4128*
4129* ______________________________________________________________________
4130* Title: MATSQR
4131* Author: MMJT
4132* Date: 5 Feb 1990
4133* Description:  This subroutine multiplies the complex matrix
4134* 'b', of logical size n x n, by itself. Result is returned in 'a'.
4135*
4136*      ARGUMENTS:
4137*            a   -  Complex*16 array to store result. (output).
4138*            b   -  Complex*16 array to be 'squared'. (input).
4139*            n   -  Logical size of matrices. (input).
4140*
4141* ______________________________________________________________________
4142*
4143      subroutine MATSQR(a, b, n)
4144      include 'DIFFaX.par'
4145*     save
4146*
4147      integer*4 n
4148      complex*16 a(MAX_L,MAX_L), b(MAX_L,MAX_L)
4149*
4150      integer*4 i, j, m
4151      complex*16 ctmp
4152*
4153      do 10 j = 1, n
4154        do 20 i = 1, n
4155          ctmp = C_ZERO
4156          do 30 m = 1, n
4157            ctmp = ctmp + b(i,m) * b(m,j)
4158   30     continue
4159          a(i,j) = ctmp
4160   20   continue
4161   10 continue
4162*
4163      return
4164      end
4165*
4166* ______________________________________________________________________
4167* Title: NMCOOR
4168* Author: MWD
4169* Date: 18 Aug 1988
4170* Description:  This subroutine multiplies the relative coordinates
4171* of each atom by 2*pi, which is the useful form for calculating phases.
4172*
4173*      ARGUMENTS:
4174*            No input arguments.
4175*
4176*      COMMON VARIABLES:
4177*            uses:  n_actual, l_n_atoms, a_pos, PI2
4178*
4179*        modifies:  a_pos
4180* ______________________________________________________________________
4181*
4182      subroutine NMCOOR()
4183      include 'DIFFaX.par'
4184      include 'DIFFaX.inc'
4185*
4186      integer*4 i, j
4187*
4188      do 10 i = 1, n_actual
4189        do 20 j = 1, l_n_atoms(i)
4190          a_pos(1,j,i) = a_pos(1,j,i) * PI2
4191          a_pos(2,j,i) = a_pos(2,j,i) * PI2
4192          a_pos(3,j,i) = a_pos(3,j,i) * PI2
4193   20   continue
4194   10 continue
4195*
4196      return
4197      end
4198*
4199* ______________________________________________________________________
4200* Title: OPTIMZ
4201* Author: MWD and MMJT
4202* Date: 8 Apr 1992; 15 Mar 1995; 24 July 1997
4203* Description:  This routine determines if any shortcuts can be taken
4204* in the calculations.
4205*
4206*      ARGUMENTS:
4207*            rootnam  -  The name of the input data file. (input).
4208*            ok       -  logical flag indicating all went well.
4209*                                                      (output).
4210*
4211*      COMMON VARIABLES:
4212*            uses:  a0, b0, d0, n_layers, l_alpha, l_actual, l_n_atoms,
4213*                   a_B, r_B11, r_B22, r_B33, r_B12, r_B23, r_B31,
4214*                   l_symmetry, l_seq, DoSymDump, SymGrpNo, lambda,
4215*                   l_cnt, CENTRO, PI, l_r, n_actual, finite_width
4216*
4217*        modifies:  there, one_B, Bs_zero, same_Bs, only_real, l_rz,
4218*                   same_layer, max_var, no_trials, tolerance, SymGrpNo,
4219*                   lambda, check_sym, has_l_mirror, theta1, theta2
4220*                   h_end, h_start, k_end, k_start, pnt_grp, same_rz
4221*                   xplcit, formfactor, all_Bs_zero,
4222*                   a_B11,a_B22,a_B33,a_B12,a_B23,a_B31
4223* ______________________________________________________________________
4224*
4225      subroutine OPTIMZ(rootnam, ok)
4226      include 'DIFFaX.par'
4227      include 'DIFFaX.inc'
4228*
4229      character*(*) rootnam
4230      logical ok
4231*
4232      character*31 sym_fnam
4233      logical EQUALB, BINPOW, GET_G
4234      integer*4 GET_SYM, i, j, j2, m, n, LENGTH
4235      real*8 HKANGL, h_val, k_val
4236      real*8 x, error, tmp, incr, z, old_lambda
4237      logical did_it(MAX_L,MAX_L)
4238*
4239* external functions
4240      external GET_SYM, LENGTH, EQUALB, BINPOW, GET_G
4241* external subroutines (Some compilers need them declared external)
4242*      external GETFNM, GET_BDS, CHK_SYM, NMCOOR, OVERLP
4243*
4244* statement function
4245* HKANGL is the angle between the vector (h_val,k_val,0) and (1,0,0)
4246      HKANGL(k_val,h_val) = atan2(k_val*sqrt(a0*b0 - d0*d0*QUARTER),
4247     |                              (h_val*a0 + k_val*d0*HALF))
4248*
4249* set up logic table for stacking transitions
4250      do 10 i = 1, n_layers
4251        do 20 j = 1, n_layers
4252          there(j,i) = l_alpha(j,i).ge.eps7
4253   20   continue
4254   10 continue
4255*
4256* see if there are any overlapping atoms
4257      call OVERLP()
4258*
4259* multiply all atom coordinates by 2*PI
4260      call NMCOOR()
4261*
4262* If calculation is to be recursive for a finite number of layers,
4263* then store the binary form of l_cnt+1 in an array for efficient
4264* matrix multiplication.
4265      ok = GET_G()
4266      if(recrsv .and. .not.inf_thick) then
4267        ok = BINPOW(l_cnt+1)
4268        if(.not.ok) then
4269          write(op,202) 'ERROR returned by BINPOW to OPTIMZ'
4270          write(op,201) 'The argument passed was l_cnt+1 = ', l_cnt+1
4271          goto 999
4272        endif
4273      endif
4274*
4275* see if Debye-Waller coefficients are same for all atoms in a layer
4276      do 30 i = 1, n_layers
4277        x = ZERO
4278        j2 = l_actual(i)
4279        m = l_n_atoms(j2)
4280        do 40 j = 1, m
4281          x = x + a_B(j,j2)
4282   40   continue
4283        x = x / dble(m)
4284        error = ZERO
4285* find absolute deviation of coefficients
4286        do 50 j = 1, m
4287          error = error + abs(a_B(j,j2) - x)
4288   50   continue
4289* get relative error
4290        if(x.ne.ZERO) error = error / (x * m)
4291        one_B(j2) = abs(error).le.eps3
4292   30 continue
4293*
4294* Check that the layer uncertainty factors are physically reasonable.
4295      do 60 i = 1, n_layers
4296        do 65 j = 1, n_layers
4297          if(there(j,i)) then
4298* check on r_B12
4299            x = r_B11(j,i)*r_B22(j,i)*a0*b0 -
4300     |          r_B12(j,i)*r_B12(j,i)*ab0*ab0
4301            if(x.lt.ZERO) then
4302              write(op,500) 'C12'
4303              write(op,501) i, j
4304              write(op,502) 'C11 and C22'
4305              ok = .false.
4306              goto 999
4307            endif
4308* check on r_B23
4309            x = r_B22(j,i)*r_B33(j,i)*b0*c0 -
4310     |          r_B23(j,i)*r_B23(j,i)*bc0*bc0
4311            if(x.lt.ZERO) then
4312              write(op,500) 'C23'
4313              write(op,501) i, j
4314              write(op,502) 'C22 and C33'
4315              ok = .false.
4316              goto 999
4317            endif
4318* check on r_B31
4319            x = r_B11(j,i)*r_B33(j,i)*a0*c0 -
4320     |          r_B31(j,i)*r_B31(j,i)*ca0*ca0
4321            if(x.lt.ZERO) then
4322              write(op,500) 'C13'
4323              write(op,501) i, j
4324              write(op,502) 'C11 and C33'
4325              ok = .false.
4326              goto 999
4327            endif
4328          endif
4329   65   continue
4330   60 continue
4331*
4332* see if stacking 'uncertainty' coefficients are same for all layers.
4333* Special flag if they are all zero.
4334      all_Bs_zero = .true.
4335      do 70 i = 1, n_layers
4336        do 80 j = 1, n_layers
4337          Bs_zero(j,i) =
4338     |         r_B11(j,i).eq.ZERO .and. r_B22(j,i).eq.ZERO .and.
4339     |         r_B33(j,i).eq.ZERO .and. r_B12(j,i).eq.ZERO .and.
4340     |         r_B23(j,i).eq.ZERO .and. r_B31(j,i).eq.ZERO
4341          all_Bs_zero = all_Bs_zero .and. Bs_zero(j,i)
4342   80   continue
4343   70 continue
4344*
4345* run through all 6 coefficients
4346* until it is clear that some are different
4347                  same_Bs = EQUALB(r_B11, a_B11)
4348      if(same_Bs) same_Bs = EQUALB(r_B22, a_B22)
4349      if(same_Bs) same_Bs = EQUALB(r_B33, a_B33)
4350      if(same_Bs) same_Bs = EQUALB(r_B12, a_B12)
4351      if(same_Bs) same_Bs = EQUALB(r_B23, a_B23)
4352      if(same_Bs) same_Bs = EQUALB(r_B31, a_B31)
4353*
4354* see if all layers are centrosymmetric
4355      only_real = .true.
4356      do 90 i = 1, n_actual
4357        only_real = only_real .and. (l_symmetry(i).eq.CENTRO)
4358   90 continue
4359*
4360* see if all z-components of the stacking vectors are the same
4361      l_rz = ZERO
4362      n = 0
4363      do 100 i = 1, n_layers
4364        do 110 j = 1, n_layers
4365          if(there(j,i)) then
4366            l_rz = l_rz + l_r(3,j,i)
4367            n = n + 1
4368          endif
4369  110   continue
4370  100 continue
4371      l_rz = l_rz / dble(n)
4372      error = ZERO
4373      do 120 i = 1, n_layers
4374        do 130 j = 1, n_layers
4375          if(there(j,i)) error = error + abs(l_r(3,j,i) - l_rz)
4376  130   continue
4377  120 continue
4378      same_rz = abs(error).le.eps4
4379*
4380* If the stacking is explicit, check to see if all the layers are
4381* the same
4382      same_layer = .false.
4383      if(xplcit) then
4384        if(l_cnt.eq.1) goto 140
4385        same_layer = .true.
4386        j = l_seq(1)
4387        i = 2
4388  150   if(l_seq(i).eq.j) then
4389          i = i + 1
4390          if(i.le.l_cnt) goto 150
4391        else
4392          same_layer = .false.
4393        endif
4394*
4395* Check if any of the layer transitions have non-zero probability
4396* initialize flags so that we do not produce duplicate error messages
4397        do 160 i = 1, n_layers
4398          do 170 j = 1, n_layers
4399            did_it(j,i) = .false.
4400  170     continue
4401  160   continue
4402*
4403* now check for legal transitions
4404        write(op,400)
4405        do 180 n = 1, l_cnt-1
4406          i = l_seq(n)
4407          j = l_seq(n+1)
4408          if(.not.there(j,i)) then
4409            ok = .false.
4410            if(.not.did_it(j,i)) then
4411              did_it(j,i) = .true.
4412              write(op,401) j, i
4413            endif
4414          endif
4415  180   continue
4416  140   continue
4417      endif
4418      if(.not.ok) goto 999
4419*
4420* Pre-compute a pseudo-Lorentzian form factor for lateral
4421* planar widths.
4422*
4423      if(finite_width) then
4424* ideally, FFACT_SIZE is odd
4425        m = FFACT_SIZE/2
4426* incr contains the correct increment size such that the first and
4427* last array elements contain zero.
4428        incr = (THREE*N_SIGMAS*N_SIGMAS + ONE) / (TWO*N_SIGMAS*m)
4429        ffact_scale = incr
4430        z = ONE + N_SIGMAS*N_SIGMAS
4431        tmp = ONE / (z*z)
4432* put the peak in the middle (if FFACT_SIZE is odd)
4433        formfactor(m+1) = ONE
4434        do 190, n = 1, m
4435          z = n*incr
4436          if(z.le.dble(N_SIGMAS)) then
4437* Lorentzian part
4438            x = ONE / (ONE + z*z)
4439          else
4440* Linear part
4441             x = tmp*(THREE*N_SIGMAS*N_SIGMAS + ONE - TWO*N_SIGMAS*z)
4442          endif
4443          if(m+n+1.le.FFACT_SIZE) formfactor(m+n+1) = x
4444          if(m-n+1.gt.0)          formfactor(m-n+1) = x
4445  190   continue
4446* compute half width in reciprocal Angstroms for use in computing the
4447* subtle impact of a-b shape broadening on 00l reflections
4448        tmp = Wa*sin(PI-cell_gamma)
4449        ffwdth = sqrt(ONE/(tmp*tmp) + ONE/(Wb*Wb))
4450      endif
4451*
4452* get diffraction symmetry
4453*
4454* establish some bounds
4455      no_trials = 25
4456      max_var = ZERO
4457* save lambda, HKL_LIM (called by THRESH), may change it.
4458      old_lambda = lambda
4459*
4460      call THRESH(ok)
4461      if(.not.ok) goto 999
4462*
4463      if(SymGrpNo.eq.11) then
4464        write(op,202)
4465     |      'Axial integration only selected.'
4466      else
4467        write(op,202)
4468     |      'Evaluating point group symmetry of diffraction data. . .'
4469      endif
4470*      if(DoSymDump) then
4471*        call GETFNM(rootnam, sym_fnam, 'sym', ok)
4472*        if(.not.ok) then
4473*          write(op,202) 'OPTIMZ: ERROR in creating symmetry dumpfile.'
4474*          goto 999
4475*        endif
4476*        if(sy.ne.op)
4477*     |      open(unit = sy, file = sym_fnam, status = 'new')
4478*        write(op,204) 'Writing symmetry data to file ''',
4479*     |                 sym_fnam(1:LENGTH(sym_fnam)),'''. . .'
4480*        write(sy,303) '''', rootnam(1:LENGTH(rootnam)),''''
4481*        write(sy,203) no_trials
4482*        write(sy,206) tiny_inty
4483*      endif
4484*
4485      check_sym = .false.
4486      if(SymGrpNo.eq.UNKNOWN) then
4487        SymGrpNo = GET_SYM(ok)
4488        if(.not.ok) goto 999
4489        write(op,200) 'Diffraction point symmetry is ',pnt_grp
4490        if(SymGrpNo.ne.1) then
4491          write(op,201) '  to within a tolerance of one part in ',
4492     |                  nint(ONE / tolerance)
4493        endif
4494      else
4495        check_sym = .true.
4496        call CHK_SYM(ok)
4497        if(.not.ok) goto 999
4498      endif
4499*
4500* restore lambda.
4501      lambda = old_lambda
4502*
4503      has_l_mirror = SymGrpNo.ne.1 .and. SymGrpNo.ne.3 .and.
4504     |               SymGrpNo.ne.5 .and. SymGrpNo.ne.6 .and.
4505     |               SymGrpNo.ne.11
4506*
4507      if(DoSymDump) then
4508        write(sy,202) ' '
4509        write(sy,204)
4510     | 'The diffraction data fits the point group symmetry ''',
4511     | pnt_grp(1:LENGTH(pnt_grp)),''''
4512        if(SymGrpNo.ne.1 .and. SymGrpNo.ne.11) then
4513          if(max_var.gt.eps6 .and. max_var.le.eps1) then
4514            write(sy,201) '  with a tolerance of one part in ',
4515     |                  nint(ONE / max_var)
4516          else if(max_var.gt.eps1) then
4517            write(sy,205) '  with a tolerance of one part in ',
4518     |                  ONE / max_var
4519          else
4520            write(sy,202)
4521     |      '  with a tolerance better than one part in a million.'
4522          endif
4523        else
4524          write(sy,202)
4525     |'By definition, all diffraction data has a center of symmetry'
4526          write(sy,202) 'thus, there is no need to test for inversion.'
4527        endif
4528* close file, unless the output was to the default device
4529        if(sy.ne.op) close (unit = sy)
4530      endif
4531* establish integration limits and weighting factors
4532      call GET_BDS()
4533* compute angles of scanning vectors relative to 1 0 0
4534      theta1 = HKANGL(k_start, h_start)
4535      theta2 = HKANGL(k_end, h_end)
4536* resolve ambiguity in the case of -1 point symmetry
4537      if(SymGrpNo.eq.1 .or. SymGrpNo.eq.11) then
4538        theta1 = -PI
4539        theta2 =  PI
4540      endif
4541  999 return
4542  200 format(1x, 2a)
4543  201 format(1x, a, i6)
4544  202 format(1x, a)
4545  203 format(1x, 'Number of trials per symmetry element = ', i4)
4546  204 format(1x, 3a)
4547  205 format(1x, a, f3.1)
4548  206 format(1x, "Threshold intensity = ", g22.6)
4549  302 format(1x, 'SYMMETRY EVALUATIONS FOR DATA IN FILE ''', a, '''')
4550  303 format(1x, 'SYMMETRY EVALUATIONS FOR DATA IN FILE ', 3a)
4551  400 format(1x,'Checking for conflicts in layer stackings . . .')
4552  401 format(1x,'ERROR: Layer ',i2,' cannot stack after layer ',i2)
4553  500 format(1x,'ERROR: Non-physical interlayer uncertainty factor ',a)
4554  501 format(1x,'       for stacking from layer ', i2, ' to ', i2)
4555  502 format(1x,'       It is too large relative to ', a)
4556      end
4557*
4558* ______________________________________________________________________
4559* Title: OVERLP
4560* Authors: MMJT
4561* Date: 24 Feb 1990
4562* Description: This function compares the coordinates of atoms in
4563* adjacent layers, and searches for any overlap. If it finds any
4564* overlap, it checks the summed occupancy to check that it does
4565* not exceed 1. If OVERLP finds atoms too close it provides a warning
4566* message only.
4567*
4568*      ARGUMENTS:
4569*            No arguments are passed.
4570*
4571*      COMMON VARIABLES:
4572*            uses:  n_layers, there, l_n_atoms, l_r, a_pos
4573*
4574*        modifies:  No COMMON variables are modified
4575* ______________________________________________________________________
4576*
4577      subroutine OVERLP()
4578      include 'DIFFaX.par'
4579      include 'DIFFaX.inc'
4580*
4581      logical invert
4582      character*33 txt
4583      integer*4 i, j, m, n, nn, j2, err_no, max_err, fact, at_num
4584      integer*4 PRUNE
4585      parameter(max_err = 100)
4586      real*8 lay(3,2*MAX_A)
4587      real*8 x1, y1, z1, x2, y2, z2, sum_occ, tol, tmp, BOUNDS
4588      parameter(tol = eps1)
4589*
4590* external functions
4591      external BOUNDS, PRUNE
4592*
4593      write(op,400) 'Checking for conflicts in atom positions . . .'
4594*
4595      err_no = 0
4596      do 10 i = 1, n_layers
4597        fact = 1
4598        invert = l_symmetry(i).eq.CENTRO
4599        if(invert) fact = 2
4600        at_num = l_n_atoms(i)
4601        do 20 j2 = 1, at_num
4602          lay(1,j2) = BOUNDS(a_pos(1,j2,i))
4603          lay(2,j2) = BOUNDS(a_pos(2,j2,i))
4604* Remember, only the a and b directions are truly periodic.
4605* The scaling along c is arbitrary. Furthermore, c is perpendicular
4606* to a and b, and is not necessarily parallel to Rii, which (if it
4607* exists) would define the third cell-repeat direction. In other words,
4608* for the general case, we cannot use the BOUNDS function along c.
4609          lay(3,j2) = a_pos(3,j2,i)
4610   20   continue
4611        if(invert) then
4612          do 30 j2 = 1, at_num
4613            lay(1,at_num+j2) = BOUNDS(-a_pos(1,j2,i))
4614            lay(2,at_num+j2) = BOUNDS(-a_pos(2,j2,i))
4615            lay(3,at_num+j2) = -a_pos(3,j2,i)
4616   30     continue
4617        endif
4618        do 40 m = 1, at_num
4619          x1 = lay(1,m)
4620          y1 = lay(2,m)
4621          z1 = lay(3,m)
4622          do 50 n = m + 1, fact * at_num
4623            if(n.gt.at_num) then
4624              nn = n - at_num
4625            else
4626              nn = n
4627            endif
4628            x2 = lay(1,n)
4629            y2 = lay(2,n)
4630            z2 = lay(3,n)
4631            if(abs(x1-x2)*cell_a.gt.tol) goto 50
4632            if(abs(y1-y2)*cell_b.gt.tol) goto 50
4633            if(abs(z1-z2)*cell_c.le.tol) then
4634              sum_occ = a_occup(nn,i) + a_occup(m,i)
4635              if((sum_occ - ONE).gt.eps4) then
4636                if(n.le.at_num) then
4637                  txt = 'are too close in layer'
4638                else
4639                  txt = '(inverted) are too close in layer'
4640                endif
4641                write(op,410) 'Atom ', a_name(nn,i), a_number(nn,i),
4642     |                   ' and atom ', a_name(m,i), a_number(m,i)
4643                write(op,412) txt(1:PRUNE(txt)), i
4644                write(op,420) 'Their combined occupancy is ', sum_occ
4645                err_no = err_no + 1
4646                if(err_no.gt.max_err) goto 999
4647              endif
4648            endif
4649   50     continue
4650   40   continue
4651   10 continue
4652*
4653* now let's look at i-j layer transitions and generate a simple warning
4654* message if it seems that the layers are intertwined.
4655      do 60 i = 1, n_layers
4656        do 70 j = 1, n_layers
4657          if(there(j,i) .and. i.ne.j) then
4658            tmp = l_r(3,j,i) +
4659     |              low_atom(l_actual(j)) - high_atom(l_actual(i))
4660            if(tmp*cell_c.le.-tol) then
4661              write(op,430) 'Atoms from layer', j,
4662     |                    ' extend into layer', i
4663            endif
4664          endif
4665   70   continue
4666   60 continue
4667*
4668      return
4669  999 write(op,401) 'WARNING: Number of errors exceeds ', max_err
4670      return
4671  400 format(1x, a)
4672  401 format(1x, a, i5)
4673  410 format(1x, 'WARNING: ', 2(2a, ' (number ', i3, ')' ) )
4674  412 format(10x, a, 1x, i2)
4675  420 format(1x, a, g12.5)
4676  430 format(1x, 'WARNING: ', 2(a, i3))
4677      end
4678*
4679*
4680* ______________________________________________________________________
4681* Title: PNTINT
4682* Author: MMJT
4683* Date: 30 July 1989
4684* Gets the intensity at the point h, k, l. This differs from the
4685* subroutine POINT only in that PNTINT does not interact with the user.
4686*
4687*      ARGUMENTS:
4688*            h   -  reciprocal vector h-component. (input).
4689*            k   -  reciprocal vector k-component. (input).
4690*            l   -  reciprocal lattice vector l-component. (input).
4691*            ok  -  logical flag indicating all went well. (output).
4692*
4693*      COMMON VARIABLES:
4694*            uses:  a0, b0, c0, d0, lambda, recrsv, rad_type, X_RAY
4695*
4696*        modifies:  no COMMON variables are modified
4697*
4698*      PNTINT returns the intensity at h, k, l
4699* ______________________________________________________________________
4700*
4701      real*8 function PNTINT(h, k, l, ok)
4702      include 'DIFFaX.par'
4703      include 'DIFFaX.inc'
4704*
4705      integer*4 h, k
4706      real*8 l
4707      logical ok
4708*
4709      real*8 S, ANGLE, W4, INTENS, INTEN2, theta, x
4710      complex*16 f(MAX_L)
4711*
4712* external functions
4713      external INTENS, INTEN2
4714* external subroutines (Some compilers need them declared external)
4715*      external GET_F, XYPHSE, PRE_MAT
4716*
4717* statement functions
4718* S is the value of 1/d**2 at hkl
4719      S(h,k,l) = h*h*a0 + k*k*b0 + l*l*c0 + h*k*d0
4720* ANGLE is the Bragg angle (in radians) of the h,k,l plane
4721      ANGLE(h,k,l) = asin(HALF * lambda * sqrt(S(h,k,l)))
4722* W4 is the X-ray polarization factor
4723      W4(theta) = HALF * (ONE + (cos(TWO*theta))**2)
4724*
4725      call XYPHSE(h, k)
4726      call PRE_MAT(h, k)
4727      call GET_F(f, S(h,k,l), l)
4728      if(recrsv) then
4729        x = INTENS(f, h, k, l, ok)
4730      else
4731        x = INTEN2(f, h, k, l, ok)
4732      endif
4733      if(.not.ok) then
4734        if(recrsv) then
4735          write(op,200) 'INTENS'
4736        else
4737          write(op,200) 'INTEN2'
4738        endif
4739        write(op,201) 'h,k,l = ', h,',', k,',', l
4740      endif
4741      if(rad_type.eq.X_RAY)  x = x * W4(ANGLE(h,k,l))
4742*
4743      PNTINT = x
4744*
4745      return
4746  200 format(1x, 'ERROR returned from ', a, ' in subroutine PNTINT')
4747  201 format(1x, 2(a, i3), a, g12.5)
4748      end
4749*
4750* ______________________________________________________________________
4751* Title: POLINT
4752* Author: MMJT, adapted from Numerical Recipes Software
4753* Date: Copyright (C) 1985
4754* Description: Given arrays xa and ya, each of length n, and given
4755* a value x, this routine returns an interpolated estimate for y(x),
4756* and an error estimate dy. If P(x) is the polynomial of degree n-1,
4757* such that P(xa_i) = ya_i, i=1,....N, then the returned value
4758* y = P(x). This routine is called by APPR_F.
4759*
4760*      ARGUMENTS:
4761*            xa  -  Array of length n containing the x values
4762*                   corresponding to the n ya values. (input).
4763*            ya  -  Array of length n containing input y values. (input).
4764*            n   -  Length of the arrays entered. (input).
4765*            x   -  x position at which an interpolated value is
4766*                   required. (input).
4767*            y   -  interpolated y value. (output).
4768*            dy  -  estimate of the error in y. (output).
4769*            ok  -  logical flag indicating all went well. (output).
4770* ______________________________________________________________________
4771*
4772      subroutine POLINT(xa,ya,n,x,y,dy,ok)
4773      include 'DIFFaX.par'
4774*     save
4775*
4776      integer*4 n
4777      real*8 x, xa(n)
4778      complex*16 ya(n), dy, y
4779      logical ok
4780*
4781      integer*4 NMAX, i, m, ns
4782      parameter (NMAX = 10)
4783      real*8 dif, dift, ho, hp
4784      complex*16 c(NMAX), d(NMAX), w, den
4785*
4786      ns = 1
4787      dif = abs(x - xa(1))
4788      do 11 i = 1, n
4789        dift = abs(x - xa(i))
4790        if(dift.lt.dif) then
4791          ns = i
4792          dif = dift
4793        endif
4794        c(i) = ya(i)
4795        d(i) = ya(i)
4796   11 continue
4797      y = ya(ns)
4798      ns = ns - 1
4799      do 13 m = 1, n - 1
4800        do 12 i = 1, n - m
4801          ho = xa(i) - x
4802          hp = xa(i + m) - x
4803          w = c(i+1) - d(i)
4804          den = dcmplx(ho - hp, ZERO)
4805          if(den.eq.C_ZERO) goto 999
4806          den = w / den
4807          d(i) = hp * den
4808          c(i) = ho * den
4809   12   continue
4810        if(2*ns.lt.n-m) then
4811          dy = c(ns + 1)
4812        else
4813          dy = d(ns)
4814          ns = ns - 1
4815        endif
4816        y = y + dy
4817   13 continue
4818*
4819      return
4820  999 ok = .false.
4821      write(op,100) 'ERROR: Zero denominator in POLINT.'
4822      return
4823  100 format(1x, a)
4824      end
4825*
4826* ______________________________________________________________________
4827* Title: PRE_MAT
4828* Author: MMJT
4829* Date: 21 Mar 1990; 21 July 1997
4830* Description:  This subroutine premultiplies some of the factors
4831* needed to calculate the matrix mat ( = Identity - alphaij * Rij)
4832* at (h,k,0) which remain constant when integrating along a streak. The
4833* pre-multiplied factors are stored in mat1, and fatsWalla_hk.
4834*
4835*      ARGUMENTS:
4836*            h   -  reciprocal lattice vector h-component. (input).
4837*            k   -  reciprocal lattice vector k-component. (input).
4838*
4839*      COMMON VARIABLES:
4840*            uses:   n_layers, l_r, there, detune, a0, b0, ab0, r_B11,
4841*                    r_B22, r_B12, same_Bs, Bs_zero, PI2, l_alpha
4842*
4843*        modifies:   l_phi, mat1, fatsWalla_hk
4844* ______________________________________________________________________
4845*
4846      subroutine PRE_MAT(h, k)
4847      include 'DIFFaX.par'
4848      include 'DIFFaX.inc'
4849*
4850      integer*4 h, k
4851*
4852      real*8 dot
4853      integer*4 i, j
4854*
4855* Set up matrix that represents the sequences
4856* For the matrix inversion routines, 'mat' and 'mat1' have to be
4857* in 'i,j' format rather than the quicker 'j,i' format
4858      do 10 i = 1, n_layers
4859        do 20 j = 1, n_layers
4860          if(there(j,i)) then
4861            dot = PI2*(h*l_r(1,j,i) + k*l_r(2,j,i))
4862            l_phi(j,i) = dcmplx(cos(dot),sin(dot))
4863            if(same_Bs.or.Bs_zero(j,i)) then
4864              mat1(i,j) =  detune(j,i) * l_alpha(j,i) * l_phi(j,i)
4865            else
4866* h-k components only. l-components are handled later by GET_MAT.
4867              mat1(i,j) = detune(j,i) * l_alpha(j,i) * l_phi(j,i)
4868     |           * exp(-QUARTER*(r_B11(j,i)*a0*h*h + r_B22(j,i)*b0*k*k)
4869     |           + HALF*r_B12(j,i)*ab0*h*k )
4870            endif
4871          else
4872            mat1(i,j) = C_ZERO
4873          endif
4874   20   continue
4875   10 continue
4876*
4877* Are all the uncertainty factors identical?
4878* Here, we compute only the h-k components.
4879* The l-components are handled later by GET_MAT.
4880      if(same_Bs) then
4881        if(all_Bs_zero) then
4882* This initialization is not actually necessary, since if we are here,
4883* fatsWalla_hk will not be needed by GET_MAT. However, let's be safe.
4884          fatsWalla_hk = ONE
4885        else
4886          fatsWalla_hk =
4887     |      exp(-(QUARTER*(a_B11*a0*h*h + a_B22*b0*k*k) +
4888     |          HALF*a_B12*ab0*h*k) )
4889        endif
4890      endif
4891*
4892      return
4893      end
4894*
4895* ______________________________________________________________________
4896* Title: PRUNE
4897* Author: MMJT
4898* Date: 4 Oct 1989
4899* Description: This function determines the position of the last
4900* non-blank character in a character string.
4901*
4902*      ARGUMENTS:
4903*            line     -  Line of characters to examine. (input).
4904*
4905* PRUNE returns the position of the last non-blank character, or zero
4906* if there was an error.
4907* ______________________________________________________________________
4908*
4909      integer*4 function PRUNE(line)
4910      include 'DIFFaX.par'
4911*     save
4912*
4913      character*(*) line
4914*
4915      integer*4 lin_len, i
4916*
4917      PRUNE = 0
4918*
4919      lin_len = len(line)
4920      i = lin_len
4921   10 if(i.gt.0) then
4922        if(line(i:i).eq.' ') then
4923          if(i.gt.1) then
4924            i = i - 1
4925            goto 10
4926          endif
4927        endif
4928      endif
4929*
4930      if(i.gt.0) PRUNE = i
4931*
4932      return
4933      end
4934*
4935* ______________________________________________________________________
4936* Title: PV
4937* Author: MMJT
4938* Date: 21st Dec 1990; 7 Mar 1995; 9 June 1998
4939* Description: This subroutine performs the pseudo-Voigt
4940* instrumental broadening. Expects the standard u, v and w and gamma
4941* parameters to be available in 'common'. PV does not conserve
4942* intensity well when the broadening width is comparable to d_theta.
4943* Data at the extreme ends of the spectrum are corrupted slightly.
4944*
4945*      ARGUMENTS:
4946*            th2_low  -  lowest 2theta angle to consider. (input).
4947*
4948*      COMMON VARIABLES:
4949*            uses:  pv_u, pv_v, pv_w, pv_gamma, d_theta, spec
4950*                   NONE, PI, RAD2DEG, blurring, th2_max
4951*
4952*        modifies:  brd_spc
4953* ______________________________________________________________________
4954*
4955      subroutine PV(th2_low)
4956      include 'DIFFaX.par'
4957      include 'DIFFaX.inc'
4958*
4959      real*8 th2_low
4960*
4961      integer*4 i, j, n_low, n_high, indx
4962      real*8 th_rng, tn_th, c00, hk_inv, th0
4963      real*8 k1, k2, k3, k4, k5, pVoigt, const, tmp, speci
4964*
4965* first check the numbers
4966      if(pv_u.eq.ZERO .and. pv_v.eq.ZERO .and. pv_w.eq.ZERO) goto 990
4967      if(pv_gamma.lt.ZERO .or. pv_gamma.gt.ONE) goto 999
4968*
4969      if(th2_low.lt.ZERO .or. th2_low.ge.th2_max) then
4970        write(op,103) 'PV: Cut-off angle ', th2_low,
4971     |        ' is out of bounds. Angle reset to zero.'
4972        th2_low = ZERO
4973      endif
4974*
4975* th2_low is the angle relative to th2_min
4976* 2*d_theta is the angular step size
4977      n_low  = int(HALF*th2_low/d_theta) + 1
4978      n_high = int(HALF*(th2_max-th2_min)/d_theta) + 1
4979*
4980      c00 = FOUR * log(TWO)
4981      const = TWO * RAD2DEG * d_theta
4982      do 10 i = 1, n_high
4983        brd_spc(i) = ZERO
4984   10 continue
4985      k1 = pv_gamma*TWO/PI
4986      k2 = (ONE - pv_gamma)*sqrt(c00/PI)
4987      k3 = -c00
4988      th0 = HALF*th2_min
4989      do 20 i = n_low, n_high
4990* get tan((2theta)/2)
4991        tn_th = tan(i*d_theta + th0)
4992        tmp = (pv_u * tn_th + pv_v) * tn_th  +  pv_w
4993        if(tmp.le.ZERO) goto 995
4994        hk_inv = ONE / sqrt(tmp)
4995        tmp = (const * hk_inv)**2
4996        k4 = k1 * hk_inv * const
4997        k5 = k2 * hk_inv * const
4998        speci = spec(i)
4999        do 30 j = n_low - i, n_high - i
5000          th_rng = tmp * j*j
5001          pVoigt = (k4/(ONE+FOUR*th_rng) + k5*exp(k3*th_rng)) * speci
5002          indx = i + j
5003          brd_spc(indx) = brd_spc(indx) + pVoigt
5004   30   continue
5005   20 continue
5006      return
5007  990 write(op,100) 'pseudo-Voigt parameters are zero in PV()'
5008      write(op,101)
5009      blurring = NONE
5010      return
5011  995 write(op,102)
5012     | 'ERROR: pseudo-Voigt spread function is complex at theta = ',
5013     |  i*d_theta
5014      write(op,100) '   u, v, w parameters are illegal.'
5015      write(op,101)
5016      blurring = NONE
5017      return
5018  999 write(op,100) 'Illegal pseudo-Voigt gamma value in PV()'
5019      write(op,101)
5020      blurring = NONE
5021      return
5022  100 format(1x, a)
5023  101 format(1x, '   pseudo-Voigt instrumental broadening not added')
5024  102 format(1x, a, g12.5)
5025  103 format(1x, a, g12.5, a)
5026      end
5027*
5028* ______________________________________________________________________
5029* Title: RAN3
5030* Authors: Press, Flannery, Teukolsky and Vetterling
5031* Date: Copyright (C) 1985
5032* Returns a uniform random deviate between 0.0 and 1.0. Set 'idum'
5033* to any negative value to initialize or reinitialize the sequence.
5034* This version is modified to return real*8 values, and enforces static
5035* storage of all local variables by use of the 'save' statement
5036* (In fact 'seed' is the important variable to save, but we save all
5037* anyway).
5038*
5039*      ARGUMENTS:
5040*            idum       -  Set -ve to initialize. (input).
5041*
5042*      RAN3 returns a real random number between 0 and 1
5043* ______________________________________________________________________
5044*
5045      real*8 function RAN3(idum)
5046      implicit none
5047      save
5048*
5049      integer*4 idum
5050*
5051      real*8 big, seed, mz, fac
5052      parameter (big=4000000.0D0,seed=1618033.0D0,mz=0.0D0,fac=2.5D-7)
5053      real*8 ma(55)
5054      real*8 mj, mk
5055      integer*4 iff, ii, i, j, inext, inextp
5056*
5057      data iff /0/
5058*
5059      if(idum.lt.0 .or. iff.eq.0) then
5060        iff = 1
5061        mj = seed - iabs(idum)
5062        mj = mod(mj,big)
5063        ma(55) = mj
5064        mk = 1.0D0
5065        do 11 i = 1, 54
5066          ii = mod(21*i,55)
5067          ma(ii) = mk
5068          mk = mj - mk
5069          if(mk.lt.mz) mk = mk + big
5070          mj = ma(ii)
5071   11   continue
5072        do 13 j = 1, 4
5073          do 12 i = 1, 55
5074            ma(i) = ma(i) - ma(1 + mod(i + 30,55))
5075            if(ma(i).lt.mz) ma(i) = ma(i) + big
5076   12     continue
5077   13   continue
5078        inext = 0
5079        inextp = 31
5080        idum = 1
5081      endif
5082*
5083      inext = inext + 1
5084      if(inext.eq.56) inext = 1
5085      inextp = inextp + 1
5086      if(inextp.eq.56) inextp = 1
5087      mj = ma(inext) - ma(inextp)
5088      if(mj.lt.mz) mj = mj + big
5089      ma(inext) = mj
5090      RAN3 = mj * fac
5091*
5092      return
5093      end
5094*
5095* ______________________________________________________________________
5096* Title: SHARP
5097* Author: MMJT
5098* Date: 29 Aug 1991
5099* Description: This subroutine determines whether or not
5100* spots on a given h, k row are sharp. It does this by examining
5101* the intensity at l = 0 (or, if absent, at l = d_l) and comparing
5102* with the intensity at a nearby l-value. If the peak is sharp, there
5103* will be a large change in intensity.
5104*
5105*      ARGUMENTS:
5106*            h    -  reciprocal vector h-component. (input).
5107*            k    -  reciprocal vector k-component. (input).
5108*            d_l  -  steps in l-component of the reciprocal lattice
5109*                    vector that sharp spots are likely to be found at.
5110*                                                            (input).
5111*
5112*      COMMON VARIABLES:
5113*            uses:  a0, b0, c0, d0, lambda, PI
5114*
5115*        modifies:  no COMMON variables are modified
5116*
5117*      SHARP returns logical .true. if it thinks h, k contains sharp
5118*      spots.
5119* ______________________________________________________________________
5120*
5121      logical function SHARP(h, k, d_l)
5122      include 'DIFFaX.par'
5123      include 'DIFFaX.inc'
5124*
5125      integer*4 h, k
5126      real*8 d_l
5127*
5128      logical ok
5129      real*8 S, PNTINT, ANGLE, LL, l, i1, i2, x, theta, l_next
5130*
5131* external subroutine (Some compilers need them declared external)
5132*      external GET_F
5133* external functions
5134      external PNTINT
5135*
5136* statement functions
5137* S is the value of 1/d**2 at hkl
5138      S(h,k,l) = h*h*a0 + k*k*b0 + l*l*c0 + h*k*d0
5139* ANGLE is the Bragg angle (in radians) of the h,k,l plane
5140      ANGLE(h,k,l) = asin(HALF * lambda * sqrt(S(h,k,l)))
5141* LL is the maximum allowable l value for a given h,k and theta
5142      LL(theta,h,k) =  sqrt( ( (TWO*sin(theta)/lambda)**2
5143     |                         - h*h*a0 - k*k*b0 - h*k*d0) / c0 )
5144*
5145      SHARP = .false.
5146*
5147* get the intensity at hkl, with l = 0 initially
5148      l = ZERO
5149   10 i1 = PNTINT(h, k, l, ok)
5150      if(.not.ok) goto 100
5151* If there is an extinction at hkl, try again at l = l + d_l
5152      if(i1.lt.eps4) then
5153        l = l + d_l
5154        if(ANGLE(h,k,l).lt.HALF*PI) then
5155          goto 10
5156        else
5157          goto 100
5158        endif
5159      endif
5160*
5161* Define a spot to be sharp if intensity is halved at l = l + d_l/100
5162      theta = ANGLE(h, k, l)
5163      x = min(d_theta, HALF*th2_max-theta)
5164      l_next = LL(theta+x, h, k)
5165      i2 = PNTINT(h, k, l+l_next*eps2, ok)
5166      if(.not.ok) goto 100
5167*
5168      SHARP = i1 .gt. TWO*i2
5169*
5170  100 return
5171      end
5172*
5173* ______________________________________________________________________
5174* Title: SMUDGE
5175* Author: MMJT
5176* Date: 30 Oct 1989; 16 April 1999
5177* Description: This subroutine convolutes 'array', of length
5178* 'arrsize' with a Gaussian of standard deviation 'sigma'.
5179* NOTE: This routine does not conserve area under the convoluted curve.
5180* It is designed so that the peak heights are unchanged.
5181*
5182*      ARGUMENTS:
5183*            array    -  The name of the input array. (input).
5184*            arrsize  -  Size of the input array. (input).
5185*            sigma    -  Standard deviation of the Gaussian as a
5186*                        multiple of the array sampling step. (input).
5187*            ok       -  logical flag indicating all went well.
5188*                                                      (output).
5189* ______________________________________________________________________
5190*
5191      subroutine SMUDGE(array, arrsize, sigma, ok)
5192      include 'DIFFaX.par'
5193      include 'DIFFaX.inc'
5194*     save
5195*
5196      integer*4 arrsize
5197      real*8 array(arrsize), sigma
5198      logical ok
5199*
5200      integer*4 m, i, j
5201      real*8 tmparr(SADSIZE)
5202      real*8 k1, k2, tmp, tmp1, tmp2, gss, normalize
5203*
5204      if(sigma.eq.ZERO) return
5205*
5206      do 10 i = 1, arrsize
5207        if(array(i).gt.maxsad) array(i) = maxsad
5208        tmparr(i) = array(i)
5209        array(i) = ZERO
5210   10 continue
5211*
5212* go out to 5 standard deviations, or to the end of the spectrum
5213      m = nint(FIVE * sigma)
5214      k1 = HALF / ( sigma * sigma )
5215      if(m.gt.arrsize) m = arrsize
5216* get normalization constant. We wish peak heights to be unchanged.
5217      normalize = ONE
5218      do 20 i = 1, m
5219        normalize = normalize + TWO * exp( -k1 * dble(i*i) )
5220   20 continue
5221*
5222      if(normalize.eq.ZERO) then
5223        write(op,100) 'ERROR in SMUDGE: Zero normalization constant.'
5224        ok = .false.
5225        goto 999
5226      endif
5227      normalize = ONE / normalize
5228*
5229      do 30 i = 0, m
5230        k2 = k1*dble(i*i)
5231        gss = exp(-k2)
5232        do 40 j = 1, arrsize
5233          tmp1 = ZERO
5234          tmp2 = ZERO
5235          if((j-i).gt.0) tmp1 = tmparr(j-i)
5236          if((j+i).le.arrsize) tmp2 = tmparr(j+i)
5237          tmp = tmp1 + tmp2
5238          if(i.eq.0) tmp = HALF * tmp
5239          array(j) = array(j) + gss * tmp * normalize
5240   40   continue
5241   30 continue
5242*
5243  999 return
5244  100 format(1x, a)
5245      end
5246*
5247* ______________________________________________________________________
5248* Title: SPHCST
5249* Author: MWD
5250* Date: 18 Aug 1988
5251* Description: This subroutine determines the constants used in
5252* determining the magnitude of reciprocal lattice vectors.
5253*
5254*      ARGUMENTS:
5255*            No input arguments.
5256*
5257*      COMMON VARIABLES:
5258*            uses:  cell_a, cell_b, cell_c, cell_gamma
5259*
5260*        modifies:  a0, b0, c0, d0, ab0, bc0, ca0
5261* ______________________________________________________________________
5262*
5263      subroutine SPHCST()
5264      include 'DIFFaX.par'
5265      include 'DIFFaX.inc'
5266*
5267      a0 = ONE / (cell_a * sin(cell_gamma))**2
5268      b0 = ONE / (cell_b * sin(cell_gamma))**2
5269      c0 = ONE / (cell_c)**2
5270      d0 = -TWO * cos(cell_gamma) /
5271     |            (cell_a * cell_b * sin(cell_gamma)**2)
5272*
5273      ab0 = sqrt(a0 * b0)
5274      bc0 = sqrt(b0 * c0)
5275      ca0 = sqrt(c0 * a0)
5276*
5277      return
5278      end
5279*
5280* ______________________________________________________________________
5281* Title: THRESH
5282* Author: MMJT
5283* Date: 21 Jan 1995
5284* Samples intensities in reciprocal space to get a "feeling" for what
5285* kind of values are out there prior to testing diffraction symmetry.
5286* CHK_SYM and GET_SYM measure the fractional deviations of intensity
5287* from (potentially) symmetry-related points in reciprocal space.
5288* This method runs into problems when the intensity being measured
5289* is close to zero. Miniscule intensity variations can appear to be
5290* huge relative to zero!
5291* This function is needed in order to obtain a (crude) estimate of
5292* which intensity values are too small to worry about, even if the
5293* relative intensity variations seem to be large.
5294* This function will be of no use if there are no streaks, that is,
5295* if the crystal is perfect.
5296*
5297*      ARGUMENTS:
5298*            ok       -  logical flag indicating all went well.
5299*                                                      (output).
5300*
5301*      COMMON VARIABLES:
5302*            uses:  a0, b0, c0, d0, lambda, no_trials,
5303*                   max_angle
5304*
5305*        modifies:  ok, max_angle, h_bnd, k_bnd, l_bnd, tiny_inty
5306*
5307* ______________________________________________________________________
5308*
5309      subroutine THRESH(ok)
5310      include 'DIFFaX.par'
5311      include 'DIFFaX.inc'
5312*
5313      logical ok
5314*
5315      integer*4 i, h, k, idum
5316      real*8 RAN3, PNTINT, S, ANGLE
5317      real*8 l, tot_int
5318*
5319* external functions
5320      external RAN3, PNTINT
5321*
5322* statement functions
5323* S is the value of 1/d**2 at hkl
5324      S(h,k,l) = h*h*a0 + k*k*b0 + l*l*c0 + h*k*d0
5325      ANGLE(h,k,l) = asin(HALF * lambda * sqrt(S(h,k,l)))
5326*
5327* initialize random numbers in RAN3
5328      idum = -1
5329*
5330* First define angular range to sample. h_bnd, k_bnd, l_bnd
5331* (defined in HKL_LIM) and max_angle, are used later on
5332* in GET_SYM and CHK_SYM.
5333      max_angle = QUARTER * PI
5334      call HKL_LIM()
5335*
5336* Sample the typical intensities that are out there
5337      tot_int = ZERO
5338      do 10 i = 1, no_trials
5339   20   h = int( dble(h_bnd + 1) * RAN3(idum) )
5340        k = int( dble(k_bnd + 1) * RAN3(idum) )
5341        l = l_bnd * RAN3(idum)
5342* make sure we are not sampling at too high an angle
5343        if(TWO * ANGLE(h, k, l) .gt. max_angle) goto 20
5344* get I(h,k,l)
5345        tot_int = tot_int + PNTINT(h, k, l, ok)
5346        if(.not.ok) goto 999
5347   10 continue
5348*
5349* Estimate some suitable fraction of the average intensity. This
5350* fraction defines a baseline intensity equivalent to zero for
5351* the purposes of the symmetry testing later on.
5352*
5353      tiny_inty = tot_int * eps5 / no_trials
5354*
5355  900 return
5356  999 write(op,100) 'ERROR in intensity calculation in THRESH'
5357      write(op,200) '   at h,k,l = ', h,',', k,',', l
5358      return
5359  100 format(1x, a)
5360  200 format(1x, a, i3, a, i3, a, f7.2)
5361      end
5362*
5363* ______________________________________________________________________
5364* Title: TRMSPC
5365* Author: MMJT
5366* Date: 20 April 1989; 7 Mar 1995
5367* Description:  This function locates a suitable cut-off angle,
5368* below which we can ignore the huge intensities which may occur
5369* close to the origin when the full adaptive integration option is used.
5370* TRMSPC is called only when theta = 0 is in the integrated range.
5371*
5372*      ARGUMENTS:
5373*            th2_low  -  a 2theta cut-off angle to be found (output)
5374*
5375*      COMMON VARIABLES:
5376*            uses:  th2_min, th2_max, d_theta, spec, RAD2DEG
5377*
5378*        modifies:  spec
5379*
5380*            TRMSPC returns logical .true. if all went well.
5381* ______________________________________________________________________
5382*
5383      logical function TRMSPC(th2_low)
5384      include 'DIFFaX.par'
5385      include 'DIFFaX.inc'
5386*
5387      real*8 th2_low
5388*
5389      integer*4 i, i_min, i_max
5390*
5391      i_max = int(HALF*(th2_max - th2_min) / d_theta) + 1
5392* spec(1) corresponds to the intensity at the origin and is always zero.
5393      i = 2
5394   20 i = i + 1
5395        if(i.ge.i_max+1) then
5396          write(op,100) 'No peaks were found in spectrum.'
5397          goto 900
5398        endif
5399* locate the first minimum after the huge peak at the origin
5400        if(spec(i).le.spec(i-1))
5401     |goto 20
5402      i_min = i - 1
5403*
5404* NOTE: The absolute angle is th2_low + th2_min
5405      th2_low = i_min * d_theta
5406*
5407  900 TRMSPC = .true.
5408*
5409      return
5410  100 format(1x, a)
5411      end
5412*
5413* ______________________________________________________________________
5414* Title: TST_MIR
5415* Author: MMJT
5416* Date: 30 July 1989; 22 Feb 1995
5417* Identifies the presence of mirror planes which contain the streaky
5418* axis. The integer argument 'mir_sym' can have one of two values,
5419* 1 or 2.
5420* mir_sym = 1 is the plane containing the cell sides h and l.
5421* mir_sym = 2 is the plane containing the cell sides k and l.
5422* For rotational symmetry greater than 2, one mirror implies the
5423* presence of the other. For diffraction symmetry group No 3 (2/M),
5424* only one or the other can occur, but we must test for both.
5425* TST_MIR returns '.true.' if a mirror was found, '.false.' otherwise.
5426*
5427*      ARGUMENTS:
5428*            mir_sym  -  Plane across which we wish to test for mirror
5429*                        symmetry. Takes the values 1 or 2. (input).
5430*            idum     -  parameter used by RAN3. Is -ve if RAN3
5431*                        is to be reset. (input).
5432*            ok       -  logical flag indicating all went well.
5433*                                                      (output).
5434*
5435*      COMMON VARIABLES:
5436*            uses:  a0, b0, c0, d0, lambda, DoSymDump, no_trials,
5437*                   max_angle, PI, PI2, RAD2DEG, cell_gamma, check_sym
5438*                   h_bnd, k_bnd, l_bnd, tolerance, tiny_inty
5439*
5440*        modifies:  max_var
5441*
5442*      TST_MIR returns logical .true. if the diffraction intensities
5443*      have the mirror symmetry about the plane requested.
5444* ______________________________________________________________________
5445*
5446      logical function TST_MIR(mir_sym, idum, ok)
5447      include 'DIFFaX.par'
5448      include 'DIFFaX.inc'
5449*
5450      integer*4 mir_sym, idum
5451      logical ok
5452*
5453      logical is_good, match, eq_sides
5454      integer*4 i, h, k, h_tmp, k_tmp
5455      logical cell90, cell120
5456      real*8 RAN3, PNTINT, S, ANGLE
5457      real*8 tiny, l
5458      real*8 i_avg, tol, i1, i2, variance, rel_var
5459      parameter (tiny = FIVE * eps4)
5460*
5461* external functions
5462      external RAN3, PNTINT
5463*
5464* statement functions
5465* S is the value of 1/d**2 at hkl
5466      S(h,k,l) = h*h*a0 + k*k*b0 + l*l*c0 + h*k*d0
5467* ANGLE is the Bragg angle (in radians) of the h,k,l plane
5468      ANGLE(h,k,l) = asin(HALF * lambda * sqrt(S(h,k,l)))
5469*
5470      cell90  = abs(cell_gamma - HALF*PI) .lt. HALF*PI*tiny
5471      cell120 = abs(cell_gamma - PI2/THREE) .lt. PI2*tiny/THREE
5472      eq_sides = abs(cell_a - cell_b) .le. HALF*eps6*(cell_a + cell_b)
5473      TST_MIR = .false.
5474      is_good = .false.
5475      if(mir_sym.lt.1 .or. mir_sym.gt.3) goto 900
5476      if(DoSymDump) then
5477        write(sy,200)
5478        if(.not.cell90 .and. .not.cell120) then
5479          write(sy,230) 'cell angle = ', cell_gamma * RAD2DEG,
5480     |      ' degrees. NO HORIZONTAL MIRRORS ARE LIKELY.'
5481          goto 900
5482        endif
5483      endif
5484*
5485      if(DoSymDump) then
5486        if(mir_sym.eq.1) then
5487          write(sy,199) 'Testing for mirror about the h-l plane'
5488        else if(mir_sym.eq.2) then
5489          write(sy,199) 'Testing for mirror about the k-l plane'
5490        else if(mir_sym.eq.3) then
5491          write(sy,199) 'Testing for mirror about the h=k,l plane'
5492        endif
5493        write(sy,200)
5494        write(sy,210)
5495      endif
5496      do 10 i = 1, no_trials
5497* get usable h,k >= 0
5498   20   if(mir_sym.eq.1) then
5499          h_tmp = int( dble(h_bnd + 1) * RAN3(idum) )
5500   30     k_tmp = int( dble(k_bnd + 1) * RAN3(idum) )
5501          if(k_tmp.eq.0) goto 30
5502        else if(mir_sym.eq.2) then
5503          k_tmp = int( dble(k_bnd + 1) * RAN3(idum) )
5504   40     h_tmp = int( dble(h_bnd + 1) * RAN3(idum) )
5505          if(h_tmp.eq.0) goto 40
5506        else if(mir_sym.eq.3) then
5507          k_tmp = int( dble(k_bnd + 1) * RAN3(idum) )
5508   45     h_tmp = int( dble(h_bnd + 1) * RAN3(idum) )
5509          if(h_tmp.eq.k_tmp) goto 45
5510        endif
5511* get usable l > 0
5512   50   l = l_bnd * RAN3(idum)
5513        if(abs(l).le.eps2) goto 50
5514* make sure we are not sampling at too high an angle
5515        if(TWO * ANGLE(h_tmp, k_tmp, l) .gt. max_angle) goto 20
5516* I(h,k,l)
5517        h = h_tmp
5518        k = k_tmp
5519        i1 = PNTINT(h, k, l, ok)
5520        if(.not.ok) goto 999
5521        if(DoSymDump) write(sy,220) h, k, l, i1
5522* mirror on h-l plane
5523        if(mir_sym.eq.1) then
5524* is the cell angle equal to 90 degrees
5525          if(cell90) then
5526* I(h,-k,l), rectangular cell
5527            h =  h_tmp
5528            k = -k_tmp
5529            i2 = PNTINT(h, k, l, ok)
5530            if(.not.ok) goto 999
5531            if(DoSymDump) write(sy,220) h, k, l, i2
5532          else if(cell120) then
5533* I(h+k,-k,l), hexagonal cell
5534            h =  h_tmp + k_tmp
5535            k = -k_tmp
5536            i2 = PNTINT(h, k, l, ok)
5537            if(.not.ok) goto 999
5538            if(DoSymDump) write(sy,220) h, k, l, i2
5539          endif
5540* else mirror on k-l plane, mir = 2
5541        else if(mir_sym.eq.2) then
5542* is the cell angle equal to 90 degrees
5543          if(cell90) then
5544* I(-h,k,l), rectangular cell
5545            h = -h_tmp
5546            k =  k_tmp
5547            i2 = PNTINT(h, k, l, ok)
5548            if(.not.ok) goto 999
5549            if(DoSymDump) write(sy,220) h, k, l, i2
5550          else if(cell120) then
5551* I(-h,h+k,l), hexagonal cell
5552            h = -h_tmp
5553            k =  h_tmp + k_tmp
5554            i2 = PNTINT(h, k, l, ok)
5555            if(.not.ok) goto 999
5556            if(DoSymDump) write(sy,220) h, k, l, i2
5557          endif
5558* else mirror on hk-l plane, mir = 3
5559        else if(mir_sym.eq.3) then
5560* The following if block is redundant, and in special
5561* cases fails to print to the .sym file. mmjt 3/18/04
5562* is the cell square
5563*          if(cell90 .and. eq_sides) then
5564* I(-h,k,l), square cell
5565*            h = k_tmp
5566*            k = h_tmp
5567*            i2 = PNTINT(h, k, l, ok)
5568*            if(.not.ok) goto 999
5569*            if(DoSymDump) write(sy,220) h, k, l, i2
5570*          else if(cell120) then
5571* I(-h,h+k,l), hexagonal cell
5572*            h = k_tmp
5573*            k = h_tmp
5574*            i2 = PNTINT(h, k, l, ok)
5575*            if(.not.ok) goto 999
5576*            if(DoSymDump) write(sy,220) h, k, l, i2
5577*          endif
5578          h = k_tmp
5579          k = h_tmp
5580          i2 = PNTINT(h, k, l, ok)
5581          if(.not.ok) goto 999
5582          if(DoSymDump) write(sy,220) h, k, l, i2
5583        endif
5584* compare mirrored intensities
5585        i_avg = HALF * (i1 + i2)
5586        variance = HALF * (abs(i_avg-i1) + abs(i_avg-i2))
5587* Be careful intensities are not actually zero
5588        if(i_avg.lt.tiny_inty) then
5589          tol = tiny_inty
5590        else
5591          tol = i_avg * tolerance
5592          rel_var = variance / i_avg
5593          if(rel_var.gt.max_var) max_var = rel_var
5594        endif
5595        match = abs(i_avg-i1).lt.tol .and. abs(i_avg-i2).lt.tol
5596        is_good = (i.eq.1 .or. is_good) .and. match
5597        if(DoSymDump) then
5598          write(sy,270)
5599          write(sy,240) i_avg
5600          write(sy,360) variance, HUNDRED * variance / i_avg
5601*          if(.not.check_sym) then
5602            write(sy,260) tol
5603            if(match) then
5604              if(mir_sym.eq.1) then
5605                write(sy,250)
5606              else if(mir_sym.eq.2) then
5607                write(sy,280)
5608              else if(mir_sym.eq.3) then
5609                write(sy,285)
5610              endif
5611            else
5612              if(mir_sym.eq.1) then
5613                write(sy,290)
5614              else if(mir_sym.eq.2) then
5615                write(sy,300)
5616              else if(mir_sym.eq.3) then
5617                write(sy,305)
5618              endif
5619            endif
5620*          endif
5621          write(sy,200)
5622        endif
5623   10 continue
5624      TST_MIR = is_good
5625*
5626      if(DoSymDump) then
5627*        if(.not.check_sym) then
5628          if(mir_sym.eq.1) then
5629            if(is_good) then
5630              write(sy,310)
5631            else
5632              write(sy,320)
5633            endif
5634          else if(mir_sym.eq.2) then
5635            if(is_good) then
5636              write(sy,330)
5637            else
5638              write(sy,340)
5639            endif
5640          else if(mir_sym.eq.3) then
5641            if(is_good) then
5642              write(sy,345)
5643            else
5644              write(sy,346)
5645            endif
5646          endif
5647*        endif
5648        write(sy,200)
5649      endif
5650*
5651  900 return
5652  999 write(op,199) 'ERROR in intensity calculation in TST_MIR'
5653      write(op,350) '   at h,k,l = ', h,',', k,',', l
5654      return
5655  199 format(1x, a)
5656  200 format(' ')
5657  210 format(1x, '  h', 5x, 'k', 7x, 'l', 20x, 'Intensity')
5658  220 format(1x, i3, 3x, i3, 2x, f9.4, 5x, f22.6)
5659  230 format(1x, a, f6.2, a)
5660  240 format(6x, 'Average Intensity = ', f22.6)
5661  250 format(1x, 'Intensities are consistent with an h-l mirror plane')
5662  260 format(1x, 'Intensity tolerance = +/-', f22.6)
5663  270 format(26x, '----------------------')
5664  280 format(1x, 'Intensities are consistent with a k-l mirror plane')
5665  285 format(1x, 'Intensities are consistent with a h=k,l mirror plane')
5666  290 format(1x, 'Intensities not consistent with an h-l mirror plane')
5667  300 format(1x, 'Intensities not consistent with a k-l mirror plane')
5668  305 format(1x, 'Intensities not consistent with a h=k,l mirror plane')
5669  310 format(1x, 'THERE IS A MIRROR ABOUT THE H-L PLANE')
5670  320 format(1x, 'THERE IS NO MIRROR ABOUT THE H-L PLANE')
5671  330 format(1x, 'THERE IS A MIRROR ABOUT THE K-L PLANE')
5672  340 format(1x, 'THERE IS NO MIRROR ABOUT THE K-L PLANE')
5673  345 format(1x, 'THERE IS A MIRROR ABOUT THE H=K,L PLANE')
5674  346 format(1x, 'THERE IS NO MIRROR ABOUT THE H=K,L PLANE')
5675  350 format(1x, a, i3, a, i3, a, f7.2)
5676  360 format(1x,'  Average variation = +/-', f22.6,'  (+/-',g9.2,'%)')
5677      end
5678*
5679* ______________________________________________________________________
5680* Title: TST_ROT
5681* Author: MMJT
5682* Date: 30 July 1989; 22 Feb 1995
5683* Identifies the rotational symmetry of the diffraction pattern.
5684* The rotational symmetry to test for is passed as 'rot_sym'.
5685* The values accepted for 'rot_sym' are 2, 3 and 4. If the diffraction
5686* intensity has the symmetry requested, TST_ROT returns '.true.'. If
5687* the pattern does not have the requested symmetry, or if an illegal
5688* value was passed (i.e. rot_sym = 5), TST_ROT returns '.false.'.
5689* NOTE. A 6-fold axis is found by calling TST_ROT twice: once for
5690* rot_sym = 2 and again for rot_sym = 3.
5691*
5692*      ARGUMENTS:
5693*            rot_sym  -  Rotational symmetry to test for.
5694*                        Accepts the values 2, 3 or 4. (input).
5695*            idum     -  parameter used by RAN3. Is -ve if RAN3
5696*                        is to be reset. (input).
5697*            ok       -  logical flag indicating all went well.
5698*                                                      (output).
5699*
5700*      COMMON VARIABLES:
5701*            uses:  a0, b0, c0, d0, lambda, DoSymDump, no_trials,
5702*                   max_angle, check_sym, h_bnd, k_bnd, l_bnd
5703*                   tolerance
5704*
5705*        modifies:  max_var
5706*
5707*      TST_ROT returns logical .true. if the diffraction intensities
5708*      have the requested rotational symmetry.
5709* ______________________________________________________________________
5710*
5711      logical function TST_ROT(rot_sym, idum, ok)
5712      include 'DIFFaX.par'
5713      include 'DIFFaX.inc'
5714*
5715      integer*4 rot_sym, idum
5716      logical ok
5717*
5718      logical is_good, match
5719      integer*4 i, h, k, h_tmp, k_tmp
5720      real*8 RAN3, PNTINT, S, ANGLE
5721      real*8 l, i_avg, tol
5722      real*8 i1, i2, i3, i4, variance, rel_var
5723*
5724* external functions
5725      external RAN3, PNTINT
5726*
5727* statement functions
5728* S is the value of 1/d**2 at hkl
5729      S(h,k,l) = h*h*a0 + k*k*b0 + l*l*c0 + h*k*d0
5730      ANGLE(h,k,l) = asin(HALF * lambda * sqrt(S(h,k,l)))
5731*
5732      TST_ROT = .false.
5733      is_good = .false.
5734* Is rot valid?
5735      if(rot_sym.lt.2 .or. rot_sym.gt.4) goto 900
5736* Now test for rotational symmetry.
5737* 2-fold and 4-fold
5738* Avoid both h, k = 0. Also, avoid l = 0, since Friedel's law will
5739* create a pseudo 2-fold.
5740      if(rot_sym.eq.2 .or. rot_sym.eq.4) then
5741        if(DoSymDump) then
5742          write(sy,210)
5743          write(sy,330) 'Testing for ', rot_sym, '-fold axis'
5744          write(sy,210)
5745        endif
5746        do 10 i = 1, no_trials
5747          if(DoSymDump) write(sy,220)
5748* get usable h,k >= 0
5749   20     h_tmp = int( dble(h_bnd + 1) * RAN3(idum) )
5750          k_tmp = int( dble(k_bnd + 1) * RAN3(idum) )
5751          if(h_tmp.eq.0 .and. k_tmp.eq.0) goto 20
5752* get usable l > 0
5753   30     l = l_bnd * RAN3(idum)
5754* keep l off the l = 0 plane, else we might confuse the inversion
5755* with a 2-fold
5756          if(abs(l).le.eps2) goto 30
5757* make sure we are not sampling at too high an angle
5758          if(TWO * ANGLE(h_tmp, k_tmp, l) .gt. max_angle) goto 20
5759* I(h,k,l)
5760          h = h_tmp
5761          k = k_tmp
5762          i1 = PNTINT(h, k, l, ok)
5763          if(.not.ok) goto 999
5764          if(DoSymDump) write(sy,230) h, k, l, i1
5765* I(-h,-k,l)
5766          h = -h_tmp
5767          k = -k_tmp
5768          i2 = PNTINT(h, k, l, ok)
5769          if(.not.ok) goto 999
5770          if(DoSymDump) write(sy,230) h, k, l, i2
5771* compare 2-fold intensities
5772          if(rot_sym.eq.2) then
5773            i_avg = HALF * (i1 + i2)
5774            variance = HALF * (abs(i_avg-i1) + abs(i_avg-i2))
5775* Be careful intensities are not actually zero
5776            if(i_avg.lt.tiny_inty) then
5777              tol = tiny_inty
5778            else
5779              tol = i_avg * tolerance
5780              rel_var = variance / i_avg
5781              if(rel_var.gt.max_var) max_var = rel_var
5782            endif
5783            match = abs(i_avg-i1).lt.tol .and. abs(i_avg-i2).lt.tol
5784            is_good = (i.eq.1 .or. is_good) .and. match
5785          else
5786* I(-k,h,l)
5787            h = -k_tmp
5788            k =  h_tmp
5789            i3 = PNTINT(h, k, l, ok)
5790            if(.not.ok) goto 999
5791            if(DoSymDump) write(sy,230) h, k, l, i3
5792* I(k,-h,l)
5793            h =  k_tmp
5794            k = -h_tmp
5795            i4 = PNTINT(h, k, l, ok)
5796            if(.not.ok) goto 999
5797            if(DoSymDump) write(sy,230) h, k, l, i4
5798* compare 4-fold intensities
5799            i_avg = QUARTER * (i1 + i2 + i3 + i4)
5800            variance = QUARTER * (abs(i_avg-i1) + abs(i_avg-i2) +
5801     |                           abs(i_avg-i3) + abs(i_avg-i4))
5802* Be careful intensities are not actually zero
5803            if(i_avg.lt.tiny_inty) then
5804              tol = tiny_inty
5805            else
5806              tol = i_avg * tolerance
5807              rel_var = variance / i_avg
5808              if(rel_var.gt.max_var) max_var = rel_var
5809            endif
5810            match = abs(i_avg-i1).lt.tol .and. abs(i_avg-i2).lt.tol
5811     |        .and. abs(i_avg-i3).lt.tol .and. abs(i_avg-i4).lt.tol
5812            is_good = (i.eq.1.or.is_good) .and. match
5813          endif
5814*
5815          if(DoSymDump) then
5816            write(sy,240)
5817            write(sy,250) i_avg
5818            write(sy,260) variance, HUNDRED * variance / i_avg
5819*            if(.not.check_sym) then
5820              write(sy,270) tol
5821              if(match) then
5822                write(sy,280) rot_sym
5823              else
5824                write(sy,290) rot_sym
5825              endif
5826*            endif
5827            write(sy,210)
5828          endif
5829   10   continue
5830        TST_ROT = is_good
5831        goto 900
5832      endif
5833* 3-fold
5834* Avoid both h, k = 0.
5835      if(rot_sym.eq.3) then
5836        if(DoSymDump) then
5837          write(sy,200) rot_sym
5838          write(sy,210)
5839          write(sy,220)
5840        endif
5841        do 40 i = 1, no_trials
5842* get usable h,k >= 0
5843   50     h_tmp = int( dble(h_bnd + 1) * RAN3(idum) )
5844          k_tmp = int( dble(k_bnd + 1) * RAN3(idum) )
5845          if(h_tmp.eq.0 .and. k_tmp.eq.0) goto 50
5846* get l (l=0 is allowed)
5847          l = l_bnd * RAN3(idum)
5848* make sure we are not sampling at too high an angle
5849          if(TWO * ANGLE(h_tmp, k_tmp, l) .gt. max_angle) goto 50
5850* I(h,k,l)
5851          h = h_tmp
5852          k = k_tmp
5853          i1 = PNTINT(h, k, l, ok)
5854          if(.not.ok) goto 999
5855          if(DoSymDump) write(sy,230) h, k, l, i1
5856* I(-h-k,h,l)
5857          h = -(h_tmp + k_tmp)
5858          k = h_tmp
5859          i2 = PNTINT(h, k, l, ok)
5860          if(.not.ok) goto 999
5861          if(DoSymDump) write(sy,230) h, k, l, i2
5862* I(k,-h-k,l)
5863          h = k_tmp
5864          k = -(h_tmp + k_tmp)
5865          i3 = PNTINT(h, k, l, ok)
5866          if(DoSymDump) write(sy,230) h, k, l, i3
5867* compare intensities
5868          i_avg = (i1 + i2 + i3) / THREE
5869          variance =
5870     |      (abs(i_avg-i1) + abs(i_avg-i2) + abs(i_avg-i3)) / THREE
5871* Be careful intensities are not actually zero
5872          if(i_avg.lt.tiny_inty) then
5873            tol = tiny_inty
5874          else
5875            tol = i_avg * tolerance
5876            rel_var = variance / i_avg
5877            if(rel_var.gt.max_var) max_var = rel_var
5878          endif
5879          match = abs(i_avg-i1).lt.tol .and. abs(i_avg-i2).lt.tol
5880     |      .and. abs(i_avg-i3).lt.tol
5881          is_good = (i.eq.1 .or. is_good) .and. match
5882          if(DoSymDump) then
5883            write(sy,240)
5884            write(sy,250) i_avg
5885            write(sy,260) variance, HUNDRED * variance / i_avg
5886*            if(.not.check_sym) then
5887              write(sy,270) tol
5888              if(match) then
5889                write(sy,280) rot_sym
5890              else
5891                write(sy,290) rot_sym
5892              endif
5893*            endif
5894            write(sy,210)
5895          endif
5896   40   continue
5897        TST_ROT = is_good
5898      endif
5899*
5900  900 if(DoSymDump) then
5901*        if(.not.check_sym) then
5902          if(is_good) then
5903            write(sy,300) rot_sym
5904          else
5905            write(sy,310) rot_sym
5906          endif
5907*        endif
5908        write(sy,210)
5909      endif
5910      return
5911*
5912  999 write(op,400) 'ERROR in intensity calculation in TST_ROT'
5913      write(op,320) '   at h,k,l = ', h,',', k,',', l
5914      return
5915  200 format(1x, 'Testing for a ', i1, '-fold axis')
5916  210 format(' ')
5917  220 format(1x, '  h', 5x, 'k', 7x, 'l', 20x, 'Intensity')
5918  230 format(1x, i3, 3x, i3, 2x, f9.4, 5x, f22.6)
5919  240 format(26x, '----------------------')
5920  250 format(6x, 'Average Intensity = ', f22.6)
5921  260 format(1x, '  Average variation = +/-',f22.6,'  (+/-',g9.2,'%)')
5922  270 format(1x, 'Intensity tolerance = +/-', f22.6)
5923  280 format(1x, 'Intensities are consistent with a ', i1, '-fold')
5924  290 format(1x, 'Intensities are not consistent with a ', i1, '-fold')
5925  300 format(1x, 'INTENSITY DISTRIBUTION HAS A ', i1, '-FOLD AXIS')
5926  310 format(1x, 'INTENSITY DISTRIBUTION HAS NO ', i1, '-FOLD AXIS')
5927  320 format(1x, a, i3, a, i3, a, f7.2)
5928  330 format(1x, a, i1, a)
5929  400 format(1x, a)
5930      end
5931*
5932* ______________________________________________________________________
5933* Title: XYPHSE
5934* Author: MMJT
5935* Date: 8 June 1990
5936* Description:  This routine pre-calculates the h and k components of
5937* phases for each atom. This is called when l = 0, since h and k
5938* are held constant while l is varied along each streak.
5939*
5940*      ARGUMENTS:
5941*            h  -  reciprocal lattice vector h-component. (input).
5942*            k  -  reciprocal lattice vector k-component. (input).
5943*
5944*      COMMON VARIABLES:
5945*            uses:  n_actual, l_n_atoms, a_pos, l_actual, n_layers
5946*
5947*        modifies:  hx_ky
5948* ______________________________________________________________________
5949*
5950      subroutine XYPHSE(h, k)
5951      include 'DIFFaX.par'
5952      include 'DIFFaX.inc'
5953*
5954      integer*4 h, k
5955*
5956      integer*4 i, m
5957*
5958      do 10 m = 1, n_actual
5959        do 20 i = 1, l_n_atoms(m)
5960          hx_ky(i,m) = h*a_pos(1,i,m) + k*a_pos(2,i,m)
5961   20   continue
5962   10 continue
5963*
5964      return
5965      end
5966*
5967* ______________________________________________________________________
5968* Title: YRDSTK
5969* Author: MMJT
5970* Date: 12 Aug 1989
5971* Description: YRDSTK checks that y is a multiple of x.
5972*
5973*      ARGUMENTS:
5974*            x   -  reference value. (input).
5975*            y   -  value to be tested. (input)
5976*            ok  -  logical flag set to .false. if y is zero. (output).
5977*
5978*      YRDSTK returns logical .true. if y is a multiple of x.
5979* ______________________________________________________________________
5980*
5981      logical function YRDSTK(x, y, ok)
5982      include 'DIFFaX.par'
5983*     save
5984*
5985      real*8 x, y
5986      logical ok
5987*
5988      real*8 tmp
5989*
5990      YRDSTK = .false.
5991      if(y.eq.ZERO) goto 999
5992      if(x.eq.ZERO) then
5993* This is the first visit to YRDSTK. Set things up.
5994        x = y
5995        YRDSTK = .true.
5996      else
5997        tmp = y / x
5998        if(abs(nint(tmp)-tmp) .le. eps3*tmp) YRDSTK = .true.
5999      endif
6000*
6001      return
6002  999 ok = .false.
6003      return
6004      end
6005*
Note: See TracBrowser for help on using the repository browser.