# Changeset 1635

Ignore:
Timestamp:
Feb 4, 2015 4:38:12 PM (9 years ago)
Message:

modulation functions & constraints

Location:
trunk
Files:
6 edited

Unmodified
Removed
• ## trunk/GSASIImath.py

 r1630 def Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata): import scipy.special as sp for iwt,wt in enumerate(waveTypes):    #atom loop! if wt == 'Fourier': A = np.array([[a,b] for a,b in zip(XSSdata[:3],XSSdata[3:])]) HdotA = twopi*np.inner(A.T,SSUniq.T[:3].T) m = SSUniq.T[3] Gp = sp.jn(-m,HdotA) return np.real(Gp),np.imag(Gp) import scipy.integrate as si m = SSUniq.T[3] nh = np.zeros(1) if XSSdata.ndim > 2: nh = np.arange(XSSdata.shape[1]) M = np.where(m>0,m+nh[:,np.newaxis],m-nh[:,np.newaxis]) A = np.array(XSSdata[:3]) B = np.array(XSSdata[3:]) HdotA = (np.inner(A.T,SSUniq.T[:3].T)+SSPhi) HdotB = (np.inner(B.T,SSUniq.T[:3].T)+SSPhi) GpA = sp.jn(M[:,np.newaxis],twopi*HdotA) GpB = sp.jn(M[:,np.newaxis],twopi*HdotB)*(1.j)**M Gp = np.sum(GpA+GpB,axis=0) return np.real(Gp).T,np.imag(Gp).T def ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata): import scipy.special as sp m = SSUniq.T[3] nh = np.zeros(1) if XSSdata.ndim > 2: nh = np.arange(XSSdata.shape[1]) M = np.where(m>0,m+nh[:,np.newaxis],m-nh[:,np.newaxis]) A = np.array([[a,b] for a,b in zip(XSSdata[:3],XSSdata[3:])]) HdotA = twopi*(np.inner(SSUniq.T[:3].T,A.T)+SSPhi) Gpm = sp.jn(M[:,np.newaxis,:]-1,HdotA) Gpp = sp.jn(M[:,np.newaxis,:]+1,HdotA) if Gpm.ndim > 3: #sum over multiple harmonics Gpm = np.sum(Gpm,axis=0) Gpp = np.sum(Gpp,axis=0) dGpdk = 0.5*(Gpm+Gpp) return np.real(dGpdk),np.imag(dGpdk) def posFourier(tau,psin,pcos):
• ## trunk/GSASIIphsGUI.py

 r1633 waveSizer.Add(waveHead) if len(waveBlk): nFour = 0 for iwave,wave in enumerate(waveBlk): if waveType == 'Fourier': nFour += 1 if not iwave: CSI = G2spc.GetSSfxuinel(waveType,xyz,uij,SGData,SSGData) CSI = G2spc.GetSSfxuinel(waveType,nFour,xyz,SGData,SSGData) else: CSI = G2spc.GetSSfxuinel('Fourier',xyz,uij,SGData,SSGData) CSI = G2spc.GetSSfxuinel('Fourier',nFour,xyz,SGData,SSGData) waveName = 'Fourier' if Stype == 'Sfrac': waveSizer.Add(wx.StaticText(waveData,label=' %s  parameters: %s'%(waveName,str(names).rstrip(']').lstrip('[').replace("'",''))),0,WACV) for ival,val in enumerate(wave[0]): if CSI[Stype][0][ival] == 0: waveVal = wx.TextCtrl(waveData,value='%.4f'%(val),style=wx.TE_READONLY) waveVal.SetBackgroundColour(VERY_LIGHT_GREY) else: if any(CSI[Stype][0][ival]): waveVal = wx.TextCtrl(waveData,value='%.4f'%(val),style=wx.TE_PROCESS_ENTER) waveVal.Bind(wx.EVT_TEXT_ENTER,OnWaveVal) waveVal.Bind(wx.EVT_KILL_FOCUS,OnWaveVal) Indx[waveVal.GetId()] = [iatm,Stype,iwave,ival] else: waveVal = wx.TextCtrl(waveData,value='%.4f'%(val),style=wx.TE_READONLY) waveVal.SetBackgroundColour(VERY_LIGHT_GREY) Waves.Add(waveVal,0,WACV) if len(wave[0]) > 6 and ival == 5:
• ## trunk/GSASIIplot.py

 r1633 Toolbar.__init__(self,plotCanvas) self.plotCanvas = plotCanvas POSITION_OF_CONFIGURE_SUBPLOTS_BTN = 6 # remove one button POSITION_OF_CONFIGURE_SUBPLOTS_BTN = 6 # remove one button, nos. start at 1! self.DeleteToolByPos(POSITION_OF_CONFIGURE_SUBPLOTS_BTN) self.parent = self.GetParent()
• ## trunk/GSASIIspc.py

 r1630 import math import sys import copy import os.path as ospath npsind = lambda x: np.sin(x*np.pi/180.) npcosd = lambda x: np.cos(x*np.pi/180.) DEBUG = True DEBUG = False ################################################################################ return CSuinel[indx[1]] def GetSSfxuinel(XYZ,UIJ,SGData,SSGData): CSI = {'Sfrac':[[1,2],[1.,1.]],'Spos':[[1,2,3, 4,5,6],[1.,1.,1., 1.,1.,1.]],    #sin & cos 'Sadp':[[1,2,3,4,5,6, 7,8,9,10,11,12],[1.,1.,1.,1.,1.,1., 1.,1.,1.,1.,1.,1.]], 'Smag':[[1,2,3, 4,5,6],[1.,1.,1., 1.,1.,1.]]} deltx = np.ones((3,4))*.01 deltx[:3,:3] = np.eye((3))*.001 deltu = np.eye((6))*.0001 def GetSSfxuinel(waveType,nH,XYZ,SGData,SSGData,debug=False): def fracCrenel(tau,Toff,Twid): Tau = (tau-Toff)%1. A = np.where(Tau
