Changeset 2205
 Timestamp:
 Apr 8, 2016 3:10:04 PM (7 years ago)
 Location:
 trunk
 Files:

 22 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIphsGUI.py
r2204 r2205 2412 2412 'viewPoint':[[0.,0.,0.],[]],} 2413 2413 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])) 2414 2440 2415 2441 def OnLaue(event): … … 2434 2460 G2plt.PlotXYZ(G2frame,XY,Layers['Sadp']['Img'].T,labelX=labels[:1], 2435 2461 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) 2436 2469 2437 2470 def CellSizer(): … … 2895 2928 sadpPlot.Bind(wx.EVT_CHECKBOX,OnSadpPlot) 2896 2929 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) 2897 2934 topSizer.Add(laueSizer,0,WACV) 2898 2935 topSizer.Add(wx.StaticText(layerData,label=' Reference unit cell for all layers:'),0,WACV) … … 2979 3016 return 2980 3017 profile = G2frame.PatternTree.GetItemPyData(G2frame.PatternId)[1] 2981 dlg.Destroy()2982 3018 ctrls = '0\n0\n3\n' 2983 3019 G2pwd.StackSim(data['Layers'],ctrls,scale,background,limits,inst,profile) … … 2987 3023 test2 = np.copy(profile[3]) 2988 3024 rat = (test1test2)/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) 2989 3028 # GSASIIpath.IPyBreak() 2990 3029 G2plt.PlotPatterns(G2frame,plotType='PWDR') … … 2993 3032 data['Layers']['Sadp']['Plane'] = simCodes[1] 2994 3033 data['Layers']['Sadp']['Lmax'] = simCodes[2] 2995 #planeChoice = ['h0l','0kl','hhl','hhl',]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','hhl',] 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) 3000 3039 G2pwd.CalcStackingSADP(data['Layers']) 3001 3040 wx.CallAfter(UpdateLayerData) … … 3003 3042 def OnSeqSimulate(event): 3004 3043 3044 cellSel = ['cellA','cellB','cellC','cellG'] 3045 transSel = ['TransP','TransX','TransY','TransZ'] 3005 3046 ctrls = '' 3047 data['Layers']['seqResults'] = [] 3048 data['Layers']['seqCodes'] = [] 3006 3049 Parms = G2pwd.GetStackParms(data['Layers']) 3007 3050 dlg = G2gd.DIFFaXcontrols(G2frame,ctrls,Parms) 3008 3051 if dlg.ShowModal() == wx.ID_OK: 3009 3052 simCodes = dlg.GetSelection() 3010 data['Layers']['Multi'] = [simCodes[2:5]]3011 3053 else: 3012 3054 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) 3014 3121 3015 3122 ################################################################################ 
trunk/GSASIIplot.py
r2203 r2205 2468 2468 ################################################################################ 2469 2469 2470 def PlotXY(G2frame,XY,XY2=None,labelX=None,labelY=None,newPlot=False,Title='' ):2470 def PlotXY(G2frame,XY,XY2=None,labelX=None,labelY=None,newPlot=False,Title='',lines=False): 2471 2471 '''simple plot of xy data, used for diagnostic purposes 2472 2472 ''' … … 2511 2511 colors=['b','g','r','c','m','k'] 2512 2512 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) 2515 2518 if len(XY2): 2516 2519 for ixy,xy in enumerate(XY2): 2517 X,Y = xy2520 X,Y = XY2[ixy] 2518 2521 Plot.plot(X,Y,colors[ixy%6],picker=False) 2519 2522 if not newPlot: 
trunk/fsource/DIFFaXsubs/DIFFaX.inc
r2191 r2205 375 375 character*4 a_name(MAX_A,MAX_L), atom_l(MAX_TA) 376 376 character*12 pnt_grp 377 character*16 sfname, cfname377 * character*16 sfname, cfname 378 378 * 379 379 logical one_B(MAX_L), one_occup(MAX_L), Bs_zero(MAX_L, MAX_L), … … 416 416 common /chars1/ a_name, atom_l 417 417 common /chars2/ pnt_grp 418 common /chars3/ sfname, cfname418 * common /chars3/ sfname, cfname 419 419 * 420 420 common /logic1/ one_B, one_occup, Bs_zero, there 
trunk/fsource/DIFFaXsubs/DIFFaXsubs.for
r2197 r2205 233 233 100 format(1x, a) 234 234 end 235 **236 ** ______________________________________________________________________237 ** Title: ATOMS238 ** Authors: MMJT239 ** Date: 18 Feb 90240 ** Description: This routine lists the legal atom names accepted by241 ** DIFFaX. These names are taken from the file 'sfname'.242 **243 ** ARGUMENTS:244 ** No arguments are used.245 **246 ** COMMON VARIABLES:247 ** uses: sfname248 **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, LENGTH257 * character*4 atom_name258 * character*80 atomline259 * character*120 line260 **261 ** external function262 * external LENGTH263 **264 ** sfname has been opened once already,265 ** no need for elaborate error checking266 * 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 Hydrogen271 * 1 read(sf,'(a)') line272 * if(line(1:4).ne.'H ') goto 1273 * backspace(unit = sf, err = 500)274 **275 ** initialize atomline276 * atomline = ' '277 ** write 10 atom names to a line278 * 10 i = 0279 * 20 read(sf, '(a)', end = 100) line280 * atom_name = line(1:4)281 * write(atomline(i*7+1:i*7+7),'(3a)') '''', atom_name, ''' '282 * i = i + 1283 * if(mod(i,10).eq.0) then284 * write(op,210) atomline285 * goto 10286 * endif287 * goto 20288 * 100 continue289 **290 * close(unit = sf, err = 600)291 **292 * 999 return293 * 500 write(op,220) 'Unable to backspace scattering factor file ''',294 *  sfname(1:LENGTH(sfname)), '''.'295 * return296 * 600 write(op,220) 'Unable to close scattering factor file ''',297 *  sfname(1:LENGTH(sfname)), '''.'298 * return299 * 200 format(1x, a)300 * 210 format(2x, a)301 * 220 format(1x, 'ERROR in ATOMS: ', 3a)302 * end303 **304 235 * ______________________________________________________________________ 305 236 * Title: BINPOW … … 671 602 205 format(1x, 'Diffraction point symmetry is found to be ''',a,'''') 672 603 300 format(1x, 'Reevaluating diffraction symmetry') 673 end674 *675 * ______________________________________________________________________676 * Title: CHOICE677 * Author: MWD678 * 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 characters681 * 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 none691 *692 integer*4 n693 * character*(*) flag, list(n)694 character*(*) flag695 character*80 list(n)696 *697 integer*4 i, j1, j2, LENGTH698 *699 * external function700 external LENGTH701 *702 i = 1703 j1 = LENGTH(flag)704 *705 10 j2 = LENGTH(list(i))706 * see if the string contained in list(i) is identical to that in flag707 if(j1.eq.j2 .and. index(flag, list(i)(1:j2)).eq.1) goto 20708 i = i + 1709 if(i.le.n) goto 10710 *711 20 if(i.gt.n) i = 1712 CHOICE = i713 *714 return715 604 end 716 605 * … … 1297 1186 * 1298 1187 * ______________________________________________________________________ 1299 * Title: CNTARG1300 * Author: MMJT1301 * Date: 24 Feb 19901302 * Description: Counts the number of arguments in the character1303 * string 'line', and returns the answer in CNTARG. Legal separators1304 * 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 none1315 *1316 character*(*) line1317 *1318 logical in_arg1319 integer*4 blank, tab, comma, i, j, arg_cnt, lin_len1320 parameter (tab=9, blank=32, comma=44)1321 *1322 in_arg = .false.1323 arg_cnt = 01324 CNTARG = 11325 *1326 lin_len = len(line)1327 do 10 i = 1, lin_len1328 j = ichar(line(i:i))1329 if(j.eq.tab .or. j.eq.blank .or. j.eq.comma) then1330 in_arg = .false.1331 else1332 if(.not.in_arg) then1333 in_arg = .true.1334 arg_cnt = arg_cnt + 11335 endif1336 endif1337 10 continue1338 *1339 * There cannot be more than (lin_len+1)/2 arguments1340 if(arg_cnt.le.((lin_len+1)/2)) CNTARG = arg_cnt1341 *1342 return1343 end1344 *1345 * ______________________________________________________________________1346 1188 * Title: DETUN 1347 1189 * Author: MMJT … … 1388 1230 C 401 format(1x, g12.5) 1389 1231 end 1390 *1391 * ______________________________________________________________________1392 ** Title: DUMP1393 ** Author: MWD and MMJT1394 ** Date: 18 Aug 1988; 15 Mar 19951395 ** Description: This subroutine prints out the data read in from1396 ** the data file. It reassures the user that the data was read in1397 ** 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, cntrl1409 ** d_theta, e_sf, FWHM, inf_thick, l_actual1410 ** l_alpha, l_cnt, l_g, l_n_atoms, l_r, l_symmetry1411 ** lambda, l_seq, n_actual, n_atoms, n_layers1412 ** pnt_grp, pv_gamma, pv_u, pv_v, pv_w, r_B111413 ** r_B12, r_B22, r_B23, r_B31, r_B33, rad_type1414 ** recrsv, rndm, th2_max, th2_min1415 ** tolerance, xplcit1416 **1417 ** modifies: no COMMON variables are modified1418 ** ______________________________________________________________________1419 **1420 * subroutine DUMP(infile,ok)1421 * include 'DIFFaX.par'1422 * include 'DIFFaX.inc'1423 **1424 * character*(*) infile1425 * logical ok1426 **1427 * real*8 scale, atom_cnt(MAX_TA), cum_atom_cnt(MAX_TA), norm1428 * integer*4 i, i2, j, n, type, num_types, tot_types1429 * integer*4 LENGTH, print_width1430 * character*31 dmpfile1431 * character*80 list(5)1432 **1433 ** external function1434 * external LENGTH1435 ** external subroutine (Some compilers need them declared external)1436 ** external GETFNM1437 **1438 * call GETFNM(infile, dmpfile, 'dmp', ok)1439 * if(.not.ok) then1440 * write(op,200) 'DUMP aborted'1441 * goto 9991442 * endif1443 * if(dp.ne.op) open(unit = dp , file = dmpfile, status = 'new')1444 **1445 * i2 = 01446 * 99 write(op,200) 'Enter 1 for full atomic position dump'1447 * read(cntrl,*,err=99,end=9999) i21448 * if(CFile) write(op,'(1x,i3)') i21449 * write(op,300) 'Writing data dump to file ''',1450 *  dmpfile(1:LENGTH(dmpfile)),'''. . .'1451 ** sundry details about layers1452 * write(dp,125) 'Number of layers = ', n_layers1453 * write(dp,125) 'Number of unique layers = ', n_actual1454 * write(dp,125) 'Number of different atom types = ', n_atoms1455 ** cell dimensions1456 * write(dp,140)1457 *  ' cell_a, cell_b, cell_c, cell_gamma ',1458 *  cell_a, cell_b, cell_c, RAD2DEG * cell_gamma1459 ** inplane layer widths1460 * if(finite_width) then1461 * if(Wa.lt.inf_width) then1462 * write(dp,126) 'width along a', Wa1463 * write(dp,128) Wa/cell_a1464 * else1465 * write(dp,127) 'a'1466 * endif1467 * if(Wb.lt.inf_width) then1468 * write(dp,126) 'width along b', Wb1469 * write(dp,128) Wb/cell_b1470 * else1471 * write(dp,127) 'b'1472 * endif1473 * else1474 * write(dp,200)1475 *  'Layers are to be treated as having infinite lateral width.'1476 * endif1477 ** radiation type1478 * list(1) = 'XRAY '1479 * list(2) = 'NEUTRON '1480 * list(3) = 'ELECTRON '1481 * write(dp,100) 'radiation type = ', list(rad_type+1)1482 ** wavelength1483 * write(dp,170) 'Wavelength lambda = ', lambda1484 ** symmetry1485 * write(dp,100)1486 * 'Diffraction point group symmetry specified = ',pnt_grp1487 * if(SymGrpNo.eq.UNKNOWN) write(dp,175)1488 * 'Tolerance on intensities in symmetry evaluation = +/',1489 *  tolerance * HUNDRED, ' %'1490 ** instrumental broadening to simulate1491 * list(1) = 'NONE'1492 * list(2) = 'GAUSSIAN'1493 * list(3) = 'LORENTZIAN'1494 * list(4) = 'PSEUDOVOIGT'1495 * i = blurring + 11496 ** see if it's a pure Gaussian pseudoVoigt1497 * if(blurring.eq.PV_GSS) i = 21498 ** see if it's a pure Lorentzian pseudoVoigt1499 * if(blurring.eq.PV_LRN) i = 31500 * 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 halfmaximum = ', FWHM1504 * if(blurring.eq.PS_VGT) then1505 * write(dp,200)1506 * 'PseudoVoigt parameters: u, v, w, gamma'1507 * write(dp,180) pv_u, pv_v, pv_w, pv_gamma1508 * endif1509 * if(blurring.eq.PV_GSS) then1510 * write(dp,200) 'Gaussian parameters: u, v, w'1511 * write(dp,185) pv_u, pv_v, pv_w1512 * endif1513 * if(blurring.eq.PV_LRN) then1514 * write(dp,200) 'Lorentzian parameters: u, v, w'1515 * write(dp,185) pv_u, pv_v, pv_w1516 * endif1517 ** are we to trim out the origin?1518 * if(blurring.ne.NONE) then1519 * if(trim_origin) then1520 * write(dp,200) ' Intensity near origin will be ignored'1521 * else1522 * write(dp,200) ' Intensity near origin will be included'1523 * endif1524 * endif1525 **1526 ** equivalent layers1527 * write(dp,200) ' '1528 * do 10 i = 1, n_layers1529 * 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 continue1533 **1534 * do 11 j = 1, MAX_TA1535 * cum_atom_cnt(j) = ZERO1536 * 11 continue1537 **1538 ** layer structure1539 * write(dp,200) ' '1540 * do 30 i = 1, n_actual1541 * write(dp,130) 'fundamental LAYER ',i1542 * list(1) = 'NONE '1543 * list(2) = 'CENTROSYMMETRIC '1544 * write(dp,100) 'symmetry = ', list(l_symmetry(i)+1)1545 ** write out the layer composition1546 * write(dp,200) 'LAYER composition is:'1547 * do 14 j = 1, MAX_TA1548 * atom_cnt(j) = ZERO1549 * 14 continue1550 * num_types = ZERO1551 * scale = ONE1552 * if(l_symmetry(i).eq.CENTRO) scale = TWO1553 * do 15 j = 1, l_n_atoms(i)1554 * type = a_type(j,i)1555 * if(type.gt.num_types) num_types = type1556 * atom_cnt(type) = atom_cnt(type) + a_occup(j,i) * scale1557 * 15 continue1558 * if(num_types.gt.tot_types) tot_types = num_types1559 ** accumulate weighted atom count1560 * do 16 j = 1, MAX_TA1561 * cum_atom_cnt(j) = cum_atom_cnt(j) + atom_cnt(j) * l_g(i)1562 * 16 continue1563 * do 17 j = 1, num_types1564 * write(dp,280) atom_l(j), ':', atom_cnt(j)1565 * 17 continue1566 * write(dp,200) ' '1567 **1568 ** do we want all the details about each atom?1569 * if(i2.ne.1) goto 301570 ** 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) then1579 * write(dp,140) 'Xray scattering factor data Ai, Bi',1580 *  (x_sf(n,a_type(j,i)),n=1,9)1581 * else if(rad_type.eq.NEUTRN) then1582 * write(dp,140) 'Neutron scattering factor data',1583 *  n_sf(a_type(j,i))1584 * else if(rad_type.eq.ELECTN) then1585 * 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 * else1588 * write(op,200) 'ERROR: Illegal radiation type in DUMP'1589 * endif1590 * 20 continue1591 * write(dp,200) ' '1592 * 30 continue1593 **1594 ** Print out average composition1595 * if(.not.xplcit) then1596 * norm = ZERO1597 * do 25 j = 1, num_types1598 * norm = norm + cum_atom_cnt(j)1599 * 25 continue1600 * write(dp,200) ' '1601 * write(dp,200) 'Average crystal composition is:'1602 * do 26 j = 1, num_types1603 * write(dp,281) atom_l(j), ':', cum_atom_cnt(j) / norm1604 * 26 continue1605 * write(dp,200) ' '1606 * endif1607 **1608 ** stacking details1609 * write(dp,200) ' '1610 * write(dp,200) 'STACKING'1611 * if(recrsv) then1612 * write(dp,210) 'RECURSIVELY'1613 * if(inf_thick) then1614 * write(dp,220) 'INFINITE'1615 * else1616 * write(dp,230) l_cnt1617 * endif1618 * else if(xplcit) then1619 * if(rndm) then1620 * write(dp,240)1621 * else1622 * write(dp,250)1623 * endif1624 * write(dp,260) 'Sequence for ', l_cnt, ' layers is:'1625 * print_width = 251626 * j = 11627 * n = print_width1628 * 35 if(n.gt.l_cnt) n = l_cnt1629 * write(dp,270) (l_seq(i), i = j, n)1630 * if(n.lt.l_cnt) then1631 * j = j + print_width1632 * n = n + print_width1633 * goto 351634 * endif1635 * endif1636 **1637 ** stacking transition data1638 * 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_layers1648 * write(dp,200)1649 * ' '1650 * do 50 j = 1, n_layers1651 * 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 continue1654 * 40 continue1655 **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 > C131660 * write(dp,200)1661 * ' i j  C11 C22 C33 C12 C13 C23'1662 * write(dp,200)1663 * ' '1664 * do 60 i = 1, n_layers1665 * write(dp,200)1666 * ' '1667 * do 70 j = 1, n_layers1668 ** MMJT: 3/15/95. Ordering of B23 and B31 swapped1669 * 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 continue1672 * 60 continue1673 **1674 * if(dp.ne.op) close(unit = dp)1675 * 999 return1676 * 9999 ok = .false.1677 * write(op,200) 'DUMP aborted'1678 * return1679 * 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 * end1708 1232 * 1709 1233 * ______________________________________________________________________ … … 1767 1291 * 1768 1292 * ______________________________________________________________________ 1769 ** Title: FNDCTL1770 ** Author: MMJT1771 ** Date: 6 June 19891772 ** Description: Checks whether or not there is a control file, and1773 ** 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, CFile1780 **1781 ** modifies: cntrl1782 **1783 ** FNDCTL returns logical .true. if a control file was found in the1784 ** current directory.1785 ** ______________________________________________________________________1786 **1787 * logical function FNDCTL()1788 * include 'DIFFaX.par'1789 * include 'DIFFaX.inc'1790 **1791 * integer*4 io_err, LENGTH1792 **1793 ** external function1794 * external LENGTH1795 **1796 ** See if control file is present1797 * inquire(file = cfname, exist = CFile)1798 ** get data from ControlFile1799 * if(CFile) then1800 * cntrl = 991801 * open(unit = cntrl, file = cfname, status = 'old',1802 *  err = 999, iostat = io_err)1803 * else1804 ** get data from user entries to the screen1805 * cntrl = ip1806 * endif1807 * FNDCTL = .true.1808 **1809 * return1810 * 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_err1814 * FNDCTL = .false.1815 * return1816 * 100 format(1x, 3a)1817 * 101 format(1x, a)1818 * 102 format(1x, a, i5)1819 * end1820 **1821 ** ______________________________________________________________________1822 ** Title: GETFIL1823 ** Author: MMJT1824 ** Date: 5 June 19891825 ** Description: Checks for the presence of the structure data file and1826 ** 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 read1833 ** from the control file. (output).1834 **1835 ** COMMON VARIABLES:1836 ** uses: cntrl, sfname, CFile1837 **1838 ** modifies: no COMMON variables are modified1839 ** ______________________________________________________________________1840 **1841 * subroutine GETFIL(fname, ok, ending)1842 * include 'DIFFaX.par'1843 * include 'DIFFaX.inc'1844 **1845 * character*(*) fname1846 * logical ok, ending1847 **1848 * character*31 tmp_nam1849 * character*200 line1850 * character*80 list(1)1851 * integer*4 unit_no, LENGTH, CHOICE1852 **1853 * integer*4 i, io_err1854 **1855 ** external function1856 * external LENGTH, CHOICE1857 ** external subroutines (Some compilers need them declared external)1858 ** external GETLNE, TOUPPR1859 **1860 * unit_no = cntrl1861 ** get file name from ControlFile1862 * 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) fname1867 ** check that the string does not say 'end'. First convert to upper case1868 ** we must be sure that the string 'end' is not part of a legitimate1869 ** file name. Use CHOICE for this.1870 * tmp_nam = fname1871 * call TOUPPR(fname)1872 * list(1) = 'END '1873 * i = CHOICE(fname, list, 1)1874 * if(i.eq.1) then1875 * ending = .true.1876 * goto 9991877 * endif1878 * fname = tmp_nam1879 **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 9991898 * 991 write(op,100)1899 * 'Problems opening structure data file ''',1900 *  fname(1:LENGTH(fname)),''''1901 * write(op,103) 'IOSTAT = ',io_err1902 * ok = .false.1903 * goto 9991904 * 992 write(op,100)1905 * 'Problems opening scattering factor file ''',1906 *  sfname(1:LENGTH(sfname)),''''1907 * write(op,103) 'IOSTAT = ',io_err1908 * ok = .false.1909 * goto 9991910 * 993 write(op,104) 'ERROR reading filename in GETFIL'1911 * ok = .false.1912 * 999 return1913 * 100 format(1x, 3a)1914 * 103 format(1x, a, i5)1915 * 104 format(1x, a)1916 * end1917 **1918 ** ______________________________________________________________________1919 ** Title: GETFNM1920 ** Author: MMJT1921 ** Date: 13 Oct 19881922 ** Description: This subroutine makes sure that we do not write over1923 ** an existing file, and creates a suitable filename (name2). name21924 ** equals 'name1.append' append is supplied by the calling routine.1925 ** If 'name1.append' is the name of an existing file, then append is1926 ** set to 'append#' where #=1,2.1927 ** CLIP is the maximum length of the appendage. Thus if append1928 ** = '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 are1932 ** 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 ** save1946 **1947 * character*(*) name1, name2, append1948 * logical ok1949 **1950 * logical fexist1951 * character*31 appnd21952 * integer*4 applen, i, myclip1953 * integer*4 LENGTH1954 **1955 ** external function1956 * external LENGTH1957 ** external subroutine (Some compilers need them declared external)1958 ** external NAMER1959 **1960 * ok = .true.1961 * fexist = .true.1962 * appnd2 = append1963 * applen = LENGTH(append)1964 * myclip = CLIP1965 * if(myclip.le.0) then1966 * write(op,302) 'Illegal CLIP value ', CLIP1967 * write(op,300) 'Default CLIP value of 3 is being assumed'1968 * myclip = 31969 * endif1970 **1971 * i = 01972 * 10 call NAMER(name1,name2,appnd2)1973 * inquire(file = name2, exist = fexist)1974 * if(fexist) then1975 * i = i + 11976 * if(i.lt.10) then1977 * if(myclip.gt.10) then1978 * write(appnd2(applen+1:applen+1),200) i1979 * else1980 * write(appnd2(myclip:myclip),200) i1981 * endif1982 * else if(i.ge.10 .and. i.lt.100) then1983 * if(myclip.gt.10) then1984 * write(appnd2(applen+1:applen+2),201) i1985 * else if(myclip.gt.1) then1986 * write(appnd2(myclip1:myclip),201) i1987 * else1988 * ok = .false.1989 * goto 9991990 * endif1991 * else if(i.ge.100 .and. i.lt.1000) then1992 * if(myclip.gt.10) then1993 * write(appnd2(applen+1:applen+3),202) i1994 * else if(myclip.gt.2) then1995 * write(appnd2(myclip2:myclip),202) i1996 * else1997 * ok = .false.1998 * goto 9991999 * endif2000 * else2001 * ok = .false.2002 * goto 9992003 * endif2004 * goto 102005 * endif2006 **2007 * return2008 * 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 * return2013 * 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 * end2020 **2021 ** ______________________________________________________________________2022 1293 * Title: GETLAY 2023 1294 * Author: MMJT … … 2105 1376 * 2106 1377 * ______________________________________________________________________ 2107 ** Title: GETLNE2108 ** Author: MWD & MMJT2109 ** Date: MWD: Original version 18 Aug 1988; MMJT: This version 9 Apr 19922110 ** Description: This function returns a nonblank, leftjustified line,2111 ** stripped of any comments, from unit 'unit_no'. Any error generates2112 ** a return to the passed linenumber, presumably an errorhandling2113 ** routine in the calling procedure. Multiple comments, and nested2114 ** braces on one line can be handled.2115 **2116 ** ARGUMENTS:2117 ** unit_no  Logical unit number that data is to2118 ** 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 routine2122 ** to return to if an error occurs. (input).2123 ** ______________________________________________________________________2124 **2125 * subroutine GETLNE(unit_no, line, *)2126 * include 'DIFFaX.par'2127 ** save2128 **2129 * integer*4 unit_no2130 * character*(*) line2131 **2132 * character*1 tab, blank2133 ** nominal line length is 200. Leave generous room for overflow.2134 * character*250 rawline2135 * integer*4 i, j, n, offset, lin_len2136 * integer*4 l_br, r_br, first_left, last_right2137 **2138 * tab = char(9)2139 * blank = char(32)2140 * lin_len = len(line)2141 **2142 * 1 read(unit_no, '(a)', err = 200, end = 200) rawline2143 * line = rawline2144 **2145 ** First check if user has gone beyond lin_len characters2146 * do 2 i = len(rawline), lin_len + 1, 12147 * if(rawline(i:i).ne.blank) goto 9982148 * 2 continue2149 **2150 ** strip out all tabs. Replace with blanks.2151 * 3 n = index(line,tab)2152 * if(n.eq.0) goto 42153 * write(line(n:n),'(a)') blank2154 * goto 32155 * 4 continue2156 **2157 ** Search for matching braces. Keep looping until all braces are found2158 * 10 first_left = 02159 * last_right = 02160 * l_br = 02161 * r_br = 02162 * i = 02163 * 20 i = i + 12164 * if(i.gt.lin_len) goto 302165 * if(line(i:i).eq.'{') then2166 * if(first_left.eq.0) first_left = i2167 * l_br = l_br + 12168 * else if(line(i:i).eq.'}') then2169 * last_right = i2170 * r_br = r_br + 12171 * endif2172 * if(r_br.gt.l_br) goto 9992173 * if(l_br.eq.0 .or. l_br.gt.r_br) goto 202174 * 30 continue2175 * if(r_br.ne.l_br) goto 9992176 **2177 ** if no more comments, then left justify the line2178 * if(l_br.eq.0) goto 502179 ** else, remove comments, and shunt data to the left2180 * offset = last_right  first_left2181 * i = first_left2182 * do 60 j = last_right + 1, lin_len2183 * line(i:i) = line(j:j)2184 * i = i + 12185 * 60 continue2186 ** blank out the end of the modified string2187 * do 70 i = lin_lenoffset, lin_len2188 * line(i:i) = blank2189 * 70 continue2190 ** repeat until no more braces are found2191 * goto 102192 **2193 ** left adjust line2194 * 50 do 80 i = 1, lin_len2195 * if((line(i:i).ne.blank)) then2196 * do 90 j = 1, lin_leni+12197 * line(j:j) = line(i+j1:i+j1)2198 * 90 continue2199 * do 100 j = lin_leni+2,lin_len2200 * line(j:j) = blank2201 * 100 continue2202 * goto 1102203 * endif2204 * 80 continue2205 ** If, after all this, the line is empty, repeat with next next line.2206 * 110 if(line(1:1).eq.blank) goto 12207 **2208 * return2209 **2210 * 998 write(op,401) 'ERROR: Data line exceeds ',lin_len, ' characters.'2211 * return 12212 * 999 write(op,400) 'ERROR: Unbalanced braces in data file.'2213 * 200 return 12214 **2215 * 400 format(1x, a)2216 * 401 format(1x, a, i3, a)2217 * end2218 **2219 ** ______________________________________________________________________2220 1378 * Title: GETSAD 2221 1379 * Author: MMJT … … 3986 3144 * 3987 3145 * ______________________________________________________________________ 3988 ** Title: GOINTR3989 ** Author: MMJT3990 ** Date: 23 Oct 19893991 ** Description: This subroutine sets up the parameters necessary to3992 ** 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, CFile3999 **4000 ** modifies: no COMMON variables are modified4001 ** ______________________________________________________________________4002 **4003 * subroutine GOINTR(ok)4004 * include 'DIFFaX.par'4005 * include 'DIFFaX.inc'4006 **4007 * logical ok4008 **4009 * integer*4 i4010 * real*8 AGLQ16, GLQ164011 **4012 ** external functions4013 * external AGLQ16, GLQ164014 ** external subroutine (Some compilers need them declared external)4015 ** external INTEGR4016 **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) i4021 * if(CFile) write(op,101) i4022 * if(i.eq.1) then4023 * call INTEGR(AGLQ16,ok)4024 * else4025 * call INTEGR(GLQ16,ok)4026 * endif4027 * if(.not.ok) goto 9994028 * i = 14029 * 2 write(op,100) 'Enter 1 to integrate another interval.'4030 * read(cntrl,*,err=2) i4031 * if(CFile) write(op,101) i4032 * if(i.eq.1) goto 14033 **4034 * 999 return4035 * 100 format(1x, a)4036 * 101 format(1x, i3)4037 * end4038 **4039 ** ______________________________________________________________________4040 ** Title: GOSADP4041 ** Author: MMJT4042 ** Date: 23 Oct 1989; 1st Mar 20004043 ** Description: This subroutine sets up a file whose name is given4044 ** by 'outfile', which is derived from the input filename 'infile'.4045 ** The selected area diffraction pattern (SADP) is calculated and is4046 ** 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, CFile4056 **4057 ** modifies: loglin, brightness4058 ** ______________________________________________________________________4059 **4060 * subroutine GOSADP(infile,outfile,ok)4061 * include 'DIFFaX.par'4062 * include 'DIFFaX.inc'4063 **4064 * character*(*) infile, outfile4065 * logical ok4066 **4067 * integer*4 i_adapt, i_plane, io_err, hk_lim4068 * real*8 AGLQ16, GLQ16, l_upper4069 **4070 ** external functions4071 * external AGLQ16, GLQ164072 ** external subroutines (Some compilers need them declared external)4073 ** external STREAK, GETFNM, GETSAD, WRTSAD4074 **4075 ** Get a suitable file name, and open file.4076 * call GETFNM(infile,outfile,'sadp',ok)4077 * if(.not.ok) goto 9994078 ** 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 output4082 **  ,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_adapt4089 * if(CFile) write(op,101) i_adapt4090 **4091 * 1 write(op,100) 'Choose a plane in reciprocal space to view.'4092 * write(op,100) ' the laxis 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_plane4095 * if(CFile) write(op,101) i_plane4096 * if(i_plane.lt.1 .or. i_plane.gt.4) then4097 * if(CFile) then4098 * write(op,100) 'Illegal choice. Default is k = 0.'4099 * i_plane = 14100 * else4101 * write(op,100) 'Illegal choice. Choose again.'4102 * goto 14103 * endif4104 * endif4105 **4106 ** Get upper bounds. The SADP will be square, so we only need one bound.4107 ** Choose laxis, since this is common to all 4 views and makes scaling4108 ** 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_upper4111 * if(CFile) write(op,104) l_upper4112 **4113 ** 8bit images or 16bit images?4114 * 3456 write(op,100) 'Choose the bitdepth of the image.'4115 * write(op,100) ' 8:  8bits 16:  16bits'4116 * read(cntrl,*,err=3456) bitdepth4117 * if(CFile) write(op,101) bitdepth4118 * if(bitdepth.ne.8 .and. bitdepth.ne.16) then4119 * write(op,100) 'Illegal bitdepth.'4120 * if(CFile) then4121 * write(op,100) 'Using 8bit as the default'4122 * bitdepth = 84123 * else4124 * write(op,100) 'Reenter. . .'4125 * goto 34564126 * endif4127 * endif4128 * if(bitdepth.eq.16) then4129 ** Bypass issue of signed or unsigned format. Use range 0  32767 only.4130 * write(op,100)4131 *  'File will be saved in unsigned 16bit (bigendian) format.'4132 * maxsad = FIFTEENBITS  ONE4133 * else4134 * write(op,100) 'File will be saved in unsigned 8bit format.'4135 * maxsad = EIGHTBITS  ONE4136 * endif4137 **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) loglin4142 * if(CFile) write(op,101) loglin4143 * if(loglin.lt.0 .or. loglin.gt.1) then4144 * write(op,100) 'Illegal intensity scaling type.'4145 * if(CFile) then4146 * write(op,100) 'Using linear as the default'4147 * loglin = 14148 * else4149 * write(op,100) 'Reenter. . .'4150 * goto 24151 * endif4152 * endif4153 **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) then4157 * 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 * else4161 * 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 * endif4165 * read(cntrl,*,err=3) brightness4166 * if(CFile) write(op,104) brightness4167 * if(brightness.le.ZERO) then4168 * write(op,100) 'Illegal value for brightness. Must be positive'4169 * if(CFile) then4170 * if(loglin.eq.0) then4171 * write(op,100) 'Using default of 1 for logarithmic scaling'4172 * brightness = ONE4173 * else4174 * write(op,100) 'Using default of 10 for linear scaling'4175 * brightness = TEN4176 * endif4177 * else4178 * if(loglin.eq.0) then4179 * write(op,100) '1 is a good default for logarithmic scaling'4180 * else4181 * write(op,100) '10 is a good default for linear scaling'4182 * endif4183 * write(op,100) 'Reenter. . .'4184 * goto 34185 * endif4186 * endif4187 **4188 ** Generate intensity data for the SAD pattern.4189 * if(i_adapt.eq.1) then4190 * call GETSAD(AGLQ16, i_plane, l_upper, hk_lim, infile, ok)4191 * else4192 * call GETSAD( GLQ16, i_plane, l_upper, hk_lim, infile, ok)4193 * endif4194 **4195 * if(ok) call WRTSAD(outfile, i_plane, l_upper, hk_lim, ok)4196 **4197 * 999 return4198 * 990 write(op,102) 'Problems opening output file ', outfile4199 * write(op,103) 'IOSTAT = ', io_err4200 * ok = .false.4201 * return4202 * 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 * end4208 **4209 ** ______________________________________________________________________4210 ** Title: GOSPEC4211 ** Author: MMJT4212 ** Date: 23 Oct 19894213 ** Description: This subroutine sets up a file whose name is given4214 ** by 'outfile', which is derived from the input filename 'infile'.4215 ** The powder pattern data is then written to 'outfile', which is closed4216 ** 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_VGT4226 ** PV_GSS, PV_LRN, cntrl4227 **4228 ** modifies: full_brd4229 ** ______________________________________________________________________4230 **4231 * subroutine GOSPEC(infile,outfile,ok)4232 * include 'DIFFaX.par'4233 * include 'DIFFaX.inc'4234 **4235 * character*(*) infile, outfile4236 * logical ok4237 **4238 * logical GETSPC, TRMSPC, RDRNGE4239 * integer*4 io_err4240 * real*8 AGLQ16, GLQ16, cut_off4241 **4242 ** external functions4243 * external RDRNGE, GETSPC, AGLQ16, GLQ16, TRMSPC4244 ** external subroutines (Some compilers need them declared external)4245 ** external GETFNM, PV, GAUSSN, LORNZN, WRTSPC4246 **4247 * write(op,100) ' '4248 * write(op,100) 'CALCULATING POWDER DIFFRACTION PATTERN. . .'4249 **4250 ** get angular range and step4251 * ok = RDRNGE()4252 * if(.not.ok) goto 9994253 ** create a new filename to save spectrum data to4254 * call GETFNM(infile, outfile, 'spc', ok)4255 * if(.not.ok) goto 9994256 * if(sp.ne.op) open(unit = sp, file = outfile, status = 'new',4257 *  err = 990, iostat = io_err)4258 * full_brd = 14259 * 3456 write(op,100) 'Enter 1 for adaptive quadrature on broad peaks.'4260 * read(cntrl,*,err=3456) full_brd4261 * if(CFile) write(op,101) full_brd4262 * if(full_brd.eq.1) then4263 * ok = GETSPC(AGLQ16, infile)4264 * else4265 * ok = GETSPC(GLQ16, infile)4266 * endif4267 ** suppress the huge peak near the origin if required4268 * cut_off = ZERO4269 * if(ok.and.(th2_min.eq.ZERO).and.trim_origin) ok = TRMSPC(cut_off)4270 * if(ok) then4271 * if(blurring.eq.GAUSS) then4272 * write(op,104) 'Gaussian'4273 * call GAUSSN(cut_off)4274 * else if(blurring.eq.LORENZ) then4275 * write(op,104) 'Lorentzian'4276 * call LORNZN(cut_off)4277 * else if(blurring.eq.PS_VGT) then4278 * write(op,104) 'pseudoVoigt'4279 * call PV(cut_off)4280 * else if(blurring.eq.PV_GSS) then4281 * write(op,104) 'Gaussian'4282 * call PV(cut_off)4283 * else if(blurring.eq.PV_LRN) then4284 * write(op,104) 'Lorentzian'4285 * call PV(cut_off)4286 * else if(blurring.ne.NONE) then4287 * write(op,100)4288 *  'Instrumental broadening type is undefined in GOSPEC.'4289 * endif4290 * endif4291 * if(ok) call WRTSPC(outfile, ok)4292 **4293 * 999 return4294 * 990 write(op,102) 'Problems opening output file ', outfile4295 * write(op,103) 'IOSTAT = ', io_err4296 * ok = .false.4297 * return4298 * 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 * end4304 **4305 ** ______________________________________________________________________4306 ** Title: GOSTRK4307 ** Author: MMJT4308 ** Date: 23 Oct 19894309 ** Description: This subroutine sets up a file whose name is given4310 ** by 'outfile', which is derived from the input filename 'infile'.4311 ** The streak intensity is then written to 'outfile', which is closed4312 ** 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, CFile4322 **4323 ** modifies: no COMMON variables are modified4324 ** ______________________________________________________________________4325 **4326 * subroutine GOSTRK(infile,outfile,ok)4327 * include 'DIFFaX.par'4328 * include 'DIFFaX.inc'4329 **4330 * character*(*) infile, outfile4331 * logical ok4332 **4333 * integer*4 i, io_err4334 * real*8 AGLQ16, GLQ164335 **4336 ** external functions4337 * external AGLQ16, GLQ164338 ** external subroutines (Some compilers need them declared external)4339 ** external STREAK, GETFNM4340 **4341 * call GETFNM(infile, outfile, 'str', ok)4342 * if(.not.ok) goto 9994343 * 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) i4349 * if(CFile) write(op,101) i4350 ** select integration function4351 * if(i.eq.1) then4352 * call STREAK(AGLQ16, outfile, ok)4353 * else4354 * call STREAK( GLQ16, outfile, ok)4355 * endif4356 **4357 * 999 return4358 * 990 write(op,102) 'Problems opening output file ', outfile4359 * write(op,103) 'IOSTAT = ', io_err4360 * ok = .false.4361 * return4362 * 100 format(1x, a)4363 * 101 format(1x, i3)4364 * 102 format(1x, 2a)4365 * 103 format(1x, a, i5)4366 * end4367 **4368 ** ______________________________________________________________________4369 3146 * Title: HKL_LIM 4370 3147 * Author: MMJT … … 4733 3510 return 4734 3511 end 4735 *4736 * ______________________________________________________________________4737 ** Title: LAYCNT4738 ** Author: MMJT4739 ** Date: 3 Oct 19894740 ** Description: Counts the number of integer arguments in the character4741 ** string 'line', and returns the answer in LAYCNT. Legal separators4742 ** are, blanks, tabs, and commas. If LAYCNT encounters anything other4743 ** than a digit (0  9), or one of the three legal separators, LAYCNT4744 ** is set to 1. This routine is similar to CNTARG, accept it will only4745 ** 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 none4756 **4757 * character*(*) line4758 **4759 * logical in_num, legit4760 * integer*4 blank, tab, comma, nought, nine, i, j, num_cnt, lin_len4761 * parameter (tab=9, blank=32, comma=44, nought=48, nine=57)4762 **4763 * in_num = .false.4764 * legit = .true.4765 * num_cnt = 04766 * LAYCNT = 14767 **4768 * lin_len = len(line)4769 * do 10 i = 1, lin_len4770 * j = ichar(line(i:i))4771 * if(j.eq.tab .or. j.eq.blank .or. j.eq.comma) then4772 * in_num = .false.4773 * else if(j.ge.nought .and. j.le.nine) then4774 * if(.not.in_num) then4775 * in_num = .true.4776 * num_cnt = num_cnt + 14777 * endif4778 * else4779 ** illegal character4780 * legit = .false.4781 * endif4782 * 10 continue4783 **4784 * if(legit) LAYCNT = num_cnt4785 **4786 * return4787 * end4788 3512 * 4789 3513 * ______________________________________________________________________ … … 5441 4165 * 5442 4166 * ______________________________________________________________________ 5443 ** Title: NAMER5444 ** Author: MMJT5445 ** Date: 13 Oct 19885446 ** Description: This subroutine creates the appropriate filenames for5447 ** 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 within5450 ** the name (ie. if the data file is called 'Structure.dat'). It deletes5451 ** 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 GETFNM5454 **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 ** save5464 **5465 * character*(*) name1, name2, append5466 **5467 * integer*4 applen, namlen, i, idot5468 * integer*4 LENGTH5469 **5470 ** external function5471 * external LENGTH5472 **5473 ** get length of the string holding the filename5474 * 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) then5481 * if(namlen.lt.MAX_NAM) then5482 * idot = namlen + 15483 * else5484 * idot = MAX_NAM5485 * namlen = MAX_NAM5486 * endif5487 * endif5488 ** truncate root filename so that the appendage will appear correctly5489 * if(idot+applen.gt.namlen) idot = namlen  applen5490 **5491 * do 10 i = 1, idot15492 * write(name2(i:i),'(a)') name1(i:i)5493 * 10 continue5494 **5495 * write(name2(idot:idot),'(a)') '.'5496 **5497 * do 20 i = 1, applen5498 * write(name2((idot+i):(idot+i)),'(a)') append(i:i)5499 * 20 continue5500 **5501 * do 30 i = idot+applen+1, namlen5502 * write(name2(i:i),'(a)') ' '5503 * 30 continue5504 **5505 * return5506 * end5507 **5508 ** ______________________________________________________________________5509 4167 * Title: NMCOOR 5510 4168 * Author: MWD … … 5540 4198 * 5541 4199 * ______________________________________________________________________ 5542 ** Title: NXTARG5543 ** Author: MMJT5544 ** Date: 21 JUL 19975545 ** Description: Removes the next argument from the character string5546 ** '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 none5559 * save5560 **5561 * character*(*) line, arg5562 **5563 * character*200 tmpline5564 * integer*4 blank, tab, comma, i, j, i1, i2, lin_len5565 * parameter (tab=9, blank=32, comma=44)5566 **5567 * NXTARG = 05568 **5569 * lin_len = len(line)5570 **5571 * i = 05572 ** Strip leading blanks5573 * 10 i = i + 15574 * if(i.gt.lin_len) goto 9975575 * j = ichar(line(i:i))5576 * if(j.eq.tab .or. j.eq.blank .or. j.eq.comma) goto 105577 * i1 = i5578 **5579 ** Read off the nonblank characters5580 * 20 i = i + 15581 * if(i.gt.lin_len) goto 305582 * j = ichar(line(i:i))5583 * if(j.eq.tab .or. j.eq.blank .or. j.eq.comma) then5584 * NXTARG = 15585 * goto 305586 * else5587 * goto 205588 * endif5589 * 30 continue5590 **5591 * i2 = i5592 * if(i2.lt.i1) goto 9985593 * 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 * return5598 * 997 NXTARG = 35599 * return5600 * 998 NXTARG = 25601 * return5602 * 999 NXTARG = 15603 * return5604 * end5605 **5606 ** ______________________________________________________________________5607 4200 * Title: OPTIMZ 5608 4201 * Author: MWD and MMJT … … 6086 4679 * 6087 4680 * ______________________________________________________________________ 6088 ** Title: PARENT6089 ** Author: MMJT6090 ** Date: 10 Aug 19896091 ** Description: This routine detects whether or not a character6092 ** string contains any parentheses, '(' and ')'. If it does, they are6093 ** removed. PARENT is used by RDPRMS to detect whether or not the6094 ** 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 balanced6101 ** pair of parentheses were found, and 1 if there was an error.6102 ** ______________________________________________________________________6103 **6104 * integer*4 function PARENT(line)6105 * implicit none6106 **6107 * character*(*) line6108 **6109 * integer*4 i, j6110 **6111 * PARENT = 06112 * i = index(line,'(')6113 * j = index(line,')')6114 ** no parentheses6115 * if(i.eq.0 .and. j.eq.0) goto 2006116 ** right hand parenthesis occurred before the left hand parenthesis!6117 * if(j.lt.i) goto 9996118 ** only one parenthesis6119 * if(i.ne.j .and. (i.eq.0 .or. j.eq.0)) goto 9996120 ** blot out the two parentheses found6121 * write(line(i:i),'(a)') ' '6122 * write(line(j:j),'(a)') ' '6123 ** check there are not any more6124 * i = index(line,'(')6125 * j = index(line,')')6126 * if(i.eq.0 .and. j.eq.0) then6127 ** there were no more parentheses. Data is probably ok.6128 * PARENT = 16129 * else6130 ** there were more parentheses. Data is suspect.6131 * goto 9996132 * endif6133 **6134 * 200 return6135 * 999 PARENT = 16136 * return6137 * end6138 **6139 ** ______________________________________________________________________6140 4681 * Title: PNTINT 6141 4682 * Author: MMJT … … 6208 4749 * 6209 4750 * ______________________________________________________________________ 6210 ** Title: POINT6211 ** Author: MWD and MMJT6212 ** Date: 18 Aug 19886213 ** Description: This subroutine prompts the user for h, k, l values6214 ** 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, wavefn6222 **6223 ** modifies: no COMMON variables are modified6224 ** ______________________________________________________________________6225 **6226 * subroutine POINT(ok)6227 * include 'DIFFaX.par'6228 * include 'DIFFaX.inc'6229 **6230 * logical ok6231 **6232 * logical GET_S, GET_S26233 * integer*4 h, k, i6234 * character chr6235 * real*8 l, INTENS, INTEN2, Shkl, x, ANGLE, SS, W4, theta, Q26236 * complex*16 f(MAX_L), s(MAX_L)6237 **6238 ** external functions6239 * external INTENS, INTEN2, GET_S, GET_S26240 ** external subroutines (Some compilers need them declared external)6241 ** external XYPHSE, PRE_MAT, GET_MAT, GET_F6242 **6243 ** statement functions6244 ** SS is the value of 1/d**2 at hkl6245 * SS(h,k,l) = h*h*a0 + k*k*b0 + l*l*c0 + h*k*d06246 ** ANGLE is the Bragg angle (in radians) of the h,k,l plane6247 * ANGLE(h,k,l) = asin(HALF * lambda * sqrt(SS(h,k,l)))6248 ** W4 is the Xray polarization factor6249 * 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, l6255 * if(CFile) write(op,401) h, k, l6256 ** check angles are meaningful6257 * Q2 = FOUR / (lambda**2)6258 * Shkl = SS(h,k,l)6259 ** Check that SS(h,k,l) is legal6260 * if(Shkl.lt.ZERO) then6261 * write(op,403) 'ERROR: In POINT(), 1/d**2 = ', Shkl6262 * ok = .false.6263 * goto 9996264 * endif6265 * if(Shkl.gt.Q2) then6266 * write(op,402) h, k, l,' exceeds 180 degree scattering angle!'6267 * write(op,400) 'Reenter. . .'6268 * goto 106269 * endif6270 ** make sure we are not going to blow up at the origin6271 * if(h.eq.0 .and. k.eq.0 .and. rad_type.eq.ELECTN) then6272 * if(Shkl.le.eps4) then6273 * write(op,400)6274 *  'Cannot integrate across the origin for electron radiation'6275 * write(op,400) 'Reenter. . .'6276 * goto 106277 * endif6278 * endif6279 ** get intensity6280 * call XYPHSE(h, k)6281 * call PRE_MAT(h, k)6282 * call GET_F(f, Shkl, l)6283 * if(recrsv) then6284 * x = INTENS(f, h, k, l, ok)6285 * if(.not.ok) goto 9996286 ** set up mat again to recall GET_S (or GET_S2)6287 * call PRE_MAT(h, k)6288 * call GET_MAT(h, k, l)6289 * if(inf_thick) then6290 ** initialize s to f, since mat is (ident  T)6291 * do 20 i = 1, n_layers6292 * s(i) =  f(i)6293 * 20 continue6294 * ok = GET_S(f, s, h, k, l)6295 * else6296 ** s is initialized within GET_S2, which then calls GET_S6297 * ok = GET_S2(f, s, h, k, l)6298 * endif6299 * else6300 * x = INTEN2(f, h, k, l, ok)6301 * endif6302 * if(.not.ok) goto 9996303 ** for diagnostics write out the s values as well6304 * 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_layers6312 * write(op,405) i, f(i)6313 * 30 continue6314 * write(op,400) ' '6315 * if(recrsv) then6316 * write(op,400) 'Average scattered wavefunctions'6317 * do 40 i = 1, n_layers6318 * write(op,407) i, s(i)6319 * 40 continue6320 * else6321 * write(op,408) 'Crystal wavefunction', wavefn6322 * endif6323 * write(op,400) ' '6324 * write(op,406) 'Intensity at ', h, k, l, ' = ', x6325 * write(op,400) 'Another point? (y/n)'6326 * read(cntrl,'(a)') chr6327 * if(chr(1:1).eq.'y' .or. chr(1:1).eq.'Y') goto 16328 **6329 * 999 return6330 * 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 * end6340 **6341 ** ______________________________________________________________________6342 4751 * Title: POLINT 6343 4752 * Author: MMJT, adapted from Numerical Recipes Software … … 6618 5027 * 6619 5028 * ______________________________________________________________________ 6620 ** Title: RDBLUR6621 ** Authors: MMJT6622 ** Date: 4 Oct 19896623 ** Description: This function reads the type of instrumental6624 ** broadening to be simulated, and the associated parameters.6625 ** Legal types are 'NONE ', 'PSEUDOVOIGT ', 'GAUSSIAN '6626 ** and 'LORENTZIAN '. PSEUDOVOIGT requires 4 parameters,6627 ** GAUSSIAN requires a standard deviation in 2theta degrees, and6628 ** LORENTZIAN requires a halfwidth in 2theta degrees.6629 ** If the keyword 'TRIM' is at the end of the data line, intensity6630 ** near the origin will be ignored when the instrumental broadening6631 ** is added.6632 **6633 ** ARGUMENTS:6634 ** unit_no  Logical unit number that the data file is to6635 ** be read from. (input).6636 ** messge  A short message indicating what went wrong6637 ** (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_origin6643 **6644 ** RDBLUR returns logical .true. if a legal function for6645 ** simulating instrumental broadening was entered.6646 **6647 ** Legal types are:6648 ** 'NONE'6649 ** 'PSEUDOVOIGT' (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_no6659 * character*(*) messge6660 **6661 * character*200 line, dummy6662 * character*80 list(4)6663 * integer*4 CHOICE, LENGTH, CNTARG, i, iflag, arg_num, lin_len6664 **6665 ** external functions6666 * external CHOICE, LENGTH, CNTARG6667 ** external subroutines (Some compilers need them declared external)6668 ** external GETLNE, WRTLNE, TOUPPR6669 **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 angle6679 ** intensities for the broadened spectrum. If present, remove the6680 ** keyword 'TRIM', and set flag.6681 * i = index(line, ' TRIM')6682 * if(i.eq.0) then6683 * trim_origin = .false.6684 * else6685 * trim_origin = .true.6686 * write(line(i+1:i+4), '(a)') ' '6687 * endif6688 **6689 * list(1) = 'NONE '6690 * list(2) = 'GAUSSIAN '6691 * list(3) = 'LORENTZIAN '6692 * list(4) = 'PSEUDOVOIGT '6693 * iflag = CHOICE(line, list, 4)6694 **6695 * if(iflag.lt.1 .or. iflag.gt.4) then6696 * write(op,100)6697 *  'Illegal instrumental broadening type. Legal types are:'6698 * do 10 i = 1, 46699 * write(op,101) ' ', list(i)6700 * 10 continue6701 * messge = 'Instrumental broadening incorrectly specified.$'6702 * goto 9996703 * endif6704 **6705 * if(iflag.eq.1) blurring = NONE6706 **6707 * if(iflag.eq.2) then6708 **6709 ** check the number of arguments6710 * dummy = line(LENGTH('GAUSSIAN')+2:lin_len)6711 * arg_num = CNTARG(dummy)6712 * if(arg_num.ne.1 .and. arg_num.ne.3) then6713 * messge =6714 *  'Illegal number of GAUSSIAN arguments. FWHM, or u,v,w.$'6715 * goto 9996716 * endif6717 **6718 * if(arg_num.eq.1) then6719 * blurring = GAUSS6720 **6721 * messge = 'Problems reading Gaussian FWHM.$'6722 ** for strict FORTRAN77 compliance6723 * write(scrtch, *, err=990) dummy6724 * rewind(scrtch, err=990)6725 * read(scrtch, *, err=999) FWHM6726 * rewind(scrtch, err=990)6727 ** for more lenient, efficient, compilers6728 ** read(dummy, *, err=999) FWHM6729 **6730 * if(FWHM.lt.ZERO) then6731 * messge = 'Gaussian FWHM is ve.$'6732 * goto 9996733 * endif6734 **6735 * if(FWHM.lt.eps7) then6736 * write(op,102) 'Gaussian FWHM is zero.'6737 * write(op,102)6738 *  ' No instrumental broadening will be simulated.'6739 * blurring = NONE6740 * endif6741 **6742 * else if(arg_num.eq.3) then6743 * blurring = PV_GSS6744 **6745 * messge = 'Problems reading Gaussian u, v, and w factors.$'6746 ** for strict FORTRAN77 compliance6747 * write(scrtch, *, err=990) dummy6748 * rewind(scrtch, err=990)6749 * read(scrtch, *, err=999) pv_u, pv_v, pv_w6750 * rewind(scrtch, err=990)6751 ** for more lenient, efficient, compilers6752 ** read(dummy, *, err=999) pv_u, pv_v, pv_w6753 **6754 * pv_gamma = ZERO6755 **6756 * if(pv_w.lt.ZERO) then6757 * write(op,102)6758 *  'WARNING: Gaussian wfactor is less than zero.'6759 * write(op,102)6760 *  ' This may cause the broadening function to fail.'6761 * endif6762 **6763 * if(abs(pv_u).lt.eps7 .and. abs(pv_v).lt.eps7 .and.6764 *  abs(pv_w).lt.eps7) then6765 * write(op,102) 'Gaussian u, v, w factors are zero.'6766 * write(op,102) 'No instrumental broadening will be applied.'6767 * blurring = NONE6768 * endif6769 **6770 ** are the parameters equivalent to a Gaussian with6771 ** constant standard deviation?6772 * if(abs(pv_u).le.eps7 .and. abs(pv_v).le.eps7 .and.6773 *  pv_w.gt.ZERO) then6774 * FWHM = sqrt(pv_w)6775 * write(op,102)6776 *  'Gaussian factors imply a constant FWHM ', FWHM6777 * blurring = GAUSS6778 * endif6779 **6780 * endif6781 **6782 * endif6783 **6784 * if(iflag.eq.3) then6785 **6786 ** check the number of arguments6787 * dummy = line(LENGTH('LORENTZIAN')+2:lin_len)6788 * arg_num = CNTARG(dummy)6789 * if(arg_num.ne.1 .and. arg_num.ne.3) then6790 * messge =6791 *  'Illegal number of LORENTZIAN arguments. FWHM, or u,v,w.$'6792 * goto 9996793 * endif6794 **6795 * if(arg_num.eq.1) then6796 * blurring = LORENZ6797 **6798 * messge = 'Problems reading Lorentzian width.$'6799 ** for strict FORTRAN77 compliance6800 * write(scrtch, *, err=990) dummy6801 * rewind(scrtch, err=990)6802 * read(scrtch, *, err=999) FWHM6803 * rewind(scrtch, err=990)6804 ** for more lenient, efficient, compilers6805 ** read(dummy, *, err=999) FWHM6806 **6807 * if(FWHM.lt.ZERO) then6808 * messge =6809 *  'Negative value for Lorentzian FWHM. Must be positive.$'6810 * goto 9996811 * endif6812 **6813 * if(FWHM.lt.eps7) then6814 * write(op,102) 'Lorentzian FWHM is zero.'6815 * write(op,102)6816 *  ' No instrumental broadening will be simulated.'6817 * blurring = NONE6818 * endif6819 **6820 * else if(arg_num.eq.3) then6821 * blurring = PV_LRN6822 **6823 * messge = 'Problems reading Lorentzian u, v, and w factors.$'6824 ** for strict FORTRAN77 compliance6825 * write(scrtch, *, err=990) dummy6826 * rewind(scrtch, err=990)6827 * read(scrtch, *, err=999) pv_u, pv_v, pv_w6828 * rewind(scrtch, err=990)6829 ** for more lenient, efficient, compilers6830 ** read(dummy, *, err=999) pv_u, pv_v, pv_w6831 **6832 * pv_gamma = ONE6833 **6834 * if(pv_w.lt.ZERO) then6835 * write(op,102)6836 *  'WARNING: Lorentzian wfactor is less than zero.'6837 * write(op,102)6838 *  ' This may cause the broadening function to fail.'6839 * endif6840 **6841 * if(pv_u.lt.eps7 .and. pv_v.lt.eps7 .and.6842 *  abs(pv_w).lt.eps7) then6843 * write(op,102) 'Lorentzian u, v, w factors are zero.'6844 * write(op,102) 'No instrumental broadening will be applied.'6845 * blurring = NONE6846 * endif6847 **6848 ** are the parameters equivalent to a Gaussian with6849 ** constant standard deviation?6850 * if(pv_u.le.eps7 .and. pv_v.le.eps7 .and. pv_w.gt.ZERO) then6851 * FWHM = sqrt(pv_w)6852 * write(op,103)6853 *  'Lorentzian factors imply a constant FWHM ', FWHM6854 * blurring = LORENZ6855 * endif6856 **6857 * endif6858 **6859 * endif6860 **6861 * 500 continue6862 **6863 * if(iflag.eq.4) then6864 * blurring = PS_VGT6865 **6866 ** check the number of arguments6867 * dummy = line(LENGTH('PSEUDOVOIGT')+2:len(line))6868 * arg_num = CNTARG(dummy)6869 * if(arg_num.ne.4) then6870 * messge = 'Exactly 4 pseudoVoigt arguments are required.$'6871 * goto 9996872 * endif6873 **6874 * messge =6875 *  'Problems reading pseudoVoigt u, v, w, gamma factors.$'6876 **6877 ** for strict FORTRAN77 compliance6878 * write(scrtch, *, err=990) dummy6879 * rewind(scrtch, err=990)6880 * read(scrtch, *, err=999) pv_u, pv_v, pv_w, pv_gamma6881 * rewind(scrtch, err=990)6882 ** for more lenient, efficient, compilers6883 ** read(dummy, *, err=999) pv_u, pv_v, pv_w, pv_gamma6884 **6885 * if(pv_w.lt.ZERO) then6886 * write(op,102)6887 *  'WARNING: pseudoVoigt w factor is less than zero.'6888 * write(op,102)6889 *  ' This may cause the pseudoVoigt function to fail.'6890 * endif6891 **6892 * if(abs(pv_u).lt.eps7 .and. abs(pv_v).lt.eps7 .and.6893 *  abs(pv_w).lt.eps7) then6894 * write(op,102) 'PseudoVoigt u, v, w factors are zero.'6895 * write(op,102) 'No instrumental broadening will be simulated.'6896 * blurring = NONE6897 * endif6898 **6899 * if(pv_gamma.lt.ZERO) then6900 * messge =6901 *  'PseudoVoigt gamma factor is ve (must be between 0 and 1).$'6902 * goto 9996903 * endif6904 **6905 * if(pv_gamma.gt.ONE) then6906 * messge = 'PseudoVoigt gamma factor is greater than 1.$'6907 * goto 9996908 * endif6909 **6910 ** are the pseudoVoigt parameters equivalent to a Gaussian with6911 ** 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) then6914 * FWHM = sqrt(pv_w)6915 * write(op,102)6916 *  'PseudoVoigt factors are equivalent to a Gaussian'6917 * write(op,103) ' with FWHM ', FWHM6918 * blurring = GAUSS6919 * endif6920 **6921 ** are the pseudoVoigt parameters equivalent to a Lorentzian with6922 ** constant halfwidth?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) then6925 * FWHM = sqrt(pv_w)6926 * write(op,102)6927 *  'PseudoVoigt factors are equivalent to a Lorentzian'6928 * write(op,103) ' with FWHM ', FWHM6929 * blurring = LORENZ6930 * endif6931 **6932 * endif6933 **6934 * RDBLUR = .true.6935 * return6936 * 990 messge = 'Problems using scratch file in RDBLUR.$'6937 * 999 call WRTLNE(line)6938 * return6939 * 100 format(1x, 'ERROR: ', a)6940 * 101 format(1x, 2a)6941 * 102 format(1x, a)6942 * 103 format(1x, a, g12.5)6943 * end6944 **6945 ** ______________________________________________________________________6946 ** Title: RDCELL6947 ** Authors: MMJT6948 ** Date: 4 Oct 19896949 ** Description: This function reads the layer mesh dimensions from6950 ** the second data line of the data file. It checks that the header6951 ** 'STRUCTURAL' is present.6952 **6953 ** ARGUMENTS:6954 ** unit_no  Logical unit number that the data file is to6955 ** be read from. (input).6956 ** messge  A short message indicating what went wrong6957 ** (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_gamma6962 **6963 ** RDCELL returns logical .true. if acceptable values6964 ** 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_no6972 * character*(*) messge6973 **6974 * character*200 line6975 **6976 ** external subroutines (Some compilers need them declared external)6977 ** external GETLNE, TOUPPR, WRTLNE6978 **6979 * RDCELL = .false.6980 **6981 ** first check that 'STRUCTURAL' header is present6982 * 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 9996987 **6988 * messge = 'a, b, c, and gamma are improperly specified.$'6989 * call GETLNE(unit_no, line, *999)6990 **6991 ** for strict FORTRAN77 compliance6992 * write(scrtch, *, err=990) line6993 * rewind(scrtch, err=990)6994 * read(scrtch, *, err=999) cell_a, cell_b, cell_c, cell_gamma6995 * rewind(scrtch, err=990)6996 ** for more lenient, efficient, compilers6997 ** read(line, *, err=999) cell_a, cell_b, cell_c, cell_gamma6998 **6999 * if(cell_a.le.ZERO) then7000 * messge = 'Negative (or zero) value for cell_a.$'7001 * goto 9997002 * endif7003 **7004 * if(cell_b.le.ZERO) then7005 * messge = 'Negative (or zero) value for cell_b.$'7006 * goto 9997007 * endif7008 **7009 * if(cell_c.le.ZERO) then7010 * messge = 'Negative (or zero) value for cell_c.$'7011 * goto 9997012 * endif7013 **7014 * if(cell_gamma.le.ZERO) then7015 * messge = 'Negative (or zero) value for cell_gamma.$'7016 * goto 9997017 * endif7018 **7019 * if(cell_gamma.ge.ONE_EIGHTY) then7020 * messge = 'Cell_gamma is greater than 180 degrees.$'7021 * goto 9997022 * endif7023 **7024 * cell_gamma = cell_gamma * DEG2RAD7025 * RDCELL = .true.7026 * return7027 * 990 messge = 'Problems using scratch file in RDCELL.$'7028 * 999 call WRTLNE(line)7029 * return7030 * end7031 **7032 ** ______________________________________________________________________7033 ** Title: RDFILE7034 ** Author: MMJT7035 ** Date: 3 Oct 1989; 3 Mar 19957036 ** Description: This function reads in the data from the userdefined7037 ** data file. The data file is in NaurBackus form, and the user may7038 ** insert '{' and '}' delimiters to comment the data. The data file7039 ** is not read in its raw form. The subroutine GETLNE is used7040 ** to read each line, stripped of its comments, and leftjustified.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 ** save7051 **7052 * character*(*) fname7053 **7054 * character*200 messge7055 * logical goodfile7056 * logical RDRADN, RDWAVL, RDBLUR, RDCELL, RDPTGP7057 * logical RDNLAY, RDWDTH, RDLAYR, RDSTAK, RDTRNS7058 * integer*4 LENGTH, unit_no7059 **7060 ** external functions7061 * external RDRADN, RDWAVL, RDBLUR, RDCELL, RDPTGP7062 * external RDNLAY, RDWDTH, RDLAYR, RDSTAK, RDTRNS7063 * external LENGTH7064 **7065 * RDFILE = .false.7066 * unit_no = df7067 * 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) then7082 * write(op,100) '''', fname(1:LENGTH(fname)),''' read in.'7083 * else7084 * write(op,200) messge(1:index(messge,'$')  1)7085 * endif7086 **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 * return7093 **7094 * 900 write(op,200) messge(1:index(messge,'$')  1)7095 * return7096 * 100 format(1x, 3a)7097 * 101 format(3a)7098 * 200 format(1x, 'ERROR: ', a)7099 * end7100 **7101 ** ______________________________________________________________________7102 ** Title: RDLAYR7103 ** Authors: MMJT7104 ** Date: 7 June 19907105 ** Description: This function reads in the layer atomic coordinates.7106 ** Each set of data must be prefaced with a layer number, and a7107 ** statement of whether or not the layer is centrosymmetric. If the7108 ** layer is centrosymmetric, only the asymmetric coordinates need7109 ** be entered (with appropriate fractional occupancies for special7110 ** positions.7111 ** If a layer has the same structure as a previously defined layer7112 ** the layer can be defined as 'LAYER i = j', where j is a previously7113 ** defined layer.7114 ** The format assumed for atom coordinates is:7115 **7116 ** 'atom name', number, x, y, z, DebyeWaller factor, occupancy7117 **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 to7123 ** be read from. (input).7124 ** messge  A short message indicating what went wrong7125 ** (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, NONE7130 **7131 ** modifies: l_symmetry, a_name, a_number, a_pos, a_B, a_occup,7132 ** l_n_atoms, n_actual7133 **7134 ** RDLAYR returns logical .true. if all the layer data7135 ** 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_no7143 * character*(*) messge7144 **7145 * character*200 line, tmpline, arg7146 * character*80 list(2)7147 * integer*4 CHOICE, LENGTH, NXTARG, m, n, i, j, i2, iflag7148 * real*8 RDNMBR, tmp7149 * logical*4 ok7150 **7151 ** external functions7152 * external CHOICE, LENGTH, NXTARG, RDNMBR7153 ** external subroutines (Some compilers need them declared external)7154 ** external GETLNE, TOUPPR, WRTLNE7155 **7156 * RDLAYR = .false.7157 **7158 * do 5 i = 1, n_layers7159 * high_atom(i) = ZERO7160 * low_atom(i) = ZERO7161 * 5 continue7162 **7163 * m = 07164 * messge = 'Unexpected EOF before LAYER 1.$'7165 * call GETLNE(unit_no, line, *999)7166 * call TOUPPR(line)7167 **7168 * do 30 i = 1, n_layers7169 * write(messge,100) 'LAYER header not seen for LAYER ',i,'.$'7170 * if(index(line, 'LAYER').ne.1) goto 9997171 ** read layer number7172 * 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 compliance7176 * write(scrtch, *, err=990) tmpline7177 * rewind(scrtch, err=990)7178 * read(scrtch, *, err=999) j7179 * rewind(scrtch, err=990)7180 ** for more lenient, efficient, compilers7181 ** read(tmpline, *, err=999) j7182 **7183 ** layers must be entered in sequence7184 * if(j.ne.i) goto 9997185 **7186 * j = index(line, '=')7187 ** see if this layer is equivalenced to another7188 * if(j.ne.0) then7189 **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 compliance7194 * write(scrtch, *, err=990) tmpline7195 * rewind(scrtch, err=990)7196 * read(scrtch, *, err=999) i27197 * rewind(scrtch, err=990)7198 ** for more lenient, efficient, compilers7199 ** read(tmpline, *, err=999) i27200 **7201 ** layer i2 must be defined already7202 * if(i2.ge.i .or. i2.lt.1) goto 9997203 * 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 * else7210 * m = m + 17211 * l_actual(i) = m7212 **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 9997225 **7226 * if(iflag.eq.1) then7227 * l_symmetry(m) = NONE7228 * else7229 * l_symmetry(m) = CENTRO7230 * endif7231 **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 9997235 **7236 * j = 07237 * 20 j = j + 17238 * if(j.gt.MAX_A) then7239 * write(op,101) 'ERROR: Too many atoms on LAYER ',i,'.'7240 * write(op,102) ' Maximum allowed = ', MAX_A7241 * messge = ' $'7242 * goto 9997243 * endif7244 **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 identifier7251 * n = NXTARG(tmpline,arg)7252 * if(n.le.0) goto 9997253 * read(arg,'(i3)') a_number(j,m)7254 **7255 ** read x_relative position7256 * n = NXTARG(tmpline,arg)7257 * if(n.le.0) goto 9997258 * a_pos(1,j,m) = RDNMBR(arg,ok)7259 * if(.not.ok) goto 9997260 **7261 ** read y_relative position7262 * n = NXTARG(tmpline,arg)7263 * if(n.le.0) goto 9997264 * a_pos(2,j,m) = RDNMBR(arg,ok)7265 * if(.not.ok) goto 9997266 **7267 ** read z_relative position7268 * n = NXTARG(tmpline,arg)7269 * if(n.le.0) goto 9997270 * a_pos(3,j,m) = RDNMBR(arg,ok)7271 * if(.not.ok) goto 9997272 **7273 ** read DebyeWaller factor7274 * n = NXTARG(tmpline,arg)7275 * if(n.le.0) goto 9997276 * a_B(j,m) = RDNMBR(arg,ok)7277 * if(.not.ok) goto 9997278 **7279 ** read occupancy7280 * n = NXTARG(tmpline,arg)7281 * if(n.le.0) goto 9997282 * a_occup(j,m) = RDNMBR(arg,ok)7283 * if(.not.ok) goto 9997284 **7285 ** Check values7286 * if(a_B(j,m).lt.ZERO) then7287 * messge = 'Negative DebyeWaller factor.$'7288 * goto 9997289 * endif7290 **7291 * if(a_occup(j,m).lt.ZERO) then7292 * messge = 'Negative atom occupancy.$'7293 * goto 9997294 * endif7295 **7296 * if(a_occup(j,m).gt.ONE) then7297 * messge = 'Occupancy greater than 1.$'7298 * goto 9997299 * endif7300 **7301 ** get extreme upper and lower atom positions for errorchecking later on7302 * tmp = a_pos(3,j,m)7303 * if(tmp.gt.high_atom(m)) high_atom(m) = tmp7304 * if(tmp.lt.low_atom(m)) low_atom(m) = tmp7305 **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 in7309 ** but first we must convert to uppercase for the test. If there are7310 ** more atoms, then we must convert 'line' to its original format.7311 * tmpline = line7312 * call TOUPPR(line)7313 ** check data is in correct sequence7314 * messge = 'STACKING data is missing.$'7315 * if(index(line,'PARAMETERS').eq.1) goto 9997316 * if(index(line,'LAYER').ne.17317 *  .and. index(line,'STACKING').ne.1) then7318 * line = tmpline7319 * goto 207320 * endif7321 * l_n_atoms(m) = j7322 * endif7323 * 30 continue7324 * n_actual = m7325 **7326 ** make sure we have genuine lowest and highest atoms in each layer.7327 * do 40 i = 1, n_actual7328 * if(l_symmetry(i).eq.CENTRO) then7329 * if(low_atom(i).gt.high_atom(i)) then7330 * high_atom(i) = low_atom(i)7331 * else if(low_atom(i).lt.high_atom(i)) then7332 * low_atom(i) = high_atom(i)7333 * endif7334 * endif7335 * 40 continue7336 **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 return7343 * 990 messge = 'Problems using scratch file in RDLAYR.$'7344 * 999 call WRTLNE(line)7345 * return7346 * 100 format(a, i2, a)7347 * 101 format(1x, a, i2, a)7348 * 102 format(1x, a, i4)7349 * 103 format(3a)7350 * end7351 **7352 ** ______________________________________________________________________7353 ** Title: RDNLAY7354 ** Authors: MMJT7355 ** Date: 4 Oct 19897356 ** Description: This function reads the first line of the data file7357 ** to extract the number of layers in the structure.7358 **7359 ** ARGUMENTS:7360 ** unit_no  Logical unit number that the data file is to7361 ** be read from. (input).7362 ** messge  A short message indicating what went wrong7363 ** (if anything) during the datafile read.7364 ** messge is terminated by a token '$'. (output).7365 **7366 ** COMMON VARIABLES:7367 ** modifies: n_layers7368 **7369 ** RDNLAY returns logical .true. if an acceptable value7370 ** 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_no7378 * character*(*) messge7379 **7380 * character*200 line7381 **7382 ** external subroutines (Some compilers need them declared external)7383 ** external GETLNE, WRTLNE7384 **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 compliance7392 * write(scrtch, *, err=990) line7393 * rewind(scrtch, err=990)7394 * read(scrtch, *, err=999) n_layers7395 * rewind(scrtch, err=990)7396 ** for more lenient, efficient, compilers7397 ** read(line, *, err=999) n_layers7398 **7399 * if(n_layers.eq.0) then7400 * messge = 'Zero number of layers entered.$'7401 * goto 9997402 * else if(n_layers.lt.0) then7403 * messge = 'Negative number of layers entered.$'7404 * goto 9997405 * else if(n_layers.gt.MAX_L) then7406 * write(messge,100)7407 *  'Too many layers: number should not exceed ',MAX_L,'.$'7408 * goto 9997409 * endif7410 **7411 * RDNLAY = .true.7412 * return7413 * 990 messge = 'Problems using scratch file in RDNLAY.$'7414 * 999 call WRTLNE(line)7415 * return7416 * 100 format(a, i2, a)7417 * end7418 **7419 ** ______________________________________________________________________7420 ** Title: RDNMBR7421 ** Author: MMJT7422 ** Date: 21 January 20057423 ** Description: Reads the character string 'numberstring' to extract a7424 ** floating point number. The character string can be either a true7425 ** 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 number7432 ** ______________________________________________________________________7433 **7434 * real*8 function RDNMBR(numberstring,ok)7435 * include 'DIFFaX.par'7436 * save7437 **7438 * character*(*) numberstring7439 * logical ok7440 **7441 * character*200 arg7442 * integer*4 numerator, denominator, i, j, lin_len7443 * integer*4 NXTARG7444 **7445 ** external functions7446 * external NXTARG7447 **7448 * ok = .true.7449 * RDNMBR = ZERO7450 **7451 * lin_len = len(numberstring)7452 * j = index(numberstring, '/')7453 * if(j.eq.0) then7454 ** It was already a number7455 * read(numberstring, *) RDNMBR7456 * else if(j.lt.0 .or. j.ge.lin_len) then7457 ** An error happened7458 * goto 9997459 * else7460 ** 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 9997465 * read(arg, *) numerator7466 * i = NXTARG(numberstring,arg)7467 * if(i.le.0) goto 9997468 * read(arg, *) denominator7469 * if(denominator.le.ZERO) goto 9997470 * RDNMBR = dble(numerator) / dble(denominator)7471 * endif7472 **7473 * 100 return7474 * 999 ok = .false.7475 * return7476 * end7477 **7478 ** ______________________________________________________________________7479 ** Title: RDPTGP7480 ** Author: MMJT7481 ** Date: 4 Oct 19897482 ** Description: This function reads the user's estimate of the7483 ** point group symmetry of the diffraction data.7484 **7485 ** ARGUMENTS:7486 ** unit_no  Logical unit number that the data file is to7487 ** be read from. (input).7488 ** messge  A short message indicating what went wrong7489 ** (if anything) during the datafile read.7490 ** messge is terminated by a token '$'. (output).7491 **7492 ** COMMON VARIABLES:7493 ** modifies: pnt_grp, SymGrpNo, tolerance7494 **7495 ** RDPTGP returns logical .true. if a legal point group7496 ** 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_no7519 * character*(*) messge7520 **7521 * character*200 line, oldline7522 * character*80 list(12)7523 * integer*4 CHOICE, LENGTH, PRUNE, i, iflag7524 **7525 ** external functions7526 * external CHOICE, LENGTH, PRUNE7527 ** external subroutines (Some compilers need them declared external)7528 ** external GETLNE, TOUPPR, WRTLNE7529 **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) then7550 * write(op,100) 'Unrecognized symmetry type. Legal types are:'7551 * do 10 i = 1, 127552 * write(op,101) ' ',list(i)7553 * 10 continue7554 * goto 9997555 * endif7556 **7557 * pnt_grp = list(iflag)7558 **7559 ** set default tolerance on intensity variability7560 * tolerance = ONE7561 **7562 * if(iflag.ge.1 .and. iflag.le.11) then7563 * SymGrpNo = iflag7564 * else if(iflag.eq.12) then7565 ** a value of UNKNOWN = 1 alerts OPTIMZ that user entered 'UNKNOWN'7566 * SymGrpNo = UNKNOWN7567 ** Strip off the symmetry point group text string.7568 * oldline = line7569 * 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) then7573 * messge = 'Problems reading symmetry tolerance parameter.$'7574 * goto 9997575 * else if(i.gt.1) then7576 ** There was something on the line, try to read tolerance7577 * messge = 'Problems reading symmetry tolerance parameter.$'7578 ** for strict FORTRAN77 compliance7579 * write(scrtch, *, err=990) line7580 * rewind(scrtch, err=990)7581 * read(scrtch, *, err=999) tolerance7582 * rewind(scrtch, err=990)7583 ** for more lenient, efficient, compilers7584 ** read(line, *, err=999) tolerance7585 **7586 * if(tolerance.le.ZERO) then7587 * messge = 'Negative (or zero) value for tolerance.$'7588 * goto 9997589 * endif7590 **7591 * if(tolerance.lt.eps2) then7592 * tolerance = eps27593 * write(op,102)7594 * write(op,103) tolerance7595 * endif7596 **7597 * endif7598 * endif7599 **7600 ** convert from a percentage7601 * tolerance = tolerance * eps27602 **7603 * RDPTGP = .true.7604 * return7605 * 990 messge = 'Problems using scratch file in RDPTGP.$'7606 * 999 call WRTLNE(line)7607 * return7608 * 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 * end7613 **7614 ** ______________________________________________________________________7615 ** Title: RDRADN7616 ** Authors: MMJT7617 ** Date: 15 Feb 19907618 ** Description: This function reads the radiation type. The choices are7619 ** XRAY, 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 to7624 ** be read from. (input).7625 ** messge  A short message indicating what went wrong7626 ** (if anything) during the datafile read.7627 ** messge is terminated by a token '$'. (output).7628 **7629 ** COMMON VARIABLES:7630 ** uses: ELECTN, NEUTRN, X_RAY7631 **7632 ** modifies: rad_type7633 **7634 ** RDRADN returns logical .true. if an acceptable7635 ** 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_no7643 * character*(*) messge7644 **7645 * character*200 line7646 * character*80 list(3)7647 * integer*4 CHOICE, iflag7648 **7649 ** external function7650 * external CHOICE7651 ** external subroutines (Some compilers need them declared external)7652 ** external GETLNE, TOUPPR, WRTLNE7653 **7654 * RDRADN = .false.7655 **7656 ** first check that 'INSTRUMENTAL' header is present7657 * 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 9997663 **7664 * messge = 'radiation type not specified.$'7665 * call GETLNE(unit_no, line, *999)7666 * call TOUPPR(line)7667 * list(1) = 'XRAY '7668 * list(2) = 'NEUTRON '7669 * list(3) = 'ELECTRON '7670 * iflag = CHOICE(line, list, 3)7671 **7672 * if(iflag.eq.1) then7673 * rad_type = X_RAY7674 * else if(iflag.eq.2) then7675 * rad_type = NEUTRN7676 * else if(iflag.eq.3) then7677 * rad_type = ELECTN7678 * else7679 * messge ='Invalid radiation type (XRAY, NEUTRON or ELECTRON).$'7680 * goto 9997681 * endif7682 **7683 * RDRADN = .true.7684 * return7685 * 999 call WRTLNE(line)7686 * return7687 * end7688 **7689 ** ______________________________________________________________________7690 ** Title: RDRNGE7691 ** Authors: MMJT7692 ** Date: 15 Feb 19907693 ** Description: This function reads the angular range for the spectrum.7694 **7695 ** ARGUMENTS:7696 ** unit_no  Logical unit number that the data file is to7697 ** be read from. (input).7698 **7699 ** COMMON VARIABLES:7700 ** uses: DEG2RAD7701 **7702 ** modifies: th2_min, th2_max, d_theta7703 **7704 ** RDRNGE returns logical .true. if an acceptable angular7705 ** 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, WRTLNE7714 **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_theta7720 * if(CFile) write(op,105) th2_min, th2_max, d_theta7721 **7722 * if(th2_min.lt.ZERO) then7723 * write(op,110) '2theta min is negative.'7724 * goto 9997725 * endif7726 **7727 * if(th2_min.ge.ONE_EIGHTY) then7728 * write(op,110) '2theta min exceeds 180 degrees.'7729 * goto 9997730 * endif7731 **7732 * if(th2_max.le.ZERO) then7733 * write(op,110) 'Negative (or zero) value for 2theta min.'7734 * goto 9997735 * endif7736 **7737 * if(th2_max.gt.ONE_EIGHTY) then7738 * write(op,110) '2theta max exceeds 180 degrees.'7739 * goto 9997740 * endif7741 **7742 ** if th2_max = 180, reduce it slightly to keep us out of trouble7743 * if(th2_max.eq.ONE_EIGHTY) th2_max = th2_max  eps47744 **7745 * if(d_theta.le.ZERO) then7746 * write(op,110) 'Negative (or zero) value for 2theta increment.'7747 * goto 9997748 * endif7749 **7750 * if(th2_min.ge.th2_max) then7751 * write(op,110) '2theta min is larger than 2theta max.'7752 * goto 9997753 * endif7754 **7755 * if((th2_max  th2_min).lt.d_theta) then7756 * write(op,110)7757 *  '2theta increment is larger than 2theta max  2theta min.'7758 * goto 9997759 * endif7760 **7761 * th2_min = th2_min * DEG2RAD7762 * th2_max = th2_max * DEG2RAD7763 * d_theta = HALF * DEG2RAD * d_theta7764 **7765 * RDRNGE = .true.7766 * return7767 **7768 * 999 return7769 * 100 format(1x, a)7770 * 105 format(1x, g11.4, 2(3x, g11.4))7771 * 110 format(1x, 'ERROR: ', a)7772 * end7773 **7774 ** ______________________________________________________________________7775 ** Title: RDSTAK7776 ** Author: MMJT7777 ** Date: 4 Oct 19897778 ** Description: This function reads the STACKING input of the datafile.7779 ** The stacking is either 'EXPLICIT', whereupon DIFFaX expects the7780 ** user to enter an explicit sequence of layers, or 'RECURSIVE'. If7781 ** 'EXPLICIT' was specified, the user may state that the sequence is7782 ** 'RANDOM', whereby DIFFaX will generate an explicit sequence of7783 ** layers weighted by the stacking probabilities specified later on7784 ** under 'PARAMETERS'. The user can specify a crystal thickness (in7785 ** terms of the number of layers) under the 'RECURSIVE' option.7786 **7787 ** ARGUMENTS:7788 ** unit_no  Logical unit number that the data file is to7789 ** be read from. (input).7790 ** messge  A short message indicating what went wrong7791 ** (if anything) during the datafile read.7792 ** messge is terminated by a token '$'. (output).7793 **7794 ** COMMON VARIABLES:7795 ** uses: n_layers7796 **7797 ** modifies: recrsv, xplcit, rndm, l_cnt, l_seq, inf_thick7798 **7799 ** RDSTAK returns logical .true. if the stacking7800 ** 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_no7808 * character*(*) messge7809 **7810 * character*200 line, tmpline7811 * character*80 list(2)7812 * integer*4 CHOICE, LAYCNT, LENGTH, i, j, iflag7813 **7814 ** external functions7815 * external CHOICE, LAYCNT, LENGTH7816 ** external subroutines (Some compilers need them declared external)7817 ** external GETLNE, TOUPPR, WRTLNE7818 **7819 * RDSTAK = .false.7820 **7821 * call GETLNE(unit_no, line, *999)7822 * call TOUPPR(line)7823 * if(index(line,'STACKING ').ne.1) then7824 * write(messge,103) 'The ', n_layers,7825 *  ' layers are read in. ''STACKING'' section not found.$'7826 * goto 9997827 * endif7828 **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) then7838 * write(op,100) 'Illegal STACKING keyword. Legal types are:'7839 * do 10 i = 1, 27840 * write(op,101) ' ',list(i)7841 * 10 continue7842 * messge = ' $'7843 * goto 9997844 * endif7845 **7846 ** set flags to indicate stacking type7847 * recrsv = .false.7848 * xplcit = .false.7849 * rndm = .false.7850 * inf_thick = .false.7851 ** initialize l_cnt7852 * l_cnt = 07853 * if(iflag.eq.1) then7854 * xplcit = .true.7855 * else if(iflag.eq.2) then7856 * recrsv = .true.7857 * endif7858 **7859 * if(recrsv) then7860 * 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) then7865 * inf_thick = .false.7866 * messge =7867 *  'Problems reading the number of RECURSIVE layers.$'7868 ** for strict FORTRAN77 compliance7869 * write(scrtch, *, err=990) line7870 * rewind(scrtch, err=990)7871 * read(scrtch, *, err=999) l_cnt7872 * rewind(scrtch, err=990)7873 ** for more lenient, efficient, compilers7874 ** read(line, *, err=999) l_cnt7875 ** check l_cnt7876 * if(l_cnt.eq.0) then7877 * messge = 'The number of RECURSIVE layers cannot be zero.$'7878 * else if(l_cnt.lt.0) then7879 * messge =7880 *  'The number of RECURSIVE layers cannot be negative.$'7881 * else if(l_cnt.gt.RCSV_MAX) then7882 * write(op,104)7883 *  'WARNING: Number of RECURSIVE layers ', l_cnt,7884 *  ' exceeds the maximum ', RCSV_MAX7885 * write(op,105) 'An INFINITE thickness is assumed.'7886 ** reset l_cnt and proceed assuming infinite thickness7887 * l_cnt = 07888 * inf_thick = .true.7889 * endif7890 * else7891 * inf_thick = .true.7892 * endif7893 **7894 * else if(xplcit) then7895 ** read in layer sequence7896 * messge = 'Problem reading EXPLICIT layer sequencing.$'7897 * call GETLNE(unit_no, line, *999)7898 * call TOUPPR(line)7899 * rndm = index(line,'RANDOM ').ne.07900 ** catch possible error  user requested INFINITE7901 * if(.not.rndm) then7902 * if(index(line,'INFINITE').ne.0) then7903 * write(op,107) 'Maximum number of layers allowed is ',XP_MAX7904 * messge = 'Too many layers for an EXPLICIT sequence.$'7905 * goto 9997906 * endif7907 * endif7908 **7909 * if(rndm) then7910 ** 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 compliance7917 * write(scrtch, *, err=990) tmpline7918 * rewind(scrtch, err=990)7919 * read(scrtch, *, err=999) l_cnt7920 * rewind(scrtch, err=990)7921 ** for more lenient, efficient, compilers7922 ** read(tmpline, *, err=999) l_cnt7923 **7924 * else7925 * l_cnt = 17926 * messge =7927 *  'Problems reading EXPLICIT layer sequence in RDSTAK.$'7928 * 20 i = LAYCNT(line)7929 **7930 * if(i.ge.0 .and. i.le.40) then7931 **7932 ** for strict FORTRAN77 compliance7933 * write(scrtch, *, err=990) line7934 * rewind(scrtch, err=990)7935 * read(scrtch, *, err=999) (l_seq(j), j=l_cnt,l_cnt+i1)7936 * rewind(scrtch, err=990)7937 ** for more lenient, efficient, compilers7938 ** read(line, *, err=999) (l_seq(j), j=l_cnt,l_cnt+i1)7939 **7940 * l_cnt = l_cnt + i7941 * else7942 * goto 9997943 * endif7944 **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 207949 * l_cnt = l_cnt  17950 **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 * endif7957 **7958 ** check we didn't enter too many layers7959 * if(l_cnt.gt.XP_MAX) then7960 * write(op,104) 'WARNING: No. of EXPLICIT layers, ',l_cnt,7961 *  ' exceeds maximum of ',XP_MAX7962 * write(op,106) XP_MAX, ' EXPLICIT layers assumed'7963 * l_cnt = XP_MAX7964 * endif7965 **7966 ** check that all the layer numbers are within bounds7967 * if(.not.rndm) then7968 * do 30 i = 1, l_cnt7969 * if(l_seq(i).gt.n_layers) then7970 * write(messge,102) 'Illegal layer number,', l_seq(i),7971 *  '. Number must not exceed', n_layers, '.$'7972 * goto 9807973 * else if(l_seq(i).lt.1) then7974 * write(messge,103) 'Illegal layer number,', l_seq(i),7975 *  '. Number must be greater than 0.$'7976 * goto 9807977 * endif7978 * 30 continue7979 * endif7980 * endif7981 **7982 * RDSTAK = .true.7983 * 980 return7984 * 990 messge = 'Problems using scratch file in RDSTAK.$'7985 * 999 call WRTLNE(line)7986 * return7987 * 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 * end7996 **7997 ** ______________________________________________________________________7998 ** Title: RDTRNS7999 ** Authors: MMJT8000 ** Date: 4 Oct 1989; 15 Mar 1995; 29th Aug 20008001 ** Description: This function reads the layer stacking probabilities8002 ** and stacking vectors.8003 **8004 ** ARGUMENTS:8005 ** unit_no  Logical unit number that the data file is to8006 ** be read from. (input).8007 ** messge  A short message indicating what went wrong8008 ** (if anything) during the datafile read.8009 ** messge is terminated by a token '$'. (output).8010 **8011 ** COMMON VARIABLES:8012 ** uses: n_layers, recrsv, rndm, xplcit8013 **8014 ** modifies: l_alpha, l_r, r_B11, r_B22, r_B33, r_B12, r_B23,8015 ** r_B318016 **8017 ** RDTRNS returns logical .true. if all the stacking8018 ** 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_no8026 * character*(*) messge8027 **8028 * character*200 line, tmpline, arg8029 * integer*4 PARENT, NXTARG, i, j, n, tempfact8030 * real*8 sum, RDNMBR8031 * logical*4 ok8032 **8033 ** external function8034 * external PARENT, NXTARG, RDNMBR8035 ** external subroutines (Some compilers need them declared external)8036 ** external GETLNE, TOUPPR, WRTLNE8037 **8038 * RDTRNS = .false.8039 **8040 * call GETLNE(unit_no, line, *999)8041 * call TOUPPR(line)8042 * if(index(line,'TRANSITIONS ').ne.1) then8043 * messge = '''TRANSITIONS'' section not found.$'8044 * goto 9998045 * endif8046 **8047 * do 10 i = 1, n_layers8048 * do 20 j = 1, n_layers8049 * 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 there8053 ** are FatsWaller (DebyeWaller for layers!) factors8054 * write(messge,100)8055 *  'Confusing use of parentheses in data for vector ',i,j,'.$'8056 * tempfact = PARENT(line)8057 * if(tempfact.lt.0) goto 9998058 ** initialize8059 * l_r(1,j,i) = ZERO8060 * l_r(2,j,i) = ZERO8061 * l_r(3,j,i) = ZERO8062 *8063 * r_B11(j,i) = ZERO8064 * r_B22(j,i) = ZERO8065 * r_B33(j,i) = ZERO8066 * r_B12(j,i) = ZERO8067 * r_B31(j,i) = ZERO8068 * r_B23(j,i) = ZERO8069 **8070 * tmpline = line8071 **8072 * if(tempfact.eq.1) then8073 * write(messge,100)8074 *  'Expected alpha_ij, R_ij, (dR_ij) for vector ',i,j,'.$'8075 * else8076 * write(messge,100)8077 *  'Expected alpha_ij, R_ij data for vector ',i,j,'.$'8078 * endif8079 **8080 ** read stacking probability8081 * n = NXTARG(tmpline,arg)8082 * if(n.le.0) goto 9998083 * l_alpha(j,i) = RDNMBR(arg,ok)8084 * if(.not.ok) goto 9998085 ** If zero, there is no need to read further8086 * if(dabs(l_alpha(j,i)).lt.eps6) goto 208087 **8088 ** read stacking xvector8089 * n = NXTARG(tmpline,arg)8090 * if(n.le.0) goto 9998091 * l_r(1,j,i) = RDNMBR(arg,ok)8092 * if(.not.ok) goto 9998093 **8094 ** read stacking yvector8095 * n = NXTARG(tmpline,arg)8096 * if(n.le.0) goto 9998097 * l_r(2,j,i) = RDNMBR(arg,ok)8098 * if(.not.ok) goto 9998099 **8100 ** read stacking zvector8101 * n = NXTARG(tmpline,arg)8102 * if(n.le.0) goto 9998103 * l_r(3,j,i) = RDNMBR(arg,ok)8104 * if(.not.ok) goto 9998105 **8106 ** Read temperature factors, if present8107 * if(tempfact.eq.1) then8108 **8109 ** read B11 "FatsWaller"8110 * n = NXTARG(tmpline,arg)8111 * if(n.le.0) goto 9998112 * r_B11(j,i) = RDNMBR(arg,ok)8113 * if(.not.ok) goto 9998114 **8115 ** read B22 "FatsWaller"8116 * n = NXTARG(tmpline,arg)8117 * if(n.le.0) goto 9998118 * r_B22(j,i) = RDNMBR(arg,ok)8119 * if(.not.ok) goto 9998120 **8121 ** read B33 "FatsWaller"8122 * n = NXTARG(tmpline,arg)8123 * if(n.le.0) goto 9998124 * r_B33(j,i) = RDNMBR(arg,ok)8125 * if(.not.ok) goto 9998126 **8127 ** read B12 "FatsWaller"8128 * n = NXTARG(tmpline,arg)8129 * if(n.le.0) goto 9998130 * r_B12(j,i) = RDNMBR(arg,ok)8131 * if(.not.ok) goto 9998132 **8133 ** read B31 "FatsWaller"8134 * n = NXTARG(tmpline,arg)8135 * if(n.le.0) goto 9998136 * r_B31(j,i) = RDNMBR(arg,ok)8137 * if(.not.ok) goto 9998138 **8139 ** read B23 "FatsWaller"8140 * n = NXTARG(tmpline,arg)8141 * if(n.le.0) goto 9998142 * r_B23(j,i) = RDNMBR(arg,ok)8143 * if(.not.ok) goto 9998144 **8145 * endif8146 **8147 ** stacking probabilities should not be negative.8148 * if(l_alpha(j,i).lt.ZERO) then8149 * write(messge,100) 'alpha',i,j,' is negative.$'8150 * goto 9998151 * endif8152 ** the Rz(i,i) vectors should not be zero or negative.8153 ** if an ii transition exists8154 * if(j.eq.i) then8155 * if(l_alpha(i,i).ne.ZERO) then8156 * if(l_r(3,i,i).le.ZERO) then8157 * write(messge,100)8158 * 'Invalid negative (or zero) zcomponent for Rz',i,i,'.$'8159 * goto 9998160 * endif8161 * endif8162 * endif8163 * 20 continue8164 * 10 continue8165 **8166 ** Check stacking probabilities.8167 ** Generate error message if they do not sum to 1, and we are using8168 ** a RECURSIVE stacking description. If EXPLICIT stacking, then we8169 ** do not use the probabilities (but we still assume they are there8170 ** to be read) but generate a warning.8171 **8172 * do 30 i = 1, n_layers8173 * sum = ZERO8174 * do 40 j = 1, n_layers8175 * sum = sum + l_alpha(j,i)8176 * 40 continue8177 * if(abs(sum  ONE).gt.eps6) then8178 * write(messge,101)8179 *  'Stacking probabilities from LAYER ',i,' do not sum to 1.$'8180 * if(recrsv.or.rndm) then8181 * goto 9808182 * else if(xplcit .and. .not.rndm) then8183 * write(op,102) messge(1:index(messge, '$')  1)8184 * endif8185 * endif8186 * 30 continue8187 **8188 * RDTRNS = .true.8189 * 980 return8190 * 999 call WRTLNE(line)8191 * return8192 * 100 format(a, '(', i2, ',', i2, ')', a)8193 * 101 format(a, i2, a)8194 * 102 format(1x, 'WARNING: ', a)8195 * end8196 **8197 ** ______________________________________________________________________8198 ** Title: RDWAVL8199 ** Authors: MMJT8200 ** Date: 4 Oct 19898201 ** Description: This function reads the radiation wavelength.8202 **8203 ** ARGUMENTS:8204 ** unit_no  Logical unit number that the data file is to8205 ** be read from. (input).8206 ** messge  A short message indicating what went wrong8207 ** (if anything) during the datafile read.8208 ** messge is terminated by a token '$'. (output).8209 **8210 ** COMMON VARIABLES:8211 ** modifies: lambda8212 **8213 ** RDWAVL returns logical .true. if an acceptable value8214 ** 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_no8222 * character*(*) messge8223 **8224 * character*200 line8225 **8226 ** external subroutines (Some compilers need them declared external)8227 ** external GETLNE, WRTLNE8228 **8229 * RDWAVL = .false.8230 **8231 * messge = 'lambda is specified improperly.$'8232 * call GETLNE(unit_no, line, *999)8233 **8234 ** for strict FORTRAN77 compliance8235 * write(scrtch, *, err=990) line8236 * rewind(scrtch, err=990)8237 * read(scrtch, *, err=999) lambda8238 * rewind(scrtch, err=990)8239 **8240 ** for more lenient, efficient compilers8241 ** read(line, *, err=999) lambda8242 **8243 * if(lambda.le.ZERO) then8244 * messge = 'lambda is negative (or zero).$'8245 * goto 9998246 * endif8247 **8248 * RDWAVL = .true.8249 * return8250 * 990 messge = 'Problems using scratch file in RDWAVL.$'8251 * 999 call WRTLNE(line)8252 * return8253 * end8254 **8255 ** ______________________________________________________________________8256 ** Title: RDWDTH8257 ** Authors: MMJT8258 ** Date: 3 Mar 19958259 ** Description: This function reads the lateral layer dimensions8260 ** from the data file. First we check that the data on the line does8261 ** not refer to the layer atomic positions. If it does, assume infinite8262 ** widths, backspace the file, and return to let RDLAYR do its job.8263 ** Unlike the other read functions, RDWDTH does not return an error in8264 ** this instance, so as to maintain compatibility with data files8265 ** written for v 1.7xx which did not treat lateral layer width8266 ** broadening.8267 **8268 ** ARGUMENTS:8269 ** unit_no  Logical unit number that the data file is to8270 ** be read from. (input).8271 ** messge  A short message indicating what went wrong8272 ** (if anything) during the datafile read.8273 ** messge is terminated by a token '$'. (output).8274 **8275 ** COMMON VARIABLES:8276 ** modifies: finite_width, Wa, Wb8277 **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_no8286 * character*(*) messge8287 **8288 * character*200 line8289 **8290 * integer*4 n, CNTARG8291 **8292 ** external function8293 * external CNTARG8294 ** external subroutines (Some compilers need them declared external)8295 ** external GETLNE, TOUPPR, WRTLNE8296 **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 so8301 ** 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) then8306 * finite_width = .false.8307 * write(messge,100)8308 *  'Problems ''backspacing'' unit ',unit_no, '.$'8309 ** backspace file so that RDLAYR will work properly8310 * backspace(unit = unit_no, err = 990)8311 * goto 9008312 * endif8313 **8314 ** assume that layer width data is specified8315 * if(index(line,'INFINITE').eq.1) then8316 * finite_width = .false.8317 * goto 9008318 * else8319 * finite_width = .true.8320 * endif8321 **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 compliance8327 * write(scrtch, *, err=990) line8328 * rewind(scrtch, err=990)8329 **8330 * if(n.eq.1) then8331 ** for strict FORTRAN77 compliance8332 * read(scrtch, *, err=999) Wa8333 * rewind(scrtch, err=990)8334 ** for more lenient, efficient compilers8335 ** read(line, *, err=999) Wa8336 * Wb = Wa8337 * else if(n.eq.2) then8338 ** for strict FORTRAN77 compliance8339 * read(scrtch, *, err=999) Wa, Wb8340 * rewind(scrtch, err=990)8341 ** for more lenient, efficient compilers8342 ** read(line, *, err=999) Wa, Wb8343 * else8344 * goto 9998345 * endif8346 **8347 * if(finite_width) then8348 * if(Wa.gt.inf_width .and. Wb.gt.inf_width) then8349 * finite_width = .false.8350 * write(op,101) 'Layers will be treated as if infinitely wide.'8351 * endif8352 * if(Wa.le.ZERO .or. Wb.le.ZERO) then8353 * write(op,101)'Illegal values for layer characteristic widths'8354 * goto 9998355 * endif8356 * endif8357 **8358 * 900 RDWDTH = .true.8359 * return8360 * 990 messge = 'Problems using scratch file in RDWDTH.$'8361 * 999 call WRTLNE(line)8362 * return8363 * 100 format(a, i2, a)8364 * 101 format(1x, a)8365 * end8366 **8367 ** ______________________________________________________________________8368 5029 * Title: RAN3 8369 5030 * Authors: Press, Flannery, Teukolsky and Vetterling … … 8433 5094 * 8434 5095 * ______________________________________________________________________ 8435 * Title: SALUTE8436 * Author: MMJT8437 * Date: 19th May, 20108438 * 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 * save8447 *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 return8473 1 format(1x, a51)8474 end8475 *8476 * ______________________________________________________________________8477 ** Title: SFC8478 ** Authors: MWD and MMJT8479 ** Date: 11 Feb 90; 7 Mar 19958480 ** Description: This routine reads the atomic scattering factor constants8481 ** from the file 'sfname' (usually given as "data.sfc"). The scattering8482 ** factors are used to generate the layer form factors. SFC returns8483 ** .false. if there are too many distinct atom types (ie. more than8484 ** 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_type8491 ** NEUTRN, X_RAY, rad_type, sfname8492 **8493 ** modifies: n_atoms, x_sf, n_sf, e_sf, a_type8494 **8495 ** SFC returns logical .true. if it read the scattering factor data8496 ** safely.8497 ** ______________________________________________________________________8498 **8499 * logical function SFC()8500 * include 'DIFFaX.par'8501 * include 'DIFFaX.inc'8502 **8503 * logical ok, new_atom, our_atom, done8504 * logical list(MAX_TA+1)8505 * integer*4 i, j, n, m, LENGTH8506 * character*4 name8507 * character*120 line8508 * real*4 tmp8509 **8510 ** external subroutine (Some compilers need them declared external)8511 ** external ATOMS8512 ** external function8513 * external LENGTH8514 **8515 * write(op,301) 'Reading scattering factor datafile ''',8516 *  sfname(1:LENGTH(sfname)),'''. . .'8517 **8518 * SFC = .false.8519 * do 10 i = 1, MAX_TA8520 * list(i) = .false.8521 * atom_l(i) = ' '8522 * 10 continue8523 * m = 08524 ** determine all distinct atom types8525 * do 30 i = 1, n_actual8526 * 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 already8530 * n = 08531 * 50 n = n + 18532 * new_atom = name.ne.atom_l(n)8533 * if(new_atom .and. n.lt.m) goto 508534 * if(new_atom) then8535 * m = m + 18536 * if(m.gt.MAX_TA) goto 2208537 * atom_l(m) = name8538 * a_type(j,i) = m8539 * else8540 * a_type(j,i) = n8541 * endif8542 * 40 continue8543 * 30 continue8544 * n_atoms = m8545 ** now find data for each atom type in file8546 ** pass through file only once8547 * 60 read(sf, 300, end=90, err=200) line8548 * name = line(1:4)8549 * call TOUPPR(name)8550 ** see if this is one of the distinct atoms8551 * i = 08552 * 70 i = i + 18553 * our_atom = name.eq.atom_l(i)8554 * if(.not.our_atom .and. i.lt.n_atoms) goto 708555 ** have we read all that we need to know?8556 * done = .true.8557 * do 80 n = 1, n_atoms8558 * done = done.and.list(n)8559 * 80 continue8560 ** If we're done, close file.8561 * if(done) goto 908562 * if(our_atom .and. .not.list(i)) then8563 ** mark this atom's data as read in8564 * list(i) = .true.8565 ** and read it in8566 * if(rad_type.eq.X_RAY) then8567 * read(line,310,err=210) (x_sf(j,i), j=1, 9)8568 * else if(rad_type.eq.NEUTRN) then8569 * read(line,320,err=210) n_sf(i)8570 * else8571 * read(line,310,err=210) (x_sf(j,i), j=1, 9), tmp, e_sf(i)8572 * endif8573 * endif8574 * goto 608575 * 90 close(sf,err=240)8576 ** see if all the data for each atom has been read in8577 * ok = .true.8578 * do 100 i = 1, n_atoms8579 * if(.not.list(i)) then8580 * 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 * endif8585 * 100 continue8586 * if(ok) then8587 * SFC = .true.8588 * write(op,400) 'Scattering factor data read in.'8589 * endif8590 * return8591 * 200 write(op,301) 'Scattering factor file ''',8592 *  sfname(1:LENGTH(sfname)), ''' defective.'8593 * close(sf,err=240)8594 * return8595 * 210 write(op,400) 'ERROR reading scattering factor data.'8596 * close(sf,err=240)8597 * return8598 * 220 write(op,402) 'There are too many types of atoms in layer ', i8599 * write(op,400) 'Atoms recorded so far are:'8600 * do 110 j = 1, MAX_TA8601 * write(op,403) ' Atom type ', j, ' ', atom_l(j)8602 * 110 continue8603 * write(op,402) ' Maximum number of types allowed is ', MAX_TA8604 * return8605 * 240 write(op,301) 'Unable to close scattering factor file ''',8606 *  sfname(1:LENGTH(sfname)), '''.'8607 * return8608 * 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 * end8618 **8619 ** ______________________________________________________________________8620 5096 * Title: SHARP 8621 5097 * Author: MMJT … … 8803 5279 * 8804 5280 * ______________________________________________________________________ 8805 ** Title: STREAK8806 ** Author: MWD & MMJT8807 ** 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 file8810 ** strkfile.8811 **8812 ** ARGUMENTS:8813 ** FN  Function name passed by reference. The8814 ** choice is between GLQ16 (nonadaptive8815 ** GaussLegendre integration), and AGLQ168816 ** (adaptive GaussLegendre integration).8817 ** (input).8818 ** strkfile  The name of the file that streak data is8819 ** 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_RAY8826 **8827 ** modifies: no COMMON variables are modified8828 ** ______________________________________________________________________8829 **8830 * subroutine STREAK(FN, strkfile, ok)8831 * include 'DIFFaX.par'8832 * include 'DIFFaX.inc'8833 **8834 * real*8 FN8835 * character*(*) strkfile8836 * logical ok8837 **8838 * logical its_hot8839 * integer*4 h, k, LENGTH, i, i_step8840 * real*8 l, theta, x, ANGLE, S, Q2, l0, l18841 * real*8 dl, W4, intervals8842 * parameter (intervals = TWENTY)8843 **8844 ** external functions8845 * external FN, LENGTH8846 ** external subroutines (Some compilers need them declared external)8847 ** external XYPHSE, PRE_MAT, GET_F8848 **8849 ** statement functions8850 ** S is the value of 1/d**2 at hkl8851 * S(h,k,l) = h*h*a0 + k*k*b0 + l*l*c0 + h*k*d08852 ** ANGLE is the Bragg angle (in radians) of the h,k,l plane8853 * ANGLE(h,k,l) = asin(HALF * lambda * sqrt(S(h,k,l)))8854 ** W4 is the Xray polarization factor8855 * 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, dl8860 * if(CFile) write(op,401) h, k, l0, l1, dl8861 ** check input8862 * if(l1.eq.l0) then8863 * write(op,400) 'Illegal input: l0 equals l1'8864 * goto 9998865 * else if(dl.eq.ZERO) then8866 * write(op,400) 'Illegal zero value of dl entered'8867 * write(op,402)'A value of ',(l1l0)/(FIVE*HUNDRED),' is assumed'8868 * else if(l1.gt.l0 .and. dl.lt.ZERO) then8869 * write(op,400) 'l1 is greater than l0. +ve dl assumed'8870 * dl = dl8871 * else if(l1.lt.l0 .and. dl.gt.ZERO) then8872 * write(op,400) 'l0 is greater than l1. ve dl assumed'8873 * dl = dl8874 * endif8875 ** 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.ELECTN8878 * if(its_hot) then8879 * write(op,400) 'Cannot scan the origin with electron radiation'8880 * write(op,400) 'Origin will be skipped.'8881 * endif8882 ** check angles are meaningful8883 * if(S(h,k,l0).gt.Q2 .or. S(h,k,l1).gt.Q2) then8884 * 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 108889 * endif8890 **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 = 108899 * i = 08900 * do 20 l = l0, l1, dl8901 ** If we are dealing with electrons, make sure we avoid the origin8902 * if(its_hot .and. l*(l+dl).le.ZERO) then8903 * x = ZERO8904 * goto 308905 * endif8906 * i = i + 18907 * if(mod(i,i_step).eq.0) write(op,405) 'Reached l = ',l8908 * x = FN(h,k,l,l+dl,ok)8909 * if(.not.ok) goto 9998910 ** note: since this is streak data, only the X_RAY input needs8911 ** 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), x8914 * 20 continue8915 * if(sk.ne.op) close(sk,err=110)8916 * write(op,404) 'Streak data file, ''',8917 *  strkfile(1:LENGTH(strkfile)),''' written to disk.'8918 * return8919 * 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 * return8923 * 110 write(op,404) 'Unable to close streak data file ''',8924 *  strkfile(1:LENGTH(strkfile)),''''8925 * return8926 * 999 write(op,405) 'ERROR encountered in streak integration at l = ',l8927 * return8928 * 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 * end8936 **8937 ** ______________________________________________________________________8938 5281 * Title: THRESH 8939 5282 * Author: MMJT … … 9019 5362 * 9020 5363 * ______________________________________________________________________ 9021 ** Title: TOUPPR9022 ** Author: MMJT9023 ** Date: 3 August 19899024 ** Description: Converts alphabetic characters in 'line' to uppercase.9025 ** Works for ascii characters, since it assumes alphabetic characters9026 ** are contiguous.9027 ** The neat trick of using9028 ** lower_case = ichar(chr).and.959029 ** will not always work in ansi77 fortran, as some compilers evaluate9030 ** the right hand side as a logical .true. or .false. rather than a9031 ** bitwise 'anded' integer.9032 ** The alternative9033 ** lower_case = ichar(chr).iand.959034 ** is equally undesirable since the operator .iand. is a military9035 ** extension, and is not ansi77.9036 **9037 ** ARGUMENTS:9038 ** line  Input line of characters to be converted to9039 ** upper case. The converted string is returned9040 ** in 'line'.9041 ** ______________________________________________________________________9042 **9043 * subroutine TOUPPR(line)9044 **9045 * character*(*) line9046 **9047 * integer*4 a, z, offset, i, itmp, lin_len9048 **9049 * a = ichar('a')9050 * z = ichar('z')9051 * offset = a  ichar('A')9052 * lin_len = len(line)9053 **9054 * i = 19055 * 10 itmp = ichar(line(i:i))9056 * if(itmp.ge.a .and. itmp.le.z) line(i:i) = char(itmp  offset)9057 * i = i + 19058 * if(i.le.lin_len)9059 * goto 109060 **9061 * return9062 * end9063 **9064 ** ______________________________________________________________________9065 5364 * Title: TRMSPC 9066 5365 * Author: MMJT … … 9632 5931 * 9633 5932 * ______________________________________________________________________ 9634 ** Title: WRTLNE9635 ** Author: MMJT9636 ** Date: 4 Oct 19899637 ** Description: This subroutine writes out the contents of9638 ** 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 ** save9648 **9649 * character*(*) line9650 **9651 * integer*4 lin_len, i9652 **9653 * lin_len = len(line)9654 * i = lin_len9655 * 10 if(i.gt.0) then9656 * if(line(i:i).eq.' ') then9657 * if(i.gt.1) then9658 * i = i  19659 * goto 109660 * endif9661 * endif9662 * endif9663 **9664 * write(op,100) 'The data on the current line was read as:'9665 * write(op,101) ' ''', line(1:i), ''''9666 **9667 * return9668 * 100 format(1x, a)9669 * 101 format(1x, 3a)9670 * end9671 **9672 ** ______________________________________________________________________9673 ** Title: WRTSAD9674 ** Author: MMJT9675 ** Date: 29 Oct 1989; 1st Mar 2000; 19th May 20109676 ** Description: This subroutine writes the selected area diffraction9677 ** pattern (SADP) data to a binary file.9678 ** Binary file output has been problematical because some compilers add9679 ** recordlength specifiers at the beginning and ending of each line.9680 ** Sometimes the specifiers are 2byte, or 4byte, and even 0byte!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 binary9685 ** selected area diffraction data is to be9686 ** written to. (input).9687 ** view  Choice of beam direction. (input).9688 ** 1 = normal to the plane k = 09689 ** 2 = normal to the plane h = 09690 ** 3 = normal to the plane h = k9691 ** 4 = normal to the plane h = k9692 ** 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_mirror9699 **9700 ** modifies: spec9701 ** ______________________________________________________________________9702 **9703 * subroutine WRTSAD(outfile, view, l_upper, hk_lim, ok)9704 * include 'DIFFaX.par'9705 * include 'DIFFaX.inc'9706 **9707 * character*(*) outfile9708 * integer*4 view, hk_lim9709 * real*8 l_upper9710 * logical ok9711 **9712 * integer*4 i, j, p1, p2, LENGTH9713 * real*8 rowint(SADSIZE), incr, sigma, x9714 **9715 ** external function9716 * external LENGTH9717 **9718 ** external subroutine (Some compilers need them declared external)9719 ** external SMUDGE9720 **9721 * write(op,102) 'Writing SADP data to file ''',9722 *  outfile(1:LENGTH(outfile)), '''. . .'9723 **9724 ** set standard deviation9725 * sigma = ZERO9726 **9727 ** Establish scaling in reciprocal space, the number of pixels per unit9728 ** first get number of pixels per l9729 * incr = dble(SADSIZE/2) / l_upper9730 * if(view.eq.1) then9731 ** number of pixels per unit h9732 * incr = incr * sqrt(a0 / c0)9733 * else if(view.eq.2) then9734 ** number of pixels per unit k9735 * incr = incr * sqrt(b0 / c0)9736 * else if(view.eq.3) then9737 ** number of pixels per unit along (h = k)9738 * incr = incr * sqrt((a0 + b0 + d0) / c0)9739 * else if(view.eq.4) then9740 ** number of pixels per unit along (h = k)9741 * incr = incr * sqrt((a0 + b0  d0) / c0)9742 * endif9743 **9744 * do 10 j = sadblock  1, 0, 19745 * do 20 i = 1, SADSIZE9746 * rowint(i) = ZERO9747 * 20 continue9748 * do 30 i = 0, hk_lim9749 * p1 = SADSIZE/2 + nint(i*incr) + 19750 * p2 = SADSIZE/2  nint(i*incr) + 19751 **9752 ** cycle if we have gone out of bounds9753 * if(p1.gt.SADSIZE .or. p1.lt.0) then9754 * goto 309755 * endif9756 * if(p2.gt.SADSIZE .or. p2.lt.0) then9757 * goto 309758 * endif9759 **9760 * x = spec(i*sadblock + j + 1) * scaleint9761 ** handle saturated pixels9762 * if(x.gt.maxsad) x = maxsad9763 * rowint(p1) = x9764 * if(has_l_mirror) then9765 * rowint(p2) = x9766 * else9767 * x = spec(i*sadblock + sadblock  j) * scaleint9768 ** handle saturated pixels9769 * if(x.gt.maxsad) x = maxsad9770 * rowint(p2) = x9771 * endif9772 * 30 continue9773 * call SMUDGE(rowint, SADSIZE, sigma, ok)9774 * if(.not.ok) goto 9999775 * if(bitdepth.eq.8) then9776 * write(sa,err=900) (char(nint(rowint(i))), i = 1, SADSIZE)9777 * else9778 ** Known bug: Writing 16bit binary is problematic on Intel processors9779 * 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 * endif9784 * 10 continue9785 **9786 ** Repeat, in reverse for the bottom half if data had a mirror.9787 * if(has_l_mirror) then9788 * do 40 j = 1, sadblock  19789 * do 50 i = 1, SADSIZE9790 * rowint(i) = ZERO9791 * 50 continue9792 * do 60 i = 0, hk_lim9793 * p1 = SADSIZE/2 + nint(i*incr) + 19794 * p2 = SADSIZE/2  nint(i*incr) + 19795 **9796 ** cycle if we have gone out of bounds9797 * if(p1.gt.SADSIZE .or. p1.lt.0) then9798 * goto 609799 * endif9800 * if(p2.gt.SADSIZE .or. p2.lt.0) then9801 * goto 609802 * endif9803 **9804 * x = spec(i*sadblock + j + 1) * scaleint9805 ** handle saturated pixels9806 * if(x.gt.maxsad) x = maxsad9807 * rowint(p1) = x9808 * rowint(p2) = x9809 * 60 continue9810 * call SMUDGE(rowint, SADSIZE, sigma, ok)9811 * if(.not.ok) goto 9999812 * if(bitdepth.eq.8) then9813 * write(sa,err=900) (char(nint(rowint(i))), i = 1, SADSIZE)9814 * else9815 ** Known bug: Writing 16bit binary is problematic on Intel processors9816 * 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 * endif9821 * 40 continue9822 ** write a blank last line to make the array SADSIZE x SADSIZE square9823 * if(bitdepth.eq.8) then9824 * write(sa,err=900) (char(0), i = 1, SADSIZE)9825 * else9826 * write(sa,err=900) (char(0), i = 1, 2*SADSIZE)9827 * endif9828 * endif9829 **9830 * if(sa.ne.op) close(unit = sa, err = 990)9831 **9832 ** MMJT 19th May 20109833 * write(op,103) SADSIZE, ' x ', SADSIZE, ' pixels: ',9834 *  bitdepth, ' bits deep.'9835 * write(op,100) ' (Caution: Because of a compilerdependent issue'9836 * write(op,100) ' the true image width in pixels is found by'9837 * if(bitdepth.eq.8) then9838 * write(op,100) ' dividing the file size in bytes by 256.)'9839 * else9840 * write(op,100) ' dividing the file size in bytes by 512.'9841 * write(op,100) ' The byteorder is bigendian.)'9842 * endif9843 **9844 * return9845 * 900 write(op,100) 'ERROR: problems writing to binary SADP file.'9846 * ok = .false.9847 * return9848 * 990 write(op,100) 'ERROR: problems closing binary SADP file.'9849 * ok = .false.9850 * return9851 * 999 write(op,100) 'ERROR: SMUDGE returned error to WRTSAD.'9852 * write(op,101) i, j9853 * return9854 * 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 * end9859 ** ______________________________________________________________________9860 ** Title: WRTSPC9861 ** Author: MWD and MMJT9862 ** Date: 18 Aug 1988; 7 Mar 19959863 ** 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, blurring9872 ** NONE, RAD2DEG9873 **9874 ** modifies: no COMMON variables are modified9875 ** ______________________________________________________________________9876 **9877 * subroutine WRTSPC(spcfile, ok)9878 * include 'DIFFaX.par'9879 * include 'DIFFaX.inc'9880 **9881 * character*(*) spcfile9882 * logical ok9883 **9884 * character*1 tab9885 * integer*4 i, n_low, n_high, LENGTH9886 * real*8 theta9887 **9888 ** external function9889 * external LENGTH9890 **9891 * n_low = 19892 * n_high = int(HALF*(th2_maxth2_min)/d_theta) + 19893 **9894 * tab = char(9)9895 * theta = th2_min * RAD2DEG9896 * 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) + 19900 * if(blurring.eq.NONE) then9901 * do 10 i = n_low, n_high9902 * write(sp,101,err=50) theta, tab, spec(i)9903 * theta = theta + TWO * RAD2DEG * d_theta9904 * 10 continue9905 * else9906 * do 20 i = n_low, n_high9907 * write(sp,100,err=50) theta, tab, spec(i), tab, brd_spc(i)9908 * theta = theta + TWO * RAD2DEG * d_theta9909 * 20 continue9910 * endif9911 **9912 * if(sp.ne.op) close(sp,err=60)9913 * write(op,202) 'Spectrum written.'9914 * write(op,202)9915 * return9916 * 50 write(op,202) 'Unable to write to spectrum file.'9917 * if(sp.ne.op) close(sp,err=60)9918 * ok = .false.9919 * return9920 * 60 write(op,202) 'Unable to close spectrum file.'9921 * ok = .false.9922 * return9923 * 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 * end9929 **9930 ** ______________________________________________________________________9931 5933 * Title: XYPHSE 9932 5934 * Author: MMJT
Note: See TracChangeset
for help on using the changeset viewer.