Changeset 4533


Ignore:
Timestamp:
Jul 23, 2020 12:28:51 PM (3 years ago)
Author:
vondreele
Message:

fix math error in incomm mag str fctr stuff $ expand equations

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIImath.py

    r4526 r4533  
    14371437    fxn of gTau points; NB: this allows only 1 mag. wave fxn.
    14381438    '''
    1439     Am = np.array(MSSdata[3:]).T   #atoms x waves x cos pos mods
    1440     Bm = np.array(MSSdata[:3]).T  #...sin pos mods
    1441     nWaves = Am.shape[1]
     1439    Am = np.array(MSSdata[3:]).T[:,0,:]   #atoms x cos pos mods; only 1 wave used
     1440    Bm = np.array(MSSdata[:3]).T[:,0,:]   #...sin pos mods
    14421441    SGMT = np.array([ops[0] for ops in SGData['SGOps']])        #not .T!!
    14431442    Sinv = np.array([nl.inv(ops[0]) for ops in SSGData['SSGOps']])
    … …  
    14531452        SGMT = np.vstack((SGMT,SGMT))
    14541453        Sinv = np.vstack((Sinv,Sinv))
    1455         SGT = np.vstack((SGT,SGT+.5))%1.
     1454        SGT = np.vstack((SGT,SGT+np.array([0.,0.,0.,.5])))%1.
    14561455    mst = Sinv[:,3,:3]
    14571456    epsinv = Sinv[:,3,3]
    14581457    phi = np.inner(XYZ,modQ).T
    14591458    TA = np.sum(mst[nxs,:,:]*(XYZ-SGT[:,:3][nxs,:,:]),axis=-1).T
    1460     tauT =  TA[nxs,:,:] + epsinv[nxs,:,nxs]*(glTau[:,nxs,nxs]-SGT[:,3][nxs,:,nxs]+phi[nxs,:,:])
    1461     modind = np.arange(nWaves)+1.
    1462     phase = modind[:,nxs,nxs]*tauT     #Nops,Natm,Nwave
    1463     psin = np.sin(twopi*phase)
    1464     pcos = np.cos(twopi*phase)
    1465     MmodA = np.sum(Am[nxs,nxs,:,:,:]*pcos[:,:,:,nxs,nxs],axis=3)/2.    #cos term
    1466     MmodB = np.sum(Bm[nxs,nxs,:,:,:]*psin[:,:,:,nxs,nxs],axis=3)/2.    #sin term
    1467     MmodA = np.sum(SGMT[nxs,:,nxs,:,:]*MmodA[:,:,:,nxs,:],axis=-1)*SGData['SpnFlp'][nxs,:,nxs,nxs]
    1468     MmodB = np.sum(SGMT[nxs,:,nxs,:,:]*MmodB[:,:,:,nxs,:],axis=-1)*SGData['SpnFlp'][nxs,:,nxs,nxs]
    1469     return MmodA,MmodB    #Ntau,Nops,Natm,Mxyz; cos & sin parts; sum matches drawn atom moments
    1470        
    1471 def MagMod2(m,XYZ,modQ,MSSdata,SGData,SSGData):
    1472     '''
    1473     this needs to make magnetic moment modulations & magnitudes as
    1474     fxn of gTau points; NB: this allows only 1 mag. wave fxn.
    1475     '''
    1476     Am = np.array(MSSdata[3:]).T[:,0,:]   #atoms x cos pos mods; only 1 wave
    1477     Bm = np.array(MSSdata[:3]).T[:,0,:]  #...sin pos mods
    1478     SGMT = np.array([ops[0] for ops in SGData['SGOps']])        #not .T!!
    1479     SSGMT = np.array([ops[0] for ops in SSGData['SSGOps']])        #not .T!!
    1480     Sinv = np.array([nl.inv(ops[0]) for ops in SSGData['SSGOps']])
    1481     SGT = np.array([ops[1] for ops in SSGData['SSGOps']])
    1482     if SGData['SGInv']:
    1483         SGMT = np.vstack((SGMT,-SGMT))
    1484         SSGMT = np.vstack((SSGMT,-SSGMT))
    1485         Sinv = np.vstack((Sinv,-Sinv))
    1486         SGT = np.vstack((SGT,-SGT))
    1487     SGMT = np.vstack([SGMT for cen in SGData['SGCen']])
    1488     SSGMT = np.vstack([SSGMT for cen in SGData['SGCen']])
    1489     Sinv = np.vstack([Sinv for cen in SGData['SGCen']])
    1490     SGT = np.vstack([SGT+cen for cen in SSGData['SSGCen']])%1.
    1491     if SGData['SGGray']:
    1492         SGMT = np.vstack((SGMT,SGMT))
    1493         SSGMT = np.vstack((SSGMT,SSGMT))
    1494         Sinv = np.vstack((Sinv,Sinv))
    1495         SGT = np.vstack((SGT,SGT+.5))%1.
    1496     epsinv = Sinv[:,3,3]
    1497     phi = np.inner(XYZ,modQ).T
    1498     TA = phi+(epsinv*(np.inner(modQ,SGT[:,:3])-SGT[:,3]))[:,nxs]    #Nops,Natm
    1499     phase = phi+(np.inner(modQ,SGT[:,:3])-SGT[:,3])[:,nxs]
    1500    
    1501     pcos = np.cos(-twopi*m[:,nxs,nxs]*phase[nxs,:,:])      #Nref,Nops,Natm
    1502     psin = np.sin(-twopi*m[:,nxs,nxs]*phase[nxs,:,:])
    1503     MmodA = TA[nxs,:,:,nxs]*(Am[nxs,nxs,:,:]*pcos[:,:,:,nxs]-Bm[nxs,nxs,:,:]*psin[:,:,:,nxs])/2.    #Nref,Nops,Natm,Mxyz
    1504     MmodB = TA[nxs,:,:,nxs]*(Am[nxs,nxs,:,:]*psin[:,:,:,nxs]+Bm[nxs,nxs,:,:]*pcos[:,:,:,nxs])/2.    #Nref,Nops,Natm,Mxyz
    1505     MmodA = np.sum(SGMT[nxs,:,nxs,:,:]*MmodA[:,:,:,nxs,:],axis=-1)*SGData['SpnFlp'][nxs,:,nxs,nxs]
    1506     MmodB = np.sum(SGMT[nxs,:,nxs,:,:]*MmodB[:,:,:,nxs,:],axis=-1)*SGData['SpnFlp'][nxs,:,nxs,nxs]
    1507     return MmodA,MmodB    #Nref,Nops,Natm,Mxyz; cos & sin parts
    1508        
     1459    phase =  TA[nxs,:,:] + epsinv[nxs,:,nxs]*(glTau[:,nxs,nxs]-SGT[:,3][nxs,:,nxs]+phi[nxs,:,:])
     1460    pcos = np.cos(twopi*phase*epsinv[nxs,:,nxs])      #Ntau,Nops,Natm
     1461    psin = np.sin(twopi*phase*epsinv[nxs,:,nxs])
     1462   
     1463    MmodAp = (Am[nxs,nxs,:,:]*pcos[:,:,:,nxs]+epsinv[nxs,:,nxs,nxs]*Bm[nxs,nxs,:,:]*psin[:,:,:,nxs])/2.    #Ntau,Nops,Natm,Mxyz
     1464    MmodBp = (Am[nxs,nxs,:,:]*psin[:,:,:,nxs]+epsinv[nxs,:,nxs,nxs]*Bm[nxs,nxs,:,:]*pcos[:,:,:,nxs])/2.
     1465    MmodAm = (Am[nxs,nxs,:,:]*pcos[:,:,:,nxs]+epsinv[nxs,:,nxs,nxs]*Bm[nxs,nxs,:,:]*psin[:,:,:,nxs])/2.
     1466    MmodBm = (Am[nxs,nxs,:,:]*psin[:,:,:,nxs]+epsinv[nxs,:,nxs,nxs]*Bm[nxs,nxs,:,:]*pcos[:,:,:,nxs])/2.
     1467    MmodAp = np.sum(SGMT[nxs,:,nxs,:,:]*MmodAp[:,:,:,nxs,:],axis=-1)*SGData['MagMom'][nxs,:,nxs,nxs]
     1468    MmodBp = np.sum(SGMT[nxs,:,nxs,:,:]*MmodBp[:,:,:,nxs,:],axis=-1)*SGData['MagMom'][nxs,:,nxs,nxs]
     1469    MmodAm = np.sum(SGMT[nxs,:,nxs,:,:]*MmodAm[:,:,:,nxs,:],axis=-1)*SGData['MagMom'][nxs,:,nxs,nxs]
     1470    MmodBm = np.sum(SGMT[nxs,:,nxs,:,:]*MmodBm[:,:,:,nxs,:],axis=-1)*SGData['MagMom'][nxs,:,nxs,nxs]
     1471    return MmodAp,MmodBp,MmodAm,MmodBm    #Ntau,Nops,Natm,Mxyz; cos & sin parts
     1472               
    15091473def Modulation(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
    15101474    '''

  • trunk/GSASIIstrMath.py

    r4531 r4533  
    15161516        if SGData['SGGray']:
    15171517            mXYZ = np.hstack((mXYZ,mXYZ))
    1518 # this done in wrong place - should be
    1519         MmodA,MmodB = G2mth.MagMod(glTau,mXYZ,modQ,MSSdata,SGData,SSGData)  #Ntau,Nops,Natm,Mxyz cos,sim parts sum matches drawing
    1520         Mmod = MmodA+MmodB
     1518
     1519        MmodAp,MmodBp,MmodAm,MmodBm = G2mth.MagMod(glTau,mXYZ,modQ,MSSdata,SGData,SSGData)  #Ntau,Nops,Natm,Mxyz cos,sim parts sum matches drawing
     1520#        MmodAp,MmodBp,MmodAm,MmodBm = G2mth.MagMod2(mXYZ,modQ,MSSdata,SGData,SSGData)  #Nops,Natm,Mxyz cos,sim parts sum matches drawing
    15211521       
    15221522        if not SGData['SGGray']:    #for fixed Mx,My,Mz
    … …  
    15271527            GSdata = SGData['MagMom'][nxs,:,nxs]*GSdata   #flip vectors according to spin flip * det(opM)
    15281528            GSdata = np.swapaxes(GSdata,0,1)    #Nop,Natm,Mxyz
    1529             Mmod += GSdata[nxs,:,:,:]
    1530 
    1531 #        Kdata = np.inner(Mmod,uAmat)          #Ntau,Nops,Natm,Mxyz
    1532 #        Kmean = np.mean(np.sqrt(np.sum(Kdata**2,axis=0)),axis=0)
    1533 #        Kdata /= Kmean     #Ntau,Nops,Natm,Mxyz  Cartesian unit vectors along mag moments
    15341529       
    15351530    FF = np.zeros(len(Tdata))
    … …  
    16081603            eM = (HM/np.sqrt(np.sum(HM**2,axis=0))).T    # normalize  HP  Nref,hkl=Unit vectors || Q
    16091604
    1610             fam0 = 0.
    1611             fbm0 = 0.
     1605#            fam0 = 0.
     1606#            fbm0 = 0.
    16121607            if not SGData['SGGray']:     #correct -fixed Mx,My,Mz contribution             
    16131608                fam0 = TMcorr[:,nxs,:,nxs]*GSdata[nxs,:,:,:]*cosm[:,:,:,nxs]    #Nref,Nops,Natm,Mxyz
    16141609                fbm0 = TMcorr[:,nxs,:,nxs]*GSdata[nxs,:,:,:]*sinm[:,:,:,nxs]
    16151610#  calc mag. structure factors; Nref,Ntau,Nops,Natm,Mxyz                           
    1616             fams = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(H[3,i]!=0,(MmodA*cosm[i,nxs,:,:,nxs]+   
    1617                 H[3,i]*MmodB*sinm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])          #Nref,Ntau,Nops,Natm,Mxyz
     1611#            fams = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(H[3,i]!=0,(MmodA*cosm[i,nxs,:,:,nxs]+   
     1612#                H[3,i]*MmodB*sinm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])          #Nref,Ntau,Nops,Natm,Mxyz
     1613#                       
     1614#            fbms = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(H[3,i]!=0,(MmodA*sinm[i,nxs,:,:,nxs]+   
     1615#                H[3,i]*MmodB*cosm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])          #Nref,Ntau,Nops,Natm,Mxyz
     1616#
     1617#            if not SGData['SGGray']:           
     1618#                fams += fam0[:,nxs,:,:,:]
     1619#                fbms += fbm0[:,nxs,:,:,:]
     1620               
     1621            fams = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(H[3,i]!=0,np.where(H[3,i]>0,
     1622                (MmodAp*cosm[i,nxs,:,:,nxs]+MmodBp*sinm[i,nxs,:,:,nxs]),
     1623                (MmodAm*cosm[i,nxs,:,:,nxs]+MmodBm*sinm[i,nxs,:,:,nxs])),0.) for i in range(mRef)])  #Nref,Nops,Natm,Mxyz
    16181624                       
    1619             fbms = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(H[3,i]!=0,(MmodA*sinm[i,nxs,:,:,nxs]+   
    1620                 H[3,i]*MmodB*cosm[i,nxs,:,:,nxs]),0.) for i in range(mRef)])          #Nref,Ntau,Nops,Natm,Mxyz
     1625            fbms = TMcorr[:,nxs,nxs,:,nxs]*np.array([np.where(H[3,i]!=0,np.where(H[3,i]>0,
     1626                (MmodAp*sinm[i,nxs,:,:,nxs]+MmodBp*cosm[i,nxs,:,:,nxs]),
     1627                (MmodAm*sinm[i,nxs,:,:,nxs]+MmodBm*cosm[i,nxs,:,:,nxs])),0.) for i in range(mRef)])  #Nref,Nops,Natm,Mxyz
    16211628
    16221629            if not SGData['SGGray']:           
    … …  
    16241631                fbms += fbm0[:,nxs,:,:,:]
    16251632               
    1626 # this is the best, but not exactly right
    16271633#sum ops & atms                               
    1628             fasm = np.sum(np.sum(fams,axis=-2),axis=-2)    #Nref,Ntau,Mxyz; sum ops & atoms
     1634            fasm = np.sum(np.sum(fams,axis=-2),axis=-2)    #Nref,Mxyz; sum ops & atoms
    16291635            fbsm = np.sum(np.sum(fbms,axis=-2),axis=-2)
    16301636#put into cartesian space
    … …  
    16351641            eDotFb = np.sum(eM[:,nxs,:]*fbcm,axis=-1)
    16361642#intensity Halpern & Johnson
    1637 #            fass = np.sum((facm-eM[:,nxs,:]*eDotFa[:,:,nxs])**2,axis=-1)
    1638 #            fbss = np.sum((fbcm-eM[:,nxs,:]*eDotFb[:,:,nxs])**2,axis=-1)
    1639             fass = np.sum(facm**2,axis=-1)-eDotFa**2
    1640             fbss = np.sum(fbcm**2,axis=-1)-eDotFb**2
     1643            fass = np.sum((facm-eM[:,nxs,:]*eDotFa[:,:,nxs])**2,axis=-1)
     1644            fbss = np.sum((fbcm-eM[:,nxs,:]*eDotFb[:,:,nxs])**2,axis=-1)
    16411645## #do integration           
    16421646            fas = np.sum(fass*glWt[nxs,:],axis=1)
Note: See TracChangeset for help on using the changeset viewer.