Changeset 3372


Ignore:
Timestamp:
May 6, 2018 1:42:01 PM (5 years ago)
Author:
vondreele
Message:

disable printing of debug info for binary location & config controls
fix bug in image plotting wrt linescan
better (but not best) mag moment derivatives

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIpath.py

    r3303 r3372  
    740740
    741741BinaryPathLoaded = False
    742 def SetBinaryPath(printInfo=True):
     742def SetBinaryPath(printInfo=False):
    743743    '''
    744744    Add location of GSAS-II shared libraries (binaries: .so or .pyd files) to path

  • trunk/GSASIIplot.py

    r3366 r3372  
    60396039        except TypeError:
    60406040            return
    6041         if Data['linescan']:
     6041        if Data['linescan'][0]:
    60426042            Page.xlim1 = Plot1.get_xlim()
    60436043        new,plotNum,Page,Plot,lim = G2frame.G2plotNB.FindPlotTab('2D Powder Image','mpl',newImage=False)

  • trunk/GSASIIstrMath.py

    r3353 r3372  
    770770#    print 'sf time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
    771771   
    772 def MagStructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
    773     ''' Compute neutron magnetic structure factors for all h,k,l for phase
    774     puts the result, F^2, in each ref[8] in refList
    775     operates on blocks of 100 reflections for speed
    776     input:
    777    
    778     :param dict refDict: where
    779         'RefList' list where each ref = h,k,l,it,d,...
    780         'FF' dict of form factors - filed in below
    781     :param np.array G:      reciprocal metric tensor
    782     :param str pfx:    phase id string
    783     :param dict SGData: space group info. dictionary output from SpcGroup
    784     :param dict calcControls:
    785     :param dict ParmDict:
    786 
    787     '''       
    788     ast = np.sqrt(np.diag(G))
    789     GS = G/np.outer(ast,ast)
    790     uAmat = G2lat.Gmat2AB(GS)[0]
    791     Mast = twopisq*np.multiply.outer(ast,ast)
    792     SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
    793     SGT = np.array([ops[1] for ops in SGData['SGOps']])
    794     Ncen = len(SGData['SGCen'])
    795     Nops = len(SGMT)*Ncen
    796     if not SGData['SGFixed']:
    797         Nops *= (1+SGData['SGInv'])
    798     MFtables = calcControls['MFtables']
    799     Bmat = G2lat.Gmat2AB(G)[1]
    800     TwinLaw = np.ones(1)
    801 #    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
    802 #    TwDict = refDict.get('TwDict',{})           
    803 #    if 'S' in calcControls[hfx+'histType']:
    804 #        NTL = calcControls[phfx+'NTL']
    805 #        NM = calcControls[phfx+'TwinNMN']+1
    806 #        TwinLaw = calcControls[phfx+'TwinLaw']
    807 #        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
    808 #        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
    809     Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
    810         GetAtomFXU(pfx,calcControls,parmDict)
    811     if not Xdata.size:          #no atoms in phase!
    812         return
    813     Mag = np.sqrt(np.array([np.inner(mag,np.inner(mag,GS)) for mag in Gdata.T]))
    814     Gdata = np.inner(Gdata.T,SGMT).T            #apply sym. ops.
    815     if SGData['SGInv'] and not SGData['SGFixed']:
    816         Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
    817     Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering--> [Mxyz,nops,natms]
    818     Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip * det(opM)
    819     Mag = np.tile(Mag[:,nxs],Nops).T  #make Mag same length as Gdata
    820     Gdata = Gdata/Mag       #normalze mag. moments
    821     Gdata = np.inner(Gdata.T,uAmat).T*np.sqrt(nl.det(GS))
    822     Uij = np.array(G2lat.U6toUij(Uijdata))
    823     bij = Mast*Uij.T
    824     blkSize = 100       #no. of reflections in a block - size seems optimal
    825     nRef = refDict['RefList'].shape[0]
    826     SQ = 1./(2.*refDict['RefList'].T[4])**2
    827     refDict['FF']['El'] = list(MFtables.keys())
    828     refDict['FF']['MF'] = np.zeros((nRef,len(MFtables)))
    829     for iel,El in enumerate(refDict['FF']['El']):
    830         refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
    831 #reflection processing begins here - big arrays!
    832     iBeg = 0
    833     while iBeg < nRef:
    834         iFin = min(iBeg+blkSize,nRef)
    835         refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
    836         H = refl.T[:3]                          #array(blkSize,3)
    837 #        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,nTwins,3) or (blkSize,3)
    838 #        TwMask = np.any(H,axis=-1)
    839 #        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
    840 #            for ir in range(blkSize):
    841 #                iref = ir+iBeg
    842 #                if iref in TwDict:
    843 #                    for i in TwDict[iref]:
    844 #                        for n in range(NTL):
    845 #                            H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
    846 #            TwMask = np.any(H,axis=-1)
    847         SQ = 1./(2.*refl.T[4])**2               #array(blkSize)
    848         SQfactor = 4.0*SQ*twopisq               #ditto prev.
    849         Uniq = np.inner(H.T,SGMT)
    850         Phi = np.inner(H.T,SGT)
    851         phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
    852         biso = -SQfactor*Uisodata[:,nxs]
    853         Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
    854         HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
    855         Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
    856         Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
    857         MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
    858         TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm
    859         if SGData['SGInv']:
    860             if not SGData['SGFixed']:
    861                 mphase = np.hstack((phase,-phase))  #OK
    862             else:
    863                 mphase = phase
    864         else:
    865             mphase = phase                    #
    866         mphase = np.array([mphase+twopi*np.inner(cen,H.T)[:,nxs,nxs] for cen in SGData['SGCen']])
    867         mphase = np.concatenate(mphase,axis=1)              #Nref,full Nop,Natm
    868         sinm = np.sin(mphase)                               #ditto - match magstrfc.for
    869         cosm = np.cos(mphase)                               #ditto
    870         HM = np.inner(Bmat.T,H.T)                             #put into cartesian space
    871         HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Gdata = MAGS & HM = UVEC in magstrfc.for both OK
    872         eDotK = np.sum(HM[:,:,nxs,nxs]*Gdata[:,nxs,:,:],axis=0)
    873         Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Gdata[:,nxs,:,:] #xyz,Nref,Nop,Natm = BPM in magstrfc.for OK
    874         fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
    875         fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
    876         fams = np.sum(np.sum(fam,axis=-1),axis=-1)                          #xyz,Nref
    877         fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)                          #ditto
    878         refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)
    879         refl.T[7] = np.copy(refl.T[9])               
    880         refl.T[10] = 0.0    #atan2d(fbs[0],fas[0]) - what is phase for mag refl?
    881 #        if 'P' in calcControls[hfx+'histType']:     #PXC, PNC & PNT: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
    882 #            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) #add fam**2 & fbm**2 here   
    883 #            refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
    884 #        else:                                       #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2
    885 #            if len(TwinLaw) > 1:
    886 #                refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
    887 #                refl.T[7] = np.sum(TwinFr*TwMask*np.sum(fas,axis=0)**2,axis=-1)+   \
    888 #                    np.sum(TwinFr*TwMask*np.sum(fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
    889 #                refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f" & use primary twin
    890 #            else:   # checked correct!!
    891 #                refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
    892 #                refl.T[7] = np.copy(refl.T[9])               
    893 #                refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
    894 ##                refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
    895         iBeg += blkSize
    896 #    print 'sf time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
    897 
    898772def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
    899773    '''Compute structure factor derivatives on blocks of reflections - for powders/nontwins only
    … …  
    1055929    return dFdvDict
    1056930   
    1057 def StructureFactorDervMag(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
     931def MagStructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
     932    ''' Compute neutron magnetic structure factors for all h,k,l for phase
     933    puts the result, F^2, in each ref[8] in refList
     934    operates on blocks of 100 reflections for speed
     935    input:
     936   
     937    :param dict refDict: where
     938        'RefList' list where each ref = h,k,l,it,d,...
     939        'FF' dict of form factors - filed in below
     940    :param np.array G:      reciprocal metric tensor
     941    :param str pfx:    phase id string
     942    :param dict SGData: space group info. dictionary output from SpcGroup
     943    :param dict calcControls:
     944    :param dict ParmDict:
     945
     946    '''       
     947    g = nl.inv(G)
     948    ast = np.sqrt(np.diag(G))
     949    ainv = np.sqrt(np.diag(g))
     950    GS = G/np.outer(ast,ast)
     951    Ginv = g/np.outer(ainv,ainv)
     952    uAmat = G2lat.Gmat2AB(GS)[0]
     953    Mast = twopisq*np.multiply.outer(ast,ast)
     954    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     955    SGT = np.array([ops[1] for ops in SGData['SGOps']])
     956    Ncen = len(SGData['SGCen'])
     957    Nops = len(SGMT)*Ncen
     958    if not SGData['SGFixed']:
     959        Nops *= (1+SGData['SGInv'])
     960    MFtables = calcControls['MFtables']
     961    Bmat = G2lat.Gmat2AB(G)[1]
     962    TwinLaw = np.ones(1)
     963#    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
     964#    TwDict = refDict.get('TwDict',{})           
     965#    if 'S' in calcControls[hfx+'histType']:
     966#        NTL = calcControls[phfx+'NTL']
     967#        NM = calcControls[phfx+'TwinNMN']+1
     968#        TwinLaw = calcControls[phfx+'TwinLaw']
     969#        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
     970#        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
     971    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
     972        GetAtomFXU(pfx,calcControls,parmDict)
     973    if not Xdata.size:          #no atoms in phase!
     974        return
     975    Mag = np.array([np.sqrt(np.inner(mag,np.inner(mag,Ginv))) for mag in Gdata.T])
     976    Gdata = np.inner(Gdata.T,SGMT).T            #apply sym. ops.
     977    if SGData['SGInv'] and not SGData['SGFixed']:
     978        Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
     979    Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering--> [Mxyz,nops,natms]
     980    Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip * det(opM)
     981    Mag = np.tile(Mag[:,nxs],Nops).T  #make Mag same length as Gdata
     982    Kdata = np.inner(Gdata.T,uAmat).T*np.sqrt(nl.det(Ginv))/Mag     #Cartesian unit vectors
     983    Uij = np.array(G2lat.U6toUij(Uijdata))
     984    bij = Mast*Uij.T
     985    blkSize = 100       #no. of reflections in a block - size seems optimal
     986    nRef = refDict['RefList'].shape[0]
     987    SQ = 1./(2.*refDict['RefList'].T[4])**2
     988    refDict['FF']['El'] = list(MFtables.keys())
     989    refDict['FF']['MF'] = np.zeros((nRef,len(MFtables)))
     990    for iel,El in enumerate(refDict['FF']['El']):
     991        refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
     992#reflection processing begins here - big arrays!
     993    iBeg = 0
     994    while iBeg < nRef:
     995        iFin = min(iBeg+blkSize,nRef)
     996        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
     997        H = refl.T[:3].T                          #array(blkSize,3)
     998#        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,nTwins,3) or (blkSize,3)
     999#        TwMask = np.any(H,axis=-1)
     1000#        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
     1001#            for ir in range(blkSize):
     1002#                iref = ir+iBeg
     1003#                if iref in TwDict:
     1004#                    for i in TwDict[iref]:
     1005#                        for n in range(NTL):
     1006#                            H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
     1007#            TwMask = np.any(H,axis=-1)
     1008        SQ = 1./(2.*refl.T[4])**2               #array(blkSize)
     1009        SQfactor = 4.0*SQ*twopisq               #ditto prev.
     1010        Uniq = np.inner(H,SGMT)
     1011        Phi = np.inner(H,SGT)
     1012        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
     1013        biso = -SQfactor*Uisodata[:,nxs]
     1014        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
     1015        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
     1016        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
     1017        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
     1018        MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
     1019        TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm
     1020        if SGData['SGInv']:
     1021            if not SGData['SGFixed']:
     1022                mphase = np.hstack((phase,-phase))  #OK
     1023            else:
     1024                mphase = phase
     1025        else:
     1026            mphase = phase                    #
     1027        mphase = np.array([mphase+twopi*np.inner(cen,H)[:,nxs,nxs] for cen in SGData['SGCen']])
     1028        mphase = np.concatenate(mphase,axis=1)              #Nref,full Nop,Natm
     1029        sinm = np.sin(mphase)                               #ditto - match magstrfc.for
     1030        cosm = np.cos(mphase)                               #ditto
     1031        HM = np.inner(Bmat.T,H)                             #put into cartesian space
     1032        HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Kdata = MAGS & HM = UVEC in magstrfc.for both OK
     1033        eDotK = np.sum(HM[:,:,nxs,nxs]*Kdata[:,nxs,:,:],axis=0)
     1034        Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Kdata[:,nxs,:,:] #xyz,Nref,Nop,Natm = BPM in magstrfc.for OK
     1035        fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
     1036        fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
     1037        fams = np.sum(np.sum(fam,axis=-1),axis=-1)                          #xyz,Nref
     1038        fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)                          #ditto
     1039        refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)
     1040        refl.T[7] = np.copy(refl.T[9])               
     1041        refl.T[10] = 0.0    #atan2d(fbs[0],fas[0]) - what is phase for mag refl?
     1042#        if 'P' in calcControls[hfx+'histType']:     #PXC, PNC & PNT: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
     1043#            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) #add fam**2 & fbm**2 here   
     1044#            refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
     1045#        else:                                       #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2
     1046#            if len(TwinLaw) > 1:
     1047#                refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
     1048#                refl.T[7] = np.sum(TwinFr*TwMask*np.sum(fas,axis=0)**2,axis=-1)+   \
     1049#                    np.sum(TwinFr*TwMask*np.sum(fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
     1050#                refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f" & use primary twin
     1051#            else:   # checked correct!!
     1052#                refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
     1053#                refl.T[7] = np.copy(refl.T[9])               
     1054#                refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
     1055##                refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
     1056        iBeg += blkSize
     1057#    print 'sf time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
     1058
     1059def MagStructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
    10581060    '''Compute magnetic structure factor derivatives on blocks of reflections - for powders/nontwins only
    10591061    input:
    … …  
    10721074    '''
    10731075    #TODO: fix mag math - moments parallel to crystal axes
     1076    g = nl.inv(G)
    10741077    ast = np.sqrt(np.diag(G))
     1078    ainv = np.sqrt(np.diag(g))
    10751079    GS = G/np.outer(ast,ast)
    1076     uAmat,uBmat = G2lat.Gmat2AB(GS)
     1080    Ginv = g/np.outer(ainv,ainv)
     1081    uAmat = G2lat.Gmat2AB(GS)[0]
    10771082    Mast = twopisq*np.multiply.outer(ast,ast)
    10781083    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
    … …  
    10891094        return {}
    10901095    mSize = len(Mdata)
    1091     Mag = np.sqrt(np.array([np.inner(mag,np.inner(mag,GS)) for mag in Gdata.T]))
    1092     dGdM = np.array([np.inner(mag,GS) for mag in Gdata.T]).T/Mag
    1093    
     1096    Mag = np.array([np.sqrt(np.inner(mag,np.inner(mag,Ginv))) for mag in Gdata.T])
     1097    dMdm = np.inner(Gdata.T,Ginv).T/Mag
    10941098    Gdata = np.inner(Gdata.T,SGMT).T            #apply sym. ops.
     1099    dG = np.inner(np.eye(3),SGMT).T
    10951100    if SGData['SGInv'] and not SGData['SGFixed']:
    10961101        Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
     1102        dG = np.hstack((dG,-dG))
    10971103    Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering
     1104    dG = np.hstack([dG for icen in range(Ncen)])
    10981105    Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip
    10991106    Mag = np.tile(Mag[:,nxs],Nops).T  #make Mag same length as Gdata
    1100     Gdata = Gdata/Mag       #normalze mag. moments
    1101     Gdata = np.inner(Gdata.T,uAmat).T*np.sqrt(nl.det(GS))       #make unit vectors in Cartesian space
    1102 #    dGdm = (1.-Gdata**2)                        #1/Mag removed - canceled out in dqmx=sum(dqdm*dGdm)
     1107    VGi = np.sqrt(nl.det(Ginv))
     1108    Kdata = np.inner(Gdata.T,uAmat).T*VGi/Mag       #make unit vectors in Cartesian space
     1109    dG = np.inner(dG.T,uAmat).T*VGi
     1110    dkdG = (np.inner(np.ones(3),uAmat)*VGi)[:,nxs,nxs]/Mag[nxs,:,:]
     1111    dkdm = dkdG-Kdata*dMdm[:,nxs,:]/Mag[nxs,:,:]
    11031112    dFdMx = np.zeros((nRef,mSize,3))
    11041113    Uij = np.array(G2lat.U6toUij(Uijdata))
    … …  
    11491158        cosm = np.cos(mphase)                               #ditto
    11501159        HM = np.inner(Bmat.T,H)                             #put into cartesian space
    1151         HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #unit vector for H
    1152         eDotK = np.sum(HM[:,:,nxs,nxs]*Gdata[:,nxs,:,:],axis=0)
    1153         Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Gdata[:,nxs,:,:] #Mxyz,Nref,Nop,Natm = BPM in magstrfc.for OK
     1160        HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #unit cartesian vector for H
     1161        eDotK = np.sum(HM[:,:,nxs,nxs]*Kdata[:,nxs,:,:],axis=0)
     1162        deDotK = HM.T[nxs,:,nxs,:]*dG.T[:,nxs,:,:]   #Mxyz,Nref,Nops
     1163        Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Kdata[:,nxs,:,:] #Mxyz,Nref,Nop,Natm = BPM in magstrfc.for OK
     1164        dqdk = np.sum(HM[:,:,nxs,nxs]*deDotK-dG.T[:,nxs,:,:],axis=3)     #Nref,Nops,Mxyyz
    11541165        NQ = np.where(np.abs(Q)>0.,1./np.abs(Q),0.)     #this sort of works esp for 1 axis moments
    1155 #        NQ2 = np.where(np.abs(Q)>0.,1./np.sqrt(np.sum(Q**2,axis=0)),0.)
    1156 #        dqdm = np.array([np.outer(hm,hm)-np.eye(3) for hm in HM.T]).T   #Mxyz,Mxyz,Nref (3x3 matrix)
    1157 #        dqmx = dqdm[:,:,:,nxs,nxs]*dGdm[:,nxs,nxs,:,:]
    1158 #        dqmx2 = np.sum(dqmx,axis=1)   #matrix * vector = vector
    1159 #        dqmx1 = np.swapaxes(np.swapaxes(np.inner(dqdm.T,dGdm.T),0,1),2,3)
    1160         dmx = (NQ*Q*dGdM[:,nxs,nxs,:])        #just use dpdM term & ignore dqdM term(small)
    1161 #        dmx = (NQ*Q*dGdM[:,nxs,:,:]-Q*dqmx2)/Mag
     1166#        dqdk = (HM*HM)-1.
     1167        dqdm = dqdk[:,:,:,nxs]*dkdm[:,nxs,:,:]
     1168        dmx = Q*dMdm[:,nxs,nxs,:]
     1169        dmx += dqdm*Mag[nxs,nxs,:,:]
    11621170       
    11631171        fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #Mxyz,Nref,Nop,Natm
    … …  
    14951503            sinm = np.sin(mphase)                               #ditto - match magstrfc.for
    14961504            cosm = np.cos(mphase)                               #ditto
    1497             HM = np.inner(Bmat.T,H)                             #put into cartesian space
     1505            HM = np.inner(Bmat,H)                             #put into cartesian space
    14981506            HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Gdata = MAGS & HM = UVEC in magstrfc.for both OK
    14991507            eDotK = np.sum(HM[:,:,nxs,nxs]*Gdata[:,nxs,:,:],axis=0)
    … …  
    38323840            else:
    38333841                if Phase['General']['Type'] == 'magnetic':
    3834                     dFdvDict = StructureFactorDervMag(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
     3842                    dFdvDict = MagStructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
    38353843                else:
    38363844                    dFdvDict = StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
    … …  
    41094117        else:   #correct!!
    41104118            if Phase['General']['Type'] == 'magnetic':  #is this going to work for single crystal mag data?
    4111                 dFdvDict = StructureFactorDervMag(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
     4119                dFdvDict = MagStructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
    41124120            else:
    41134121                dFdvDict = StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
Note: See TracChangeset for help on using the changeset viewer.