Changeset 303


Ignore:
Timestamp:
Jun 20, 2011 2:37:07 PM (10 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.