Changeset 4139


Ignore:
Timestamp:
Sep 6, 2019 3:35:55 PM (2 years ago)
Author:
vondreele
Message:

more mag. str. fctr. stuff

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r4138 r4139  
    14681468    return np.array([cosHA,sinHA])             # 2 x refBlk x SGops x atoms
    14691469   
    1470 #def MagModulation(H,HP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt):
    1471 #    '''
    1472 #    H: array nRefBlk x ops X hklt
    1473 #    HP: array nRefBlk x ops X hklt proj to hkl
    1474 #    nWaves: list number of waves for frac, pos, uij & mag
    1475 #    Fmod: array 2 x atoms x waves    (sin,cos terms)
    1476 #    Xmod: array atoms X 3 X ngl
    1477 #    Umod: array atoms x 3x3 x ngl
    1478 #    glTau,glWt: arrays Gauss-Lorentzian pos & wts
    1479 #    '''
    1480 #   
    1481 #    if nWaves[2]:       #uij (adp) waves
    1482 #        if len(HP.shape) > 2:
    1483 #            HbH = np.exp(-np.sum(HP[:,:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
    1484 #        else:
    1485 #            HbH = np.exp(-np.sum(HP[:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
    1486 #    else:
    1487 #        HbH = 1.0
    1488 #    HdotX = np.inner(HP,Xmod)                   #refBlk x ops x atoms X ngl
    1489 #    if len(H.shape) > 2:
    1490 #        D = H[:,:,3:]*glTau[nxs,nxs,:]              #m*e*tau; refBlk x ops X ngl
    1491 #        HdotXD = twopi*(HdotX+D[:,:,nxs,:])
    1492 #    else:
    1493 #        D = H[:,3:]*glTau[nxs,:]              #m*e*tau; refBlk x ops X ngl
    1494 #        HdotXD = twopi*(HdotX+D[:,nxs,:])
    1495 #    M = np.swapaxes(Mmod,1,2)
    1496 #    cosHA = np.sum(M[nxs,nxs,:,:,:]*(Fmod*HbH*np.cos(HdotXD)[:,:,:,nxs,:]*glWt),axis=-1)       #real part; refBlk X ops x atoms; sum for G-L integration
    1497 #    sinHA = np.sum(M[nxs,nxs,:,:,:]*(Fmod*HbH*np.sin(HdotXD)[:,:,:,nxs,:]*glWt),axis=-1)       #imag part; ditto
    1498 #    return np.array([cosHA,sinHA])             # 2 x refBlk x SGops x atoms
    1499 #
     1470def MagModulation(H,HP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt):
     1471    '''
     1472    H: array nRefBlk x ops X hklt
     1473    HP: array nRefBlk x ops X hklt proj to hkl
     1474    nWaves: list number of waves for frac, pos, uij & mag
     1475    Fmod: array 2 x atoms x waves    (sin,cos terms)
     1476    Xmod: array atoms X 3 X ngl
     1477    Umod: array atoms x 3x3 x ngl
     1478    Mmod: array atoms x 3x3 x ngl
     1479    glTau,glWt: arrays Gauss-Lorentzian pos & wts
     1480    '''
     1481   
     1482    if nWaves[2]:       #uij (adp) waves
     1483        if len(HP.shape) > 2:
     1484            HbH = np.exp(-np.sum(HP[:,:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
     1485        else:
     1486            HbH = np.exp(-np.sum(HP[:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
     1487    else:
     1488        HbH = 1.0
     1489    HdotX = np.inner(HP,Xmod)                   #refBlk x ops x atoms X ngl
     1490    if len(H.shape) > 2:
     1491        D = H[:,:,3:]*glTau[nxs,nxs,:]              #m*e*tau; refBlk x ops X ngl
     1492        HdotXD = twopi*(HdotX+D[:,:,nxs,:])
     1493    else:
     1494        D = H[:,3:]*glTau[nxs,:]              #m*e*tau; refBlk x ops X ngl
     1495        HdotXD = twopi*(HdotX+D[:,nxs,:])
     1496    M = np.swapaxes(Mmod,1,2)
     1497    cosHA = np.sum(M[nxs,nxs,:,:,:]*(Fmod*HbH*np.cos(HdotXD)[:,:,:,nxs,:]*glWt),axis=-1)       #real part; refBlk X ops x atoms; sum for G-L integration
     1498    sinHA = np.sum(M[nxs,nxs,:,:,:]*(Fmod*HbH*np.sin(HdotXD)[:,:,:,nxs,:]*glWt),axis=-1)       #imag part; ditto
     1499    return np.array([cosHA,sinHA])             # 2 x refBlk x SGops x atoms
     1500
    15001501def ModulationTw(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
    15011502    '''
  • trunk/GSASIIstrMath.py

    r4138 r4139  
    14821482    phfx = pfx.split(':')[0]+hfx
    14831483    ast = np.sqrt(np.diag(G))
    1484 #    GS = G/np.outer(ast,ast)
    1485 #    uAmat = G2lat.Gmat2AB(GS)[0]
     1484    g = nl.inv(G)
     1485    ainv = np.sqrt(np.diag(g))
     1486    GS = G/np.outer(ast,ast)
     1487    uAmat = G2lat.Gmat2AB(GS)[0]
     1488    Ginv = g/np.outer(ainv,ainv)
    14861489    Mast = twopisq*np.multiply.outer(ast,ast)   
    14871490    SGInv = SGData['SGInv']
    14881491    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
    1489     Ncen = len(SGData['SGCen'])
     1492#    Ncen = len(SGData['SGCen'])
    14901493    Nops = len(SGMT)*(1+SGData['SGInv'])
    14911494    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
     
    15111514        mXYZ = np.array([[xyz[0] for xyz in list(G2spc.GenAtom(xyz,SGData,All=True,Move=True))] for xyz in (Xdata+dXdata).T])%1. #Natn,Nop,xyz
    15121515        MmodA,MmodB = G2mth.MagMod(glTau,mXYZ,modQ,MSSdata,SGData,SSGData)  #Ntau,Nops,Natm,Mxyz sum matches drawing
     1516        Mmod = MmodA+MmodB
    15131517       
    15141518        if not SGData['SGGray']:    #for fixed Mx,My,Mz
     
    15191523            GSdata = SGData['MagMom'][nxs,:,nxs]*GSdata   #flip vectors according to spin flip * det(opM)
    15201524            GSdata = np.swapaxes(GSdata,0,1)    #Nop,Natm,Mxyz
     1525            Mmod += GSdata[nxs,:,:,:]
     1526           
     1527        Mag = np.array([np.sqrt(np.inner(mag,np.inner(mag,Ginv))) for mag in Mmod])     #Ntau,Nop,Natm
     1528        Kdata = np.inner(Mmod,uAmat)     #Cartesian vectors; Mxyz,Natm,Nop,Ntau
     1529        Kmean = np.sqrt(np.sum(Kdata**2,axis=-1))
     1530        Kdata /= Kmean[:,:,:,nxs]      #Ntau,Nop,Natm,Mxyz
     1531
    15211532
    15221533    FF = np.zeros(len(Tdata))
     
    15951606#                     
    15961607            HM = np.inner(Bmat,HP.T)                            #put into cartesian space
    1597             eM = HM/np.sqrt(np.sum(HM**2,axis=0))               #& normalize
    1598 #for fixed moments --> m=0 reflections
     1608            eM = HM/np.sqrt(np.sum(HM**2,axis=0)).T               #& normalize; Nref,Mhkl
     1609            eDotK = np.sum(eM[:,nxs,nxs,:]*Kdata[:,nxs,:,:],axis=0)
     1610           
    15991611            fam0 = 0.
    16001612            fbm0 = 0.
     
    16111623           
    16121624            if not SGData['SGGray']:
    1613                 fams *= 0.5
    1614                 fbms *= 0.5
     1625                fams *= .5
     1626                fbms *= .5
    16151627                fams += fam0[:,nxs,:,:,:]
    16161628                fbms += fbm0[:,nxs,:,:,:]
    1617                        
     1629               
     1630#   do GL integral before making mag F's - less correct; gives F's where there shouldn't be any                     
     1631           
     1632            faqs = np.sum(fams**2,axis=-1)*(1.-eDotK**2)      #mag intensity calc F^2-(e.F)^2
     1633            fbqs = np.sum(fbms**2,axis=-1)*(1.-eDotK**2)
     1634           
     1635            fass = np.sum(glWt[nxs,:,nxs,nxs]*faqs,axis=1)         #GL integration
     1636            fbss = np.sum(glWt[nxs,:,nxs,nxs]*fbqs,axis=1)
     1637           
     1638            fas = np.sum(np.sum(fass,axis=-1),axis=-1)      #Nref,Ntau,Mxyz; sum ops & atoms
     1639            fbs = np.sum(np.sum(fbss,axis=-1),axis=-1)
     1640
     1641# do GL integral after making mag. F's - better but not right               
    16181642            famqs = np.sum(np.sum(fams,axis=-2),axis=-2)      #Nref,Ntau,Mxyz; sum ops & atoms
    16191643            fbmqs = np.sum(np.sum(fbms,axis=-2),axis=-2)
     
    16311655            fbss = np.sum(fbmqs**2,axis=-1)*(1.-np.sum(eM.T[:,nxs,:]*fbmcs,axis=-1)**2)
    16321656           
    1633             fas = np.sum(glWt*fass,axis=1)
     1657            fas = np.sum(glWt*fass,axis=1)          #GL integration
    16341658            fbs = np.sum(glWt*fbss,axis=1)
    16351659           
Note: See TracChangeset for help on using the changeset viewer.