Changeset 1966
- Timestamp:
- Aug 28, 2015 3:02:12 PM (7 years ago)
- Location:
- trunk
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIgrid.py
r1958 r1966 2831 2831 Hmax = np.array([int(np.max(refList.T[0])),int(np.max(refList.T[1])),int(np.max(refList.T[2]))]) 2832 2832 Vpoint = [int(np.mean(refList.T[0])),int(np.mean(refList.T[1])),int(np.mean(refList.T[2]))] 2833 controls = {'Type' : 'Fosq','Iscale' : False,'HKLmax' : Hmax,'HKLmin' : Hmin, 2833 controls = {'Type' : 'Fosq','Iscale' : False,'HKLmax' : Hmax,'HKLmin' : Hmin,'Zone':False, 2834 2834 'FoMax' : FoMax,'Scale' : 1.0,'Drawing':{'viewPoint':[Vpoint,[]],'default':Vpoint[:], 2835 'backColor':[0,0,0],'depthFog':False,'Zclip':10.0,'cameraPos':10.,'Zstep':0.05, 2836 'Scale':1.0,'oldxy':[],'viewDir':[ 1,0,0]},'Super':Super,'SuperVec':SuperVec}2835 'backColor':[0,0,0],'depthFog':False,'Zclip':10.0,'cameraPos':10.,'Zstep':0.05,'viewUp':[0,1,0], 2836 'Scale':1.0,'oldxy':[],'viewDir':[0,0,1]},'Super':Super,'SuperVec':SuperVec} 2837 2837 G2plt.Plot3DSngl(G2frame,newPlot=True,Data=controls,hklRef=refList,Title=phaseName) 2838 2838 … … 2861 2861 Hmax = np.array([int(np.max(refList.T[0])),int(np.max(refList.T[1])),int(np.max(refList.T[2]))]) 2862 2862 Vpoint = [int(np.mean(refList.T[0])),int(np.mean(refList.T[1])),int(np.mean(refList.T[2]))] 2863 controls = {'Type' : 'Fosq','Iscale' : False,'HKLmax' : Hmax,'HKLmin' : Hmin, 2863 controls = {'Type' : 'Fosq','Iscale' : False,'HKLmax' : Hmax,'HKLmin' : Hmin,'Zone':False, 2864 2864 'FoMax' : FoMax,'Scale' : 1.0,'Drawing':{'viewPoint':[Vpoint,[]],'default':Vpoint[:], 2865 'backColor':[0,0,0],'depthFog':False,'Zclip':10.0,'cameraPos':10.,'Zstep':0.05, 2865 'backColor':[0,0,0],'depthFog':False,'Zclip':10.0,'cameraPos':10.,'Zstep':0.05,'viewUp':[0,1,0], 2866 2866 'Scale':1.0,'oldxy':[],'viewDir':[1,0,0]},'Super':Super,'SuperVec':SuperVec} 2867 2867 G2plt.Plot3DSngl(G2frame,newPlot=True,Data=controls,hklRef=refList,Title=phaseName) -
trunk/GSASIImath.py
r1958 r1966 924 924 ################################################################################ 925 925 926 def Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM ):926 def Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM,Mast): 927 927 import scipy.special as sp 928 928 import scipy.integrate as si 929 930 def cosMod(tau,H,A,B,S): 931 uTau = posFourier(tau,A,B,S) 932 SdotU = twopi*(np.inner(uTau.T,H[:3].T))+H[3]*tau 933 return list(np.cos(SdotU))[0] 934 935 def sinMod(tau,H,A,B,S): 936 uTau = posFourier(tau,A,B,S) 937 SdotU = twopi*(np.inner(uTau.T,H[:3]))+H[3]*tau 938 return list(np.sin(SdotU))[0] 939 940 def expModInt(H,A,B,S): 941 intReal = np.array([[si.romberg(cosMod,0,1,args=(h,a,b,s),rtol=1.e-1) for a,b in zip(A,B)] for h,s in zip(H,S)]) 942 intImag = np.array([[si.romberg(sinMod,0,1,args=(h,a,b,s),rtol=1.e-1) for a,b in zip(A,B)] for h,s in zip(H,S)]) 943 # return np.sqrt(intReal**2+intImag**2) 944 return intReal,intImag 945 929 946 Smult,TauT = SStauM # both atoms x SGops 930 947 m = SSUniq.T[3] 931 A = np.array(XSSdata[:3]) #sin mods x waves x atoms 932 B = np.array(XSSdata[3:]) #cos mods... 933 if XSSdata.ndim > 2: 934 nh = np.arange(XSSdata.shape[1]) #0..max harmonic order-1 935 M = np.where(m>0,m+nh[:,np.newaxis],m-nh[:,np.newaxis]) # waves x SGops 936 HdotA = (np.inner(A.T,SSUniq.T[:3].T)) #atoms X waves x SGops 937 HdotB = (np.inner(B.T,SSUniq.T[:3].T)) # ditto 938 Hdot = np.sqrt(HdotA**2+HdotB**2) 939 Hphi = np.sum(np.arctan2(HdotB,HdotA)*M,axis=1) 940 GpA = np.sum(sp.jn(-M,twopi*Hdot),axis=1) # sum on waves 948 Ax = np.array(XSSdata[:3]).T #atoms x waves x sin pos mods 949 Bx = np.array(XSSdata[3:]).T #...cos pos mods 950 Af = np.array(FSSdata[0]).T #sin frac mods x waves x atoms 951 Bf = np.array(FSSdata[1]).T #cos frac mods... 952 Ab = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T #atoms x waves x sin Uij mods 953 Bb = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T #...cos Uij mods 954 # GSASIIpath.IPyBreak() 955 if Ax.ndim > 2: 956 GpA = np.array(expModInt(SSUniq,Ax,Bx,Smult)) 941 957 else: 942 HdotA = np.inner(A.T,SSUniq.T[:3].T) #atoms X SGops 943 HdotB = np.inner(B.T,SSUniq.T[:3].T) # ditto 944 Hdot = np.sqrt(HdotA**2+HdotB**2) 945 Hphi = np.arctan2(HdotB,HdotA) 946 GpA = sp.jn(-m,twopi*Hdot) # atoms x SGops 947 948 # GSASIIpath.IPyBreak() 949 return GpA.T,Hphi.T # SGops x atoms 950 951 def ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM): 958 GpA = np.array(expModInt(SSUniq,Ax[:,np.newaxis,:],Bx[:,np.newaxis,:],Smult)) 959 return GpA # SGops x atoms 960 961 def ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM,Mast): 952 962 import scipy.special as sp 953 963 Smult,TauT = SStauM # both atoms x SGops … … 1063 1073 return drawAtoms,Fade 1064 1074 1075 def BessJn(nmax,x): 1076 ''' compute Bessel function J(n,x) from scipy routine & recurrance relation 1077 returns sequence of J(n,x) for n in range [-nmax...0...nmax] 1078 1079 :param integer nmax: maximul order for Jn(x) 1080 :param float x: argument for Jn(x) 1081 1082 :returns numpy array: [J(-nmax,x)...J(0,x)...J(nmax,x)] 1083 1084 ''' 1085 import scipy.special as sp 1086 bessJn = np.zeros(2*nmax+1) 1087 bessJn[nmax] = sp.j0(x) 1088 bessJn[nmax+1] = sp.j1(x) 1089 bessJn[nmax-1] = -bessJn[nmax+1] 1090 for i in range(2,nmax+1): 1091 bessJn[i+nmax] = 2*(i-1)*bessJn[nmax+i-1]/x-bessJn[nmax+i-2] 1092 bessJn[nmax-i] = bessJn[i+nmax]*(-1)**i 1093 return bessJn 1094 1095 def BessIn(nmax,x): 1096 ''' compute modified Bessel function I(n,x) from scipy routines & recurrance relation 1097 returns sequence of I(n,x) for n in range [-nmax...0...nmax] 1098 1099 :param integer nmax: maximul order for In(x) 1100 :param float x: argument for In(x) 1101 1102 :returns numpy array: [I(-nmax,x)...I(0,x)...I(nmax,x)] 1103 1104 ''' 1105 import scipy.special as sp 1106 bessIn = np.zeros(2*nmax+1) 1107 bessIn[nmax] = sp.i0(x) 1108 bessIn[nmax+1] = sp.i1(x) 1109 bessIn[nmax-1] = bessIn[nmax+1] 1110 for i in range(2,nmax+1): 1111 bessIn[i+nmax] = bessIn[nmax+i-2]-2*(i-1)*bessIn[nmax+i-1]/x 1112 bessIn[nmax-i] = bessIn[i+nmax] 1113 return bessIn 1114 1065 1115 1066 1116 ################################################################################ -
trunk/GSASIIplot.py
r1958 r1966 572 572 global ifBox 573 573 Choice = {'F':'Fo','S':'Fosq','U':'Unit','D':'dFsq','W':'dFsq/sig'} 574 viewChoice = {'H':[[1,0,0],[0,0,1]],'K':[[0,1,0],[1,0,0]],'L':[[0,0,1],[0,1,0]]} 574 575 try: 575 576 keyCode = event.GetKeyCode() … … 579 580 except AttributeError: #if from OnKeyBox above 580 581 key = str(event.key).upper() 581 if key in ['C']: 582 if key in ['C','H','K','L']: 583 if key == 'C': 584 Data['Zone'] = False 585 key = 'L' 582 586 drawingData['viewPoint'][0] = drawingData['default'] 583 drawingData['viewDir'] = [0,0,1] 587 drawingData['viewDir'] = viewChoice[key][0] 588 drawingData['viewUp'] = viewChoice[key][1] 584 589 drawingData['oldxy'] = [] 585 V0 = np.array([0,0,1]) 586 V = np.inner(Amat,V0) 587 V /= np.sqrt(np.sum(V**2)) 588 A = np.arccos(np.sum(V*V0)) 589 Q = G2mth.AV2Q(A,[0,1,0]) 590 drawingData['Quaternion'] = Q 591 Q = drawingData['Quaternion'] 590 # V0 = np.array(viewChoice[key][0]) 591 # V = np.inner(Amat,V0) 592 # V /= np.sqrt(np.sum(V**2)) 593 # A = np.arccos(np.sum(V*V0)) 594 # Q = G2mth.AV2Q(A,viewChoice[key][1]) 595 drawingData['Quaternion'] = [-1,0,0,0] #the old stuff always gave this... 596 elif key in 'Z': 597 Data['Zone'] = not Data['Zone'] 592 598 elif key in 'B': 593 599 ifBox = not ifBox … … 883 889 cPos = drawingData['cameraPos'] 884 890 Zclip = drawingData['Zclip']*cPos/20. 891 if Data['Zone']: 892 Zclip = 0.1 885 893 Q = drawingData['Quaternion'] 886 894 Tx,Ty,Tz = drawingData['viewPoint'][0] … … 901 909 glViewport(0,0,VS[0],VS[1]) 902 910 gluPerspective(20.,aspect,cPos-Zclip,cPos+Zclip) 903 gluLookAt(0,0,cPos,0,0,0,0,1,0) 911 vDir = drawingData['viewDir'] 912 vUp = drawingData['viewUp'] 913 # gluLookAt(0,0,cPos,0,0,0,0,1,0) 914 gluLookAt(vDir[0]*cPos,vDir[1]*cPos,vDir[2]*cPos, 0,0,0, vUp[0],vUp[1],vUp[2]) 904 915 SetLights() 905 916 … … 936 947 Page.SetFocus() 937 948 Page.Choice = None 938 choice = [' save as/key:','jpeg','tiff','bmp','c: recenter to default','b: toggle box ','+: increase scale', 939 '-: decrease scale','f: Fobs','s: Fobs**2','u: unit','d: Fo-Fc','w: DF/sig','i: toggle intensity scaling'] 949 choice = [' save as/key:','jpeg','tiff','bmp','h: view down h','k: view down k','l: view down l', 950 'z: zero zone toggle','c: reset to default','b: toggle box ','+: increase scale','-: decrease scale', 951 'f: Fobs','s: Fobs**2','u: unit','d: Fo-Fc','w: DF/sig','i: toggle intensity scaling'] 940 952 cb = wx.ComboBox(G2frame.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,choices=choice) 941 953 cb.Bind(wx.EVT_COMBOBOX, OnKeyBox) -
trunk/GSASIIpwdGUI.py
r1959 r1966 3121 3121 Hmax = np.array([int(np.max(refList.T[0])),int(np.max(refList.T[1])),int(np.max(refList.T[2]))]) 3122 3122 Vpoint = [int(np.mean(refList.T[0])),int(np.mean(refList.T[1])),int(np.mean(refList.T[2]))] 3123 controls = {'Type':'Fosq','Iscale':False,'HKLmax':Hmax,'HKLmin':Hmin, 3123 controls = {'Type':'Fosq','Iscale':False,'HKLmax':Hmax,'HKLmin':Hmin,'Zone':False, 3124 3124 'FoMax' : FoMax,'Scale' : 1.0,'Drawing':{'viewPoint':[Vpoint,[]],'default':Vpoint[:], 3125 'backColor':[0,0,0],'depthFog':False,'Zclip':10.0,'cameraPos':10.,'Zstep':0.05, 3126 'Scale':1.0,'oldxy':[],'viewDir':[ 1,0,0]},'Super':Super,'SuperVec':SuperVec}3125 'backColor':[0,0,0],'depthFog':False,'Zclip':10.0,'cameraPos':10.,'Zstep':0.05,'viewUp':[0,1,0], 3126 'Scale':1.0,'oldxy':[],'viewDir':[0,0,1]},'Super':Super,'SuperVec':SuperVec} 3127 3127 G2plt.Plot3DSngl(G2frame,newPlot=True,Data=controls,hklRef=refList,Title=phaseName) 3128 3128 -
trunk/GSASIIstrMath.py
r1958 r1966 814 814 if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict: 815 815 Flack = 1.-2.*parmDict[phfx+'Flack'] 816 time0 = time.time() 817 nref = len(refDict['RefList'])/100 816 818 for iref,refl in enumerate(refDict['RefList']): 817 819 if 'T' in calcControls[hfx+'histType']: … … 915 917 dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \ 916 918 2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T 919 print ' derivative time %.4f\r'%(time.time()-time0) 917 920 #loop over atoms - each dict entry is list of derivatives for all the reflections 918 921 if nTwin > 1: … … 998 1001 dat = G2el.getFFvalues(FFtables,0.) 999 1002 refDict['FF']['El'] = dat.keys() 1000 refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat))) 1003 refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat))) 1004 time0 = time.time() 1005 nref = len(refDict['RefList'])/100 1001 1006 for iref,refl in enumerate(refDict['RefList']): 1002 1007 if 'NT' in calcControls[hfx+'histType']: … … 1028 1033 Phi = np.hstack((Phi,-Phi)) 1029 1034 SSPhi = np.hstack((SSPhi,-SSPhi)) 1030 GfpuA ,Hphi = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM)1035 GfpuA = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM,Mast) 1031 1036 phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+Phi[:,np.newaxis]) 1032 1037 sinp = np.sin(phase) … … 1040 1045 fb = np.array([(FF+FP-Bab)*sinp*Tcorr,FPP*cosp*Tcorr]) 1041 1046 fa *= GfpuA 1042 fb *= GfpuA 1047 fb *= GfpuA 1043 1048 fas = np.real(np.sum(np.sum(fa,axis=1),axis=1)) 1044 1049 fbs = np.real(np.sum(np.sum(fb,axis=1),axis=1)) … … 1048 1053 refl[7+im] = np.sum(fasq)+np.sum(fbsq) 1049 1054 refl[10+im] = atan2d(fbs[0],fas[0]) 1055 if not iref%nref: print 'ref no. %d time %.4f\r'%(iref,time.time()-time0), 1056 print '\nref no. %d time %.4f\r'%(iref,time.time()-time0) 1050 1057 # GSASIIpath.IPyBreak() 1051 1058 … … 1094 1101 Phi = np.inner(H[:3],SGT) 1095 1102 SSPhi = np.inner(H,SSGT) 1096 GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM )1097 dGAdk,dGBdk = G2mth.ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM )1103 GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM,Mast) 1104 dGAdk,dGBdk = G2mth.ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM,Mast) 1098 1105 #need ModulationDerv here dGAdXsin, etc 1099 1106 phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+Phi[np.newaxis,:])
Note: See TracChangeset
for help on using the changeset viewer.