Changeset 495


Ignore:
Timestamp:
Feb 24, 2012 2:28:18 PM (10 years ago)
Author:
vondreele
Message:

back to mine

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIstruct.py

    r494 r495  
    11#GSASIIstructure - structure computation routines
    2 
    32########### SVN repository information ###################
    4 
    53# $Date$
    6 
    74# $Author$
    8 
    95# $Revision$
    10 
    116# $URL$
    12 
    137# $Id$
    14 
    158########### SVN repository information ###################
    16 
    17 MaxThreads = 1
    18 
    19 MaxProcess = 2
    20 
    21 
    22 
    239import sys
    24 
    2510import os
    26 
    2711import os.path as ospath
    28 
    2912import time
    30 
    3113import math
    32 
    3314import cPickle
    34 
    3515import multiprocessing as mp
    36 
    3716import numpy as np
    38 
    3917import numpy.linalg as nl
    40 
    4118import scipy.optimize as so
    42 
    4319import GSASIIpath
    44 
    4520import GSASIIElem as G2el
    46 
    4721import GSASIIlattice as G2lat
    48 
    4922import GSASIIspc as G2spc
    50 
    5123import GSASIIpwd as G2pwd
    52 
    5324import GSASIImapvars as G2mv
    54 
    5525import GSASIImath as G2mth
    5626
    57 #import buggerwx
    58 
    59 try:
    60 
    61     import mkl
    62 
    63 except:
    64 
    65     print "MKL module not present"
    66 
    67 
    68 
    69 
    70 
    7127sind = lambda x: np.sin(x*np.pi/180.)
    72 
    7328cosd = lambda x: np.cos(x*np.pi/180.)
    74 
    7529tand = lambda x: np.tan(x*np.pi/180.)
    76 
    7730asind = lambda x: 180.*np.arcsin(x)/np.pi
    78 
    7931acosd = lambda x: 180.*np.arccos(x)/np.pi
    80 
    8132atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
    8233
    8334
    84 
    85 
    86 
    8735def GetControls(GPXfile):
    88 
    8936    ''' Returns dictionary of control items found in GSASII gpx file
    90 
    9137    input:
    92 
    9338        GPXfile = .gpx full file name
    94 
    9539    return:
    96 
    9740        Controls = dictionary of control items
    98 
    9941    '''
    100 
    101     Controls = {'deriv type':'analytical','min dM/M':0.0001,'shift factor':1.}
    102 
     42    Controls = {'deriv type':'analytic Hessian','max cyc':3,'max Hprocess':1,
     43        'max Rprocess':1,'min dM/M':0.0001,'shift factor':1.}
    10344    file = open(GPXfile,'rb')
    104 
    10545    while True:
    106 
    10746        try:
    108 
    10947            data = cPickle.load(file)
    110 
    11148        except EOFError:
    112 
    11349            break
    114 
    11550        datum = data[0]
    116 
    11751        if datum[0] == 'Controls':
    118 
    11952            Controls.update(datum[1])
    120 
    12153    file.close()
    122 
    12354    return Controls
    124 
    125    
    126 
     55   
    12756def GetConstraints(GPXfile):
    128 
    12957    constList = []
    130 
    13158    file = open(GPXfile,'rb')
    132 
    13359    while True:
    134 
    13560        try:
    136 
    13761            data = cPickle.load(file)
    138 
    13962        except EOFError:
    140 
    14163            break
    142 
    14364        datum = data[0]
    144 
    14565        if datum[0] == 'Constraints':
    146 
    14766            constDict = datum[1]
    148 
    14967            for item in constDict:
    150 
    15168                constList += constDict[item]
    152 
    15369    file.close()
    154 
    15570    constDict = []
    156 
    15771    constFlag = []
    158 
    15972    fixedList = []
    160 
    16173    for item in constList:
    162 
    16374        if item[-2]:
    164 
    16575            fixedList.append(str(item[-2]))
    166 
    16776        else:
    168 
    16977            fixedList.append('0')
    170 
    17178        if item[-1]:
    172 
    17379            constFlag.append(['VARY'])
    174 
    17580        else:
    176 
    17781            constFlag.append([])
    178 
    17982        itemDict = {}
    180 
    18183        for term in item[:-2]:
    182 
    18384            itemDict[term[1]] = term[0]
    184 
    18585        constDict.append(itemDict)
    186 
    18786    return constDict,constFlag,fixedList
    188 
    189    
    190 
     87   
    19188def GetPhaseNames(GPXfile):
    192 
    19389    ''' Returns a list of phase names found under 'Phases' in GSASII gpx file
    194 
    19590    input:
    196 
    19791        GPXfile = gpx full file name
    198 
    19992    return:
    200 
    20193        PhaseNames = list of phase names
    202 
    20394    '''
    204 
    20595    file = open(GPXfile,'rb')
    206 
    20796    PhaseNames = []
    208 
    20997    while True:
    210 
    21198        try:
    212 
    21399            data = cPickle.load(file)
    214 
    215100        except EOFError:
    216 
    217101            break
    218 
    219102        datum = data[0]
    220 
    221103        if 'Phases' == datum[0]:
    222 
    223104            for datus in data[1:]:
    224 
    225105                PhaseNames.append(datus[0])
    226 
    227106    file.close()
    228 
    229107    return PhaseNames
    230108
    231 
    232 
    233109def GetAllPhaseData(GPXfile,PhaseName):
    234 
    235110    ''' Returns the entire dictionary for PhaseName from GSASII gpx file
    236 
    237111    input:
    238 
    239112        GPXfile = gpx full file name
    240 
    241113        PhaseName = phase name
    242 
    243114    return:
    244 
    245115        phase dictionary
    246 
    247116    '''       
    248 
    249117    file = open(GPXfile,'rb')
    250 
    251118    General = {}
    252 
    253119    Atoms = []
    254 
    255120    while True:
    256 
    257121        try:
    258 
    259122            data = cPickle.load(file)
    260 
    261123        except EOFError:
    262 
    263124            break
    264 
    265125        datum = data[0]
    266 
    267126        if 'Phases' == datum[0]:
    268 
    269127            for datus in data[1:]:
    270 
    271128                if datus[0] == PhaseName:
    272 
    273129                    break
    274 
    275130    file.close()
    276 
    277131    return datus[1]
    278 
    279    
    280 
     132   
    281133def GetHistograms(GPXfile,hNames):
    282 
    283134    """ Returns a dictionary of histograms found in GSASII gpx file
    284 
    285135    input:
    286 
    287136        GPXfile = .gpx full file name
    288 
    289137        hNames = list of histogram names
    290 
    291138    return:
    292 
    293139        Histograms = dictionary of histograms (types = PWDR & HKLF)
    294 
    295140    """
    296 
    297141    file = open(GPXfile,'rb')
    298 
    299142    Histograms = {}
    300 
    301143    while True:
    302 
    303144        try:
    304 
    305145            data = cPickle.load(file)
    306 
    307146        except EOFError:
    308 
    309147            break
    310 
    311148        datum = data[0]
    312 
    313149        hist = datum[0]
    314 
    315150        if hist in hNames:
    316 
    317151            if 'PWDR' in hist[:4]:
    318 
    319152                PWDRdata = {}
    320 
    321153                PWDRdata['Data'] = datum[1][1]          #powder data arrays
    322 
    323154                PWDRdata[data[2][0]] = data[2][1]       #Limits
    324 
    325155                PWDRdata[data[3][0]] = data[3][1]       #Background
    326 
    327156                PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
    328 
    329157                PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
    330 
    331158                try:
    332 
    333159                    PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
    334 
    335160                except IndexError:
    336 
    337161                    PWDRdata['Reflection lists'] = {}
    338 
    339    
    340 
     162   
    341163                Histograms[hist] = PWDRdata
    342 
    343164            elif 'HKLF' in hist[:4]:
    344 
    345165                HKLFdata = []
    346 
    347166                datum = data[0]
    348 
    349167                HKLFdata = datum[1:][0]
    350 
    351168                Histograms[hist] = HKLFdata           
    352 
    353169    file.close()
    354 
    355170    return Histograms
    356 
    357    
    358 
     171   
    359172def GetHistogramNames(GPXfile,hType):
    360 
    361173    """ Returns a list of histogram names found in GSASII gpx file
    362 
    363174    input:
    364 
    365175        GPXfile = .gpx full file name
    366 
    367176        hType = list ['PWDR','HKLF']
    368 
    369177    return:
    370 
    371178        HistogramNames = list of histogram names (types = PWDR & HKLF)
    372 
    373179    """
    374 
    375180    file = open(GPXfile,'rb')
    376 
    377181    HistogramNames = []
    378 
    379182    while True:
    380 
    381183        try:
    382 
    383184            data = cPickle.load(file)
    384 
    385185        except EOFError:
    386 
    387186            break
    388 
    389187        datum = data[0]
    390 
    391188        if datum[0][:4] in hType:
    392 
    393189            HistogramNames.append(datum[0])
    394 
    395190    file.close()
    396 
    397191    return HistogramNames
    398 
    399    
    400 
     192   
    401193def GetUsedHistogramsAndPhases(GPXfile):
    402 
    403194    ''' Returns all histograms that are found in any phase
    404 
    405195    and any phase that uses a histogram
    406 
    407196    input:
    408 
    409197        GPXfile = .gpx full file name
    410 
    411198    return:
    412 
    413199        Histograms = dictionary of histograms as {name:data,...}
    414 
    415200        Phases = dictionary of phases that use histograms
    416 
    417201    '''
    418 
    419202    phaseNames = GetPhaseNames(GPXfile)
    420 
    421203    histoList = GetHistogramNames(GPXfile,['PWDR','HKLF'])
    422 
    423204    allHistograms = GetHistograms(GPXfile,histoList)
    424 
    425205    phaseData = {}
    426 
    427206    for name in phaseNames:
    428 
    429207        phaseData[name] =  GetAllPhaseData(GPXfile,name)
    430 
    431208    Histograms = {}
    432 
    433209    Phases = {}
    434 
    435210    for phase in phaseData:
    436 
    437211        Phase = phaseData[phase]
    438 
    439212        if Phase['Histograms']:
    440 
    441213            if phase not in Phases:
    442 
    443214                pId = phaseNames.index(phase)
    444 
    445215                Phase['pId'] = pId
    446 
    447216                Phases[phase] = Phase
    448 
    449217            for hist in Phase['Histograms']:
    450 
    451218                if hist not in Histograms:
    452 
    453219                    Histograms[hist] = allHistograms[hist]
    454 
    455220                    #future restraint, etc. histograms here           
    456 
    457221                    hId = histoList.index(hist)
    458 
    459222                    Histograms[hist]['hId'] = hId
    460 
    461223    return Histograms,Phases
    462 
    463    
    464 
     224   
     225def getBackupName2(GPXfile,makeBack=True):      #not work correctly
     226    GPXpath,GPXname = ospath.split(GPXfile)
     227    if GPXpath == '': GPXpath = '.'
     228    Name = ospath.splitext(GPXname)[0]
     229    files = os.listdir(GPXpath)
     230    last = 0
     231    for name in files:
     232        name = name.split('.')
     233        if len(name) >= 3 and name[0] == Name and 'bak' in name[-2]:
     234            if makeBack:
     235                last = max(last,int(name[-2].strip('bak'))+1)
     236            else:
     237                last = max(last,int(name[-2].strip('bak')))
     238    GPXback = ospath.join(GPXpath,GPXname.rstrip('.'.join(name[-2:]))+'.bak'+str(last)+'.gpx')
     239    return GPXback
     240
     241def getBackupName(GPXfile,makeBack):       #recovered old one
     242    GPXpath,GPXname = ospath.split(GPXfile)
     243    if GPXpath == '': GPXpath = '.'
     244    Name = ospath.splitext(GPXname)[0]
     245    files = os.listdir(GPXpath)
     246    last = 0
     247    for name in files:
     248        name = name.split('.')
     249        if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
     250            if makeBack:
     251                last = max(last,int(name[1].strip('bak'))+1)
     252            else:
     253                last = max(last,int(name[1].strip('bak')))
     254    GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
     255    return GPXback   
     256       
    465257def GPXBackup(GPXfile,makeBack=True):
    466 
    467258    import distutils.file_util as dfu
    468 
    469     GPXpath,GPXname = ospath.split(GPXfile)
    470 
    471     if GPXpath == '': GPXpath = '.'
    472 
    473     Name = ospath.splitext(GPXname)[0]
    474 
    475     files = os.listdir(GPXpath)
    476 
    477     last = 0
    478 
    479     for name in files:
    480 
    481         name = name.split('.')
    482 
    483         if len(name) >= 3 and name[0] == Name and 'bak' in name[-2]:
    484 
    485             if makeBack:
    486 
    487                 last = max(last,int(name[-2].strip('bak'))+1)
    488 
     259    GPXback = getBackupName(GPXfile,makeBack)
     260    dfu.copy_file(GPXfile,GPXback)
     261    return GPXback
     262
     263def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,CovData,makeBack=True):
     264    ''' Updates gpxfile from all histograms that are found in any phase
     265    and any phase that used a histogram
     266    input:
     267        GPXfile = .gpx full file name
     268        Histograms = dictionary of histograms as {name:data,...}
     269        Phases = dictionary of phases that use histograms
     270        CovData = dictionary of refined variables, varyList, & covariance matrix
     271        makeBack = True if new backup of .gpx file is to be made; else use the last one made
     272    '''
     273                       
     274    GPXback = GPXBackup(GPXfile,makeBack)
     275    print '\n',135*'-'
     276    print 'Read from file:',GPXback
     277    print 'Save to file  :',GPXfile
     278    infile = open(GPXback,'rb')
     279    outfile = open(GPXfile,'wb')
     280    while True:
     281        try:
     282            data = cPickle.load(infile)
     283        except EOFError:
     284            break
     285        datum = data[0]
     286#        print 'read: ',datum[0]
     287        if datum[0] == 'Phases':
     288            for iphase in range(len(data)):
     289                if data[iphase][0] in Phases:
     290                    phaseName = data[iphase][0]
     291                    data[iphase][1].update(Phases[phaseName])
     292        elif datum[0] == 'Covariance':
     293            data[0][1] = CovData
     294        try:
     295            histogram = Histograms[datum[0]]
     296#            print 'found ',datum[0]
     297            data[0][1][1] = histogram['Data']
     298            for datus in data[1:]:
     299#                print '    read: ',datus[0]
     300                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
     301                    datus[1] = histogram[datus[0]]
     302        except KeyError:
     303            pass
     304                               
     305        cPickle.dump(data,outfile,1)
     306    infile.close()
     307    outfile.close()
     308    print 'GPX file save successful'
     309   
     310def SetSeqResult(GPXfile,Histograms,SeqResult):
     311    GPXback = GPXBackup(GPXfile)
     312    print '\n',135*'-'
     313    print 'Read from file:',GPXback
     314    print 'Save to file  :',GPXfile
     315    infile = open(GPXback,'rb')
     316    outfile = open(GPXfile,'wb')
     317    while True:
     318        try:
     319            data = cPickle.load(infile)
     320        except EOFError:
     321            break
     322        datum = data[0]
     323        if datum[0] == 'Sequental results':
     324            data[0][1] = SeqResult
     325        try:
     326            histogram = Histograms[datum[0]]
     327            data[0][1][1] = histogram['Data']
     328            for datus in data[1:]:
     329                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
     330                    datus[1] = histogram[datus[0]]
     331        except KeyError:
     332            pass
     333                               
     334        cPickle.dump(data,outfile,1)
     335    infile.close()
     336    outfile.close()
     337    print 'GPX file save successful'
     338                       
     339def GetPWDRdata(GPXfile,PWDRname):
     340    ''' Returns powder data from GSASII gpx file
     341    input:
     342        GPXfile = .gpx full file name
     343        PWDRname = powder histogram name as obtained from GetHistogramNames
     344    return:
     345        PWDRdata = powder data dictionary with:
     346            Data - powder data arrays, Limits, Instrument Parameters, Sample Parameters
     347       
     348    '''
     349    file = open(GPXfile,'rb')
     350    PWDRdata = {}
     351    while True:
     352        try:
     353            data = cPickle.load(file)
     354        except EOFError:
     355            break
     356        datum = data[0]
     357        if datum[0] == PWDRname:
     358            PWDRdata['Data'] = datum[1][1]          #powder data arrays
     359            PWDRdata[data[2][0]] = data[2][1]       #Limits
     360            PWDRdata[data[3][0]] = data[3][1]       #Background
     361            PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
     362            PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
     363            try:
     364                PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
     365            except IndexError:
     366                PWDRdata['Reflection lists'] = {}
     367    file.close()
     368    return PWDRdata
     369   
     370def GetHKLFdata(GPXfile,HKLFname):
     371    ''' Returns single crystal data from GSASII gpx file
     372    input:
     373        GPXfile = .gpx full file name
     374        HKLFname = single crystal histogram name as obtained from GetHistogramNames
     375    return:
     376        HKLFdata = single crystal data list of reflections: for each reflection:
     377            HKLF = [np.array([h,k,l]),FoSq,sigFoSq,FcSq,Fcp,Fcpp,phase]
     378    '''
     379    file = open(GPXfile,'rb')
     380    HKLFdata = []
     381    while True:
     382        try:
     383            data = cPickle.load(file)
     384        except EOFError:
     385            break
     386        datum = data[0]
     387        if datum[0] == HKLFname:
     388            HKLFdata = datum[1:][0]
     389    file.close()
     390    return HKLFdata
     391   
     392def ShowBanner():
     393    print 80*'*'
     394    print '   General Structure Analysis System-II Crystal Structure Refinement'
     395    print '              by Robert B. Von Dreele & Brian H. Toby'
     396    print '                Argonne National Laboratory(C), 2010'
     397    print ' This product includes software developed by the UChicago Argonne, LLC,'
     398    print '            as Operator of Argonne National Laboratory.'
     399    print 80*'*','\n'
     400
     401def ShowControls(Controls):
     402    print ' Least squares controls:'
     403    print ' Refinement type: ',Controls['deriv type']
     404    if 'Hessian' in Controls['deriv type']:
     405        print ' Maximum number of cycles:',Controls['max cyc']
     406        print ' Maximum histogram processes:',Controls['max Hprocess']
     407        print ' Maximum reflection processes:',Controls['max Rprocess']
     408    else:
     409        print ' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
     410    print ' Initial shift factor: ','%.3f'%(Controls['shift factor'])
     411   
     412def GetFFtable(General):
     413    ''' returns a dictionary of form factor data for atom types found in General
     414    input:
     415        General = dictionary of phase info.; includes AtomTypes
     416    return:
     417        FFtable = dictionary of form factor data; key is atom type
     418    '''
     419    atomTypes = General['AtomTypes']
     420    FFtable = {}
     421    for El in atomTypes:
     422        FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])
     423        for item in FFs:
     424            if item['Symbol'] == El.upper():
     425                FFtable[El] = item
     426    return FFtable
     427   
     428def GetBLtable(General):
     429    ''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General
     430    input:
     431        General = dictionary of phase info.; includes AtomTypes & Isotopes
     432    return:
     433        BLtable = dictionary of scattering length data; key is atom type
     434    '''
     435    atomTypes = General['AtomTypes']
     436    BLtable = {}
     437    isotopes = General['Isotopes']
     438    isotope = General['Isotope']
     439    for El in atomTypes:
     440        BLtable[El] = [isotope[El],isotopes[El][isotope[El]]]
     441    return BLtable
     442       
     443def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
     444    if SGLaue in ['-1','2/m','mmm']:
     445        return                      #no Pawley symmetry required constraints
     446    for i,varyI in enumerate(pawleyVary):
     447        refI = int(varyI.split(':')[-1])
     448        ih,ik,il = PawleyRef[refI][:3]
     449        for varyJ in pawleyVary[0:i]:
     450            refJ = int(varyJ.split(':')[-1])
     451            jh,jk,jl = PawleyRef[refJ][:3]
     452            if SGLaue in ['4/m','4/mmm']:
     453                isum = ih**2+ik**2
     454                jsum = jh**2+jk**2
     455                if abs(il) == abs(jl) and isum == jsum:
     456                    G2mv.StoreEquivalence(varyJ,(varyI,))
     457            elif SGLaue in ['3R','3mR']:
     458                isum = ih**2+ik**2+il**2
     459                jsum = jh**2+jk**2*jl**2
     460                isum2 = ih*ik+ih*il+ik*il
     461                jsum2 = jh*jk+jh*jl+jk*jl
     462                if isum == jsum and isum2 == jsum2:
     463                    G2mv.StoreEquivalence(varyJ,(varyI,))
     464            elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
     465                isum = ih**2+ik**2+ih*ik
     466                jsum = jh**2+jk**2+jh*jk
     467                if abs(il) == abs(jl) and isum == jsum:
     468                    G2mv.StoreEquivalence(varyJ,(varyI,))
     469            elif SGLaue in ['m3','m3m']:
     470                isum = ih**2+ik**2+il**2
     471                jsum = jh**2+jk**2+jl**2
     472                if isum == jsum:
     473                    G2mv.StoreEquivalence(varyJ,(varyI,))
     474                   
     475def cellVary(pfx,SGData):
     476    if SGData['SGLaue'] in ['-1',]:
     477        return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
     478    elif SGData['SGLaue'] in ['2/m',]:
     479        if SGData['SGUniq'] == 'a':
     480            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
     481        elif SGData['SGUniq'] == 'b':
     482            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
     483        else:
     484            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
     485    elif SGData['SGLaue'] in ['mmm',]:
     486        return [pfx+'A0',pfx+'A1',pfx+'A2']
     487    elif SGData['SGLaue'] in ['4/m','4/mmm']:
     488        return [pfx+'A0',pfx+'A2']
     489    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
     490        return [pfx+'A0',pfx+'A2']
     491    elif SGData['SGLaue'] in ['3R', '3mR']:
     492        return [pfx+'A0',pfx+'A3']                       
     493    elif SGData['SGLaue'] in ['m3m','m3']:
     494        return [pfx+'A0',]
     495                   
     496def GetPhaseData(PhaseData,Print=True):
     497           
     498    def PrintFFtable(FFtable):
     499        print '\n X-ray scattering factors:'
     500        print '   Symbol     fa                                      fb                                      fc'
     501        print 99*'-'
     502        for Ename in FFtable:
     503            ffdata = FFtable[Ename]
     504            fa = ffdata['fa']
     505            fb = ffdata['fb']
     506            print ' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' %  \
     507                (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
     508               
     509    def PrintBLtable(BLtable):
     510        print '\n Neutron scattering factors:'
     511        print '   Symbol   isotope       mass       b       resonant terms'
     512        print 99*'-'
     513        for Ename in BLtable:
     514            bldata = BLtable[Ename]
     515            isotope = bldata[0]
     516            mass = bldata[1][0]
     517            blen = bldata[1][1]
     518            bres = []
     519            if len(bldata[1]) > 2:
     520                bres = bldata[1][2:]
     521            line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
     522            for item in bres:
     523                line += '%10.5g'%(item)
     524            print line
     525               
     526    def PrintAtoms(General,Atoms):
     527        print '\n Atoms:'
     528        line = '   name    type  refine?   x         y         z    '+ \
     529            '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
     530        if General['Type'] == 'magnetic':
     531            line += '   Mx     My     Mz'
     532        elif General['Type'] == 'macromolecular':
     533            line = ' res no  residue  chain '+line
     534        print line
     535        if General['Type'] == 'nuclear':
     536            print 135*'-'
     537            for i,at in enumerate(Atoms):
     538                line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \
     539                    '%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9])
     540                if at[9] == 'I':
     541                    line += '%8.4f'%(at[10])+48*' '
     542                else:
     543                    line += 8*' '
     544                    for j in range(6):
     545                        line += '%8.4f'%(at[11+j])
     546                print line
     547       
     548    def PrintTexture(textureData):
     549        topstr = '\n Spherical harmonics texture: Order:' + \
     550            str(textureData['Order'])
     551        if textureData['Order']:
     552            print topstr+' Refine? '+str(textureData['SH Coeff'][0])
     553        else:
     554            print topstr
     555            return
     556        names = ['omega','chi','phi']
     557        line = '\n'
     558        for name in names:
     559            line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
     560        print line
     561        print '\n Texture coefficients:'
     562        ptlbls = ' names :'
     563        ptstr =  ' values:'
     564        SHcoeff = textureData['SH Coeff'][1]
     565        for item in SHcoeff:
     566            ptlbls += '%12s'%(item)
     567            ptstr += '%12.4f'%(SHcoeff[item])
     568        print ptlbls
     569        print ptstr   
     570       
     571    if Print: print ' Phases:'
     572    phaseVary = []
     573    phaseDict = {}
     574    phaseConstr = {}
     575    pawleyLookup = {}
     576    FFtables = {}                   #scattering factors - xrays
     577    BLtables = {}                   # neutrons
     578    Natoms = {}
     579    AtMults = {}
     580    AtIA = {}
     581    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
     582    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
     583    for name in PhaseData:
     584        General = PhaseData[name]['General']
     585        pId = PhaseData[name]['pId']
     586        pfx = str(pId)+'::'
     587        FFtable = GetFFtable(General)
     588        BLtable = GetBLtable(General)
     589        FFtables.update(FFtable)
     590        BLtables.update(BLtable)
     591        Atoms = PhaseData[name]['Atoms']
     592        try:
     593            PawleyRef = PhaseData[name]['Pawley ref']
     594        except KeyError:
     595            PawleyRef = []
     596        SGData = General['SGData']
     597        SGtext = G2spc.SGPrint(SGData)
     598        cell = General['Cell']
     599        A = G2lat.cell2A(cell[1:7])
     600        phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5]})
     601        if cell[0]:
     602            phaseVary += cellVary(pfx,SGData)
     603        Natoms[pfx] = 0
     604        if Atoms:
     605            if General['Type'] == 'nuclear':
     606                Natoms[pfx] = len(Atoms)
     607                for i,at in enumerate(Atoms):
     608                    phaseDict.update({pfx+'Atype:'+str(i):at[1],pfx+'Afrac:'+str(i):at[6],pfx+'Amul:'+str(i):at[8],
     609                        pfx+'Ax:'+str(i):at[3],pfx+'Ay:'+str(i):at[4],pfx+'Az:'+str(i):at[5],
     610                        pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0.,         #refined shifts for x,y,z
     611                        pfx+'AI/A:'+str(i):at[9],})
     612                    if at[9] == 'I':
     613                        phaseDict[pfx+'AUiso:'+str(i)] = at[10]
     614                    else:
     615                        phaseDict.update({pfx+'AU11:'+str(i):at[11],pfx+'AU22:'+str(i):at[12],pfx+'AU33:'+str(i):at[13],
     616                            pfx+'AU12:'+str(i):at[14],pfx+'AU13:'+str(i):at[15],pfx+'AU23:'+str(i):at[16]})
     617                    if 'F' in at[2]:
     618                        phaseVary.append(pfx+'Afrac:'+str(i))
     619                    if 'X' in at[2]:
     620                        xId,xCoef = G2spc.GetCSxinel(at[7])
     621                        delnames = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
     622                        for j in range(3):
     623                            if xId[j] > 0:                               
     624                                phaseVary.append(delnames[j])
     625                                for k in range(j):
     626                                    if xId[j] == xId[k]:
     627                                        G2mv.StoreEquivalence(delnames[k],((delnames[j],xCoef[j]),))
     628                    if 'U' in at[2]:
     629                        if at[9] == 'I':
     630                            phaseVary.append(pfx+'AUiso:'+str(i))
     631                        else:
     632                            uId,uCoef = G2spc.GetCSuinel(at[7])[:2]
     633                            names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
     634                                pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
     635                            for j in range(6):
     636                                if uId[j] > 0:                               
     637                                    phaseVary.append(names[j])
     638                                    for k in range(j):
     639                                        if uId[j] == uId[k]:
     640                                            G2mv.StoreEquivalence(names[k],((names[j],uCoef[j]),))
     641#            elif General['Type'] == 'magnetic':
     642#            elif General['Type'] == 'macromolecular':
     643
     644               
     645            textureData = General['SH Texture']
     646            if textureData['Order']:
     647                phaseDict[pfx+'SHorder'] = textureData['Order']
     648                phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
     649                for name in ['omega','chi','phi']:
     650                    phaseDict[pfx+'SH '+name] = textureData['Sample '+name][1]
     651                    if textureData['Sample '+name][0]:
     652                        phaseVary.append(pfx+'SH '+name)
     653                for name in textureData['SH Coeff'][1]:
     654                    phaseDict[pfx+name] = textureData['SH Coeff'][1][name]
     655                    if textureData['SH Coeff'][0]:
     656                        phaseVary.append(pfx+name)
     657               
     658            if Print:
     659                print '\n Phase name: ',General['Name']
     660                print 135*'-'
     661                PrintFFtable(FFtable)
     662                PrintBLtable(BLtable)
     663                print ''
     664                for line in SGtext: print line
     665                PrintAtoms(General,Atoms)
     666                print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
     667                    ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
     668                    '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
     669                PrintTexture(textureData)
     670                   
     671        elif PawleyRef:
     672            pawleyVary = []
     673            for i,refl in enumerate(PawleyRef):
     674                phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
     675                pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
     676                if refl[5]:
     677                    pawleyVary.append(pfx+'PWLref:'+str(i))
     678            GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary)      #does G2mv.StoreEquivalence
     679            phaseVary += pawleyVary
     680               
     681    return Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
     682   
     683def cellFill(pfx,SGData,parmDict,sigDict):
     684    if SGData['SGLaue'] in ['-1',]:
     685        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
     686            parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
     687        sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
     688            sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
     689    elif SGData['SGLaue'] in ['2/m',]:
     690        if SGData['SGUniq'] == 'a':
     691            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
     692                parmDict[pfx+'A3'],0,0]
     693            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
     694                sigDict[pfx+'A3'],0,0]
     695        elif SGData['SGUniq'] == 'b':
     696            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
     697                0,parmDict[pfx+'A4'],0]
     698            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
     699                0,sigDict[pfx+'A4'],0]
     700        else:
     701            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
     702                0,0,parmDict[pfx+'A5']]
     703            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
     704                0,0,sigDict[pfx+'A5']]
     705    elif SGData['SGLaue'] in ['mmm',]:
     706        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
     707        sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
     708    elif SGData['SGLaue'] in ['4/m','4/mmm']:
     709        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
     710        sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
     711    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
     712        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
     713            parmDict[pfx+'A0'],0,0]
     714        sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
     715    elif SGData['SGLaue'] in ['3R', '3mR']:
     716        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
     717            parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
     718        sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
     719    elif SGData['SGLaue'] in ['m3m','m3']:
     720        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
     721        sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
     722    return A,sigA
     723       
     724def getCellEsd(pfx,SGData,A,covData):
     725    dpr = 180./np.pi
     726    rVsq = G2lat.calc_rVsq(A)
     727    G,g = G2lat.A2Gmat(A)       #get recip. & real metric tensors
     728    cell = np.array(G2lat.Gmat2cell(g))   #real cell
     729    cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
     730    scos = cosd(cellst[3:6])
     731    ssin = sind(cellst[3:6])
     732    scot = scos/ssin
     733    rcos = cosd(cell[3:6])
     734    rsin = sind(cell[3:6])
     735    rcot = rcos/rsin
     736    RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
     737    varyList = covData['varyList']
     738    covMatrix = covData['covMatrix']
     739    vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
     740    Ax = np.array(A)
     741    Ax[3:] /= 2.
     742    drVdA = np.array([Ax[1]*Ax[2]-Ax[5]**2,Ax[0]*Ax[2]-Ax[4]**2,Ax[0]*Ax[1]-Ax[3]**2,
     743        Ax[4]*Ax[5]-Ax[2]*Ax[3],Ax[3]*Ax[5]-Ax[1]*Ax[4],Ax[3]*Ax[4]-Ax[0]*Ax[5]])
     744    srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
     745    Vol = 1/np.sqrt(rVsq)
     746    sigVol = Vol**3*np.sqrt(srcvlsq)/2.
     747    R123 = Ax[0]*Ax[1]*Ax[2]
     748    dsasdg = np.zeros((3,6))
     749    dadg = np.zeros((6,6))
     750    for i0 in range(3):         #0  1   2
     751        i1 = (i0+1)%3           #1  2   0
     752        i2 = (i1+1)%3           #2  0   1
     753        i3 = 5-i2               #3  5   4
     754        i4 = 5-i1               #4  3   5
     755        i5 = 5-i0               #5  4   3
     756        dsasdg[i0][i1] = 0.5*scot[i0]*scos[i0]/Ax[i1]
     757        dsasdg[i0][i2] = 0.5*scot[i0]*scos[i0]/Ax[i2]
     758        dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
     759        denmsq = Ax[i0]*(R123-Ax[i1]*Ax[i4]**2-Ax[i2]*Ax[i3]**2+(Ax[i4]*Ax[i3])**2)
     760        denom = np.sqrt(denmsq)
     761        dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq*(R123-0.5*Ax[i1]*Ax[i4]**2-0.5*Ax[i2]*Ax[i3]**2)
     762        dadg[i5][i1] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i2]-Ax[i0]*Ax[i4]**2)
     763        dadg[i5][i2] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i1]-Ax[i0]*Ax[i3]**2)
     764        dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i2]*Ax[i3]-Ax[i3]*Ax[i4]**2)
     765        dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i1]*Ax[i4]-Ax[i3]**2*Ax[i4])
     766        dadg[i5][i5] = -Ax[i0]/denom
     767    for i0 in range(3):
     768        i1 = (i0+1)%3
     769        i2 = (i1+1)%3
     770        i3 = 5-i2
     771        for ij in range(6):
     772            dadg[i0][ij] = cell[i0]*(rcot[i2]*dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
     773            if ij == i0:
     774                dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
     775            dadg[i3][ij] = -dadg[i3][ij]*rsin[2-i0]*dpr
     776    sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
     777    var = np.diag(sigMat)
     778    CS = np.where(var>0.,np.sqrt(var),0.)
     779    cellSig = [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]  #exchange sig(alp) & sig(gam) to get in right order
     780    return cellSig           
     781   
     782def SetPhaseData(parmDict,sigDict,Phases,covData):
     783   
     784    def PrintAtomsAndSig(General,Atoms,atomsSig):
     785        print '\n Atoms:'
     786        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
     787        if General['Type'] == 'magnetic':
     788            line += '   Mx     My     Mz'
     789        elif General['Type'] == 'macromolecular':
     790            line = ' res no  residue  chain '+line
     791        print line
     792        if General['Type'] == 'nuclear':
     793            print 135*'-'
     794            fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
     795                11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
     796            noFXsig = {3:[10*' ','%10s'],4:[10*' ','%10s'],5:[10*' ','%10s'],6:[8*' ','%8s']}
     797            for i,at in enumerate(Atoms):
     798                name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
     799                valstr = ' values:'
     800                sigstr = ' sig   :'
     801                for ind in [3,4,5,6]:
     802                    sigind = str(i)+':'+str(ind)
     803                    valstr += fmt[ind]%(at[ind])                   
     804                    if sigind in atomsSig:
     805                        sigstr += fmt[ind]%(atomsSig[sigind])
     806                    else:
     807                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
     808                if at[9] == 'I':
     809                    valstr += fmt[10]%(at[10])
     810                    if str(i)+':10' in atomsSig:
     811                        sigstr += fmt[10]%(atomsSig[str(i)+':10'])
     812                    else:
     813                        sigstr += 8*' '
     814                else:
     815                    valstr += 8*' '
     816                    sigstr += 8*' '
     817                    for ind in [11,12,13,14,15,16]:
     818                        sigind = str(i)+':'+str(ind)
     819                        valstr += fmt[ind]%(at[ind])
     820                        if sigind in atomsSig:                       
     821                            sigstr += fmt[ind]%(atomsSig[sigind])
     822                        else:
     823                            sigstr += 8*' '
     824                print name
     825                print valstr
     826                print sigstr
     827               
     828    def PrintSHtextureAndSig(textureData,SHtextureSig):
     829        print '\n Spherical harmonics texture: Order:' + str(textureData['Order'])
     830        names = ['omega','chi','phi']
     831        namstr = '  names :'
     832        ptstr =  '  values:'
     833        sigstr = '  esds  :'
     834        for name in names:
     835            namstr += '%12s'%(name)
     836            ptstr += '%12.3f'%(textureData['Sample '+name][1])
     837            if 'Sample '+name in SHtextureSig:
     838                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
    489839            else:
    490 
    491                 last = max(last,int(name[-2].strip('bak')))
    492 
    493     GPXback = ospath.join(GPXpath,GPXname.rstrip('.'.join(name[-2:]))+'.bak'+str(last)+'.gpx')
    494 
    495     dfu.copy_file(GPXfile,GPXback)
    496 
    497     return GPXback
    498 
    499        
    500 
    501 def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,CovData,makeBack=True):
    502 
    503     ''' Updates gpxfile from all histograms that are found in any phase
    504 
    505     and any phase that used a histogram
    506 
     840                sigstr += 12*' '
     841        print namstr
     842        print ptstr
     843        print sigstr
     844        print '\n Texture coefficients:'
     845        namstr = '  names :'
     846        ptstr =  '  values:'
     847        sigstr = '  esds  :'
     848        SHcoeff = textureData['SH Coeff'][1]
     849        for name in SHcoeff:
     850            namstr += '%12s'%(name)
     851            ptstr += '%12.3f'%(SHcoeff[name])
     852            if name in SHtextureSig:
     853                sigstr += '%12.3f'%(SHtextureSig[name])
     854            else:
     855                sigstr += 12*' '
     856        print namstr
     857        print ptstr
     858        print sigstr
     859       
     860           
     861    print '\n Phases:'
     862    for phase in Phases:
     863        print ' Result for phase: ',phase
     864        Phase = Phases[phase]
     865        General = Phase['General']
     866        SGData = General['SGData']
     867        Atoms = Phase['Atoms']
     868        cell = General['Cell']
     869        pId = Phase['pId']
     870        pfx = str(pId)+'::'
     871        if cell[0]:
     872            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
     873            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
     874            print ' Reciprocal metric tensor: '
     875            ptfmt = "%15.9f"
     876            names = ['A11','A22','A33','A12','A13','A23']
     877            namstr = '  names :'
     878            ptstr =  '  values:'
     879            sigstr = '  esds  :'
     880            for name,a,siga in zip(names,A,sigA):
     881                namstr += '%15s'%(name)
     882                ptstr += ptfmt%(a)
     883                if siga:
     884                    sigstr += ptfmt%(siga)
     885                else:
     886                    sigstr += 15*' '
     887            print namstr
     888            print ptstr
     889            print sigstr
     890            cell[1:7] = G2lat.A2cell(A)
     891            cell[7] = G2lat.calc_V(A)
     892            print ' New unit cell:'
     893            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
     894            names = ['a','b','c','alpha','beta','gamma','Volume']
     895            namstr = '  names :'
     896            ptstr =  '  values:'
     897            sigstr = '  esds  :'
     898            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
     899                namstr += '%12s'%(name)
     900                ptstr += fmt%(a)
     901                if siga:
     902                    sigstr += fmt%(siga)
     903                else:
     904                    sigstr += 12*' '
     905            print namstr
     906            print ptstr
     907            print sigstr
     908           
     909        if 'Pawley' in Phase['General']['Type']:
     910            pawleyRef = Phase['Pawley ref']
     911            for i,refl in enumerate(pawleyRef):
     912                key = pfx+'PWLref:'+str(i)
     913                refl[6] = abs(parmDict[key])        #suppress negative Fsq
     914                if key in sigDict:
     915                    refl[7] = sigDict[key]
     916                else:
     917                    refl[7] = 0
     918        else:
     919            atomsSig = {}
     920            if General['Type'] == 'nuclear':
     921                for i,at in enumerate(Atoms):
     922                    names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
     923                        10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
     924                        14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
     925                    for ind in [3,4,5,6]:
     926                        at[ind] = parmDict[names[ind]]
     927                        if ind in [3,4,5]:
     928                            name = names[ind].replace('A','dA')
     929                        else:
     930                            name = names[ind]
     931                        if name in sigDict:
     932                            atomsSig[str(i)+':'+str(ind)] = sigDict[name]
     933                    if at[9] == 'I':
     934                        at[10] = parmDict[names[10]]
     935                        if names[10] in sigDict:
     936                            atomsSig[str(i)+':10'] = sigDict[names[10]]
     937                    else:
     938                        for ind in [11,12,13,14,15,16]:
     939                            at[ind] = parmDict[names[ind]]
     940                            if names[ind] in sigDict:
     941                                atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
     942            PrintAtomsAndSig(General,Atoms,atomsSig)
     943       
     944        textureData = General['SH Texture']   
     945        if textureData['Order']:
     946            SHtextureSig = {}
     947            for name in ['omega','chi','phi']:
     948                aname = pfx+'SH '+name
     949                textureData['Sample '+name][1] = parmDict[aname]
     950                if aname in sigDict:
     951                    SHtextureSig['Sample '+name] = sigDict[aname]
     952            for name in textureData['SH Coeff'][1]:
     953                aname = pfx+name
     954                textureData['SH Coeff'][1][name] = parmDict[aname]
     955                if aname in sigDict:
     956                    SHtextureSig[name] = sigDict[aname]
     957            PrintSHtextureAndSig(textureData,SHtextureSig)
     958
     959def GetHistogramPhaseData(Phases,Histograms,Print=True):
     960   
     961    def PrintSize(hapData):
     962        if hapData[0] in ['isotropic','uniaxial']:
     963            line = '\n Size model    : %9s'%(hapData[0])
     964            line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
     965            if hapData[0] == 'uniaxial':
     966                line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
     967            print line
     968        else:
     969            print '\n Size model    : %s'%(hapData[0])
     970            Snames = ['S11','S22','S33','S12','S13','S23']
     971            ptlbls = ' names :'
     972            ptstr =  ' values:'
     973            varstr = ' refine:'
     974            for i,name in enumerate(Snames):
     975                ptlbls += '%12s' % (name)
     976                ptstr += '%12.6f' % (hapData[4][i])
     977                varstr += '%12s' % (str(hapData[5][i]))
     978            print ptlbls
     979            print ptstr
     980            print varstr
     981       
     982    def PrintMuStrain(hapData,SGData):
     983        if hapData[0] in ['isotropic','uniaxial']:
     984            line = '\n Mustrain model: %9s'%(hapData[0])
     985            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
     986            if hapData[0] == 'uniaxial':
     987                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
     988            print line
     989        else:
     990            print '\n Mustrain model: %s'%(hapData[0])
     991            Snames = G2spc.MustrainNames(SGData)
     992            ptlbls = ' names :'
     993            ptstr =  ' values:'
     994            varstr = ' refine:'
     995            for i,name in enumerate(Snames):
     996                ptlbls += '%12s' % (name)
     997                ptstr += '%12.6f' % (hapData[4][i])
     998                varstr += '%12s' % (str(hapData[5][i]))
     999            print ptlbls
     1000            print ptstr
     1001            print varstr
     1002
     1003    def PrintHStrain(hapData,SGData):
     1004        print '\n Hydrostatic/elastic strain: '
     1005        Hsnames = G2spc.HStrainNames(SGData)
     1006        ptlbls = ' names :'
     1007        ptstr =  ' values:'
     1008        varstr = ' refine:'
     1009        for i,name in enumerate(Hsnames):
     1010            ptlbls += '%12s' % (name)
     1011            ptstr += '%12.6f' % (hapData[0][i])
     1012            varstr += '%12s' % (str(hapData[1][i]))
     1013        print ptlbls
     1014        print ptstr
     1015        print varstr
     1016
     1017    def PrintSHPO(hapData):
     1018        print '\n Spherical harmonics preferred orientation: Order:' + \
     1019            str(hapData[4])+' Refine? '+str(hapData[2])
     1020        ptlbls = ' names :'
     1021        ptstr =  ' values:'
     1022        for item in hapData[5]:
     1023            ptlbls += '%12s'%(item)
     1024            ptstr += '%12.3f'%(hapData[5][item])
     1025        print ptlbls
     1026        print ptstr
     1027   
     1028    hapDict = {}
     1029    hapVary = []
     1030    controlDict = {}
     1031    poType = {}
     1032    poAxes = {}
     1033    spAxes = {}
     1034    spType = {}
     1035   
     1036    for phase in Phases:
     1037        HistoPhase = Phases[phase]['Histograms']
     1038        SGData = Phases[phase]['General']['SGData']
     1039        cell = Phases[phase]['General']['Cell'][1:7]
     1040        A = G2lat.cell2A(cell)
     1041        pId = Phases[phase]['pId']
     1042        histoList = HistoPhase.keys()
     1043        histoList.sort()
     1044        for histogram in histoList:
     1045            try:
     1046                Histogram = Histograms[histogram]
     1047            except KeyError:                       
     1048                #skip if histogram not included e.g. in a sequential refinement
     1049                continue
     1050            hapData = HistoPhase[histogram]
     1051            hId = Histogram['hId']
     1052            limits = Histogram['Limits'][1]
     1053            inst = Histogram['Instrument Parameters']
     1054            inst = dict(zip(inst[3],inst[1]))
     1055            Zero = inst['Zero']
     1056            if 'C' in inst['Type']:
     1057                try:
     1058                    wave = inst['Lam']
     1059                except KeyError:
     1060                    wave = inst['Lam1']
     1061                dmin = wave/(2.0*sind(limits[1]/2.0))
     1062            pfx = str(pId)+':'+str(hId)+':'
     1063            for item in ['Scale','Extinction']:
     1064                hapDict[pfx+item] = hapData[item][0]
     1065                if hapData[item][1]:
     1066                    hapVary.append(pfx+item)
     1067            names = G2spc.HStrainNames(SGData)
     1068            for i,name in enumerate(names):
     1069                hapDict[pfx+name] = hapData['HStrain'][0][i]
     1070                if hapData['HStrain'][1][i]:
     1071                    hapVary.append(pfx+name)
     1072            controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
     1073            if hapData['Pref.Ori.'][0] == 'MD':
     1074                hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
     1075                controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
     1076                if hapData['Pref.Ori.'][2]:
     1077                    hapVary.append(pfx+'MD')
     1078            else:                           #'SH' spherical harmonics
     1079                controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
     1080                controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
     1081                for item in hapData['Pref.Ori.'][5]:
     1082                    hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
     1083                    if hapData['Pref.Ori.'][2]:
     1084                        hapVary.append(pfx+item)
     1085            for item in ['Mustrain','Size']:
     1086                controlDict[pfx+item+'Type'] = hapData[item][0]
     1087                if hapData[item][0] in ['isotropic','uniaxial']:
     1088                    hapDict[pfx+item+':i'] = hapData[item][1][0]
     1089                    if hapData[item][2][0]:
     1090                        hapVary.append(pfx+item+':i')
     1091                    if hapData[item][0] == 'uniaxial':
     1092                        controlDict[pfx+item+'Axis'] = hapData[item][3]
     1093                        hapDict[pfx+item+':a'] = hapData[item][1][1]
     1094                        if hapData[item][2][1]:
     1095                            hapVary.append(pfx+item+':a')
     1096                else:       #generalized for mustrain or ellipsoidal for size
     1097                    if item == 'Mustrain':
     1098                        names = G2spc.MustrainNames(SGData)
     1099                        pwrs = []
     1100                        for name in names:
     1101                            h,k,l = name[1:]
     1102                            pwrs.append([int(h),int(k),int(l)])
     1103                        controlDict[pfx+'MuPwrs'] = pwrs
     1104                    for i in range(len(hapData[item][4])):
     1105                        sfx = ':'+str(i)
     1106                        hapDict[pfx+item+sfx] = hapData[item][4][i]
     1107                        if hapData[item][5][i]:
     1108                            hapVary.append(pfx+item+sfx)
     1109                           
     1110            if Print:
     1111                print '\n Phase: ',phase,' in histogram: ',histogram
     1112                print 135*'-'
     1113                print ' Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
     1114                print ' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
     1115                if hapData['Pref.Ori.'][0] == 'MD':
     1116                    Ax = hapData['Pref.Ori.'][3]
     1117                    print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
     1118                        ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
     1119                else: #'SH' for spherical harmonics
     1120                    PrintSHPO(hapData['Pref.Ori.'])
     1121                PrintSize(hapData['Size'])
     1122                PrintMuStrain(hapData['Mustrain'],SGData)
     1123                PrintHStrain(hapData['HStrain'],SGData)
     1124            HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
     1125            refList = []
     1126            for h,k,l,d in HKLd:
     1127                ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
     1128                if ext:
     1129                    continue
     1130                if 'C' in inst['Type']:
     1131                    pos = 2.0*asind(wave/(2.0*d))+Zero
     1132                    if limits[0] < pos < limits[1]:
     1133                        refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0,{}])
     1134                        #last item should contain form factor values by atom type
     1135                else:
     1136                    raise ValueError
     1137            Histogram['Reflection Lists'][phase] = refList
     1138    return hapVary,hapDict,controlDict
     1139   
     1140def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True):
     1141   
     1142    def PrintSizeAndSig(hapData,sizeSig):
     1143        line = '\n Size model:     %9s'%(hapData[0])
     1144        refine = False
     1145        if hapData[0] in ['isotropic','uniaxial']:
     1146            line += ' equatorial:%12.3f'%(hapData[1][0])
     1147            if sizeSig[0][0]:
     1148                line += ', sig: %8.3f'%(sizeSig[0][0])
     1149                refine = True
     1150            if hapData[0] == 'uniaxial':
     1151                line += ' axial:%12.3f'%(hapData[1][1])
     1152                if sizeSig[0][1]:
     1153                    refine = True
     1154                    line += ', sig: %8.3f'%(sizeSig[0][1])
     1155            if refine:
     1156                print line
     1157        else:
     1158            Snames = ['S11','S22','S33','S12','S13','S23']
     1159            ptlbls = ' name  :'
     1160            ptstr =  ' value :'
     1161            sigstr = ' sig   :'
     1162            for i,name in enumerate(Snames):
     1163                ptlbls += '%12s' % (name)
     1164                ptstr += '%12.6f' % (hapData[4][i])
     1165                if sizeSig[1][i]:
     1166                    refine = True
     1167                    sigstr += '%12.6f' % (sizeSig[1][i])
     1168                else:
     1169                    sigstr += 12*' '
     1170            if refine:
     1171                print line
     1172                print ptlbls
     1173                print ptstr
     1174                print sigstr
     1175       
     1176    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
     1177        line = '\n Mustrain model: %9s'%(hapData[0])
     1178        refine = False
     1179        if hapData[0] in ['isotropic','uniaxial']:
     1180            line += ' equatorial:%12.1f'%(hapData[1][0])
     1181            if mustrainSig[0][0]:
     1182                line += ', sig: %8.1f'%(mustrainSig[0][0])
     1183                refine = True
     1184            if hapData[0] == 'uniaxial':
     1185                line += ' axial:%12.1f'%(hapData[1][1])
     1186                if mustrainSig[0][1]:
     1187                     line += ', sig: %8.1f'%(mustrainSig[0][1])
     1188            if refine:
     1189                print line
     1190        else:
     1191            Snames = G2spc.MustrainNames(SGData)
     1192            ptlbls = ' name  :'
     1193            ptstr =  ' value :'
     1194            sigstr = ' sig   :'
     1195            for i,name in enumerate(Snames):
     1196                ptlbls += '%12s' % (name)
     1197                ptstr += '%12.6f' % (hapData[4][i])
     1198                if mustrainSig[1][i]:
     1199                    refine = True
     1200                    sigstr += '%12.6f' % (mustrainSig[1][i])
     1201                else:
     1202                    sigstr += 12*' '
     1203            if refine:
     1204                print line
     1205                print ptlbls
     1206                print ptstr
     1207                print sigstr
     1208           
     1209    def PrintHStrainAndSig(hapData,strainSig,SGData):
     1210        Hsnames = G2spc.HStrainNames(SGData)
     1211        ptlbls = ' name  :'
     1212        ptstr =  ' value :'
     1213        sigstr = ' sig   :'
     1214        refine = False
     1215        for i,name in enumerate(Hsnames):
     1216            ptlbls += '%12s' % (name)
     1217            ptstr += '%12.6g' % (hapData[0][i])
     1218            if name in strainSig:
     1219                refine = True
     1220                sigstr += '%12.6g' % (strainSig[name])
     1221            else:
     1222                sigstr += 12*' '
     1223        if refine:
     1224            print '\n Hydrostatic/elastic strain: '
     1225            print ptlbls
     1226            print ptstr
     1227            print sigstr
     1228       
     1229    def PrintSHPOAndSig(hapData,POsig):
     1230        print '\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
     1231        ptlbls = ' names :'
     1232        ptstr =  ' values:'
     1233        sigstr = ' sig   :'
     1234        for item in hapData[5]:
     1235            ptlbls += '%12s'%(item)
     1236            ptstr += '%12.3f'%(hapData[5][item])
     1237            if item in POsig:
     1238                sigstr += '%12.3f'%(POsig[item])
     1239            else:
     1240                sigstr += 12*' '
     1241        print ptlbls
     1242        print ptstr
     1243        print sigstr
     1244   
     1245    for phase in Phases:
     1246        HistoPhase = Phases[phase]['Histograms']
     1247        SGData = Phases[phase]['General']['SGData']
     1248        pId = Phases[phase]['pId']
     1249        histoList = HistoPhase.keys()
     1250        histoList.sort()
     1251        for histogram in histoList:
     1252            try:
     1253                Histogram = Histograms[histogram]
     1254            except KeyError:                       
     1255                #skip if histogram not included e.g. in a sequential refinement
     1256                continue
     1257            print '\n Phase: ',phase,' in histogram: ',histogram
     1258            print 130*'-'
     1259            hapData = HistoPhase[histogram]
     1260            hId = Histogram['hId']
     1261            pfx = str(pId)+':'+str(hId)+':'
     1262            print ' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
     1263                %(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
     1264           
     1265            PhFrExtPOSig = {}
     1266            for item in ['Scale','Extinction']:
     1267                hapData[item][0] = parmDict[pfx+item]
     1268                if pfx+item in sigDict:
     1269                    PhFrExtPOSig[item] = sigDict[pfx+item]
     1270            if hapData['Pref.Ori.'][0] == 'MD':
     1271                hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
     1272                if pfx+'MD' in sigDict:
     1273                    PhFrExtPOSig['MD'] = sigDict[pfx+'MD']
     1274            else:                           #'SH' spherical harmonics
     1275                for item in hapData['Pref.Ori.'][5]:
     1276                    hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
     1277                    if pfx+item in sigDict:
     1278                        PhFrExtPOSig[item] = sigDict[pfx+item]
     1279            if Print:
     1280                if 'Scale' in PhFrExtPOSig:
     1281                    print ' Phase fraction  : %10.4f, sig %10.4f'%(hapData['Scale'][0],PhFrExtPOSig['Scale'])
     1282                if 'Extinction' in PhFrExtPOSig:
     1283                    print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction'])
     1284                if hapData['Pref.Ori.'][0] == 'MD':
     1285                    if 'MD' in PhFrExtPOSig:
     1286                        print ' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD'])
     1287                else:
     1288                    PrintSHPOAndSig(hapData['Pref.Ori.'],PhFrExtPOSig)
     1289            SizeMuStrSig = {'Mustrain':[[0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
     1290                'Size':[[0,0],[0 for i in range(len(hapData['Size'][4]))]],
     1291                'HStrain':{}}                 
     1292            for item in ['Mustrain','Size']:
     1293                if hapData[item][0] in ['isotropic','uniaxial']:                   
     1294                    hapData[item][1][0] = parmDict[pfx+item+':i']
     1295                    if item == 'Size':
     1296                        hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
     1297                    if pfx+item+':i' in sigDict:
     1298                        SizeMuStrSig[item][0][0] = sigDict[pfx+item+':i']
     1299                    if hapData[item][0] == 'uniaxial':
     1300                        hapData[item][1][1] = parmDict[pfx+item+':a']
     1301                        if item == 'Size':
     1302                            hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
     1303                        if pfx+item+':a' in sigDict:
     1304                            SizeMuStrSig[item][0][1] = sigDict[pfx+item+':a']
     1305                else:       #generalized for mustrain or ellipsoidal for size
     1306                    for i in range(len(hapData[item][4])):
     1307                        sfx = ':'+str(i)
     1308                        hapData[item][4][i] = parmDict[pfx+item+sfx]
     1309                        if pfx+item+sfx in sigDict:
     1310                            SizeMuStrSig[item][1][i] = sigDict[pfx+item+sfx]
     1311            names = G2spc.HStrainNames(SGData)
     1312            for i,name in enumerate(names):
     1313                hapData['HStrain'][0][i] = parmDict[pfx+name]
     1314                if pfx+name in sigDict:
     1315                    SizeMuStrSig['HStrain'][name] = sigDict[pfx+name]
     1316            if Print:
     1317                PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size'])
     1318                PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData)
     1319                PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig['HStrain'],SGData)
     1320   
     1321def GetHistogramData(Histograms,Print=True):
     1322   
     1323    def GetBackgroundParms(hId,Background):
     1324        Back = Background[0]
     1325        Debye = Background[1]
     1326        bakType,bakFlag = Back[:2]
     1327        backVals = Back[3:]
     1328        backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
     1329        backDict = dict(zip(backNames,backVals))
     1330        backVary = []
     1331        if bakFlag:
     1332            backVary = backNames
     1333        backDict[':'+str(hId)+':nDebye'] = Debye['nDebye']
     1334        debyeDict = {}
     1335        debyeList = []
     1336        for i in range(Debye['nDebye']):
     1337            debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
     1338            debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
     1339            debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
     1340        debyeVary = []
     1341        for item in debyeList:
     1342            if item[1]:
     1343                debyeVary.append(item[0])
     1344        backDict.update(debyeDict)
     1345        backVary += debyeVary   
     1346        return bakType,backDict,backVary           
     1347       
     1348    def GetInstParms(hId,Inst):
     1349        insVals,insFlags,insNames = Inst[1:4]
     1350        dataType = insVals[0]
     1351        instDict = {}
     1352        insVary = []
     1353        pfx = ':'+str(hId)+':'
     1354        for i,flag in enumerate(insFlags):
     1355            insName = pfx+insNames[i]
     1356            instDict[insName] = insVals[i]
     1357            if flag:
     1358                insVary.append(insName)
     1359        instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001)
     1360        instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001)
     1361        instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
     1362        return dataType,instDict,insVary
     1363       
     1364    def GetSampleParms(hId,Sample):
     1365        sampVary = []
     1366        hfx = ':'+str(hId)+':'       
     1367        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
     1368            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
     1369        Type = Sample['Type']
     1370        if 'Bragg' in Type:             #Bragg-Brentano
     1371            for item in ['Scale','Shift','Transparency']:       #surface roughness?, diffuse scattering?
     1372                sampDict[hfx+item] = Sample[item][0]
     1373                if Sample[item][1]:
     1374                    sampVary.append(hfx+item)
     1375        elif 'Debye' in Type:        #Debye-Scherrer
     1376            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
     1377                sampDict[hfx+item] = Sample[item][0]
     1378                if Sample[item][1]:
     1379                    sampVary.append(hfx+item)
     1380        return Type,sampDict,sampVary
     1381       
     1382    def PrintBackground(Background):
     1383        Back = Background[0]
     1384        Debye = Background[1]
     1385        print '\n Background function: ',Back[0],' Refine?',bool(Back[1])
     1386        line = ' Coefficients: '
     1387        for i,back in enumerate(Back[3:]):
     1388            line += '%10.3f'%(back)
     1389            if i and not i%10:
     1390                line += '\n'+15*' '
     1391        print line
     1392        if Debye['nDebye']:
     1393            print '\n Debye diffuse scattering coefficients'
     1394            parms = ['DebyeA','DebyeR','DebyeU']
     1395            line = ' names :'
     1396            for parm in parms:
     1397                line += '%16s'%(parm)
     1398            print line
     1399            for j,term in enumerate(Debye['debyeTerms']):
     1400                line = ' term'+'%2d'%(j)+':'
     1401                for i in range(3):
     1402                    line += '%10.4g %5s'%(term[2*i],bool(term[2*i+1]))                   
     1403                print line
     1404       
     1405    def PrintInstParms(Inst):
     1406        print '\n Instrument Parameters:'
     1407        ptlbls = ' name  :'
     1408        ptstr =  ' value :'
     1409        varstr = ' refine:'
     1410        instNames = Inst[3][1:]
     1411        for i,name in enumerate(instNames):
     1412            ptlbls += '%12s' % (name)
     1413            ptstr += '%12.6f' % (Inst[1][i+1])
     1414            if name in ['Lam1','Lam2','Azimuth']:
     1415                varstr += 12*' '
     1416            else:
     1417                varstr += '%12s' % (str(bool(Inst[2][i+1])))
     1418        print ptlbls
     1419        print ptstr
     1420        print varstr
     1421       
     1422    def PrintSampleParms(Sample):
     1423        print '\n Sample Parameters:'
     1424        print ' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
     1425            (Sample['Omega'],Sample['Chi'],Sample['Phi'])
     1426        ptlbls = ' name  :'
     1427        ptstr =  ' value :'
     1428        varstr = ' refine:'
     1429        if 'Bragg' in Sample['Type']:
     1430            for item in ['Scale','Shift','Transparency']:
     1431                ptlbls += '%14s'%(item)
     1432                ptstr += '%14.4f'%(Sample[item][0])
     1433                varstr += '%14s'%(str(bool(Sample[item][1])))
     1434           
     1435        elif 'Debye' in Type:        #Debye-Scherrer
     1436            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
     1437                ptlbls += '%14s'%(item)
     1438                ptstr += '%14.4f'%(Sample[item][0])
     1439                varstr += '%14s'%(str(bool(Sample[item][1])))
     1440
     1441        print ptlbls
     1442        print ptstr
     1443        print varstr
     1444       
     1445
     1446    histDict = {}
     1447    histVary = []
     1448    controlDict = {}
     1449    histoList = Histograms.keys()
     1450    histoList.sort()
     1451    for histogram in histoList:
     1452        Histogram = Histograms[histogram]
     1453        hId = Histogram['hId']
     1454        pfx = ':'+str(hId)+':'
     1455        controlDict[pfx+'Limits'] = Histogram['Limits'][1]
     1456       
     1457        Background = Histogram['Background']
     1458        Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
     1459        controlDict[pfx+'bakType'] = Type
     1460        histDict.update(bakDict)
     1461        histVary += bakVary
     1462       
     1463        Inst = Histogram['Instrument Parameters']
     1464        Type,instDict,insVary = GetInstParms(hId,Inst)
     1465        controlDict[pfx+'histType'] = Type
     1466        if pfx+'Lam1' in instDict:
     1467            controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
     1468        else:
     1469            controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
     1470        histDict.update(instDict)
     1471        histVary += insVary
     1472       
     1473        Sample = Histogram['Sample Parameters']
     1474        Type,sampDict,sampVary = GetSampleParms(hId,Sample)
     1475        controlDict[pfx+'instType'] = Type
     1476        histDict.update(sampDict)
     1477        histVary += sampVary
     1478
     1479        if Print:
     1480            print '\n Histogram: ',histogram,' histogram Id: ',hId
     1481            print 135*'-'
     1482            Units = {'C':' deg','T':' msec'}
     1483            units = Units[controlDict[pfx+'histType'][2]]
     1484            Limits = controlDict[pfx+'Limits']
     1485            print ' Instrument type: ',Sample['Type']
     1486            print ' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)     
     1487            PrintSampleParms(Sample)
     1488            PrintInstParms(Inst)
     1489            PrintBackground(Background)
     1490       
     1491    return histVary,histDict,controlDict
     1492   
     1493def SetHistogramData(parmDict,sigDict,Histograms,Print=True):
     1494   
     1495    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
     1496        Back = Background[0]
     1497        Debye = Background[1]
     1498        lenBack = len(Back[3:])
     1499        backSig = [0 for i in range(lenBack+3*Debye['nDebye'])]
     1500        for i in range(lenBack):
     1501            Back[3+i] = parmDict[pfx+'Back:'+str(i)]
     1502            if pfx+'Back:'+str(i) in sigDict:
     1503                backSig[i] = sigDict[pfx+'Back:'+str(i)]
     1504        if Debye['nDebye']:
     1505            for i in range(Debye['nDebye']):
     1506                names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
     1507                for j,name in enumerate(names):
     1508                    Debye['debyeTerms'][i][2*j] = parmDict[name]
     1509                    if name in sigDict:
     1510                        backSig[lenBack+3*i+j] = sigDict[name]           
     1511        return backSig
     1512       
     1513    def SetInstParms(pfx,Inst,parmDict,sigDict):
     1514        insVals,insFlags,insNames = Inst[1:4]
     1515        instSig = [0 for i in range(len(insVals))]
     1516        for i,flag in enumerate(insFlags):
     1517            insName = pfx+insNames[i]
     1518            insVals[i] = parmDict[insName]
     1519            if insName in sigDict:
     1520                instSig[i] = sigDict[insName]
     1521        return instSig
     1522       
     1523    def SetSampleParms(pfx,Sample,parmDict,sigDict):
     1524        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
     1525            sampSig = [0 for i in range(3)]
     1526            for i,item in enumerate(['Scale','Shift','Transparency']):       #surface roughness?, diffuse scattering?
     1527                Sample[item][0] = parmDict[pfx+item]
     1528                if pfx+item in sigDict:
     1529                    sampSig[i] = sigDict[pfx+item]
     1530        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
     1531            sampSig = [0 for i in range(4)]
     1532            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
     1533                Sample[item][0] = parmDict[pfx+item]
     1534                if pfx+item in sigDict:
     1535                    sampSig[i] = sigDict[pfx+item]
     1536        return sampSig
     1537       
     1538    def PrintBackgroundSig(Background,backSig):
     1539        Back = Background[0]
     1540        Debye = Background[1]
     1541        lenBack = len(Back[3:])
     1542        valstr = ' value : '
     1543        sigstr = ' sig   : '
     1544        refine = False
     1545        for i,back in enumerate(Back[3:]):
     1546            valstr += '%10.4g'%(back)
     1547            if Back[1]:
     1548                refine = True
     1549                sigstr += '%10.4g'%(backSig[i])
     1550            else:
     1551                sigstr += 10*' '
     1552        if refine:
     1553            print '\n Background function: ',Back[0]
     1554            print valstr
     1555            print sigstr
     1556        if Debye['nDebye']:
     1557            ifAny = False
     1558            ptfmt = "%12.5f"
     1559            names =  ' names :'
     1560            ptstr =  ' values:'
     1561            sigstr = ' esds  :'
     1562            for item in sigDict:
     1563                if 'Debye' in item:
     1564                    ifAny = True
     1565                    names += '%12s'%(item)
     1566                    ptstr += ptfmt%(parmDict[item])
     1567                    sigstr += ptfmt%(sigDict[item])
     1568            if ifAny:
     1569                print '\n Debye diffuse scattering coefficients'
     1570                print names
     1571                print ptstr
     1572                print sigstr
     1573       
     1574    def PrintInstParmsSig(Inst,instSig):
     1575        ptlbls = ' names :'
     1576        ptstr =  ' value :'
     1577        sigstr = ' sig   :'
     1578        instNames = Inst[3][1:]
     1579        refine = False
     1580        for i,name in enumerate(instNames):
     1581            ptlbls += '%12s' % (name)
     1582            ptstr += '%12.6f' % (Inst[1][i+1])
     1583            if instSig[i+1]:
     1584                refine = True
     1585                sigstr += '%12.6f' % (instSig[i+1])
     1586            else:
     1587                sigstr += 12*' '
     1588        if refine:
     1589            print '\n Instrument Parameters:'
     1590            print ptlbls
     1591            print ptstr
     1592            print sigstr
     1593       
     1594    def PrintSampleParmsSig(Sample,sampleSig):
     1595        ptlbls = ' names :'
     1596        ptstr =  ' values:'
     1597        sigstr = ' sig   :'
     1598        refine = False
     1599        if 'Bragg' in Sample['Type']:
     1600            for i,item in enumerate(['Scale','Shift','Transparency']):
     1601                ptlbls += '%14s'%(item)
     1602                ptstr += '%14.4f'%(Sample[item][0])
     1603                if sampleSig[i]:
     1604                    refine = True
     1605                    sigstr += '%14.4f'%(sampleSig[i])
     1606                else:
     1607                    sigstr += 14*' '
     1608           
     1609        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
     1610            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
     1611                ptlbls += '%14s'%(item)
     1612                ptstr += '%14.4f'%(Sample[item][0])
     1613                if sampleSig[i]:
     1614                    refine = True
     1615                    sigstr += '%14.4f'%(sampleSig[i])
     1616                else:
     1617                    sigstr += 14*' '
     1618
     1619        if refine:
     1620            print '\n Sample Parameters:'
     1621            print ptlbls
     1622            print ptstr
     1623            print sigstr
     1624       
     1625    histoList = Histograms.keys()
     1626    histoList.sort()
     1627    for histogram in histoList:
     1628        if 'PWDR' in histogram:
     1629            Histogram = Histograms[histogram]
     1630            hId = Histogram['hId']
     1631            pfx = ':'+str(hId)+':'
     1632            Background = Histogram['Background']
     1633            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
     1634           
     1635            Inst = Histogram['Instrument Parameters']
     1636            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
     1637       
     1638            Sample = Histogram['Sample Parameters']
     1639            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
     1640
     1641            print '\n Histogram: ',histogram,' histogram Id: ',hId
     1642            print 135*'-'
     1643            print ' Final refinement wRp = %.2f%% on %d observations in this histogram'%(Histogram['wRp'],Histogram['Nobs'])
     1644            if Print:
     1645                print ' Instrument type: ',Sample['Type']
     1646                PrintSampleParmsSig(Sample,sampSig)
     1647                PrintInstParmsSig(Inst,instSig)
     1648                PrintBackgroundSig(Background,backSig)
     1649
     1650def GetAtomFXU(pfx,calcControls,parmDict):
     1651    Natoms = calcControls['Natoms'][pfx]
     1652    Tdata = Natoms*[' ',]
     1653    Mdata = np.zeros(Natoms)
     1654    IAdata = Natoms*[' ',]
     1655    Fdata = np.zeros(Natoms)
     1656    FFdata = []
     1657    BLdata = []
     1658    Xdata = np.zeros((3,Natoms))
     1659    dXdata = np.zeros((3,Natoms))
     1660    Uisodata = np.zeros(Natoms)
     1661    Uijdata = np.zeros((6,Natoms))
     1662    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
     1663        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
     1664        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
     1665        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
     1666        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
     1667    for iatm in range(Natoms):
     1668        for key in keys:
     1669            parm = pfx+key+str(iatm)
     1670            if parm in parmDict:
     1671                keys[key][iatm] = parmDict[parm]
     1672    return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
     1673   
     1674def getFFvalues(FFtables,SQ):
     1675    FFvals = {}
     1676    for El in FFtables:
     1677        FFvals[El] = G2el.ScatFac(FFtables[El],SQ)[0]
     1678    return FFvals
     1679   
     1680def getBLvalues(BLtables):
     1681    BLvals = {}
     1682    for El in BLtables:
     1683        BLvals[El] = BLtables[El][1]
     1684    return BLvals
     1685       
     1686def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
     1687    ''' Compute structure factors for all h,k,l for phase
    5071688    input:
    508 
    509         GPXfile = .gpx full file name
    510 
    511         Histograms = dictionary of histograms as {name:data,...}
    512 
    513         Phases = dictionary of phases that use histograms
    514 
    515         CovData = dictionary of refined variables, varyList, & covariance matrix
    516 
    517         makeBack = True if new backup of .gpx file is to be made; else use the last one made
    518 
     1689        refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
     1690        G:      reciprocal metric tensor
     1691        pfx:    phase id string
     1692        SGData: space group info. dictionary output from SpcGroup
     1693        calcControls:
     1694        ParmDict:
     1695    puts result F^2 in each ref[8] in refList
     1696    '''       
     1697    twopi = 2.0*np.pi
     1698    twopisq = 2.0*np.pi**2
     1699    ast = np.sqrt(np.diag(G))
     1700    Mast = twopisq*np.multiply.outer(ast,ast)
     1701    FFtables = calcControls['FFtables']
     1702    BLtables = calcControls['BLtables']
     1703    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     1704    FF = np.zeros(len(Tdata))
     1705    if 'N' in parmDict[hfx+'Type']:
     1706        FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
     1707    else:
     1708        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
     1709        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
     1710    maxPos = len(SGData['SGOps'])
     1711    Uij = np.array(G2lat.U6toUij(Uijdata))
     1712    bij = Mast*Uij.T
     1713    for refl in refList:
     1714        fbs = np.array([0,0])
     1715        H = refl[:3]
     1716        SQ = 1./(2.*refl[4])**2
     1717        if not len(refl[-1]):                #no form factors
     1718            if 'N' in parmDict[hfx+'Type']:
     1719                refl[-1] = getBLvalues(BLtables)
     1720            else:       #'X'
     1721                refl[-1] = getFFvalues(FFtables,SQ)
     1722        for i,El in enumerate(Tdata):           
     1723            FF[i] = refl[-1][El]           
     1724        SQfactor = 4.0*SQ*twopisq
     1725        Uniq = refl[11]
     1726        phi = refl[12]
     1727        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
     1728        sinp = np.sin(phase)
     1729        cosp = np.cos(phase)
     1730        occ = Mdata*Fdata/len(Uniq)
     1731        biso = -SQfactor*Uisodata
     1732        Tiso = np.where(biso<1.,np.exp(biso),1.0)
     1733        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
     1734        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
     1735        Tcorr = Tiso*Tuij
     1736        fa = np.array([(FF+FP)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr])
     1737        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
     1738        if not SGData['SGInv']:
     1739            fb = np.array([(FF+FP)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
     1740            fbs = np.sum(np.sum(fb,axis=1),axis=1)
     1741        fasq = fas**2
     1742        fbsq = fbs**2        #imaginary
     1743        refl[9] = np.sum(fasq)+np.sum(fbsq)
     1744        refl[10] = atan2d(fbs[0],fas[0])
     1745    return refList
     1746   
     1747def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
     1748    twopi = 2.0*np.pi
     1749    twopisq = 2.0*np.pi**2
     1750    ast = np.sqrt(np.diag(G))
     1751    Mast = twopisq*np.multiply.outer(ast,ast)
     1752    FFtables = calcControls['FFtables']
     1753    BLtables = calcControls['BLtables']
     1754    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     1755    FF = np.zeros(len(Tdata))
     1756    if 'N' in parmDict[hfx+'Type']:
     1757        FP = 0.
     1758        FPP = 0.
     1759    else:
     1760        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
     1761        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
     1762    maxPos = len(SGData['SGOps'])       
     1763    Uij = np.array(G2lat.U6toUij(Uijdata))
     1764    bij = Mast*Uij.T
     1765    dFdvDict = {}
     1766    dFdfr = np.zeros((len(refList),len(Mdata)))
     1767    dFdx = np.zeros((len(refList),len(Mdata),3))
     1768    dFdui = np.zeros((len(refList),len(Mdata)))
     1769    dFdua = np.zeros((len(refList),len(Mdata),6))
     1770    for iref,refl in enumerate(refList):
     1771        H = np.array(refl[:3])
     1772        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
     1773        for i,El in enumerate(Tdata):           
     1774            FF[i] = refl[-1][El]           
     1775        SQfactor = 8.0*SQ*np.pi**2
     1776        Uniq = refl[11]
     1777        phi = refl[12]
     1778        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
     1779        sinp = np.sin(phase)
     1780        cosp = np.cos(phase)
     1781        occ = Mdata*Fdata/len(Uniq)
     1782        biso = -SQfactor*Uisodata
     1783        Tiso = np.where(biso<1.,np.exp(biso),1.0)
     1784#        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
     1785        HbH = -np.inner(H,np.inner(bij,H))
     1786        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
     1787        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
     1788        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
     1789        Tcorr = Tiso*Tuij
     1790        fot = (FF+FP)*occ*Tcorr
     1791        fotp = FPP*occ*Tcorr
     1792        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
     1793        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
     1794       
     1795        fas = np.sum(np.sum(fa,axis=1),axis=1)
     1796        fbs = np.sum(np.sum(fb,axis=1),axis=1)
     1797        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
     1798        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
     1799        #sum below is over Uniq
     1800        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
     1801        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
     1802        dfadui = np.sum(-SQfactor*fa,axis=2)
     1803        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
     1804        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for     
     1805        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
     1806        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
     1807        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
     1808        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
     1809        if not SGData['SGInv']:
     1810            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #problem here if occ=0 for some atom
     1811            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)         
     1812            dfbdui = np.sum(-SQfactor*fb,axis=2)
     1813            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
     1814            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
     1815            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
     1816            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
     1817            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
     1818        #loop over atoms - each dict entry is list of derivatives for all the reflections
     1819    for i in range(len(Mdata)):     
     1820        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
     1821        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
     1822        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
     1823        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
     1824        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
     1825        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
     1826        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
     1827        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
     1828        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
     1829        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
     1830        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
     1831    return dFdvDict
     1832       
     1833def Dict2Values(parmdict, varylist):
     1834    '''Use before call to leastsq to setup list of values for the parameters
     1835    in parmdict, as selected by key in varylist'''
     1836    return [parmdict[key] for key in varylist]
     1837   
     1838def Values2Dict(parmdict, varylist, values):
     1839    ''' Use after call to leastsq to update the parameter dictionary with
     1840    values corresponding to keys in varylist'''
     1841    parmdict.update(zip(varylist,values))
     1842   
     1843def GetNewCellParms(parmDict,varyList):
     1844    newCellDict = {}
     1845    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],['A'+str(i) for i in range(6)]))
     1846    for item in varyList:
     1847        keys = item.split(':')
     1848        if keys[2] in Ddict:
     1849            key = keys[0]+'::'+Ddict[keys[2]]
     1850            parm = keys[0]+'::'+keys[2]
     1851            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
     1852    return newCellDict
     1853   
     1854def ApplyXYZshifts(parmDict,varyList):
     1855    ''' takes atom x,y,z shift and applies it to corresponding atom x,y,z value
     1856        input:
     1857            parmDict - parameter dictionary
     1858            varyList - list of variables
     1859        returns:
     1860            newAtomDict - dictitemionary of new atomic coordinate names & values;
     1861                key is parameter shift name
    5191862    '''
    520 
    521                        
    522 
    523     GPXback = GPXBackup(GPXfile,makeBack)
    524 
    525     print '\n',135*'-'
    526 
    527     print 'Read from file:',GPXback
    528 
    529     print 'Save to file  :',GPXfile
    530 
    531     infile = open(GPXback,'rb')
    532 
    533     outfile = open(GPXfile,'wb')
     1863    newAtomDict = {}
     1864    for item in parmDict:
     1865        if 'dA' in item:
     1866            parm = ''.join(item.split('d'))
     1867            parmDict[parm] += parmDict[item]
     1868            newAtomDict[item] = [parm,parmDict[parm]]
     1869    return newAtomDict
     1870   
     1871def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
     1872    IFCoup = 'Bragg' in calcControls[hfx+'instType']
     1873    odfCor = 1.0
     1874    H = refl[:3]
     1875    cell = G2lat.Gmat2cell(g)
     1876    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
     1877    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
     1878    phi,beta = G2lat.CrsAng(H,cell,SGData)
     1879    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
     1880    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
     1881    for item in SHnames:
     1882        L,M,N = eval(item.strip('C'))
     1883        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
     1884        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
     1885        Lnorm = G2lat.Lnorm(L)
     1886        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
     1887    return odfCor
     1888   
     1889def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
     1890    FORPI = 12.5663706143592
     1891    IFCoup = 'Bragg' in calcControls[hfx+'instType']
     1892    odfCor = 1.0
     1893    dFdODF = {}
     1894    dFdSA = [0,0,0]
     1895    H = refl[:3]
     1896    cell = G2lat.Gmat2cell(g)
     1897    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
     1898    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
     1899    phi,beta = G2lat.CrsAng(H,cell,SGData)
     1900    psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
     1901    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
     1902    for item in SHnames:
     1903        L,M,N = eval(item.strip('C'))
     1904        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
     1905        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
     1906        Lnorm = G2lat.Lnorm(L)
     1907        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
     1908        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
     1909        for i in range(3):
     1910            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
     1911    return odfCor,dFdODF,dFdSA
     1912   
     1913def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
     1914    odfCor = 1.0
     1915    H = refl[:3]
     1916    cell = G2lat.Gmat2cell(g)
     1917    Sangl = [0.,0.,0.]
     1918    if 'Bragg' in calcControls[hfx+'instType']:
     1919        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
     1920        IFCoup = True
     1921    else:
     1922        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
     1923        IFCoup = False
     1924    phi,beta = G2lat.CrsAng(H,cell,SGData)
     1925    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
     1926    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
     1927    for item in SHnames:
     1928        L,N = eval(item.strip('C'))
     1929        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
     1930        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
     1931    return odfCor
     1932   
     1933def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
     1934    FORPI = 12.5663706143592
     1935    odfCor = 1.0
     1936    dFdODF = {}
     1937    H = refl[:3]
     1938    cell = G2lat.Gmat2cell(g)
     1939    Sangl = [0.,0.,0.]
     1940    if 'Bragg' in calcControls[hfx+'instType']:
     1941        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
     1942        IFCoup = True
     1943    else:
     1944        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
     1945        IFCoup = False
     1946    phi,beta = G2lat.CrsAng(H,cell,SGData)
     1947    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
     1948    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
     1949    for item in SHnames:
     1950        L,N = eval(item.strip('C'))
     1951        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
     1952        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
     1953        dFdODF[phfx+item] = Kcsl*Lnorm
     1954    return odfCor,dFdODF
     1955   
     1956def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
     1957    POcorr = 1.0
     1958    if calcControls[phfx+'poType'] == 'MD':
     1959        MD = parmDict[phfx+'MD']
     1960        if MD != 1.0:
     1961            MDAxis = calcControls[phfx+'MDAxis']
     1962            sumMD = 0
     1963            for H in refl[11]:           
     1964                cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
     1965                A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
     1966                sumMD += A**3
     1967            POcorr = sumMD/len(refl[11])
     1968    else:   #spherical harmonics
     1969        if calcControls[pfx+'SHord']:
     1970            POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
     1971    return POcorr
     1972   
     1973def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
     1974    POcorr = 1.0
     1975    POderv = {}
     1976    if calcControls[phfx+'poType'] == 'MD':
     1977        MD = parmDict[phfx+'MD']
     1978        MDAxis = calcControls[phfx+'MDAxis']
     1979        sumMD = 0
     1980        sumdMD = 0
     1981        for H in refl[11]:           
     1982            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
     1983            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
     1984            sumMD += A**3
     1985            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
     1986        POcorr = sumMD/len(refl[11])
     1987        POderv[phfx+'MD'] = sumdMD/len(refl[11])
     1988    else:   #spherical harmonics
     1989        if calcControls[pfx+'SHord']:
     1990            POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
     1991    return POcorr,POderv
     1992   
     1993def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
     1994    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
     1995    if 'X' in parmDict[hfx+'Type']:
     1996        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
     1997    Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
     1998    if pfx+'SHorder' in parmDict:
     1999        Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
     2000    refl[13] = Icorr       
     2001   
     2002def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
     2003    dIdsh = 1./parmDict[hfx+'Scale']
     2004    dIdsp = 1./parmDict[phfx+'Scale']
     2005    if 'X' in parmDict[hfx+'Type']:
     2006        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
     2007        dIdPola /= pola
     2008    else:       #'N'
     2009        dIdPola = 0.0
     2010    POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
     2011    for iPO in dIdPO:
     2012        dIdPO[iPO] /= POcorr
     2013    dFdODF = {}
     2014    dFdSA = [0,0,0]
     2015    if pfx+'SHorder' in parmDict:
     2016        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
     2017        for iSH in dFdODF:
     2018            dFdODF[iSH] /= odfCor
     2019        for i in range(3):
     2020            dFdSA[i] /= odfCor
     2021    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA
     2022       
     2023def GetSampleGam(refl,wave,G,GB,phfx,calcControls,parmDict):
     2024    costh = cosd(refl[5]/2.)
     2025    #crystallite size
     2026    if calcControls[phfx+'SizeType'] == 'isotropic':
     2027        gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:i']*costh)
     2028    elif calcControls[phfx+'SizeType'] == 'uniaxial':
     2029        H = np.array(refl[:3])
     2030        P = np.array(calcControls[phfx+'SizeAxis'])
     2031        cosP,sinP = G2lat.CosSinAngle(H,P,G)
     2032        gam = (1.8*wave/np.pi)/(parmDict[phfx+'Size:i']*parmDict[phfx+'Size:a']*costh)
     2033        gam *= np.sqrt((sinP*parmDict[phfx+'Size:a'])**2+(cosP*parmDict[phfx+'Size:i'])**2)
     2034    else:           #ellipsoidal crystallites
     2035        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
     2036        H = np.array(refl[:3])
     2037        lenR = G2pwd.ellipseSize(H,Sij,GB)
     2038        gam = 1.8*wave/(np.pi*costh*lenR)
     2039    #microstrain               
     2040    if calcControls[phfx+'MustrainType'] == 'isotropic':
     2041        gam += 0.018*parmDict[phfx+'Mustrain:i']*tand(refl[5]/2.)/np.pi
     2042    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
     2043        H = np.array(refl[:3])
     2044        P = np.array(calcControls[phfx+'MustrainAxis'])
     2045        cosP,sinP = G2lat.CosSinAngle(H,P,G)
     2046        Si = parmDict[phfx+'Mustrain:i']
     2047        Sa = parmDict[phfx+'Mustrain:a']
     2048        gam += 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
     2049    else:       #generalized - P.W. Stephens model
     2050        pwrs = calcControls[phfx+'MuPwrs']
     2051        sum = 0
     2052        for i,pwr in enumerate(pwrs):
     2053            sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
     2054        gam += 0.018*refl[4]**2*tand(refl[5]/2.)*sum           
     2055    return gam
     2056       
     2057def GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict):
     2058    gamDict = {}
     2059    costh = cosd(refl[5]/2.)
     2060    tanth = tand(refl[5]/2.)
     2061    #crystallite size derivatives
     2062    if calcControls[phfx+'SizeType'] == 'isotropic':
     2063        gamDict[phfx+'Size:i'] = -1.80*wave/(np.pi*costh)
     2064    elif calcControls[phfx+'SizeType'] == 'uniaxial':
     2065        H = np.array(refl[:3])
     2066        P = np.array(calcControls[phfx+'SizeAxis'])
     2067        cosP,sinP = G2lat.CosSinAngle(H,P,G)
     2068        Si = parmDict[phfx+'Size:i']
     2069        Sa = parmDict[phfx+'Size:a']
     2070        gami = (1.8*wave/np.pi)/(Si*Sa)
     2071        sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
     2072        gam = gami*sqtrm/costh           
     2073        gamDict[phfx+'Size:i'] = gami*Si*cosP**2/(sqtrm*costh)-gam/Si
     2074        gamDict[phfx+'Size:a'] = gami*Sa*sinP**2/(sqtrm*costh)-gam/Sa         
     2075    else:           #ellipsoidal crystallites
     2076        const = 1.8*wave/(np.pi*costh)
     2077        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
     2078        H = np.array(refl[:3])
     2079        R,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
     2080        for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
     2081            gamDict[item] = -(const/R**2)*dRdS[i]
     2082    #microstrain derivatives               
     2083    if calcControls[phfx+'MustrainType'] == 'isotropic':
     2084        gamDict[phfx+'Mustrain:i'] =  0.018*tanth/np.pi           
     2085    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
     2086        H = np.array(refl[:3])
     2087        P = np.array(calcControls[phfx+'MustrainAxis'])
     2088        cosP,sinP = G2lat.CosSinAngle(H,P,G)
     2089        Si = parmDict[phfx+'Mustrain:i']
     2090        Sa = parmDict[phfx+'Mustrain:a']
     2091        gami = 0.018*Si*Sa*tanth/np.pi
     2092        sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
     2093        gam = gami/sqtrm
     2094        gamDict[phfx+'Mustrain:i'] = gam/Si-gami*Si*cosP**2/sqtrm**3
     2095        gamDict[phfx+'Mustrain:a'] = gam/Sa-gami*Sa*sinP**2/sqtrm**3
     2096    else:       #generalized - P.W. Stephens model
     2097        pwrs = calcControls[phfx+'MuPwrs']
     2098        const = 0.018*refl[4]**2*tanth
     2099        for i,pwr in enumerate(pwrs):
     2100            gamDict[phfx+'Mustrain:'+str(i)] = const*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
     2101    return gamDict
     2102       
     2103def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
     2104    h,k,l = refl[:3]
     2105    dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
     2106    d = np.sqrt(dsq)
     2107
     2108    refl[4] = d
     2109    pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
     2110    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
     2111    if 'Bragg' in calcControls[hfx+'instType']:
     2112        pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
     2113            parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
     2114    else:               #Debye-Scherrer - simple but maybe not right
     2115        pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
     2116    return pos
     2117
     2118def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
     2119    dpr = 180./np.pi
     2120    h,k,l = refl[:3]
     2121    dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
     2122    dst = np.sqrt(dstsq)
     2123    pos = refl[5]-parmDict[hfx+'Zero']
     2124    const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
     2125    dpdw = const*dst
     2126    dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
     2127    dpdA *= const*wave/(2.0*dst)
     2128    dpdZ = 1.0
     2129    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
     2130    if 'Bragg' in calcControls[hfx+'instType']:
     2131        dpdSh = -4.*const*cosd(pos/2.0)
     2132        dpdTr = -const*sind(pos)*100.0
     2133        return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
     2134    else:               #Debye-Scherrer - simple but maybe not right
     2135        dpdXd = -const*cosd(pos)
     2136        dpdYd = -const*sind(pos)
     2137        return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
     2138           
     2139def GetHStrainShift(refl,SGData,phfx,parmDict):
     2140    laue = SGData['SGLaue']
     2141    uniq = SGData['SGUniq']
     2142    h,k,l = refl[:3]
     2143    if laue in ['m3','m3m']:
     2144        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
     2145            refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
     2146    elif laue in ['6/m','6/mmm','3m1','31m','3']:
     2147        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
     2148    elif laue in ['3R','3mR']:
     2149        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
     2150    elif laue in ['4/m','4/mmm']:
     2151        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
     2152    elif laue in ['mmm']:
     2153        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
     2154    elif laue in ['2/m']:
     2155        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
     2156        if uniq == 'a':
     2157            Dij += parmDict[phfx+'D23']*k*l
     2158        elif uniq == 'b':
     2159            Dij += parmDict[phfx+'D13']*h*l
     2160        elif uniq == 'c':
     2161            Dij += parmDict[phfx+'D12']*h*k
     2162    else:
     2163        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
     2164            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
     2165    return Dij*refl[4]**2*tand(refl[5]/2.0)
     2166           
     2167def GetHStrainShiftDerv(refl,SGData,phfx):
     2168    laue = SGData['SGLaue']
     2169    uniq = SGData['SGUniq']
     2170    h,k,l = refl[:3]
     2171    if laue in ['m3','m3m']:
     2172        dDijDict = {phfx+'D11':h**2+k**2+l**2,
     2173            phfx+'eA':((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
     2174    elif laue in ['6/m','6/mmm','3m1','31m','3']:
     2175        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
     2176    elif laue in ['3R','3mR']:
     2177        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
     2178    elif laue in ['4/m','4/mmm']:
     2179        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
     2180    elif laue in ['mmm']:
     2181        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
     2182    elif laue in ['2/m']:
     2183        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
     2184        if uniq == 'a':
     2185            dDijDict[phfx+'D23'] = k*l
     2186        elif uniq == 'b':
     2187            dDijDict[phfx+'D13'] = h*l
     2188        elif uniq == 'c':
     2189            dDijDict[phfx+'D12'] = h*k
     2190            names.append()
     2191    else:
     2192        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
     2193            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
     2194    for item in dDijDict:
     2195        dDijDict[item] *= refl[4]**2*tand(refl[5]/2.0)
     2196    return dDijDict
     2197   
     2198def GetFprime(controlDict,Histograms):
     2199    FFtables = controlDict['FFtables']
     2200    if not FFtables:
     2201        return
     2202    histoList = Histograms.keys()
     2203    histoList.sort()
     2204    for histogram in histoList:
     2205        if 'PWDR' in histogram[:4]:
     2206            Histogram = Histograms[histogram]
     2207            hId = Histogram['hId']
     2208            hfx = ':%d:'%(hId)
     2209            keV = controlDict[hfx+'keV']
     2210            for El in FFtables:
     2211                Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
     2212                FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
     2213                FFtables[El][hfx+'FP'] = FP
     2214                FFtables[El][hfx+'FPP'] = FPP               
     2215           
     2216def GetFobsSq(Histograms,Phases,parmDict,calcControls):
     2217    histoList = Histograms.keys()
     2218    histoList.sort()
     2219    for histogram in histoList:
     2220        if 'PWDR' in histogram[:4]:
     2221            Histogram = Histograms[histogram]
     2222            hId = Histogram['hId']
     2223            hfx = ':%d:'%(hId)
     2224            Limits = calcControls[hfx+'Limits']
     2225            shl = max(parmDict[hfx+'SH/L'],0.002)
     2226            Ka2 = False
     2227            kRatio = 0.0
     2228            if hfx+'Lam1' in parmDict.keys():
     2229                Ka2 = True
     2230                lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
     2231                kRatio = parmDict[hfx+'I(L2)/I(L1)']
     2232            x,y,w,yc,yb,yd = Histogram['Data']
     2233            ymb = np.array(y-yb)
     2234            ycmb = np.array(yc-yb)
     2235            ratio = np.where(ycmb!=0.,ymb/ycmb,0.0)           
     2236            xB = np.searchsorted(x,Limits[0])
     2237            xF = np.searchsorted(x,Limits[1])
     2238            refLists = Histogram['Reflection Lists']
     2239            for phase in refLists:
     2240                Phase = Phases[phase]
     2241                pId = Phase['pId']
     2242                phfx = '%d:%d:'%(pId,hId)
     2243                refList = refLists[phase]
     2244                sumFo = 0.0
     2245                sumdF = 0.0
     2246                sumFosq = 0.0
     2247                sumdFsq = 0.0
     2248                for refl in refList:
     2249                    if 'C' in calcControls[hfx+'histType']:
     2250                        yp = np.zeros_like(yb)
     2251                        Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
     2252                        iBeg = np.searchsorted(x[xB:xF],refl[5]-fmin)
     2253                        iFin = np.searchsorted(x[xB:xF],refl[5]+fmax)
     2254                        iFin2 = iFin
     2255                        yp[iBeg:iFin] = refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here                           
     2256                        if Ka2:
     2257                            pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
     2258                            Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
     2259                            iBeg2 = np.searchsorted(x,pos2-fmin)
     2260                            iFin2 = np.searchsorted(x,pos2+fmax)
     2261                            yp[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])        #and here
     2262                        refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[13]*(1.+kRatio)),0.0))
     2263                    elif 'T' in calcControls[hfx+'histType']:
     2264                        print 'TOF Undefined at present'
     2265                        raise Exception    #no TOF yet
     2266                    Fo = np.sqrt(np.abs(refl[8]))
     2267                    Fc = np.sqrt(np.abs(refl[9]))
     2268                    sumFo += Fo
     2269                    sumFosq += refl[8]**2
     2270                    sumdF += np.abs(Fo-Fc)
     2271                    sumdFsq += (refl[8]-refl[9])**2
     2272                Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
     2273                Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
     2274                Histogram[phfx+'Nref'] = len(refList)
     2275               
     2276def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
     2277   
     2278    def GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
     2279        U = parmDict[hfx+'U']
     2280        V = parmDict[hfx+'V']
     2281        W = parmDict[hfx+'W']
     2282        X = parmDict[hfx+'X']
     2283        Y = parmDict[hfx+'Y']
     2284        tanPos = tand(refl[5]/2.0)
     2285        sig = U*tanPos**2+V*tanPos+W        #save peak sigma
     2286        sig = max(0.001,sig)
     2287        gam = X/cosd(refl[5]/2.0)+Y*tanPos+GetSampleGam(refl,wave,G,GB,phfx,calcControls,parmDict) #save peak gamma
     2288        gam = max(0.001,gam)
     2289        return sig,gam
     2290               
     2291    hId = Histogram['hId']
     2292    hfx = ':%d:'%(hId)
     2293    bakType = calcControls[hfx+'bakType']
     2294    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
     2295    yc = np.zeros_like(yb)
     2296       
     2297    if 'C' in calcControls[hfx+'histType']:   
     2298        shl = max(parmDict[hfx+'SH/L'],0.002)
     2299        Ka2 = False
     2300        if hfx+'Lam1' in parmDict.keys():
     2301            wave = parmDict[hfx+'Lam1']
     2302            Ka2 = True
     2303            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
     2304            kRatio = parmDict[hfx+'I(L2)/I(L1)']
     2305        else:
     2306            wave = parmDict[hfx+'Lam']
     2307    else:
     2308        print 'TOF Undefined at present'
     2309        raise ValueError
     2310    for phase in Histogram['Reflection Lists']:
     2311        refList = Histogram['Reflection Lists'][phase]
     2312        Phase = Phases[phase]
     2313        pId = Phase['pId']
     2314        pfx = '%d::'%(pId)
     2315        phfx = '%d:%d:'%(pId,hId)
     2316        hfx = ':%d:'%(hId)
     2317        SGData = Phase['General']['SGData']
     2318        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
     2319        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
     2320        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
     2321        Vst = np.sqrt(nl.det(G))    #V*
     2322        if 'Pawley' not in Phase['General']['Type']:
     2323            refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
     2324        for refl in refList:
     2325            if 'C' in calcControls[hfx+'histType']:
     2326                h,k,l = refl[:3]
     2327                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
     2328                Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
     2329                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
     2330                refl[6:8] = GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)    #peak sig & gam
     2331                GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)    #puts corrections in refl[13]
     2332                refl[13] *= Vst*Lorenz
     2333                if 'Pawley' in Phase['General']['Type']:
     2334                    try:
     2335                        refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
     2336                    except KeyError:
     2337#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
     2338                        continue
     2339                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
     2340                iBeg = np.searchsorted(x,refl[5]-fmin)
     2341                iFin = np.searchsorted(x,refl[5]+fmax)
     2342                if not iBeg+iFin:       #peak below low limit - skip peak
     2343                    continue
     2344                elif not iBeg-iFin:     #peak above high limit - done
     2345                    break
     2346                yc[iBeg:iFin] += refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
     2347                if Ka2:
     2348                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
     2349                    Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
     2350                    iBeg = np.searchsorted(x,pos2-fmin)
     2351                    iFin = np.searchsorted(x,pos2+fmax)
     2352                    if not iBeg+iFin:       #peak below low limit - skip peak
     2353                        continue
     2354                    elif not iBeg-iFin:     #peak above high limit - done
     2355                        return yc,yb
     2356                    yc[iBeg:iFin] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
     2357            elif 'T' in calcControls[hfx+'histType']:
     2358                print 'TOF Undefined at present'
     2359                raise Exception    #no TOF yet
     2360    return yc,yb
     2361   
     2362def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
     2363   
     2364    def cellVaryDerv(pfx,SGData,dpdA):
     2365        if SGData['SGLaue'] in ['-1',]:
     2366            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
     2367                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
     2368        elif SGData['SGLaue'] in ['2/m',]:
     2369            if SGData['SGUniq'] == 'a':
     2370                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
     2371            elif SGData['SGUniq'] == 'b':
     2372                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
     2373            else:
     2374                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
     2375        elif SGData['SGLaue'] in ['mmm',]:
     2376            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
     2377        elif SGData['SGLaue'] in ['4/m','4/mmm']:
     2378#            return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]]
     2379            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
     2380        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
     2381#            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]]
     2382            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
     2383        elif SGData['SGLaue'] in ['3R', '3mR']:
     2384            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
     2385        elif SGData['SGLaue'] in ['m3m','m3']:
     2386#            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]]
     2387            return [[pfx+'A0',dpdA[0]]]
     2388    # create a list of dependent variables and set up a dictionary to hold their derivatives
     2389    dependentVars = G2mv.GetDependentVars()
     2390    depDerivDict = {}
     2391    for j in dependentVars:
     2392        depDerivDict[j] = np.zeros(shape=(len(x)))
     2393    #print 'dependent vars',dependentVars
     2394    lenX = len(x)               
     2395    hId = Histogram['hId']
     2396    hfx = ':%d:'%(hId)
     2397    bakType = calcControls[hfx+'bakType']
     2398    dMdv = np.zeros(shape=(len(varylist),len(x)))
     2399    dMdb,dMddb = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
     2400    if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
     2401        bBpos =varylist.index(hfx+'Back:0')
     2402        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
     2403    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
     2404    for name in varylist:
     2405        if 'Debye' in name:
     2406            id = int(name.split(':')[-1])
     2407            parm = name[:int(name.rindex(':'))]
     2408            ip = names.index(parm)
     2409            dMdv[varylist.index(name)] = dMddb[3*id+ip]
     2410    if 'C' in calcControls[hfx+'histType']:   
     2411        dx = x[1]-x[0]
     2412        shl = max(parmDict[hfx+'SH/L'],0.002)
     2413        Ka2 = False
     2414        if hfx+'Lam1' in parmDict.keys():
     2415            wave = parmDict[hfx+'Lam1']
     2416            Ka2 = True
     2417            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
     2418            kRatio = parmDict[hfx+'I(L2)/I(L1)']
     2419        else:
     2420            wave = parmDict[hfx+'Lam']
     2421    else:
     2422        print 'TOF Undefined at present'
     2423        raise ValueError
     2424    for phase in Histogram['Reflection Lists']:
     2425        refList = Histogram['Reflection Lists'][phase]
     2426        Phase = Phases[phase]
     2427        SGData = Phase['General']['SGData']
     2428        pId = Phase['pId']
     2429        pfx = '%d::'%(pId)
     2430        phfx = '%d:%d:'%(pId,hId)
     2431        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
     2432        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
     2433        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
     2434        if 'Pawley' not in Phase['General']['Type']:
     2435            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
     2436        for iref,refl in enumerate(refList):
     2437            if 'C' in calcControls[hfx+'histType']:        #CW powder
     2438                h,k,l = refl[:3]
     2439                dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
     2440                if 'Pawley' in Phase['General']['Type']:
     2441                    try:
     2442                        refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
     2443                    except KeyError:
     2444#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
     2445                        continue
     2446                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
     2447                iBeg = np.searchsorted(x,refl[5]-fmin)
     2448                iFin = np.searchsorted(x,refl[5]+fmax)
     2449                if not iBeg+iFin:       #peak below low limit - skip peak
     2450                    continue
     2451                elif not iBeg-iFin:     #peak above high limit - done
     2452                    break
     2453                pos = refl[5]
     2454                tanth = tand(pos/2.0)
     2455                costh = cosd(pos/2.0)
     2456                lenBF = iFin-iBeg
     2457                dMdpk = np.zeros(shape=(6,lenBF))
     2458                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
     2459                for i in range(1,5):
     2460                    dMdpk[i] += 100.*dx*refl[13]*refl[9]*dMdipk[i]
     2461                dMdpk[0] += 100.*dx*refl[13]*refl[9]*dMdipk[0]
     2462                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
     2463                if Ka2:
     2464                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
     2465                    kdelt = int((pos2-refl[5])/dx)               
     2466                    iBeg2 = min(lenX,iBeg+kdelt)
     2467                    iFin2 = min(lenX,iFin+kdelt)
     2468                    if iBeg2-iFin2:
     2469                        lenBF2 = iFin2-iBeg2
     2470                        dMdpk2 = np.zeros(shape=(6,lenBF2))
     2471                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
     2472                        for i in range(1,5):
     2473                            dMdpk2[i] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[i]
     2474                        dMdpk2[0] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[0]
     2475                        dMdpk2[5] = 100.*dx*refl[13]*dMdipk2[0]
     2476                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
     2477                if 'Pawley' in Phase['General']['Type']:
     2478                    try:
     2479                        idx = varylist.index(pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]))
     2480                        dMdv[idx][iBeg:iFin] = dervDict['int']/refl[9]
     2481                        if Ka2:
     2482                            dMdv[idx][iBeg2:iFin2] = dervDict2['int']/refl[9]
     2483                        # Assuming Pawley variables not in constraints
     2484                    except ValueError:
     2485                        pass
     2486                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
     2487                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
     2488                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
     2489                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
     2490                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
     2491                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
     2492                    hfx+'DisplaceY':[dpdY,'pos'],}
     2493                for name in names:
     2494                    item = names[name]
     2495                    if name in varylist:
     2496                        dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
     2497                        if Ka2:
     2498                            dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
     2499                    elif name in dependentVars:
     2500                        if Ka2:
     2501                            depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
     2502                        depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
     2503                for iPO in dIdPO:
     2504                    if iPO in varylist:
     2505                        dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
     2506                        if Ka2:
     2507                            dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
     2508                    elif iPO in dependentVars:
     2509                        depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
     2510                        if Ka2:
     2511                            depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
     2512                for i,name in enumerate(['omega','chi','phi']):
     2513                    aname = pfx+'SH '+name
     2514                    if aname in varylist:
     2515                        dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
     2516                        if Ka2:
     2517                            dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
     2518                    elif aname in dependentVars:
     2519                        depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
     2520                        if Ka2:
     2521                            depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
     2522                for iSH in dFdODF:
     2523                    if iSH in varylist:
     2524                        dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
     2525                        if Ka2:
     2526                            dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
     2527                    elif iSH in dependentVars:
     2528                        depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
     2529                        if Ka2:
     2530                            depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
     2531                cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
     2532                for name,dpdA in cellDervNames:
     2533                    if name in varylist:
     2534                        dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
     2535                        if Ka2:
     2536                            dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
     2537                    elif name in dependentVars:
     2538                        depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
     2539                        if Ka2:
     2540                            depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
     2541                dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
     2542                for name in dDijDict:
     2543                    if name in varylist:
     2544                        dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
     2545                        if Ka2:
     2546                            dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
     2547                    elif name in dependentVars:
     2548                        depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
     2549                        if Ka2:
     2550                            depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
     2551                gamDict = GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
     2552                for name in gamDict:
     2553                    if name in varylist:
     2554                        dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
     2555                        if Ka2:
     2556                            dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
     2557                    elif name in dependentVars:
     2558                        depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
     2559                        if Ka2:
     2560                            depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
     2561                                               
     2562            elif 'T' in calcControls[hfx+'histType']:
     2563                print 'TOF Undefined at present'
     2564                raise Exception    #no TOF yet
     2565            #do atom derivatives -  for F,X & U so far             
     2566            corr = dervDict['int']/refl[9]
     2567            if Ka2:
     2568                corr2 = dervDict2['int']/refl[9]
     2569            for name in varylist+dependentVars:
     2570                try:
     2571                    aname = name.split(pfx)[1][:2]
     2572                    if aname not in ['Af','dA','AU']: continue # skip anything not an atom param
     2573                except IndexError:
     2574                    continue
     2575                if name in varylist:
     2576                    dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
     2577                    if Ka2:
     2578                        dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
     2579                elif name in dependentVars:
     2580                    depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
     2581                    if Ka2:
     2582                        depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
     2583    # now process derivatives in constraints
     2584    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
     2585    return dMdv
     2586
     2587def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
     2588    parmdict.update(zip(varylist,values))
     2589    G2mv.Dict2Map(parmdict,varylist)
     2590    Histograms,Phases = HistoPhases
     2591    nvar = len(varylist)
     2592    dMdv = np.empty(0)
     2593    histoList = Histograms.keys()
     2594    histoList.sort()
     2595    for histogram in histoList:
     2596        if 'PWDR' in histogram[:4]:
     2597            Histogram = Histograms[histogram]
     2598            hId = Histogram['hId']
     2599            hfx = ':%d:'%(hId)
     2600            Limits = calcControls[hfx+'Limits']
     2601            x,y,w,yc,yb,yd = Histogram['Data']
     2602            xB = np.searchsorted(x,Limits[0])
     2603            xF = np.searchsorted(x,Limits[1])
     2604            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
     2605                varylist,Histogram,Phases,calcControls,pawleyLookup)
     2606            if len(dMdv):
     2607                dMdv = np.concatenate((dMdv.T,dMdvh.T)).T
     2608            else:
     2609                dMdv = dMdvh
     2610    return dMdv
     2611
     2612def ComputePowderHessian(args):
     2613    Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup = args
     2614    hId = Histogram['hId']
     2615    hfx = ':%d:'%(hId)
     2616    Limits = calcControls[hfx+'Limits']
     2617    x,y,w,yc,yb,yd = Histogram['Data']
     2618    dy = y-yc
     2619    xB = np.searchsorted(x,Limits[0])
     2620    xF = np.searchsorted(x,Limits[1])
     2621    dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(
     2622        parmdict,x[xB:xF],
     2623        varylist,Histogram,Phases,calcControls,pawleyLookup)
     2624    Vec = np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1)
     2625    Hess = np.inner(dMdvh,dMdvh)
     2626    return Vec,Hess
     2627
     2628#def HessRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
     2629#    parmdict.update(zip(varylist,values))
     2630#    G2mv.Dict2Map(parmdict,varylist)
     2631#    Histograms,Phases = HistoPhases
     2632#    nvar = len(varylist)
     2633#    Hess = np.empty(0)
     2634#    histoList = Histograms.keys()
     2635#    histoList.sort()
     2636#    for histogram in histoList:
     2637#        if 'PWDR' in histogram[:4]:
     2638#            Histogram = Histograms[histogram]
     2639#            hId = Histogram['hId']
     2640#            hfx = ':%d:'%(hId)
     2641#            Limits = calcControls[hfx+'Limits']
     2642#            x,y,w,yc,yb,yd = Histogram['Data']
     2643#            dy = y-yc
     2644#            xB = np.searchsorted(x,Limits[0])
     2645#            xF = np.searchsorted(x,Limits[1])
     2646#            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
     2647#                varylist,Histogram,Phases,calcControls,pawleyLookup)
     2648#            if dlg:
     2649#                dlg.Update(Histogram['wRp'],newmsg='Hessian for histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0]
     2650#            if len(Hess):
     2651#                Vec += np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1)
     2652#                Hess += np.inner(dMdvh,dMdvh)
     2653#            else:
     2654#                Vec = np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1)
     2655#                Hess = np.inner(dMdvh,dMdvh)
     2656#    return Vec,Hess
     2657
     2658def HessRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
     2659    parmdict.update(zip(varylist,values))
     2660    G2mv.Dict2Map(parmdict,varylist)
     2661    Histograms,Phases = HistoPhases
     2662    nvar = len(varylist)
     2663    HessSum = np.zeros((nvar,nvar))
     2664    VecSum = np.zeros(nvar)
     2665    histoList = Histograms.keys()
     2666    histoList.sort()
     2667    argList = []
     2668    MaxProcess = calcControls['max Hprocess']
     2669    for histogram in histoList:
     2670        if 'PWDR' in histogram[:4]:
     2671            Histogram = Histograms[histogram]
     2672            argList.append([
     2673                Histogram,parmdict,varylist,Phases,
     2674                calcControls,pawleyLookup])
     2675    if MaxProcess > 1:
     2676        mpPool = mp.Pool(processes=MaxProcess)
     2677        results = mpPool.map(ComputePowderHessian,argList)
     2678        for Vec,Hess in results:
     2679            VecSum += Vec
     2680            HessSum += Hess
     2681    else:
     2682        for arg in argList:
     2683            Vec,Hess = ComputePowderHessian(arg)
     2684            VecSum += Vec
     2685            HessSum += Hess
     2686    return VecSum,HessSum
     2687
     2688def ComputePowderProfile(args):
     2689    Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup = args
     2690    hId = Histogram['hId']
     2691    hfx = ':%d:'%(hId)
     2692    x,y,w,yc,yb,yd = Histogram['Data']
     2693    Limits = calcControls[hfx+'Limits']
     2694    xB = np.searchsorted(x,Limits[0])
     2695    xF = np.searchsorted(x,Limits[1])
     2696    yc,yb = getPowderProfile(parmdict,x[xB:xF],
     2697                            varylist,Histogram,Phases,calcControls,
     2698                            pawleyLookup)
     2699    return xB,xF,yc,yb,Histogram['Reflection Lists']
     2700
     2701#def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
     2702#    parmdict.update(zip(varylist,values))
     2703#    Values2Dict(parmdict, varylist, values)
     2704#    G2mv.Dict2Map(parmdict,varylist)
     2705#    Histograms,Phases = HistoPhases
     2706#    M = np.empty(0)
     2707#    sumwYo = 0
     2708#    Nobs = 0
     2709#    histoList = Histograms.keys()
     2710#    histoList.sort()
     2711#    for histogram in histoList:
     2712#        if 'PWDR' in histogram[:4]:
     2713#            Histogram = Histograms[histogram]
     2714#            hId = Histogram['hId']
     2715#            hfx = ':%d:'%(hId)
     2716#            Limits = calcControls[hfx+'Limits']
     2717#            x,y,w,yc,yb,yd = Histogram['Data']
     2718#            yc *= 0.0                           #zero full calcd profiles
     2719#            yb *= 0.0
     2720#            yd *= 0.0
     2721#            xB = np.searchsorted(x,Limits[0])
     2722#            xF = np.searchsorted(x,Limits[1])
     2723#            Histogram['Nobs'] = xF-xB
     2724#            Nobs += Histogram['Nobs']
     2725#            Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
     2726#            sumwYo += Histogram['sumwYo']
     2727#            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
     2728#                varylist,Histogram,Phases,calcControls,pawleyLookup)
     2729#            yc[xB:xF] += yb[xB:xF]
     2730#            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
     2731#            Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
     2732#            wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF])
     2733#            Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
     2734#            if dlg:
     2735#                dlg.Update(Histogram['wRp'],newmsg='For histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0]
     2736#            M = np.concatenate((M,wdy))
     2737#    Histograms['sumwYo'] = sumwYo
     2738#    Histograms['Nobs'] = Nobs
     2739#    Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)
     2740#    if dlg:
     2741#        GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0]
     2742#        if not GoOn:
     2743#            parmDict['saved values'] = values
     2744#            raise Exception         #Abort!!
     2745#    return M
     2746#   
     2747def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
     2748    parmdict.update(zip(varylist,values))
     2749    Values2Dict(parmdict, varylist, values)
     2750    G2mv.Dict2Map(parmdict,varylist)
     2751    Histograms,Phases = HistoPhases
     2752    M = np.empty(0)
     2753    sumwYo = 0
     2754    Nobs = 0
     2755    histoList = Histograms.keys()
     2756    histoList.sort()
     2757    argList = []
     2758    MaxProcess = calcControls['max Hprocess']
     2759    for histogram in histoList:
     2760        if 'PWDR' in histogram[:4]:
     2761            Histogram = Histograms[histogram]
     2762            argList.append(
     2763                [Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup]
     2764                )
     2765    if MaxProcess > 1:
     2766        mpPool = mp.Pool(processes=MaxProcess)
     2767        results = mpPool.map(ComputePowderProfile,argList)
     2768        for arg,res in zip(argList,results):
     2769            xB,xF,ycSect,ybSect,RL = res
     2770            Histogram = arg[0]
     2771            Histogram['Reflection Lists'] = RL
     2772            x,y,w,yc,yb,yd = Histogram['Data']
     2773            yc *= 0.0                           #zero full calcd profiles
     2774            yb *= 0.0
     2775            yd *= 0.0
     2776            Histogram['Nobs'] = xF-xB
     2777            Nobs += Histogram['Nobs']
     2778            Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
     2779            sumwYo += Histogram['sumwYo']
     2780           
     2781            yc[xB:xF] = ycSect
     2782            yb[xB:xF] = ybSect
     2783            yc[xB:xF] += yb[xB:xF]
     2784            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
     2785            Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
     2786            wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF])
     2787            Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
     2788            M = np.concatenate((M,wdy))
     2789    else:
     2790        for arg in argList:
     2791            xB,xF,ycSect,ybSect,RL = ComputePowderProfile(arg)
     2792            Histogram = arg[0]
     2793            hId = Histogram['hId']
     2794            x,y,w,yc,yb,yd = Histogram['Data']
     2795            yc *= 0.0                           #zero full calcd profiles
     2796            yb *= 0.0
     2797            yd *= 0.0
     2798            Histogram['Nobs'] = xF-xB
     2799            Nobs += Histogram['Nobs']
     2800            Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
     2801            sumwYo += Histogram['sumwYo']
     2802           
     2803            yc[xB:xF] = ycSect
     2804            yb[xB:xF] = ybSect
     2805            yc[xB:xF] += yb[xB:xF]
     2806            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
     2807            Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
     2808            wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF])
     2809            Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
     2810            if dlg:
     2811                dlg.Update(Histogram['wRp'],newmsg='For histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0]
     2812            M = np.concatenate((M,wdy))
     2813
     2814    Histograms['sumwYo'] = sumwYo
     2815    Histograms['Nobs'] = Nobs
     2816    Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)
     2817    if dlg:
     2818        GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0]
     2819        if not GoOn:
     2820            parmDict['saved values'] = values
     2821            raise Exception         #Abort!!
     2822    return M       
     2823                   
     2824def Refine(GPXfile,dlg):
     2825    import pytexture as ptx
     2826    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
     2827   
     2828    ShowBanner()
     2829    varyList = []
     2830    parmDict = {}
     2831    G2mv.InitVars()   
     2832    Controls = GetControls(GPXfile)
     2833    ShowControls(Controls)
     2834    calcControls = {}
     2835    calcControls.update(Controls)           
     2836    constrDict,constrFlag,fixedList = GetConstraints(GPXfile)
     2837    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
     2838    if not Phases:
     2839        print ' *** ERROR - you have no histograms to refine! ***'
     2840        print ' *** Refine aborted ***'
     2841        raise Exception
     2842    if not Histograms:
     2843        print ' *** ERROR - you have no data to refine with! ***'
     2844        print ' *** Refine aborted ***'
     2845        raise Exception       
     2846    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases)
     2847    calcControls['Natoms'] = Natoms
     2848    calcControls['FFtables'] = FFtables
     2849    calcControls['BLtables'] = BLtables
     2850    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
     2851    calcControls.update(controlDict)
     2852    histVary,histDict,controlDict = GetHistogramData(Histograms)
     2853    calcControls.update(controlDict)
     2854    varyList = phaseVary+hapVary+histVary
     2855    parmDict.update(phaseDict)
     2856    parmDict.update(hapDict)
     2857    parmDict.update(histDict)
     2858    GetFprime(calcControls,Histograms)
     2859    # do constraint processing
     2860    try:
     2861        groups,parmlist = G2mv.GroupConstraints(constrDict)
     2862        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList)
     2863    except:
     2864        print ' *** ERROR - your constraints are internally inconsistent ***'
     2865        raise Exception(' *** Refine aborted ***')
     2866    # check to see which generated parameters are fully varied
     2867    msg = G2mv.SetVaryFlags(varyList)
     2868    if msg:
     2869        print ' *** ERROR - you have not set the refine flags for constraints consistently! ***'
     2870        print msg
     2871        raise Exception(' *** Refine aborted ***')
     2872    G2mv.Map2Dict(parmDict,varyList)
     2873#    print G2mv.VarRemapShow(varyList)
    5342874
    5352875    while True:
    536 
     2876        begin = time.time()
     2877        values =  np.array(Dict2Values(parmDict, varyList))
     2878        Ftol = Controls['min dM/M']
     2879        Factor = Controls['shift factor']
     2880        maxCyc = Controls['max cyc']
     2881        if 'Jacobian' in Controls['deriv type']:           
     2882            result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
     2883                ftol=Ftol,col_deriv=True,factor=Factor,
     2884                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
     2885            ncyc = int(result[2]['nfev']/2)
     2886        elif 'Hessian' in Controls['deriv type']:
     2887            result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
     2888                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
     2889            ncyc = result[2]['num cyc']+1                           
     2890        else:           #'numeric'
     2891            result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
     2892                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
     2893            ncyc = int(result[2]['nfev']/len(varyList))
     2894#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
     2895#        for item in table: print item,table[item]               #useful debug - are things shifting?
     2896        runtime = time.time()-begin
     2897        chisq = np.sum(result[2]['fvec']**2)
     2898        Values2Dict(parmDict, varyList, result[0])
     2899        G2mv.Dict2Map(parmDict,varyList)
     2900       
     2901        Rwp = np.sqrt(chisq/Histograms['sumwYo'])*100.      #to %
     2902        GOF = chisq/(Histograms['Nobs']-len(varyList))
     2903        print '\n Refinement results:'
     2904        print 135*'-'
     2905        print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
     2906        print ' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
     2907        print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
    5372908        try:
    538 
    539             data = cPickle.load(infile)
    540 
    541         except EOFError:
    542 
    543             break
    544 
    545         datum = data[0]
    546 
    547 #        print 'read: ',datum[0]
    548 
    549         if datum[0] == 'Phases':
    550 
    551             for iphase in range(len(data)):
    552 
    553                 if data[iphase][0] in Phases:
    554 
    555                     phaseName = data[iphase][0]
    556 
    557                     data[iphase][1].update(Phases[phaseName])
    558 
    559         elif datum[0] == 'Covariance':
    560 
    561             data[0][1] = CovData
    562 
    563         try:
    564 
    565             histogram = Histograms[datum[0]]
    566 
    567 #            print 'found ',datum[0]
    568 
    569             data[0][1][1] = histogram['Data']
    570 
    571             for datus in data[1:]:
    572 
    573 #                print '    read: ',datus[0]
    574 
    575                 if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
    576 
    577                     datus[1] = histogram[datus[0]]
    578 
    579         except KeyError:
    580 
    581             pass
    582 
    583                                
    584 
    585         cPickle.dump(data,outfile,1)
    586 
    587     infile.close()
    588 
    589     outfile.close()
    590 
    591     print 'GPX file save successful'
    592 
    593    
    594 
    595 def SetSeqResult(GPXfile,Histograms,SeqResult):
    596 
    597     GPXback = GPXBackup(GPXfile)
    598 
    599     print '\n',135*'-'
    600 
    601     print 'Read from file:',GPXback
    602 
    603     print 'Save to file  :',GPXfile
    604 
    605     infile = open(GPXback,'rb')
    606 
    607     outfile = open(GPXfile,'wb')
    608 
    609     while True:
    610 
    611         try:
    612 
    613             data = cPickle.load(infile)
    614 
    615         except EOFError:
    616 
    617             break
    618 
    619         datum = data[0]
    620 
    621         if datum[0] == 'Sequental results':
    622 
    623             data[0][1] = SeqResult
    624 
    625         try:
    626 
    627             histogram = Histograms[datum[0]]
    628 
    629             data[0][1][1] = histogram['Data']
    630 
    631             for datus in data[1:]:
    632 
    633                 if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
    634 
    635                     datus[1] = histogram[datus[0]]
    636 
    637         except KeyError:
    638 
    639             pass
    640 
    641                                
    642 
    643         cPickle.dump(data,outfile,1)
    644 
    645     infile.close()
    646 
    647     outfile.close()
    648 
    649     print 'GPX file save successful'
    650 
    651                        
    652 
    653 def GetPWDRdata(GPXfile,PWDRname):
    654 
    655     ''' Returns powder data from GSASII gpx file
    656 
    657     input:
    658 
    659         GPXfile = .gpx full file name
    660 
    661         PWDRname = powder histogram name as obtained from GetHistogramNames
    662 
    663     return:
    664 
    665         PWDRdata = powder data dictionary with:
    666 
    667             Data - powder data arrays, Limits, Instrument Parameters, Sample Parameters
    668 
    669        
    670 
    671     '''
    672 
    673     file = open(GPXfile,'rb')
    674 
    675     PWDRdata = {}
    676 
    677     while True:
    678 
    679         try:
    680 
    681             data = cPickle.load(file)
    682 
    683         except EOFError:
    684 
    685             break
    686 
    687         datum = data[0]
    688 
    689         if datum[0] == PWDRname:
    690 
    691             PWDRdata['Data'] = datum[1][1]          #powder data arrays
    692 
    693             PWDRdata[data[2][0]] = data[2][1]       #Limits
    694 
    695             PWDRdata[data[3][0]] = data[3][1]       #Background
    696 
    697             PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
    698 
    699             PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
    700 
     2909            covMatrix = result[1]*GOF
     2910            sig = np.sqrt(np.diag(covMatrix))
     2911            if np.any(np.isnan(sig)):
     2912                print '*** Least squares aborted - some invalid esds possible ***'
     2913#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
     2914#            for item in table: print item,table[item]               #useful debug - are things shifting?
     2915            break                   #refinement succeeded - finish up!
     2916        except TypeError:          #result[1] is None on singular matrix
     2917            print '**** Refinement failed - singular matrix ****'
     2918            if 'Hessian' in Controls['deriv type']:
     2919                for i in result[2]['psing'].reverse():
     2920                        print 'Removing parameter: ',varyList[i]
     2921                        del(varyList[i])                   
     2922            else:
     2923                Ipvt = result[2]['ipvt']
     2924                for i,ipvt in enumerate(Ipvt):
     2925                    if not np.sum(result[2]['fjac'],axis=1)[i]:
     2926                        print 'Removing parameter: ',varyList[ipvt-1]
     2927                        del(varyList[ipvt-1])
     2928                        break
     2929
     2930#    print 'dependentParmList: ',G2mv.dependentParmList
     2931#    print 'arrayList: ',G2mv.arrayList
     2932#    print 'invarrayList: ',G2mv.invarrayList
     2933#    print 'indParmList: ',G2mv.indParmList
     2934#    print 'fixedDict: ',G2mv.fixedDict
     2935#    print 'test1'
     2936    GetFobsSq(Histograms,Phases,parmDict,calcControls)
     2937#    print 'test2'
     2938    sigDict = dict(zip(varyList,sig))
     2939    newCellDict = GetNewCellParms(parmDict,varyList)
     2940    newAtomDict = ApplyXYZshifts(parmDict,varyList)
     2941    covData = {'variables':result[0],'varyList':varyList,'sig':sig,
     2942        'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
     2943    # add the uncertainties into the esd dictionary (sigDict)
     2944    sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
     2945    SetPhaseData(parmDict,sigDict,Phases,covData)
     2946    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
     2947    SetHistogramData(parmDict,sigDict,Histograms)
     2948    G2mv.PrintIndependentVars(parmDict,varyList,sigDict)
     2949    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData)
     2950   
     2951#for testing purposes!!!
     2952#    import cPickle
     2953#    file = open('structTestdata.dat','wb')
     2954#    cPickle.dump(parmDict,file,1)
     2955#    cPickle.dump(varyList,file,1)
     2956#    for histogram in Histograms:
     2957#        if 'PWDR' in histogram[:4]:
     2958#            Histogram = Histograms[histogram]
     2959#    cPickle.dump(Histogram,file,1)
     2960#    cPickle.dump(Phases,file,1)
     2961#    cPickle.dump(calcControls,file,1)
     2962#    cPickle.dump(pawleyLookup,file,1)
     2963#    file.close()
     2964
     2965    if dlg:
     2966        return Rwp
     2967
     2968def SeqRefine(GPXfile,dlg):
     2969    import pytexture as ptx
     2970    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
     2971   
     2972    ShowBanner()
     2973    print ' Sequential Refinement'
     2974    G2mv.InitVars()   
     2975    Controls = GetControls(GPXfile)
     2976    ShowControls(Controls)           
     2977    constrDict,constrFlag,fixedList = GetConstraints(GPXfile)
     2978    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
     2979    if not Phases:
     2980        print ' *** ERROR - you have no histograms to refine! ***'
     2981        print ' *** Refine aborted ***'
     2982        raise Exception
     2983    if not Histograms:
     2984        print ' *** ERROR - you have no data to refine with! ***'
     2985        print ' *** Refine aborted ***'
     2986        raise Exception
     2987    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,False)
     2988    if 'Seq Data' in Controls:
     2989        histNames = Controls['Seq Data']
     2990    else:
     2991        histNames = GetHistogramNames(GPXfile,['PWDR',])
     2992    if 'Reverse Seq' in Controls:
     2993        if Controls['Reverse Seq']:
     2994            histNames.reverse()
     2995    SeqResult = {'histNames':histNames}
     2996    makeBack = True
     2997    for ihst,histogram in enumerate(histNames):
     2998        ifPrint = False
     2999        if dlg:
     3000            dlg.SetTitle('Residual for histogram '+str(ihst))
     3001        calcControls = {}
     3002        calcControls['Natoms'] = Natoms
     3003        calcControls['FFtables'] = FFtables
     3004        calcControls['BLtables'] = BLtables
     3005        varyList = []
     3006        parmDict = {}
     3007        Histo = {histogram:Histograms[histogram],}
     3008        hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histo,False)
     3009        calcControls.update(controlDict)
     3010        histVary,histDict,controlDict = GetHistogramData(Histo,False)
     3011        calcControls.update(controlDict)
     3012        varyList = phaseVary+hapVary+histVary
     3013        if not ihst:
     3014            saveVaryList = varyList[:]
     3015            for i,item in enumerate(saveVaryList):
     3016                items = item.split(':')
     3017                if items[1]:
     3018                    items[1] = ''
     3019                item = ':'.join(items)
     3020                saveVaryList[i] = item
     3021            SeqResult['varyList'] = saveVaryList
     3022        else:
     3023            newVaryList = varyList[:]
     3024            for i,item in enumerate(newVaryList):
     3025                items = item.split(':')
     3026                if items[1]:
     3027                    items[1] = ''
     3028                item = ':'.join(items)
     3029                newVaryList[i] = item
     3030            if newVaryList != SeqResult['varyList']:
     3031                print newVaryList
     3032                print SeqResult['varyList']
     3033                print '**** ERROR - variable list for this histogram does not match previous'
     3034                raise Exception
     3035        parmDict.update(phaseDict)
     3036        parmDict.update(hapDict)
     3037        parmDict.update(histDict)
     3038        GetFprime(calcControls,Histo)
     3039        constrDict,constrFlag,fixedList = G2mv.InputParse([])        #constraints go here?
     3040        groups,parmlist = G2mv.GroupConstraints(constrDict)
     3041        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList)
     3042        G2mv.Map2Dict(parmDict,varyList)
     3043   
     3044        while True:
     3045            begin = time.time()
     3046            values =  np.array(Dict2Values(parmDict, varyList))
     3047            Ftol = Controls['min dM/M']
     3048            Factor = Controls['shift factor']
     3049            maxCyc = Controls['max cyc']
     3050
     3051            if 'Jacobian' in Controls['deriv type']:           
     3052                result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
     3053                    ftol=Ftol,col_deriv=True,factor=Factor,
     3054                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
     3055                ncyc = int(result[2]['nfev']/2)
     3056            elif 'Hessian' in Controls['deriv type']:
     3057                result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
     3058                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
     3059                ncyc = result[2]['num cyc']+1                           
     3060            else:           #'numeric'
     3061                result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
     3062                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
     3063                ncyc = int(result[2]['nfev']/len(varyList))
     3064
     3065
     3066
     3067            runtime = time.time()-begin
     3068            chisq = np.sum(result[2]['fvec']**2)
     3069            Values2Dict(parmDict, varyList, result[0])
     3070            G2mv.Dict2Map(parmDict,varyList)
     3071           
     3072            Rwp = np.sqrt(chisq/Histo['sumwYo'])*100.      #to %
     3073            GOF = chisq/(Histo['Nobs']-len(varyList))
     3074            print '\n Refinement results for histogram: v'+histogram
     3075            print 135*'-'
     3076            print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histo['Nobs'],' Number of parameters: ',len(varyList)
     3077            print ' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
     3078            print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
    7013079            try:
    702 
    703                 PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
    704 
    705             except IndexError:
    706 
    707                 PWDRdata['Reflection lists'] = {}
    708 
    709     file.close()
    710 
    711     return PWDRdata
    712 
    713    
    714 
    715 def GetHKLFdata(GPXfile,HKLFname):
    716 
    717     ''' Returns single crystal data from GSASII gpx file
    718 
    719     input:
    720 
    721         GPXfile = .gpx full file name
    722 
    723         HKLFname = single crystal histogram name as obtained from GetHistogramNames
    724 
    725     return:
    726 
    727         HKLFdata = single crystal data list of reflections: for each reflection:
    728 
    729             HKLF = [np.array([h,k,l]),FoSq,sigFoSq,FcSq,Fcp,Fcpp,phase]
    730 
    731     '''
    732 
    733     file = open(GPXfile,'rb')
    734 
    735     HKLFdata = []
    736 
    737     while True:
    738 
    739         try:
    740 
    741             data = cPickle.load(file)
    742 
    743         except EOFError:
    744 
    745             break
    746 
    747         datum = data[0]
    748 
    749         if datum[0] == HKLFname:
    750 
    751             HKLFdata = datum[1:][0]
    752 
    753     file.close()
    754 
    755     return HKLFdata
    756 
    757    
    758 
    759 def ShowBanner():
    760 
    761     print 80*'*'
    762 
    763     print '   General Structure Analysis System-II Crystal Structure Refinement'
    764 
    765     print '     by Robert B. Von Dreele, Argonne National Laboratory(C), 2010'
    766 
    767     print ' This product includes software developed by the UChicago Argonne, LLC,'
    768 
    769     print '            as Operator of Argonne National Laboratory.'
    770 
    771     print 80*'*','\n'
    772 
    773 
    774 
    775 def ShowControls(Controls):
    776 
    777     print ' Least squares controls:'
    778 
    779     print ' Derivative type: ',Controls['deriv type']
    780 
    781     print ' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
    782 
    783     print ' Initial shift factor: ','%.3f'%(Controls['shift factor'])
    784 
    785    
    786 
    787 def GetFFtable(General):
    788 
    789     ''' returns a dictionary of form factor data for atom types found in General
    790 
    791     input:
    792 
    793         General = dictionary of phase info.; includes AtomTypes
    794 
    795     return:
    796 
    797         FFtable = dictionary of form factor data; key is atom type
    798 
    799     '''
    800 
    801     atomTypes = General['AtomTypes']
    802 
    803     FFtable = {}
    804 
    805     for El in atomTypes:
    806 
    807         FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])
    808 
    809         for item in FFs:
    810 
    811             if item['Symbol'] == El.upper():
    812 
    813                 FFtable[El] = item
    814 
    815     return FFtable
    816 
    817    
    818 
    819 def GetBLtable(General):
    820 
    821     ''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General
    822 
    823     input:
    824 
    825         General = dictionary of phase info.; includes AtomTypes & Isotopes
    826 
    827     return:
    828 
    829         BLtable = dictionary of scattering length data; key is atom type
    830 
    831     '''
    832 
    833     atomTypes = General['AtomTypes']
    834 
    835     BLtable = {}
    836 
    837     isotopes = General['Isotopes']
    838 
    839     isotope = General['Isotope']
    840 
    841     for El in atomTypes:
    842 
    843         BLtable[El] = [isotope[El],isotopes[El][isotope[El]]]
    844 
    845     return BLtable
    846 
    847        
    848 
    849 def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
    850 
    851     if SGLaue in ['-1','2/m','mmm']:
    852 
    853         return                      #no Pawley symmetry required constraints
    854 
    855     for i,varyI in enumerate(pawleyVary):
    856 
    857         refI = int(varyI.split(':')[-1])
    858 
    859         ih,ik,il = PawleyRef[refI][:3]
    860 
    861         for varyJ in pawleyVary[0:i]:
    862 
    863             refJ = int(varyJ.split(':')[-1])
    864 
    865             jh,jk,jl = PawleyRef[refJ][:3]
    866 
    867             if SGLaue in ['4/m','4/mmm']:
    868 
    869                 isum = ih**2+ik**2
    870 
    871                 jsum = jh**2+jk**2
    872 
    873                 if abs(il) == abs(jl) and isum == jsum:
    874 
    875                     G2mv.StoreEquivalence(varyJ,(varyI,))
    876 
    877             elif SGLaue in ['3R','3mR']:
    878 
    879                 isum = ih**2+ik**2+il**2
    880 
    881                 jsum = jh**2+jk**2*jl**2
    882 
    883                 isum2 = ih*ik+ih*il+ik*il
    884 
    885                 jsum2 = jh*jk+jh*jl+jk*jl
    886 
    887                 if isum == jsum and isum2 == jsum2:
    888 
    889                     G2mv.StoreEquivalence(varyJ,(varyI,))
    890 
    891             elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
    892 
    893                 isum = ih**2+ik**2+ih*ik
    894 
    895                 jsum = jh**2+jk**2+jh*jk
    896 
    897                 if abs(il) == abs(jl) and isum == jsum:
    898 
    899                     G2mv.StoreEquivalence(varyJ,(varyI,))
    900 
    901             elif SGLaue in ['m3','m3m']:
    902 
    903                 isum = ih**2+ik**2+il**2
    904 
    905                 jsum = jh**2+jk**2+jl**2
    906 
    907                 if isum == jsum:
    908 
    909                     G2mv.StoreEquivalence(varyJ,(varyI,))
    910 
    911                    
    912 
    913 def cellVary(pfx,SGData):
    914 
    915     if SGData['SGLaue'] in ['-1',]:
    916 
    917         return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
    918 
    919     elif SGData['SGLaue'] in ['2/m',]:
    920 
    921         if SGData['SGUniq'] == 'a':
    922 
    923             return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
    924 
    925         elif SGData['SGUniq'] == 'b':
    926 
    927             return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
    928 
    929         else:
    930 
    931             return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
    932 
    933     elif SGData['SGLaue'] in ['mmm',]:
    934 
    935         return [pfx+'A0',pfx+'A1',pfx+'A2']
    936 
    937     elif SGData['SGLaue'] in ['4/m','4/mmm']:
    938 
    939         return [pfx+'A0',pfx+'A2']
    940 
    941     elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
    942 
    943         return [pfx+'A0',pfx+'A2']
    944 
    945     elif SGData['SGLaue'] in ['3R', '3mR']:
    946 
    947         return [pfx+'A0',pfx+'A3']                       
    948 
    949     elif SGData['SGLaue'] in ['m3m','m3']:
    950 
    951         return [pfx+'A0',]
    952 
    953                    
    954 
    955 def GetPhaseData(PhaseData,Print=True):
     3080                covMatrix = result[1]*GOF
     3081                sig = np.sqrt(np.diag(covMatrix))
     3082                if np.any(np.isnan(sig)):
     3083                    print '*** Least squares aborted - some invalid esds possible ***'
     3084                    ifPrint = True
     3085                break                   #refinement succeeded - finish up!
     3086            except TypeError:          #result[1] is None on singular matrix
     3087                print '**** Refinement failed - singular matrix ****'
     3088                if 'Hessian' in Controls['deriv type']:
     3089                    for i in result[2]['psing'].reverse():
     3090                            print 'Removing parameter: ',varyList[i]
     3091                            del(varyList[i])                   
     3092                else:
     3093                    Ipvt = result[2]['ipvt']
     3094                    for i,ipvt in enumerate(Ipvt):
     3095                        if not np.sum(result[2]['fjac'],axis=1)[i]:
     3096                            print 'Removing parameter: ',varyList[ipvt-1]
     3097                            del(varyList[ipvt-1])
     3098                            break
     3099   
     3100        GetFobsSq(Histo,Phases,parmDict,calcControls)
     3101        sigDict = dict(zip(varyList,sig))
     3102        newCellDict = GetNewCellParms(parmDict,varyList)
     3103        newAtomDict = ApplyXYZshifts(parmDict,varyList)
     3104        covData = {'variables':result[0],'varyList':varyList,'sig':sig,
     3105            'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
     3106        SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint)
     3107        SetHistogramData(parmDict,sigDict,Histo,ifPrint)
     3108        SeqResult[histogram] = covData
     3109        SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,covData,makeBack)
     3110        makeBack = False
     3111    SetSeqResult(GPXfile,Histograms,SeqResult)
     3112
     3113def DistAngle(DisAglCtls,DisAglData):
     3114    import numpy.ma as ma
     3115   
     3116    def ShowBanner(name):
     3117        print 80*'*'
     3118        print '   Interatomic Distances and Angles for phase '+name
     3119        print 80*'*','\n'
     3120
     3121    ShowBanner(DisAglCtls['Name'])
     3122    SGData = DisAglData['SGData']
     3123    SGtext = G2spc.SGPrint(SGData)
     3124    for line in SGtext: print line
     3125    Cell = DisAglData['Cell']
     3126   
     3127    Amat,Bmat = G2lat.cell2AB(Cell[:6])
     3128    covData = {}
     3129    if 'covData' in DisAglData:   
     3130        covData = DisAglData['covData']
     3131        covMatrix = covData['covMatrix']
     3132        varyList = covData['varyList']
     3133        pfx = str(DisAglData['pId'])+'::'
     3134        A = G2lat.cell2A(Cell[:6])
     3135        cellSig = getCellEsd(pfx,SGData,A,covData)
     3136        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
     3137        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
     3138        line = '\n Unit cell:'
     3139        for name,vals in zip(names,valEsd):
     3140            line += name+vals 
     3141        print line
     3142    else:
     3143        print '\n Unit cell: a = ','%.5f'%(Cell[0]),' b = ','%.5f'%(Cell[1]),' c = ','%.5f'%(Cell[2]), \
     3144            ' alpha = ','%.3f'%(Cell[3]),' beta = ','%.3f'%(Cell[4]),' gamma = ', \
     3145            '%.3f'%(Cell[5]),' volume = ','%.3f'%(Cell[6])
     3146    Factor = DisAglCtls['Factors']
     3147    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
     3148    Units = np.array([                   #is there a nicer way to make this?
     3149        [-1,-1,-1],[-1,-1,0],[-1,-1,1],[-1,0,-1],[-1,0,0],[-1,0,1],[-1,1,-1],[-1,1,0],[-1,1,1],
     3150        [0,-1,-1],[0,-1,0],[0,-1,1],[0,0,-1],[0,0,0],[0,0,1],[0,1,-1],[0,1,0],[0,1,1],
     3151        [1,-1,-1],[1,-1,0],[1,-1,1],[1,0,-1],[1,0,0],[1,0,1],[1,1,-1],[1,1,0],[1,1,1]])
     3152    origAtoms = DisAglData['OrigAtoms']
     3153    targAtoms = DisAglData['TargAtoms']
     3154    for Oatom in origAtoms:
     3155        OxyzNames = ''
     3156        IndBlist = []
     3157        Dist = []
     3158        Vect = []
     3159        VectA = []
     3160        angles = []
     3161        for Tatom in targAtoms:
     3162            Xvcov = []
     3163            TxyzNames = ''
     3164            if 'covData' in DisAglData:
     3165                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
     3166                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
     3167                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
     3168            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
     3169            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
     3170            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
     3171            for Txyz,Top,Tunit in result:
     3172                Dx = (Txyz-np.array(Oatom[3:6]))+Units
     3173                dx = np.inner(Amat,Dx)
     3174                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
     3175                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
     3176                if np.any(IndB):
     3177                    for indb in IndB:
     3178                        for i in range(len(indb)):
     3179                            if str(dx.T[indb][i]) not in IndBlist:
     3180                                IndBlist.append(str(dx.T[indb][i]))
     3181                                unit = Units[indb][i]
     3182                                tunit = '[%2d%2d%2d]'%(unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
     3183                                pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
     3184                                sig = 0.0
     3185                                if len(Xvcov):
     3186                                    sig = np.sqrt(np.inner(pdpx,np.inner(Xvcov,pdpx)))
     3187                                Dist.append([Oatom[1],Tatom[1],tunit,Top,ma.getdata(dist[indb])[i],sig])
     3188                                if (Dist[-1][-1]-AsumR) <= 0.:
     3189                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
     3190                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
     3191                                else:
     3192                                    Vect.append([0.,0.,0.])
     3193                                    VectA.append([])
     3194        Vect = np.array(Vect)
     3195        angles = np.zeros((len(Vect),len(Vect)))
     3196        angsig = np.zeros((len(Vect),len(Vect)))
     3197        for i,veca in enumerate(Vect):
     3198            if np.any(veca):
     3199                for j,vecb in enumerate(Vect):
     3200                    if np.any(vecb):
     3201                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
     3202        line = ''
     3203        for i,x in enumerate(Oatom[3:6]):
     3204            if len(Xvcov):
     3205                line += '%12s'%(G2mth.ValEsd(x,np.sqrt(Xvcov[i][i]),True))
     3206            else:
     3207                line += '%12.5f'%(x)
     3208        print '\n Distances & angles for ',Oatom[1],' at ',line
     3209        print 80*'*'
     3210        line = ''
     3211        for dist in Dist[:-1]:
     3212            line += '%12s'%(dist[1].center(12))
     3213        print '  To       cell +(sym. op.)      dist.  ',line
     3214        for i,dist in enumerate(Dist):
     3215            line = ''
     3216            for j,angle in enumerate(angles[i][0:i]):
     3217                sig = angsig[i][j]
     3218                if angle:
     3219                    if sig:
     3220                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
     3221                    else:
     3222                        val = '%.3f'%(angle)
     3223                        line += '%12s'%(val.center(12))
     3224                else:
     3225                    line += 12*' '
     3226            if dist[5]:            #sig exists!
     3227                val = G2mth.ValEsd(dist[4],dist[5])
     3228            else:
     3229                val = '%8.4f'%(dist[4])
     3230            print '  %8s%10s+(%4d) %12s'%(dist[1].ljust(8),dist[2].ljust(10),dist[3],val.center(12)),line
     3231
     3232def DisAglTor(DATData):
     3233    SGData = DATData['SGData']
     3234    Cell = DATData['Cell']
     3235   
     3236    Amat,Bmat = G2lat.cell2AB(Cell[:6])
     3237    covData = {}
     3238    pfx = ''
     3239    if 'covData' in DATData:   
     3240        covData = DATData['covData']
     3241        covMatrix = covData['covMatrix']
     3242        varyList = covData['varyList']
     3243        pfx = str(DATData['pId'])+'::'
     3244    Datoms = []
     3245    Oatoms = []
     3246    for i,atom in enumerate(DATData['Datoms']):
     3247        symop = atom[-1].split('+')
     3248        if len(symop) == 1:
     3249            symop.append('0,0,0')       
     3250        symop[0] = int(symop[0])
     3251        symop[1] = eval(symop[1])
     3252        atom.append(symop)
     3253        Datoms.append(atom)
     3254        oatom = DATData['Oatoms'][i]
     3255        names = ['','','']
     3256        if pfx:
     3257            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
     3258        oatom += [names,]
     3259        Oatoms.append(oatom)
     3260    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
     3261    if DATData['Natoms'] == 4:  #torsion
     3262        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
     3263        print ' Torsion angle for '+DATData['Name']+' atom sequence: ',atmSeq,'=',G2mth.ValEsd(Tors,sig)
     3264        print ' NB: Atom sequence determined by selection order'
     3265        return      # done with torsion
     3266    elif DATData['Natoms'] == 3:  #angle
     3267        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
     3268        print ' Angle in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Ang,sig)
     3269        print ' NB: Atom sequence determined by selection order'
     3270    else:   #2 atoms - distance
     3271        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
     3272        print ' Distance in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Dist,sig)
     3273               
     3274def BestPlane(PlaneData):
     3275
     3276    def ShowBanner(name):
     3277        print 80*'*'
     3278        print '   Best plane result for phase '+name
     3279        print 80*'*','\n'
     3280
     3281    ShowBanner(PlaneData['Name'])
     3282
     3283    Cell = PlaneData['Cell']   
     3284    Amat,Bmat = G2lat.cell2AB(Cell[:6])       
     3285    Atoms = PlaneData['Atoms']
     3286    sumXYZ = np.zeros(3)
     3287    XYZ = []
     3288    Natoms = len(Atoms)
     3289    for atom in Atoms:
     3290        xyz = np.array(atom[3:6])
     3291        XYZ.append(xyz)
     3292        sumXYZ += xyz
     3293    sumXYZ /= Natoms
     3294    XYZ = np.array(XYZ)-sumXYZ
     3295    XYZ = np.inner(Amat,XYZ).T
     3296    Zmat = np.zeros((3,3))
     3297    for i,xyz in enumerate(XYZ):
     3298        Zmat += np.outer(xyz.T,xyz)
     3299    print ' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2])
     3300    Evec,Emat = nl.eig(Zmat)
     3301    Evec = np.sqrt(Evec)/(Natoms-3)
     3302    Order = np.argsort(Evec)
     3303    XYZ = np.inner(XYZ,Emat.T).T
     3304    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
     3305    print ' Atoms in Cartesian best plane coordinates:'
     3306    print ' Name         X         Y         Z'
     3307    for i,xyz in enumerate(XYZ):
     3308        print ' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2])
     3309    print '\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]])   
    9563310
    9573311           
    958 
    959     def PrintFFtable(FFtable):
    960 
    961         print '\n X-ray scattering factors:'
    962 
    963         print '   Symbol     fa                                      fb                                      fc'
    964 
    965         print 99*'-'
    966 
    967         for Ename in FFtable:
    968 
    969             ffdata = FFtable[Ename]
    970 
    971             fa = ffdata['fa']
    972 
    973             fb = ffdata['fb']
    974 
    975             print ' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' %  \
    976 
    977                 (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
    978 
    979                
    980 
    981     def PrintBLtable(BLtable):
    982 
    983         print '\n Neutron scattering factors:'
    984 
    985         print '   Symbol   isotope       mass       b       resonant terms'
    986 
    987         print 99*'-'
    988 
    989         for Ename in BLtable:
    990 
    991             bldata = BLtable[Ename]
    992 
    993             isotope = bldata[0]
    994 
    995             mass = bldata[1][0]
    996 
    997             blen = bldata[1][1]
    998 
    999             bres = []
    1000 
    1001             if len(bldata[1]) > 2:
    1002 
    1003                 bres = bldata[1][2:]
    1004 
    1005             line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
    1006 
    1007             for item in bres:
    1008 
    1009                 line += '%10.5g'%(item)
    1010 
    1011             print line
    1012 
    1013                
    1014 
    1015     def PrintAtoms(General,Atoms):
    1016 
    1017         print '\n Atoms:'
    1018 
    1019         line = '   name    type  refine?   x         y         z    '+ \
    1020 
    1021             '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
    1022 
    1023         if General['Type'] == 'magnetic':
    1024 
    1025             line += '   Mx     My     Mz'
    1026 
    1027         elif General['Type'] == 'macromolecular':
    1028 
    1029             line = ' res no  residue  chain '+line
    1030 
    1031         print line
    1032 
    1033         if General['Type'] == 'nuclear':
    1034 
    1035             print 135*'-'
    1036 
    1037             for i,at in enumerate(Atoms):
    1038 
    1039                 line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \
    1040 
    1041                     '%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9])
    1042 
    1043                 if at[9] == 'I':
    1044 
    1045                     line += '%8.4f'%(at[10])+48*' '
    1046 
    1047                 else:
    1048 
    1049                     line += 8*' '
    1050 
    1051                     for j in range(6):
    1052 
    1053                         line += '%8.4f'%(at[11+j])
    1054 
    1055                 print line
    1056 
    1057        
    1058 
    1059     def PrintTexture(textureData):
    1060 
    1061         topstr = '\n Spherical harmonics texture: Order:' + \
    1062 
    1063             str(textureData['Order'])
    1064 
    1065         if textureData['Order']:
    1066 
    1067             print topstr+' Refine? '+str(textureData['SH Coeff'][0])
    1068 
    1069         else:
    1070 
    1071             print topstr
    1072 
    1073             return
    1074 
    1075         names = ['omega','chi','phi']
    1076 
    1077         line = '\n'
    1078 
    1079         for name in names:
    1080 
    1081             line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
    1082 
    1083         print line
    1084 
    1085         print '\n Texture coefficients:'
    1086 
    1087         ptlbls = ' names :'
    1088 
    1089         ptstr =  ' values:'
    1090 
    1091         SHcoeff = textureData['SH Coeff'][1]
    1092 
    1093         for item in SHcoeff:
    1094 
    1095             ptlbls += '%12s'%(item)
    1096 
    1097             ptstr += '%12.4f'%(SHcoeff[item])
    1098 
    1099         print ptlbls
    1100 
    1101         print ptstr   
    1102 
    1103        
    1104 
    1105     if Print: print ' Phases:'
    1106 
    1107     phaseVary = []
    1108 
    1109     phaseDict = {}
    1110 
    1111     phaseConstr = {}
    1112 
    1113     pawleyLookup = {}
    1114 
    1115     FFtables = {}                   #scattering factors - xrays
    1116 
    1117     BLtables = {}                   # neutrons
    1118 
    1119     Natoms = {}
    1120 
    1121     AtMults = {}
    1122 
    1123     AtIA = {}
    1124 
    1125     shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
    1126 
    1127     SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
    1128 
    1129     for name in PhaseData:
    1130 
    1131         General = PhaseData[name]['General']
    1132 
    1133         pId = PhaseData[name]['pId']
    1134 
    1135         pfx = str(pId)+'::'
    1136 
    1137         FFtable = GetFFtable(General)
    1138 
    1139         BLtable = GetBLtable(General)
    1140 
    1141         FFtables.update(FFtable)
    1142 
    1143         BLtables.update(BLtable)
    1144 
    1145         Atoms = PhaseData[name]['Atoms']
    1146 
    1147         try:
    1148 
    1149             PawleyRef = PhaseData[name]['Pawley ref']
    1150 
    1151         except KeyError:
    1152 
    1153             PawleyRef = []
    1154 
    1155         SGData = General['SGData']
    1156 
    1157         SGtext = G2spc.SGPrint(SGData)
    1158 
    1159         cell = General['Cell']
    1160 
    1161         A = G2lat.cell2A(cell[1:7])
    1162 
    1163         phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5]})
    1164 
    1165         if cell[0]:
    1166 
    1167             phaseVary += cellVary(pfx,SGData)
    1168 
    1169         Natoms[pfx] = 0
    1170 
    1171         if Atoms:
    1172 
    1173             if General['Type'] == 'nuclear':
    1174 
    1175                 Natoms[pfx] = len(Atoms)
    1176 
    1177                 for i,at in enumerate(Atoms):
    1178 
    1179                     phaseDict.update({pfx+'Atype:'+str(i):at[1],pfx+'Afrac:'+str(i):at[6],pfx+'Amul:'+str(i):at[8],
    1180 
    1181                         pfx+'Ax:'+str(i):at[3],pfx+'Ay:'+str(i):at[4],pfx+'Az:'+str(i):at[5],
    1182 
    1183                         pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0.,         #refined shifts for x,y,z
    1184 
    1185                         pfx+'AI/A:'+str(i):at[9],})
    1186 
    1187                     if at[9] == 'I':
    1188 
    1189                         phaseDict[pfx+'AUiso:'+str(i)] = at[10]
    1190 
    1191                     else:
    1192 
    1193                         phaseDict.update({pfx+'AU11:'+str(i):at[11],pfx+'AU22:'+str(i):at[12],pfx+'AU33:'+str(i):at[13],
    1194 
    1195                             pfx+'AU12:'+str(i):at[14],pfx+'AU13:'+str(i):at[15],pfx+'AU23:'+str(i):at[16]})
    1196 
    1197                     if 'F' in at[2]:
    1198 
    1199                         phaseVary.append(pfx+'Afrac:'+str(i))
    1200 
    1201                     if 'X' in at[2]:
    1202 
    1203                         xId,xCoef = G2spc.GetCSxinel(at[7])
    1204 
    1205                         delnames = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
    1206 
    1207                         for j in range(3):
    1208 
    1209                             if xId[j] > 0:                               
    1210 
    1211                                 phaseVary.append(delnames[j])
    1212 
    1213                                 for k in range(j):
    1214 
    1215                                     if xId[j] == xId[k]:
    1216 
    1217                                         G2mv.StoreEquivalence(delnames[k],((delnames[j],xCoef[j]),))
    1218 
    1219                     if 'U' in at[2]:
    1220 
    1221                         if at[9] == 'I':
    1222 
    1223                             phaseVary.append(pfx+'AUiso:'+str(i))
    1224 
    1225                         else:
    1226 
    1227                             uId,uCoef = G2spc.GetCSuinel(at[7])[:2]
    1228 
    1229                             names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
    1230 
    1231                                 pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
    1232 
    1233                             for j in range(6):
    1234 
    1235                                 if uId[j] > 0:                               
    1236 
    1237                                     phaseVary.append(names[j])
    1238 
    1239                                     for k in range(j):
    1240 
    1241                                         if uId[j] == uId[k]:
    1242 
    1243                                             G2mv.StoreEquivalence(names[k],((names[j],uCoef[j]),))
    1244 
    1245 #            elif General['Type'] == 'magnetic':
    1246 
    1247 #            elif General['Type'] == 'macromolecular':
    1248 
    1249 
    1250 
    1251                
    1252 
    1253             if 'SH Texture' in General:
    1254 
    1255                 textureData = General['SH Texture']
    1256 
    1257                 phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
    1258 
    1259                 phaseDict[pfx+'SHorder'] = textureData['Order']
    1260 
    1261                 for name in ['omega','chi','phi']:
    1262 
    1263                     phaseDict[pfx+'SH '+name] = textureData['Sample '+name][1]
    1264 
    1265                     if textureData['Sample '+name][0]:
    1266 
    1267                         phaseVary.append(pfx+'SH '+name)
    1268 
    1269                 for name in textureData['SH Coeff'][1]:
    1270 
    1271                     phaseDict[pfx+name] = textureData['SH Coeff'][1][name]
    1272 
    1273                     if textureData['SH Coeff'][0]:
    1274 
    1275                         phaseVary.append(pfx+name)
    1276 
    1277                
    1278 
    1279             if Print:
    1280 
    1281                 print '\n Phase name: ',General['Name']
    1282 
    1283                 print 135*'-'
    1284 
    1285                 PrintFFtable(FFtable)
    1286 
    1287                 PrintBLtable(BLtable)
    1288 
    1289                 print ''
    1290 
    1291                 for line in SGtext: print line
    1292 
    1293                 PrintAtoms(General,Atoms)
    1294 
    1295                 print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
    1296 
    1297                     ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
    1298 
    1299                     '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
    1300 
    1301                 if 'SH Texture' in General:
    1302 
    1303                     PrintTexture(textureData)
    1304 
    1305                    
    1306 
    1307         elif PawleyRef:
    1308 
    1309             pawleyVary = []
    1310 
    1311             for i,refl in enumerate(PawleyRef):
    1312 
    1313                 phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
    1314 
    1315                 pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
    1316 
    1317                 if refl[5]:
    1318 
    1319                     pawleyVary.append(pfx+'PWLref:'+str(i))
    1320 
    1321             GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary)      #does G2mv.StoreEquivalence
    1322 
    1323             phaseVary += pawleyVary
    1324 
    1325                
    1326 
    1327     return Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
    1328 
    1329    
    1330 
    1331 def cellFill(pfx,SGData,parmDict,sigDict):
    1332 
    1333     if SGData['SGLaue'] in ['-1',]:
    1334 
    1335         A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
    1336 
    1337             parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
    1338 
    1339         sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
    1340 
    1341             sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
    1342 
    1343     elif SGData['SGLaue'] in ['2/m',]:
    1344 
    1345         if SGData['SGUniq'] == 'a':
    1346 
    1347             A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
    1348 
    1349                 parmDict[pfx+'A3'],0,0]
    1350 
    1351             sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
    1352 
    1353                 sigDict[pfx+'A3'],0,0]
    1354 
    1355         elif SGData['SGUniq'] == 'b':
    1356 
    1357             A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
    1358 
    1359                 0,parmDict[pfx+'A4'],0]
    1360 
    1361             sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
    1362 
    1363                 0,sigDict[pfx+'A4'],0]
    1364 
    1365         else:
    1366 
    1367             A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
    1368 
    1369                 0,0,parmDict[pfx+'A5']]
    1370 
    1371             sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
    1372 
    1373                 0,0,sigDict[pfx+'A5']]
    1374 
    1375     elif SGData['SGLaue'] in ['mmm',]:
    1376 
    1377         A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
    1378 
    1379         sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
    1380 
    1381     elif SGData['SGLaue'] in ['4/m','4/mmm']:
    1382 
    1383         A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
    1384 
    1385         sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
    1386 
    1387     elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
    1388 
    1389         A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
    1390 
    1391             parmDict[pfx+'A0'],0,0]
    1392 
    1393         sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
    1394 
    1395     elif SGData['SGLaue'] in ['3R', '3mR']:
    1396 
    1397         A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
    1398 
    1399             parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
    1400 
    1401         sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
    1402 
    1403     elif SGData['SGLaue'] in ['m3m','m3']:
    1404 
    1405         A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
    1406 
    1407         sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
    1408 
    1409     return A,sigA
    1410 
    1411        
    1412 
    1413 def getCellEsd(pfx,SGData,A,covData):
    1414 
    1415     dpr = 180./np.pi
    1416 
    1417     rVsq = G2lat.calc_rVsq(A)
    1418 
    1419     G,g = G2lat.A2Gmat(A)       #get recip. & real metric tensors
    1420 
    1421     cell = np.array(G2lat.Gmat2cell(g))   #real cell
    1422 
    1423     cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
    1424 
    1425     scos = cosd(cellst[3:6])
    1426 
    1427     ssin = sind(cellst[3:6])
    1428 
    1429     scot = scos/ssin
    1430 
    1431     rcos = cosd(cell[3:6])
    1432 
    1433     rsin = sind(cell[3:6])
    1434 
    1435     rcot = rcos/rsin
    1436 
    1437     RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
    1438 
    1439     varyList = covData['varyList']
    1440 
    1441     covMatrix = covData['covMatrix']
    1442 
    1443     vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
    1444 
    1445     Ax = np.array(A)
    1446 
    1447     Ax[3:] /= 2.
    1448 
    1449     drVdA = np.array([Ax[1]*Ax[2]-Ax[5]**2,Ax[0]*Ax[2]-Ax[4]**2,Ax[0]*Ax[1]-Ax[3]**2,
    1450 
    1451         Ax[4]*Ax[5]-Ax[2]*Ax[3],Ax[3]*Ax[5]-Ax[1]*Ax[4],Ax[3]*Ax[4]-Ax[0]*Ax[5]])
    1452 
    1453     srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
    1454 
    1455     Vol = 1/np.sqrt(rVsq)
    1456 
    1457     sigVol = Vol**3*np.sqrt(srcvlsq)/2.
    1458 
    1459     R123 = Ax[0]*Ax[1]*Ax[2]
    1460 
    1461     dsasdg = np.zeros((3,6))
    1462 
    1463     dadg = np.zeros((6,6))
    1464 
    1465     for i0 in range(3):         #0  1   2
    1466 
    1467         i1 = (i0+1)%3           #1  2   0
    1468 
    1469         i2 = (i1+1)%3           #2  0   1
    1470 
    1471         i3 = 5-i2               #3  5   4
    1472 
    1473         i4 = 5-i1               #4  3   5
    1474 
    1475         i5 = 5-i0               #5  4   3
    1476 
    1477         dsasdg[i0][i1] = 0.5*scot[i0]*scos[i0]/Ax[i1]
    1478 
    1479         dsasdg[i0][i2] = 0.5*scot[i0]*scos[i0]/Ax[i2]
    1480 
    1481         dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
    1482 
    1483         denmsq = Ax[i0]*(R123-Ax[i1]*Ax[i4]**2-Ax[i2]*Ax[i3]**2+(Ax[i4]*Ax[i3])**2)
    1484 
    1485         denom = np.sqrt(denmsq)
    1486 
    1487         dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq*(R123-0.5*Ax[i1]*Ax[i4]**2-0.5*Ax[i2]*Ax[i3]**2)
    1488 
    1489         dadg[i5][i1] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i2]-Ax[i0]*Ax[i4]**2)
    1490 
    1491         dadg[i5][i2] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i1]-Ax[i0]*Ax[i3]**2)
    1492 
    1493         dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i2]*Ax[i3]-Ax[i3]*Ax[i4]**2)
    1494 
    1495         dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i1]*Ax[i4]-Ax[i3]**2*Ax[i4])
    1496 
    1497         dadg[i5][i5] = -Ax[i0]/denom
    1498 
    1499     for i0 in range(3):
    1500 
    1501         i1 = (i0+1)%3
    1502 
    1503         i2 = (i1+1)%3
    1504 
    1505         i3 = 5-i2
    1506 
    1507         for ij in range(6):
    1508 
    1509             dadg[i0][ij] = cell[i0]*(rcot[i2]*dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
    1510 
    1511             if ij == i0:
    1512 
    1513                 dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
    1514 
    1515             dadg[i3][ij] = -dadg[i3][ij]*rsin[2-i0]*dpr
    1516 
    1517     sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
    1518 
    1519     var = np.diag(sigMat)
    1520 
    1521     CS = np.where(var>0.,np.sqrt(var),0.)
    1522 
    1523     cellSig = [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]  #exchange sig(alp) & sig(gam) to get in right order
    1524 
    1525     return cellSig           
    1526 
    1527    
    1528 
    1529 def SetPhaseData(parmDict,sigDict,Phases,covData):
    1530 
    1531    
    1532 
    1533     def PrintAtomsAndSig(General,Atoms,atomsSig):
    1534 
    1535         print '\n Atoms:'
    1536 
    1537         line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
    1538 
    1539         if General['Type'] == 'magnetic':
    1540 
    1541             line += '   Mx     My     Mz'
    1542 
    1543         elif General['Type'] == 'macromolecular':
    1544 
    1545             line = ' res no  residue  chain '+line
    1546 
    1547         print line
    1548 
    1549         if General['Type'] == 'nuclear':
    1550 
    1551             print 135*'-'
    1552 
    1553             fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
    1554 
    1555                 11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
    1556 
    1557             noFXsig = {3:[10*' ','%10s'],4:[10*' ','%10s'],5:[10*' ','%10s'],6:[8*' ','%8s']}
    1558 
    1559             for i,at in enumerate(Atoms):
    1560 
    1561                 name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
    1562 
    1563                 valstr = ' values:'
    1564 
    1565                 sigstr = ' sig   :'
    1566 
    1567                 for ind in [3,4,5,6]:
    1568 
    1569                     sigind = str(i)+':'+str(ind)
    1570 
    1571                     valstr += fmt[ind]%(at[ind])                   
    1572 
    1573                     if sigind in atomsSig:
    1574 
    1575                         sigstr += fmt[ind]%(atomsSig[sigind])
    1576 
    1577                     else:
    1578 
    1579                         sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
    1580 
    1581                 if at[9] == 'I':
    1582 
    1583                     valstr += fmt[10]%(at[10])
    1584 
    1585                     if str(i)+':10' in atomsSig:
    1586 
    1587                         sigstr += fmt[10]%(atomsSig[str(i)+':10'])
    1588 
    1589                     else:
    1590 
    1591                         sigstr += 8*' '
    1592 
    1593                 else:
    1594 
    1595                     valstr += 8*' '
    1596 
    1597                     sigstr += 8*' '
    1598 
    1599                     for ind in [11,12,13,14,15,16]:
    1600 
    1601                         sigind = str(i)+':'+str(ind)
    1602 
    1603                         valstr += fmt[ind]%(at[ind])
    1604 
    1605                         if sigind in atomsSig:                       
    1606 
    1607                             sigstr += fmt[ind]%(atomsSig[sigind])
    1608 
    1609                         else:
    1610 
    1611                             sigstr += 8*' '
    1612 
    1613                 print name
    1614 
    1615                 print valstr
    1616 
    1617                 print sigstr
    1618 
    1619                
    1620 
    1621     def PrintSHtextureAndSig(textureData,SHtextureSig):
    1622 
    1623         print '\n Spherical harmonics texture: Order:' + str(textureData['Order'])
    1624 
    1625         names = ['omega','chi','phi']
    1626 
    1627         namstr = '  names :'
    1628 
    1629         ptstr =  '  values:'
    1630 
    1631         sigstr = '  esds  :'
    1632 
    1633         for name in names:
    1634 
    1635             namstr += '%12s'%(name)
    1636 
    1637             ptstr += '%12.3f'%(textureData['Sample '+name][1])
    1638 
    1639             if 'Sample '+name in SHtextureSig:
    1640 
    1641                 sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
    1642 
    1643             else:
    1644 
    1645                 sigstr += 12*' '
    1646 
    1647         print namstr
    1648 
    1649         print ptstr
    1650 
    1651         print sigstr
    1652 
    1653         print '\n Texture coefficients:'
    1654 
    1655         namstr = '  names :'
    1656 
    1657         ptstr =  '  values:'
    1658 
    1659         sigstr = '  esds  :'
    1660 
    1661         SHcoeff = textureData['SH Coeff'][1]
    1662 
    1663         for name in SHcoeff:
    1664 
    1665             namstr += '%12s'%(name)
    1666 
    1667             ptstr += '%12.3f'%(SHcoeff[name])
    1668 
    1669             if name in SHtextureSig:
    1670 
    1671                 sigstr += '%12.3f'%(SHtextureSig[name])
    1672 
    1673             else:
    1674 
    1675                 sigstr += 12*' '
    1676 
    1677         print namstr
    1678 
    1679         print ptstr
    1680 
    1681         print sigstr
    1682 
    1683        
    1684 
    1685            
    1686 
    1687     print '\n Phases:'
    1688 
    1689     for phase in Phases:
    1690 
    1691         print ' Result for phase: ',phase
    1692 
    1693         Phase = Phases[phase]
    1694 
    1695         General = Phase['General']
    1696 
    1697         SGData = General['SGData']
    1698 
    1699         Atoms = Phase['Atoms']
    1700 
    1701         cell = General['Cell']
    1702 
    1703         pId = Phase['pId']
    1704 
    1705         pfx = str(pId)+'::'
    1706 
    1707         if cell[0]:
    1708 
    1709             A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
    1710 
    1711             cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
    1712 
    1713             print ' Reciprocal metric tensor: '
    1714 
    1715             ptfmt = "%15.9f"
    1716 
    1717             names = ['A11','A22','A33','A12','A13','A23']
    1718 
    1719             namstr = '  names :'
    1720 
    1721             ptstr =  '  values:'
    1722 
    1723             sigstr = '  esds  :'
    1724 
    1725             for name,a,siga in zip(names,A,sigA):
    1726 
    1727                 namstr += '%15s'%(name)
    1728 
    1729                 ptstr += ptfmt%(a)
    1730 
    1731                 if siga:
    1732 
    1733                     sigstr += ptfmt%(siga)
    1734 
    1735                 else:
    1736 
    1737                     sigstr += 15*' '
    1738 
    1739             print namstr
    1740 
    1741             print ptstr
    1742 
    1743             print sigstr
    1744 
    1745             cell[1:7] = G2lat.A2cell(A)
    1746 
    1747             cell[7] = G2lat.calc_V(A)
    1748 
    1749             print ' New unit cell:'
    1750 
    1751             ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
    1752 
    1753             names = ['a','b','c','alpha','beta','gamma','Volume']
    1754 
    1755             namstr = '  names :'
    1756 
    1757             ptstr =  '  values:'
    1758 
    1759             sigstr = '  esds  :'
    1760 
    1761             for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
    1762 
    1763                 namstr += '%12s'%(name)
    1764 
    1765                 ptstr += fmt%(a)
    1766 
    1767                 if siga:
    1768 
    1769                     sigstr += fmt%(siga)
    1770 
    1771                 else:
    1772 
    1773                     sigstr += 12*' '
    1774 
    1775             print namstr
    1776 
    1777             print ptstr
    1778 
    1779             print sigstr
    1780 
    1781            
    1782 
    1783         if 'Pawley' in Phase['General']['Type']:
    1784 
    1785             pawleyRef = Phase['Pawley ref']
    1786 
    1787             for i,refl in enumerate(pawleyRef):
    1788 
    1789                 key = pfx+'PWLref:'+str(i)
    1790 
    1791                 refl[6] = abs(parmDict[key])        #suppress negative Fsq
    1792 
    1793                 if key in sigDict:
    1794 
    1795                     refl[7] = sigDict[key]
    1796 
    1797                 else:
    1798 
    1799                     refl[7] = 0
    1800 
    1801         else:
    1802 
    1803             atomsSig = {}
    1804 
    1805             if General['Type'] == 'nuclear':
    1806 
    1807                 for i,at in enumerate(Atoms):
    1808 
    1809                     names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
    1810 
    1811                         10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
    1812 
    1813                         14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
    1814 
    1815                     for ind in [3,4,5,6]:
    1816 
    1817                         at[ind] = parmDict[names[ind]]
    1818 
    1819                         if ind in [3,4,5]:
    1820 
    1821                             name = names[ind].replace('A','dA')
    1822 
    1823                         else:
    1824 
    1825                             name = names[ind]
    1826 
    1827                         if name in sigDict:
    1828 
    1829                             atomsSig[str(i)+':'+str(ind)] = sigDict[name]
    1830 
    1831                     if at[9] == 'I':
    1832 
    1833                         at[10] = parmDict[names[10]]
    1834 
    1835                         if names[10] in sigDict:
    1836 
    1837                             atomsSig[str(i)+':10'] = sigDict[names[10]]
    1838 
    1839                     else:
    1840 
    1841                         for ind in [11,12,13,14,15,16]:
    1842 
    1843                             at[ind] = parmDict[names[ind]]
    1844 
    1845                             if names[ind] in sigDict:
    1846 
    1847                                 atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
    1848 
    1849             PrintAtomsAndSig(General,Atoms,atomsSig)
    1850 
    1851        
    1852 
    1853         if 'SH Texture' in General:
    1854 
    1855             textureData = General['SH Texture']   
    1856 
    1857             if textureData['Order']:
    1858 
    1859                 SHtextureSig = {}
    1860 
    1861                 for name in ['omega','chi','phi']:
    1862 
    1863                     aname = pfx+'SH '+name
    1864 
    1865                     textureData['Sample '+name][1] = parmDict[aname]
    1866 
    1867                     if aname in sigDict:
    1868 
    1869                         SHtextureSig['Sample '+name] = sigDict[aname]
    1870 
    1871                 for name in textureData['SH Coeff'][1]:
    1872 
    1873                     aname = pfx+name
    1874 
    1875                     textureData['SH Coeff'][1][name] = parmDict[aname]
    1876 
    1877                     if aname in sigDict:
    1878 
    1879                         SHtextureSig[name] = sigDict[aname]
    1880 
    1881                 PrintSHtextureAndSig(textureData,SHtextureSig)
    1882 
    1883 
    1884 
    1885 def GetHistogramPhaseData(Phases,Histograms,Print=True):
    1886 
    1887    
    1888 
    1889     def PrintSize(hapData):
    1890 
    1891         if hapData[0] in ['isotropic','uniaxial']:
    1892 
    1893             line = '\n Size model    : %9s'%(hapData[0])
    1894 
    1895             line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
    1896 
    1897             if hapData[0] == 'uniaxial':
    1898 
    1899                 line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
    1900 
    1901             print line
    1902 
    1903         else:
    1904 
    1905             print '\n Size model    : %s'%(hapData[0])
    1906 
    1907             Snames = ['S11','S22','S33','S12','S13','S23']
    1908 
    1909             ptlbls = ' names :'
    1910 
    1911             ptstr =  ' values:'
    1912 
    1913             varstr = ' refine:'
    1914 
    1915             for i,name in enumerate(Snames):
    1916 
    1917                 ptlbls += '%12s' % (name)
    1918 
    1919                 ptstr += '%12.6f' % (hapData[4][i])
    1920 
    1921                 varstr += '%12s' % (str(hapData[5][i]))
    1922 
    1923             print ptlbls
    1924 
    1925             print ptstr
    1926 
    1927             print varstr
    1928 
    1929        
    1930 
    1931     def PrintMuStrain(hapData,SGData):
    1932 
    1933         if hapData[0] in ['isotropic','uniaxial']:
    1934 
    1935             line = '\n Mustrain model: %9s'%(hapData[0])
    1936 
    1937             line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
    1938 
    1939             if hapData[0] == 'uniaxial':
    1940 
    1941                 line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
    1942 
    1943             print line
    1944 
    1945         else:
    1946 
    1947             print '\n Mustrain model: %s'%(hapData[0])
    1948 
    1949             Snames = G2spc.MustrainNames(SGData)
    1950 
    1951             ptlbls = ' names :'
    1952 
    1953             ptstr =  ' values:'
    1954 
    1955             varstr = ' refine:'
    1956 
    1957             for i,name in enumerate(Snames):
    1958 
    1959                 ptlbls += '%12s' % (name)
    1960 
    1961                 ptstr += '%12.6f' % (hapData[4][i])
    1962 
    1963                 varstr += '%12s' % (str(hapData[5][i]))
    1964 
    1965             print ptlbls
    1966 
    1967             print ptstr
    1968 
    1969             print varstr
    1970 
    1971 
    1972 
    1973     def PrintHStrain(hapData,SGData):
    1974 
    1975         print '\n Hydrostatic/elastic strain: '
    1976 
    1977         Hsnames = G2spc.HStrainNames(SGData)
    1978 
    1979         ptlbls = ' names :'
    1980 
    1981         ptstr =  ' values:'
    1982 
    1983         varstr = ' refine:'
    1984 
    1985         for i,name in enumerate(Hsnames):
    1986 
    1987             ptlbls += '%12s' % (name)
    1988 
    1989             ptstr += '%12.6f' % (hapData[0][i])
    1990 
    1991             varstr += '%12s' % (str(hapData[1][i]))
    1992 
    1993         print ptlbls
    1994 
    1995         print ptstr
    1996 
    1997         print varstr
    1998 
    1999 
    2000 
    2001     def PrintSHPO(hapData):
    2002 
    2003         print '\n Spherical harmonics preferred orientation: Order:' + \
    2004 
    2005             str(hapData[4])+' Refine? '+str(hapData[2])
    2006 
    2007         ptlbls = ' names :'
    2008 
    2009         ptstr =  ' values:'
    2010 
    2011         for item in hapData[5]:
    2012 
    2013             ptlbls += '%12s'%(item)
    2014 
    2015             ptstr += '%12.3f'%(hapData[5][item])
    2016 
    2017         print ptlbls
    2018 
    2019         print ptstr
    2020 
    2021    
    2022 
    2023     hapDict = {}
    2024 
    2025     hapVary = []
    2026 
    2027     controlDict = {}
    2028 
    2029     poType = {}
    2030 
    2031     poAxes = {}
    2032 
    2033     spAxes = {}
    2034 
    2035     spType = {}
    2036 
    2037    
    2038 
    2039     for phase in Phases:
    2040 
    2041         HistoPhase = Phases[phase]['Histograms']
    2042 
    2043         SGData = Phases[phase]['General']['SGData']
    2044 
    2045         cell = Phases[phase]['General']['Cell'][1:7]
    2046 
    2047         A = G2lat.cell2A(cell)
    2048 
    2049         pId = Phases[phase]['pId']
    2050 
    2051         histoList = HistoPhase.keys()
    2052 
    2053         histoList.sort()
    2054 
    2055         for histogram in histoList:
    2056 
    2057             try:
    2058 
    2059                 Histogram = Histograms[histogram]
    2060 
    2061             except KeyError:                       
    2062 
    2063                 #skip if histogram not included e.g. in a sequential refinement
    2064 
    2065                 continue
    2066 
    2067             hapData = HistoPhase[histogram]
    2068 
    2069             hId = Histogram['hId']
    2070 
    2071             limits = Histogram['Limits'][1]
    2072 
    2073             inst = Histogram['Instrument Parameters']
    2074 
    2075             inst = dict(zip(inst[3],inst[1]))
    2076 
    2077             Zero = inst['Zero']
    2078 
    2079             if 'C' in inst['Type']:
    2080 
    2081                 try:
    2082 
    2083                     wave = inst['Lam']
    2084 
    2085                 except KeyError:
    2086 
    2087                     wave = inst['Lam1']
    2088 
    2089                 dmin = wave/(2.0*sind(limits[1]/2.0))
    2090 
    2091             pfx = str(pId)+':'+str(hId)+':'
    2092 
    2093             for item in ['Scale','Extinction']:
    2094 
    2095                 hapDict[pfx+item] = hapData[item][0]
    2096 
    2097                 if hapData[item][1]:
    2098 
    2099                     hapVary.append(pfx+item)
    2100 
    2101             names = G2spc.HStrainNames(SGData)
    2102 
    2103             for i,name in enumerate(names):
    2104 
    2105                 hapDict[pfx+name] = hapData['HStrain'][0][i]
    2106 
    2107                 if hapData['HStrain'][1][i]:
    2108 
    2109                     hapVary.append(pfx+name)
    2110 
    2111             controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
    2112 
    2113             if hapData['Pref.Ori.'][0] == 'MD':
    2114 
    2115                 hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
    2116 
    2117                 controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
    2118 
    2119                 if hapData['Pref.Ori.'][2]:
    2120 
    2121                     hapVary.append(pfx+'MD')
    2122 
    2123             else:                           #'SH' spherical harmonics
    2124 
    2125                 controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
    2126 
    2127                 controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
    2128 
    2129                 for item in hapData['Pref.Ori.'][5]:
    2130 
    2131                     hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
    2132 
    2133                     if hapData['Pref.Ori.'][2]:
    2134 
    2135                         hapVary.append(pfx+item)
    2136 
    2137             for item in ['Mustrain','Size']:
    2138 
    2139                 controlDict[pfx+item+'Type'] = hapData[item][0]
    2140 
    2141                 if hapData[item][0] in ['isotropic','uniaxial']:
    2142 
    2143                     hapDict[pfx+item+':i'] = hapData[item][1][0]
    2144 
    2145                     if hapData[item][2][0]:
    2146 
    2147                         hapVary.append(pfx+item+':i')
    2148 
    2149                     if hapData[item][0] == 'uniaxial':
    2150 
    2151                         controlDict[pfx+item+'Axis'] = hapData[item][3]
    2152 
    2153                         hapDict[pfx+item+':a'] = hapData[item][1][1]
    2154 
    2155                         if hapData[item][2][1]:
    2156 
    2157                             hapVary.append(pfx+item+':a')
    2158 
    2159                 else:       #generalized for mustrain or ellipsoidal for size
    2160 
    2161                     if item == 'Mustrain':
    2162 
    2163                         names = G2spc.MustrainNames(SGData)
    2164 
    2165                         pwrs = []
    2166 
    2167                         for name in names:
    2168 
    2169                             h,k,l = name[1:]
    2170 
    2171                             pwrs.append([int(h),int(k),int(l)])
    2172 
    2173                         controlDict[pfx+'MuPwrs'] = pwrs
    2174 
    2175                     for i in range(len(hapData[item][4])):
    2176 
    2177                         sfx = ':'+str(i)
    2178 
    2179                         hapDict[pfx+item+sfx] = hapData[item][4][i]
    2180 
    2181                         if hapData[item][5][i]:
    2182 
    2183                             hapVary.append(pfx+item+sfx)
    2184 
    2185                            
    2186 
    2187             if Print:
    2188 
    2189                 print '\n Phase: ',phase,' in histogram: ',histogram
    2190 
    2191                 print 135*'-'
    2192 
    2193                 print ' Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
    2194 
    2195                 print ' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
    2196 
    2197                 if hapData['Pref.Ori.'][0] == 'MD':
    2198 
    2199                     Ax = hapData['Pref.Ori.'][3]
    2200 
    2201                     print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
    2202 
    2203                         ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
    2204 
    2205                 else: #'SH' for spherical harmonics
    2206 
    2207                     PrintSHPO(hapData['Pref.Ori.'])
    2208 
    2209                 PrintSize(hapData['Size'])
    2210 
    2211                 PrintMuStrain(hapData['Mustrain'],SGData)
    2212 
    2213                 PrintHStrain(hapData['HStrain'],SGData)
    2214 
    2215             HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
    2216 
    2217             refList = []
    2218 
    2219             for h,k,l,d in HKLd:
    2220 
    2221                 ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
    2222 
    2223                 if ext:
    2224 
    2225                     continue
    2226 
    2227                 if 'C' in inst['Type']:
    2228 
    2229                     pos = 2.0*asind(wave/(2.0*d))+Zero
    2230 
    2231                     if limits[0] < pos < limits[1]:
    2232 
    2233                         refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0])
    2234 
    2235                 else:
    2236 
    2237                     raise ValueError
    2238 
    2239             Histogram['Reflection Lists'][phase] = refList
    2240 
    2241     return hapVary,hapDict,controlDict
    2242 
    2243    
    2244 
    2245 def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True):
    2246 
    2247    
    2248 
    2249     def PrintSizeAndSig(hapData,sizeSig):
    2250 
    2251         line = '\n Size model:     %9s'%(hapData[0])
    2252 
    2253         refine = False
    2254 
    2255         if hapData[0] in ['isotropic','uniaxial']:
    2256 
    2257             line += ' equatorial:%12.3f'%(hapData[1][0])
    2258 
    2259             if sizeSig[0][0]:
    2260 
    2261                 line += ', sig: %8.3f'%(sizeSig[0][0])
    2262 
    2263                 refine = True
    2264 
    2265             if hapData[0] == 'uniaxial':
    2266 
    2267                 line += ' axial:%12.3f'%(hapData[1][1])
    2268 
    2269                 if sizeSig[0][1]:
    2270 
    2271                     refine = True
    2272 
    2273                     line += ', sig: %8.3f'%(sizeSig[0][1])
    2274 
    2275             if refine:
    2276 
    2277                 print line
    2278 
    2279         else:
    2280 
    2281             Snames = ['S11','S22','S33','S12','S13','S23']
    2282 
    2283             ptlbls = ' name  :'
    2284 
    2285             ptstr =  ' value :'
    2286 
    2287             sigstr = ' sig   :'
    2288 
    2289             for i,name in enumerate(Snames):
    2290 
    2291                 ptlbls += '%12s' % (name)
    2292 
    2293                 ptstr += '%12.6f' % (hapData[4][i])
    2294 
    2295                 if sizeSig[1][i]:
    2296 
    2297                     refine = True
    2298 
    2299                     sigstr += '%12.6f' % (sizeSig[1][i])
    2300 
    2301                 else:
    2302 
    2303                     sigstr += 12*' '
    2304 
    2305             if refine:
    2306 
    2307                 print line
    2308 
    2309                 print ptlbls
    2310 
    2311                 print ptstr
    2312 
    2313                 print sigstr
    2314 
    2315        
    2316 
    2317     def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
    2318 
    2319         line = '\n Mustrain model: %9s'%(hapData[0])
    2320 
    2321         refine = False
    2322 
    2323         if hapData[0] in ['isotropic','uniaxial']:
    2324 
    2325             line += ' equatorial:%12.1f'%(hapData[1][0])
    2326 
    2327             if mustrainSig[0][0]:
    2328 
    2329                 line += ', sig: %8.1f'%(mustrainSig[0][0])
    2330 
    2331                 refine = True
    2332 
    2333             if hapData[0] == 'uniaxial':
    2334 
    2335                 line += ' axial:%12.1f'%(hapData[1][1])
    2336 
    2337                 if mustrainSig[0][1]:
    2338 
    2339                      line += ', sig: %8.1f'%(mustrainSig[0][1])
    2340 
    2341             if refine:
    2342 
    2343                 print line
    2344 
    2345         else:
    2346 
    2347             Snames = G2spc.MustrainNames(SGData)
    2348 
    2349             ptlbls = ' name  :'
    2350 
    2351             ptstr =  ' value :'
    2352 
    2353             sigstr = ' sig   :'
    2354 
    2355             for i,name in enumerate(Snames):
    2356 
    2357                 ptlbls += '%12s' % (name)
    2358 
    2359                 ptstr += '%12.6f' % (hapData[4][i])
    2360 
    2361                 if mustrainSig[1][i]:
    2362 
    2363                     refine = True
    2364 
    2365                     sigstr += '%12.6f' % (mustrainSig[1][i])
    2366 
    2367                 else:
    2368 
    2369                     sigstr += 12*' '
    2370 
    2371             if refine:
    2372 
    2373                 print line
    2374 
    2375                 print ptlbls
    2376 
    2377                 print ptstr
    2378 
    2379                 print sigstr
    2380 
    2381            
    2382 
    2383     def PrintHStrainAndSig(hapData,strainSig,SGData):
    2384 
    2385         Hsnames = G2spc.HStrainNames(SGData)
    2386 
    2387         ptlbls = ' name  :'
    2388 
    2389         ptstr =  ' value :'
    2390 
    2391         sigstr = ' sig   :'
    2392 
    2393         refine = False
    2394 
    2395         for i,name in enumerate(Hsnames):
    2396 
    2397             ptlbls += '%12s' % (name)
    2398 
    2399             ptstr += '%12.6g' % (hapData[0][i])
    2400 
    2401             if name in strainSig:
    2402 
    2403                 refine = True
    2404 
    2405                 sigstr += '%12.6g' % (strainSig[name])
    2406 
    2407             else:
    2408 
    2409                 sigstr += 12*' '
    2410 
    2411         if refine:
    2412 
    2413             print '\n Hydrostatic/elastic strain: '
    2414 
    2415             print ptlbls
    2416 
    2417             print ptstr
    2418 
    2419             print sigstr
    2420 
    2421        
    2422 
    2423     def PrintSHPOAndSig(hapData,POsig):
    2424 
    2425         print '\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
    2426 
    2427         ptlbls = ' names :'
    2428 
    2429         ptstr =  ' values:'
    2430 
    2431         sigstr = ' sig   :'
    2432 
    2433         for item in hapData[5]:
    2434 
    2435             ptlbls += '%12s'%(item)
    2436 
    2437             ptstr += '%12.3f'%(hapData[5][item])
    2438 
    2439             if item in POsig:
    2440 
    2441                 sigstr += '%12.3f'%(POsig[item])
    2442 
    2443             else:
    2444 
    2445                 sigstr += 12*' '
    2446 
    2447         print ptlbls
    2448 
    2449         print ptstr
    2450 
    2451         print sigstr
    2452 
    2453    
    2454 
    2455     for phase in Phases:
    2456 
    2457         HistoPhase = Phases[phase]['Histograms']
    2458 
    2459         SGData = Phases[phase]['General']['SGData']
    2460 
    2461         pId = Phases[phase]['pId']
    2462 
    2463         histoList = HistoPhase.keys()
    2464 
    2465         histoList.sort()
    2466 
    2467         for histogram in histoList:
    2468 
    2469             try:
    2470 
    2471                 Histogram = Histograms[histogram]
    2472 
    2473             except KeyError:                       
    2474 
    2475                 #skip if histogram not included e.g. in a sequential refinement
    2476 
    2477                 continue
    2478 
    2479             print '\n Phase: ',phase,' in histogram: ',histogram
    2480 
    2481             print 130*'-'
    2482 
    2483             hapData = HistoPhase[histogram]
    2484 
    2485             hId = Histogram['hId']
    2486 
    2487             pfx = str(pId)+':'+str(hId)+':'
    2488 
    2489             print ' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
    2490 
    2491                 %(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
    2492 
    2493            
    2494 
    2495             PhFrExtPOSig = {}
    2496 
    2497             for item in ['Scale','Extinction']:
    2498 
    2499                 hapData[item][0] = parmDict[pfx+item]
    2500 
    2501                 if pfx+item in sigDict:
    2502 
    2503                     PhFrExtPOSig[item] = sigDict[pfx+item]
    2504 
    2505             if hapData['Pref.Ori.'][0] == 'MD':
    2506 
    2507                 hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
    2508 
    2509                 if pfx+'MD' in sigDict:
    2510 
    2511                     PhFrExtPOSig['MD'] = sigDict[pfx+'MD']
    2512 
    2513             else:                           #'SH' spherical harmonics
    2514 
    2515                 for item in hapData['Pref.Ori.'][5]:
    2516 
    2517                     hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
    2518 
    2519                     if pfx+item in sigDict:
    2520 
    2521                         PhFrExtPOSig[item] = sigDict[pfx+item]
    2522 
    2523             if Print:
    2524 
    2525                 if 'Scale' in PhFrExtPOSig:
    2526 
    2527                     print ' Phase fraction  : %10.4f, sig %10.4f'%(hapData['Scale'][0],PhFrExtPOSig['Scale'])
    2528 
    2529                 if 'Extinction' in PhFrExtPOSig:
    2530 
    2531                     print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction'])
    2532 
    2533                 if hapData['Pref.Ori.'][0] == 'MD':
    2534 
    2535                     if 'MD' in PhFrExtPOSig:
    2536 
    2537                         print ' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD'])
    2538 
    2539                 else:
    2540 
    2541                     PrintSHPOAndSig(hapData['Pref.Ori.'],PhFrExtPOSig)
    2542 
    2543             SizeMuStrSig = {'Mustrain':[[0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
    2544 
    2545                 'Size':[[0,0],[0 for i in range(len(hapData['Size'][4]))]],
    2546 
    2547                 'HStrain':{}}                 
    2548 
    2549             for item in ['Mustrain','Size']:
    2550 
    2551                 if hapData[item][0] in ['isotropic','uniaxial']:                   
    2552 
    2553                     hapData[item][1][0] = parmDict[pfx+item+':i']
    2554 
    2555                     if item == 'Size':
    2556 
    2557                         hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
    2558 
    2559                     if pfx+item+':i' in sigDict:
    2560 
    2561                         SizeMuStrSig[item][0][0] = sigDict[pfx+item+':i']
    2562 
    2563                     if hapData[item][0] == 'uniaxial':
    2564 
    2565                         hapData[item][1][1] = parmDict[pfx+item+':a']
    2566 
    2567                         if item == 'Size':
    2568 
    2569                             hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
    2570 
    2571                         if pfx+item+':a' in sigDict:
    2572 
    2573                             SizeMuStrSig[item][0][1] = sigDict[pfx+item+':a']
    2574 
    2575                 else:       #generalized for mustrain or ellipsoidal for size
    2576 
    2577                     for i in range(len(hapData[item][4])):
    2578 
    2579                         sfx = ':'+str(i)
    2580 
    2581                         hapData[item][4][i] = parmDict[pfx+item+sfx]
    2582 
    2583                         if pfx+item+sfx in sigDict:
    2584 
    2585                             SizeMuStrSig[item][1][i] = sigDict[pfx+item+sfx]
    2586 
    2587             names = G2spc.HStrainNames(SGData)
    2588 
    2589             for i,name in enumerate(names):
    2590 
    2591                 hapData['HStrain'][0][i] = parmDict[pfx+name]
    2592 
    2593                 if pfx+name in sigDict:
    2594 
    2595                     SizeMuStrSig['HStrain'][name] = sigDict[pfx+name]
    2596 
    2597             if Print:
    2598 
    2599                 PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size'])
    2600 
    2601                 PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData)
    2602 
    2603                 PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig['HStrain'],SGData)
    2604 
    2605    
    2606 
    2607 def GetHistogramData(Histograms,Print=True):
    2608 
    2609    
    2610 
    2611     def GetBackgroundParms(hId,Background):
    2612 
    2613         Back = Background[0]
    2614 
    2615         Debye = Background[1]
    2616 
    2617         bakType,bakFlag = Back[:2]
    2618 
    2619         backVals = Back[3:]
    2620 
    2621         backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
    2622 
    2623         backDict = dict(zip(backNames,backVals))
    2624 
    2625         backVary = []
    2626 
    2627         if bakFlag:
    2628 
    2629             backVary = backNames
    2630 
    2631         backDict[':'+str(hId)+':nDebye'] = Debye['nDebye']
    2632 
    2633         debyeDict = {}
    2634 
    2635         debyeList = []
    2636 
    2637         for i in range(Debye['nDebye']):
    2638 
    2639             debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
    2640 
    2641             debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
    2642 
    2643             debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
    2644 
    2645         debyeVary = []
    2646 
    2647         for item in debyeList:
    2648 
    2649             if item[1]:
    2650 
    2651                 debyeVary.append(item[0])
    2652 
    2653         backDict.update(debyeDict)
    2654 
    2655         backVary += debyeVary   
    2656 
    2657         return bakType,backDict,backVary           
    2658 
    2659        
    2660 
    2661     def GetInstParms(hId,Inst):
    2662 
    2663         insVals,insFlags,insNames = Inst[1:4]
    2664 
    2665         dataType = insVals[0]
    2666 
    2667         instDict = {}
    2668 
    2669         insVary = []
    2670 
    2671         pfx = ':'+str(hId)+':'
    2672 
    2673         for i,flag in enumerate(insFlags):
    2674 
    2675             insName = pfx+insNames[i]
    2676 
    2677             instDict[insName] = insVals[i]
    2678 
    2679             if flag:
    2680 
    2681                 insVary.append(insName)
    2682 
    2683         instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001)
    2684 
    2685         instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001)
    2686 
    2687         instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
    2688 
    2689         return dataType,instDict,insVary
    2690 
    2691        
    2692 
    2693     def GetSampleParms(hId,Sample):
    2694 
    2695         sampVary = []
    2696 
    2697         hfx = ':'+str(hId)+':'       
    2698 
    2699         sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
    2700 
    2701             hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
    2702 
    2703         Type = Sample['Type']
    2704 
    2705         if 'Bragg' in Type:             #Bragg-Brentano
    2706 
    2707             for item in ['Scale','Shift','Transparency']:       #surface roughness?, diffuse scattering?
    2708 
    2709                 sampDict[hfx+item] = Sample[item][0]
    2710 
    2711                 if Sample[item][1]:
    2712 
    2713                     sampVary.append(hfx+item)
    2714 
    2715         elif 'Debye' in Type:        #Debye-Scherrer
    2716 
    2717             for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
    2718 
    2719                 sampDict[hfx+item] = Sample[item][0]
    2720 
    2721                 if Sample[item][1]:
    2722 
    2723                     sampVary.append(hfx+item)
    2724 
    2725         return Type,sampDict,sampVary
    2726 
    2727        
    2728 
    2729     def PrintBackground(Background):
    2730 
    2731         Back = Background[0]
    2732 
    2733         Debye = Background[1]
    2734 
    2735         print '\n Background function: ',Back[0],' Refine?',bool(Back[1])
    2736 
    2737         line = ' Coefficients: '
    2738 
    2739         for i,back in enumerate(Back[3:]):
    2740 
    2741             line += '%10.3f'%(back)
    2742 
    2743             if i and not i%10:
    2744 
    2745                 line += '\n'+15*' '
    2746 
    2747         print line
    2748 
    2749         if Debye['nDebye']:
    2750 
    2751             print '\n Debye diffuse scattering coefficients'
    2752 
    2753             parms = ['DebyeA','DebyeR','DebyeU']
    2754 
    2755             line = ' names :'
    2756 
    2757             for parm in parms:
    2758 
    2759                 line += '%16s'%(parm)
    2760 
    2761             print line
    2762 
    2763             for j,term in enumerate(Debye['debyeTerms']):
    2764 
    2765                 line = ' term'+'%2d'%(j)+':'
    2766 
    2767                 for i in range(3):
    2768 
    2769                     line += '%10.4g %5s'%(term[2*i],bool(term[2*i+1]))                   
    2770 
    2771                 print line
    2772 
    2773        
    2774 
    2775     def PrintInstParms(Inst):
    2776 
    2777         print '\n Instrument Parameters:'
    2778 
    2779         ptlbls = ' name  :'
    2780 
    2781         ptstr =  ' value :'
    2782 
    2783         varstr = ' refine:'
    2784 
    2785         instNames = Inst[3][1:]
    2786 
    2787         for i,name in enumerate(instNames):
    2788 
    2789             ptlbls += '%12s' % (name)
    2790 
    2791             ptstr += '%12.6f' % (Inst[1][i+1])
    2792 
    2793             if name in ['Lam1','Lam2','Azimuth']:
    2794 
    2795                 varstr += 12*' '
    2796 
    2797             else:
    2798 
    2799                 varstr += '%12s' % (str(bool(Inst[2][i+1])))
    2800 
    2801         print ptlbls
    2802 
    2803         print ptstr
    2804 
    2805         print varstr
    2806 
    2807        
    2808 
    2809     def PrintSampleParms(Sample):
    2810 
    2811         print '\n Sample Parameters:'
    2812 
    2813         print ' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
    2814 
    2815             (Sample['Omega'],Sample['Chi'],Sample['Phi'])
    2816 
    2817         ptlbls = ' name  :'
    2818 
    2819         ptstr =  ' value :'
    2820 
    2821         varstr = ' refine:'
    2822 
    2823         if 'Bragg' in Sample['Type']:
    2824 
    2825             for item in ['Scale','Shift','Transparency']:
    2826 
    2827                 ptlbls += '%14s'%(item)
    2828 
    2829                 ptstr += '%14.4f'%(Sample[item][0])
    2830 
    2831                 varstr += '%14s'%(str(bool(Sample[item][1])))
    2832 
    2833            
    2834 
    2835         elif 'Debye' in Type:        #Debye-Scherrer
    2836 
    2837             for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
    2838 
    2839                 ptlbls += '%14s'%(item)
    2840 
    2841                 ptstr += '%14.4f'%(Sample[item][0])
    2842 
    2843                 varstr += '%14s'%(str(bool(Sample[item][1])))
    2844 
    2845 
    2846 
    2847         print ptlbls
    2848 
    2849         print ptstr
    2850 
    2851         print varstr
    2852 
    2853        
    2854 
    2855 
    2856 
    2857     histDict = {}
    2858 
    2859     histVary = []
    2860 
    2861     controlDict = {}
    2862 
    2863     histoList = Histograms.keys()
    2864 
    2865     histoList.sort()
    2866 
    2867     for histogram in histoList:
    2868 
    2869         Histogram = Histograms[histogram]
    2870 
    2871         hId = Histogram['hId']
    2872 
    2873         pfx = ':'+str(hId)+':'
    2874 
    2875         controlDict[pfx+'Limits'] = Histogram['Limits'][1]
    2876 
    2877        
    2878 
    2879         Background = Histogram['Background']
    2880 
    2881         Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
    2882 
    2883         controlDict[pfx+'bakType'] = Type
    2884 
    2885         histDict.update(bakDict)
    2886 
    2887         histVary += bakVary
    2888 
    2889        
    2890 
    2891         Inst = Histogram['Instrument Parameters']
    2892 
    2893         Type,instDict,insVary = GetInstParms(hId,Inst)
    2894 
    2895         controlDict[pfx+'histType'] = Type
    2896 
    2897         if pfx+'Lam1' in instDict:
    2898 
    2899             controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
    2900 
    2901         else:
    2902 
    2903             controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
    2904 
    2905         histDict.update(instDict)
    2906 
    2907         histVary += insVary
    2908 
    2909        
    2910 
    2911         Sample = Histogram['Sample Parameters']
    2912 
    2913         Type,sampDict,sampVary = GetSampleParms(hId,Sample)
    2914 
    2915         controlDict[pfx+'instType'] = Type
    2916 
    2917         histDict.update(sampDict)
    2918 
    2919         histVary += sampVary
    2920 
    2921 
    2922 
    2923         if Print:
    2924 
    2925             print '\n Histogram: ',histogram,' histogram Id: ',hId
    2926 
    2927             print 135*'-'
    2928 
    2929             Units = {'C':' deg','T':' msec'}
    2930 
    2931             units = Units[controlDict[pfx+'histType'][2]]
    2932 
    2933             Limits = controlDict[pfx+'Limits']
    2934 
    2935             print ' Instrument type: ',Sample['Type']
    2936 
    2937             print ' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)     
    2938 
    2939             PrintSampleParms(Sample)
    2940 
    2941             PrintInstParms(Inst)
    2942 
    2943             PrintBackground(Background)
    2944 
    2945        
    2946 
    2947     return histVary,histDict,controlDict
    2948 
    2949    
    2950 
    2951 def SetHistogramData(parmDict,sigDict,Histograms,Print=True):
    2952 
    2953    
    2954 
    2955     def SetBackgroundParms(pfx,Background,parmDict,sigDict):
    2956 
    2957         Back = Background[0]
    2958 
    2959         Debye = Background[1]
    2960 
    2961         lenBack = len(Back[3:])
    2962 
    2963         backSig = [0 for i in range(lenBack+3*Debye['nDebye'])]
    2964 
    2965         for i in range(lenBack):
    2966 
    2967             Back[3+i] = parmDict[pfx+'Back:'+str(i)]
    2968 
    2969             if pfx+'Back:'+str(i) in sigDict:
    2970 
    2971                 backSig[i] = sigDict[pfx+'Back:'+str(i)]
    2972 
    2973         if Debye['nDebye']:
    2974 
    2975             for i in range(Debye['nDebye']):
    2976 
    2977                 names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
    2978 
    2979                 for j,name in enumerate(names):
    2980 
    2981                     Debye['debyeTerms'][i][2*j] = parmDict[name]
    2982 
    2983                     if name in sigDict: