Changeset 798 for trunk/GSASIIpwd.py


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

more peak fitting work including TOF

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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:
Note: See TracChangeset for help on using the changeset viewer.