Changeset 1484


Ignore:
Timestamp:
Sep 3, 2014 10:34:01 AM (9 years ago)
Author:
vondreele
Message:

change lower limit of LGmix to 0 from 0.1 - recompile Win-32 & Win-64 fortran sources
fix ellipseSizeDeriv - didn't do it correctly before
fix mustrain & size fxns & derivatives for CW (minor) & TOF major; all checked against numerical derivatives for both lorentzian & gaussian contributions as well as mixing coeff.

Location:
trunk
Files:
19 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIpwd.py

    r1474 r1484  
    788788    dRdS = np.zeros(6)
    789789    for i in range(6):
    790         dSij = Sij[:]
    791         dSij[i] += delt
    792         dRdS[i] = (ellipseSize(H,dSij,GB)-lenR)/delt
     790        Sij[i] -= delt
     791        lenM = ellipseSize(H,Sij,GB)
     792        Sij[i] += 2.*delt
     793        lenP = ellipseSize(H,Sij,GB)
     794        Sij[i] -= delt
     795        dRdS[i] = (lenP-lenM)/(2.*delt)
    793796    return lenR,dRdS
    794797

  • trunk/GSASIIstrMath.py

    r1482 r1484  
    11591159def GetSampleSigGam(refl,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
    11601160    'Needs a doc string'
    1161     if 'C' in calcControls[hfx+'histType']:
     1161    if 'C' in calcControls[hfx+'histType']:     #All checked & OK
    11621162        costh = cosd(refl[5]/2.)
    11631163        #crystallite size
    … …  
    11911191                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
    11921192            Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*np.sqrt(Sum)/np.pi
    1193     elif 'T' in calcControls[hfx+'histType']:
     1193    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
    11941194        #crystallite size
    1195         if calcControls[phfx+'SizeType'] == 'isotropic':
    1196             Sgam = 1.e-4*parmDict[hfx+'difC']*parmDict[phfx+'Size;i']
    1197         elif calcControls[phfx+'SizeType'] == 'uniaxial':
     1195        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
     1196            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4]**2/parmDict[phfx+'Size;i']
     1197        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
    11981198            H = np.array(refl[:3])
    11991199            P = np.array(calcControls[phfx+'SizeAxis'])
    12001200            cosP,sinP = G2lat.CosSinAngle(H,P,G)
    1201             Sgam = 1.e-4*parmDict[hfx+'difC']*(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
    1202             Sgam /= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
    1203         else:           #ellipsoidal crystallites
     1201            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
     1202            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
     1203        else:           #ellipsoidal crystallites   #OK
    12041204            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
    12051205            H = np.array(refl[:3])
    12061206            lenR = G2pwd.ellipseSize(H,Sij,GB)
    1207             Sgam = 1.e-4*parmDict[hfx+'difC']*(refl[4]**2*lenR)
     1207            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4]**2/lenR
    12081208        #microstrain               
    1209         if calcControls[phfx+'MustrainType'] == 'isotropic':
     1209        if calcControls[phfx+'MustrainType'] == 'isotropic':    #OK
    12101210            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4]*parmDict[phfx+'Mustrain;i']
    1211         elif calcControls[phfx+'MustrainType'] == 'uniaxial':
     1211        elif calcControls[phfx+'MustrainType'] == 'uniaxial':   #OK
    12121212            H = np.array(refl[:3])
    12131213            P = np.array(calcControls[phfx+'MustrainAxis'])
    … …  
    12161216            Sa = parmDict[phfx+'Mustrain;a']
    12171217            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
    1218         else:       #generalized - P.W. Stephens model
     1218        else:       #generalized - P.W. Stephens model  OK
     1219            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
    12191220            Sum = 0
    1220             Strms = G2spc.MustrainCoeff(refl[:3],SGData)
    12211221            for i,strm in enumerate(Strms):
    12221222                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
    1223             Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4]**2*Sum
     1223            Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4]**3
    12241224           
    12251225    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
    … …  
    12321232    gamDict = {}
    12331233    sigDict = {}
    1234     if 'C' in calcControls[hfx+'histType']:
     1234    if 'C' in calcControls[hfx+'histType']:         #All checked & OK
    12351235        costh = cosd(refl[5]/2.)
    12361236        tanth = tand(refl[5]/2.)
    … …  
    12461246            Si = parmDict[phfx+'Size;i']
    12471247            Sa = parmDict[phfx+'Size;a']
    1248             gami = (1.8*wave/np.pi)/(Si*Sa)
     1248            gami = 1.8*wave/(costh*np.pi*Si*Sa)
    12491249            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
    12501250            Sgam = gami*sqtrm
    1251             gam = Sgam/costh
    1252             dsi = (gami*Si*cosP**2/(sqtrm*costh)-gam/Si)
    1253             dsa = (gami*Sa*sinP**2/(sqtrm*costh)-gam/Sa)
     1251            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
     1252            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
    12541253            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
    12551254            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
    … …  
    13021301        gamDict[phfx+'Mustrain;mx'] = Mgam
    13031302        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
    1304     else:   #'T'OF
    1305         if calcControls[phfx+'SizeType'] == 'isotropic':
    1306             Sgam = 1.e-4*parmDict[hfx+'difC']*parmDict[phfx+'Size;i']
    1307             gamDict[phfx+'Size;i'] = 1.e-4*parmDict[hfx+'difC']*parmDict[phfx+'Size;mx']
    1308             sigDict[phfx+'Size;i'] = 2.e-4*parmDict[hfx+'difC']*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
    1309         elif calcControls[phfx+'SizeType'] == 'uniaxial':
    1310             const = 1.e-4*refl[4]*parmDict[hfx+'difC']
     1303    else:   #'T'OF - All checked & OK
     1304        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
     1305            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4]**2/parmDict[phfx+'Size;i']
     1306            gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
     1307            sigDict[phfx+'Size;i'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])**2/(ateln2*parmDict[phfx+'Size;i'])
     1308        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
     1309            const = 1.e-4*parmDict[hfx+'difC']*refl[4]**2
    13111310            H = np.array(refl[:3])
    13121311            P = np.array(calcControls[phfx+'SizeAxis'])
    … …  
    13141313            Si = parmDict[phfx+'Size;i']
    13151314            Sa = parmDict[phfx+'Size;a']
    1316             gami = const*(Si*Sa)
     1315            gami = const/(Si*Sa)
    13171316            sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
    1318             Sgam = gami/sqtrm
    1319             dsi = -gami*Si*cosP**2/sqtrm**3
    1320             dsa = -gami*Sa*sinP**2/sqtrm**3
    1321             gamDict[phfx+'Size;i'] = const*parmDict[phfx+'Size;mx']*Sa/8.
    1322             gamDict[phfx+'Size;a'] = const*parmDict[phfx+'Size;mx']*Si/8.
     1317            Sgam = gami*sqtrm
     1318            dsi = gami*Si*cosP**2/sqtrm-Sgam/Si
     1319            dsa = gami*Sa*sinP**2/sqtrm-Sgam/Sa
     1320            gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
     1321            gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
    13231322            sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
    13241323            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
    1325         else:           #ellipsoidal crystallites
    1326             const = 1.e-4*parmDict[hfx+'difC']
     1324        else:           #OK  ellipsoidal crystallites
     1325            const = 1.e-4*parmDict[hfx+'difC']*refl[4]**2
    13271326            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
    13281327            H = np.array(refl[:3])
    … …  
    13321331                gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
    13331332                sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
    1334         gamDict[phfx+'Size;mx'] = Sgam
    1335         sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
     1333        gamDict[phfx+'Size;mx'] = Sgam  #OK
     1334        sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2  #OK
    13361335               
    13371336        #microstrain derivatives               
    13381337        if calcControls[phfx+'MustrainType'] == 'isotropic':
    13391338            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4]*parmDict[phfx+'Mustrain;i']
    1340             gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']
    1341             sigDict[phfx+'Mustrain;i'] =  2.e-6*refl[4]*parmDict[hfx+'difC']*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
     1339            gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
     1340            sigDict[phfx+'Mustrain;i'] =  2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])**2/(ateln2*parmDict[phfx+'Mustrain;i'])       
    13421341        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
    13431342            H = np.array(refl[:3])
    … …  
    13461345            Si = parmDict[phfx+'Mustrain;i']
    13471346            Sa = parmDict[phfx+'Mustrain;a']
    1348             gami = 1.e-6*parmDict[hfx+'difC']*Si*Sa
     1347            gami = 1.e-6*parmDict[hfx+'difC']*refl[4]*Si*Sa
    13491348            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
    13501349            Mgam = gami/sqtrm
    13511350            dsi = -gami*Si*cosP**2/sqtrm**3
    13521351            dsa = -gami*Sa*sinP**2/sqtrm**3
    1353             gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']*refl[4]
    1354             gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']*refl[4]
    1355             sigDict[phfx+'Mustrain;i'] = 2*refl[4]*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
    1356             sigDict[phfx+'Mustrain;a'] = 2*refl[4]*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
    1357         else:       #generalized - P.W. Stephens model
     1352            gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
     1353            gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
     1354            sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
     1355            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
     1356        else:       #generalized - P.W. Stephens model OK
    13581357            pwrs = calcControls[phfx+'MuPwrs']
    13591358            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
    1360             const = 1.e-6*parmDict[hfx+'difC']*refl[4]**2
    1361             sum = 0
     1359            const = 1.e-6*parmDict[hfx+'difC']*refl[4]**3
     1360            Sum = 0
    13621361            for i,strm in enumerate(Strms):
    1363                 sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
    1364                 gamDict[phfx+'Mustrain:'+str(i)] = const*strm*parmDict[phfx+'Mustrain;mx']
    1365                 sigDict[phfx+'Mustrain:'+str(i)] = \
    1366                     2.*const*term*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
    1367             Mgam = const*sum
     1362                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
     1363                gamDict[phfx+'Mustrain:'+str(i)] = strm*parmDict[phfx+'Mustrain;mx']/2.
     1364                sigDict[phfx+'Mustrain:'+str(i)] = strm*(1.-parmDict[phfx+'Mustrain;mx'])**2
     1365            Mgam = const*np.sqrt(Sum)
    13681366            for i in range(len(Strms)):
    1369                 sigDict[phfx+'Mustrain:'+str(i)] *= Mgam       
     1367                gamDict[phfx+'Mustrain:'+str(i)] *= Mgam/Sum
     1368                sigDict[phfx+'Mustrain:'+str(i)] *= const**2/ateln2       
    13701369        gamDict[phfx+'Mustrain;mx'] = Mgam
    13711370        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
    … …  
    15901589        Ssig,Sgam = GetSampleSigGam(refl,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
    15911590        sig += Ssig
    1592         sig = max(0.001,sig)
    15931591        gam += Sgam
    1594         gam = max(0.001,gam)
    15951592        return sig,gam
    15961593       

  • trunk/fsource/powsubs/epsvoigt.for

    r829 r1484  
    4949!CODE:
    5050
    51       SQSGT = MAX(SQRT(SIG),0.01)
    52       GAM = MAX(GAM,0.1)
     51      SQSGT = MAX(SQRT(SIG),0.00001)
     52      GAM = MAX(GAM,0.00001)
    5353      FWHG = STOFW*SQSGT
    5454      PGL = FWHG**5
Note: See TracChangeset for help on using the changeset viewer.