Changeset 2205


Ignore:
Timestamp:
Apr 8, 2016 3:10:04 PM (7 years ago)
Author:
vondreele
Message:

implement sequential stacking fault calculations with plotting
fix PlotXY

Location:
trunk
Files:
22 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIphsGUI.py

    r2204 r2205  
    24122412            'viewPoint':[[0.,0.,0.],[]],}
    24132413        Indx = {}
     2414       
     2415        def SetCell(laue,cell):
     2416            if laue in ['-3','-3m','6/m','6/mmm','4/m','4/mmm']:                   
     2417                cell[4] = cell[5] = 90.
     2418                cell[6] = 120.
     2419                if laue in ['4/m','4/mmm']:
     2420                    cell[6] = 90.
     2421                if ObjId == 0:
     2422                    cell[1] = cell[2] = value
     2423                    Obj.SetValue("%.5f"%(cell[1]))
     2424                else:
     2425                    cell[3] = value
     2426                    Obj.SetValue("%.5f"%(cell[3]))
     2427            elif laue in ['mmm']:
     2428                cell[ObjId+1] = value
     2429                cell[4] = cell[5] = cell[6] = 90.
     2430                Obj.SetValue("%.5f"%(cell[ObjId+1]))
     2431            elif laue in ['2/m','-1']:
     2432                cell[4] = cell[5] = 90.
     2433                if ObjId != 3:
     2434                    cell[ObjId+1] = value
     2435                    Obj.SetValue("%.5f"%(cell[ObjId+1]))
     2436                else:
     2437                    cell[6] = value
     2438                    Obj.SetValue("%.3f"%(cell[6]))
     2439            cell[7] = G2lat.calc_V(G2lat.cell2A(cell[1:7]))
    24142440
    24152441        def OnLaue(event):
    … …  
    24342460            G2plt.PlotXYZ(G2frame,XY,Layers['Sadp']['Img'].T,labelX=labels[:-1],
    24352461                labelY=labels[-1],newPlot=False,Title=Layers['Sadp']['Plane'])
     2462               
     2463        def OnSeqPlot(event):
     2464            seqPlot.SetValue(False)
     2465            resultXY,resultXY2 = Layers['seqResults']
     2466            pName = Layers['seqCodes'][0]
     2467            G2plt.PlotXY(G2frame,resultXY,XY2=resultXY2,labelX=r'$\mathsf{2\theta}$',
     2468            labelY='Intensity',newPlot=True,Title='Sequential simulations on '+pName,lines=False)
    24362469           
    24372470        def CellSizer():
    … …  
    28952928            sadpPlot.Bind(wx.EVT_CHECKBOX,OnSadpPlot)
    28962929            laueSizer.Add(sadpPlot,0,WACV)
     2930        if 'seqResults' in Layers:
     2931            seqPlot = wx.CheckBox(layerData,label=' Plot sequential result?')
     2932            seqPlot.Bind(wx.EVT_CHECKBOX,OnSeqPlot)
     2933            laueSizer.Add(seqPlot,0,WACV)
    28972934        topSizer.Add(laueSizer,0,WACV)
    28982935        topSizer.Add(wx.StaticText(layerData,label=' Reference unit cell for all layers:'),0,WACV)
    … …  
    29793016                return           
    29803017            profile = G2frame.PatternTree.GetItemPyData(G2frame.PatternId)[1]
    2981             dlg.Destroy()       
    29823018            ctrls = '0\n0\n3\n'
    29833019            G2pwd.StackSim(data['Layers'],ctrls,scale,background,limits,inst,profile)
    … …  
    29873023            test2 = np.copy(profile[3])
    29883024            rat = (test1-test2)/test1
     3025            XY = np.vstack((profile[0],rat))
     3026            G2plt.PlotXY(G2frame,[XY,],XY2=[],labelX=r'$\mathsf{2\theta}$',
     3027                labelY='ratio',newPlot=True,Title='DIFFaX vs GSASII',lines=True)
    29893028#            GSASIIpath.IPyBreak()
    29903029            G2plt.PlotPatterns(G2frame,plotType='PWDR')
    … …  
    29933032            data['Layers']['Sadp']['Plane'] = simCodes[1]
    29943033            data['Layers']['Sadp']['Lmax'] = simCodes[2]
    2995 #            planeChoice = ['h0l','0kl','hhl','h-hl',]
    2996 #            lmaxChoice = [str(i+1) for i in range(6)]
    2997 #            ctrls = '0\n0\n4\n1\n%d\n1\n16\n1\n%d\n0\nend\n'%    \
    2998 #                (planeChoice.index(simCodes[1])+1,lmaxChoice.index(simCodes[2])+1)
    2999 #            G2pwd.StackSim(data['Layers'],ctrls)
     3034            planeChoice = ['h0l','0kl','hhl','h-hl',]
     3035            lmaxChoice = [str(i+1) for i in range(6)]
     3036            ctrls = '0\n0\n4\n1\n%d\n1\n16\n1\n%d\n0\nend\n'%    \
     3037                (planeChoice.index(simCodes[1])+1,lmaxChoice.index(simCodes[2])+1)
     3038            G2pwd.StackSim(data['Layers'],ctrls)
    30003039            G2pwd.CalcStackingSADP(data['Layers'])
    30013040        wx.CallAfter(UpdateLayerData)
    … …  
    30033042    def OnSeqSimulate(event):
    30043043       
     3044        cellSel = ['cellA','cellB','cellC','cellG']
     3045        transSel = ['TransP','TransX','TransY','TransZ']
    30053046        ctrls = ''
     3047        data['Layers']['seqResults'] = []
     3048        data['Layers']['seqCodes'] = []
    30063049        Parms = G2pwd.GetStackParms(data['Layers'])
    30073050        dlg = G2gd.DIFFaXcontrols(G2frame,ctrls,Parms)
    30083051        if dlg.ShowModal() == wx.ID_OK:
    30093052            simCodes = dlg.GetSelection()
    3010             data['Layers']['Multi'] = [simCodes[2:5]]
    30113053        else:
    30123054            return
    3013         print 'do sequence of simulations on...',parm,parmRange,parmStep
     3055        UseList = []
     3056        for item in data['Histograms']:
     3057            if 'PWDR' in item:
     3058                UseList.append(item)
     3059        if not UseList:
     3060            wx.MessageBox('No PWDR data for this phase to simulate',caption='Data error',style=wx.ICON_EXCLAMATION)
     3061            return
     3062        dlg = wx.SingleChoiceDialog(G2frame,'Data to simulate','Select',UseList)
     3063        if dlg.ShowModal() == wx.ID_OK:
     3064            sel = dlg.GetSelection()
     3065            HistName = UseList[sel]
     3066        else:
     3067            return
     3068        dlg.Destroy()
     3069        PWDR = data['Histograms'][HistName]
     3070        G2frame.PatternId = G2gd.GetPatternTreeItemId(G2frame,G2frame.root,HistName)
     3071        sample = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(
     3072            G2frame,G2frame.PatternId, 'Sample Parameters'))
     3073        scale = sample['Scale'][0]
     3074        background = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(
     3075            G2frame,G2frame.PatternId, 'Background'))       
     3076        limits = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(
     3077            G2frame,G2frame.PatternId, 'Limits'))[1]
     3078        inst = G2frame.PatternTree.GetItemPyData(
     3079            G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Instrument Parameters'))[0]
     3080        if 'T' in inst['Type'][0]:
     3081            wx.MessageBox("Can't simulate neutron TOF patterns yet",caption='Data error',style=wx.ICON_EXCLAMATION)
     3082            return           
     3083        profile = np.copy(G2frame.PatternTree.GetItemPyData(G2frame.PatternId)[1])
     3084        resultXY2 = []
     3085        resultXY = [np.vstack((profile[0],profile[1])),]
     3086        data['Layers']['selInst'] = simCodes[1]
     3087        data['Layers']['seqCodes'] = simCodes[2:]
     3088        Layers = copy.deepcopy(data['Layers'])
     3089        pName = simCodes[2]
     3090        BegFin = simCodes[3]
     3091        nSteps = simCodes[4]
     3092        laue = Layers['Laue']
     3093        dStep = (BegFin[1]-BegFin[0])/nSteps
     3094        vals = np.linspace(BegFin[0],BegFin[1],nSteps+1,True)
     3095        for val in vals:
     3096            if 'cell' in pName:
     3097                cellId = cellSel.index(pName)
     3098                cell = Layers['Cell']
     3099                cell[cellId+1] = val
     3100                if laue in ['-3','-3m','6/m','6/mmm','4/m','4/mmm']:                   
     3101                    cell[2] = cell[1]
     3102                cell[7] = G2lat.calc_V(G2lat.cell2A(cell[1:7]))
     3103                Layers['Cell'] = cell
     3104            elif 'Trans' in pName:
     3105                names = pName.split(';')
     3106                transId = transSel.index(names[0])
     3107                iY = int(names[1])
     3108                iX = int(names[2])
     3109                Trans = Layers['Transitions'][iY]
     3110                if not transId:     #i.e. probability
     3111                    osum = 1.-Trans[iX][0]
     3112                    nsum = 1.-val
     3113                    for i in range(len(Trans)):
     3114                        if i != iX:
     3115                            Trans[i][0] *= (nsum/osum)
     3116                Trans[iX][transId] = val
     3117            G2pwd.CalcStackingPWDR(Layers,scale,background,limits,inst,profile)
     3118            resultXY2.append([np.vstack((profile[0],profile[3])),][0])
     3119        data['Layers']['seqResults'] = [resultXY,resultXY2]
     3120        wx.CallAfter(UpdateLayerData)
    30143121       
    30153122################################################################################

  • trunk/GSASIIplot.py

    r2203 r2205  
    24682468################################################################################
    24692469           
    2470 def PlotXY(G2frame,XY,XY2=None,labelX=None,labelY=None,newPlot=False,Title=''):
     2470def PlotXY(G2frame,XY,XY2=None,labelX=None,labelY=None,newPlot=False,Title='',lines=False):
    24712471    '''simple plot of xy data, used for diagnostic purposes
    24722472    '''
    … …  
    25112511    colors=['b','g','r','c','m','k']
    25122512    for ixy,xy in enumerate(XY):
    2513         X,Y = xy
    2514         Plot.plot(X,Y,colors[ixy%6]+'+',picker=False)
     2513        X,Y = XY[ixy]
     2514        if lines:
     2515            Plot.plot(X,Y,colors[ixy%6],picker=False)
     2516        else:
     2517            Plot.plot(X,Y,colors[ixy%6]+'+',picker=False)
    25152518    if len(XY2):
    25162519        for ixy,xy in enumerate(XY2):
    2517             X,Y = xy
     2520            X,Y = XY2[ixy]
    25182521            Plot.plot(X,Y,colors[ixy%6],picker=False)
    25192522    if not newPlot:

  • trunk/fsource/DIFFaXsubs/DIFFaX.inc

    r2191 r2205  
    375375      character*4 a_name(MAX_A,MAX_L), atom_l(MAX_TA)
    376376      character*12 pnt_grp
    377       character*16 sfname, cfname
     377*      character*16 sfname, cfname
    378378*
    379379      logical one_B(MAX_L), one_occup(MAX_L), Bs_zero(MAX_L, MAX_L),
    … …  
    416416      common /chars1/ a_name, atom_l
    417417      common /chars2/ pnt_grp
    418       common /chars3/ sfname, cfname
     418*      common /chars3/ sfname, cfname
    419419*
    420420      common /logic1/ one_B, one_occup, Bs_zero, there

  • trunk/fsource/DIFFaXsubs/DIFFaXsubs.for

    r2197 r2205  
    233233  100 format(1x, a)
    234234      end
    235 **
    236 ** ______________________________________________________________________
    237 ** Title: ATOMS
    238 ** Authors: MMJT
    239 ** Date: 18 Feb 90
    240 ** Description: This routine lists the legal atom names accepted by
    241 ** DIFFaX. These names are taken from the file 'sfname'.
    242 **
    243 **      ARGUMENTS:
    244 **           No arguments are used.
    245 **
    246 **      COMMON VARIABLES:
    247 **            uses:  sfname
    248 **
    249 **        modifies:  No COMMON variables are modified.
    250 ** ______________________________________________________________________
    251 **
    252 *      subroutine ATOMS()
    253 *      include 'DIFFaX.par'
    254 *      include 'DIFFaX.inc'
    255 **
    256 *      integer*4 i, LENGTH
    257 *      character*4 atom_name
    258 *      character*80 atomline
    259 *      character*120 line
    260 **
    261 ** external function
    262 *      external LENGTH
    263 **
    264 ** sfname has been opened once already,
    265 ** no need for elaborate error checking
    266 *      open(unit = sf, file = sfname, status = 'old', err = 999)
    267 **
    268 *      write(op,200) 'Legal atom types are:'
    269 **
    270 ** skip first few lines which contain no data, until we come to Hydrogen
    271 *    1 read(sf,'(a)') line
    272 *        if(line(1:4).ne.'H   ') goto 1
    273 *      backspace(unit = sf, err = 500)
    274 **
    275 ** initialize atomline
    276 *      atomline = ' '
    277 ** write 10 atom names to a line
    278 *   10 i = 0
    279 *   20 read(sf, '(a)', end = 100) line
    280 *        atom_name = line(1:4)
    281 *        write(atomline(i*7+1:i*7+7),'(3a)') '''', atom_name, ''' '
    282 *        i = i + 1
    283 *        if(mod(i,10).eq.0) then
    284 *          write(op,210) atomline
    285 *          goto 10
    286 *        endif
    287 *        goto 20
    288 *  100 continue
    289 **
    290 *      close(unit = sf, err = 600)
    291 **
    292 *  999 return
    293 *  500 write(op,220) 'Unable to backspace scattering factor file ''',
    294 *     |    sfname(1:LENGTH(sfname)), '''.'
    295 *      return
    296 *  600 write(op,220) 'Unable to close scattering factor file ''',
    297 *     |    sfname(1:LENGTH(sfname)), '''.'
    298 *      return
    299 *  200 format(1x, a)
    300 *  210 format(2x, a)
    301 *  220 format(1x, 'ERROR in ATOMS: ', 3a)
    302 *      end
    303 **
    304235* ______________________________________________________________________
    305236* Title: BINPOW
    … …  
    671602  205 format(1x, 'Diffraction point symmetry is found to be ''',a,'''')
    672603  300 format(1x, 'Re-evaluating diffraction symmetry')
    673       end
    674 *
    675 * ______________________________________________________________________
    676 * Title: CHOICE
    677 * Author: MWD
    678 * Date: 18 Aug 1988 (modified by MMJT 6 August 1989)
    679 * Description: This function compares a string with an array of keys,
    680 * and returns the index of the key that matches the first characters
    681 * of the string. CHOICE returns -1 if an error is found.
    682 *
    683 *      ARGUMENTS:
    684 *            flag  -  Character string to be matched. (input).
    685 *            list  -  Array of character strings. (input).
    686 *            n     -  Length of the array 'list'. (input).
    687 * ______________________________________________________________________
    688 *
    689       integer*4 function CHOICE(flag, list, n)
    690       implicit none
    691 *
    692       integer*4 n
    693 *      character*(*) flag, list(n)
    694       character*(*) flag
    695       character*80 list(n)
    696 *
    697       integer*4 i, j1, j2, LENGTH
    698 *
    699 * external function
    700       external LENGTH
    701 *
    702       i = 1
    703       j1 = LENGTH(flag)
    704 *
    705    10 j2 = LENGTH(list(i))
    706 * see if the string contained in list(i) is identical to that in flag
    707       if(j1.eq.j2 .and. index(flag, list(i)(1:j2)).eq.1) goto 20
    708         i = i + 1
    709         if(i.le.n) goto 10
    710 *
    711    20 if(i.gt.n) i = -1
    712       CHOICE = i
    713 *
    714       return
    715604      end
    716605*
    … …  
    12971186*
    12981187* ______________________________________________________________________
    1299 * Title: CNTARG
    1300 * Author: MMJT
    1301 * Date: 24 Feb 1990
    1302 * Description: Counts the number of arguments in the character
    1303 * string 'line', and returns the answer in CNTARG. Legal separators
    1304 * are, blanks, tabs, and commas.
    1305 *
    1306 *      ARGUMENTS:
    1307 *            line  -  Line of characters to be parsed. (input).
    1308 *
    1309 *      CNTARG returns the number of integer arguments in the line.
    1310 *      Returns -1 if an error occurred.
    1311 * ______________________________________________________________________
    1312 *
    1313       integer*4 function CNTARG(line)
    1314       implicit none
    1315 *
    1316       character*(*) line
    1317 *
    1318       logical in_arg
    1319       integer*4 blank, tab, comma, i, j, arg_cnt, lin_len
    1320       parameter (tab=9, blank=32, comma=44)
    1321 *
    1322       in_arg = .false.
    1323       arg_cnt = 0
    1324       CNTARG = -1
    1325 *
    1326       lin_len = len(line)
    1327       do 10 i = 1, lin_len
    1328         j = ichar(line(i:i))
    1329         if(j.eq.tab .or. j.eq.blank .or. j.eq.comma) then
    1330           in_arg = .false.
    1331         else
    1332           if(.not.in_arg) then
    1333             in_arg = .true.
    1334             arg_cnt = arg_cnt + 1
    1335           endif
    1336         endif
    1337    10 continue
    1338 *
    1339 * There cannot be more than (lin_len+1)/2 arguments
    1340       if(arg_cnt.le.((lin_len+1)/2)) CNTARG = arg_cnt
    1341 *
    1342       return
    1343       end
    1344 *
    1345 * ______________________________________________________________________
    13461188* Title: DETUN
    13471189* Author: MMJT
    … …  
    13881230C  401 format(1x, g12.5)
    13891231      end
    1390 *
    1391 * ______________________________________________________________________
    1392 ** Title: DUMP
    1393 ** Author: MWD and MMJT
    1394 ** Date: 18 Aug 1988; 15 Mar 1995
    1395 ** Description: This subroutine prints out the data read in from
    1396 ** the data file. It reassures the user that the data was read in
    1397 ** correctly.
    1398 **
    1399 **      ARGUMENTS:
    1400 **            infile  -  The name of the input data file. (input).
    1401 **            ok      -  logical flag indicating all went well.
    1402 **                                                      (output).
    1403 **
    1404 **      COMMON VARIABLES:
    1405 **            uses:      CFile, GAUSS, LORENZ, NEUTRN, PI2, PS_VGT,
    1406 **                       RAD2DEG, SymGrpNo, X_RAY, a_B, a_name,
    1407 **                       a_number, a_occup, a_pos, a_type, blurring,
    1408 **                       cell_a, cell_b, cell_c, cell_gamma, cntrl
    1409 **                       d_theta, e_sf, FWHM, inf_thick, l_actual
    1410 **                       l_alpha, l_cnt, l_g, l_n_atoms, l_r, l_symmetry
    1411 **                       lambda, l_seq, n_actual, n_atoms, n_layers
    1412 **                       pnt_grp, pv_gamma, pv_u, pv_v, pv_w, r_B11
    1413 **                       r_B12, r_B22, r_B23, r_B31, r_B33, rad_type
    1414 **                       recrsv, rndm, th2_max, th2_min
    1415 **                       tolerance, xplcit
    1416 **
    1417 **        modifies:      no COMMON variables are modified
    1418 ** ______________________________________________________________________
    1419 **
    1420 *      subroutine DUMP(infile,ok)
    1421 *      include 'DIFFaX.par'
    1422 *      include 'DIFFaX.inc'
    1423 **
    1424 *      character*(*) infile
    1425 *      logical ok
    1426 **
    1427 *      real*8 scale, atom_cnt(MAX_TA), cum_atom_cnt(MAX_TA), norm
    1428 *      integer*4 i, i2, j, n, type, num_types, tot_types
    1429 *      integer*4 LENGTH, print_width
    1430 *      character*31 dmpfile
    1431 *      character*80 list(5)
    1432 **
    1433 ** external function
    1434 *      external LENGTH
    1435 ** external subroutine (Some compilers need them declared external)
    1436 **      external GETFNM
    1437 **
    1438 *      call GETFNM(infile, dmpfile, 'dmp', ok)
    1439 *      if(.not.ok) then
    1440 *        write(op,200) 'DUMP aborted'
    1441 *        goto 999
    1442 *      endif
    1443 *      if(dp.ne.op) open(unit = dp , file = dmpfile, status = 'new')
    1444 **
    1445 *      i2 = 0
    1446 *   99 write(op,200) 'Enter 1 for full atomic position dump'
    1447 *      read(cntrl,*,err=99,end=9999) i2
    1448 *      if(CFile) write(op,'(1x,i3)') i2
    1449 *      write(op,300) 'Writing data dump to file ''',
    1450 *     |               dmpfile(1:LENGTH(dmpfile)),'''. . .'
    1451 ** sundry details about layers
    1452 *      write(dp,125) 'Number of layers = ', n_layers
    1453 *      write(dp,125) 'Number of unique layers = ', n_actual
    1454 *      write(dp,125) 'Number of different atom types = ', n_atoms
    1455 ** cell dimensions
    1456 *      write(dp,140)
    1457 *     |    ' cell_a,      cell_b,      cell_c,      cell_gamma ',
    1458 *     |      cell_a, cell_b, cell_c, RAD2DEG * cell_gamma
    1459 ** in-plane layer widths
    1460 *      if(finite_width) then
    1461 *        if(Wa.lt.inf_width) then
    1462 *          write(dp,126) 'width along a', Wa
    1463 *          write(dp,128) Wa/cell_a
    1464 *        else
    1465 *          write(dp,127) 'a'
    1466 *        endif
    1467 *        if(Wb.lt.inf_width) then
    1468 *          write(dp,126) 'width along b', Wb
    1469 *          write(dp,128) Wb/cell_b
    1470 *        else
    1471 *          write(dp,127) 'b'
    1472 *        endif
    1473 *      else
    1474 *        write(dp,200)
    1475 *     |    'Layers are to be treated as having infinite lateral width.'
    1476 *      endif
    1477 ** radiation type
    1478 *      list(1) = 'X-RAY '
    1479 *      list(2) = 'NEUTRON '
    1480 *      list(3) = 'ELECTRON '
    1481 *      write(dp,100) 'radiation type = ', list(rad_type+1)
    1482 ** wavelength
    1483 *      write(dp,170) 'Wavelength lambda = ', lambda
    1484 ** symmetry
    1485 *      write(dp,100)
    1486 *     |'Diffraction point group symmetry specified = ',pnt_grp
    1487 *      if(SymGrpNo.eq.UNKNOWN) write(dp,175)
    1488 *     |'Tolerance on intensities in symmetry evaluation = +/-',
    1489 *     |                    tolerance * HUNDRED, ' %'
    1490 ** instrumental broadening to simulate
    1491 *      list(1) = 'NONE'
    1492 *      list(2) = 'GAUSSIAN'
    1493 *      list(3) = 'LORENTZIAN'
    1494 *      list(4) = 'PSEUDO-VOIGT'
    1495 *      i = blurring + 1
    1496 ** see if it's a pure Gaussian pseudo-Voigt
    1497 *      if(blurring.eq.PV_GSS) i = 2
    1498 ** see if it's a pure Lorentzian pseudo-Voigt
    1499 *      if(blurring.eq.PV_LRN) i = 3
    1500 *      write(dp,100)
    1501 *     |  'Instrumental broadening to be simulated = ', list(i)
    1502 *      if(blurring.eq.GAUSS .or. blurring.eq.LORENZ) write(dp,170)
    1503 *     |'Full width half-maximum = ', FWHM
    1504 *      if(blurring.eq.PS_VGT) then
    1505 *        write(dp,200)
    1506 *     |'Pseudo-Voigt parameters:   u,       v,       w,     gamma'
    1507 *        write(dp,180) pv_u, pv_v, pv_w, pv_gamma
    1508 *      endif
    1509 *      if(blurring.eq.PV_GSS) then
    1510 *        write(dp,200) 'Gaussian parameters:      u,       v,       w'
    1511 *        write(dp,185) pv_u, pv_v, pv_w
    1512 *      endif
    1513 *      if(blurring.eq.PV_LRN) then
    1514 *        write(dp,200) 'Lorentzian parameters:    u,       v,       w'
    1515 *        write(dp,185) pv_u, pv_v, pv_w
    1516 *      endif
    1517 ** are we to trim out the origin?
    1518 *      if(blurring.ne.NONE) then
    1519 *        if(trim_origin) then
    1520 *          write(dp,200) '  Intensity near origin will be ignored'
    1521 *        else
    1522 *          write(dp,200) '  Intensity near origin will be included'
    1523 *        endif
    1524 *      endif
    1525 **
    1526 ** equivalent layers
    1527 *      write(dp,200) ' '
    1528 *      do 10 i = 1, n_layers
    1529 *        write(dp,110) 'LAYER ',i,
    1530 *     |            ' is equivalent to fundamental LAYER ',l_actual(i)
    1531 *        write(dp,170) '   Existence probability = ', l_g(i)
    1532 *   10 continue
    1533 **
    1534 *        do 11 j = 1, MAX_TA
    1535 *          cum_atom_cnt(j) = ZERO
    1536 *   11   continue
    1537 **
    1538 ** layer structure
    1539 *      write(dp,200) ' '
    1540 *      do 30 i = 1, n_actual
    1541 *        write(dp,130) 'fundamental LAYER ',i
    1542 *        list(1) = 'NONE '
    1543 *        list(2) = 'CENTROSYMMETRIC '
    1544 *        write(dp,100) 'symmetry = ', list(l_symmetry(i)+1)
    1545 ** write out the layer composition
    1546 *        write(dp,200) 'LAYER composition is:'
    1547 *        do 14 j = 1, MAX_TA
    1548 *          atom_cnt(j) = ZERO
    1549 *   14   continue
    1550 *        num_types = ZERO
    1551 *        scale = ONE
    1552 *        if(l_symmetry(i).eq.CENTRO) scale = TWO
    1553 *        do 15 j = 1, l_n_atoms(i)
    1554 *          type = a_type(j,i)
    1555 *          if(type.gt.num_types) num_types = type
    1556 *          atom_cnt(type) = atom_cnt(type) + a_occup(j,i) * scale
    1557 *   15   continue
    1558 *        if(num_types.gt.tot_types) tot_types = num_types
    1559 ** accumulate weighted atom count
    1560 *        do 16 j = 1, MAX_TA
    1561 *          cum_atom_cnt(j) = cum_atom_cnt(j) + atom_cnt(j) * l_g(i)
    1562 *   16   continue
    1563 *        do 17 j = 1, num_types
    1564 *          write(dp,280) atom_l(j), ':', atom_cnt(j)
    1565 *   17   continue
    1566 *        write(dp,200) ' '
    1567 **
    1568 ** do we want all the details about each atom?
    1569 *        if(i2.ne.1) goto 30
    1570 ** yes, go for it.
    1571 *        do 20 j = 1, l_n_atoms(i)
    1572 *            write(dp,200) ' '
    1573 *          write(dp,120) 'ATOM number', a_number(j,i), ' in layer',i,
    1574 *     |     ' is ''', a_name(j,i),''''
    1575 *          write(dp,121)
    1576 *          write(dp,122) a_pos(1,j,i), a_pos(2,j,i), a_pos(3,j,i),
    1577 *     |      a_B(j,i), a_occup(j,i)
    1578 *          if(rad_type.eq.X_RAY) then
    1579 *            write(dp,140) 'X-ray scattering factor data Ai, Bi',
    1580 *     |                         (x_sf(n,a_type(j,i)),n=1,9)
    1581 *          else if(rad_type.eq.NEUTRN) then
    1582 *            write(dp,140) 'Neutron scattering factor data',
    1583 *     |                         n_sf(a_type(j,i))
    1584 *          else if(rad_type.eq.ELECTN) then
    1585 *            write(dp,140)'Electron scattering factor data Ai, Bi and Z',
    1586 *     |               (x_sf(n,a_type(j,i)),n=1,9), e_sf(a_type(j,i))
    1587 *          else
    1588 *            write(op,200) 'ERROR: Illegal radiation type in DUMP'
    1589 *          endif
    1590 *   20   continue
    1591 *        write(dp,200) ' '
    1592 *   30 continue
    1593 **
    1594 ** Print out average composition
    1595 *      if(.not.xplcit) then
    1596 *        norm = ZERO
    1597 *        do 25 j = 1, num_types
    1598 *          norm = norm + cum_atom_cnt(j)
    1599 *   25   continue
    1600 *        write(dp,200) ' '
    1601 *        write(dp,200) 'Average crystal composition is:'
    1602 *        do 26 j = 1, num_types
    1603 *          write(dp,281) atom_l(j), ':', cum_atom_cnt(j) / norm
    1604 *   26   continue
    1605 *        write(dp,200) ' '
    1606 *      endif
    1607 **
    1608 ** stacking details
    1609 *      write(dp,200) ' '
    1610 *      write(dp,200) 'STACKING'
    1611 *      if(recrsv) then
    1612 *        write(dp,210) 'RECURSIVELY'
    1613 *        if(inf_thick) then
    1614 *          write(dp,220) 'INFINITE'
    1615 *        else
    1616 *          write(dp,230) l_cnt
    1617 *        endif
    1618 *      else if(xplcit) then
    1619 *        if(rndm) then
    1620 *          write(dp,240)
    1621 *        else
    1622 *          write(dp,250)
    1623 *        endif
    1624 *        write(dp,260) 'Sequence for ', l_cnt, ' layers is:'
    1625 *        print_width = 25
    1626 *        j = 1
    1627 *        n = print_width
    1628 *   35   if(n.gt.l_cnt) n = l_cnt
    1629 *        write(dp,270) (l_seq(i), i = j, n)
    1630 *        if(n.lt.l_cnt) then
    1631 *          j = j + print_width
    1632 *          n = n + print_width
    1633 *          goto 35
    1634 *        endif
    1635 *      endif
    1636 **
    1637 ** stacking transition data
    1638 *      write(dp,200) ' '
    1639 *      write(dp,200) ' '
    1640 *      write(dp,200) 'TRANSITIONS'
    1641 *      write(dp,200)
    1642 *     |      'Layer stacking probabilities and stacking vectors.'
    1643 *      write(dp,200)
    1644 *     |'  i  j |   alpha_ij      Rx_ij        Ry_ij        Rz_ij'
    1645 *      write(dp,200)
    1646 *     |'  -----|-------------------------------------------------'
    1647 *      do 40 i = 1, n_layers
    1648 *        write(dp,200)
    1649 *     |'       |'
    1650 *        do 50 j = 1, n_layers
    1651 *          write(dp,150) i, j, l_alpha(j,i),
    1652 *     |                   l_r(1,j,i), l_r(2,j,i), l_r(3,j,i)
    1653 *   50   continue
    1654 *   40 continue
    1655 **
    1656 *      write(dp,200) ' '
    1657 *      write(dp,200)
    1658 *     |      'Anisotropic Layer stacking ''uncertainty'' factors.'
    1659 ** MMJT: 3/15/95. Ordering of B23 and B31 swapped. B31 -> C13
    1660 *      write(dp,200)
    1661 *     |'  i  j |     C11      C22      C33      C12      C13      C23'
    1662 *      write(dp,200)
    1663 *     |'  -----|-------------------------------------------------------'
    1664 *      do 60 i = 1, n_layers
    1665 *        write(dp,200)
    1666 *     |'       |'
    1667 *        do 70 j = 1, n_layers
    1668 ** MMJT: 3/15/95. Ordering of B23 and B31 swapped
    1669 *          write(dp,151) i, j, r_B11(j,i), r_B22(j,i), r_B33(j,i),
    1670 *     |                        r_B12(j,i), r_B31(j,i), r_B23(j,i)
    1671 *   70   continue
    1672 *   60 continue
    1673 **
    1674 *      if(dp.ne.op) close(unit = dp)
    1675 *  999 return
    1676 * 9999 ok = .false.
    1677 *      write(op,200) 'DUMP aborted'
    1678 *      return
    1679 *  100 format(1x, 2a)
    1680 *  110 format(1x, a, i4, a, i4)
    1681 *  120 format(1x, a, i3, a, i2, 3a)
    1682 *  121 format(4x,'x_rel',8x,'y_rel',8x,'z_rel',10x,'DW',10x,'Occ')
    1683 *  122 format(1x, 5g13.5)
    1684 *  125 format(1x, a40, i4)
    1685 *  126 format(1x, 'layer characteristic ', a, ' = ', f9.2, ' Angstroms')
    1686 *  127 format(1x, 'layer characteristic ', a, ' = INFINITY Angstroms')
    1687 *  128 format(1x, '   which is equivalent to ', f9.2, ' unit cells')
    1688 *  130 format(1x, a, i4)
    1689 *  140 format(1x, a, / 5g13.5, / 4g13.5, i4)
    1690 *  150 format(1x, 2i3, ' |', 4(1x, g12.5))
    1691 *  151 format(1x, 2i3, ' |', 6(1x, f8.3))
    1692 *  170 format(1x, a, g12.5)
    1693 *  175 format(1x, a, g12.5, a)
    1694 *  180 format(21x, 4(2x, f7.3))
    1695 *  185 format(21x, 3(2x, f7.3))
    1696 *  200 format(1x, a)
    1697 *  210 format(1x, 'Stacking is to be treated ', a, ' by DIFFaX.')
    1698 *  220 format(1x, 'Number of layers along the fault direction is ', a)
    1699 *  230 format(1x, 'There are ',i5 ,' layers along the fault direction')
    1700 *  240 format(1x, 'Sequencing is defined RANDOMLY by DIFFaX.')
    1701 *  250 format(1x, 'Sequencing is defined EXPLICITLY by the user.')
    1702 *  260 format(1x, a, i4, a)
    1703 *  270 format(1x, 30i3)
    1704 *  280 format(23x, 2a, f7.2)
    1705 *  281 format(23x, 2a, f7.4)
    1706 *  300 format(1x, 3a)
    1707 *      end
    17081232*
    17091233* ______________________________________________________________________
    … …  
    17671291*
    17681292* ______________________________________________________________________
    1769 ** Title: FNDCTL
    1770 ** Author: MMJT
    1771 ** Date: 6 June 1989
    1772 ** Description: Checks whether or not there is a control file, and
    1773 ** sets the i/o unit specifier 'cntrl' to 99.
    1774 **
    1775 **      ARGUMENTS:
    1776 **           No arguments are used. All data is in 'COMMON'.
    1777 **
    1778 **      COMMON VARIABLES:
    1779 **            uses:   cfname, CFile
    1780 **
    1781 **        modifies:   cntrl
    1782 **
    1783 **      FNDCTL returns logical .true. if a control file was found in the
    1784 **      current directory.
    1785 ** ______________________________________________________________________
    1786 **
    1787 *      logical function FNDCTL()
    1788 *      include 'DIFFaX.par'
    1789 *      include 'DIFFaX.inc'
    1790 **
    1791 *      integer*4 io_err, LENGTH
    1792 **
    1793 ** external function
    1794 *      external LENGTH
    1795 **
    1796 ** See if control file is present
    1797 *      inquire(file = cfname, exist = CFile)
    1798 ** get data from ControlFile
    1799 *      if(CFile) then
    1800 *        cntrl = 99
    1801 *        open(unit = cntrl, file = cfname, status = 'old',
    1802 *     |      err = 999, iostat = io_err)
    1803 *      else
    1804 ** get data from user entries to the screen
    1805 *        cntrl = ip
    1806 *      endif
    1807 *      FNDCTL = .true.
    1808 **
    1809 *      return
    1810 *  999 write(op,100)
    1811 *     |'A control file named ''', cfname(1:LENGTH(cfname)), ''''
    1812 *      write(op,101) 'exists, but there were problems opening it.'
    1813 *      write(op,102) 'iostat = ', io_err
    1814 *      FNDCTL = .false.
    1815 *      return
    1816 *  100 format(1x, 3a)
    1817 *  101 format(1x, a)
    1818 *  102 format(1x, a, i5)
    1819 *      end
    1820 **
    1821 ** ______________________________________________________________________
    1822 ** Title: GETFIL
    1823 ** Author: MMJT
    1824 ** Date: 5 June 1989
    1825 ** Description: Checks for the presence of the structure data file and
    1826 ** the structure factor data file in the current directory.
    1827 **
    1828 **      ARGUMENTS:
    1829 **            fname   -  The name of the input data file. (output).
    1830 **            ok      -  logical flag indicating all went well.
    1831 **                                                      (output).
    1832 **            ending  -  returns logical .true. if 'end' was read
    1833 **                       from the control file. (output).
    1834 **
    1835 **      COMMON VARIABLES:
    1836 **            uses:       cntrl, sfname, CFile
    1837 **
    1838 **        modifies:       no COMMON variables are modified
    1839 ** ______________________________________________________________________
    1840 **
    1841 *      subroutine GETFIL(fname, ok, ending)
    1842 *      include 'DIFFaX.par'
    1843 *      include 'DIFFaX.inc'
    1844 **
    1845 *      character*(*) fname
    1846 *      logical ok, ending
    1847 **
    1848 *      character*31 tmp_nam
    1849 *      character*200 line
    1850 *      character*80 list(1)
    1851 *      integer*4 unit_no, LENGTH, CHOICE
    1852 **
    1853 *      integer*4 i, io_err
    1854 **
    1855 ** external function
    1856 *      external LENGTH, CHOICE
    1857 ** external subroutines (Some compilers need them declared external)
    1858 **      external GETLNE, TOUPPR
    1859 **
    1860 *      unit_no = cntrl
    1861 ** get file name from ControlFile
    1862 *      write(op,104) 'Enter name of structure file.'
    1863 ** use GETLNE to strip off any leading blanks etc...
    1864 *      call GETLNE(unit_no, line, *993)
    1865 *      write(fname,'(a)') line(1:len(fname))
    1866 *      if(CFile) write(op,104) fname
    1867 ** check that the string does not say 'end'. First convert to upper case
    1868 ** we must be sure that the string 'end' is not part of a legitimate
    1869 ** file name. Use CHOICE for this.
    1870 *      tmp_nam = fname
    1871 *      call TOUPPR(fname)
    1872 *      list(1) = 'END '
    1873 *      i = CHOICE(fname, list, 1)
    1874 *      if(i.eq.1) then
    1875 *        ending = .true.
    1876 *        goto 999
    1877 *      endif
    1878 *      fname = tmp_nam
    1879 **
    1880 ** Open structure data file.
    1881 *      write(op,100) 'Looking for structure data file ''',
    1882 *     |            fname(1:LENGTH(fname)),''''
    1883 *      open(unit = df, file = fname, status = 'old',
    1884 *     |            err = 991, iostat = io_err)
    1885 *      write(op,100) 'Opening structure file ''',
    1886 *     |       fname(1:LENGTH(fname)),''''
    1887 **
    1888 **Open scattering factor data file.
    1889 *      write(op,100)
    1890 *     | 'Looking for scattering factor data file ''',
    1891 *     |            sfname(1:LENGTH(sfname)),''''
    1892 *      open(unit = sf, file = sfname, status = 'old',
    1893 *     |            err = 992, iostat = io_err)
    1894 *      write(op,100) 'Opening scattering factor data file ''',
    1895 *     |       sfname(1:LENGTH(sfname)),''''
    1896 **
    1897 *      goto 999
    1898 *  991 write(op,100)
    1899 *     |'Problems opening structure data file ''',
    1900 *     |            fname(1:LENGTH(fname)),''''
    1901 *      write(op,103) 'IOSTAT = ',io_err
    1902 *      ok = .false.
    1903 *      goto 999
    1904 *  992 write(op,100)
    1905 *     |'Problems opening scattering factor file ''',
    1906 *     |            sfname(1:LENGTH(sfname)),''''
    1907 *      write(op,103) 'IOSTAT = ',io_err
    1908 *      ok = .false.
    1909 *      goto 999
    1910 *  993 write(op,104) 'ERROR reading filename in GETFIL'
    1911 *      ok = .false.
    1912 *  999 return
    1913 *  100 format(1x, 3a)
    1914 *  103 format(1x, a, i5)
    1915 *  104 format(1x, a)
    1916 *      end
    1917 **
    1918 ** ______________________________________________________________________
    1919 ** Title: GETFNM
    1920 ** Author: MMJT
    1921 ** Date: 13 Oct 1988
    1922 ** Description: This subroutine makes sure that we do not write over
    1923 ** an existing file, and creates a suitable filename (name2). name2
    1924 ** equals 'name1.append' append is supplied by the calling routine.
    1925 ** If 'name1.append' is the name of an existing file, then append is
    1926 ** set to 'append#' where #=1,2.
    1927 ** CLIP is the maximum length of the appendage. Thus if append
    1928 ** = 'fred', and CLIP is 3, then the appendages become 'fre',
    1929 ** 'fr1', 'fr2', ...'f10', 'f11' ... up to 'f99'.
    1930 ** If CLIP is greater than 10, there is no restriction.
    1931 ** NOTE: It is up to the user to ensure that filename lengths are
    1932 ** compatible with the operating system that DIFFaX is being run on.
    1933 ** 'CLIP' and 'MAX_NAM' are set in 'DIFFaX.par'
    1934 **
    1935 **      ARGUMENTS:
    1936 **            name1   -  The name of the input data file. (input).
    1937 **            name2   -  A derivative filename. (output).
    1938 **            append  -  A token to be appended to name2. (input).
    1939 **            ok      -  logical flag signalling all went well.
    1940 **                                                      (output).
    1941 ** ______________________________________________________________________
    1942 **
    1943 *      subroutine GETFNM(name1,name2,append,ok)
    1944 *      include 'DIFFaX.par'
    1945 **     save
    1946 **
    1947 *      character*(*) name1, name2, append
    1948 *      logical ok
    1949 **
    1950 *      logical fexist
    1951 *      character*31 appnd2
    1952 *      integer*4 applen, i, myclip
    1953 *      integer*4 LENGTH
    1954 **
    1955 ** external function
    1956 *      external LENGTH
    1957 ** external subroutine (Some compilers need them declared external)
    1958 **      external NAMER
    1959 **
    1960 *      ok = .true.
    1961 *      fexist = .true.
    1962 *      appnd2 = append
    1963 *      applen = LENGTH(append)
    1964 *      myclip = CLIP
    1965 *      if(myclip.le.0) then
    1966 *        write(op,302) 'Illegal CLIP value ', CLIP
    1967 *        write(op,300) 'Default CLIP value of 3 is being assumed'
    1968 *        myclip = 3
    1969 *      endif
    1970 **
    1971 *      i = 0
    1972 *   10 call NAMER(name1,name2,appnd2)
    1973 *      inquire(file = name2, exist = fexist)
    1974 *      if(fexist) then
    1975 *        i = i + 1
    1976 *        if(i.lt.10) then
    1977 *          if(myclip.gt.10) then
    1978 *            write(appnd2(applen+1:applen+1),200) i
    1979 *          else
    1980 *            write(appnd2(myclip:myclip),200) i
    1981 *          endif
    1982 *        else if(i.ge.10 .and. i.lt.100) then
    1983 *          if(myclip.gt.10) then
    1984 *            write(appnd2(applen+1:applen+2),201) i
    1985 *          else if(myclip.gt.1) then
    1986 *            write(appnd2(myclip-1:myclip),201) i
    1987 *          else
    1988 *            ok = .false.
    1989 *            goto 999
    1990 *          endif
    1991 *        else if(i.ge.100 .and. i.lt.1000) then
    1992 *          if(myclip.gt.10) then
    1993 *            write(appnd2(applen+1:applen+3),202) i
    1994 *          else if(myclip.gt.2) then
    1995 *            write(appnd2(myclip-2:myclip),202) i
    1996 *          else
    1997 *            ok = .false.
    1998 *            goto 999
    1999 *          endif
    2000 *        else
    2001 *          ok = .false.
    2002 *          goto 999
    2003 *        endif
    2004 *        goto 10
    2005 *      endif
    2006 **
    2007 *      return
    2008 *  999 write(op,300) 'No new files can be created. . .'
    2009 *      write(op,300) '      . . . some old files should be deleted.'
    2010 *      write(op,301)
    2011 *     |      'Last file created was ''',name2(1:LENGTH(name2)),''''
    2012 *      return
    2013 *  200 format(i1)
    2014 *  201 format(i2)
    2015 *  202 format(i3)
    2016 *  300 format(1x, a)
    2017 *  301 format(1x, 3a)
    2018 *  302 format(1x, a, i5)
    2019 *      end
    2020 **
    2021 ** ______________________________________________________________________
    20221293* Title: GETLAY
    20231294* Author: MMJT
    … …  
    21051376*
    21061377* ______________________________________________________________________
    2107 ** Title: GETLNE
    2108 ** Author: MWD  & MMJT
    2109 ** Date: MWD: Original version 18 Aug 1988; MMJT: This version 9 Apr 1992
    2110 ** Description: This function returns a non-blank, left-justified line,
    2111 ** stripped of any comments, from unit 'unit_no'. Any error generates
    2112 ** a return to the passed linenumber, presumably an error-handling
    2113 ** routine in the calling procedure. Multiple comments, and nested
    2114 ** braces on one line can be handled.
    2115 **
    2116 **      ARGUMENTS:
    2117 **            unit_no  -  Logical unit number that data is to
    2118 **                        be read from. (input).
    2119 **            line     -  Line of data read as a character string.
    2120 **                                                          (output).
    2121 **            *        -  The label of a line in the calling routine
    2122 **                        to return to if an error occurs. (input).
    2123 ** ______________________________________________________________________
    2124 **
    2125 *      subroutine GETLNE(unit_no, line, *)
    2126 *      include 'DIFFaX.par'
    2127 **     save
    2128 **
    2129 *      integer*4 unit_no
    2130 *      character*(*) line
    2131 **
    2132 *      character*1 tab, blank
    2133 ** nominal line length is 200. Leave generous room for overflow.
    2134 *      character*250 rawline
    2135 *      integer*4 i, j, n, offset, lin_len
    2136 *      integer*4 l_br, r_br, first_left, last_right
    2137 **
    2138 *      tab = char(9)
    2139 *      blank = char(32)
    2140 *      lin_len = len(line)
    2141 **
    2142 *    1 read(unit_no, '(a)', err = 200, end = 200) rawline
    2143 *      line = rawline
    2144 **
    2145 ** First check if user has gone beyond lin_len characters
    2146 *      do 2 i = len(rawline), lin_len + 1, -1
    2147 *        if(rawline(i:i).ne.blank) goto 998
    2148 *    2 continue
    2149 **
    2150 ** strip out all tabs. Replace with blanks.
    2151 *    3 n = index(line,tab)
    2152 *      if(n.eq.0) goto 4
    2153 *        write(line(n:n),'(a)') blank
    2154 *        goto 3
    2155 *    4 continue
    2156 **
    2157 ** Search for matching braces. Keep looping until all braces are found
    2158 *   10   first_left = 0
    2159 *        last_right = 0
    2160 *        l_br = 0
    2161 *        r_br = 0
    2162 *        i = 0
    2163 *   20   i = i + 1
    2164 *          if(i.gt.lin_len) goto 30
    2165 *          if(line(i:i).eq.'{') then
    2166 *            if(first_left.eq.0) first_left = i
    2167 *            l_br = l_br + 1
    2168 *          else if(line(i:i).eq.'}') then
    2169 *            last_right = i
    2170 *            r_br = r_br + 1
    2171 *          endif
    2172 *          if(r_br.gt.l_br) goto 999
    2173 *          if(l_br.eq.0 .or. l_br.gt.r_br) goto 20
    2174 *   30   continue
    2175 *        if(r_br.ne.l_br) goto 999
    2176 **
    2177 ** if no more comments, then left justify the line
    2178 *        if(l_br.eq.0) goto 50
    2179 ** else, remove comments, and shunt data to the left
    2180 *        offset = last_right - first_left
    2181 *        i = first_left
    2182 *        do 60 j = last_right + 1, lin_len
    2183 *          line(i:i) = line(j:j)
    2184 *          i = i + 1
    2185 *   60   continue
    2186 ** blank out the end of the modified string
    2187 *        do 70 i = lin_len-offset, lin_len
    2188 *          line(i:i) = blank
    2189 *   70   continue
    2190 ** repeat until no more braces are found
    2191 *      goto 10
    2192 **
    2193 ** left adjust line
    2194 *   50 do 80 i = 1, lin_len
    2195 *        if((line(i:i).ne.blank)) then
    2196 *          do 90 j = 1, lin_len-i+1
    2197 *            line(j:j) = line(i+j-1:i+j-1)
    2198 *   90     continue
    2199 *          do 100 j = lin_len-i+2,lin_len
    2200 *            line(j:j) = blank
    2201 *  100     continue
    2202 *          goto 110
    2203 *        endif
    2204 *   80 continue
    2205 ** If, after all this, the line is empty, repeat with next next line.
    2206 *  110 if(line(1:1).eq.blank) goto 1
    2207 **
    2208 *      return
    2209 **
    2210 *  998 write(op,401) 'ERROR: Data line exceeds ',lin_len, ' characters.'
    2211 *      return 1
    2212 *  999 write(op,400) 'ERROR: Unbalanced braces in data file.'
    2213 *  200 return 1
    2214 **
    2215 *  400 format(1x, a)
    2216 *  401 format(1x, a, i3, a)
    2217 *      end
    2218 **
    2219 ** ______________________________________________________________________
    22201378* Title: GETSAD
    22211379* Author: MMJT
    … …  
    39863144*
    39873145* ______________________________________________________________________
    3988 ** Title: GOINTR
    3989 ** Author: MMJT
    3990 ** Date: 23 Oct 1989
    3991 ** Description: This subroutine sets up the parameters necessary to
    3992 ** integrate over an interval of reciprocal space.
    3993 **
    3994 **      ARGUMENTS:
    3995 **            ok  -  logical flag indicating all went well. (output).
    3996 **
    3997 **      COMMON VARIABLES:
    3998 **            uses:  cntrl, CFile
    3999 **
    4000 **        modifies:  no COMMON variables are modified
    4001 ** ______________________________________________________________________
    4002 **
    4003 *      subroutine GOINTR(ok)
    4004 *      include 'DIFFaX.par'
    4005 *      include 'DIFFaX.inc'
    4006 **
    4007 *      logical ok
    4008 **
    4009 *      integer*4 i
    4010 *      real*8 AGLQ16, GLQ16
    4011 **
    4012 ** external functions
    4013 *      external AGLQ16, GLQ16
    4014 ** external subroutine (Some compilers need them declared external)
    4015 **      external INTEGR
    4016 **
    4017 *    1 write(op,100) ' '
    4018 *      write(op,100) 'INTEGRATING INTENSITY IN AN INTERVAL ALONG l. . .'
    4019 *      write(op,100) 'Enter 1 for adaptive quadrature.'
    4020 *      read(cntrl,*,err=1) i
    4021 *      if(CFile) write(op,101) i
    4022 *      if(i.eq.1) then
    4023 *        call INTEGR(AGLQ16,ok)
    4024 *      else
    4025 *        call INTEGR(GLQ16,ok)
    4026 *      endif
    4027 *      if(.not.ok) goto 999
    4028 *      i = 1
    4029 *    2 write(op,100) 'Enter 1 to integrate another interval.'
    4030 *      read(cntrl,*,err=2) i
    4031 *      if(CFile) write(op,101) i
    4032 *      if(i.eq.1) goto 1
    4033 **
    4034 *  999 return
    4035 *  100 format(1x, a)
    4036 *  101 format(1x, i3)
    4037 *      end
    4038 **
    4039 ** ______________________________________________________________________
    4040 ** Title: GOSADP
    4041 ** Author: MMJT
    4042 ** Date: 23 Oct 1989; 1st Mar 2000
    4043 ** Description: This subroutine sets up a file whose name is given
    4044 ** by 'outfile', which is derived from the input filename 'infile'.
    4045 ** The selected area diffraction pattern (SADP) is calculated and is
    4046 ** then written in binary format to 'outfile'.
    4047 **
    4048 **      ARGUMENTS:
    4049 **            infile   -  name of input file (input).
    4050 **            outfile  -  name of output file (output).
    4051 **            ok       -  logical flag indicating all went well.
    4052 **                                                      (output).
    4053 **
    4054 **      COMMON VARIABLES:
    4055 **            uses:  cntrl, CFile
    4056 **
    4057 **        modifies:  loglin, brightness
    4058 ** ______________________________________________________________________
    4059 **
    4060 *      subroutine GOSADP(infile,outfile,ok)
    4061 *      include 'DIFFaX.par'
    4062 *      include 'DIFFaX.inc'
    4063 **
    4064 *      character*(*) infile, outfile
    4065 *      logical ok
    4066 **
    4067 *      integer*4 i_adapt, i_plane, io_err, hk_lim
    4068 *      real*8 AGLQ16, GLQ16, l_upper
    4069 **
    4070 ** external functions
    4071 *      external AGLQ16, GLQ16
    4072 ** external subroutines (Some compilers need them declared external)
    4073 **      external STREAK, GETFNM, GETSAD, WRTSAD
    4074 **
    4075 ** Get a suitable file name, and open file.
    4076 *      call GETFNM(infile,outfile,'sadp',ok)
    4077 *      if(.not.ok) goto 999
    4078 ** Open unformatted for binary write.
    4079 *      if(sa.ne.op) open(unit = sa, file = outfile, status = 'new',
    4080 *     |       form = 'unformatted', err = 990, iostat = io_err)
    4081 ** some compilers need the following added for true binary output
    4082 **     |       ,recordtype = 'stream')
    4083 **
    4084 *      write(op,100) ' '
    4085 *      write(op,100)'CALCULATING SELECTED AREA DIFFRACTION PATTERN. . .'
    4086 **
    4087 * 1234 write(op,100) 'Enter 1 for adaptive quadrature.'
    4088 *      read(cntrl,*,err=1234) i_adapt
    4089 *      if(CFile) write(op,101) i_adapt
    4090 **
    4091 *    1 write(op,100) 'Choose a plane in reciprocal space to view.'
    4092 *      write(op,100) '       the l-axis is included by default.'
    4093 *      write(op,100) '1: k = 0.   2: h = 0.   3: h = k.   4: h = -k.'
    4094 *      read(cntrl,*,err=1) i_plane
    4095 *      if(CFile) write(op,101) i_plane
    4096 *      if(i_plane.lt.1 .or. i_plane.gt.4) then
    4097 *        if(CFile) then
    4098 *          write(op,100) 'Illegal choice. Default is k = 0.'
    4099 *          i_plane = 1
    4100 *        else
    4101 *          write(op,100) 'Illegal choice. Choose again.'
    4102 *          goto 1
    4103 *        endif
    4104 *      endif
    4105 **
    4106 ** Get upper bounds. The SADP will be square, so we only need one bound.
    4107 ** Choose l-axis, since this is common to all 4 views and makes scaling
    4108 ** simpler if more than one SAD pattern is requested.
    4109 * 2345 write(op,100) 'Enter maximum value of l.'
    4110 *      read(cntrl,*,err=2345) l_upper
    4111 *      if(CFile) write(op,104) l_upper
    4112 **
    4113 ** 8-bit images or 16-bit images?
    4114 * 3456 write(op,100) 'Choose the bit-depth of the image.'
    4115 *      write(op,100) '   8: - 8-bits  16: - 16-bits'
    4116 *      read(cntrl,*,err=3456) bitdepth
    4117 *      if(CFile) write(op,101) bitdepth
    4118 *      if(bitdepth.ne.8 .and. bitdepth.ne.16) then
    4119 *        write(op,100) 'Illegal bit-depth.'
    4120 *        if(CFile) then
    4121 *          write(op,100) 'Using 8-bit as the default'
    4122 *          bitdepth = 8
    4123 *        else
    4124 *          write(op,100) 'Re-enter. . .'
    4125 *          goto 3456
    4126 *        endif
    4127 *      endif
    4128 *      if(bitdepth.eq.16) then
    4129 ** Bypass issue of signed or unsigned format. Use range 0 - 32767 only.
    4130 *        write(op,100)
    4131 *     | 'File will be saved in unsigned 16-bit (big-endian) format.'
    4132 *        maxsad = FIFTEENBITS - ONE
    4133 *      else
    4134 *        write(op,100) 'File will be saved in unsigned 8-bit format.'
    4135 *        maxsad = EIGHTBITS - ONE
    4136 *      endif
    4137 **
    4138 ** Logarithmic or linear intensity scaling?
    4139 *    2 write(op,100) 'Choose intensity scaling'
    4140 *      write(op,100) '   0: - Logarithmic  1: - Linear'
    4141 *      read(cntrl,*,err=2) loglin
    4142 *      if(CFile) write(op,101) loglin
    4143 *      if(loglin.lt.0 .or. loglin.gt.1) then
    4144 *        write(op,100) 'Illegal intensity scaling type.'
    4145 *        if(CFile) then
    4146 *          write(op,100) 'Using linear as the default'
    4147 *          loglin = 1
    4148 *        else
    4149 *          write(op,100) 'Re-enter. . .'
    4150 *          goto 2
    4151 *        endif
    4152 *      endif
    4153 **
    4154 ** By how much do we wish to saturate?
    4155 *    3 write(op,100) 'Enter a brightness (must be +ve)'
    4156 *      if(bitdepth.eq.16) then
    4157 *        write(op,100)'  1 - scales intensities to the range 0 - 32767'
    4158 *        write(op,100)' 10 - scales intensities to the range 0 - 327670,'
    4159 *        write(op,100)'      but values above 65535 will be saturated'
    4160 *      else
    4161 *        write(op,100)'  1 - scales intensities to the range 0 - 255'
    4162 *        write(op,100)' 10 - scales intensities to the range 0 - 2550,'
    4163 *        write(op,100)'      but values above 255 will be saturated'
    4164 *      endif
    4165 *      read(cntrl,*,err=3) brightness
    4166 *      if(CFile) write(op,104) brightness
    4167 *      if(brightness.le.ZERO) then
    4168 *        write(op,100) 'Illegal value for brightness. Must be positive'
    4169 *        if(CFile) then
    4170 *          if(loglin.eq.0) then
    4171 *            write(op,100) 'Using default of 1 for logarithmic scaling'
    4172 *            brightness = ONE
    4173 *          else
    4174 *            write(op,100) 'Using default of 10 for linear scaling'
    4175 *            brightness = TEN
    4176 *          endif
    4177 *        else
    4178 *          if(loglin.eq.0) then
    4179 *            write(op,100) '1 is a good default for logarithmic scaling'
    4180 *          else
    4181 *            write(op,100) '10 is a good default for linear scaling'
    4182 *          endif
    4183 *          write(op,100) 'Re-enter. . .'
    4184 *          goto 3
    4185 *        endif
    4186 *      endif
    4187 **
    4188 ** Generate intensity data for the SAD pattern.
    4189 *      if(i_adapt.eq.1) then
    4190 *        call GETSAD(AGLQ16, i_plane, l_upper, hk_lim, infile, ok)
    4191 *      else
    4192 *        call GETSAD( GLQ16, i_plane, l_upper, hk_lim, infile, ok)
    4193 *      endif
    4194 **
    4195 *      if(ok) call WRTSAD(outfile, i_plane, l_upper, hk_lim, ok)
    4196 **
    4197 *  999 return
    4198 *  990 write(op,102) 'Problems opening output file ', outfile
    4199 *      write(op,103) 'IOSTAT = ', io_err
    4200 *      ok = .false.
    4201 *      return
    4202 *  100 format(1x, a)
    4203 *  101 format(1x, i3)
    4204 *  102 format(1x, 2a)
    4205 *  103 format(1x, a, i5)
    4206 *  104 format(1x, g12.5)
    4207 *      end
    4208 **
    4209 ** ______________________________________________________________________
    4210 ** Title: GOSPEC
    4211 ** Author: MMJT
    4212 ** Date: 23 Oct 1989
    4213 ** Description: This subroutine sets up a file whose name is given
    4214 ** by 'outfile', which is derived from the input filename 'infile'.
    4215 ** The powder pattern data is then written to 'outfile', which is closed
    4216 ** by the subroutine 'WRTSPC'.
    4217 **
    4218 **      ARGUMENTS:
    4219 **            infile   -  name of input file (input)
    4220 **            outfile  -  name of output file (output)
    4221 **            ok       -  logical flag indicating all went well.
    4222 **                                                      (output).
    4223 **
    4224 **      COMMON VARIABLES:
    4225 **            uses:  blurring, CFile, GAUSS, LORENZ, NONE, PS_VGT
    4226 **                   PV_GSS, PV_LRN, cntrl
    4227 **
    4228 **        modifies:  full_brd
    4229 ** ______________________________________________________________________
    4230 **
    4231 *      subroutine GOSPEC(infile,outfile,ok)
    4232 *      include 'DIFFaX.par'
    4233 *      include 'DIFFaX.inc'
    4234 **
    4235 *      character*(*) infile, outfile
    4236 *      logical ok
    4237 **
    4238 *      logical GETSPC, TRMSPC, RDRNGE
    4239 *      integer*4 io_err
    4240 *      real*8 AGLQ16, GLQ16, cut_off
    4241 **
    4242 ** external functions
    4243 *      external RDRNGE, GETSPC, AGLQ16, GLQ16, TRMSPC
    4244 ** external subroutines (Some compilers need them declared external)
    4245 **      external GETFNM, PV, GAUSSN, LORNZN, WRTSPC
    4246 **
    4247 *      write(op,100) ' '
    4248 *      write(op,100) 'CALCULATING POWDER DIFFRACTION PATTERN. . .'
    4249 **
    4250 ** get angular range and step
    4251 *      ok = RDRNGE()
    4252 *      if(.not.ok) goto 999
    4253 ** create a new filename to save spectrum data to
    4254 *      call GETFNM(infile, outfile, 'spc', ok)
    4255 *      if(.not.ok) goto 999
    4256 *      if(sp.ne.op) open(unit = sp, file = outfile, status = 'new',
    4257 *     |            err = 990, iostat = io_err)
    4258 *      full_brd = 1
    4259 * 3456 write(op,100) 'Enter 1 for adaptive quadrature on broad peaks.'
    4260 *      read(cntrl,*,err=3456) full_brd
    4261 *      if(CFile) write(op,101) full_brd
    4262 *      if(full_brd.eq.1) then
    4263 *        ok = GETSPC(AGLQ16, infile)
    4264 *      else
    4265 *        ok = GETSPC(GLQ16, infile)
    4266 *      endif
    4267 ** suppress the huge peak near the origin if required
    4268 *      cut_off = ZERO
    4269 *      if(ok.and.(th2_min.eq.ZERO).and.trim_origin) ok = TRMSPC(cut_off)
    4270 *      if(ok) then
    4271 *        if(blurring.eq.GAUSS) then
    4272 *          write(op,104) 'Gaussian'
    4273 *          call GAUSSN(cut_off)
    4274 *        else if(blurring.eq.LORENZ) then
    4275 *          write(op,104) 'Lorentzian'
    4276 *          call LORNZN(cut_off)
    4277 *        else if(blurring.eq.PS_VGT) then
    4278 *          write(op,104) 'pseudo-Voigt'
    4279 *          call PV(cut_off)
    4280 *        else if(blurring.eq.PV_GSS) then
    4281 *          write(op,104) 'Gaussian'
    4282 *          call PV(cut_off)
    4283 *        else if(blurring.eq.PV_LRN) then
    4284 *          write(op,104) 'Lorentzian'
    4285 *          call PV(cut_off)
    4286 *        else if(blurring.ne.NONE) then
    4287 *          write(op,100)
    4288 *     |         'Instrumental broadening type is undefined in GOSPEC.'
    4289 *        endif
    4290 *      endif
    4291 *      if(ok) call WRTSPC(outfile, ok)
    4292 **
    4293 *  999 return
    4294 *  990 write(op,102) 'Problems opening output file ', outfile
    4295 *      write(op,103) 'IOSTAT = ', io_err
    4296 *      ok = .false.
    4297 *      return
    4298 *  100 format(1x, a)
    4299 *  101 format(1x, i3)
    4300 *  102 format(1x, 2a)
    4301 *  103 format(1x, a, i5)
    4302 *  104 format(1x, 'Adding ', a, ' instrumental broadening . . .')
    4303 *      end
    4304 **
    4305 ** ______________________________________________________________________
    4306 ** Title: GOSTRK
    4307 ** Author: MMJT
    4308 ** Date: 23 Oct 1989
    4309 ** Description: This subroutine sets up a file whose name is given
    4310 ** by 'outfile', which is derived from the input filename 'infile'.
    4311 ** The streak intensity is then written to 'outfile', which is closed
    4312 ** by the subroutine 'STREAK'.
    4313 **
    4314 **      ARGUMENTS:
    4315 **            infile   -  name of input file. (input).
    4316 **            outfile  -  name of output file. (output).
    4317 **            ok       -  logical flag indicating all went well.
    4318 **                                                       (output).
    4319 **
    4320 **      COMMON VARIABLES:
    4321 **            uses:  cntrl, CFile
    4322 **
    4323 **        modifies:  no COMMON variables are modified
    4324 ** ______________________________________________________________________
    4325 **
    4326 *      subroutine GOSTRK(infile,outfile,ok)
    4327 *      include 'DIFFaX.par'
    4328 *      include 'DIFFaX.inc'
    4329 **
    4330 *      character*(*) infile, outfile
    4331 *      logical ok
    4332 **
    4333 *      integer*4 i, io_err
    4334 *      real*8 AGLQ16, GLQ16
    4335 **
    4336 ** external functions
    4337 *      external AGLQ16, GLQ16
    4338 ** external subroutines (Some compilers need them declared external)
    4339 **      external STREAK, GETFNM
    4340 **
    4341 *      call GETFNM(infile, outfile, 'str', ok)
    4342 *      if(.not.ok) goto 999
    4343 *      if(sk.ne.op) open(unit = sk, file = outfile, status = 'new',
    4344 *     |                              err = 990, iostat = io_err)
    4345 * 4567 write(op,100) ' '
    4346 *      write(op,100) 'CALCULATING INTENSITY ALONG A STREAK. . .'
    4347 *      write(op,100) 'Enter 1 for adaptive quadrature'
    4348 *      read(cntrl,*,err=4567) i
    4349 *      if(CFile) write(op,101) i
    4350 ** select integration function
    4351 *      if(i.eq.1) then
    4352 *        call STREAK(AGLQ16, outfile, ok)
    4353 *      else
    4354 *        call STREAK( GLQ16, outfile, ok)
    4355 *      endif
    4356 **
    4357 *  999 return
    4358 *  990 write(op,102) 'Problems opening output file ', outfile
    4359 *      write(op,103) 'IOSTAT = ', io_err
    4360 *      ok = .false.
    4361 *      return
    4362 *  100 format(1x, a)
    4363 *  101 format(1x, i3)
    4364 *  102 format(1x, 2a)
    4365 *  103 format(1x, a, i5)
    4366 *      end
    4367 **
    4368 ** ______________________________________________________________________
    43693146* Title: HKL_LIM
    43703147* Author: MMJT
    … …  
    47333510      return
    47343511      end
    4735 *
    4736 * ______________________________________________________________________
    4737 ** Title: LAYCNT
    4738 ** Author: MMJT
    4739 ** Date: 3 Oct 1989
    4740 ** Description: Counts the number of integer arguments in the character
    4741 ** string 'line', and returns the answer in LAYCNT. Legal separators
    4742 ** are, blanks, tabs, and commas. If LAYCNT encounters anything other
    4743 ** than a digit (0 - 9), or one of the three legal separators, LAYCNT
    4744 ** is set to -1. This routine is similar to CNTARG, accept it will only
    4745 ** accept integers as arguments.
    4746 **
    4747 **      ARGUMENTS:
    4748 **            line  -  Line of characters to be parsed. (input).
    4749 **
    4750 **      LAYCNT returns the number of integer arguments in the line.
    4751 **      Returns -1 if an illegal character was detected.
    4752 ** ______________________________________________________________________
    4753 **
    4754 *      integer*4 function LAYCNT(line)
    4755 *      implicit none
    4756 **
    4757 *      character*(*) line
    4758 **
    4759 *      logical in_num, legit
    4760 *      integer*4 blank, tab, comma, nought, nine, i, j, num_cnt, lin_len
    4761 *      parameter (tab=9, blank=32, comma=44, nought=48, nine=57)
    4762 **
    4763 *      in_num = .false.
    4764 *      legit = .true.
    4765 *      num_cnt = 0
    4766 *      LAYCNT = -1
    4767 **
    4768 *      lin_len = len(line)
    4769 *      do 10 i = 1, lin_len
    4770 *        j = ichar(line(i:i))
    4771 *        if(j.eq.tab .or. j.eq.blank .or. j.eq.comma) then
    4772 *          in_num = .false.
    4773 *        else if(j.ge.nought .and. j.le.nine) then
    4774 *          if(.not.in_num) then
    4775 *            in_num = .true.
    4776 *            num_cnt = num_cnt + 1
    4777 *          endif
    4778 *        else
    4779 ** illegal character
    4780 *          legit = .false.
    4781 *        endif
    4782 *   10 continue
    4783 **
    4784 *      if(legit) LAYCNT = num_cnt
    4785 **
    4786 *      return
    4787 *      end
    47883512*
    47893513* ______________________________________________________________________
    … …  
    54414165*
    54424166* ______________________________________________________________________
    5443 ** Title: NAMER
    5444 ** Author: MMJT
    5445 ** Date: 13 Oct 1988
    5446 ** Description: This subroutine creates the appropriate filenames for
    5447 ** the various files needed.
    5448 ** It generates the filename by reading the structure data file name.
    5449 ** It scans the name to see if a period ('.') is already within
    5450 ** the name (ie. if the data file is called 'Structure.dat'). It deletes
    5451 ** the characters after the period and appends whatever is in append.
    5452 ** 'name1' and 'name2' are assumed to be the same physical length.
    5453 ** This subroutine is called by GETFNM
    5454 **
    5455 **      ARGUMENTS:
    5456 **            name1   -  The name of the input data file. (input).
    5457 **            name2   -  A derivative filename. (output).
    5458 **            append  -  A token to be appended to name2. (input).
    5459 ** ______________________________________________________________________
    5460 **
    5461 *      subroutine NAMER(name1,name2,append)
    5462 *      include 'DIFFaX.par'
    5463 **     save
    5464 **
    5465 *      character*(*) name1, name2, append
    5466 **
    5467 *      integer*4 applen, namlen, i, idot
    5468 *      integer*4 LENGTH
    5469 **
    5470 ** external function
    5471 *      external LENGTH
    5472 **
    5473 ** get length of the string holding the filename
    5474 *      namlen = len(name1)
    5475 *      name2  = ' '
    5476 *      applen = LENGTH(append)
    5477 **
    5478 *      idot = index(name1,'.')
    5479 *      if(idot.eq.0) idot = index(name1,' ')
    5480 *      if(idot.eq.0) then
    5481 *        if(namlen.lt.MAX_NAM) then
    5482 *          idot = namlen + 1
    5483 *        else
    5484 *          idot = MAX_NAM
    5485 *          namlen = MAX_NAM
    5486 *        endif
    5487 *      endif
    5488 ** truncate root filename so that the appendage will appear correctly
    5489 *      if(idot+applen.gt.namlen) idot = namlen - applen
    5490 **
    5491 *      do 10 i = 1, idot-1
    5492 *        write(name2(i:i),'(a)') name1(i:i)
    5493 *   10 continue
    5494 **
    5495 *      write(name2(idot:idot),'(a)') '.'
    5496 **
    5497 *      do 20 i = 1, applen
    5498 *        write(name2((idot+i):(idot+i)),'(a)') append(i:i)
    5499 *   20 continue
    5500 **
    5501 *      do 30 i = idot+applen+1, namlen
    5502 *        write(name2(i:i),'(a)') ' '
    5503 *   30 continue
    5504 **
    5505 *      return
    5506 *      end
    5507 **
    5508 ** ______________________________________________________________________
    55094167* Title: NMCOOR
    55104168* Author: MWD
    … …  
    55404198*
    55414199* ______________________________________________________________________
    5542 ** Title: NXTARG
    5543 ** Author: MMJT
    5544 ** Date: 21 JUL 1997
    5545 ** Description: Removes the next argument from the character string
    5546 ** 'line', and returns that argument in the character string 'arg'.
    5547 **
    5548 **      ARGUMENTS:
    5549 **            line  -  Line of characters to be parsed. (input).
    5550 **            arg   -  First argument found in 'line'.  (output).
    5551 **
    5552 **      Returns  1 if an argument was found.
    5553 **      Returns  0 if no argument was found.
    5554 **      Returns <0 if an error occurred.
    5555 ** ______________________________________________________________________
    5556 **
    5557 *      integer*4 function NXTARG(line,arg)
    5558 *      implicit none
    5559 *      save
    5560 **
    5561 *      character*(*) line, arg
    5562 **
    5563 *      character*200 tmpline
    5564 *      integer*4 blank, tab, comma, i, j, i1, i2, lin_len
    5565 *      parameter (tab=9, blank=32, comma=44)
    5566 **
    5567 *      NXTARG = 0
    5568 **
    5569 *      lin_len = len(line)
    5570 **
    5571 *      i = 0
    5572 ** Strip leading blanks
    5573 *   10 i = i + 1
    5574 *        if(i.gt.lin_len) goto 997
    5575 *        j = ichar(line(i:i))
    5576 *        if(j.eq.tab .or. j.eq.blank .or. j.eq.comma) goto 10
    5577 *      i1 = i
    5578 **
    5579 ** Read off the non-blank characters
    5580 *   20 i = i + 1
    5581 *        if(i.gt.lin_len) goto 30
    5582 *        j = ichar(line(i:i))
    5583 *        if(j.eq.tab .or. j.eq.blank .or. j.eq.comma) then
    5584 *          NXTARG = 1
    5585 *          goto 30
    5586 *        else
    5587 *          goto 20
    5588 *        endif
    5589 *   30 continue
    5590 **
    5591 *      i2 = i
    5592 *      if(i2.lt.i1) goto 998
    5593 *      write(arg, '(a)', err=999) line(i1:i2)
    5594 *      write(tmpline, '(a)', err=999) line(i2:lin_len)
    5595 *      write(line, '(a)', err=999) tmpline(1:lin_len)
    5596 **
    5597 *      return
    5598 *  997 NXTARG = -3
    5599 *      return
    5600 *  998 NXTARG = -2
    5601 *      return
    5602 *  999 NXTARG = -1
    5603 *      return
    5604 *      end
    5605 **
    5606 ** ______________________________________________________________________
    56074200* Title: OPTIMZ
    56084201* Author: MWD and MMJT
    … …  
    60864679*
    60874680* ______________________________________________________________________
    6088 ** Title: PARENT
    6089 ** Author: MMJT
    6090 ** Date: 10 Aug 1989
    6091 ** Description: This routine detects whether or not a character
    6092 ** string contains any parentheses, '(' and ')'. If it does, they are
    6093 ** removed. PARENT is used by RDPRMS to detect whether or not the
    6094 ** anisotropic temperature factors are being used.
    6095 **
    6096 **      ARGUMENTS:
    6097 **            line  -  Line of characters to parse. (input).
    6098 **                  -  Line stripped of any parentheses. (output).
    6099 **
    6100 **      PARENT returns 0 if no parentheses were found, 1 if a balanced
    6101 **      pair of parentheses were found, and -1 if there was an error.
    6102 ** ______________________________________________________________________
    6103 **
    6104 *      integer*4 function PARENT(line)
    6105 *      implicit none
    6106 **
    6107 *      character*(*) line
    6108 **
    6109 *      integer*4 i, j
    6110 **
    6111 *      PARENT = 0
    6112 *      i = index(line,'(')
    6113 *      j = index(line,')')
    6114 ** no parentheses
    6115 *      if(i.eq.0 .and. j.eq.0) goto 200
    6116 ** right hand parenthesis occurred before the left hand parenthesis!
    6117 *      if(j.lt.i) goto 999
    6118 ** only one parenthesis
    6119 *      if(i.ne.j .and. (i.eq.0 .or. j.eq.0)) goto 999
    6120 ** blot out the two parentheses found
    6121 *      write(line(i:i),'(a)') ' '
    6122 *      write(line(j:j),'(a)') ' '
    6123 ** check there are not any more
    6124 *      i = index(line,'(')
    6125 *      j = index(line,')')
    6126 *      if(i.eq.0 .and. j.eq.0) then
    6127 ** there were no more parentheses. Data is probably ok.
    6128 *        PARENT = 1
    6129 *      else
    6130 ** there were more parentheses. Data is suspect.
    6131 *        goto 999
    6132 *      endif
    6133 **
    6134 *  200 return
    6135 *  999 PARENT = -1
    6136 *      return
    6137 *      end
    6138 **
    6139 ** ______________________________________________________________________
    61404681* Title: PNTINT
    61414682* Author: MMJT
    … …  
    62084749*
    62094750* ______________________________________________________________________
    6210 ** Title: POINT
    6211 ** Author: MWD and MMJT
    6212 ** Date: 18 Aug 1988
    6213 ** Description:  This subroutine prompts the user for h, k, l values
    6214 ** and displays the intensity at that point.
    6215 **
    6216 **      ARGUMENTS:
    6217 **            ok  -  logical flag indicating all went well. (output).
    6218 **
    6219 **      COMMON VARIABLES:
    6220 **            uses:  a0, b0, c0, d0, lambda, cntrl, CFile, xplcit,
    6221 **                   recrsv, rad_type, n_layers, inf_thick, wavefn
    6222 **
    6223 **        modifies:  no COMMON variables are modified
    6224 ** ______________________________________________________________________
    6225 **
    6226 *      subroutine POINT(ok)
    6227 *      include 'DIFFaX.par'
    6228 *      include 'DIFFaX.inc'
    6229 **
    6230 *      logical ok
    6231 **
    6232 *      logical GET_S, GET_S2
    6233 *      integer*4 h, k, i
    6234 *      character chr
    6235 *      real*8 l, INTENS, INTEN2, Shkl, x, ANGLE, SS, W4, theta, Q2
    6236 *      complex*16 f(MAX_L), s(MAX_L)
    6237 **
    6238 ** external functions
    6239 *      external INTENS, INTEN2, GET_S, GET_S2
    6240 ** external subroutines (Some compilers need them declared external)
    6241 **      external XYPHSE, PRE_MAT, GET_MAT, GET_F
    6242 **
    6243 ** statement functions
    6244 ** SS is the value of 1/d**2 at hkl
    6245 *      SS(h,k,l) = h*h*a0 + k*k*b0 + l*l*c0 + h*k*d0
    6246 ** ANGLE is the Bragg angle (in radians) of the h,k,l plane
    6247 *      ANGLE(h,k,l) = asin(HALF * lambda * sqrt(SS(h,k,l)))
    6248 ** W4 is the X-ray polarization factor
    6249 *      W4(theta) = HALF * (ONE + (cos(TWO*theta))**2)
    6250 **
    6251 *    1 write(op,400) ' '
    6252 *      write(op,400) 'CALCULATING INTENSITY AT A POINT. . .'
    6253 *   10 write(op,400) 'Enter h, k, l'
    6254 *      read(cntrl,*,err=10)  h, k, l
    6255 *      if(CFile) write(op,401) h, k, l
    6256 ** check angles are meaningful
    6257 *      Q2 = FOUR / (lambda**2)
    6258 *      Shkl = SS(h,k,l)
    6259 ** Check that SS(h,k,l) is legal
    6260 *      if(Shkl.lt.ZERO) then
    6261 *        write(op,403) 'ERROR: In POINT(), 1/d**2 = ', Shkl
    6262 *        ok = .false.
    6263 *        goto 999
    6264 *      endif
    6265 *      if(Shkl.gt.Q2) then
    6266 *        write(op,402) h, k, l,' exceeds 180 degree scattering angle!'
    6267 *        write(op,400) 'Re-enter. . .'
    6268 *        goto 10
    6269 *      endif
    6270 ** make sure we are not going to blow up at the origin
    6271 *      if(h.eq.0 .and. k.eq.0 .and. rad_type.eq.ELECTN) then
    6272 *        if(Shkl.le.eps4) then
    6273 *          write(op,400)
    6274 *     |   'Cannot integrate across the origin for electron radiation'
    6275 *          write(op,400) 'Re-enter. . .'
    6276 *          goto 10
    6277 *        endif
    6278 *      endif
    6279 ** get intensity
    6280 *      call XYPHSE(h, k)
    6281 *      call PRE_MAT(h, k)
    6282 *      call GET_F(f, Shkl, l)
    6283 *      if(recrsv) then
    6284 *        x = INTENS(f, h, k, l, ok)
    6285 *        if(.not.ok) goto 999
    6286 ** set up mat again to re-call GET_S (or GET_S2)
    6287 *        call PRE_MAT(h, k)
    6288 *        call GET_MAT(h, k, l)
    6289 *        if(inf_thick) then
    6290 ** initialize s to -f, since mat is -(ident - T)
    6291 *          do 20 i = 1, n_layers
    6292 *            s(i) = - f(i)
    6293 *   20     continue
    6294 *          ok = GET_S(f, s, h, k, l)
    6295 *        else
    6296 ** s is initialized within GET_S2, which then calls GET_S
    6297 *          ok = GET_S2(f, s, h, k, l)
    6298 *        endif
    6299 *      else
    6300 *        x = INTEN2(f, h, k, l, ok)
    6301 *      endif
    6302 *      if(.not.ok) goto 999
    6303 ** for diagnostics write out the s values as well
    6304 *      if(rad_type.eq.X_RAY)  x = x * W4(ANGLE(h,k,l))
    6305 *      write(op,404)
    6306 *     |     '2theta = ', RAD2DEG * TWO * ANGLE(h,k,l), ' degrees'
    6307 *      if(Shkl .gt. ZERO) write(op,403) 'd = ', ONE / sqrt(Shkl)
    6308 *      if(Shkl .ge. ZERO) write(op,403) '1/d = ', sqrt(Shkl)
    6309 *      write(op,400) ' '
    6310 *      write(op,400) 'Layer scattering factors'
    6311 *      do 30 i = 1, n_layers
    6312 *        write(op,405) i, f(i)
    6313 *   30 continue
    6314 *      write(op,400) ' '
    6315 *      if(recrsv) then
    6316 *        write(op,400) 'Average scattered wavefunctions'
    6317 *        do 40 i = 1, n_layers
    6318 *          write(op,407) i, s(i)
    6319 *   40   continue
    6320 *      else
    6321 *        write(op,408) 'Crystal wavefunction', wavefn
    6322 *      endif
    6323 *      write(op,400) ' '
    6324 *      write(op,406) 'Intensity at ', h, k, l, ' = ', x
    6325 *      write(op,400) 'Another point? (y/n)'
    6326 *      read(cntrl,'(a)') chr
    6327 *      if(chr(1:1).eq.'y' .or. chr(1:1).eq.'Y') goto 1
    6328 **
    6329 *  999 return
    6330 *  400 format(1x, a)
    6331 *  401 format(1x, 2i3, g12.5)
    6332 *  402 format(1x, 2i3, g12.5, a)
    6333 *  403 format(1x, a, g12.5)
    6334 *  404 format(1x, a, g12.5, a)
    6335 *  405 format(1x, 'f(', i2, ') = (', g14.7, ',', g14.7, ')')
    6336 *  406 format(1x, a, 2i3, g12.5, a, g14.7)
    6337 *  407 format(1x, 'psi(', i2,') = (', g14.7, ',', g14.7, ')')
    6338 *  408 format(1x, a, ' psi = (', g14.7, ',', g14.7, ')')
    6339 *      end
    6340 **
    6341 ** ______________________________________________________________________
    63424751* Title: POLINT
    63434752* Author: MMJT, adapted from Numerical Recipes Software
    … …  
    66185027*
    66195028* ______________________________________________________________________
    6620 ** Title: RDBLUR
    6621 ** Authors: MMJT
    6622 ** Date: 4 Oct 1989
    6623 ** Description: This function reads the type of instrumental
    6624 ** broadening to be simulated, and the associated parameters.
    6625 ** Legal types are 'NONE ', 'PSEUDO-VOIGT ', 'GAUSSIAN '
    6626 ** and 'LORENTZIAN '. PSEUDO-VOIGT requires 4 parameters,
    6627 ** GAUSSIAN requires a standard deviation in 2theta degrees, and
    6628 ** LORENTZIAN requires a halfwidth in 2theta degrees.
    6629 ** If the keyword 'TRIM' is at the end of the data line, intensity
    6630 ** near the origin will be ignored when the instrumental broadening
    6631 ** is added.
    6632 **
    6633 **      ARGUMENTS:
    6634 **            unit_no  -  Logical unit number that the data file is to
    6635 **                        be read from. (input).
    6636 **            messge   -  A short message indicating what went wrong
    6637 **                        (if anything) during the datafile read.
    6638 **                        messge is terminated by a token '$'. (output).
    6639 **
    6640 **      COMMON VARIABLES:
    6641 **        modifies:  blurring, pv_u, pv_v, pv_w, pv_gamma,
    6642 **                   FWHM, GAUSS, LORENZ, NONE, PS_VGT, trim_origin
    6643 **
    6644 **      RDBLUR returns logical .true. if a legal function for
    6645 **      simulating instrumental broadening was entered.
    6646 **
    6647 **      Legal types are:
    6648 **             'NONE'
    6649 **             'PSEUDO-VOIGT' (with u, v, w and gamma parameters)
    6650 **             'GAUSSIAN' (with a standard deviation in degrees)
    6651 **             'LORENTZIAN' (with a half width in degrees)
    6652 ** ______________________________________________________________________
    6653 **
    6654 *      logical function RDBLUR(unit_no, messge)
    6655 *      include 'DIFFaX.par'
    6656 *      include 'DIFFaX.inc'
    6657 **
    6658 *      integer*4 unit_no
    6659 *      character*(*) messge
    6660 **
    6661 *      character*200 line, dummy
    6662 *      character*80 list(4)
    6663 *      integer*4 CHOICE, LENGTH, CNTARG, i, iflag, arg_num, lin_len
    6664 **
    6665 ** external functions
    6666 *      external CHOICE, LENGTH, CNTARG
    6667 ** external subroutines (Some compilers need them declared external)
    6668 **      external GETLNE, WRTLNE, TOUPPR
    6669 **
    6670 *      RDBLUR = .false.
    6671 **
    6672 **
    6673 *      messge = 'Instrumental broadening parameters were not specified.$'
    6674 *      call GETLNE(unit_no, line, *999)
    6675 *      call TOUPPR(line)
    6676 *      lin_len = len(line)
    6677 **
    6678 ** First, check whether or not user wants to trim the low angle
    6679 ** intensities for the broadened spectrum. If present, remove the
    6680 ** keyword 'TRIM', and set flag.
    6681 *      i = index(line, ' TRIM')
    6682 *      if(i.eq.0) then
    6683 *        trim_origin = .false.
    6684 *      else
    6685 *        trim_origin = .true.
    6686 *        write(line(i+1:i+4), '(a)') '    '
    6687 *      endif
    6688 **
    6689 *      list(1) = 'NONE '
    6690 *      list(2) = 'GAUSSIAN '
    6691 *      list(3) = 'LORENTZIAN '
    6692 *      list(4) = 'PSEUDO-VOIGT '
    6693 *      iflag   =  CHOICE(line, list, 4)
    6694 **
    6695 *      if(iflag.lt.1 .or. iflag.gt.4) then
    6696 *        write(op,100)
    6697 *     |  'Illegal instrumental broadening type. Legal types are:'
    6698 *        do 10 i = 1, 4
    6699 *          write(op,101) '               ', list(i)
    6700 *   10   continue
    6701 *        messge = 'Instrumental broadening incorrectly specified.$'
    6702 *        goto 999
    6703 *      endif
    6704 **
    6705 *      if(iflag.eq.1) blurring = NONE
    6706 **
    6707 *      if(iflag.eq.2) then
    6708 **
    6709 ** check the number of arguments
    6710 *        dummy = line(LENGTH('GAUSSIAN')+2:lin_len)
    6711 *        arg_num = CNTARG(dummy)
    6712 *        if(arg_num.ne.1 .and. arg_num.ne.3) then
    6713 *          messge =
    6714 *     |       'Illegal number of GAUSSIAN arguments. FWHM, or u,v,w.$'
    6715 *          goto 999
    6716 *        endif
    6717 **
    6718 *        if(arg_num.eq.1) then
    6719 *          blurring = GAUSS
    6720 **
    6721 *          messge = 'Problems reading Gaussian FWHM.$'
    6722 ** for strict FORTRAN77 compliance
    6723 *          write(scrtch, *, err=990) dummy
    6724 *          rewind(scrtch, err=990)
    6725 *          read(scrtch, *, err=999) FWHM
    6726 *          rewind(scrtch, err=990)
    6727 ** for more lenient, efficient, compilers
    6728 **          read(dummy, *, err=999) FWHM
    6729 **
    6730 *          if(FWHM.lt.ZERO) then
    6731 *            messge = 'Gaussian FWHM is -ve.$'
    6732 *            goto 999
    6733 *          endif
    6734 **
    6735 *          if(FWHM.lt.eps7) then
    6736 *            write(op,102) 'Gaussian FWHM is zero.'
    6737 *            write(op,102)
    6738 *     |        ' No instrumental broadening will be simulated.'
    6739 *            blurring = NONE
    6740 *          endif
    6741 **
    6742 *        else if(arg_num.eq.3) then
    6743 *          blurring = PV_GSS
    6744 **
    6745 *          messge = 'Problems reading Gaussian u, v, and w factors.$'
    6746 ** for strict FORTRAN77 compliance
    6747 *          write(scrtch, *, err=990) dummy
    6748 *          rewind(scrtch, err=990)
    6749 *          read(scrtch, *, err=999) pv_u, pv_v, pv_w
    6750 *          rewind(scrtch, err=990)
    6751 ** for more lenient, efficient, compilers
    6752 **          read(dummy, *, err=999) pv_u, pv_v, pv_w
    6753 **
    6754 *          pv_gamma = ZERO
    6755 **
    6756 *          if(pv_w.lt.ZERO) then
    6757 *            write(op,102)
    6758 *     |        'WARNING: Gaussian w-factor is less than zero.'
    6759 *            write(op,102)
    6760 *     |        '  This may cause the broadening function to fail.'
    6761 *          endif
    6762 **
    6763 *          if(abs(pv_u).lt.eps7 .and. abs(pv_v).lt.eps7 .and.
    6764 *     |       abs(pv_w).lt.eps7) then
    6765 *            write(op,102) 'Gaussian u, v, w factors are zero.'
    6766 *            write(op,102) 'No instrumental broadening will be applied.'
    6767 *            blurring = NONE
    6768 *          endif
    6769 **
    6770 ** are the parameters equivalent to a Gaussian with
    6771 ** constant standard deviation?
    6772 *          if(abs(pv_u).le.eps7 .and. abs(pv_v).le.eps7 .and.
    6773 *     |                                    pv_w.gt.ZERO) then
    6774 *            FWHM = sqrt(pv_w)
    6775 *            write(op,102)
    6776 *     |       'Gaussian factors imply a constant FWHM ', FWHM
    6777 *            blurring = GAUSS
    6778 *          endif
    6779 **
    6780 *        endif
    6781 **
    6782 *      endif
    6783 **
    6784 *      if(iflag.eq.3) then
    6785 **
    6786 ** check the number of arguments
    6787 *        dummy = line(LENGTH('LORENTZIAN')+2:lin_len)
    6788 *        arg_num = CNTARG(dummy)
    6789 *        if(arg_num.ne.1 .and. arg_num.ne.3) then
    6790 *          messge =
    6791 *     |     'Illegal number of LORENTZIAN arguments. FWHM, or u,v,w.$'
    6792 *          goto 999
    6793 *        endif
    6794 **
    6795 *        if(arg_num.eq.1) then
    6796 *          blurring = LORENZ
    6797 **
    6798 *          messge = 'Problems reading Lorentzian width.$'
    6799 ** for strict FORTRAN77 compliance
    6800 *          write(scrtch, *, err=990) dummy
    6801 *          rewind(scrtch, err=990)
    6802 *          read(scrtch, *, err=999) FWHM
    6803 *          rewind(scrtch, err=990)
    6804 ** for more lenient, efficient, compilers
    6805 **          read(dummy, *, err=999) FWHM
    6806 **
    6807 *          if(FWHM.lt.ZERO) then
    6808 *            messge =
    6809 *     |        'Negative value for Lorentzian FWHM. Must be positive.$'
    6810 *            goto 999
    6811 *          endif
    6812 **
    6813 *          if(FWHM.lt.eps7) then
    6814 *            write(op,102) 'Lorentzian FWHM is zero.'
    6815 *            write(op,102)
    6816 *     |      ' No instrumental broadening will be simulated.'
    6817 *            blurring = NONE
    6818 *          endif
    6819 **
    6820 *        else if(arg_num.eq.3) then
    6821 *          blurring = PV_LRN
    6822 **
    6823 *          messge = 'Problems reading Lorentzian u, v, and w factors.$'
    6824 ** for strict FORTRAN77 compliance
    6825 *          write(scrtch, *, err=990) dummy
    6826 *          rewind(scrtch, err=990)
    6827 *          read(scrtch, *, err=999) pv_u, pv_v, pv_w
    6828 *          rewind(scrtch, err=990)
    6829 ** for more lenient, efficient, compilers
    6830 **          read(dummy, *, err=999) pv_u, pv_v, pv_w
    6831 **
    6832 *          pv_gamma = ONE
    6833 **
    6834 *          if(pv_w.lt.ZERO) then
    6835 *            write(op,102)
    6836 *     |        'WARNING: Lorentzian w-factor is less than zero.'
    6837 *            write(op,102)
    6838 *     |        '  This may cause the broadening function to fail.'
    6839 *          endif
    6840 **
    6841 *          if(pv_u.lt.eps7 .and. pv_v.lt.eps7 .and.
    6842 *     |                         abs(pv_w).lt.eps7) then
    6843 *            write(op,102) 'Lorentzian u, v, w factors are zero.'
    6844 *            write(op,102) 'No instrumental broadening will be applied.'
    6845 *            blurring = NONE
    6846 *          endif
    6847 **
    6848 ** are the parameters equivalent to a Gaussian with
    6849 ** constant standard deviation?
    6850 *          if(pv_u.le.eps7 .and. pv_v.le.eps7 .and. pv_w.gt.ZERO) then
    6851 *            FWHM = sqrt(pv_w)
    6852 *            write(op,103)
    6853 *     |       'Lorentzian factors imply a constant FWHM ', FWHM
    6854 *            blurring = LORENZ
    6855 *          endif
    6856 **
    6857 *        endif
    6858 **
    6859 *      endif
    6860 **
    6861 *  500 continue
    6862 **
    6863 *      if(iflag.eq.4) then
    6864 *        blurring = PS_VGT
    6865 **
    6866 ** check the number of arguments
    6867 *        dummy = line(LENGTH('PSEUDO-VOIGT')+2:len(line))
    6868 *        arg_num = CNTARG(dummy)
    6869 *        if(arg_num.ne.4) then
    6870 *          messge = 'Exactly 4 pseudo-Voigt arguments are required.$'
    6871 *          goto 999
    6872 *        endif
    6873 **
    6874 *        messge =
    6875 *     |      'Problems reading pseudo-Voigt u, v, w, gamma factors.$'
    6876 **
    6877 ** for strict FORTRAN77 compliance
    6878 *        write(scrtch, *, err=990) dummy
    6879 *        rewind(scrtch, err=990)
    6880 *        read(scrtch, *, err=999) pv_u, pv_v, pv_w, pv_gamma
    6881 *        rewind(scrtch, err=990)
    6882 ** for more lenient, efficient, compilers
    6883 **        read(dummy, *, err=999) pv_u, pv_v, pv_w, pv_gamma
    6884 **
    6885 *        if(pv_w.lt.ZERO) then
    6886 *          write(op,102)
    6887 *     |        'WARNING: pseudo-Voigt w factor is less than zero.'
    6888 *          write(op,102)
    6889 *     |        '  This may cause the pseudo-Voigt function to fail.'
    6890 *        endif
    6891 **
    6892 *        if(abs(pv_u).lt.eps7 .and. abs(pv_v).lt.eps7 .and.
    6893 *     |                         abs(pv_w).lt.eps7) then
    6894 *          write(op,102) 'Pseudo-Voigt u, v, w factors are zero.'
    6895 *          write(op,102) 'No instrumental broadening will be simulated.'
    6896 *          blurring = NONE
    6897 *        endif
    6898 **
    6899 *        if(pv_gamma.lt.ZERO) then
    6900 *          messge =
    6901 *     |   'Pseudo-Voigt gamma factor is -ve (must be between 0 and 1).$'
    6902 *          goto 999
    6903 *        endif
    6904 **
    6905 *        if(pv_gamma.gt.ONE) then
    6906 *          messge = 'Pseudo-Voigt gamma factor is greater than 1.$'
    6907 *          goto 999
    6908 *        endif
    6909 **
    6910 ** are the pseudo-Voigt parameters equivalent to a Gaussian with
    6911 ** constant standard deviation?
    6912 *        if(abs(pv_u).le.eps7 .and. abs(pv_v).le.eps7 .and.
    6913 *     |     pv_w.gt.ZERO .and. pv_gamma.le.eps4) then
    6914 *          FWHM = sqrt(pv_w)
    6915 *          write(op,102)
    6916 *     |        'Pseudo-Voigt factors are equivalent to a Gaussian'
    6917 *          write(op,103) ' with FWHM ', FWHM
    6918 *          blurring = GAUSS
    6919 *        endif
    6920 **
    6921 ** are the pseudo-Voigt parameters equivalent to a Lorentzian with
    6922 ** constant half-width?
    6923 *        if(abs(pv_u).le.eps7 .and. abs(pv_v).le.eps7 .and.
    6924 *     |     pv_w.gt.ZERO .and. abs(pv_gamma - ONE).le.eps4) then
    6925 *          FWHM = sqrt(pv_w)
    6926 *          write(op,102)
    6927 *     |        'Pseudo-Voigt factors are equivalent to a Lorentzian'
    6928 *          write(op,103) ' with FWHM ', FWHM
    6929 *          blurring = LORENZ
    6930 *        endif
    6931 **
    6932 *      endif
    6933 **
    6934 *      RDBLUR = .true.
    6935 *      return
    6936 *  990 messge = 'Problems using scratch file in RDBLUR.$'
    6937 *  999 call WRTLNE(line)
    6938 *      return
    6939 *  100 format(1x, 'ERROR: ', a)
    6940 *  101 format(1x, 2a)
    6941 *  102 format(1x, a)
    6942 *  103 format(1x, a, g12.5)
    6943 *      end
    6944 **
    6945 ** ______________________________________________________________________
    6946 ** Title: RDCELL
    6947 ** Authors: MMJT
    6948 ** Date: 4 Oct 1989
    6949 ** Description: This function reads the layer mesh dimensions from
    6950 ** the second data line of the data file. It checks that the header
    6951 ** 'STRUCTURAL' is present.
    6952 **
    6953 **      ARGUMENTS:
    6954 **            unit_no  -  Logical unit number that the data file is to
    6955 **                        be read from. (input).
    6956 **            messge   -  A short message indicating what went wrong
    6957 **                        (if anything) during the datafile read.
    6958 **                        messge is terminated by a token '$'. (output).
    6959 **
    6960 **      COMMON VARIABLES:
    6961 **        modifies:  cell_a, cell_b, cell_c, cell_gamma
    6962 **
    6963 **      RDCELL returns logical .true. if acceptable values
    6964 **      for the layer mesh dimensions were read.
    6965 ** ______________________________________________________________________
    6966 **
    6967 *      logical function RDCELL(unit_no, messge)
    6968 *      include 'DIFFaX.par'
    6969 *      include 'DIFFaX.inc'
    6970 **
    6971 *      integer*4 unit_no
    6972 *      character*(*) messge
    6973 **
    6974 *      character*200 line
    6975 **
    6976 ** external subroutines (Some compilers need them declared external)
    6977 **      external GETLNE, TOUPPR, WRTLNE
    6978 **
    6979 *      RDCELL = .false.
    6980 **
    6981 ** first check that 'STRUCTURAL' header is present
    6982 *      messge = '''STRUCTURAL'' header could not be read.$'
    6983 *      call GETLNE(unit_no, line, *999)
    6984 *      call TOUPPR(line)
    6985 *      messge = '''STRUCTURAL'' header not found.$'
    6986 *      if(index(line,'STRUCTURAL ').eq.0) goto 999
    6987 **
    6988 *      messge = 'a, b, c, and gamma are improperly specified.$'
    6989 *      call GETLNE(unit_no, line, *999)
    6990 **
    6991 ** for strict FORTRAN77 compliance
    6992 *      write(scrtch, *, err=990) line
    6993 *      rewind(scrtch, err=990)
    6994 *      read(scrtch, *, err=999) cell_a, cell_b, cell_c, cell_gamma
    6995 *      rewind(scrtch, err=990)
    6996 ** for more lenient, efficient, compilers
    6997 **      read(line, *, err=999) cell_a, cell_b, cell_c, cell_gamma
    6998 **
    6999 *      if(cell_a.le.ZERO) then
    7000 *        messge = 'Negative (or zero) value for cell_a.$'
    7001 *        goto 999
    7002 *      endif
    7003 **
    7004 *      if(cell_b.le.ZERO) then
    7005 *        messge = 'Negative (or zero) value for cell_b.$'
    7006 *        goto 999
    7007 *      endif
    7008 **
    7009 *      if(cell_c.le.ZERO) then
    7010 *        messge = 'Negative (or zero) value for cell_c.$'
    7011 *        goto 999
    7012 *      endif
    7013 **
    7014 *      if(cell_gamma.le.ZERO) then
    7015 *        messge = 'Negative (or zero) value for cell_gamma.$'
    7016 *        goto 999
    7017 *      endif
    7018 **
    7019 *      if(cell_gamma.ge.ONE_EIGHTY) then
    7020 *        messge = 'Cell_gamma is greater than 180 degrees.$'
    7021 *        goto 999
    7022 *      endif
    7023 **
    7024 *      cell_gamma = cell_gamma * DEG2RAD
    7025 *      RDCELL = .true.
    7026 *      return
    7027 *  990 messge = 'Problems using scratch file in RDCELL.$'
    7028 *  999 call WRTLNE(line)
    7029 *      return
    7030 *      end
    7031 **
    7032 ** ______________________________________________________________________
    7033 ** Title: RDFILE
    7034 ** Author: MMJT
    7035 ** Date: 3 Oct 1989; 3 Mar 1995
    7036 ** Description: This function reads in the data from the user-defined
    7037 ** data file. The data file is in Naur-Backus form, and the user may
    7038 ** insert '{' and '}' delimiters to comment the data. The data file
    7039 ** is not read in its raw form. The subroutine GETLNE is used
    7040 ** to read each line, stripped of its comments, and left-justified.
    7041 **
    7042 **      ARGUMENTS:
    7043 **           fname  -  The name of the input data file to read. (input).
    7044 **
    7045 **      RDFILE returns logical .true. if the file read went well.
    7046 ** ______________________________________________________________________
    7047 **
    7048 *      logical function RDFILE(fname)
    7049 *      include 'DIFFaX.par'
    7050 **     save
    7051 **
    7052 *      character*(*) fname
    7053 **
    7054 *      character*200 messge
    7055 *      logical goodfile
    7056 *      logical RDRADN, RDWAVL, RDBLUR, RDCELL, RDPTGP
    7057 *      logical RDNLAY, RDWDTH, RDLAYR, RDSTAK, RDTRNS
    7058 *      integer*4 LENGTH, unit_no
    7059 **
    7060 ** external functions
    7061 *      external RDRADN, RDWAVL, RDBLUR, RDCELL, RDPTGP
    7062 *      external RDNLAY, RDWDTH, RDLAYR, RDSTAK, RDTRNS
    7063 *      external LENGTH
    7064 **
    7065 *      RDFILE = .false.
    7066 *      unit_no = df
    7067 *      write(op,100) 'Reading ''', fname(1:LENGTH(fname)),''''
    7068 **
    7069 ** now read the file, one datum at a time.
    7070 *                   goodfile = RDRADN(unit_no, messge)
    7071 *      if(goodfile) goodfile = RDWAVL(unit_no, messge)
    7072 *      if(goodfile) goodfile = RDBLUR(unit_no, messge)
    7073 *      if(goodfile) goodfile = RDCELL(unit_no, messge)
    7074 *      if(goodfile) goodfile = RDPTGP(unit_no, messge)
    7075 *      if(goodfile) goodfile = RDNLAY(unit_no, messge)
    7076 *      if(goodfile) goodfile = RDWDTH(unit_no, messge)
    7077 *      if(goodfile) goodfile = RDLAYR(unit_no, messge)
    7078 *      if(goodfile) goodfile = RDSTAK(unit_no, messge)
    7079 *      if(goodfile) goodfile = RDTRNS(unit_no, messge)
    7080 **
    7081 *      if(goodfile) then
    7082 *        write(op,100) '''', fname(1:LENGTH(fname)),''' read in.'
    7083 *      else
    7084 *        write(op,200) messge(1:index(messge,'$') - 1)
    7085 *      endif
    7086 **
    7087 *      write(messge,101)
    7088 *     |      'Problems closing ''', fname(1:LENGTH(fname)),'''$'
    7089 *      close(df,err=900)
    7090 ** if all went well, set RDFILE to .true.
    7091 *      if(goodfile) RDFILE = .true.
    7092 *      return
    7093 **
    7094 *  900 write(op,200) messge(1:index(messge,'$') - 1)
    7095 *      return
    7096 *  100 format(1x, 3a)
    7097 *  101 format(3a)
    7098 *  200 format(1x, 'ERROR: ', a)
    7099 *      end
    7100 **
    7101 ** ______________________________________________________________________
    7102 ** Title: RDLAYR
    7103 ** Authors: MMJT
    7104 ** Date: 7 June 1990
    7105 ** Description: This function reads in the layer atomic coordinates.
    7106 ** Each set of data must be prefaced with a layer number, and a
    7107 ** statement of whether or not the layer is centrosymmetric. If the
    7108 ** layer is centrosymmetric, only the asymmetric coordinates need
    7109 ** be entered (with appropriate fractional occupancies for special
    7110 ** positions.
    7111 ** If a layer has the same structure as a previously defined layer
    7112 ** the layer can be defined as 'LAYER i = j', where j is a previously
    7113 ** defined layer.
    7114 ** The format assumed for atom coordinates is:
    7115 **
    7116 **    'atom name', number, x, y, z, Debye-Waller factor, occupancy
    7117 **
    7118 ** 'atom name' must occupy 4 characters, eg. 'Si4+', 'O 1-' etc...
    7119 ** The allowed atom names can be found in the file 'data.sfc'.
    7120 **
    7121 **      ARGUMENTS:
    7122 **            unit_no  -  Logical unit number that the data file is to
    7123 **                        be read from. (input).
    7124 **            messge   -  A short message indicating what went wrong
    7125 **                        (if anything) during the datafile read.
    7126 **                        messge is terminated by a token '$'. (output).
    7127 **
    7128 **      COMMON VARIABLES:
    7129 **            uses:  n_layers, l_actual, CENTRO, NONE
    7130 **
    7131 **        modifies:  l_symmetry, a_name, a_number, a_pos, a_B, a_occup,
    7132 **                   l_n_atoms, n_actual
    7133 **
    7134 **      RDLAYR returns logical .true. if all the layer data
    7135 **      was read successfully.
    7136 ** ______________________________________________________________________
    7137 **
    7138 *      logical function RDLAYR(unit_no, messge)
    7139 *      include 'DIFFaX.par'
    7140 *      include 'DIFFaX.inc'
    7141 **
    7142 *      integer*4 unit_no
    7143 *      character*(*) messge
    7144 **
    7145 *      character*200 line, tmpline, arg
    7146 *      character*80 list(2)
    7147 *      integer*4 CHOICE, LENGTH, NXTARG, m, n, i, j, i2, iflag
    7148 *      real*8 RDNMBR, tmp
    7149 *      logical*4 ok
    7150 **
    7151 ** external functions
    7152 *      external CHOICE, LENGTH, NXTARG, RDNMBR
    7153 ** external subroutines (Some compilers need them declared external)
    7154 **      external GETLNE, TOUPPR, WRTLNE
    7155 **
    7156 *      RDLAYR = .false.
    7157 **
    7158 *      do 5 i = 1, n_layers
    7159 *        high_atom(i) = ZERO
    7160 *        low_atom(i)  = ZERO
    7161 *    5 continue
    7162 **
    7163 *      m = 0
    7164 *      messge = 'Unexpected EOF before LAYER 1.$'
    7165 *      call GETLNE(unit_no, line, *999)
    7166 *      call TOUPPR(line)
    7167 **
    7168 *      do 30 i = 1, n_layers
    7169 *        write(messge,100) 'LAYER header not seen for LAYER ',i,'.$'
    7170 *        if(index(line, 'LAYER').ne.1) goto 999
    7171 ** read layer number
    7172 *        write(messge,100) 'LAYER ',i,' number sequence invalid.$'
    7173 *        write(tmpline, '(a)', err=999)
    7174 *     |       line(LENGTH('LAYER')+1:len(line))
    7175 ** for strict FORTRAN77 compliance
    7176 *        write(scrtch, *, err=990) tmpline
    7177 *        rewind(scrtch, err=990)
    7178 *        read(scrtch, *, err=999) j
    7179 *        rewind(scrtch, err=990)
    7180 ** for more lenient, efficient, compilers
    7181 **        read(tmpline, *, err=999) j
    7182 **
    7183 ** layers must be entered in sequence
    7184 *        if(j.ne.i) goto 999
    7185 **
    7186 *        j = index(line, '=')
    7187 ** see if this layer is equivalenced to another
    7188 *        if(j.ne.0) then
    7189 **
    7190 *          write(messge,100)
    7191 *     |           'Invalid equivalence number for LAYER ',i,'.$'
    7192 *          write(tmpline, '(a)', err=999) line(j+1:len(line))
    7193 ** for strict FORTRAN77 compliance
    7194 *          write(scrtch, *, err=990) tmpline
    7195 *          rewind(scrtch, err=990)
    7196 *          read(scrtch, *, err=999) i2
    7197 *          rewind(scrtch, err=990)
    7198 ** for more lenient, efficient, compilers
    7199 **        read(tmpline, *, err=999) i2
    7200 **
    7201 ** layer i2 must be defined already
    7202 *          if(i2.ge.i .or. i2.lt.1) goto 999
    7203 *          l_actual(i) = l_actual(i2)
    7204 **
    7205 *          write(messge,100) 'Unexpected EOF in LAYER ',i,'.$'
    7206 *          call GETLNE(unit_no, line, *999)
    7207 *          call TOUPPR(line)
    7208 **
    7209 *        else
    7210 *          m = m + 1
    7211 *          l_actual(i) = m
    7212 **
    7213 *          write(messge,100)
    7214 *     |           'No symmetry statement for LAYER ',i,'.$'
    7215 *          call GETLNE(unit_no, line, *999)
    7216 *          call TOUPPR(line)
    7217 **
    7218 *          list(1) = 'NONE '
    7219 *          list(2) = 'CENTROSYMMETRIC '
    7220 *          iflag = CHOICE(line, list, 2)
    7221 **
    7222 *          write(messge,100)
    7223 *     |           'Invalid symmetry statement for LAYER ',i,'.$'
    7224 *          if(iflag.lt.1 .or. iflag.gt.2) goto 999
    7225 **
    7226 *          if(iflag.eq.1) then
    7227 *            l_symmetry(m) = NONE
    7228 *          else
    7229 *            l_symmetry(m) = CENTRO
    7230 *          endif
    7231 **
    7232 *          write(messge,100) 'No atoms specified for LAYER ',i,'.$'
    7233 *          call GETLNE(unit_no, line, *999)
    7234 *          if(index(line,'LAYER').eq.1) goto 999
    7235 **
    7236 *          j = 0
    7237 *   20       j = j + 1
    7238 *            if(j.gt.MAX_A) then
    7239 *              write(op,101) 'ERROR: Too many atoms on LAYER ',i,'.'
    7240 *              write(op,102) '       Maximum allowed = ', MAX_A
    7241 *              messge = ' $'
    7242 *              goto 999
    7243 *            endif
    7244 **
    7245 *            a_name(j,m) = line(1:4)
    7246 **
    7247 *            write(messge,103)
    7248 *     |        'Problems reading data for atom ''', a_name(j,m), '''.$'
    7249 *            write(tmpline, '(a)', err=999) line(5:len(line))
    7250 ** read atom identifier
    7251 *            n = NXTARG(tmpline,arg)
    7252 *            if(n.le.0) goto 999
    7253 *            read(arg,'(i3)') a_number(j,m)
    7254 **
    7255 ** read x_relative position
    7256 *            n = NXTARG(tmpline,arg)
    7257 *            if(n.le.0) goto 999
    7258 *            a_pos(1,j,m) = RDNMBR(arg,ok)
    7259 *            if(.not.ok) goto 999
    7260 **
    7261 ** read y_relative position
    7262 *            n = NXTARG(tmpline,arg)
    7263 *            if(n.le.0) goto 999
    7264 *            a_pos(2,j,m) = RDNMBR(arg,ok)
    7265 *            if(.not.ok) goto 999
    7266 **
    7267 ** read z_relative position
    7268 *            n = NXTARG(tmpline,arg)
    7269 *            if(n.le.0) goto 999
    7270 *            a_pos(3,j,m) = RDNMBR(arg,ok)
    7271 *            if(.not.ok) goto 999
    7272 **
    7273 ** read Debye-Waller factor
    7274 *            n = NXTARG(tmpline,arg)
    7275 *            if(n.le.0) goto 999
    7276 *            a_B(j,m) = RDNMBR(arg,ok)
    7277 *            if(.not.ok) goto 999
    7278 **
    7279 ** read occupancy
    7280 *            n = NXTARG(tmpline,arg)
    7281 *            if(n.le.0) goto 999
    7282 *            a_occup(j,m) = RDNMBR(arg,ok)
    7283 *            if(.not.ok) goto 999
    7284 **
    7285 ** Check values
    7286 *            if(a_B(j,m).lt.ZERO) then
    7287 *              messge = 'Negative Debye-Waller factor.$'
    7288 *              goto 999
    7289 *            endif
    7290 **
    7291 *            if(a_occup(j,m).lt.ZERO) then
    7292 *              messge = 'Negative atom occupancy.$'
    7293 *              goto 999
    7294 *            endif
    7295 **
    7296 *            if(a_occup(j,m).gt.ONE) then
    7297 *              messge = 'Occupancy greater than 1.$'
    7298 *              goto 999
    7299 *            endif
    7300 **
    7301 ** get extreme upper and lower atom positions for error-checking later on
    7302 *            tmp = a_pos(3,j,m)
    7303 *            if(tmp.gt.high_atom(m)) high_atom(m) = tmp
    7304 *            if(tmp.lt.low_atom(m))   low_atom(m) = tmp
    7305 **
    7306 *            write(messge,103)'Unexpected EOF after atom ',line(1:4),'$'
    7307 *            call GETLNE(unit_no, line, *999)
    7308 ** see if all atomic specifications for this layer have been read in
    7309 ** but first we must convert to uppercase for the test. If there are
    7310 ** more atoms, then we must convert 'line' to its original format.
    7311 *            tmpline = line
    7312 *            call TOUPPR(line)
    7313 ** check data is in correct sequence
    7314 *            messge = 'STACKING data is missing.$'
    7315 *            if(index(line,'PARAMETERS').eq.1) goto 999
    7316 *            if(index(line,'LAYER').ne.1
    7317 *     |             .and. index(line,'STACKING').ne.1) then
    7318 *              line = tmpline
    7319 *              goto 20
    7320 *            endif
    7321 *            l_n_atoms(m) = j
    7322 *        endif
    7323 *   30 continue
    7324 *      n_actual = m
    7325 **
    7326 ** make sure we have genuine lowest and highest atoms in each layer.
    7327 *      do 40 i = 1, n_actual
    7328 *        if(l_symmetry(i).eq.CENTRO) then
    7329 *          if(-low_atom(i).gt.high_atom(i)) then
    7330 *            high_atom(i) = -low_atom(i)
    7331 *          else if(-low_atom(i).lt.high_atom(i)) then
    7332 *            low_atom(i) = -high_atom(i)
    7333 *          endif
    7334 *        endif
    7335 *   40 continue
    7336 **
    7337 ** If we have reached here, then we have read one line too many.
    7338 *      write(messge,100) 'Problems ''backspacing'' unit ',unit_no, '.$'
    7339 *      backspace(unit = unit_no, err = 980)
    7340 **
    7341 *      RDLAYR = .true.
    7342 *  980 return
    7343 *  990 messge = 'Problems using scratch file in RDLAYR.$'
    7344 *  999 call WRTLNE(line)
    7345 *      return
    7346 *  100 format(a, i2, a)
    7347 *  101 format(1x, a, i2, a)
    7348 *  102 format(1x, a, i4)
    7349 *  103 format(3a)
    7350 *      end
    7351 **
    7352 ** ______________________________________________________________________
    7353 ** Title: RDNLAY
    7354 ** Authors: MMJT
    7355 ** Date: 4 Oct 1989
    7356 ** Description: This function reads the first line of the data file
    7357 ** to extract the number of layers in the structure.
    7358 **
    7359 **      ARGUMENTS:
    7360 **            unit_no  -  Logical unit number that the data file is to
    7361 **                        be read from. (input).
    7362 **            messge   -  A short message indicating what went wrong
    7363 **                        (if anything) during the datafile read.
    7364 **                        messge is terminated by a token '$'. (output).
    7365 **
    7366 **      COMMON VARIABLES:
    7367 **        modifies:  n_layers
    7368 **
    7369 **      RDNLAY returns logical .true. if an acceptable value
    7370 **      for the number of layers was read.
    7371 ** ______________________________________________________________________
    7372 **
    7373 *      logical function RDNLAY(unit_no, messge)
    7374 *      include 'DIFFaX.par'
    7375 *      include 'DIFFaX.inc'
    7376 **
    7377 *      integer*4 unit_no
    7378 *      character*(*) messge
    7379 **
    7380 *      character*200 line
    7381 **
    7382 ** external subroutines (Some compilers need them declared external)
    7383 **      external GETLNE, WRTLNE
    7384 **
    7385 *      RDNLAY = .false.
    7386 **
    7387 *      messge = 'No lines in file.$'
    7388 *      call GETLNE(unit_no, line, *999)
    7389 **
    7390 *      messge = 'Problems reading number of layers.$'
    7391 ** for strict FORTRAN77 compliance
    7392 *      write(scrtch, *, err=990) line
    7393 *      rewind(scrtch, err=990)
    7394 *      read(scrtch, *, err=999) n_layers
    7395 *      rewind(scrtch, err=990)
    7396 ** for more lenient, efficient, compilers
    7397 **      read(line, *, err=999) n_layers
    7398 **
    7399 *      if(n_layers.eq.0) then
    7400 *        messge = 'Zero number of layers entered.$'
    7401 *        goto 999
    7402 *      else if(n_layers.lt.0) then
    7403 *        messge = 'Negative number of layers entered.$'
    7404 *        goto 999
    7405 *      else if(n_layers.gt.MAX_L) then
    7406 *        write(messge,100)
    7407 *     |  'Too many layers: number should not exceed ',MAX_L,'.$'
    7408 *        goto 999
    7409 *      endif
    7410 **
    7411 *      RDNLAY = .true.
    7412 *      return
    7413 *  990 messge = 'Problems using scratch file in RDNLAY.$'
    7414 *  999 call WRTLNE(line)
    7415 *      return
    7416 *  100 format(a, i2, a)
    7417 *      end
    7418 **
    7419 ** ______________________________________________________________________
    7420 ** Title: RDNMBR
    7421 ** Author: MMJT
    7422 ** Date: 21 January 2005
    7423 ** Description: Reads the character string 'numberstring' to extract a
    7424 ** floating point number. The character string can be either a true
    7425 ** number, such as '0.3333', or can be expressed as a ratio, '1/3'.
    7426 **
    7427 **      ARGUMENTS:
    7428 **            arg   -  Character string containing numeric data.(input).
    7429 **            ok    -  True of all went well.(output).
    7430 **
    7431 **      Returns the floating point number
    7432 ** ______________________________________________________________________
    7433 **
    7434 *      real*8 function RDNMBR(numberstring,ok)
    7435 *      include 'DIFFaX.par'
    7436 *      save
    7437 **
    7438 *      character*(*) numberstring
    7439 *      logical ok
    7440 **
    7441 *      character*200 arg
    7442 *      integer*4 numerator, denominator, i, j, lin_len
    7443 *      integer*4 NXTARG
    7444 **
    7445 ** external functions
    7446 *      external NXTARG
    7447 **
    7448 *      ok = .true.
    7449 *      RDNMBR = ZERO
    7450 **
    7451 *      lin_len = len(numberstring)
    7452 *      j = index(numberstring, '/')
    7453 *      if(j.eq.0) then
    7454 ** It was already a number
    7455 *        read(numberstring, *) RDNMBR
    7456 *      else if(j.lt.0 .or. j.ge.lin_len) then
    7457 ** An error happened
    7458 *        goto 999
    7459 *      else
    7460 ** It was expressed as a ratio. Blank out the slash, '/'.
    7461 *        numberstring(j:j) = ' '
    7462 ** 'numberstring' now contains two arguments, numerator and denominator.
    7463 *        i = NXTARG(numberstring,arg)
    7464 *        if(i.le.0) goto 999
    7465 *        read(arg, *) numerator
    7466 *        i = NXTARG(numberstring,arg)
    7467 *        if(i.le.0) goto 999
    7468 *        read(arg, *) denominator
    7469 *        if(denominator.le.ZERO) goto 999
    7470 *        RDNMBR = dble(numerator) / dble(denominator)
    7471 *      endif
    7472 **
    7473 *  100 return
    7474 *  999 ok = .false.
    7475 *      return
    7476 *      end
    7477 **
    7478 ** ______________________________________________________________________
    7479 ** Title: RDPTGP
    7480 ** Author: MMJT
    7481 ** Date: 4 Oct 1989
    7482 ** Description: This function reads the user's estimate of the
    7483 ** point group symmetry of the diffraction data.
    7484 **
    7485 **      ARGUMENTS:
    7486 **            unit_no  -  Logical unit number that the data file is to
    7487 **                        be read from. (input).
    7488 **            messge   -  A short message indicating what went wrong
    7489 **                        (if anything) during the datafile read.
    7490 **                        messge is terminated by a token '$'. (output).
    7491 **
    7492 **      COMMON VARIABLES:
    7493 **        modifies:  pnt_grp, SymGrpNo, tolerance
    7494 **
    7495 **      RDPTGP returns logical .true. if a legal point group
    7496 **      was entered, or if 'unknown' (with a search tolerance)
    7497 **      was specified.
    7498 **
    7499 **      Legal entries are:
    7500 **                     '-1'
    7501 **                     '2/M(1)'    (diad along streaks)
    7502 **                     '2/M(2)'    (mirror contains streaks)
    7503 **                     'MMM'
    7504 **                     '-3'
    7505 **                     '-3M'
    7506 **                     '4/M'
    7507 **                     '4/MMM'
    7508 **                     '6/M'
    7509 **                     '6/MMM'
    7510 **                     'AXIAL'     (for integration along 00l only)
    7511 **                     'UNKNOWN'   (followed by an optional tolerance)
    7512 ** ______________________________________________________________________
    7513 **
    7514 *      logical function RDPTGP(unit_no, messge)
    7515 *      include 'DIFFaX.par'
    7516 *      include 'DIFFaX.inc'
    7517 **
    7518 *      integer*4 unit_no
    7519 *      character*(*) messge
    7520 **
    7521 *      character*200 line, oldline
    7522 *      character*80 list(12)
    7523 *      integer*4 CHOICE, LENGTH, PRUNE, i, iflag
    7524 **
    7525 ** external functions
    7526 *      external CHOICE, LENGTH, PRUNE
    7527 ** external subroutines (Some compilers need them declared external)
    7528 **      external GETLNE, TOUPPR, WRTLNE
    7529 **
    7530 *      RDPTGP = .false.
    7531 **
    7532 *      messge = 'Illegal diffraction point group symmetry.$'
    7533 *      call GETLNE(unit_no, line, *999)
    7534 *      call TOUPPR(line)
    7535 **
    7536 *      list(1)  = '-1 '
    7537 *      list(2)  = '2/M(1) '
    7538 *      list(3)  = '2/M(2) '
    7539 *      list(4)  = 'MMM '
    7540 *      list(5)  = '-3 '
    7541 *      list(6)  = '-3M '
    7542 *      list(7)  = '4/M '
    7543 *      list(8)  = '4/MMM '
    7544 *      list(9)  = '6/M '
    7545 *      list(10) = '6/MMM '
    7546 *      list(11) = 'AXIAL '
    7547 *      list(12) = 'UNKNOWN '
    7548 *      iflag = CHOICE(line, list, 12)
    7549 *      if(iflag.lt.1 .or. iflag.gt.12) then
    7550 *        write(op,100) 'Unrecognized symmetry type. Legal types are:'
    7551 *        do 10 i = 1, 12
    7552 *          write(op,101) '               ',list(i)
    7553 *   10   continue
    7554 *        goto 999
    7555 *      endif
    7556 **
    7557 *      pnt_grp = list(iflag)
    7558 **
    7559 ** set default tolerance on intensity variability
    7560 *      tolerance = ONE
    7561 **
    7562 *      if(iflag.ge.1 .and. iflag.le.11) then
    7563 *        SymGrpNo = iflag
    7564 *      else if(iflag.eq.12) then
    7565 ** a value of UNKNOWN = -1 alerts OPTIMZ that user entered 'UNKNOWN'
    7566 *        SymGrpNo = UNKNOWN
    7567 ** Strip off the symmetry point group text string.
    7568 *        oldline = line
    7569 *        write(line,'(a)') oldline(LENGTH(list(iflag))+1:len(oldline))
    7570 ** If the remainder of the line is blank, no tolerance was specified.
    7571 *        i = PRUNE(line)
    7572 *        if(i.le.0) then
    7573 *          messge = 'Problems reading symmetry tolerance parameter.$'
    7574 *          goto 999
    7575 *        else if(i.gt.1) then
    7576 ** There was something on the line, try to read tolerance
    7577 *          messge = 'Problems reading symmetry tolerance parameter.$'
    7578 ** for strict FORTRAN77 compliance
    7579 *          write(scrtch, *, err=990) line
    7580 *          rewind(scrtch, err=990)
    7581 *          read(scrtch, *, err=999) tolerance
    7582 *          rewind(scrtch, err=990)
    7583 ** for more lenient, efficient, compilers
    7584 **          read(line, *, err=999) tolerance
    7585 **
    7586 *          if(tolerance.le.ZERO) then
    7587 *            messge = 'Negative (or zero) value for tolerance.$'
    7588 *            goto 999
    7589 *          endif
    7590 **
    7591 *          if(tolerance.lt.eps2) then
    7592 *            tolerance = eps2
    7593 *            write(op,102)
    7594 *            write(op,103) tolerance
    7595 *          endif
    7596 **
    7597 *        endif
    7598 *      endif
    7599 **
    7600 ** convert from a percentage
    7601 *      tolerance = tolerance * eps2
    7602 **
    7603 *      RDPTGP = .true.
    7604 *      return
    7605 *  990 messge = 'Problems using scratch file in RDPTGP.$'
    7606 *  999 call WRTLNE(line)
    7607 *      return
    7608 *  100 format(1x, 'ERROR: ', a)
    7609 *  101 format(1x, 2a)
    7610 *  102 format(1x, 'WARNING: Tolerance for symmetry tests is too small.')
    7611 *  103 format(1x, '         Tolerance is reset to ', g12.5, ' percent')
    7612 *      end
    7613 **
    7614 ** ______________________________________________________________________
    7615 ** Title: RDRADN
    7616 ** Authors: MMJT
    7617 ** Date: 15 Feb 1990
    7618 ** Description: This function reads the radiation type. The choices are
    7619 ** X-RAY, NEUTRON and ELECTRON. Checks that the keword 'INSTRUMENTAL'
    7620 ** is present.
    7621 **
    7622 **      ARGUMENTS:
    7623 **            unit_no  -  Logical unit number that the data file is to
    7624 **                        be read from. (input).
    7625 **            messge   -  A short message indicating what went wrong
    7626 **                        (if anything) during the datafile read.
    7627 **                        messge is terminated by a token '$'. (output).
    7628 **
    7629 **      COMMON VARIABLES:
    7630 **            uses:  ELECTN, NEUTRN, X_RAY
    7631 **
    7632 **        modifies:  rad_type
    7633 **
    7634 **      RDRADN returns logical .true. if an acceptable
    7635 **      radiation type was read.
    7636 ** ______________________________________________________________________
    7637 **
    7638 *      logical function RDRADN(unit_no, messge)
    7639 *      include 'DIFFaX.par'
    7640 *      include 'DIFFaX.inc'
    7641 **
    7642 *      integer*4 unit_no
    7643 *      character*(*) messge
    7644 **
    7645 *      character*200 line
    7646 *      character*80 list(3)
    7647 *      integer*4 CHOICE, iflag
    7648 **
    7649 ** external function
    7650 *      external CHOICE
    7651 ** external subroutines (Some compilers need them declared external)
    7652 **      external GETLNE, TOUPPR, WRTLNE
    7653 **
    7654 *      RDRADN = .false.
    7655 **
    7656 ** first check that 'INSTRUMENTAL' header is present
    7657 *      messge = '''INSTRUMENTAL'' header could not be read.$'
    7658 *      call GETLNE(unit_no, line, *999)
    7659 *      call TOUPPR(line)
    7660 *      messge = '''INSTRUMENTAL'' header not found.$'
    7661 *      iflag  =index(line,'INSTRUMENTAL ')
    7662 *      if(iflag.eq.0) goto 999
    7663 **
    7664 *      messge = 'radiation type not specified.$'
    7665 *      call GETLNE(unit_no, line, *999)
    7666 *      call TOUPPR(line)
    7667 *      list(1) = 'X-RAY '
    7668 *      list(2) = 'NEUTRON '
    7669 *      list(3) = 'ELECTRON '
    7670 *      iflag = CHOICE(line, list, 3)
    7671 **
    7672 *      if(iflag.eq.1) then
    7673 *        rad_type = X_RAY
    7674 *      else if(iflag.eq.2) then
    7675 *        rad_type = NEUTRN
    7676 *      else if(iflag.eq.3) then
    7677 *        rad_type = ELECTN
    7678 *      else
    7679 *        messge ='Invalid radiation type (X-RAY, NEUTRON or ELECTRON).$'
    7680 *        goto 999
    7681 *      endif
    7682 **
    7683 *      RDRADN = .true.
    7684 *      return
    7685 *  999 call WRTLNE(line)
    7686 *      return
    7687 *      end
    7688 **
    7689 ** ______________________________________________________________________
    7690 ** Title: RDRNGE
    7691 ** Authors: MMJT
    7692 ** Date: 15 Feb 1990
    7693 ** Description: This function reads the angular range for the spectrum.
    7694 **
    7695 **      ARGUMENTS:
    7696 **            unit_no  -  Logical unit number that the data file is to
    7697 **                        be read from. (input).
    7698 **
    7699 **      COMMON VARIABLES:
    7700 **            uses:  DEG2RAD
    7701 **
    7702 **        modifies:  th2_min, th2_max, d_theta
    7703 **
    7704 **      RDRNGE returns logical .true. if an acceptable angular
    7705 **      range was read.
    7706 ** ______________________________________________________________________
    7707 **
    7708 *      logical function RDRNGE()
    7709 *      include 'DIFFaX.par'
    7710 *      include 'DIFFaX.inc'
    7711 **
    7712 ** external subroutines (Some compilers need them declared external)
    7713 **      external GETLNE, WRTLNE
    7714 **
    7715 *      RDRNGE = .false.
    7716 **
    7717 *  123 write(op,100) 'Enter angular range:'
    7718 *      write(op,100) '  2theta min, 2theta max, 2theta increment.'
    7719 *      read(cntrl,*,err=123) th2_min, th2_max, d_theta
    7720 *      if(CFile) write(op,105) th2_min, th2_max, d_theta
    7721 **
    7722 *      if(th2_min.lt.ZERO) then
    7723 *        write(op,110) '2theta min is negative.'
    7724 *        goto 999
    7725 *      endif
    7726 **
    7727 *      if(th2_min.ge.ONE_EIGHTY) then
    7728 *        write(op,110) '2theta min exceeds 180 degrees.'
    7729 *        goto 999
    7730 *      endif
    7731 **
    7732 *      if(th2_max.le.ZERO) then
    7733 *        write(op,110) 'Negative (or zero) value for 2theta min.'
    7734 *        goto 999
    7735 *      endif
    7736 **
    7737 *      if(th2_max.gt.ONE_EIGHTY) then
    7738 *        write(op,110) '2theta max exceeds 180 degrees.'
    7739 *        goto 999
    7740 *      endif
    7741 **
    7742 ** if th2_max = 180, reduce it slightly to keep us out of trouble
    7743 *      if(th2_max.eq.ONE_EIGHTY) th2_max = th2_max - eps4
    7744 **
    7745 *      if(d_theta.le.ZERO) then
    7746 *        write(op,110) 'Negative (or zero) value for 2theta increment.'
    7747 *        goto 999
    7748 *      endif
    7749 **
    7750 *      if(th2_min.ge.th2_max) then
    7751 *        write(op,110) '2theta min is larger than 2theta max.'
    7752 *        goto 999
    7753 *      endif
    7754 **
    7755 *      if((th2_max - th2_min).lt.d_theta) then
    7756 *        write(op,110)
    7757 *     |  '2theta increment is larger than 2theta max - 2theta min.'
    7758 *        goto 999
    7759 *      endif
    7760 **
    7761 *      th2_min = th2_min * DEG2RAD
    7762 *      th2_max = th2_max * DEG2RAD
    7763 *      d_theta = HALF * DEG2RAD * d_theta
    7764 **
    7765 *      RDRNGE = .true.
    7766 *      return
    7767 **
    7768 *  999 return
    7769 *  100 format(1x, a)
    7770 *  105 format(1x, g11.4, 2(3x, g11.4))
    7771 *  110 format(1x, 'ERROR: ', a)
    7772 *      end
    7773 **
    7774 ** ______________________________________________________________________
    7775 ** Title: RDSTAK
    7776 ** Author: MMJT
    7777 ** Date: 4 Oct 1989
    7778 ** Description: This function reads the STACKING input of the datafile.
    7779 ** The stacking is either 'EXPLICIT', whereupon DIFFaX expects the
    7780 ** user to enter an explicit sequence of layers, or 'RECURSIVE'. If
    7781 ** 'EXPLICIT' was specified, the user may state that the sequence is
    7782 ** 'RANDOM', whereby DIFFaX will generate an explicit sequence of
    7783 ** layers weighted by the stacking probabilities specified later on
    7784 ** under 'PARAMETERS'. The user can specify a crystal thickness (in
    7785 ** terms of the number of layers) under the 'RECURSIVE' option.
    7786 **
    7787 **      ARGUMENTS:
    7788 **            unit_no  -  Logical unit number that the data file is to
    7789 **                        be read from. (input).
    7790 **            messge   -  A short message indicating what went wrong
    7791 **                        (if anything) during the datafile read.
    7792 **                        messge is terminated by a token '$'. (output).
    7793 **
    7794 **      COMMON VARIABLES:
    7795 **            uses:  n_layers
    7796 **
    7797 **        modifies:  recrsv, xplcit, rndm, l_cnt, l_seq, inf_thick
    7798 **
    7799 **      RDSTAK returns logical .true. if the stacking
    7800 **      method was properly specified.
    7801 ** ______________________________________________________________________
    7802 **
    7803 *      logical function RDSTAK(unit_no, messge)
    7804 *      include 'DIFFaX.par'
    7805 *      include 'DIFFaX.inc'
    7806 **
    7807 *      integer*4 unit_no
    7808 *      character*(*) messge
    7809 **
    7810 *      character*200 line, tmpline
    7811 *      character*80 list(2)
    7812 *      integer*4 CHOICE, LAYCNT, LENGTH, i, j, iflag
    7813 **
    7814 ** external functions
    7815 *      external CHOICE, LAYCNT, LENGTH
    7816 ** external subroutines (Some compilers need them declared external)
    7817 **      external GETLNE, TOUPPR, WRTLNE
    7818 **
    7819 *      RDSTAK = .false.
    7820 **
    7821 *      call GETLNE(unit_no, line, *999)
    7822 *      call TOUPPR(line)
    7823 *      if(index(line,'STACKING ').ne.1) then
    7824 *        write(messge,103) 'The ', n_layers,
    7825 *     |   ' layers are read in. ''STACKING'' section not found.$'
    7826 *        goto 999
    7827 *      endif
    7828 **
    7829 *      messge = 'Illegal STACKING description.$'
    7830 *      call GETLNE(unit_no, line, *999)
    7831 *      call TOUPPR(line)
    7832 **
    7833 *      list(1) = 'EXPLICIT '
    7834 *      list(2) = 'RECURSIVE '
    7835 *      iflag = CHOICE(line, list, 2)
    7836 **
    7837 *      if(iflag.lt.1 .or. iflag.gt.2) then
    7838 *        write(op,100) 'Illegal STACKING keyword. Legal types are:'
    7839 *        do 10 i = 1, 2
    7840 *          write(op,101) '               ',list(i)
    7841 *   10   continue
    7842 *        messge = ' $'
    7843 *        goto 999
    7844 *      endif
    7845 **
    7846 ** set flags to indicate stacking type
    7847 *      recrsv    = .false.
    7848 *      xplcit    = .false.
    7849 *      rndm      = .false.
    7850 *      inf_thick = .false.
    7851 ** initialize l_cnt
    7852 *      l_cnt = 0
    7853 *      if(iflag.eq.1) then
    7854 *        xplcit = .true.
    7855 *      else if(iflag.eq.2) then
    7856 *        recrsv = .true.
    7857 *      endif
    7858 **
    7859 *      if(recrsv) then
    7860 *        messge = 'Problem reading crystal thickness for recursion.$'
    7861 *        call GETLNE(unit_no, line, *999)
    7862 *        call TOUPPR(line)
    7863 **
    7864 *        if(index(line,'INFINITE').eq.0) then
    7865 *          inf_thick = .false.
    7866 *          messge =
    7867 *     |       'Problems reading the number of RECURSIVE layers.$'
    7868 ** for strict FORTRAN77 compliance
    7869 *          write(scrtch, *, err=990) line
    7870 *          rewind(scrtch, err=990)
    7871 *          read(scrtch, *, err=999) l_cnt
    7872 *          rewind(scrtch, err=990)
    7873 ** for more lenient, efficient, compilers
    7874 **          read(line, *, err=999) l_cnt
    7875 ** check l_cnt
    7876 *          if(l_cnt.eq.0) then
    7877 *            messge = 'The number of RECURSIVE layers cannot be zero.$'
    7878 *          else if(l_cnt.lt.0) then
    7879 *            messge =
    7880 *     |           'The number of RECURSIVE layers cannot be negative.$'
    7881 *          else if(l_cnt.gt.RCSV_MAX) then
    7882 *            write(op,104)
    7883 *     |       'WARNING: Number of RECURSIVE layers ', l_cnt,
    7884 *     |       ' exceeds the maximum ', RCSV_MAX
    7885 *            write(op,105) 'An INFINITE thickness is assumed.'
    7886 ** re-set l_cnt and proceed assuming infinite thickness
    7887 *            l_cnt = 0
    7888 *            inf_thick = .true.
    7889 *          endif
    7890 *        else
    7891 *          inf_thick = .true.
    7892 *        endif
    7893 **
    7894 *      else if(xplcit) then
    7895 ** read in layer sequence
    7896 *        messge = 'Problem reading EXPLICIT layer sequencing.$'
    7897 *        call GETLNE(unit_no, line, *999)
    7898 *        call TOUPPR(line)
    7899 *        rndm = index(line,'RANDOM ').ne.0
    7900 ** catch possible error - user requested INFINITE
    7901 *        if(.not.rndm) then
    7902 *          if(index(line,'INFINITE').ne.0) then
    7903 *            write(op,107) 'Maximum number of layers allowed is ',XP_MAX
    7904 *            messge = 'Too many layers for an EXPLICIT sequence.$'
    7905 *            goto 999
    7906 *          endif
    7907 *        endif
    7908 **
    7909 *        if(rndm) then
    7910 ** the user asked for a random layer sequencing.
    7911 **
    7912 *          messge =
    7913 *     |       'Problems reading the number of random EXPLICIT layers.$'
    7914 *          write(tmpline, '(a)', err=999)
    7915 *     |        line(LENGTH('RANDOM')+2:len(line))
    7916 ** for strict FORTRAN77 compliance
    7917 *          write(scrtch, *, err=990) tmpline
    7918 *          rewind(scrtch, err=990)
    7919 *          read(scrtch, *, err=999) l_cnt
    7920 *          rewind(scrtch, err=990)
    7921 ** for more lenient, efficient, compilers
    7922 **          read(tmpline, *, err=999) l_cnt
    7923 **
    7924 *        else
    7925 *          l_cnt = 1
    7926 *          messge =
    7927 *     |       'Problems reading EXPLICIT layer sequence in RDSTAK.$'
    7928 *   20     i = LAYCNT(line)
    7929 **
    7930 *          if(i.ge.0 .and. i.le.40) then
    7931 **
    7932 ** for strict FORTRAN77 compliance
    7933 *            write(scrtch, *, err=990) line
    7934 *            rewind(scrtch, err=990)
    7935 *            read(scrtch, *, err=999) (l_seq(j), j=l_cnt,l_cnt+i-1)
    7936 *            rewind(scrtch, err=990)
    7937 ** for more lenient, efficient, compilers
    7938 **            read(line, *, err=999) (l_seq(j), j=l_cnt,l_cnt+i-1)
    7939 **
    7940 *            l_cnt = l_cnt + i
    7941 *          else
    7942 *            goto 999
    7943 *          endif
    7944 **
    7945 *          call GETLNE(unit_no, line, *999)
    7946 *          call TOUPPR(line)
    7947 ** check we have not run into the 'TRANSITIONS' section yet.
    7948 *          if(index(line,'TRANSITIONS').ne.1) goto 20
    7949 *          l_cnt = l_cnt - 1
    7950 **
    7951 ** If we are here, we have gone one line too far. Back up.
    7952 *          write(messge,103)
    7953 *     |       'Problems ''backspacing'' unit ', unit_no, '.$'
    7954 *          backspace(unit = unit_no, err = 980)
    7955 **
    7956 *        endif
    7957 **
    7958 ** check we didn't enter too many layers
    7959 *        if(l_cnt.gt.XP_MAX) then
    7960 *          write(op,104) 'WARNING: No. of EXPLICIT layers, ',l_cnt,
    7961 *     |      ' exceeds maximum of ',XP_MAX
    7962 *          write(op,106) XP_MAX, ' EXPLICIT layers assumed'
    7963 *          l_cnt = XP_MAX
    7964 *        endif
    7965 **
    7966 ** check that all the layer numbers are within bounds
    7967 *        if(.not.rndm) then
    7968 *          do 30 i = 1, l_cnt
    7969 *            if(l_seq(i).gt.n_layers) then
    7970 *              write(messge,102) 'Illegal layer number,', l_seq(i),
    7971 *     |         '. Number must not exceed', n_layers, '.$'
    7972 *              goto 980
    7973 *            else if(l_seq(i).lt.1) then
    7974 *              write(messge,103) 'Illegal layer number,', l_seq(i),
    7975 *     |         '. Number must be greater than 0.$'
    7976 *              goto 980
    7977 *            endif
    7978 *   30     continue
    7979 *        endif
    7980 *      endif
    7981 **
    7982 *      RDSTAK = .true.
    7983 *  980 return
    7984 *  990 messge = 'Problems using scratch file in RDSTAK.$'
    7985 *  999 call WRTLNE(line)
    7986 *      return
    7987 *  100 format(1x, 'ERROR: ', a)
    7988 *  101 format(1x, 2a)
    7989 *  102 format(2(a, i4), a)
    7990 *  103 format(a, i2, a)
    7991 *  104 format(1x, 2(a, i4))
    7992 *  105 format(1x, a)
    7993 *  106 format(1x, i4, a)
    7994 *  107 format(1x, a, i5)
    7995 *      end
    7996 **
    7997 ** ______________________________________________________________________
    7998 ** Title: RDTRNS
    7999 ** Authors: MMJT
    8000 ** Date: 4 Oct 1989; 15 Mar 1995; 29th Aug 2000
    8001 ** Description: This function reads the layer stacking probabilities
    8002 ** and stacking vectors.
    8003 **
    8004 **      ARGUMENTS:
    8005 **            unit_no  -  Logical unit number that the data file is to
    8006 **                        be read from. (input).
    8007 **            messge   -  A short message indicating what went wrong
    8008 **                        (if anything) during the datafile read.
    8009 **                        messge is terminated by a token '$'. (output).
    8010 **
    8011 **      COMMON VARIABLES:
    8012 **            uses:  n_layers, recrsv, rndm, xplcit
    8013 **
    8014 **        modifies:  l_alpha, l_r, r_B11, r_B22, r_B33, r_B12, r_B23,
    8015 **                   r_B31
    8016 **
    8017 **      RDTRNS returns logical .true. if all the stacking
    8018 **      parameters were properly specified.
    8019 ** ______________________________________________________________________
    8020 **
    8021 *      logical function RDTRNS(unit_no, messge)
    8022 *      include 'DIFFaX.par'
    8023 *      include 'DIFFaX.inc'
    8024 **
    8025 *      integer*4 unit_no
    8026 *      character*(*) messge
    8027 **
    8028 *      character*200 line, tmpline, arg
    8029 *      integer*4 PARENT, NXTARG, i, j, n, tempfact
    8030 *      real*8 sum, RDNMBR
    8031 *      logical*4 ok
    8032 **
    8033 ** external function
    8034 *      external PARENT, NXTARG, RDNMBR
    8035 ** external subroutines (Some compilers need them declared external)
    8036 **      external GETLNE, TOUPPR, WRTLNE
    8037 **
    8038 *      RDTRNS = .false.
    8039 **
    8040 *      call GETLNE(unit_no, line, *999)
    8041 *      call TOUPPR(line)
    8042 *      if(index(line,'TRANSITIONS ').ne.1) then
    8043 *        messge = '''TRANSITIONS'' section not found.$'
    8044 *        goto 999
    8045 *      endif
    8046 **
    8047 *      do 10 i = 1, n_layers
    8048 *        do 20 j = 1, n_layers
    8049 *          write(messge,100)
    8050 *     |     'Unexpected EOF at vector ',i,j,' in ''TRANSITIONS''.$'
    8051 *          call GETLNE(unit_no, line, *999)
    8052 ** Track down if there are parentheses in this line. If so then there
    8053 ** are FatsWaller (Debye-Waller for layers!) factors
    8054 *          write(messge,100)
    8055 *     |      'Confusing use of parentheses in data for vector ',i,j,'.$'
    8056 *          tempfact = PARENT(line)
    8057 *          if(tempfact.lt.0) goto 999
    8058 ** initialize
    8059 *          l_r(1,j,i) = ZERO
    8060 *          l_r(2,j,i) = ZERO
    8061 *          l_r(3,j,i) = ZERO
    8062 *
    8063 *          r_B11(j,i) = ZERO
    8064 *          r_B22(j,i) = ZERO
    8065 *          r_B33(j,i) = ZERO
    8066 *          r_B12(j,i) = ZERO
    8067 *          r_B31(j,i) = ZERO
    8068 *          r_B23(j,i) = ZERO
    8069 **
    8070 *          tmpline = line
    8071 **
    8072 *          if(tempfact.eq.1) then
    8073 *            write(messge,100)
    8074 *     |     'Expected alpha_ij, R_ij, (dR_ij) for vector ',i,j,'.$'
    8075 *          else
    8076 *            write(messge,100)
    8077 *     |     'Expected alpha_ij, R_ij data for vector ',i,j,'.$'
    8078 *          endif
    8079 **
    8080 ** read stacking probability
    8081 *          n = NXTARG(tmpline,arg)
    8082 *          if(n.le.0) goto 999
    8083 *          l_alpha(j,i) = RDNMBR(arg,ok)
    8084 *          if(.not.ok) goto 999
    8085 ** If zero, there is no need to read further
    8086 *          if(dabs(l_alpha(j,i)).lt.eps6) goto 20
    8087 **           
    8088 ** read stacking x-vector
    8089 *          n = NXTARG(tmpline,arg)
    8090 *          if(n.le.0) goto 999
    8091 *          l_r(1,j,i) = RDNMBR(arg,ok)
    8092 *          if(.not.ok) goto 999
    8093 **           
    8094 ** read stacking y-vector
    8095 *          n = NXTARG(tmpline,arg)
    8096 *          if(n.le.0) goto 999
    8097 *          l_r(2,j,i) = RDNMBR(arg,ok)
    8098 *          if(.not.ok) goto 999
    8099 **           
    8100 ** read stacking z-vector
    8101 *          n = NXTARG(tmpline,arg)
    8102 *          if(n.le.0) goto 999
    8103 *          l_r(3,j,i) = RDNMBR(arg,ok)
    8104 *          if(.not.ok) goto 999
    8105 **
    8106 ** Read temperature factors, if present
    8107 *          if(tempfact.eq.1) then
    8108 **           
    8109 ** read B11 "Fats-Waller"
    8110 *            n = NXTARG(tmpline,arg)
    8111 *            if(n.le.0) goto 999
    8112 *            r_B11(j,i) = RDNMBR(arg,ok)
    8113 *            if(.not.ok) goto 999
    8114 **           
    8115 ** read B22 "Fats-Waller"
    8116 *            n = NXTARG(tmpline,arg)
    8117 *            if(n.le.0) goto 999
    8118 *            r_B22(j,i) = RDNMBR(arg,ok)
    8119 *            if(.not.ok) goto 999
    8120 **           
    8121 ** read B33 "Fats-Waller"
    8122 *            n = NXTARG(tmpline,arg)
    8123 *            if(n.le.0) goto 999
    8124 *            r_B33(j,i) = RDNMBR(arg,ok)
    8125 *            if(.not.ok) goto 999
    8126 **           
    8127 ** read B12 "Fats-Waller"
    8128 *            n = NXTARG(tmpline,arg)
    8129 *            if(n.le.0) goto 999
    8130 *            r_B12(j,i) = RDNMBR(arg,ok)
    8131 *            if(.not.ok) goto 999
    8132 **           
    8133 ** read B31 "Fats-Waller"
    8134 *            n = NXTARG(tmpline,arg)
    8135 *            if(n.le.0) goto 999
    8136 *            r_B31(j,i) = RDNMBR(arg,ok)
    8137 *            if(.not.ok) goto 999
    8138 **           
    8139 ** read B23 "Fats-Waller"
    8140 *            n = NXTARG(tmpline,arg)
    8141 *            if(n.le.0) goto 999
    8142 *            r_B23(j,i) = RDNMBR(arg,ok)
    8143 *            if(.not.ok) goto 999
    8144 **           
    8145 *          endif
    8146 **
    8147 ** stacking probabilities should not be negative.
    8148 *          if(l_alpha(j,i).lt.ZERO) then
    8149 *            write(messge,100) 'alpha',i,j,' is negative.$'
    8150 *            goto 999
    8151 *          endif
    8152 ** the Rz(i,i) vectors should not be zero or negative.
    8153 ** if an i-i transition exists
    8154 *          if(j.eq.i) then
    8155 *            if(l_alpha(i,i).ne.ZERO) then
    8156 *              if(l_r(3,i,i).le.ZERO) then
    8157 *                write(messge,100)
    8158 *     |'Invalid negative (or zero) z-component for Rz',i,i,'.$'
    8159 *                goto 999
    8160 *              endif
    8161 *            endif
    8162 *          endif
    8163 *   20   continue
    8164 *   10 continue
    8165 **
    8166 ** Check stacking probabilities.
    8167 ** Generate error message if they do not sum to 1, and we are using
    8168 ** a RECURSIVE stacking description. If EXPLICIT stacking, then we
    8169 ** do not use the probabilities (but we still assume they are there
    8170 ** to be read) but generate a warning.
    8171 **
    8172 *      do 30 i = 1, n_layers
    8173 *        sum = ZERO
    8174 *        do 40 j = 1, n_layers
    8175 *          sum = sum + l_alpha(j,i)
    8176 *   40   continue
    8177 *        if(abs(sum - ONE).gt.eps6) then
    8178 *          write(messge,101)
    8179 *     |     'Stacking probabilities from LAYER ',i,' do not sum to 1.$'
    8180 *          if(recrsv.or.rndm) then
    8181 *            goto 980
    8182 *          else if(xplcit .and. .not.rndm) then
    8183 *            write(op,102) messge(1:index(messge, '$') - 1)
    8184 *          endif
    8185 *        endif
    8186 *   30 continue
    8187 **
    8188 *      RDTRNS = .true.
    8189 *  980 return
    8190 *  999 call WRTLNE(line)
    8191 *      return
    8192 *  100 format(a, '(', i2, ',', i2, ')', a)
    8193 *  101 format(a, i2, a)
    8194 *  102 format(1x, 'WARNING: ', a)
    8195 *      end
    8196 **
    8197 ** ______________________________________________________________________
    8198 ** Title: RDWAVL
    8199 ** Authors: MMJT
    8200 ** Date: 4 Oct 1989
    8201 ** Description: This function reads the radiation wavelength.
    8202 **
    8203 **      ARGUMENTS:
    8204 **            unit_no  -  Logical unit number that the data file is to
    8205 **                        be read from. (input).
    8206 **            messge   -  A short message indicating what went wrong
    8207 **                        (if anything) during the datafile read.
    8208 **                        messge is terminated by a token '$'. (output).
    8209 **
    8210 **      COMMON VARIABLES:
    8211 **        modifies:  lambda
    8212 **
    8213 **      RDWAVL returns logical .true. if an acceptable value
    8214 **      for the radiation wavelength was read.
    8215 ** ______________________________________________________________________
    8216 **
    8217 *      logical function RDWAVL(unit_no, messge)
    8218 *      include 'DIFFaX.par'
    8219 *      include 'DIFFaX.inc'
    8220 **
    8221 *      integer*4 unit_no
    8222 *      character*(*) messge
    8223 **
    8224 *      character*200 line
    8225 **
    8226 ** external subroutines (Some compilers need them declared external)
    8227 **      external GETLNE, WRTLNE
    8228 **
    8229 *      RDWAVL = .false.
    8230 **
    8231 *      messge = 'lambda is specified improperly.$'
    8232 *      call GETLNE(unit_no, line, *999)
    8233 **
    8234 ** for strict FORTRAN77 compliance
    8235 *      write(scrtch, *, err=990) line
    8236 *      rewind(scrtch, err=990)
    8237 *      read(scrtch, *, err=999) lambda
    8238 *      rewind(scrtch, err=990)
    8239 **
    8240 ** for more lenient, efficient compilers
    8241 **      read(line, *, err=999) lambda
    8242 **
    8243 *      if(lambda.le.ZERO) then
    8244 *        messge = 'lambda is negative (or zero).$'
    8245 *        goto 999
    8246 *      endif
    8247 **
    8248 *      RDWAVL = .true.
    8249 *      return
    8250 *  990 messge = 'Problems using scratch file in RDWAVL.$'
    8251 *  999 call WRTLNE(line)
    8252 *      return
    8253 *      end
    8254 **
    8255 ** ______________________________________________________________________
    8256 ** Title: RDWDTH
    8257 ** Authors: MMJT
    8258 ** Date: 3 Mar 1995
    8259 ** Description: This function reads the lateral layer dimensions
    8260 ** from the data file. First we check that the data on the line does
    8261 ** not refer to the layer atomic positions. If it does, assume infinite
    8262 ** widths, backspace the file, and return to let RDLAYR do its job.
    8263 ** Unlike the other read functions, RDWDTH does not return an error in
    8264 ** this instance, so as to maintain compatibility with data files
    8265 ** written for v 1.7xx which did not treat lateral layer width
    8266 ** broadening.
    8267 **
    8268 **      ARGUMENTS:
    8269 **            unit_no  -  Logical unit number that the data file is to
    8270 **                        be read from. (input).
    8271 **            messge   -  A short message indicating what went wrong
    8272 **                        (if anything) during the datafile read.
    8273 **                        messge is terminated by a token '$'. (output).
    8274 **
    8275 **      COMMON VARIABLES:
    8276 **        modifies:  finite_width, Wa, Wb
    8277 **
    8278 **      RDWAVL returns logical .true. if no errors were encountered.
    8279 ** ______________________________________________________________________
    8280 **
    8281 *      logical function RDWDTH(unit_no, messge)
    8282 *      include 'DIFFaX.par'
    8283 *      include 'DIFFaX.inc'
    8284 **
    8285 *      integer*4 unit_no
    8286 *      character*(*) messge
    8287 **
    8288 *      character*200 line
    8289 **
    8290 *      integer*4 n, CNTARG
    8291 **
    8292 ** external function
    8293 *      external CNTARG
    8294 ** external subroutines (Some compilers need them declared external)
    8295 **      external GETLNE, TOUPPR, WRTLNE
    8296 **
    8297 *      RDWDTH = .false.
    8298 **
    8299 ** Does the next line belong to the layer structural data? If so,
    8300 ** assume infinite layer widths. No error message is generated so
    8301 ** that there is compatibility with data files for versions 1.7xx.
    8302 *      messge = 'layer width, or layer data, are specified improperly.$'
    8303 *      call GETLNE(unit_no, line, *999)
    8304 *      call TOUPPR(line)
    8305 *      if(index(line, 'LAYER').eq.1) then
    8306 *        finite_width = .false.
    8307 *        write(messge,100)
    8308 *     |   'Problems ''backspacing'' unit ',unit_no, '.$'
    8309 ** backspace file so that RDLAYR will work properly
    8310 *        backspace(unit = unit_no, err = 990)
    8311 *        goto 900
    8312 *      endif
    8313 **
    8314 ** assume that layer width data is specified
    8315 *      if(index(line,'INFINITE').eq.1) then
    8316 *        finite_width = .false.
    8317 *        goto 900
    8318 *      else
    8319 *        finite_width = .true.
    8320 *      endif
    8321 **
    8322 ** are both Wa and Wb specified, or just a generic width?
    8323 *      messge = 'layer characteristic widths are specified improperly.$'
    8324 *      n = CNTARG(line)
    8325 **
    8326 ** for strict FORTRAN77 compliance
    8327 *      write(scrtch, *, err=990) line
    8328 *      rewind(scrtch, err=990)
    8329 **
    8330 *      if(n.eq.1) then
    8331 ** for strict FORTRAN77 compliance
    8332 *        read(scrtch, *, err=999) Wa
    8333 *        rewind(scrtch, err=990)
    8334 ** for more lenient, efficient compilers
    8335 **        read(line, *, err=999) Wa
    8336 *        Wb = Wa
    8337 *      else if(n.eq.2) then
    8338 ** for strict FORTRAN77 compliance
    8339 *        read(scrtch, *, err=999) Wa, Wb
    8340 *        rewind(scrtch, err=990)
    8341 ** for more lenient, efficient compilers
    8342 **        read(line, *, err=999) Wa, Wb
    8343 *      else
    8344 *        goto 999
    8345 *      endif
    8346 **
    8347 *      if(finite_width) then
    8348 *        if(Wa.gt.inf_width .and. Wb.gt.inf_width) then
    8349 *          finite_width = .false.
    8350 *          write(op,101) 'Layers will be treated as if infinitely wide.'
    8351 *        endif
    8352 *        if(Wa.le.ZERO .or. Wb.le.ZERO) then
    8353 *          write(op,101)'Illegal values for layer characteristic widths'
    8354 *          goto 999
    8355 *        endif
    8356 *      endif
    8357 **
    8358 *  900 RDWDTH = .true.
    8359 *      return
    8360 *  990 messge = 'Problems using scratch file in RDWDTH.$'
    8361 *  999 call WRTLNE(line)
    8362 *      return
    8363 *  100 format(a, i2, a)
    8364 *  101 format(1x, a)
    8365 *      end
    8366 **
    8367 ** ______________________________________________________________________
    83685029* Title: RAN3
    83695030* Authors: Press, Flannery, Teukolsky and Vetterling
    … …  
    84335094*
    84345095* ______________________________________________________________________
    8435 * Title: SALUTE
    8436 * Author: MMJT
    8437 * Date: 19th May, 2010
    8438 * Draws the DIFFaX flag. Saluting is optional.
    8439 *
    8440 *      ARGUMENTS:
    8441 *           No arguments are needed.
    8442 * ______________________________________________________________________
    8443 *
    8444       subroutine SALUTE()
    8445       include 'DIFFaX.par'
    8446 *     save
    8447 *
    8448       write(op,1) '***************************************************'
    8449       write(op,1) '***************************************************'
    8450       write(op,1) '*                                                 *'
    8451       write(op,1) '*  DDDD     II   FFFFFF   FFFFFF           X   X  *'
    8452       write(op,1) '*  D   D    II   F        F        aaaa     X X   *'
    8453       write(op,1) '*  D    D   II   FFFF     FFFF    a    a     X    *'
    8454       write(op,1) '*  D   D    II   F        F       a   aa    X X   *'
    8455       write(op,1) '*  DDDD     II   F        F        aaa a   X   X  *'
    8456       write(op,1) '*                                                 *'
    8457       write(op,1) '***************************************************'
    8458       write(op,1) '****************** DIFFaX v1.813 ******************'
    8459       write(op,1) '***************************************************'
    8460       write(op,1) '***************** 19th May,  2010 *****************'
    8461       write(op,1) '***************************************************'
    8462       write(op,1) '*                                                 *'
    8463       write(op,1) '*   A computer program for calculating            *'
    8464       write(op,1) '*   Diffraction Intensity From Faulted Crystals   *'
    8465       write(op,1) '*                                                 *'
    8466       write(op,1) '*   Authors: Michael M. J. Treacy                 *'
    8467       write(op,1) '*            Michael W. Deem                      *'
    8468       write(op,1) '*                                                 *'
    8469       write(op,1) '***************************************************'
    8470       write(op,1)
    8471 *
    8472       return
    8473     1 format(1x, a51)
    8474       end
    8475 *
    8476 * ______________________________________________________________________
    8477 ** Title: SFC
    8478 ** Authors: MWD and MMJT
    8479 ** Date: 11 Feb 90; 7 Mar 1995
    8480 ** Description: This routine reads the atomic scattering factor constants
    8481 ** from the file 'sfname' (usually given as "data.sfc"). The scattering
    8482 ** factors are used to generate the layer form factors. SFC returns
    8483 ** .false. if there are too many distinct atom types (ie. more than
    8484 ** MAX_TA). The value of MAX_TA can be reset in 'DIFFaX.par'.
    8485 **
    8486 **      ARGUMENTS:
    8487 **           No arguments are used. All data is in 'COMMON'.
    8488 **
    8489 **      COMMON VARIABLES:
    8490 **            uses:  n_actual, l_n_atoms, a_name, n_atoms, rad_type
    8491 **                   NEUTRN, X_RAY, rad_type, sfname
    8492 **
    8493 **        modifies:  n_atoms, x_sf, n_sf, e_sf, a_type
    8494 **
    8495 **      SFC returns logical .true. if it read the scattering factor data
    8496 **      safely.
    8497 ** ______________________________________________________________________
    8498 **
    8499 *      logical function SFC()
    8500 *      include 'DIFFaX.par'
    8501 *      include 'DIFFaX.inc'
    8502 **
    8503 *      logical ok, new_atom, our_atom, done
    8504 *      logical list(MAX_TA+1)
    8505 *      integer*4 i, j, n, m, LENGTH
    8506 *      character*4 name
    8507 *      character*120 line
    8508 *      real*4 tmp
    8509 **
    8510 ** external subroutine (Some compilers need them declared external)
    8511 **      external ATOMS
    8512 ** external function
    8513 *      external LENGTH
    8514 **
    8515 *      write(op,301) 'Reading scattering factor datafile ''',
    8516 *     |                   sfname(1:LENGTH(sfname)),'''. . .'
    8517 **
    8518 *      SFC = .false.
    8519 *      do 10 i = 1, MAX_TA
    8520 *        list(i) = .false.
    8521 *        atom_l(i) = ' '
    8522 *   10 continue
    8523 *      m = 0
    8524 ** determine all distinct atom types
    8525 *      do 30 i = 1, n_actual
    8526 *        do 40 j = 1, l_n_atoms(i)
    8527 *          name = a_name(j,i)
    8528 *          call TOUPPR(name)
    8529 ** see if this name has been seen already
    8530 *          n = 0
    8531 *   50     n = n + 1
    8532 *            new_atom = name.ne.atom_l(n)
    8533 *            if(new_atom .and. n.lt.m) goto 50
    8534 *          if(new_atom) then
    8535 *            m = m + 1
    8536 *            if(m.gt.MAX_TA) goto 220
    8537 *            atom_l(m) = name
    8538 *            a_type(j,i) = m
    8539 *          else
    8540 *            a_type(j,i) = n
    8541 *          endif
    8542 *   40   continue
    8543 *   30 continue
    8544 *      n_atoms = m
    8545 ** now find data for each atom type in file
    8546 ** pass through file only once
    8547 *   60 read(sf, 300, end=90, err=200) line
    8548 *        name = line(1:4)
    8549 *        call TOUPPR(name)
    8550 ** see if this is one of the distinct atoms
    8551 *        i = 0
    8552 *   70   i = i + 1
    8553 *          our_atom = name.eq.atom_l(i)
    8554 *        if(.not.our_atom .and. i.lt.n_atoms) goto 70
    8555 ** have we read all that we need to know?
    8556 *        done = .true.
    8557 *        do 80 n = 1, n_atoms
    8558 *          done = done.and.list(n)
    8559 *   80   continue
    8560 ** If we're done, close file.
    8561 *        if(done) goto 90
    8562 *        if(our_atom .and. .not.list(i)) then
    8563 ** mark this atom's data as read in
    8564 *          list(i) = .true.
    8565 ** and read it in
    8566 *          if(rad_type.eq.X_RAY) then
    8567 *            read(line,310,err=210) (x_sf(j,i), j=1, 9)
    8568 *          else if(rad_type.eq.NEUTRN) then
    8569 *            read(line,320,err=210) n_sf(i)
    8570 *          else
    8571 *            read(line,310,err=210) (x_sf(j,i), j=1, 9), tmp, e_sf(i)
    8572 *          endif
    8573 *        endif
    8574 *      goto 60
    8575 *   90 close(sf,err=240)
    8576 ** see if all the data for each atom has been read in
    8577 *      ok = .true.
    8578 *      do 100 i = 1, n_atoms
    8579 *        if(.not.list(i)) then
    8580 *          ok = .false.
    8581 *          write(op,330) 'ERROR: Data for atom ''', atom_l(i),
    8582 *     |     ''' not found in file ''', sfname(1:LENGTH(sfname)), ''''
    8583 *          call ATOMS()
    8584 *        endif
    8585 *  100 continue
    8586 *      if(ok) then
    8587 *        SFC = .true.
    8588 *        write(op,400) 'Scattering factor data read in.'
    8589 *      endif
    8590 *      return
    8591 *  200 write(op,301) 'Scattering factor file ''',
    8592 *     |    sfname(1:LENGTH(sfname)), ''' defective.'
    8593 *      close(sf,err=240)
    8594 *      return
    8595 *  210 write(op,400) 'ERROR reading scattering factor data.'
    8596 *      close(sf,err=240)
    8597 *      return
    8598 *  220 write(op,402) 'There are too many types of atoms in layer ', i
    8599 *      write(op,400) 'Atoms recorded so far are:'
    8600 *      do 110 j = 1, MAX_TA
    8601 *        write(op,403) '       Atom type ', j, '   ', atom_l(j)
    8602 *  110 continue
    8603 *      write(op,402) '  Maximum number of types allowed is ', MAX_TA
    8604 *      return
    8605 *  240 write(op,301) 'Unable to close scattering factor file ''',
    8606 *     |    sfname(1:LENGTH(sfname)), '''.'
    8607 *      return
    8608 *  300 format(a)
    8609 *  301 format(1x, 3a)
    8610 *  310 format(t5, 10f11.6, i3)
    8611 *  320 format(t104, f11.6)
    8612 *  330 format(1x, 5a)
    8613 *  400 format(1x, a)
    8614 *  401 format(1x, 2a)
    8615 *  402 format(1x, a, i2)
    8616 *  403 format(1x, a, i2, 2a)
    8617 *      end
    8618 **
    8619 ** ______________________________________________________________________
    86205096* Title: SHARP
    86215097* Author: MMJT
    … …  
    88035279*
    88045280* ______________________________________________________________________
    8805 ** Title: STREAK
    8806 ** Author: MWD & MMJT
    8807 ** Date: 18 Aug 1988, revised 22 June 1989.
    8808 ** Description:  This routine outputs the integrated (ie. averaged)
    8809 ** intensitites from h,k,l0 to h,k,l1, in steps of dl, to the file
    8810 ** strkfile.
    8811 **
    8812 **      ARGUMENTS:
    8813 **            FN        -  Function name passed by reference. The
    8814 **                         choice is between GLQ16 (non-adaptive
    8815 **                         Gauss-Legendre integration), and AGLQ16
    8816 **                         (adaptive Gauss-Legendre integration).
    8817 **                                                            (input).
    8818 **            strkfile  -  The name of the file that streak data is
    8819 **                         to be output to. (input).
    8820 **            ok        -  logical flag indicating all went well.
    8821 **                                                      (output).
    8822 **
    8823 **      COMMON VARIABLES:
    8824 **            uses:  a0, b0, c0, d0, lambda, cntrl, CFile, xplcit,
    8825 **                   rad_type, X_RAY
    8826 **
    8827 **        modifies:  no COMMON variables are modified
    8828 ** ______________________________________________________________________
    8829 **
    8830 *      subroutine STREAK(FN, strkfile, ok)
    8831 *      include 'DIFFaX.par'
    8832 *      include 'DIFFaX.inc'
    8833 **
    8834 *      real*8 FN
    8835 *      character*(*) strkfile
    8836 *      logical ok
    8837 **
    8838 *      logical its_hot
    8839 *      integer*4 h, k, LENGTH, i, i_step
    8840 *      real*8 l, theta, x, ANGLE, S, Q2, l0, l1
    8841 *      real*8 dl, W4, intervals
    8842 *      parameter (intervals = TWENTY)
    8843 **
    8844 ** external functions
    8845 *      external FN, LENGTH
    8846 ** external subroutines (Some compilers need them declared external)
    8847 **      external XYPHSE, PRE_MAT, GET_F
    8848 **
    8849 ** statement functions
    8850 ** S is the value of 1/d**2 at hkl
    8851 *      S(h,k,l) = h*h*a0 + k*k*b0 + l*l*c0 + h*k*d0
    8852 ** ANGLE is the Bragg angle (in radians) of the h,k,l plane
    8853 *      ANGLE(h,k,l) = asin(HALF * lambda * sqrt(S(h,k,l)))
    8854 ** W4 is the X-ray polarization factor
    8855 *      W4(theta) = HALF * (ONE + (cos(TWO*theta))**2)
    8856 **
    8857 *      Q2 = FOUR / (lambda**2)
    8858 *   10 write(op,400) 'Enter h, k, l0, l1, delta l'
    8859 *      read(cntrl,*,err=10) h, k, l0, l1, dl
    8860 *      if(CFile) write(op,401) h, k, l0, l1, dl
    8861 ** check input
    8862 *      if(l1.eq.l0) then
    8863 *        write(op,400) 'Illegal input: l0 equals l1'
    8864 *        goto 999
    8865 *      else if(dl.eq.ZERO) then
    8866 *        write(op,400) 'Illegal zero value of dl entered'
    8867 *        write(op,402)'A value of ',(l1-l0)/(FIVE*HUNDRED),' is assumed'
    8868 *      else if(l1.gt.l0 .and. dl.lt.ZERO) then
    8869 *        write(op,400) 'l1 is greater than l0. +ve dl assumed'
    8870 *        dl = -dl
    8871 *      else if(l1.lt.l0 .and. dl.gt.ZERO) then
    8872 *        write(op,400) 'l0 is greater than l1. -ve dl assumed'
    8873 *        dl = -dl
    8874 *      endif
    8875 ** The origin may be hotter than hell! Let's check first.
    8876 *      its_hot = (h.eq.0 .and. k.eq.0) .and.
    8877 *     |         l0*l1.le.ZERO .and. rad_type.eq.ELECTN
    8878 *      if(its_hot) then
    8879 *        write(op,400) 'Cannot scan the origin with electron radiation'
    8880 *        write(op,400) 'Origin will be skipped.'
    8881 *      endif
    8882 ** check angles are meaningful
    8883 *      if(S(h,k,l0).gt.Q2 .or. S(h,k,l1).gt.Q2) then
    8884 *        if(S(h,k,l0).gt.Q2) write(op,403) h, k, l0,
    8885 *     |            ' exceeds 180 degree scattering angle!'
    8886 *        if(S(h,k,l1).gt.Q2) write(op,403) h, k, l1,
    8887 *     |            ' exceeds 180 degree scattering angle!'
    8888 *        goto 10
    8889 *      endif
    8890 **
    8891 *      write(op,404) 'Writing streak data to file ''',
    8892 *     |      strkfile(1:LENGTH(strkfile)),'''. . .'
    8893 *      call XYPHSE(h, k)
    8894 *      call PRE_MAT(h, k)
    8895 **
    8896 *      i_step = nint( (l1 - l0) / (dl * intervals) )
    8897 ** Immune system to the rescue!
    8898 *      if(i_step.le.0) i_step = 10
    8899 *      i = 0
    8900 *      do 20 l = l0, l1, dl
    8901 ** If we are dealing with electrons, make sure we avoid the origin
    8902 *        if(its_hot .and. l*(l+dl).le.ZERO) then
    8903 *          x = ZERO
    8904 *          goto 30
    8905 *        endif
    8906 *        i = i + 1
    8907 *        if(mod(i,i_step).eq.0) write(op,405) 'Reached l = ',l
    8908 *        x = FN(h,k,l,l+dl,ok)
    8909 *        if(.not.ok) goto 999
    8910 ** note: since this is streak data, only the X_RAY input needs
    8911 ** correcting for polarization.
    8912 *        if(rad_type.eq.X_RAY)  x = x * W4(ANGLE(h,k,l+HALF*dl))
    8913 *   30   write(sk,406,err=100) l, char(9), x
    8914 *   20 continue
    8915 *      if(sk.ne.op) close(sk,err=110)
    8916 *      write(op,404) 'Streak data file, ''',
    8917 *     |      strkfile(1:LENGTH(strkfile)),''' written to disk.'
    8918 *      return
    8919 *  100 write(op,404) 'ERROR writing to streak data file ''',
    8920 *     |      strkfile(1:LENGTH(strkfile)),''''
    8921 *      if(sk.ne.op) close(sk,err=110)
    8922 *      return
    8923 *  110 write(op,404) 'Unable to close streak data file ''',
    8924 *     |      strkfile(1:LENGTH(strkfile)),''''
    8925 *      return
    8926 *  999 write(op,405) 'ERROR encountered in streak integration at l = ',l
    8927 *      return
    8928 *  400 format(1x, a)
    8929 *  401 format(1x, 2i3, 3f10.5)
    8930 *  402 format(1x, a, f10.5, a)
    8931 *  403 format(1x, 2i3, f10.5, a)
    8932 *  404 format(1x, 3a)
    8933 *  405 format(1x, a, f10.5)
    8934 *  406 format(1x, e12.5, a, e14.6)
    8935 *      end
    8936 **
    8937 ** ______________________________________________________________________
    89385281* Title: THRESH
    89395282* Author: MMJT
    … …  
    90195362*
    90205363* ______________________________________________________________________
    9021 ** Title: TOUPPR
    9022 ** Author: MMJT
    9023 ** Date: 3 August 1989
    9024 ** Description: Converts alphabetic characters in 'line' to uppercase.
    9025 ** Works for ascii characters, since it assumes alphabetic characters
    9026 ** are contiguous.
    9027 ** The neat trick of using
    9028 **                 lower_case = ichar(chr).and.95
    9029 ** will not always work in ansi77 fortran, as some compilers evaluate
    9030 ** the right hand side as a logical .true. or .false. rather than a
    9031 ** bitwise 'anded' integer.
    9032 ** The alternative
    9033 **                 lower_case = ichar(chr).iand.95
    9034 ** is equally undesirable since the operator .iand. is a military
    9035 ** extension, and is not ansi77.
    9036 **
    9037 **      ARGUMENTS:
    9038 **            line  -  Input line of characters to be converted to
    9039 **                     upper case. The converted string is returned
    9040 **                     in 'line'.
    9041 ** ______________________________________________________________________
    9042 **
    9043 *      subroutine TOUPPR(line)
    9044 **
    9045 *      character*(*) line
    9046 **
    9047 *      integer*4 a, z, offset, i, itmp, lin_len
    9048 **
    9049 *      a = ichar('a')
    9050 *      z = ichar('z')
    9051 *      offset = a - ichar('A')
    9052 *      lin_len = len(line)
    9053 **
    9054 *      i = 1
    9055 *   10 itmp = ichar(line(i:i))
    9056 *        if(itmp.ge.a .and. itmp.le.z) line(i:i) = char(itmp - offset)
    9057 *        i = i + 1
    9058 *        if(i.le.lin_len)
    9059 *     |goto 10
    9060 **
    9061 *      return
    9062 *      end
    9063 **
    9064 ** ______________________________________________________________________
    90655364* Title: TRMSPC
    90665365* Author: MMJT
    … …  
    96325931*
    96335932* ______________________________________________________________________
    9634 ** Title: WRTLNE
    9635 ** Author: MMJT
    9636 ** Date: 4 Oct 1989
    9637 ** Description: This subroutine writes out the contents of
    9638 ** the character string 'line' with an appropriate preamble.
    9639 ** This routine strips off trailing blanks.
    9640 **
    9641 **      ARGUMENTS:
    9642 **            line     -  Line of characters to write out. (input).
    9643 ** ______________________________________________________________________
    9644 **
    9645 *      subroutine WRTLNE(line)
    9646 *      include 'DIFFaX.par'
    9647 **     save
    9648 **
    9649 *      character*(*) line
    9650 **
    9651 *      integer*4 lin_len, i
    9652 **
    9653 *      lin_len = len(line)
    9654 *      i = lin_len
    9655 *   10 if(i.gt.0) then
    9656 *        if(line(i:i).eq.' ') then
    9657 *          if(i.gt.1) then
    9658 *            i = i - 1
    9659 *            goto 10
    9660 *          endif
    9661 *        endif
    9662 *      endif
    9663 **
    9664 *      write(op,100) 'The data on the current line was read as:'
    9665 *      write(op,101) '  ''', line(1:i), ''''
    9666 **
    9667 *      return
    9668 *  100 format(1x, a)
    9669 *  101 format(1x, 3a)
    9670 *      end
    9671 **
    9672 ** ______________________________________________________________________
    9673 ** Title: WRTSAD
    9674 ** Author: MMJT
    9675 ** Date: 29 Oct 1989; 1st Mar 2000; 19th May 2010
    9676 ** Description: This subroutine writes the selected area diffraction
    9677 ** pattern (SADP) data to a binary file.
    9678 ** Binary file output has been problematical because some compilers add
    9679 ** record-length specifiers at the beginning and ending of each line.
    9680 ** Sometimes the specifiers are 2-byte, or 4-byte, and even 0-byte!
    9681 ** A statement to screen specifies how to compute the true row width.
    9682 **
    9683 **      ARGUMENTS:
    9684 **            outfile  -  The name of the file that the binary
    9685 **                        selected area diffraction data is to be
    9686 **                        written to. (input).
    9687 **            view     -  Choice of beam direction. (input).
    9688 **                              1  =  normal to the plane k = 0
    9689 **                              2  =  normal to the plane h = 0
    9690 **                              3  =  normal to the plane h = k
    9691 **                              4  =  normal to the plane h = -k
    9692 **            l_upper  -  Upper limit of l. (input).
    9693 **            hk_lim   -  Upper limit of h (or k). (input).
    9694 **            ok       -  logical flag indicating all went well.
    9695 **                                                      (output).
    9696 **
    9697 **      COMMON VARIABLES:
    9698 **            uses:  a0, b0, c0, d0, sadblock, scaleint, has_l_mirror
    9699 **
    9700 **        modifies:  spec
    9701 ** ______________________________________________________________________
    9702 **
    9703 *      subroutine WRTSAD(outfile, view, l_upper, hk_lim, ok)
    9704 *      include 'DIFFaX.par'
    9705 *      include 'DIFFaX.inc'
    9706 **
    9707 *      character*(*) outfile
    9708 *      integer*4 view, hk_lim
    9709 *      real*8 l_upper
    9710 *      logical ok
    9711 **
    9712 *      integer*4 i, j, p1, p2, LENGTH
    9713 *      real*8 rowint(SADSIZE), incr, sigma, x
    9714 **
    9715 ** external function
    9716 *      external LENGTH
    9717 **
    9718 ** external subroutine (Some compilers need them declared external)
    9719 **      external SMUDGE
    9720 **
    9721 *      write(op,102) 'Writing SADP data to file ''',
    9722 *     |  outfile(1:LENGTH(outfile)), '''. . .'
    9723 **
    9724 ** set standard deviation
    9725 *      sigma = ZERO
    9726 **
    9727 ** Establish scaling in reciprocal space, the number of pixels per unit
    9728 ** first get number of pixels per l
    9729 *      incr = dble(SADSIZE/2) / l_upper
    9730 *      if(view.eq.1) then
    9731 ** number of pixels per unit h
    9732 *        incr = incr * sqrt(a0 / c0)
    9733 *      else if(view.eq.2) then
    9734 ** number of pixels per unit k
    9735 *        incr = incr * sqrt(b0 / c0)
    9736 *      else if(view.eq.3) then
    9737 ** number of pixels per unit along (h = k)
    9738 *        incr = incr * sqrt((a0 + b0 + d0) / c0)
    9739 *      else if(view.eq.4) then
    9740 ** number of pixels per unit along (h = -k)
    9741 *        incr = incr * sqrt((a0 + b0 - d0) / c0)
    9742 *      endif
    9743 **
    9744 *      do 10 j = sadblock - 1, 0, -1
    9745 *        do 20 i = 1, SADSIZE
    9746 *          rowint(i) = ZERO
    9747 *   20   continue
    9748 *        do 30 i = 0, hk_lim
    9749 *          p1 = SADSIZE/2 + nint(i*incr) + 1
    9750 *          p2 = SADSIZE/2 - nint(i*incr) + 1
    9751 **
    9752 ** cycle if we have gone out of bounds
    9753 *          if(p1.gt.SADSIZE .or. p1.lt.0) then
    9754 *            goto 30
    9755 *          endif
    9756 *          if(p2.gt.SADSIZE .or. p2.lt.0) then
    9757 *            goto 30
    9758 *          endif
    9759 **
    9760 *          x = spec(i*sadblock + j + 1) * scaleint
    9761 ** handle saturated pixels
    9762 *          if(x.gt.maxsad) x = maxsad
    9763 *          rowint(p1) = x
    9764 *          if(has_l_mirror) then
    9765 *            rowint(p2) = x
    9766 *          else
    9767 *            x = spec(i*sadblock + sadblock - j) * scaleint
    9768 ** handle saturated pixels
    9769 *            if(x.gt.maxsad) x = maxsad
    9770 *            rowint(p2) = x
    9771 *          endif
    9772 *   30   continue
    9773 *        call SMUDGE(rowint, SADSIZE, sigma, ok)
    9774 *        if(.not.ok) goto 999
    9775 *        if(bitdepth.eq.8) then
    9776 *          write(sa,err=900) (char(nint(rowint(i))), i = 1, SADSIZE)
    9777 *        else
    9778 ** Known bug: Writing 16-bit binary is problematic on Intel processors
    9779 *            write(sa,err=900)
    9780 *     | (char(int(rowint(i)/256)),
    9781 *     |  char(int(rowint(i)-(int(rowint(i)/256))*256)),i=1,SADSIZE)
    9782 **            write(sa,err=900) (iidnnt(rowint(i)), i = 1, SADSIZE)
    9783 *        endif
    9784 *   10 continue
    9785 **
    9786 ** Repeat, in reverse for the bottom half if data had a mirror.
    9787 *      if(has_l_mirror) then
    9788 *        do 40 j = 1, sadblock - 1
    9789 *          do 50 i = 1, SADSIZE
    9790 *            rowint(i) = ZERO
    9791 *   50     continue
    9792 *          do 60 i = 0, hk_lim
    9793 *            p1 = SADSIZE/2 + nint(i*incr) + 1
    9794 *            p2 = SADSIZE/2 - nint(i*incr) + 1
    9795 **
    9796 ** cycle if we have gone out of bounds
    9797 *            if(p1.gt.SADSIZE .or. p1.lt.0) then
    9798 *              goto 60
    9799 *            endif
    9800 *            if(p2.gt.SADSIZE .or. p2.lt.0) then
    9801 *              goto 60
    9802 *            endif
    9803 **
    9804 *            x = spec(i*sadblock + j + 1) * scaleint
    9805 ** handle saturated pixels
    9806 *            if(x.gt.maxsad) x = maxsad
    9807 *            rowint(p1) = x
    9808 *            rowint(p2) = x
    9809 *   60     continue
    9810 *          call SMUDGE(rowint, SADSIZE, sigma, ok)
    9811 *          if(.not.ok) goto 999
    9812 *          if(bitdepth.eq.8) then
    9813 *            write(sa,err=900) (char(nint(rowint(i))), i = 1, SADSIZE)
    9814 *          else
    9815 ** Known bug: Writing 16-bit binary is problematic on Intel processors
    9816 *            write(sa,err=900)
    9817 *     | (char(int(rowint(i)/256)),
    9818 *     |  char(int(rowint(i)-(int(rowint(i)/256))*256)),i=1,SADSIZE)
    9819 **            write(sa,err=900) (iidnnt(rowint(i)), i = 1, SADSIZE)
    9820 *          endif
    9821 *   40   continue
    9822 ** write a blank last line to make the array SADSIZE x SADSIZE square
    9823 *        if(bitdepth.eq.8) then
    9824 *          write(sa,err=900) (char(0), i = 1, SADSIZE)
    9825 *        else
    9826 *          write(sa,err=900) (char(0), i = 1, 2*SADSIZE)
    9827 *        endif
    9828 *      endif
    9829 **
    9830 *      if(sa.ne.op) close(unit = sa, err = 990)
    9831 **
    9832 ** MMJT 19th May 2010
    9833 *      write(op,103) SADSIZE, ' x ', SADSIZE, ' pixels: ',
    9834 *     |              bitdepth, ' bits deep.'
    9835 *      write(op,100) ' (Caution: Because of a compiler-dependent issue'
    9836 *      write(op,100) '  the true image width in pixels is found by'
    9837 *      if(bitdepth.eq.8) then
    9838 *        write(op,100) '  dividing the file size in bytes by 256.)'
    9839 *      else
    9840 *        write(op,100) '  dividing the file size in bytes by 512.'
    9841 *        write(op,100) '  The byte-order is big-endian.)'
    9842 *      endif
    9843 **
    9844 *      return
    9845 *  900 write(op,100) 'ERROR: problems writing to binary SADP file.'
    9846 *      ok = .false.
    9847 *      return
    9848 *  990 write(op,100) 'ERROR: problems closing binary SADP file.'
    9849 *      ok = .false.
    9850 *      return
    9851 *  999 write(op,100) 'ERROR: SMUDGE returned error to WRTSAD.'
    9852 *      write(op,101) i, j
    9853 *      return
    9854 *  100 format(1x, a)
    9855 *  101 format(5x, 'with local variables i = ', i5, ' j = ', i5)
    9856 *  102 format(1x, 3a)
    9857 *  103 format(1x, i4, a, i4, a, i2, a)
    9858 *      end
    9859 ** ______________________________________________________________________
    9860 ** Title: WRTSPC
    9861 ** Author: MWD and MMJT
    9862 ** Date: 18 Aug 1988; 7 Mar 1995
    9863 ** Description:  This routine writes the spectrum arrays to file.
    9864 **
    9865 **      ARGUMENTS:
    9866 **            spcfile  -  The name of the output data file. (input).
    9867 **            ok       -  logical flag indicating all went well.
    9868 **                                                      (output).
    9869 **
    9870 **      COMMON VARIABLES:
    9871 **            uses:  th2_min, th2_max, d_theta, spec, brd_spc, blurring
    9872 **                   NONE, RAD2DEG
    9873 **
    9874 **        modifies:  no COMMON variables are modified
    9875 ** ______________________________________________________________________
    9876 **
    9877 *      subroutine WRTSPC(spcfile, ok)
    9878 *      include 'DIFFaX.par'
    9879 *      include 'DIFFaX.inc'
    9880 **
    9881 *      character*(*) spcfile
    9882 *      logical ok
    9883 **
    9884 *      character*1 tab
    9885 *      integer*4 i, n_low, n_high, LENGTH
    9886 *      real*8 theta
    9887 **
    9888 ** external function
    9889 *      external LENGTH
    9890 **
    9891 *      n_low  = 1
    9892 *      n_high = int(HALF*(th2_max-th2_min)/d_theta) + 1
    9893 **
    9894 *      tab = char(9)
    9895 *      theta = th2_min * RAD2DEG
    9896 *      write(op,200) 'Writing spectrum data to file ''',
    9897 *     |      spcfile(1:LENGTH(spcfile)), '''. . .'
    9898 **      do 10 i = int(HALF*th2_min / d_theta) + 1,
    9899 **     |                  int(HALF*th2_max / d_theta) + 1
    9900 *      if(blurring.eq.NONE) then
    9901 *        do 10 i = n_low, n_high
    9902 *          write(sp,101,err=50) theta, tab, spec(i)
    9903 *          theta = theta + TWO * RAD2DEG * d_theta
    9904 *   10   continue
    9905 *      else
    9906 *        do 20 i = n_low, n_high
    9907 *          write(sp,100,err=50) theta, tab, spec(i), tab, brd_spc(i)
    9908 *          theta = theta + TWO * RAD2DEG * d_theta
    9909 *   20   continue
    9910 *      endif
    9911 **
    9912 *      if(sp.ne.op) close(sp,err=60)
    9913 *      write(op,202) 'Spectrum written.'
    9914 *      write(op,202)
    9915 *      return
    9916 *   50 write(op,202) 'Unable to write to spectrum file.'
    9917 *      if(sp.ne.op) close(sp,err=60)
    9918 *      ok = .false.
    9919 *      return
    9920 *   60 write(op,202) 'Unable to close spectrum file.'
    9921 *      ok = .false.
    9922 *      return
    9923 *  100 format(1x, e12.5, 2(a, g13.6))
    9924 *  101 format(1x, e12.5, a, g13.6)
    9925 *  200 format(1x, 3a)
    9926 *  201 format(1x, a, i2, 2a)
    9927 *  202 format(1x, a)
    9928 *      end
    9929 **
    9930 ** ______________________________________________________________________
    99315933* Title: XYPHSE
    99325934* Author: MMJT
Note: See TracChangeset for help on using the changeset viewer.