Changeset 303


Ignore:
Timestamp:
Jun 20, 2011 2:37:07 PM (12 years ago)
Author:
vondreele
Message:

Allow reading of xye (Topas style) files
Begin implementation of spherical harmonics texture
Refactor indexing; replace cell refinement from peak positions

Location:
trunk
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASII.py

    r295 r303  
    319319        self.CheckNotebook()
    320320        dlg = wx.FileDialog(self, 'Choose files', '.', '',
    321             'GSAS fxye files (*.fxye)|*.fxye|GSAS fxy files (*.fxy)|*.fxy|All files (*.*)|*.*',
     321            'GSAS fxye files (*.fxye)|*.fxye|GSAS fxy files (*.fxy)|*.fxy|Topas xye files (*.xye)|*.xye|All files (*.*)|*.*',
    322322            wx.OPEN | wx.MULTIPLE)
    323323        if self.dirname: dlg.SetDirectory(self.dirname)
    … …  
    875875            {'General':{'Name':PhaseName,'Type':'nuclear','SGData':SGData,
    876876            'Cell':[False,10.,10.,10.,90.,90.,90,1000.],
    877             'Pawley dmin':1.0},'Atoms':[],'Drawing':{},'Histograms':{},'Pawley ref':[],'Models':{}})
     877            'Pawley dmin':1.0},'Atoms':[],'Drawing':{},'Histograms':{},'Pawley ref':[],
     878            'Models':{},'SH Texture':{'Order':0,'Model':'cylindrical','Sample omega':[False,0.0],
     879            'Sample chi':[False,0.0],'Sample phi':[False,0.0],'SH Coeff':[False,{}],
     880            'SHShow':False,'PFhkl':[0,0,1]}})
    878881       
    879882    def OnDeletePhase(self,event):

  • trunk/GSASIIIO.py

    r296 r303  
    4747'''+File.readline()
    4848    dlg = wx.MessageDialog(self, Title, 'Is this the file you want?',
    49     wx.YES_NO | wx.ICON_QUESTION)
     49        wx.YES_NO | wx.ICON_QUESTION)
    5050    try:
    5151        result = dlg.ShowModal()
    … …  
    5454    if result == wx.ID_NO: return (0,0)
    5555    Temperature = 300
     56   
     57    if '.xye' in filename:      #Topas style xye file (e.g. 2-th, I, sig) - no iparm file/no BANK record
     58        dlg = wx.MessageDialog(self,'''Is this laboratory Cu Ka1/Ka2 data?
     59(No = 0.6A wavelength synchrotron data)
     60Change wavelength in Instrument Parameters if needed''','Data type?',
     61            wx.YES_NO | wx.ICON_QUESTION)
     62        try:
     63            result = dlg.ShowModal()
     64        finally:
     65            dlg.Destroy()
     66        print result
     67        if result == wx.ID_YES:
     68            Iparm = {}                                               #Assume CuKa lab data
     69            Iparm['INS   HTYPE '] = 'PXC '
     70            Iparm['INS  1 ICONS'] = '  1.540500  1.544300       0.0         0       0.7    0       0.5   '
     71            Iparm['INS  1PRCF1 '] = '    3    8      0.01                                                '
     72            Iparm['INS  1PRCF11'] = '   2.000000E+00  -2.000000E+00   5.000000E+00   0.000000E+00        '
     73            Iparm['INS  1PRCF12'] = '   0.000000E+00   0.000000E+00   0.150000E-01   0.150000E-01        '
     74        else:
     75            Iparm = {}                                               #Assume 0.6A synchrotron data
     76            Iparm['INS   HTYPE '] = 'PXC '
     77            Iparm['INS  1 ICONS'] = '  0.600000  0.000000       0.0         0      0.99    0       0.5   '
     78            Iparm['INS  1PRCF1 '] = '    3    8      0.01                                                '
     79            Iparm['INS  1PRCF11'] = '   1.000000E+00  -1.000000E+00   0.300000E+00   0.000000E+00        '
     80            Iparm['INS  1PRCF12'] = '   0.000000E+00   0.000000E+00   0.100000E-01   0.100000E-01        '
     81                       
    5682       
    57     self.IparmName = GetInstrumentFile(self,filename)
    58     if self.IparmName:
    59         Iparm = GetInstrumentData(self.IparmName)
    60     else:
    61         Iparm = {}                                               #Assume CuKa lab data if no iparm file
    62         Iparm['INS   HTYPE '] = 'PXC '
    63         Iparm['INS  1 ICONS'] = '  1.540500  1.544300       0.0         0       0.7    0       0.5   '
    64         Iparm['INS  1PRCF1 '] = '    3    8      0.01                                                '
    65         Iparm['INS  1PRCF11'] = '   2.000000E+00  -2.000000E+00   5.000000E+00   0.000000E+00        '
    66         Iparm['INS  1PRCF12'] = '   0.000000E+00   0.000000E+00   0.150000E-01   0.150000E-01        '
     83    else:                       #GSAS style fxye or fxy file (e.g. 100*2-th, I, sig)
     84        self.IparmName = GetInstrumentFile(self,filename)
     85        if self.IparmName:
     86            Iparm = GetInstrumentData(self.IparmName)
     87        else:
     88            Iparm = {}                                               #Assume CuKa lab data if no iparm file
     89            Iparm['INS   HTYPE '] = 'PXC '
     90            Iparm['INS  1 ICONS'] = '  1.540500  1.544300       0.0         0       0.7    0       0.5   '
     91            Iparm['INS  1PRCF1 '] = '    3    8      0.01                                                '
     92            Iparm['INS  1PRCF11'] = '   2.000000E+00  -2.000000E+00   5.000000E+00   0.000000E+00        '
     93            Iparm['INS  1PRCF12'] = '   0.000000E+00   0.000000E+00   0.150000E-01   0.150000E-01        '
    6794    S = 1
    6895    Banks = []
    … …  
    78105                    Banks.append(S)
    79106                    Pos.append(File.tell())
     107                elif '.xye' in filename:    #No BANK in a xye file
     108                    Banks.append('BANK 1 XYE')
     109                    Pos.append(File.tell())
     110                    break
    80111            else:
    81112                Comments.append(S[:-1])
    … …  
    283314    if 'FXYE' in Bank:
    284315        return GetFXYEdata(filename,Pos,Bank,DataType)
     316    elif ' XYE' in Bank:
     317        return GetXYEdata(filename,Pos,Bank,DataType)
    285318    elif 'FXY' in Bank:
    286319        return GetFXYdata(filename,Pos,Bank,DataType)
    … …  
    317350    N = len(x)
    318351    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
     352   
     353def GetXYEdata(filename,Pos,Bank,DataType):
     354    File = open(filename,'Ur')
     355    File.seek(Pos)
     356    x = []
     357    y = []
     358    w = []
     359    S = File.readline()
     360    while S:
     361        vals = S.split()
     362        try:
     363            x.append(float(vals[0]))
     364            f = float(vals[1])
     365            if f <= 0.0:
     366                y.append(0.0)
     367                w.append(1.0)
     368            else:
     369                y.append(float(vals[1]))
     370                w.append(1.0/float(vals[2])**2)
     371            S = File.readline()
     372        except ValueError:
     373            break
     374    File.close()
     375    N = len(x)
     376    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
     377   
    319378   
    320379def GetFXYdata(filename,Pos,Bank,DataType):
    … …  
    9551014   
    9561015def ReadEXPPhase(self,filename):
     1016    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
     1017    textureData = {'Order':0,'Model':'cylindrical','Sample omega':[False,0.0],
     1018        'Sample chi':[False,0.0],'Sample phi':[False,0.0],'SH Coeff':[False,{}],
     1019        'SHShow':False,'PFhkl':[0,0,1]}
     1020    shNcof = 0
    9571021    file = open(filename, 'Ur')
    9581022    Phase = {}
    … …  
    10031067            SpGrp = EXPphase[key][:15].strip()
    10041068            E,SGData = G2spc.SpcGroup(SpGrp)
     1069        elif 'OD    ' in key:
     1070            SHdata = EXPphase[key].split()
     1071            textureData['Order'] = int(SHdata[0])
     1072            textureData['Model'] = shModels[int(SHdata[2])]
     1073            textureData['Sample omega'] = [False,float(SHdata[6])]
     1074            textureData['Sample chi'] = [False,float(SHdata[7])]
     1075            textureData['Sample phi'] = [False,float(SHdata[8])]
     1076            shNcof = int(SHdata[1])
    10051077    Atoms = []
    10061078    if Ptype == 'nuclear':
    … …  
    10371109                Atoms.append(Atom)
    10381110    Volume = G2lat.calc_V(G2lat.cell2A(abc+angles))
     1111    if shNcof:
     1112        shCoef = {}
     1113        nRec = [i+1 for i in range((shNcof-1)/6+1)]
     1114        for irec in nRec:
     1115            ODkey = keyList[0][:6]+'OD'+'%3dA'%(irec)
     1116            indx = EXPphase[ODkey].split()
     1117            ODkey = ODkey[:-1]+'B'
     1118            vals = EXPphase[ODkey].split()
     1119            for i,val in enumerate(vals):
     1120                key = 'C(%s,%s,%s)'%(indx[3*i],indx[3*i+1],indx[3*i+2])
     1121                shCoef[key] = float(val)
     1122        textureData['SH Coeff'] = [False,shCoef]
     1123           
    10391124    Phase['General'] = {'Name':PhaseName,'Type':Ptype,'SGData':SGData,
    1040         'Cell':[False,]+abc+angles+[Volume,]}
     1125        'Cell':[False,]+abc+angles+[Volume,],'Pawley ref':[],'Models':{},'SH Texture':textureData}
    10411126    Phase['Atoms'] = Atoms
    10421127    Phase['Drawing'] = {}

  • trunk/GSASIIindex.py

    r302 r303  
    1616import pypowder as pyp              #assumes path has been amended to include correctr bin directory
    1717import GSASIIlattice as G2lat
     18import scipy.optimize as so
    1819
    1920# trig functions in degrees
    … …  
    194195def IndexPeaks(peaks,HKL):
    195196    import bisect
    196     hklds = [1000.0]                                    #make bounded list of available d-spacings
    197197    N = len(HKL)
    198198    if N == 0: return False
    199     for hkl in HKL:
    200         hklds.append(hkl[3])
    201     hklds.append(0.0)
     199    hklds = list(np.array(HKL).T[3])+[1000.0,0.0,]
    202200    hklds.sort()                                        # ascending sort - upper bound at end
    203201    hklmax = [0,0,0]
    … …  
    215213                dnew = min(dm,dp)
    216214                if dold > dnew:                             # new better - zero out old
    217                     for j in range(3):
    218                         opeak[j+4] = 0
     215                    opeak[4:7] = [0,0,0]
    219216                    opeak[8] = 0.
    220217                else:                                       # old better - do nothing
    221218                    continue               
    222219            hkl[4] = ipk
    223             for j in range(3):
    224                 peak[j+4] = hkl[j]
     220            peak[4:7] = hkl[:3]
    225221            peak[8] = hkl[3]                                # fill in d-calc
    226222    for peak in peaks:
    … …  
    235231    else:
    236232        return False
    237    
    238 def FitHKL(ibrav,peaks,A,wtP):
    239     def ShiftTest(a,b):
    240         if b < -0.1*a:
    241             b = -0.0001*a
    242         return b
    243     smin = 0.
    244     first = True
    245     for peak in peaks:
    246         if peak[2] and peak[3]:
    247             h,k,l = H = peak[4:7]
    248             Qo = 1./peak[7]**2
    249             Qc = G2lat.calc_rDsq(H,A)
    250             try:
    251                 peak[8] = 1./math.sqrt(Qc)
    252             except:
    253                 print G2lat.A2invcell(A)
    254             delt = Qo-Qc
    255             smin += delt**2
    256             dp = []
    257             if ibrav in [0,1,2]:            #m3m
    258                 dp.append(h*h+k*k+l*l)
    259             elif ibrav in [3,4]:            #R3H, P3/m & P6/mmm
    260                 dp.append(h*h+k*k+h*k)
    261                 dp.append(l*l)
    262             elif ibrav in [5,6]:            #4/mmm
    263                 dp.append(h*h+k*k)
    264                 dp.append(l*l)
    265             elif ibrav in [7,8,9,10]:       #mmm
    266                 dp.append(h*h)
    267                 dp.append(k*k)
    268                 dp.append(l*l)
    269             elif ibrav in [11,12]:          #2/m
    270                 dp.append(h*h)
    271                 dp.append(k*k)
    272                 dp.append(l*l)
    273                 dp.append(h*l)
    274             else:                           #1
    275 #    derivatives for a0*h^2+a1*k^2+a2*l^2+a3*h*k+a4*h*l+a5*k*l
    276                 dp.append(h*h)
    277                 dp.append(k*k)
    278                 dp.append(l*l)
    279                 dp.append(h*k)
    280                 dp.append(h*l)
    281                 dp.append(k*l)
    282             if first:
    283                 first = False
    284                 M = len(dp)
    285                 B = np.zeros(shape=(M,M))
    286                 V = np.zeros(shape=(M))
    287             dwt = peak[7]**wtP
    288             B,V = pyp.buildmv(delt*dwt,dwt,M,dp,B,V)
    289     if nl.det(B) > 0.0:
    290         try:
    291             b = nl.solve(B,V)
    292             B = nl.inv(B)
    293             sig = np.diag(B)
    294         except SingularMatrix:
    295             return False,0
    296         if ibrav in [0,1,2]:            #m3m
    297             A[0] += ShiftTest(A[0],b[0])
    298             A[1] = A[2] = A[0]
    299         elif ibrav in [3,4]:            #R3H, P3/m & P6/mmm
    300             A[0] += ShiftTest(A[0],b[0])
    301             A[1] = A[3] = A[0]
    302             A[2] += ShiftTest(A[2],b[1])
    303         elif ibrav in [5,6]:            #4/mmm
    304             A[0] += ShiftTest(A[0],b[0])
    305             A[1] = A[0]
    306             A[2] += ShiftTest(A[2],b[1])
    307         elif ibrav in [7,8,9,10]:       #mmm
    308             for i in range(3):
    309                 A[i] += ShiftTest(A[i],b[i])
    310         elif ibrav in [11,12]:          #2/m
    311             for i in range(3):
    312                 A[i] += ShiftTest(A[i],b[i])
    313             A[4] += ShiftTest(A[4],b[3])
    314             A[4] = min(1.4*math.sqrt(A[0]*A[2]),A[4])   #min beta star = 45
    315         else:                           #1
    316             oldV = math.sqrt(1./G2lat.calc_rVsq(A))
    317             oldA = A[:]
    318             for i in range(6):
    319                 A[i] += b[i]*0.2
    320             A[3] = min(1.1*math.sqrt(max(0,A[1]*A[2])),A[3])
    321             A[4] = min(1.1*math.sqrt(max(0,A[0]*A[2])),A[4])
    322             A[5] = min(1.1*math.sqrt(max(0,A[0]*A[1])),A[5])
    323             ratio = math.sqrt(1./G2lat.calc_rVsq(A))/oldV
    324             if 0.9 > ratio or ratio > 1.1:
    325                 A = oldA
    326 #    else:
    327 #        print B
    328 #        print V
    329 #        for peak in peaks:
    330 #            print peak
    331     return True,smin
    332        
     233
     234def FitHKL(ibrav,peaks,A,Pwr):
     235   
     236    def Values2A(ibrav,values):
     237        if ibrav in [0,1,2]:
     238            return [values[0],values[0],values[0],0,0,0]
     239        elif ibrav in [3,4]:
     240            return [values[0],values[0],values[1],values[0],0,0]
     241        elif ibrav in [5,6]:
     242            return [values[0],values[0],values[1],0,0,0]
     243        elif ibrav in [7,8,9,10]:
     244            return [values[0],values[1],values[2],0,0,0]
     245        elif ibrav in [11,12]:
     246            return [values[0],values[1],values[2],0,values[3],0]
     247        else:
     248            return values
     249           
     250    def A2values(ibrav,A):
     251        if ibrav in [0,1,2]:
     252            return [A[0],]
     253        elif ibrav in [3,4,5,6]:
     254            return [A[0],A[2]]
     255        elif ibrav in [7,8,9,10]:
     256            return [A[0],A[1],A[2]]
     257        elif ibrav in [11,12]:
     258            return [A[0],A[1],A[2],A[4]]
     259        else:
     260            return A
     261   
     262    def errFit(values,ibrav,d,H,Pwr):
     263        A = Values2A(ibrav,values)
     264        Qo = 1./d**2
     265        Qc = G2lat.calc_rDsq(H,A)
     266        return (Qo-Qc)*d**Pwr
     267   
     268    Peaks = np.array(peaks).T
     269    values = A2values(ibrav,A)   
     270    result = so.leastsq(errFit,values,args=(ibrav,Peaks[7],Peaks[4:7],Pwr),full_output=True)
     271    A = Values2A(ibrav,result[0])
     272    return True,np.sum(errFit(result[0],ibrav,Peaks[7],Peaks[4:7],Pwr)**2),A
     273           
    333274def FitHKLZ(ibrav,peaks,Z,A):
    334275    return A,Z
    … …  
    377318    dmin = getDmin(peaks)
    378319    smin = 1.0e10
    379     pwr = 6
     320    pwr = 3
    380321    maxTries = 10
    381322    if ibrav == 13:
    … …  
    386327    HKL = G2lat.GenHBravais(dmin,ibrav,A)
    387328    while IndexPeaks(peaks,HKL):
    388         Pwr = pwr - 2*(tries % 2)
     329        Pwr = pwr - (tries % 2)
    389330        HKL = []
    390331        tries += 1
    391332        osmin = smin
    392333        oldA = A
    393         OK,smin = FitHKL(ibrav,peaks,A,Pwr)
    394         for a in A[:3]:
    395             if a < 0:
    396                 A = oldA
    397                 OK = False
    398                 break
     334        OK,smin,A = FitHKL(ibrav,peaks,A,Pwr)
     335        if min(A[:3]) <= 0:
     336            A = oldA
     337            OK = False
     338            break
    399339        if OK:
    400340            HKL = G2lat.GenHBravais(dmin,ibrav,A)
    … …  
    404344        if tries > maxTries: break
    405345    if OK:
    406         OK,smin = FitHKL(ibrav,peaks,A,4)
     346        OK,smin,A = FitHKL(ibrav,peaks,A,2)
     347        Peaks = np.array(peaks).T
     348        H = Peaks[4:7]
     349        Peaks[8] = 1./np.sqrt(G2lat.calc_rDsq(H,A))
     350        peaks = Peaks.T
     351       
    407352    M20,X20 = calc_M20(peaks,HKL)
    408     return len(HKL),M20,X20
     353    return len(HKL),M20,X20,A
    409354       
    410355def findBestCell(dlg,ncMax,A,Ntries,ibrav,peaks,V1):
    … …  
    433378            Anew = ranAbyV(ibrav,amin,amax,V1)
    434379        HKL = G2lat.GenHBravais(dmin,ibrav,Anew)
    435         if len(HKL) > mHKL[ibrav] and IndexPeaks(peaks,HKL):
    436             Lhkl,M20,X20 = refinePeaks(peaks,ibrav,Anew)
     380       
     381        if IndexPeaks(peaks,HKL) and len(HKL) > mHKL[ibrav]:
     382            Lhkl,M20,X20,Anew = refinePeaks(peaks,ibrav,Anew)
    437383            Asave.append([calc_M20(peaks,HKL),Anew[:]])
    438384            if ibrav == 9:                          #C-centered orthorhombic
    439385                for i in range(2):
    440386                    Anew = rotOrthoA(Anew[:])
    441                     Lhkl,M20,X20 = refinePeaks(peaks,ibrav,Anew)
     387                    Lhkl,M20,X20,Anew = refinePeaks(peaks,ibrav,Anew)
    442388                    HKL = G2lat.GenHBravais(dmin,ibrav,Anew)
    443389                    IndexPeaks(peaks,HKL)
    … …  
    445391            elif ibrav == 11:                      #C-centered monoclinic
    446392                Anew = swapMonoA(Anew[:])
    447                 Lhkl,M20,X20 = refinePeaks(peaks,ibrav,Anew)
     393                Lhkl,M20,X20,Anew = refinePeaks(peaks,ibrav,Anew)
    448394                HKL = G2lat.GenHBravais(dmin,ibrav,Anew)
    449395                IndexPeaks(peaks,HKL)
    … …  
    461407    X = sortM20(Asave)
    462408    if X:
    463         Lhkl,M20,X20 = refinePeaks(peaks,ibrav,X[0][1])
     409        Lhkl,M20,X20,A = refinePeaks(peaks,ibrav,X[0][1])
    464410        return GoOn,Lhkl,M20,X20,X[0][1]
    465411    else:
    466         return GoOn,0,0,0,0
     412        return GoOn,0,0,0,Anew
    467413       
    468414def monoCellReduce(ibrav,A):
    … …  
    489435def DoIndexPeaks(peaks,inst,controls,bravais):
    490436   
    491     def peakDspace(peaks,A):
    492         for peak in peaks:
    493             if peak[3]:
    494                 dsq = calc_rDsq(peak[4:7],A)
    495                 if dsq > 0:
    496                     peak[8] = 1./math.sqrt(dsq)
    497         return
    498437    delt = 0.005                                     #lowest d-spacing cushion - can be fixed?
    499438    amin = 2.5
    … …  
    541480                            if ibrav > 2:
    542481                                if not N2:
    543                                     GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,0,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
     482                                    A = []
     483                                    GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,A,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
    544484                                if A:
    545485                                    GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,A[:],N1s[ibrav],ibrav,peaks,0)

  • trunk/GSASIIlattice.py

    r296 r303  
    1717tand = lambda x: np.tan(x*np.pi/180.)
    1818atand = lambda x: 180.*np.arctan(x)/np.pi
    19 atan2d = lambda y,x: 180.*np.atan2(y,x)/np.pi
     19atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
    2020cosd = lambda x: np.cos(x*np.pi/180.)
    2121acosd = lambda x: 180.*np.arccos(x)/np.pi
    … …  
    597597                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
    598598    return sortHKLd(HKL,True,True)
     599
     600#Spherical harmonics routines
     601def OdfChk(SGLaue,L,M):
     602    if not L%2 and abs(M) <= L:
     603        if SGLaue == '0':                      #cylindrical symmetry
     604            if M == 0: return True
     605        elif SGLaue == '-1':
     606            return True
     607        elif SGLaue == '2/m':
     608            if not abs(M)%2: return True
     609        elif SGLaue == 'mmm':
     610            if not abs(M)%2 and M >= 0: return True
     611        elif SGLaue == '4/m':
     612            if not abs(M)%4: return True
     613        elif SGLaue == '4/mmm':
     614            if not abs(M)%4 and M >= 0: return True
     615        elif SGLaue in ['3R','3']:
     616            if not abs(M)%3: return True
     617        elif SGLaue in ['3mR','3m1','31m']:
     618            if not abs(M)%3 and M >= 0: return True
     619        elif SGLaue == '6/m':
     620            if not abs(M)%6: return True
     621        elif SGLaue == '6/mmm':
     622            if not abs(M)%6 and M >= 0: return True
     623        elif SGLaue == 'm3':
     624            if M > 0:
     625                if L%12 == 2:
     626                    if M <= L/12: return True
     627                else:
     628                    if M <= L/12+1: return True
     629        elif SGLaue == 'm3m':
     630            if M > 0:
     631                if L%12 == 2:
     632                    if M <= L/12: return True
     633                else:
     634                    if M <= L/12+1: return True
     635    return False
     636       
     637def GenSHCoeff(SGLaue,SamSym,L):
     638    coeffNames = []
     639    for iord in [2*i+2 for i in range(L/2)]:
     640        for m in [i-iord for i in range(2*iord+1)]:
     641            if OdfChk(SamSym,iord,m):
     642                for n in [i-iord for i in range(2*iord+1)]:
     643                    if OdfChk(SGLaue,iord,n):
     644                        coeffNames.append('C(%d,%d,%d)'%(iord,m,n))
     645    return coeffNames
     646
     647def CrsAng(H,cell,SGData):
     648    a,b,c,al,be,ga = cell
     649    SQ3 = 1.732050807569
     650    H1 = np.array([1,0,0])
     651    H2 = np.array([0,1,0])
     652    H3 = np.array([0,0,1])
     653    H4 = np.array([1,1,1])
     654    G,g = cell2Gmat(cell)
     655    Laue = SGData['SGLaue']
     656    Naxis = SGData['SGUniq']
     657    DH = np.inner(H,np.inner(G,H))
     658    if Laue == '2/m':
     659        if Naxis == 'a':
     660            DR = np.inner(H1,np.inner(G,H1))
     661            DHR = np.inner(H,np.inner(G,H1))
     662        elif Naxis == 'b':
     663            DR = np.inner(H2,np.inner(G,H2))
     664            DHR = np.inner(H,np.inner(G,H2))
     665        else:
     666            DR = np.inner(H3,np.inner(G,H3))
     667            DHR = np.inner(H,np.inner(G,H3))
     668    elif Laue in ['R3','R3m']:
     669        DR = np.inner(H4,np.inner(G,H4))
     670        DRH = np.inner(H,np.inner(G,H4))
     671       
     672    else:
     673        DR = np.inner(H3,np.inner(G,H3))
     674        DHR = np.inner(H,np.inner(G,H3))
     675    DHR /= np.sqrt(DR*DH)
     676    phi = acosd(DHR)
     677    if Laue == '-1':
     678        BA = H[1]*a/(b-H[0]*cosd(ga))
     679        BB = H[0]*sind(ga)**2
     680    elif Laue == '2/m':
     681        if Naxis == 'a':
     682            BA = H[2]*b/(c-H[1]*cods(al))
     683            BB = H[1]*sind(al)**2
     684        elif Naxis == 'b':
     685            BA = H[0]*c/(a-H[2]*cosd(be))
     686            BB = H[2]*sind(be)**2
     687        else:
     688            BA = H[1]*a/(b-H[0]*cosd(ga))
     689            BB = H[0]*sind(ga)**2
     690    elif Laue in ['mmm','4/m','4/mmm']:
     691        BA = H[1]*a
     692        BB = H[0]*b
     693   
     694    elif Laue in ['3R','3mR']:
     695        BA = H[0]+H[1]-2.0*H[2]
     696        BB = SQ3*(H[0]-H[1])
     697    elif Laue in ['m3','m3m']:
     698        BA = H[1]
     699        BB = H[0]
     700    else:
     701        BA = H[0]+2.0*H[1]
     702        BB = SQ3*H[0]
     703    beta = atan2d(BA,BB)
     704    return phi,beta
     705   
     706def SamAng(Tth,Gangls,Sangl,IFCoup):
     707    if IFCoup:
     708        GSomeg = sind(Gangls[2]+Tth)
     709        GComeg = cosd(Gangls[2]+Tth)
     710    else:
     711        GSomeg = sind(Gangls[2])
     712        GComeg = cosd(Gangls[2])
     713    GSTth = sind(Tth)
     714    GCTth = cosd(Tth)     
     715    GSazm = sind(Gangls[3])
     716    GCazm = cosd(Gangls[3])
     717    GSchi = sind(Gangls[1])
     718    GCchi = cosd(Gangls[1])
     719    GSphi = sind(Gangls[0]+Sangl[2])
     720    GCphi = cosd(Gangls[0]+Sangl[2])
     721    SSomeg = sind(Sangl[0])
     722    SComeg = cosd(Sangl[0])
     723    SSchi = sind(Sangl[1])
     724    SCchi = cosd(Sangl[1])
     725    AT = -GSTth*GComeg+GCTth*GCazm*GSomeg
     726    BT = GSTth*GSomeg+GCTth*GCazm*GComeg
     727    CT = -GCTth*GSazm*GSchi
     728    DT = -GCTth*GSazm*GCchi
     729   
     730    BC1 = -AT*GSphi+(CT+BT*GCchi)*GCphi
     731    BC2 = DT-BT*GSchi
     732    BC3 = AT*GCphi+(CT+BT*GCchi)*GSphi
     733     
     734    BC = BC1*SComeg*SCchi+BC2*SComeg*SSchi-BC3*SSomeg     
     735    psi = acosd(BC)
     736     
     737    BA = -BC1*SSchi+BC2*SCchi
     738    BB = BC1*SSomeg*SCchi+BC2*SSomeg*SSchi+BC3*SComeg
     739    gam = atand2(BB,BA)
     740       
     741    return psi,gam
     742   
     743def Flnh(Start,SHCoef,phi,beta,SGData):
     744    import pytexture as ptx
     745    BOH = {
     746    'L=2':[[],[],[]],
     747    'L=4':[[0.30469720,0.36418281],[],[]],
     748    'L=6':[[-0.14104740,0.52775103],[],[]],
     749    'L=8':[[0.28646862,0.21545346,0.32826995],[],[]],
     750    'L=10':[[-0.16413497,0.33078546,0.39371345],[],[]],
     751    'L=12':[[0.26141975,0.27266871,0.03277460,0.32589402],
     752        [0.09298802,-0.23773812,0.49446631],[]],
     753    'L=14':[[-0.17557309,0.25821932,0.27709173,0.33645360],[],[]],
     754    'L=16':[[0.24370673,0.29873515,0.06447688,0.00377,0.32574495],
     755        [0.12039646,-0.25330128,0.23950998,0.40962508],[]],
     756    'L=18':[[-0.16914245,0.17017340,0.34598142,0.07433932,0.32696037],
     757        [-0.06901768,0.16006562,-0.24743528,0.47110273],[]],
     758    'L=20':[[0.23067026,0.31151832,0.09287682,0.01089683,0.00037564,0.32573563],
     759        [0.13615420,-0.25048007,0.12882081,0.28642879,0.34620433],[]],
     760    'L=22':[[-0.16109560,0.10244188,0.36285175,0.13377513,0.01314399,0.32585583],
     761        [-0.09620055,0.20244115,-0.22389483,0.17928946,0.42017231],[]],
     762    'L=24':[[0.22050742,0.31770654,0.11661736,0.02049853,0.00150861,0.00003426,0.32573505],
     763        [0.13651722,-0.21386648,0.00522051,0.33939435,0.10837396,0.32914497],
     764        [0.05378596,-0.11945819,0.16272298,-0.26449730,0.44923956]],
     765    'L=26':[[-0.15435003,0.05261630,0.35524646,0.18578869,0.03259103,0.00186197,0.32574594],
     766        [-0.11306511,0.22072681,-0.18706142,0.05439948,0.28122966,0.35634355],[]],
     767    'L=28':[[0.21225019,0.32031716,0.13604702,0.03132468,0.00362703,0.00018294,0.00000294,0.32573501],
     768        [0.13219496,-0.17206256,-0.08742608,0.32671661,0.17973107,0.02567515,0.32619598],
     769        [0.07989184,-0.16735346,0.18839770,-0.20705337,0.12926808,0.42715602]],
     770    'L=30':[[-0.14878368,0.01524973,0.33628434,0.22632587,0.05790047,0.00609812,0.00022898,0.32573594],
     771        [-0.11721726,0.20915005,-0.11723436,-0.07815329,0.31318947,0.13655742,0.33241385],
     772        [-0.04297703,0.09317876,-0.11831248,0.17355132,-0.28164031,0.42719361]],
     773    'L=32':[[0.20533892,0.32087437,0.15187897,0.04249238,0.00670516,0.00054977,0.00002018,0.00000024,0.32573501],
     774        [0.12775091,-0.13523423,-0.14935701,0.28227378,0.23670434,0.05661270,0.00469819,0.32578978],
     775        [0.09703829,-0.19373733,0.18610682,-0.14407046,0.00220535,0.26897090,0.36633402]],
     776    'L=34':[[-0.14409234,-0.01343681,0.31248977,0.25557722,0.08571889,0.01351208,0.00095792,0.00002550,0.32573508],
     777        [-0.11527834,0.18472133,-0.04403280,-0.16908618,0.27227021,0.21086614,0.04041752,0.32688152],
     778        [-0.06773139,0.14120811,-0.15835721,0.18357456,-0.19364673,0.08377174,0.43116318]]
     779    }
     780   
     781    FORPI = 12.5663706143592
     782    RSQPI = 0.5641895835478
     783    SQ2 = 1.414213562373
     784
     785    if Start:
     786        ptx.pyqlmninit()
     787        Start = False
     788    Fln = np.zeros(len(SHCoef))
     789    for i,term in enumerate(SHCoef):
     790        if abs(SHCoef[term]) > 1e-6:
     791            l,m,n = eval(term.strip('C'))
     792            lNorm = 4.*np.pi*(2.*l+1.)
     793            if SGData['SGLaue'] in ['m3','m3m']:
     794                L = [i*4 for i in range(l/4)]
     795                Kcl = 0.0
     796                for i in L:
     797                    im = i/4
     798                    pcrs = ptx.pyplmpsi(l,i,phi)
     799                    Kcl += BOH['L='+str(l)][n-1][l/2-1]*pcrs*cosd(i*beta)       
     800            else:                #all but cubic
     801                pcrs = ptx.pyplmpsi(l,n,phi)*RSQPI
     802                if not n:
     803                    pcrs /= SQ2
     804                if SGData['SGLaue'] in ['mmm','4/mmm','6/mmm']:
     805                    if n%6 == 3:
     806                        Kcl = pcrs*sind(n*beta)
     807                    else:
     808                        Kcl = pcrs*cosd(n*beta)
     809                elif SGData['SGLaue'] in ['3mR','3m1','31m']:
     810                    Kcl = pcrs*cosd(n*beta)
     811                else:
     812                    Kcl = pcrs*(cosd(n*beta)+sind(n*beta))
     813                                   
     814            Fln[i] = SHCoef[term]*Kcl*lNorm
     815    ODFln = dict(zip(SHCoef.keys(),list(zip(SHCoef.values(),Fln))))
     816    return ODFln
     817   
     818def polfcal(ODFln,SamSym,psi,gam):
     819    RSQPI = 0.5641895835478
     820    SQ2 = 1.414213562373
     821    PolVal = 1.0
     822    for term in ODFln:
     823        if ODFln[term][1] > 1.e-3:
     824            l,m,n = eval(term.strip('C'))
     825            psrs = ptx.pyplmpsi(l,m,psi)
     826            if SamSym in ['-1','2/m']:
     827                if m:
     828                    Ksl = RSQPI*psrs*(cosd(m*gam)+sind(m*gam))
     829                else:
     830                    ksl = RSQPI*psrs/SQ2
     831            else:
     832                if m:
     833                    Ksl = RSQPI*psrs*cosd(m*gam)
     834                else:
     835                    Ksl = RSQPI*psrs/SQ2
     836        PolVal += ODFln[term][1]*Ksl
     837    return PolVal
    599838   
    600839# output from uctbx computed on platform darwin on 2010-05-28

  • trunk/GSASIIphsGUI.py

    r264 r303  
    17341734            FindBonds()
    17351735            G2plt.PlotStructure(self,data)
    1736            
    1737        
    1738 
    17391736
    17401737        dataDisplay = wx.Panel(drawOptions)
    … …  
    18501847        else:
    18511848            self.dataFrame.DataMenu.Enable(G2gd.wxID_DATADELETE,False)           
    1852         generalData = data['General']
     1849        generalData = data['General']       
    18531850        SGData = generalData['SGData']
     1851        try:
     1852            textureData = generalData['SH Texture']
     1853        except KeyError:            #fix old files!
     1854            textureData = generalData['SH Texture'] = {'Order':0,'Model':'cylindrical',
     1855                'Sample omega':[False,0.0],'Sample chi':[False,0.0],'Sample phi':[False,0.0],
     1856                'SH Coeff':[False,{}],'SHShow':False,'PFhkl':[0,0,1]}
     1857        shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
     1858        SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
     1859       
     1860        shAngles = ['omega','chi','phi']
    18541861        keyList = UseList.keys()
    18551862        keyList.sort()
    18561863        Indx = {}
     1864       
     1865        def SetSHCoef():
     1866            cofNames = G2lat.GenSHCoeff(SGData['SGLaue'],SamSym[textureData['Model']],textureData['Order'])
     1867            newSHCoef = dict(zip(cofNames,np.zeros(len(cofNames))))
     1868            SHCoeff = textureData['SH Coeff'][1]
     1869            for cofName in SHCoeff:
     1870                if cofName in  cofNames:
     1871                    newSHCoef[cofName] = SHCoeff[cofName]
     1872            return newSHCoef
    18571873       
    18581874        def OnShowData(event):
    … …  
    18611877            UseList[hist]['Show'] = Obj.GetValue()
    18621878            UpdateDData()
     1879           
     1880        def OnShOrder(event):
     1881            textureData['Order'] = int(shOrder.GetValue())
     1882            textureData['SH Coeff'][1] = SetSHCoef()
     1883            UpdateDData()
     1884                       
     1885        def OnShModel(event):
     1886            textureData['Model'] = shModel.GetValue()
     1887            textureData['SH Coeff'][1] = SetSHCoef()
     1888            UpdateDData()
     1889           
     1890        def OnSHRefine(event):
     1891            textureData['SH Coeff'][0] = shRef.GetValue()
     1892           
     1893        def OnSHShow(event):
     1894            textureData['SHShow'] = shShow.GetValue()
     1895            UpdateDData()
     1896           
     1897        def OnAngRef(event):
     1898            Obj = event.GetEventObject()
     1899            textureData[angIndx[Obj.GetId()]][0] = Obj.GetValue()
     1900           
     1901        def OnAngValue(event):
     1902            Obj = event.GetEventObject()
     1903            try:
     1904                value =  float(Obj.GetValue())
     1905            except ValueError:
     1906                value = textureData[valIndx[Obj.GetId()]][1]
     1907            Obj.SetValue('%8.2f'%(value))
     1908            textureData[valIndx[Obj.GetId()]][1] = value
     1909           
     1910        def OnODFValue(event):
     1911            Obj = event.GetEventObject()
     1912            try:
     1913                value =  float(Obj.GetValue())
     1914            except ValueError:
     1915                value = textureData['SH Coeff'][1][ODFIndx[Obj.GetId()]]
     1916            Obj.SetValue('%8.3f'%(value))
     1917            textureData['SH Coeff'][1][ODFIndx[Obj.GetId()]] = value
     1918           
     1919        def OnPFValue(event):
     1920            Obj = event.GetEventObject()
     1921            try:
     1922                value =  int(Obj.GetValue())
     1923            except ValueError:
     1924                value = textureData['PFhkl'][pfIndx[Obj.GetId()]]
     1925            Obj.SetValue('%3d'%(value))
     1926            textureData['PFhkl'][pfIndx[Obj.GetId()]] = value
     1927            G2plt.PlotSphHarmTexture(self,generalData)           
    18631928       
    18641929        def OnScaleRef(event):
    … …  
    18761941            Obj.SetValue("%.4f"%(UseList[Indx[Obj.GetId()]]['Scale'][0]))          #reset in case of error
    18771942           
    1878         def OnCutoffVal(event):
    1879             Obj = event.GetEventObject()
    1880             try:
    1881                 cutoff = float(Obj.GetValue())
    1882                 if cutoff > 0:
    1883                     UseList[Indx[Obj.GetId()]]['Cutoff'] = cutoff
    1884             except ValueError:
    1885                 pass
    1886             Obj.SetValue("%.3f"%(UseList[Indx[Obj.GetId()]]['Cutoff']))          #reset in case of error
    1887 
    18881943        def OnSizeType(event):
    18891944            Obj = event.GetEventObject()
    … …  
    19291984            hist = Indx[Obj.GetId()]
    19301985            UseList[hist]['Mustrain'][0] = Obj.GetValue()
     1986            UpdateDData()
    19311987            G2plt.PlotStrain(self,data)
    1932             UpdateDData()
    19331988           
    19341989        def OnStrainRef(event):
    … …  
    20312086        dataDisplay = wx.Panel(DData)
    20322087        mainSizer = wx.BoxSizer(wx.VERTICAL)
     2088        mainSizer.Add(wx.StaticText(dataDisplay,-1,'Spherical harmonics texture data for '+PhaseName+':'),0,wx.ALIGN_CENTER_VERTICAL)
     2089        mainSizer.Add((0,5),0)
     2090        shSizer = wx.BoxSizer(wx.HORIZONTAL)
     2091        shSizer.Add(wx.StaticText(dataDisplay,-1,'Texture model: '),0,wx.ALIGN_CENTER_VERTICAL)
     2092        shModel = wx.ComboBox(dataDisplay,-1,value=textureData['Model'],choices=shModels,
     2093            style=wx.CB_READONLY|wx.CB_DROPDOWN)
     2094        shModel.Bind(wx.EVT_COMBOBOX,OnShModel)
     2095        shSizer.Add(shModel,0,wx.ALIGN_CENTER_VERTICAL)
     2096        shSizer.Add(wx.StaticText(dataDisplay,-1,'  Harmonic order: '),0,wx.ALIGN_CENTER_VERTICAL)
     2097        shOrder = wx.ComboBox(dataDisplay,-1,value=str(textureData['Order']),choices=[str(2*i) for i in range(18)],
     2098            style=wx.CB_READONLY|wx.CB_DROPDOWN)
     2099        shOrder.Bind(wx.EVT_COMBOBOX,OnShOrder)
     2100        shSizer.Add(shOrder,0,wx.ALIGN_CENTER_VERTICAL)
     2101        shSizer.Add((5,0),0)
     2102        shRef = wx.CheckBox(dataDisplay,label=' Refine texture?')
     2103        shRef.SetValue(textureData['SH Coeff'][0])
     2104        shRef.Bind(wx.EVT_CHECKBOX, OnSHRefine)
     2105        shSizer.Add(shRef,0,wx.ALIGN_CENTER_VERTICAL)
     2106        shShow = wx.CheckBox(dataDisplay,label=' Show coeff.?')
     2107        shShow.SetValue(textureData['SHShow'])
     2108        shShow.Bind(wx.EVT_CHECKBOX, OnSHShow)
     2109        shSizer.Add(shShow,0,wx.ALIGN_CENTER_VERTICAL)
     2110        mainSizer.Add(shSizer,0,0)
     2111        mainSizer.Add((0,5),0)
     2112        if textureData['SHShow']:
     2113            mainSizer.Add(wx.StaticText(dataDisplay,-1,'Spherical harmonic coefficients: '),0,wx.ALIGN_CENTER_VERTICAL)
     2114            mainSizer.Add((0,5),0)
     2115            ODFSizer = wx.FlexGridSizer(2,8,2,2)
     2116            ODFIndx = {}
     2117            ODFkeys = textureData['SH Coeff'][1].keys()
     2118            ODFkeys.sort()
     2119            for item in ODFkeys:
     2120                ODFSizer.Add(wx.StaticText(dataDisplay,-1,item),0,wx.ALIGN_CENTER_VERTICAL)
     2121                ODFval = wx.TextCtrl(dataDisplay,wx.ID_ANY,'%8.3f'%(textureData['SH Coeff'][1][item]),style=wx.TE_PROCESS_ENTER)
     2122                ODFIndx[ODFval.GetId()] = item
     2123                ODFval.Bind(wx.EVT_TEXT_ENTER,OnODFValue)
     2124                ODFval.Bind(wx.EVT_KILL_FOCUS,OnODFValue)
     2125                ODFSizer.Add(ODFval,0,wx.ALIGN_CENTER_VERTICAL)
     2126            mainSizer.Add(ODFSizer,0,wx.ALIGN_CENTER_VERTICAL)
     2127            mainSizer.Add((0,5),0)
     2128        PFSizer = wx.BoxSizer(wx.HORIZONTAL)
     2129        PFSizer.Add(wx.StaticText(dataDisplay,-1,'Display pole figure for HKL: '),0,wx.ALIGN_CENTER_VERTICAL)
     2130        pfIndx = {}
     2131        for i in range(3):
     2132            pfVal = wx.TextCtrl(dataDisplay,wx.ID_ANY,'%3d'%(textureData['PFhkl'][i]),size=(40,20),style=wx.TE_PROCESS_ENTER)
     2133            pfIndx[pfVal.GetId()] = i
     2134            pfVal.Bind(wx.EVT_TEXT_ENTER,OnPFValue)
     2135            pfVal.Bind(wx.EVT_KILL_FOCUS,OnPFValue)
     2136            PFSizer.Add(pfVal,0,wx.ALIGN_CENTER_VERTICAL)           
     2137        mainSizer.Add(PFSizer,0,wx.ALIGN_CENTER_VERTICAL)
     2138        mainSizer.Add((0,5),0)
     2139        mainSizer.Add(wx.StaticText(dataDisplay,-1,'Sample orientation angles: '),0,wx.ALIGN_CENTER_VERTICAL)
     2140        mainSizer.Add((0,5),0)
     2141        angSizer = wx.BoxSizer(wx.HORIZONTAL)
     2142        angIndx = {}
     2143        valIndx = {}
     2144        for item in ['Sample omega','Sample chi','Sample phi']:
     2145            angRef = wx.CheckBox(dataDisplay,label=item+': ')
     2146            angRef.SetValue(textureData[item][0])
     2147            angIndx[angRef.GetId()] = item
     2148            angRef.Bind(wx.EVT_CHECKBOX, OnAngRef)
     2149            angSizer.Add(angRef,0,wx.ALIGN_CENTER_VERTICAL)
     2150            angVal = wx.TextCtrl(dataDisplay,wx.ID_ANY,'%8.2f'%(textureData[item][1]),style=wx.TE_PROCESS_ENTER)
     2151            valIndx[angVal.GetId()] = item
     2152            angVal.Bind(wx.EVT_TEXT_ENTER,OnAngValue)
     2153            angVal.Bind(wx.EVT_KILL_FOCUS,OnAngValue)
     2154            angSizer.Add(angVal,0,wx.ALIGN_CENTER_VERTICAL)
     2155            angSizer.Add((5,0),0)
     2156        mainSizer.Add(angSizer,0,wx.ALIGN_CENTER_VERTICAL)
    20332157        mainSizer.Add(wx.StaticText(dataDisplay,-1,'Histogram data for '+PhaseName+':'),0,wx.ALIGN_CENTER_VERTICAL)
    20342158        for item in keyList:
    20352159            histData = UseList[item]
    2036             mainSizer.Add(wx.StaticText(dataDisplay,-1,50*'_'))               
    20372160            mainSizer.Add((5,5),0)
    20382161            showData = wx.CheckBox(dataDisplay,label=' Show '+item)
    … …  
    20592182               
    20602183            if item[:4] == 'PWDR' and UseList[item]['Show']:
    2061                 cutoffSizer = wx.BoxSizer(wx.HORIZONTAL)
    2062                 cutoffSizer.Add(wx.StaticText(dataDisplay,label=' Peak cutoff ratio: '),0,wx.ALIGN_CENTER_VERTICAL)
    2063                 cutoffVal = wx.TextCtrl(dataDisplay,wx.ID_ANY,'%.3f'%(UseList[item]['Cutoff']),
    2064                     style=wx.TE_PROCESS_ENTER)               
    2065                 Indx[cutoffVal.GetId()] = item
    2066                 cutoffVal.Bind(wx.EVT_TEXT_ENTER,OnCutoffVal)
    2067                 cutoffVal.Bind(wx.EVT_KILL_FOCUS,OnCutoffVal)
    2068                 cutoffSizer.Add(cutoffVal,0,wx.ALIGN_CENTER_VERTICAL)
    2069                 mainSizer.Add(cutoffSizer)
    20702184                mainSizer.Add((0,5),0)
    20712185                sizeSizer = wx.BoxSizer(wx.HORIZONTAL)
    … …  
    22342348
    22352349        dataDisplay.SetSizer(mainSizer)
     2350        mainSizer.FitInside(self.dataFrame)
    22362351        Size = mainSizer.Fit(self.dataFrame)
    22372352        Size[0] = max(Size[0],300)+20
    2238         Size[1] += 30                           #compensate for status bar
     2353        Size[1] += 30                        #compensate for status bar
    22392354        DData.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-10)
    22402355        dataDisplay.SetSize(Size)
    … …  
    22902405                            'Mustrain':['isotropic',[1.0,0.0],[False,False],[0,0,1],
    22912406                                NShkl*[0.01,],NShkl*[False,]],                           
    2292                             'Extinction':[0.0,False],'Cutoff':0.01}
     2407                            'Extinction':[0.0,False]}
    22932408                    data['Histograms'] = UseList
    22942409                    UpdateDData()

  • trunk/GSASIIplot.py

    r296 r303  
    2828import GSASIIlattice as G2lat
    2929import GSASIIspc as G2spc
     30import pytexture as ptx
    3031from  OpenGL.GL import *
    3132from OpenGL.GLU import *
    … …  
    154155       
    155156def PlotSngl(self,newPlot=False):
     157    '''Single crystal structure factor plotting package - displays zone of reflections as rings proportional
     158        to F, F**2, etc. as requested
     159    '''
    156160    from matplotlib.patches import Circle
    157161    global HKL,HKLF
    … …  
    279283       
    280284def PlotPatterns(self,newPlot=False):
     285    '''Powder pattern plotting package - displays single or multiple powder patterns as intensity vs
     286    2-theta or q (future TOF). Can display multiple patterns as "waterfall plots" or contour plots. Log I
     287    plotting available.
     288    '''
    281289    global HKL
    282290   
    … …  
    302310        if self.PatternTree.GetItemText(PickId) == 'Peak List':
    303311            if ind.all() != [0]:                                    #picked a data point
    304                 ins = [Parms[x] for x in ['U','V','W','X','Y','SH/L']]
    305                 if self.qPlot:                              #qplot - convert back to 2-theta
    306                     xy[0] = 2.0*asind(xy[0]*wave/(4*math.pi))
    307                 sig = ins[0]*tand(xy[0]/2.0)**2+ins[1]*tand(xy[0]/2.0)+ins[2]
    308                 gam = ins[3]/cosd(xy[0]/2.0)+ins[4]*tand(xy[0]/2.0)           
    309                 data = self.PatternTree.GetItemPyData(self.PickId)
    310                 XY = [xy[0],0, xy[1],1, sig,0, gam,0]       #default refine intensity 1st
     312                if 'C' in Parms['Type']:                            #CW data - TOF later in an elif
     313                    ins = [Parms[x] for x in ['U','V','W','X','Y']]
     314                    if self.qPlot:                              #qplot - convert back to 2-theta
     315                        xy[0] = 2.0*asind(xy[0]*wave/(4*math.pi))
     316                    sig = ins[0]*tand(xy[0]/2.0)**2+ins[1]*tand(xy[0]/2.0)+ins[2]
     317                    gam = ins[3]/cosd(xy[0]/2.0)+ins[4]*tand(xy[0]/2.0)           
     318                    data = self.PatternTree.GetItemPyData(self.PickId)
     319                    XY = [xy[0],0, xy[1],1, sig,0, gam,0]       #default refine intensity 1st
    311320                data.append(XY)
    312321                G2pdG.UpdatePeakGrid(self,data)
    … …  
    318327                LimitId = G2gd.GetPatternTreeItemId(self,PatternId, 'Limits')
    319328                data = self.PatternTree.GetItemPyData(LimitId)
    320                 if self.qPlot:                              #qplot - convert back to 2-theta
    321                     xy[0] = 2.0*asind(xy[0]*wave/(4*math.pi))
     329                if 'C' in Parms['Type']:                            #CW data - TOF later in an elif
     330                    if self.qPlot:                              #qplot - convert back to 2-theta
     331                        xy[0] = 2.0*asind(xy[0]*wave/(4*math.pi))
    322332                if mouse.button==1:
    323333                    data[1][0] = min(xy[0],data[1][1])
    … …  
    528538        else:
    529539            Choice = (' key press','l: offset left','r: offset right','d: offset down',
    530                 'u: offset up','0: reset offset','n: log(I) on','c: contour on',
     540                'u: offset up','o: reset offset','n: log(I) on','c: contour on',
    531541                'q: toggle q plot','s: toggle single plot','+: no selection')
    532542    cb = wx.ComboBox(self.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,
    … …  
    683693   
    684694def PlotISFG(self,newPlot=False,type=''):
     695    ''' PLotting package for PDF analysis; displays I(q), S(q), F(q) and G(r) as single
     696    or multiple plots with waterfall and contour plots as options
     697    '''
    685698    if not type:
    686699        type = self.G2plotNB.plotList[self.G2plotNB.nb.GetSelection()]
    … …  
    894907       
    895908def PlotXY(self,XY,newPlot=False,type=''):
    896     #simple plot of xy data
    897 
     909    '''simple plot of xy data, used for diagnostic purposes
     910    '''
    898911    def OnMotion(event):
    899912        xpos = event.xdata
    … …  
    943956
    944957def PlotPowderLines(self):
     958    ''' plotting of powder lines (i.e. no powder pattern) as sticks
     959    '''
    945960    global HKL
    946961
    … …  
    9931008
    9941009def PlotPeakWidths(self):
     1010    ''' Plotting of instrument broadening terms as function of 2-theta (future TOF)
     1011    Seen when "Instrument Parameters" chosen from powder pattern data tree
     1012    '''
    9951013    PatternId = self.PatternId
    9961014    limitID = G2gd.GetPatternTreeItemId(self,PatternId, 'Limits')
    … …  
    10841102   
    10851103def PlotStrain(self,data):
    1086     # in this instance data is for a phase
     1104    '''Plot 3D microstrain figure. In this instance data is for a phase
     1105    '''
    10871106    PatternId = self.PatternId
    10881107    generalData = data['General']
    … …  
    11611180            Plot.set_zlabel('Z')
    11621181    Page.canvas.draw()
     1182   
     1183def PlotSphHarmTexture(self,generalData):
     1184    '''Pole figure, inverse pole figure(?), 3D pole distribution and 3D inverse pole distribution(?)
     1185    plotting; Need way to select 
     1186    pole figure or pole distribution to be displayed - do in key enter menu
     1187    dict generalData contains all phase info needed
     1188    '''
     1189    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
     1190    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
     1191    SGData = generalData['SGData']
     1192    textureData = generalData['SH Texture']
     1193    print 'Texture plot'
     1194    SHData = generalData['SH Texture']
     1195    SHCoef = SHData['SH Coeff'][1]
     1196    cell = generalData['Cell'][1:7]
     1197    PH = np.array(SHData['PFhkl'])
     1198    phi,beta = G2lat.CrsAng(PH,cell,SGData)
     1199    Start = True
     1200    ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
     1201
    11631202           
    11641203def PlotExposedImage(self,newPlot=False,event=None):
     1204    '''General access module for 2D image plotting
     1205    '''
    11651206    plotNo = self.G2plotNB.nb.GetSelection()
    11661207    if self.G2plotNB.nb.GetPageText(plotNo) == '2D Powder Image':
    … …  
    11701211
    11711212def PlotImage(self,newPlot=False,event=None,newImage=True):
     1213    '''Plot of 2D detector images as contoured plot. Also plot calibration ellipses,
     1214    masks, etc.
     1215    '''
    11721216    from matplotlib.patches import Ellipse,Arc,Circle,Polygon
    11731217    import numpy.ma as ma
    … …  
    15961640       
    15971641def PlotIntegration(self,newPlot=False,event=None):
     1642    '''Plot of 2D image after image integration with 2-theta and azimuth as coordinates
     1643    '''
    15981644           
    15991645    def OnMotion(event):
    … …  
    16611707               
    16621708def PlotTRImage(self,tax,tay,taz,newPlot=False):
    1663     #a test plot routine - not normally used
     1709    '''a test plot routine - not normally used
     1710    '''
    16641711           
    16651712    def OnMotion(event):
    … …  
    17331780       
    17341781def PlotStructure(self,data):
     1782    '''Crystal structure plotting package. Can show structures as balls, sticks, lines,
     1783    thermal motion ellipsoids and polyhedra
     1784    '''
    17351785    generalData = data['General']
    17361786    cell = generalData['Cell'][1:7]

  • trunk/GSASIIpwd.py

    r301 r303  
    2323
    2424import GSASIIpath
    25 #import pypowder as pyp
    2625import GSASIIplot as G2plt
    2726import GSASIIlattice as G2lat

  • trunk/GSASIIpwdGUI.py

    r300 r303  
    790790        ibrav = bravaisSymb.index(controls[5])
    791791        dmin = G2indx.getDmin(peaks)-0.005
    792         Lhkl,M20,X20 = G2indx.refinePeaks(peaks,ibrav,A)
     792        Lhkl,M20,X20,A = G2indx.refinePeaks(peaks,ibrav,A)
    793793        controls[6:12] = G2lat.A2cell(A)
    794794        controls[12] = G2lat.calc_V(A)
    … …  
    825825            return
    826826        self.dataFrame.CopyCell.Enable(False)
     827        self.dataFrame.RefineCell.Enable(False)
    827828        OK,dmin,cells = G2indx.DoIndexPeaks(peaks,inst,controls,bravais)
    828829        if OK:
    … …  
    839840                else:
    840841                    G2plt.PlotPatterns(self)
    841         self.dataFrame.CopyCell.Enable(True)
    842         self.dataFrame.IndexPeaks.Enable(True)
    843         self.dataFrame.MakeNewPhase.Enable(True)
    844         UpdateUnitCellsGrid(self,data)
     842            self.dataFrame.CopyCell.Enable(True)
     843            self.dataFrame.IndexPeaks.Enable(True)
     844            self.dataFrame.MakeNewPhase.Enable(True)
     845            UpdateUnitCellsGrid(self,data)
    845846               
    846847    def CopyUnitCell(event):
Note: See TracChangeset for help on using the changeset viewer.