• ## trunk/GSASIIstrIO.py

 r1630 waveType = AtomSS['waveType'] phaseDict[pfx+'waveType:'+str(i)] = waveType CSI = G2spc.GetSSfxuinel(at[cx:cx+3],at[cia+2:cia+8],SGData,SSGData) for Stype in ['Sfrac','Spos','Sadp','Smag']: Waves = AtomSS[Stype] uId,uCoef = CSI[Stype] for iw,wave in enumerate(Waves): if not iw: CSI = G2spc.GetSSfxuinel(waveType,at[cx:cx+3],SGData,SSGData) else: CSI = G2spc.GetSSfxuinel('Fourier',at[cx:cx+3],SGData,SSGData) stiw = str(i)+':'+str(iw) if Stype == 'Spos': cx,ct,cs,cia = General['AtomPtrs'] print >>pFile,'\n Modulation waves' #        names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tzero','Xslope','Yslope','Zslope'], #            'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos', #            'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']} #        print >>pFile,135*'-' #        for i,at in enumerate(Atoms): #            AtomSS = at[-1]['SS1'] #            for Stype in ['Sfrac','Spos','Sadp','Smag']: #                Waves = AtomSS[Stype] #                if len(Waves): #                    print >>pFile,' atom: %s, site sym: %s, type: %s wave type: %s:'    \ #                        %(at[ct-1],at[cs],Stype,AtomSS['waveType']) #                    line = '' #                    for item in names[Stype]: #                        line += '%8s '%(item) #                    print >>pFile,line #                for wave in Waves: #                    line = '' #                    for item in wave[0]: #                        line += '%8.4f '%(item) #                    line += ' Refine? '+str(wave[1]) #                    print >>pFile,line names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tzero','Xslope','Yslope','Zslope'], 'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos', 'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']} print >>pFile,135*'-' for i,at in enumerate(Atoms): AtomSS = at[-1]['SS1'] waveType = AtomSS['waveType'] for Stype in ['Sfrac','Spos','Sadp','Smag']: Waves = AtomSS[Stype] if len(Waves): print >>pFile,' atom: %s, site sym: %s, type: %s wave type: %s:'    \ %(at[ct-1],at[cs],Stype,waveType) line = '' for iw,wave in enumerate(Waves): stiw = ':'+str(i)+':'+str(iw) namstr = '  names :' valstr = '  values:' sigstr = '  esds  :' if Stype == 'Spos': nt = 6 ot = 0 if waveType in ['Sawtooth','ZigZag'] and not iw: nt = 4 ot = 6 for j in range(nt): name = names['Spos'][j+ot] namstr += '%12s'%(name) valstr += '%12.4f'%(wave[0][j]) if name+stiw in wavesSig: sigstr += '%12.4f'%(wavesSig[name+stiw]) else: sigstr += 12*' ' elif Stype == 'Sfrac': ot = 0 if 'Crenel' in waveType and not iw: ot = 2 for j in range(2): name = names['Sfrac'][j+ot] namstr += '%12s'%(names['Sfrac'][j+ot]) valstr += '%12.4f'%(wave[0][j]) if name+stiw in wavesSig: sigstr += '%12.4f'%(wavesSig[name+stiw]) else: sigstr += 12*' ' elif Stype == 'Sadp': for j in range(12): name = names['Sadp'][j] namstr += '%10s'%(names['Sadp'][j]) valstr += '%10.6f'%(wave[0][j]) if name+stiw in wavesSig: sigstr += '%10.6f'%(wavesSig[name+stiw]) else: sigstr += 10*' ' elif Stype == 'Smag': for j in range(6): name = names['Smag'][j] namstr += '%12s'%(names['Smag'][j]) valstr += '%12.4f'%(wave[0][j]) if name+stiw in wavesSig: sigstr += '%12.4f'%(wavesSig[name+stiw]) else: sigstr += 12*' ' print >>pFile,namstr print >>pFile,valstr print >>pFile,sigstr if Stype == 'Spos': if waveType in ['ZigZag','Sawtooth'] and not iw: names = [pfx+'Tzero:'+stiw,pfx+'Xslope:'+stiw,pfx+'Yslope:'+stiw,pfx+'Zslope:'+stiw] names = ['Tzero:'+stiw,'Xslope:'+stiw,'Yslope:'+stiw,'Zslope:'+stiw] else: names = [pfx+'Xsin:'+stiw,pfx+'Ysin:'+stiw,pfx+'Zsin:'+stiw, pfx+'Xcos:'+stiw,pfx+'Ycos:'+stiw,pfx+'Zcos:'+stiw] names = ['Xsin:'+stiw,'Ysin:'+stiw,'Zsin:'+stiw, 'Xcos:'+stiw,'Ycos:'+stiw,'Zcos:'+stiw] elif Stype == 'Sadp': names = [pfx+'U11sin:'+stiw,pfx+'U22sin:'+stiw,pfx+'U33sin:'+stiw, pfx+'U12sin:'+stiw,pfx+'U13sin:'+stiw,pfx+'U23sin:'+stiw, pfx+'U11cos:'+stiw,pfx+'U22cos:'+stiw,pfx+'U33cos:'+stiw, pfx+'U12cos:'+stiw,pfx+'U13cos:'+stiw,pfx+'U23cos:'+stiw] names = ['U11sin:'+stiw,'U22sin:'+stiw,'U33sin:'+stiw, 'U12sin:'+stiw,'U13sin:'+stiw,'U23sin:'+stiw, 'U11cos:'+stiw,'U22cos:'+stiw,'U33cos:'+stiw, 'U12cos:'+stiw,'U13cos:'+stiw,'U23cos:'+stiw] elif Stype == 'Sfrac': if 'Crenel' in waveType and not iw: names = [pfx+'Fzero:'+stiw,pfx+'Fwid:'+stiw] names = ['Fzero:'+stiw,'Fwid:'+stiw] else: names = [pfx+'Fsin:'+stiw,pfx+'Fcos:'+stiw] names = ['Fsin:'+stiw,'Fcos:'+stiw] elif Stype == 'Smag': names = [pfx+'MXsin:'+stiw,pfx+'MYsin:'+stiw,pfx+'MZsin:'+stiw, pfx+'MXcos:'+stiw,pfx+'MYcos:'+stiw,pfx+'MZcos:'+stiw] names = ['MXsin:'+stiw,'MYsin:'+stiw,'MZsin:'+stiw, 'MXcos:'+stiw,'MYcos:'+stiw,'MZcos:'+stiw] for iname,name in enumerate(names): AtomSS[Stype][iw][0][iname] = parmDict[name] if name in sigDict: wavesSig[name] = sigDict[name] AtomSS[Stype][iw][0][iname] = parmDict[pfx+name] if pfx+name in sigDict: wavesSig[name] = sigDict[pfx+name] PrintAtomsAndSig(General,Atoms,atomsSig)
