Changeset 798


Ignore:
Timestamp:
Nov 14, 2012 4:37:26 PM (11 years ago)
Author:
vondreele
Message:

more peak fitting work including TOF

Location:
trunk
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASII.py

    r797 r798  
    658658                return [G2IO.makeInstDict(names,data,codes),{}]
    659659            elif 'T' in DataType:
    660                 names = ['Type','2-theta','difC','difA','Zero','alpha','beta-0','beta-1','var-inst','X','Y','Azimuth']
    661                 codes = [0,0,0,0,0,0,0,0,0,0,0,0]
     660                names = ['Type','2-theta','difC','difA','Zero','alpha','beta-0','beta-1','sig-0','sig-1','X','Y','Azimuth']
     661                codes = [0,0,0,0,0,0,0,0,0,0,0,0,0]
    662662                azm = 0.
    663663                if 'INS  1DETAZM' in Iparm:
    … …  
    673673                    data.extend([G2IO.sfloat(s[1]),G2IO.sfloat(s[2]),G2IO.sfloat(s[3])])
    674674                    s = Iparm['INS  1PRCF12'].split()
    675                     data.extend([G2IO.sfloat(s[1]),0.0,0.0,azm])
     675                    data.extend([0.0,G2IO.sfloat(s[1]),0.0,0.0,azm])
    676676                elif abs(pfType) in [3,4,5]:
    677677                    data.extend([G2IO.sfloat(s[0]),G2IO.sfloat(s[1]),G2IO.sfloat(s[2])])
    678678                    if abs(pfType) == 4:
    679                         data.extend([G2IO.sfloat(s[3]),0.0,0.0,azm])
     679                        data.extend([0.0,G2IO.sfloat(s[3]),0.0,0.0,azm])
    680680                    else:
    681681                        s = Iparm['INS  1PRCF12'].split()
    682                         data.extend([G2IO.sfloat(s[0]),0.0,0.0,azm])                       
     682                        data.extend([0.0,G2IO.sfloat(s[0]),0.0,0.0,azm])                       
    683683                Inst1 = G2IO.makeInstDict(names,data,codes)
    684684                Inst2 = {}
    … …  
    692692                        Inst2['Pdabc'].append([float(t) for t in s])
    693693                    Inst2['Pdabc'] = np.array(Inst2['Pdabc'])
    694                     Inst2['Pdabc'][3] += Inst2['Pdabc'][0]*Inst1['difC'][0] #turn 3rd col into TOF
     694                    Inst2['Pdabc'].T[3] += Inst2['Pdabc'].T[0]*Inst1['difC'][0] #turn 3rd col into TOF
    695695                if 'INS  1I ITYP' in Iparm:
    696696                    s = Iparm['INS  1I ITYP'].split()
    … …  
    866866            if 'T' in Iparm1['Type'][0]:
    867867                if not rd.clockWd and rd.GSAS:
    868                     rd.powderdata[0] *= 100.
     868                    rd.powderdata[0] *= 100.        #put back the CW centideg correction
    869869                cw = np.diff(rd.powderdata[0])
    870870                rd.powderdata[0] = rd.powderdata[0][:-1]+cw/2.
    871871                rd.powderdata[1] = rd.powderdata[1][:-1]/cw
    872                 rd.powderdata[2] = rd.powderdata[2][:-1]/cw**2
     872                rd.powderdata[2] = rd.powderdata[2][:-1]*cw**2  #1/var=w at this point
    873873                if 'Itype' in Iparm2:
    874874                    Ibeg = np.searchsorted(rd.powderdata[0],Iparm2['Tminmax'][0])
    … …  
    878878                    rd.powderdata[1] = rd.powderdata[1][Ibeg:Ifin]/YI
    879879                    var = 1./rd.powderdata[2][Ibeg:Ifin]
    880                     var += rd.powderdata[1]**2+WYI
     880                    var += WYI*rd.powderdata[1]**2
    881881                    var /= YI**2
    882882                    rd.powderdata[2] = 1./var

  • trunk/GSASIIIO.py

    r797 r798  
    13341334    defaultIparm_lbl.append('10m TOF backscattering bank')
    13351335    defaultIparms.append({
     1336        'INS   HTYPE ':'PNT',
    13361337        'INS  1 ICONS':'   5000.00      0.00      0.00',
    13371338        'INS  1BNKPAR':'    1.0000   150.000',       
    … …  
    13421343    defaultIparm_lbl.append('10m TOF 90deg bank')
    13431344    defaultIparms.append({
     1345        'INS   HTYPE ':'PNT',
    13441346        'INS  1 ICONS':'   3500.00      0.00      0.00',
    13451347        'INS  1BNKPAR':'    1.0000    90.000',       
    13461348        'INS  1PRCF1 ':'    1    8   0.01000',
    13471349        'INS  1PRCF11':'   0.000000E+00   5.000000E+00   3.000000E-02   4.000000E-03',
     1350        'INS  1PRCF12':'   0.000000E+00   8.000000E+01   0.000000E+00   0.000000E+00',       
     1351        })
     1352    defaultIparm_lbl.append('63m POWGEN 90deg bank')
     1353    defaultIparms.append({
     1354        'INS   HTYPE ':'PNT',
     1355        'INS  1 ICONS':'  22585.80      0.00      0.00',
     1356        'INS  1BNKPAR':'    1.0000    90.000',       
     1357        'INS  1PRCF1 ':'    1    8   0.01000',
     1358        'INS  1PRCF11':'   0.000000E+00   1.000000E+00   3.000000E-02   4.000000E-03',
    13481359        'INS  1PRCF12':'   0.000000E+00   8.000000E+01   0.000000E+00   0.000000E+00',       
    13491360        })

  • trunk/GSASIImath.py

    r797 r798  
    985985    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
    986986        for x in ['U','V','W','X','Y']:
    987             ins[x] = Parms[x][1]
     987            ins[x] = Parms[x][0]
    988988        if ifQ:                              #qplot - convert back to 2-theta
    989989            pos = 2.0*asind(pos*wave/(4*math.pi))
    … …  
    992992        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
    993993    else:
    994         dsp = pos/Parms['difC'][1]
     994        if ifQ:
     995            dsp = 2.*np.pi/pos
     996            pos = Parms['difC']*dsp
     997        else:
     998            dsp = pos/Parms['difC'][1]
    995999        if 'Pdabc' in Parms2:
    9961000            for x in ['var-inst','X','Y']:
    997                 ins[x] = Parms[x][1]
     1001                ins[x] = Parms[x][0]
    9981002            Pdabc = Parms2['Pdabc'].T
    9991003            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
    10001004            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
    10011005        else:
    1002             for x in ['alpha','beta-0','beta-0','var-inst','X','Y']:
    1003                 ins[x] = Parms[x][1]
     1006            for x in ['alpha','beta-0','beta-1','sig-0','sig-1','X','Y']:
     1007                ins[x] = Parms[x][0]
    10041008            alp = ins['alpha']/dsp
    1005             bet = ins['beta-0']+ins['beta-0']/dsp**4
    1006         sig = ins['var-inst']*dsp**2
     1009            bet = ins['beta-0']+ins['beta-1']/dsp**4
     1010        sig = ins['sig-0']+ins['sig-1']*dsp**2
    10071011        gam = ins['X']*dsp+ins['Y']*dsp**2
    10081012        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]

  • trunk/GSASIIplot.py

    r797 r798  
    5252npsind = lambda x: np.sin(x*np.pi/180.)
    5353npcosd = lambda x: np.cos(x*np.pi/180.)
     54nptand = lambda x: np.tan(x*np.pi/180.)
    5455npacosd = lambda x: 180.*np.arccos(x)/np.pi
    5556npasind = lambda x: 180.*np.arcsin(x)/np.pi
    … …  
    12201221           
    12211222def PlotPeakWidths(G2frame):
    1222     ''' Plotting of instrument broadening terms as function of 2-theta (future TOF)
     1223    ''' Plotting of instrument broadening terms as function of 2-theta
    12231224    Seen when "Instrument Parameters" chosen from powder pattern data tree
    12241225    '''
    12251226    sig = lambda Th,U,V,W: 1.17741*math.sqrt(U*tand(Th)**2+V*tand(Th)+W)*math.pi/18000.
    12261227    gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))*math.pi/18000.
    1227     gamFW = lambda s,g: math.exp(math.log(s**5+2.69269*s**4*g+2.42843*s**3*g**2+4.47163*s**2*g**3+0.07842*s*g**4+g**5)/5.)
     1228    gamFW = lambda s,g: np.exp(np.log(s**5+2.69269*s**4*g+2.42843*s**3*g**2+4.47163*s**2*g**3+0.07842*s*g**4+g**5)/5.)
    12281229#    gamFW2 = lambda s,g: math.sqrt(s**2+(0.4654996*g)**2)+.5345004*g  #Ubaldo Bafile - private communication
    12291230    PatternId = G2frame.PatternId
    … …  
    12361237        G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Instrument Parameters'))
    12371238    if 'C' in Parms['Type'][0]:
    1238         try:           
    1239             lam = Parms['Lam'][1]
    1240         except KeyError:
    1241             lam = Parms['Lam1'][1]
    1242         GU = Parms['U'][0]
    1243         GV = Parms['V'][0]
    1244         GW = Parms['W'][0]
    1245         LX = Parms['X'][0]
    1246         LY = Parms['Y'][0]
     1239        lam = G2mth.getWave(Parms)
    12471240    else:
    12481241        difC = Parms['difC'][0]
    1249         difA = Parms['difA'][0]
    1250         Zero = Parms['Zero'][0]
    1251         alp = Parms['alpha'][0]
    1252         bet0 = Parms['beta-0'][0]
    1253         bet1 = Parms['beta-1'][0]
    1254         sig = Parms['var-inst'][0]
    1255         LX = Parms['X'][0]
    1256         LY = Parms['Y'][0]
    12571242    peakID = G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Peak List')
    12581243    if peakID:
    … …  
    12861271            Xmin,Xmax = limits[1]
    12871272            Xmin = min(0.5,max(Xmin,1))
    1288             Xmin /= 2
    1289             Xmax /= 2
    1290             nPts = 100
    1291             delt = (Xmax-Xmin)/nPts
    1292             thetas = []
    1293             for i in range(nPts):
    1294                 thetas.append(Xmin+i*delt)
    1295             for theta in thetas:
    1296                 X.append(4.0*math.pi*sind(theta)/lam)              #q
    1297                 s = sig(theta,GU,GV,GW)
    1298                 g = gam(theta,LX,LY)
    1299                 G = gamFW(g,s)
    1300                 Y.append(s/tand(theta))
    1301                 Z.append(g/tand(theta))
    1302                 W.append(G/tand(theta))
    1303             Plot.plot(X,Y,color='r',label='Gaussian')
    1304             Plot.plot(X,Z,color='g',label='Lorentzian')
    1305             Plot.plot(X,W,color='b',label='G+L')
     1273            X = np.linspace(Xmin,Xmax,num=101,endpoint=True)
     1274            Q = 4.*np.pi*npsind(X/2.)/lam
     1275            Z = np.ones_like(X)
     1276            data = G2mth.setPeakparms(Parms,Parms2,X,Z)
     1277            s = 1.17741*np.sqrt(data[4])*np.pi/18000.
     1278            g = data[6]*np.pi/18000.
     1279            G = gamFW(g,s)
     1280            Y = s/nptand(X/2.)
     1281            Z = g/nptand(X/2.)
     1282            W = G/nptand(X/2.)
     1283            Plot.plot(Q,Y,color='r',label='Gaussian')
     1284            Plot.plot(Q,Z,color='g',label='Lorentzian')
     1285            Plot.plot(Q,W,color='b',label='G+L')
     1286           
    13061287            X = []
    13071288            Y = []
    … …  
    13311312    else:
    13321313        Plot.set_title('Instrument and sample peak coefficients')
    1333         Plot.set_xlabel(r'$TOF, \mu s$',fontsize=14)
    1334         Plot.set_ylabel(r'$\alpha, \beta, \Delta T$',fontsize=14)
     1314        Plot.set_xlabel(r'$q, \AA^{-1}$',fontsize=14)
     1315        Plot.set_ylabel(r'$\alpha, \beta, \Delta q/q, \Delta d/d$',fontsize=14)
    13351316        Xmin,Xmax = limits[1]
    1336         if 'Pabc' in Parms2:
    1337             Pabc = Parms2['Pabc']
    1338             print Pabc.shape,Pabc[0]
    1339            
    1340         else:           
    1341             T = np.linspace(Xmin,Xmax,num=101,endpoint=True)
    1342             ds = T/difC
    1343             A = alp/ds
    1344             B = bet0+bet1/ds**4
    1345             D = difA*ds**2+Zero
    1346             Plot.plot(T,A,color='r',label='Alpha')
    1347             Plot.plot(T,B,color='g',label='Beta')
    1348             Plot.plot(T,D,color='b',label='Delta-T')
     1317        T = np.linspace(Xmin,Xmax,num=101,endpoint=True)
     1318        Z = np.ones_like(T)
     1319        data = G2mth.setPeakparms(Parms,Parms2,T,Z)
     1320#       = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
     1321        ds = T/difC
     1322        Q = 2.*np.pi/ds
     1323        A = data[4]
     1324        B = data[6]
     1325        S = 1.17741*np.sqrt(data[8])/T
     1326        G = data[10]/T
     1327        Plot.plot(Q,A,color='r',label='Alpha')
     1328        Plot.plot(Q,B,color='g',label='Beta')
     1329        Plot.plot(Q,S,color='b',label='Gaussian')
     1330        Plot.plot(Q,G,color='m',label='Lorentzian')
    13491331        T = []
    13501332        A = []
    13511333        B = []
     1334        S = []
     1335        G = []
    13521336        W = []
     1337        Q = []
    13531338        V = []
    13541339        for peak in peaks:
    … …  
    13561341            A.append(peak[4])
    13571342            B.append(peak[6])
    1358             W.append(peak[0]/difC)
    1359            
    1360        
    1361         Plot.plot(T,A,'+',color='r',label='Alpha peak')
    1362         Plot.plot(T,B,'+',color='g',label='Beta peak')
    1363         Plot.plot(T,W,'+',color='k',label='T/difC')
     1343            Q.append(2.*np.pi*difC/peak[0])
     1344            S.append(1.17741*np.sqrt(peak[8])/peak[0])
     1345            G.append(peak[10]/peak[0])
     1346           
     1347       
     1348        Plot.plot(Q,A,'+',color='r',label='Alpha peak')
     1349        Plot.plot(Q,B,'+',color='g',label='Beta peak')
     1350        Plot.plot(Q,S,'+',color='b',label='Gaussian peak')
     1351        Plot.plot(Q,G,'+',color='m',label='Lorentzian peak')
    13641352        Plot.legend(loc='best')
    13651353        Page.canvas.draw()

  • trunk/GSASIIpwd.py

    r797 r798  
    196196    return sumNoAtoms/Vol
    197197           
    198 #def MultGetQ(Tth,MuT,Geometry,b=88.0,a=0.01):
    199 #    NS = 500
    200 #    Gama = np.linspace(0.,np.pi/2.,NS,False)[1:]
    201 #    Cgama = np.cos(Gama)[:,np.newaxis]
    202 #    Sgama = np.sin(Gama)[:,np.newaxis]
    203 #    CSgama = 1.0/Sgama
    204 #    Delt = Gama[1]-Gama[0]
    205 #    emc = 7.94e-26
    206 #    Navo = 6.023e23
    207 #    Cth = npcosd(Tth/2.0)
    208 #    CTth = npcosd(Tth)
    209 #    Sth = npcosd(Tth/2.0)
    210 #    STth = npsind(Tth)
    211 #    CSth = 1.0/Sth
    212 #    CSTth = 1.0/STth
    213 #    SCth = 1.0/Cth
    214 #    SCTth = 1.0/CTth
    215 #    if 'Bragg' in Geometry:
    216 #        G = 8.0*Delt*Navo*emc*Sth/((1.0-CTth**2)*(1.0-np.exp(-2.0*MuT*CSth)))
    217 #        Y1 = np.pi
    218 #        Y2 = np.pi/2.0
    219 #        Y3 = 3.*np.pi/8. #3pi/4?
    220 #        W = 1.0/(Sth+np.fabs(Sgama))
    221 #        W += np.exp(-MuT*CSth)*(2.0*np.fabs(Sgama)*np.exp(-MuT*np.fabs(CSgama))-
    222 #            (Sth+np.fabs(Sgama))*np.exp(-MuT*CSth))/(Sth**2-Sgama**2)
    223 #        Fac0 = Sth**2*Sgama**2
    224 #        X = Fac0+(Fac0+CTth)**2/2
    225 #        Y = Cgama**2*Cth**2*(1.0-Fac0-CTth)
    226 #        Z = Cgama**4*Cth**4/2.0
    227 #        E = 2.0*(1.0-a)/(b*Cgama/Cth)
    228 #        F1 = (2.0+b*(1.0+Sth*Sgama))/(b*Cth*Cgama) #trouble if < 1
    229 #        F2 = (2.0+b*(1.0-Sth*Sgama))/(b*Cth*Cgama)
    230 #        T1 = np.pi/np.sqrt(F1**2-1.0)
    231 #        T2 = np.pi/np.sqrt(F2**2-1.0)
    232 #        Y4 = T1+T2
    233 #        Y5 = F1**2*T1+F2**2*T2-np.pi*(F1+F2)
    234 #        Y6 = F1**4*T1+F2**4*T2-np.pi*(F1+F2)/2.0-np.pi*(F1**3+F2**3)
    235 #        Y7 = (T2-T1)/(F1-F2)
    236 #        YT = F2**2*T2-F1**2*T1
    237 #        Y8 = Y1+YT/(F1-F2)
    238 #        Y9 = Y2+(F2**4*T2-F1**4*T1)/(F1-F2)+Y1*((F1+F2)**2-F1*F2)
    239 #        M = (a**2*(X*Y1+Y*Y2+Z*Y3)+a*E*(X*Y4+Y*Y5+Z*Y6)+E**2*(X*Y7+Y*Y8+Z*Y9))*Cgama
    240 #       
    241 #        Q = np.sum(W*M,axis=0)
    242 #        return Q*G*NS/(NS-1.)
    243 ##
    244 ##      cos2th=2.0d*costh^2 - 1.0d
    245 ##      G= delta * 8.0d * Na * emc * sinth/(1.0d + cos2th^2)/(1.0d - exp(-2.0d*mut*cscth))
    246 ##      Y1=3.1415926d
    247 ##      Y2=Y1*0.5d
    248 ##      Y3=Y2*0.75d
    249 ##      for i=1,num_steps-1 do begin
    250 ##         cosgama=double(cos(gama[i]))
    251 ##         singama=double(sin(gama[i]))
    252 ##         cscgama=1.0d / singama
    253 ##
    254 ##         W=1.0d/(sinth+abs(singama))
    255 ##         W=W+exp(-1.0*mut*cscth)*(2.0d*abs(singama)*exp(-1.0d*mut*abs(cscgama))- $
    256 ##                 (sinth+abs(singama))*exp(-1.0d*mut*cscth))/(sinth^2-singama^2)
    257 ##
    258 ##         factor0=sinth^2*singama^2
    259 ##         X=factor0+(factor0+cos2th)^2/2.0d
    260 ##         Y=cosgama^2*(1.0d - factor0-cos2th)*costh^2
    261 ##         Z=cosgama^4/2.0d*costh^4
    262 ##         E=2.0d*(1.0-a)/b/cosgama/costh
    263 ##
    264 ##         F1=1.0d/b/cosgama*(2.0d + b*(1.0+sinth*singama))/costh
    265 ##         F2=1.0d/b/cosgama*(2.0d + b*(1.0-sinth*singama))/costh
    266 ##
    267 ##         T1=3.14159/sqrt(F1^2-1.0d)
    268 ##         T2=3.14159/sqrt(F2^2-1.0d)
    269 ##         Y4=T1+T2
    270 ##         Y5=F1^2*T1+F2^2*T2-3.14159*(F1+F2)
    271 ##         Y6=F1^4*T1+F2^4*T2-3.14159*(F1+F2)/2.0-3.14159*(F1^3+F2^3)
    272 ##         Y7=(T2-T1)/(F1-F2)
    273 ##         Y8=Y1+(F2^2*T2-F1^2*T1)/(F1-F2)
    274 ##         Y9=Y2+(F2^4*T2-F1^4*T1)/(F1-F2)+Y1*((F1+F2)^2-F1*F2)
    275 ##         M=(a^2*(X*Y1+Y*Y2+Z*Y3)+a*E*(X*Y4+Y*Y5+Z*Y6)+E^2* $
    276 ##                      (X*Y7+Y*Y8+Z*Y9))*cosgama
    277 ##
    278 ##         Q=Q+W*M
    279 ##
    280 ##      endfor
    281 ##      Q=double(num_steps)/Double(num_steps-1)*Q*G
    282 ##      end
    283 #    elif 'Cylinder' in Geometry:
    284 #        Q = 0.
    285 #    elif 'Fixed' in Geometry:   #Dwiggens & Park, Acta Cryst. A27, 264 (1971) with corrections
    286 #        EMA = np.exp(-MuT*(1.0-SCTth))
    287 #        Fac1 = (1.-EMA)/(1.0-SCTth)
    288 #        G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*Fac1)
    289 #        Fac0 = Cgama/(1-Sgama*SCTth)
    290 #        Wp = Fac0*(Fac1-(EMA-np.exp(-MuT*(CSgama-SCTth)))/(CSgama-1.0))
    291 #        Fac0 = Cgama/(1.0+Sgama*SCTth)
    292 #        Wm = Fac0*(Fac1+(np.exp(-MuT*(1.0+CSgama))-1.0)/(CSgama+1.0))
    293 #        X = (Sgama**2+CTth**2*(1.0-Sgama**2+Sgama**4))/2.0
    294 #        Y = Sgama**3*Cgama*STth*CTth
    295 #        Z = Cgama**2*(1.0+Sgama**2)*STth**2/2.0
    296 #        Fac2 = 1.0/(b*Cgama*STth)
    297 #        U = 2.0*(1.0-a)*Fac2
    298 #        V = (2.0+b*(1.0-CTth*Sgama))*Fac2
    299 #        Mp = 2.0*np.pi*(a+2.0*(1.0-a)/(2.0+b*(1.0-Sgama)))*(a*X+a*Z/2.0-U*Y+U*(X+Y*V+Z*V**2)/np.sqrt(V**2-1.0)-U*Z*V)
    300 #        V = (2.0+b*(1.0+CTth*Sgama))*Fac2
    301 #        Y = -Y
    302 #        Mm = 2.0*np.pi*(a+2.0*(1.0-a)/(2.0+b*(1.0+Sgama)))*(a*X+a*Z/2.0-U*Y+U*(X+Y*V+Z*V**2)/np.sqrt(V**2-1.0)-U*Z*V)
    303 #        Q = np.sum(Wp*Mp+Wm*Mm,axis=0)
    304 #        return Q*G*NS/(NS-1.)
    305 #    elif 'Tilting' in Geometry:
    306 #        EMA = np.exp(-MuT*SCth)
    307 #        G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*EMA)
    308 ##          Wplus = expmutsecth/(1.0d - singama*secth) + singama/mut/(1.0 -singama*secth)/(1.0-singama*secth)* $
    309 ##                                                       (Exp(-1.0*mut*cscgama) - expmutsecth)
    310 ##          Wminus = expmutsecth/(1.0d + singama*secth) - singama/mut/(1.0 +singama*secth)/(1.0+singama*secth)* $
    311 ##                                                        expmutsecth*(1.0d - expmutsecth*Exp(-1.0*mut*cscgama))
    312 #        Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)/(1.0-Sgama*SCth)*(np.exp(-MuT*CSgama)-EMA)
    313 ##        Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)**2*(np.exp(-MuT*CSgama)-EMA)
    314 #        Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)/(1.0+Sgama*SCth)*EMA*(1.0-EMA*np.exp(-MuT*CSgama))
    315 ##        Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)**2*EMA*(1.0-EMA*np.exp(-MuT*CSgama))
    316 #        X = 0.5*(Cth**2*(Cth**2*Sgama**4-4.0*Sth**2*Cgama**2)+1.0)
    317 #        Y = Cgama**2*(1.0+Cgama**2)*Cth**2*Sth**2
    318 #        Z = 0.5*Cgama**4*Sth**4
    319 ##          X = 0.5*(costh*costh*(costh*costh*singama*singama*singama*singama - $
    320 ##                           4.0*sinth*sinth*cosgama*cosgama) +1.0d)
    321 ##
    322 ##          Y = cosgama*cosgama*(1.0 + cosgama*cosgama)*costh*costh*sinth*sinth
    323 ##
    324 ##          Z= 0.5 *cosgama*cosgama*cosgama*cosgama* (sinth^4)
    325 ##
    326 #        AlP = 2.0+b*(1.0-Cth*Sgama)
    327 #        AlM = 2.0+b*(1.0+Cth*Sgama)
    328 ##          alphaplus = 2.0 + b*(1.0 - costh*singama)
    329 ##          alphaminus = 2.0 + b*(1.0 + costh*singama)
    330 #        BeP = np.sqrt(np.fabs(AlP**2-(b*Cgama*Sth)**2))
    331 #        BeM = np.sqrt(np.fabs(AlM**2-(b*Cgama*Sth)**2))
    332 ##          betaplus = Sqrt(Abs(alphaplus*alphaplus - b*b*cosgama*cosgama*sinth*sinth))
    333 ##          betaminus = Sqrt(Abs(alphaminus*alphaminus - b*b*cosgama*cosgama*sinth*sinth))
    334 #        Mp = Cgama*(np.pi*a**2*(2.0*X+Y+0.75*Z)+(2.0*np.pi*(1.0-a))*(1.0-a+a*AlP)* \
    335 #            (4.0*X/AlP/BeP+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlP/BeP-1.0)+
    336 #            2.0/b**4*AlP/BeP*AlP**2-2.0/b**4*AlP**2-Cgama**2/b/b*Sth*2))
    337 ##          Mplus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + $
    338 ##                   (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaplus)* $
    339 ##                   (4.0*x/alphaplus/betaplus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaplus/betaplus -1.0) + $
    340 ##                   2.0/(b^4)*alphaplus/betaplus*alphaplus*alphaplus - 2.0/(b^4)*alphaplus*alphaplus - $
    341 ##                   cosgama*cosgama/b/b*sinth*sinth))
    342 #        Mm =Cgama*(np.pi*a**2*(2.0*X+Y+0.75*Z)+(2.0*np.pi*(1.0-a))*(1.0-a+a*AlM)* \
    343 #            (4.0*X/AlM/BeM+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlM/BeM-1.0)+
    344 #            2.0/b**4*AlM/BeM*AlM**2-2.0/b**4*AlM**2-Cgama**2/b/b*Sth*2))
    345 ##          Mminus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + $
    346 ##                   (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaminus)* $
    347 ##                   (4.0*x/alphaminus/betaminus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaminus/betaminus -1.0) + $
    348 ##                   2.0/(b^4)*alphaminus/betaminus*alphaminus*alphaminus - 2.0/(b^4)*alphaminus*alphaminus - $
    349 ##                   cosgama*cosgama/b/b*sinth*sinth))
    350 #        Q = np.sum(Wp*Mp+Wm*Mm,axis=0)
    351 #        return Q*G*NS/(NS-1.)
    352 ##       expmutsecth = Exp(-1.0*mut*secth)
    353 ##       G= delta * 2.0 * Na * emc /(1.0+costth^2)/expmutsecth
    354 ##       for i=1, num_steps-1 do begin
    355 ##          cosgama=double(cos(gama[i]))
    356 ##          singama=double(sin(gama[i]))
    357 ##          cscgama=1.0d/singama
    358 ##
    359 ##
    360 ##     ; print, "W", min(wplus), max(wplus), min(wminus), max(wminus)
    361 ##
    362 ##
    363 ##
    364 ##
    365 ##    ;               print, a,b
    366 ##  ; print, "M", min(mplus), max(mplus), min(mminus), max(mminus)
    367 ##          Q=Q+ Wplus*Mplus + Wminus*Mminus
    368 ##      endfor
    369 ##      Q=double(num_steps)/double(num_steps-1)*Q*G
    370 ##   ;   print, min(q), max(q)
    371 ##      end
    372 
    373 #def MultiScattering(Geometry,ElList,Tth):
    374 #    BN = BD = 0.0
    375 #    Amu = 0.0
    376 #    for El in ElList:
    377 #        el = ElList[El]
    378 #        BN += el['Z']*el['FormulaNo']
    379 #        BD += el['FormulaNo']
    380 #        Amu += el['FormulaNo']*el['mu']
    381        
    382198def CalcPDF(data,inst,xydata):
    383199    auxPlot = []
    … …  
    913729                return yb+yc
    914730    else:
     731        Pdabc = parmDict['Pdabc']
    915732        difC = parmDict['difC']
    916733        alp0 = parmDict['alpha']
    917734        bet0 = parmDict['beta-0']
    918735        bet1 = parmDict['beta-1']
    919         sig0 = parmDict['var-inst']
     736        sig0 = parmDict['sig-0']
     737        sig1 = parmDict['sig-1']
    920738        X = parmDict['X']
    921739        Y = parmDict['Y']
    … …  
    931749                    alp = parmDict[alpName]
    932750                else:
    933                     alp = alp0/dsp
     751                    if len(Pdabc):
     752                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
     753                    else:
     754                        alp = alp0/dsp
    934755                betName = 'bet'+str(iPeak)
    935756                if betName in varyList:
    936757                    bet = parmDict[betName]
    937758                else:
    938                     bet = bet0+bet1/dsp**4
     759                    if len(Pdabc):
     760                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
     761                    else:
     762                        bet = bet0+bet1/dsp**4
    939763                sigName = 'sig'+str(iPeak)
    940764                if sigName in varyList:
    941765                    sig = parmDict[sigName]
    942766                else:
    943                     sig = sig0*dsp**2
     767                    sig = sig0+sig1*dsp**2
    944768                gamName = 'gam'+str(iPeak)
    945769                if gamName in varyList:
    … …  
    1080904                break
    1081905    else:
     906        Pdabc = parmDict['Pdabc']
    1082907        difC = parmDict['difC']
    1083908        alp0 = parmDict['alpha']
    1084909        bet0 = parmDict['beta-0']
    1085910        bet1 = parmDict['beta-1']
    1086         sig0 = parmDict['var-inst']
     911        sig0 = parmDict['sig-0']
     912        sig1 = parmDict['sig-1']
    1087913        X = parmDict['X']
    1088914        Y = parmDict['Y']
    … …  
    1098924                    alp = parmDict[alpName]
    1099925                else:
    1100                     alp = alp0/dsp
     926                    if len(Pdabc):
     927                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
     928                    else:
     929                        alp = alp0/dsp
    1101930                    dada0 = 1./dsp
    1102931                betName = 'bet'+str(iPeak)
    … …  
    1104933                    bet = parmDict[betName]
    1105934                else:
    1106                     bet = bet0+bet1/dsp**4
     935                    if len(Pdabc):
     936                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
     937                    else:
     938                        bet = bet0+bet1/dsp**4
    1107939                    dbdb0 = 1.
    1108940                    dbdb1 = 1/dsp**4
    … …  
    1111943                    sig = parmDict[sigName]
    1112944                else:
    1113                     sig = sig0*dsp**2
    1114                     dsds0 = dsp**2
     945                    sig = sig0+sig1*dsp**2
     946                    dsds0 = 1.
     947                    dsds1 = dsp**2
    1115948                gamName = 'gam'+str(iPeak)
    1116949                if gamName in varyList:
    … …  
    1153986                if 'beta-1' in varyList:
    1154987                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
    1155                 if 'var-inst' in varyList:
    1156                     dMdv[varyList.index('var-inst')] += dsds0*dervDict['sig']
     988                if 'sig-0' in varyList:
     989                    dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
     990                if 'sig-1' in varyList:
     991                    dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
    1157992                if 'X' in varyList:
    1158993                    dMdv[varyList.index('X')] += dsdX*dervDict['gam']
    … …  
    13051140            insVals.append(Inst[parm][1])
    13061141            if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
    1307                 'beta-0','beta-1','var-inst',] and Inst[parm][2]:
     1142                'beta-0','beta-1','sig-0','sig-1',] and Inst[parm][2]:
    13081143                    insVary.append(parm)
    13091144        instDict = dict(zip(insNames,insVals))
    … …  
    13271162                    else:
    13281163                        dsp = pos/Inst['difC'][1]
    1329                         parmDict[sigName] = parmDict['var-inst']*dsp**2
     1164                        parmDict[sigName] = parmDict['sig-0']+parmDict['sig-1']*dsp**2
    13301165                gamName = 'gam'+str(iPeak)
    13311166                if gamName not in varyList:
    … …  
    13471182        for parm in Inst:
    13481183            if parm in  ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
    1349                 'beta-0','beta-1','var-inst',]:
     1184                'beta-0','beta-1','sig-0','sig-1',]:
    13501185                ptlbls += "%s" % (parm.center(12))
    13511186                ptstr += ptfmt % (Inst[parm][1])
    … …  
    13651200            names = ['pos','int','sig','gam']
    13661201        else:
    1367             names = ['pos','int','alpha','beta','sig','gam']
     1202            names = ['pos','int','alp','bet','sig','gam']
    13681203        for i,peak in enumerate(Peaks):
    13691204            for j,name in enumerate(names):
    … …  
    13791214            names = ['pos','int','sig','gam']
    13801215        else:
    1381             names = ['pos','int','alpha','beta','sig','gam']
     1216            names = ['pos','int','alp','bet','sig','gam']
    13821217        for i,peak in enumerate(Peaks):
    13831218            pos = parmDict['pos'+str(i)]
    … …  
    13961231                        peak[2*j] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
    13971232                    else:
    1398                         peak[2*j] = parmDict['var-inst']*dsp**2
     1233                        peak[2*j] = parmDict['sig-0']+parmDict['sig-1']*dsp**2
    13991234                elif 'gam' in parName:
    14001235                    if 'C' in Inst['Type'][0]:
    … …  
    14081243            names = ['pos','int','sig','gam']
    14091244        else:
    1410             names = ['pos','int','alpha','beta','sig','gam']           
     1245            names = ['pos','int','alp','bet','sig','gam']           
    14111246        head = 13*' '
    14121247        for name in names:
    … …  
    14161251            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
    14171252        else:
    1418             ptfmt = {'pos':"%10.2f",'int':"%10.1f",'alpha':"%10.3f",'beta':"%10.5f",'sig':"%10.3f",'gam':"%10.3f"}
     1253            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%10.3f",'bet':"%10.5f",'sig':"%10.3f",'gam':"%10.3f"}
    14191254        for i,peak in enumerate(Peaks):
    14201255            ptstr =  ':'
    … …  
    14541289        Ftol = 1.0
    14551290    x,y,w,yc,yb,yd = data               #these are numpy arrays!
     1291    yc *= 0.                            #set calcd ones to zero
     1292    yb *= 0.
     1293    yd *= 0.
    14561294    xBeg = np.searchsorted(x,Limits[0])
    14571295    xFin = np.searchsorted(x,Limits[1])
    … …  
    14631301    parmDict.update(insDict)
    14641302    parmDict.update(peakDict)
     1303    parmDict['Pdabc'] = []      #dummy Pdabc
     1304    parmDict.update(Inst2)      #put in real one if there
    14651305    varyList = bakVary+insVary+peakVary
    14661306    while True:

  • trunk/GSASIIpwdGUI.py

    r796 r798  
    174174        PatternId = G2frame.PatternId
    175175        PickId = G2frame.PickId
    176         Inst = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Instrument Parameters'))[0]
     176        Inst,Inst2 = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Instrument Parameters'))
    177177        peaks = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Peak List'))
    178178        if not peaks:
    179179            G2frame.ErrorDialog('No peaks!','Nothing to do!')
    180180            return
     181        newpeaks = []
    181182        for peak in peaks:
    182             if 'C' in Inst['Type'][0]:
    183                 peak[4] = Inst['U'][1]*tand(peak[0]/2.0)**2+Inst['V'][1]*tand(peak[0]/2.0)+Inst['W'][1]
    184                 peak[6] = Inst['X'][1]/cosd(peak[0]/2.0)+Inst['Y'][1]*tand(peak[0]/2.0)
    185             else:
    186                 dsp = peak[0]/Inst['difC'][0]
    187                 peak[4] = Inst['alpha'][0]*dsp
    188                 peak[6] = Inst['beta-0'][0]+Inst['beta-1'][0]/dsp**4
    189                 peak[8] = Inst['var-inst'][0]*dsp**2
    190                 peak[10] = Inst['X'][0]*dsp**2+Inst['Y'][0]*dsp
    191         UpdatePeakGrid(G2frame,peaks)
     183            newpeaks.append(G2mth.setPeakparms(Inst,Inst2,peak[0],peak[2]))
     184        G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Peak List'),newpeaks)
     185        UpdatePeakGrid(G2frame,newpeaks)
    192186               
    193187    def RefreshPeakGrid(event):
    … …  
    292286        colLabels = ['position','refine','intensity','refine','alpha','refine',
    293287            'beta','refine','sigma','refine','gamma','refine']
    294         Types = [wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_BOOL,
    295             wg.GRID_VALUE_FLOAT+':10,1',wg.GRID_VALUE_BOOL,
     288        Types = [wg.GRID_VALUE_FLOAT+':10,1',wg.GRID_VALUE_BOOL,
     289            wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_BOOL,
    296290            wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_BOOL,
    297291            wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_BOOL,
    … …  
    653647        for key in keys:
    654648            if key in ['Type','U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
    655                 'beta-0','beta-1','var-inst','Polariz.','Lam','Azimuth','2-theta',
     649                'beta-0','beta-1','sig-0','sig-1','Polariz.','Lam','Azimuth','2-theta',
    656650                'difC','difA','Zero']:
    657651                good.append(key)
    … …  
    676670        'FeKa':[1.93597,1.93991],'CoKa':[1.78892,1.79278],'MoKa':[0.70926,0.713543],
    677671        'AgKa':[0.559363,0.563775]}
     672    Inst2 = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,
     673            G2frame.PatternId,'Instrument Parameters'))[1]
    678674       
    679675    def inst2data(inst,ref,data):
    … …  
    689685        if doAnyway or event.GetRow() == 1:
    690686            peaks = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Peak List'))
    691             if 'C' in data['Type'][0]:                                       #update powder peak parameters
    692                 for peak in peaks:
    693                     peak[4] = insVal['U']*tand(peak[0]/2.0)**2+insVal['V']*tand(peak[0]/2.0)+insVal['W']
    694                     peak[6] = insVal['X']/cosd(peak[0]/2.0)+insVal['Y']*tand(peak[0]/2.0)
    695             else:
    696                 for peak in peaks:
    697                     dsp = peak[0]/insVal['difC']
    698                     peak[4] = insVal['alpha']*dsp
    699                     peak[6] = insVal['beta-0']+insVal['beta-1']/dsp**4
    700                     peak[8] = insVal['var-inst']*dsp**2
    701                     peak[10] = insVal['X']*dsp**2+insVal['Y']*dsp
     687            newpeaks = []
     688            for peak in peaks:
     689                newpeaks.append(G2mth.setPeakparms(data,Inst2,peak[0],peak[2]))
     690            G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Peak List'),newpeaks)
    702691                   
    703692    def OnLoad(event):
    … …  
    728717                    S = File.readline()               
    729718                File.close()
     719                Inst,Inst2 = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId,'Instrument Parameters'))
    730720                inst = G2IO.makeInstDict(newItems,newVals,len(newVals)*[False,])
    731                 G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId,'Instrument Parameters'),data)
     721                G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId,'Instrument Parameters'),[inst,Inst2])
    732722                RefreshInstrumentGrid(event,doAnyway=True)          #to get peaks updated
    733723                UpdateInstrumentGrid(G2frame,data)
    … …  
    996986            instSizer.Add((5,5),0)
    997987            instSizer.Add((5,5),0)
    998             for item in ['difC','difA','Zero','alpha','beta-0','beta-1','var-inst','X','Y']:
     988            for item in ['difC','difA','Zero','alpha','beta-0','beta-1','sig-0','sig-1','X','Y']:
    999989                fmt = '%10.3f'
    1000990                if 'beta' in item:
Note: See TracChangeset for help on using the changeset viewer.