Changeset 1957


Ignore:
Timestamp:
Aug 14, 2015 3:08:55 PM (8 years ago)
Author:
vondreele
Message:

work on SS structure factors
some refactoring of SS special pos code
rearrange sf routines in G2strMath - some math error in the SS sf codes
make SS plots of structure behave properly

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIImath.py

    r1956 r1957  
    924924################################################################################
    925925   
    926 def Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata):
     926def Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM):
    927927    import scipy.special as sp
    928928    import scipy.integrate as si
     929    Smult,TauT = SStauM             # both atoms x SGops
    929930    m = SSUniq.T[3]
    930931    nh = np.zeros(1)
    931932    if XSSdata.ndim > 2:
    932         nh = np.arange(XSSdata.shape[1])       
    933     M = np.where(m>0,m+nh[:,np.newaxis],m-nh[:,np.newaxis])
    934     A = np.array(XSSdata[:3])
    935     B = np.array(XSSdata[3:])
    936     HdotA = (np.inner(A.T,SSUniq.T[:3].T)+SSPhi)
    937     HdotB = (np.inner(B.T,SSUniq.T[:3].T)+SSPhi)
    938     GpA = sp.jn(M[:,np.newaxis],twopi*HdotA)
    939     GpB = sp.jn(M[:,np.newaxis],twopi*HdotB)*(1.j)**M
    940     Gp = np.sum(GpA+GpB,axis=0)
    941     return np.real(Gp).T,np.imag(Gp).T
    942    
    943 def ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata):
     933        nh = np.arange(XSSdata.shape[1])        #0..max harmonic order-1
     934    M = np.where(m>0,m+nh[:,np.newaxis],m-nh[:,np.newaxis])     # waves x SGops
     935    A = np.array(XSSdata[:3])   #sin mods x waves x atoms
     936    B = np.array(XSSdata[3:])   #cos mods...
     937    HdotA = (np.inner(A.T,SSUniq.T[:3].T)+SSPhi+TauT[:,np.newaxis,:])    # atoms X waves x SGops
     938    HdotB = (np.inner(B.T,SSUniq.T[:3].T)+SSPhi)    #+TauT[:,np.newaxis,:])    # ditto
     939    GpA = sp.jn(M,twopi*HdotA)                      # atoms x waves x SGops
     940    GpB = sp.jn(M,twopi*HdotB)*(1.j)**M             # ditto
     941    Gp = np.sum(GpA+GpB,axis=1)                     # atoms x SGops
     942    return np.real(Gp).T,np.imag(Gp).T              # SGops x atoms
     943   
     944def ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM):
    944945    import scipy.special as sp
    945946    m = SSUniq.T[3]
    … …  
    950951    A = np.array([[a,b] for a,b in zip(XSSdata[:3],XSSdata[3:])])
    951952    HdotA = twopi*(np.inner(SSUniq.T[:3].T,A.T)+SSPhi)
    952     Gpm = sp.jn(M[:,np.newaxis,:]-1,HdotA)
    953     Gpp = sp.jn(M[:,np.newaxis,:]+1,HdotA)
     953    Gpm = sp.jn(M-1,HdotA)
     954    Gpp = sp.jn(M+1,HdotA)
    954955    if Gpm.ndim > 3: #sum over multiple harmonics
    955956        Gpm = np.sum(Gpm,axis=0)
    … …  
    958959    return np.real(dGpdk),np.imag(dGpdk)
    959960   
    960 def posFourier(tau,psin,pcos):
    961     A = np.array([ps[:,np.newaxis]*np.sin(2*np.pi*(i+1)*tau) for i,ps in enumerate(psin)])
     961def posFourier(tau,psin,pcos,smul):
     962    A = np.array([ps[:,np.newaxis]*np.sin(2*np.pi*(i+1)*tau) for i,ps in enumerate(psin)])*smul
    962963    B = np.array([pc[:,np.newaxis]*np.cos(2*np.pi*(i+1)*tau) for i,pc in enumerate(pcos)])
    963964    return np.sum(A,axis=0)+np.sum(B,axis=0)
    … …  
    988989    generalData = data['General']
    989990    SGData = generalData['SGData']
     991    SSGData = generalData['SSGData']
    990992    cx,ct,cs,cia = generalData['AtomPtrs']
    991993    drawingData = data['Drawing']
    … …  
    9991001        Spos = atom[-1]['SS1']['Spos']
    10001002        Sadp = atom[-1]['SS1']['Sadp']
    1001         wave = np.zeros(3)
    1002         uwave = np.zeros(6)
    1003         if len(Spos):
    1004             scof = []
    1005             ccof = []
    1006             for i,spos in enumerate(Spos):
    1007                 if waveType in ['Sawtooth','ZigZag'] and not i:
    1008                     Toff = spos[0][0]
    1009                     slopes = np.array(spos[0][1:])
    1010                     if waveType == 'Sawtooth':
    1011                         wave = posSawtooth(tau,Toff,slopes)
    1012                     elif waveType == 'ZigZag':
    1013                         wave = posZigZag(tau,Toff,slopes)
    1014                 else:
    1015                     scof.append(spos[0][:3])
    1016                     ccof.append(spos[0][3:])
    1017             wave += np.sum(posFourier(tau,np.array(scof),np.array(ccof)),axis=1)
    1018         if len(Sadp):
    1019             scof = []
    1020             ccof = []
    1021             for i,sadp in enumerate(Sadp):
    1022                 scof.append(sadp[0][:6])
    1023                 ccof.append(sadp[0][6:])
    1024             uwave += np.sum(posFourier(tau,np.array(scof),np.array(ccof)),axis=1)
    10251003        indx = FindAtomIndexByIDs(drawAtoms,dci,[atom[cia+8],],True)
    10261004        for ind in indx:
    10271005            drawatom = drawAtoms[ind]
    10281006            opr = drawatom[dcs-1]
     1007            sop,ssop,icent = G2spc.OpsfromStringOps(opr,SGData,SSGData)
     1008            sdet,ssdet,dtau,dT,tauT = G2spc.getTauT(tau,sop,ssop,atxyz)
     1009            smul = sdet*ssdet
     1010#            tauT *= icent
     1011            wave = np.zeros(3)
     1012            uwave = np.zeros(6)
     1013            if len(Spos):
     1014                scof = []
     1015                ccof = []
     1016                for i,spos in enumerate(Spos):
     1017                    if waveType in ['Sawtooth','ZigZag'] and not i:
     1018                        Toff = spos[0][0]
     1019                        slopes = np.array(spos[0][1:])
     1020                        if waveType == 'Sawtooth':
     1021                            wave = posSawtooth(tauT,Toff,slopes)
     1022                        elif waveType == 'ZigZag':
     1023                            wave = posZigZag(tauT,Toff,slopes)
     1024                    else:
     1025                        scof.append(spos[0][:3])
     1026                        ccof.append(spos[0][3:])
     1027                wave += np.sum(posFourier(tauT,np.array(scof),np.array(ccof),smul),axis=1)
     1028            if len(Sadp):
     1029                scof = []
     1030                ccof = []
     1031                for i,sadp in enumerate(Sadp):
     1032                    scof.append(sadp[0][:6])
     1033                    ccof.append(sadp[0][6:])
     1034                uwave += np.sum(posFourier(tauT,np.array(scof),np.array(ccof),smul),axis=1)
    10291035            if atom[cia] == 'A':                   
    10301036                X,U = G2spc.ApplyStringOps(opr,SGData,atxyz+wave,atuij+uwave)
    … …  
    20762082                    h,k,l,m = -hkl+Hmax
    20772083                    Fhkl[h,k,l,m] = F*phasem
     2084                elif 'Fcalc' in mapData['MapType']:
     2085                    F = np.where(Fcsq>0.,np.sqrt(Fcsq),0.)
     2086                    h,k,l,m = hkl+Hmax
     2087                    Fhkl[h,k,l,m] = F*phasep
     2088                    h,k,l,m = -hkl+Hmax
     2089                    Fhkl[h,k,l,m] = F*phasem                   
    20782090                elif 'delt-F' in mapData['MapType']:
    20792091                    dF = np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)

  • trunk/GSASIIphsGUI.py

    r1953 r1957  
    24992499            mapSizer.Add(refList,0,WACV)
    25002500           
    2501             mapTypes = ['Fobs','delt-F']
     2501            mapTypes = ['Fobs','Fcalc','delt-F']
    25022502            mapSizer.Add(wx.StaticText(waveData,label=' Map type: '),0,WACV)
    25032503            mapType = wx.ComboBox(waveData,-1,value=Map['MapType'],choices=mapTypes,
    … …  
    25592559        mapSig = np.std(mapData['rho'])
    25602560        print mapData['MapType']+' computed: rhomax = %.3f rhomin = %.3f sigma = %.3f'%(np.max(mapData['rho']),np.min(mapData['rho']),mapSig)
     2561        wx.CallAfter(UpdateWavesData)
    25612562           
    25622563################################################################################

  • trunk/GSASIIplot.py

    r1951 r1957  
    30313031                scof.append(spos[0][:3])
    30323032                ccof.append(spos[0][3:])
    3033         wave += G2mth.posFourier(tau,np.array(scof),np.array(ccof))
     3033        wave += G2mth.posFourier(tau,np.array(scof),np.array(ccof),-1)  #hmm, why -1 work for Na2CO3?
    30343034    if mapData['Flip']:
    30353035        Title = 'Charge flip'

  • trunk/GSASIIspc.py

    r1956 r1957  
    903903    '''
    904904    modsym,gensym = SSymbol.replace(' ','').split(')')
     905    if gensym in ['0','00','000','0000']:       #get rid of extraneous symbols
     906        gensym = ''
    905907    nfrac = modsym.count('/')
    906908    modsym = modsym.lstrip('(')
    … …  
    14311433    return CSuinel[indx[1]]
    14321434   
     1435def getTauT(tau,sop,ssop,XYZ):
     1436    ssopinv = nl.inv(ssop[0])
     1437    mst = ssopinv[3][:3]
     1438    epsinv = ssopinv[3][3]
     1439    sdet = nl.det(sop[0])
     1440    ssdet = nl.det(ssop[0])
     1441    dtau = mst*(XYZ-sop[1])-epsinv*ssop[1][3]
     1442    dT = 1.0
     1443    if np.any(dtau%.5):
     1444        dT = np.tan(np.pi*np.sum(dtau%.5))
     1445    tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3])
     1446    return sdet,ssdet,dtau,dT,tauT
     1447   
     1448def OpsfromStringOps(A,SGData,SSGData):
     1449    SGOps = SGData['SGOps']
     1450    SSGOps = SSGData['SSGOps']
     1451    Ax = A.split('+')
     1452    Ax[0] = int(Ax[0])
     1453    iC = 1
     1454    if Ax[0] < 0:
     1455        iC = -1
     1456    Ax[0] = abs(Ax[0])
     1457    nA = Ax[0]%100-1
     1458    return SGOps[nA],SSGOps[nA],iC
     1459   
    14331460def GetSSfxuinel(waveType,nH,XYZ,SGData,SSGData,debug=False):
    14341461   
    … …  
    14431470                csi[i] = parms.index(csi[i])
    14441471        return CSI
    1445    
     1472       
    14461473    def fracCrenel(tau,Toff,Twid):
    14471474        Tau = (tau-Toff[:,np.newaxis])%1.
    … …  
    14861513        for i in SdIndx:
    14871514            sop = Sop[i]
    1488             ssop = SSop[i]
    1489             ssopinv = nl.inv(ssop[0])
    1490             mst = ssopinv[3][:3]
    1491             epsinv = ssopinv[3][3]
    1492             sdet = nl.det(sop[0])
    1493             ssdet = nl.det(ssop[0])
    1494             dtau = mst*(XYZ-sop[1])-epsinv*ssop[1][3]
    1495             dT = 1.0
    1496             if np.any(dtau%.5):
    1497                 dT = np.tan(np.pi*np.sum(dtau%.5))
    1498             tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3])
     1515            ssop = SSop[i]           
     1516            sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
    14991517            fsc = np.ones(2,dtype='i')
    15001518            if 'Crenel' in waveType:
    … …  
    15511569            sop = Sop[i]
    15521570            ssop = SSop[i]
     1571            sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
    15531572            xsc = np.ones(6,dtype='i')
    1554             ssopinv = nl.inv(ssop[0])
    1555             mst = ssopinv[3][:3]
    1556             epsinv = ssopinv[3][3]
    1557             sdet = nl.det(sop[0])
    1558             ssdet = nl.det(ssop[0])
    1559             dtau = mst*(XYZ-sop[1])-epsinv*ssop[1][3]
    1560             dT = 1.0
    1561             if np.any(dtau%.5):
    1562                 dT = np.tan(np.pi*np.sum(dtau%.5))
    1563             tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3])
    15641573            if waveType == 'Fourier':
    15651574                dXT = posFourier(np.sort(tauT),nH,delt6[:3],delt6[3:])   #+np.array(XYZ)[:,np.newaxis,np.newaxis]
    … …  
    16341643            sop = Sop[i]
    16351644            ssop = SSop[i]
    1636             ssopinv = nl.inv(ssop[0])
    1637             mst = ssopinv[3][:3]
    1638             epsinv = ssopinv[3][3]
    1639             sdet = nl.det(sop[0])
    1640             ssdet = nl.det(ssop[0])
    1641             dtau = mst*(XYZ-sop[1])-epsinv*ssop[1][3]
    1642             dT = 1.0
    1643             if np.any(dtau%.5):
    1644                 dT = np.tan(np.pi*np.sum(dtau%.5))
    1645             tauT = np.inner(mst,XYZ-sop[1])+epsinv*(tau-ssop[1][3])
     1645            sdet,ssdet,dtau,dT,tauT = getTauT(tau,sop,ssop,XYZ)
    16461646            usc = np.ones(12,dtype='i')
    16471647            dUT = posFourier(tauT,nH,delt12[:6],delt12[6:])                  #Uij modulations - 6x12x49 array
    … …  
    21142114    else:
    21152115        cellA = np.zeros(3)
    2116     newX = Cen+(1-2*iC)*(np.inner(M,X)+T)+cellA
     2116    newX = Cen+(1-2*iC)*(np.inner(M,X).T+T)+cellA
    21172117    if len(Uij):
    21182118        U = Uij2U(Uij)

  • trunk/GSASIIstrMath.py

    r1955 r1957  
    629629                if parm in parmDict:
    630630                    keys[key][m][iatm] = parmDict[parm]
    631     return waveTypes,FSSdata.squeeze(),XSSdata.squeeze(),USSdata.squeeze(),MSSdata.squeeze()   
    632    
    633 def SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
    634     '''
    635     Compute super structure factors for all h,k,l,m for phase
    636     puts the result, F^2, in each ref[8+im] in refList
    637     input:
    638    
    639     :param dict refDict: where
    640         'RefList' list where each ref = h,k,l,m,d,...
    641         'FF' dict of form factors - filed in below
    642     :param np.array G:      reciprocal metric tensor
    643     :param str pfx:    phase id string
    644     :param dict SGData: space group info. dictionary output from SpcGroup
    645     :param dict calcControls:
    646     :param dict ParmDict:
    647 
    648     '''       
    649     phfx = pfx.split(':')[0]+hfx
    650     ast = np.sqrt(np.diag(G))
    651     Mast = twopisq*np.multiply.outer(ast,ast)
    652     SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
    653     SGT = np.array([ops[1] for ops in SGData['SGOps']])
    654     SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
    655     SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
    656     FFtables = calcControls['FFtables']
    657     BLtables = calcControls['BLtables']
    658     Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
    659     waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
    660     FF = np.zeros(len(Tdata))
    661     if 'NC' in calcControls[hfx+'histType']:
    662         FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
    663     else:
    664         FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
    665         FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
    666     Uij = np.array(G2lat.U6toUij(Uijdata))
    667     bij = Mast*Uij.T
    668     if not len(refDict['FF']):
    669         if 'N' in calcControls[hfx+'histType']:
    670             dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
    671         else:
    672             dat = G2el.getFFvalues(FFtables,0.)       
    673         refDict['FF']['El'] = dat.keys()
    674         refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))   
    675     for iref,refl in enumerate(refDict['RefList']):
    676         if 'NT' in calcControls[hfx+'histType']:
    677             FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl[14+im])
    678         fbs = np.array([0,0])
    679         H = refl[:4]
    680         SQ = 1./(2.*refl[4+im])**2
    681         SQfactor = 4.0*SQ*twopisq
    682         Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
    683         if not np.any(refDict['FF']['FF'][iref]):                #no form factors - 1st time thru StructureFactor
    684             if 'N' in calcControls[hfx+'histType']:
    685                 dat = G2el.getBLvalues(BLtables)
    686                 refDict['FF']['FF'][iref] = dat.values()
    687             else:       #'X'
    688                 dat = G2el.getFFvalues(FFtables,SQ)
    689                 refDict['FF']['FF'][iref] = dat.values()
    690         Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
    691         FF = refDict['FF']['FF'][iref][Tindx]
    692         Uniq = np.inner(H[:3],SGMT)
    693         SSUniq = np.inner(H,SSGMT)
    694         Phi = np.inner(H[:3],SGT)
    695         SSPhi = np.inner(H,SSGT)
    696         GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)       
    697         phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+Phi[:,np.newaxis])
    698         sinp = np.sin(phase)
    699         cosp = np.cos(phase)
    700         biso = -SQfactor*Uisodata
    701         Tiso = np.where(biso<1.,np.exp(biso),1.0)
    702         HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
    703         Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
    704         Tcorr = Tiso*Tuij*Mdata*Fdata/len(Uniq)
    705         fa = np.array([(FF+FP-Bab)*cosp*Tcorr,-FPP*sinp*Tcorr])
    706         fb = np.zeros_like(fa)
    707         if not SGData['SGInv']:
    708             fb = np.array([(FF+FP-Bab)*sinp*Tcorr,FPP*cosp*Tcorr])
    709         fa = fa*GfpuA-fb*GfpuB
    710         fb = fb*GfpuA+fa*GfpuB
    711         fas = np.real(np.sum(np.sum(fa,axis=1),axis=1))        #real
    712         fbs = np.real(np.sum(np.sum(fb,axis=1),axis=1))
    713            
    714         fasq = fas**2
    715         fbsq = fbs**2        #imaginary
    716         refl[9+im] = np.sum(fasq)+np.sum(fbsq)
    717         refl[7+im] = np.sum(fasq)+np.sum(fbsq)
    718         refl[10+im] = atan2d(fbs[0],fas[0])
     631    return waveTypes,FSSdata.squeeze(),XSSdata.squeeze(),USSdata.squeeze(),MSSdata.squeeze()
     632   
     633def GetSSTauM(SGOps,SSOps,pfx,calcControls,XData):
     634   
     635    Natoms = calcControls['Natoms'][pfx]
     636    maxSSwave = calcControls['maxSSwave'][pfx]
     637    Smult = np.zeros((Natoms,len(SGOps)))
     638    TauT = np.zeros((Natoms,len(SGOps)))
     639    for ix,xyz in enumerate(XData.T):
     640        for isym,(sop,ssop) in enumerate(zip(SGOps,SSOps)):
     641            sdet,ssdet,dtau,dT,tauT = G2spc.getTauT(0,sop,ssop,xyz)
     642            Smult[ix][isym] = sdet*ssdet
     643            TauT[ix][isym] = tauT
     644    return Smult,TauT
    719645   
    720646def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
    … …  
    1024950    return dFdvDict
    1025951   
     952def SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
     953    '''
     954    Compute super structure factors for all h,k,l,m for phase
     955    puts the result, F^2, in each ref[8+im] in refList
     956    input:
     957   
     958    :param dict refDict: where
     959        'RefList' list where each ref = h,k,l,m,d,...
     960        'FF' dict of form factors - filed in below
     961    :param np.array G:      reciprocal metric tensor
     962    :param str pfx:    phase id string
     963    :param dict SGData: space group info. dictionary output from SpcGroup
     964    :param dict calcControls:
     965    :param dict ParmDict:
     966
     967    '''       
     968    phfx = pfx.split(':')[0]+hfx
     969    ast = np.sqrt(np.diag(G))
     970    Mast = twopisq*np.multiply.outer(ast,ast)
     971   
     972    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     973    SGT = np.array([ops[1] for ops in SGData['SGOps']])
     974    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
     975    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
     976    FFtables = calcControls['FFtables']
     977    BLtables = calcControls['BLtables']
     978    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     979    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
     980    SStauM = GetSSTauM(SGData['SGOps'],SSGData['SSGOps'],pfx,calcControls,Xdata)
     981    FF = np.zeros(len(Tdata))
     982    if 'NC' in calcControls[hfx+'histType']:
     983        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
     984    else:
     985        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
     986        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
     987    Uij = np.array(G2lat.U6toUij(Uijdata))
     988    bij = Mast*Uij.T
     989    if not len(refDict['FF']):
     990        if 'N' in calcControls[hfx+'histType']:
     991            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
     992        else:
     993            dat = G2el.getFFvalues(FFtables,0.)       
     994        refDict['FF']['El'] = dat.keys()
     995        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))   
     996    for iref,refl in enumerate(refDict['RefList']):
     997        if 'NT' in calcControls[hfx+'histType']:
     998            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl[14+im])
     999        fbs = np.array([0,0])
     1000        H = refl[:4]
     1001        SQ = 1./(2.*refl[4+im])**2
     1002        SQfactor = 4.0*SQ*twopisq
     1003        Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
     1004        if not np.any(refDict['FF']['FF'][iref]):                #no form factors - 1st time thru StructureFactor
     1005            if 'N' in calcControls[hfx+'histType']:
     1006                dat = G2el.getBLvalues(BLtables)
     1007                refDict['FF']['FF'][iref] = dat.values()
     1008            else:       #'X'
     1009                dat = G2el.getFFvalues(FFtables,SQ)
     1010                refDict['FF']['FF'][iref] = dat.values()
     1011        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
     1012        FF = refDict['FF']['FF'][iref][Tindx]
     1013        Uniq = np.inner(H[:3],SGMT)
     1014        SSUniq = np.inner(H,SSGMT)
     1015        Phi = np.inner(H[:3],SGT)
     1016        SSPhi = np.inner(H,SSGT)
     1017        GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM)       
     1018        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+Phi[:,np.newaxis])
     1019        sinp = np.sin(phase)
     1020        cosp = np.cos(phase)
     1021        biso = -SQfactor*Uisodata
     1022        Tiso = np.where(biso<1.,np.exp(biso),1.0)
     1023        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
     1024        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
     1025        Tcorr = Tiso*Tuij*Mdata*Fdata/len(Uniq)
     1026        fa = np.array([(FF+FP-Bab)*cosp*Tcorr,-FPP*sinp*Tcorr])     #2 x sym x atoms
     1027        fb = np.zeros_like(fa)                                      #ditto
     1028        if not SGData['SGInv']:
     1029            fb = np.array([(FF+FP-Bab)*sinp*Tcorr,FPP*cosp*Tcorr])
     1030        fa = fa*GfpuA[np.newaxis,:,:]-fb*GfpuB[np.newaxis,:,:]
     1031        fb = fb*GfpuA[np.newaxis,:,:]+fa*GfpuB[np.newaxis,:,:]
     1032        fas = np.real(np.sum(np.sum(fa,axis=1),axis=1))        #real
     1033        fbs = np.real(np.sum(np.sum(fb,axis=1),axis=1))
     1034        fasq = fas**2
     1035        fbsq = fbs**2        #imaginary
     1036        refl[9+im] = np.sum(fasq)+np.sum(fbsq)
     1037        refl[7+im] = np.sum(fasq)+np.sum(fbsq)
     1038        refl[10+im] = atan2d(fbs[0],fas[0])
     1039   
    10261040def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):
    10271041    'Needs a doc string'
    … …  
    10381052    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
    10391053    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
     1054    SStauM = GetSSTauM(SGData['SGOps'],SSGData['SSGOps'],pfx,calcControls,Xdata)
    10401055    mSize = len(Mdata)
    10411056    FF = np.zeros(len(Tdata))
    … …  
    10671082        Phi = np.inner(H[:3],SGT)
    10681083        SSPhi = np.inner(H,SSGT)
    1069         GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)
    1070         dGAdk,dGBdk = G2mth.ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata)
     1084        GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM)
     1085        dGAdk,dGBdk = G2mth.ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM)
    10711086        #need ModulationDerv here dGAdXsin, etc 
    10721087        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+Phi[np.newaxis,:])
    … …  
    10851100        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
    10861101        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
    1087         GfpuA = np.swapaxes(GfpuA,1,2)
    1088         GfpuB = np.swapaxes(GfpuB,1,2)
    1089         fa = fa*GfpuA-fb*GfpuB
    1090         fb = fb*GfpuA+fa*GfpuB
     1102        fa = fa*GfpuA[:,:,np.newaxis]-fb*GfpuB[:,:,np.newaxis]
     1103        fb = fb*GfpuA[:,:,np.newaxis]+fa*GfpuB[:,:,np.newaxis]
    10911104       
    10921105        fas = np.sum(np.sum(fa,axis=1),axis=1)
    … …  
    10941107        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
    10951108        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
    1096         fax = fax*GfpuA-fbx*GfpuB
    1097         fbx = fbx*GfpuA+fax*GfpuB
     1109        fax = fax*GfpuA[:,:,np.newaxis]-fbx*GfpuB[:,:,np.newaxis]
     1110        fbx = fbx*GfpuA[:,:,np.newaxis]+fax*GfpuB[:,:,np.newaxis]
    10981111        #sum below is over Uniq
    10991112        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
Note: See TracChangeset for help on using the changeset viewer.