• ## trunk/GSASIIstrMath.py

 r1630 Tcorr = Tiso*Tuij*Mdata*Fdata/len(Uniq) fa = np.array([(FF+FP-Bab)*cosp*Tcorr,-FPP*sinp*Tcorr]) fb = 0. fb = np.zeros_like(fa) if not SGData['SGInv']: fb = np.array([(FF+FP-Bab)*sinp*Tcorr,FPP*cosp*Tcorr]) fa = fa*GfpuA.T-fb*GfpuB.T fb = fb*GfpuA.T+fa*GfpuB.T fas = np.sum(np.sum(fa,axis=1),axis=1)        #real fbs = np.sum(np.sum(fb,axis=1),axis=1) fa = fa*GfpuA-fb*GfpuB fb = fb*GfpuA+fa*GfpuB fas = np.real(np.sum(np.sum(fa,axis=1),axis=1))        #real fbs = np.real(np.sum(np.sum(fb,axis=1),axis=1)) fasq = fas**2 FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im]) H = np.array(refl[:4]) if H[3]: continue SQ = 1./(2.*refl[4+im])**2             # or (sin(theta)/lambda)**2 SQfactor = 8.0*SQ*np.pi**2 Phi = np.inner(H[:3],SGT) SSPhi = np.inner(H,SSGT) GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata) GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata) dGAdk,dGBdk = G2mth.ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata) #need ModulationDerv here dGAdXsin, etc phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+Phi[np.newaxis,:]) sinp = np.sin(phase) biso = -SQfactor*Uisodata Tiso = np.where(biso<1.,np.exp(biso),1.0) HbH = -np.inner(H,np.inner(bij,H)) HbH = -np.inner(H[:3],np.inner(bij,H[:3])) Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq]) Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij]) fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp]) GfpuA = np.swapaxes(GfpuA,1,2) GfpuB = np.swapaxes(GfpuB,1,2) fa = fa*GfpuA-fb*GfpuB fb = fb*GfpuA+fa*GfpuB fas = np.sum(np.sum(fa,axis=1),axis=1) fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp]) fax = fax*GfpuA-fbx*GfpuB fbx = fbx*GfpuA+fax*GfpuB #sum below is over Uniq dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1]) dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T #need dFdXsin, etc for modulations... if not SGData['SGInv']: dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1]) dFdbab[iref] += 2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T #need dFdXsin, etc for modulations... #loop over atoms - each dict entry is list of derivatives for all the reflections for i in range(len(Mdata)): dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i] dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i] #need dFdvDict[pfx+'Xsin:'+str[i]:str(m)], etc for modulations... dFdvDict[pfx+'BabA'] = dFdbab.T[0] dFdvDict[pfx+'BabU'] = dFdbab.T[1]
Note: See TracChangeset for help on using the changeset viewer.