Changeset 376


Ignore:
Timestamp:
Sep 20, 2011 1:42:17 PM (10 years ago)
Author:
vondreele
Message:

fix intensity type derivatives - scale, polarization, pref. ori., etc.

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIpwd.py

    r375 r376  
    128128       
    129129def Polarization(Pola,Tth,Azm=0.0):
    130 #   Calculate angle dependent x-ray polarization correction (not scaled correctly!)
    131 #   Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
    132 #   Azm: azimuthal angle e.g. 0.0 in plane of polarization
    133 #   Tth: 2-theta scattering angle - can be numpy array
    134 #       which (if either) of these is "right"?
    135 #    return (Pola*npcosd(Azm)**2+(1.-Pola)*npsind(Azm)**2)*npcosd(Tth)**2+ \
    136 #        Pola*npsind(Azm)**2+(1.-Pola)*npcosd(Azm)**2
     130    """   Calculate angle dependent x-ray polarization correction (not scaled correctly!)
     131    input:
     132        Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
     133        Azm: azimuthal angle e.g. 0.0 in plane of polarization
     134        Tth: 2-theta scattering angle - can be numpy array
     135            which (if either) of these is "right"?
     136    return:
     137        pola = (Pola*npcosd(Azm)**2+(1.-Pola)*npsind(Azm)**2)*npcosd(Tth)**2+ \
     138            Pola*npsind(Azm)**2+(1.-Pola)*npcosd(Azm)**2
     139        dpdPola: derivative needed for least squares
     140    """
    137141    pola = ((1.0-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+   \
    138142        (1.0-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
    139143    dpdPola = -npsind(Tth)**2*(npsind(Azm)**2-npcosd(Azm)**2)
    140     return pola,dpdPola/pola
     144    return pola,dpdPola
    141145   
    142146def Oblique(ObCoeff,Tth):
  • trunk/GSASIIstruct.py

    r375 r376  
    169169        Histograms = dictionary of histograms as {name:data,...}
    170170        Phases = dictionary of phases that use histograms
    171         CovData = [varyList,covariance matrix]
     171        CovData = dictionary of refined variables, varyList, & covariance matrix
    172172    '''
    173173                       
     
    205205                    data[iphase][1] = Phases[phaseName]
    206206        elif datum[0] == 'Covariance':
    207             varyList,covMatrix = CovData
    208             data[0][1] = {'varyList':varyList,'covariance':covMatrix}
     207            data[0][1] = CovData
    209208        try:
    210209            histogram = Histograms[datum[0]]
     
    11251124    parmdict.update(zip(varylist,values))
    11261125   
     1126def GetPrefOri(refl,G,phfx,calcControls,parmDict):
     1127    if calcControls[phfx+'poType'] == 'MD':
     1128        MD = parmDict[phfx+'MD']
     1129        MDAxis = calcControls[phfx+'MDAxis']
     1130        sumMD = 0
     1131        for H in refl[10]:           
     1132            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
     1133            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
     1134            sumMD += A**3
     1135        POcorr = sumMD/len(refl[10])
     1136    else:   #spherical harmonics
     1137        POcorr = 1.0
     1138    return POcorr
     1139   
     1140def GetPrefOriDerv(refl,G,phfx,calcControls,parmDict):
     1141    POderv = {}
     1142    if calcControls[phfx+'poType'] == 'MD':
     1143        MD = parmDict[phfx+'MD']
     1144        MDAxis = calcControls[phfx+'MDAxis']
     1145        sumMD = 0
     1146        sumdMD = 0
     1147        for H in refl[10]:           
     1148            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
     1149            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
     1150            sumMD += A**3
     1151            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
     1152        POcorr = sumMD/len(refl[10])
     1153        POderv[phfx+'MD'] = sumdMD/len(refl[10])
     1154    else:   #spherical harmonics
     1155        POcorr = 1.0
     1156    return POcorr,POderv
     1157   
     1158def GetIntensityCorr(refl,G,phfx,hfx,calcControls,parmDict):
     1159    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
     1160    Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
     1161    Icorr *= GetPrefOri(refl,G,phfx,calcControls,parmDict)       
     1162    return Icorr
     1163   
     1164def GetIntensityDerv(refl,G,phfx,hfx,calcControls,parmDict):
     1165    Icorr = GetIntensityCorr(refl,G,phfx,hfx,calcControls,parmDict)
     1166    pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
     1167    POcorr,dIdPO = GetPrefOriDerv(refl,G,phfx,calcControls,parmDict)
     1168    dIdPola *= refl[3]/pola
     1169    for iPO in dIdPO:
     1170        dIdPO[iPO] *= refl[3]/POcorr
     1171    dIdsh = refl[3]/parmDict[hfx+'Scale']
     1172    dIdsp = refl[3]/parmDict[phfx+'Scale']
     1173    return Icorr,dIdsh,dIdsp,dIdPola,dIdPO
     1174       
    11271175def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
    11281176   
     
    11581206            gam += 0.018*refl[4]**2*tand(refl[5]/2.)*sum           
    11591207        return gam
    1160        
    1161     def GetMarchDollase(refl,G,phfx,calcControls,parmDict):
    1162         MD = parmDict[phfx+'MD']
    1163         MDAxis = calcControls[phfx+'MDAxis']
    1164         sumMD = 0
    1165         for H in refl[10]:           
    1166             cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
    1167             sumMD += 1/np.sqrt((MD*cosP)**2+sinP**2/MD)**3
    1168         return sumMD/len(refl[10])
    1169        
    1170     def GetIntensityCorr(refl,G,phfx,hfx,calcControls,parmDict):
    1171         Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
    1172         Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
    1173         if calcControls[phfx+'poType'] == 'MD':
    1174             Icorr *= GetMarchDollase(refl,G,phfx,calcControls,parmDict)
    1175        
    1176         return Icorr
    11771208       
    11781209    def GetHStrainShift(refl,SGData,phfx,parmDict):
     
    12821313                    continue
    12831314                elif not iBeg-iFin:     #peak above high limit - done
    1284                     return yc,yb
     1315                    break
    12851316                yc[iBeg:iFin] += Icorr*refl[8]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
    12861317                if Ka2:
     
    12941325                        return yc,yb
    12951326                    yc[iBeg:iFin] += Icorr*refl[8]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
    1296             else:
    1297                 raise ValueError
     1327            elif 'T' in calcControls[hfx+'histType']:
     1328                print 'TOF Undefined at present'
     1329                raise ValueError    #no TOF yet
    12981330    return yc,yb   
    12991331           
     
    13411373                gamDict[phfx+'Mustrain:'+str(i)] = const*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
    13421374        return gamDict
    1343        
    1344     def GetMarchDollaseDerv(refl,G,phfx,calcControls,parmDict):
    1345         MD = parmDict[phfx+'MD']
    1346         MDAxis = calcControls[phfx+'MDAxis']
    1347         sumMD = 0
    1348         sumdMD = 0
    1349         for H in refl[10]:           
    1350             cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
    1351             A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
    1352             sumMD += A**3
    1353             sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
    1354         return sumMD/len(refl[10]),refl[8]*sumdMD
    1355        
    1356     def GetIntensityDerv(refl,G,phfx,hfx,calcControls,parmDict):
    1357         Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
    1358         pola,dpdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
    1359         dIdpola = refl[8]*Icorr*dpdPola/refl[3]
    1360         Icorr *= pola
    1361         MDcorr,dIdMD = GetMarchDollaseDerv(refl,G,phfx,calcControls,parmDict)
    1362         Icorr *= MDcorr
    1363         dIdsh = Icorr/parmDict[hfx+'Scale']
    1364         dIdsp = Icorr/parmDict[phfx+'Scale']
    1365         dIdMD *= Icorr
    1366         return Icorr,dIdsh,dIdsp,dIdpola,dIdMD
    13671375       
    13681376    def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
     
    14761484            sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)])
    14771485        for refl in refList:
    1478             if 'C' in calcControls[hfx+'histType']:
     1486            if 'C' in calcControls[hfx+'histType']:        #CW powder
    14791487                h,k,l = refl[:3]
    14801488                iref = pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]
    1481                 Icorr,dIdsh,dIdsp,dIdpola,dIdMD = GetIntensityDerv(refl,G,phfx,hfx,calcControls,parmDict)
    1482                 hkl = refl[:3]
     1489                Icorr,dIdsh,dIdsp,dIdpola,dIdPO = GetIntensityDerv(refl,G,phfx,hfx,calcControls,parmDict)
     1490                if 'Pawley' in Phase['General']['Type']:
     1491                    try:
     1492                        refl[8] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
     1493                    except KeyError:
     1494#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
     1495                        continue
     1496                else:
     1497                    raise ValueError       #wants strctrfacr deriv calc here
     1498                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
     1499                iBeg = np.searchsorted(x,refl[5]-fmin)
     1500                iFin = np.searchsorted(x,refl[5]+fmax)
     1501                if not iBeg+iFin:       #peak below low limit - skip peak
     1502                    continue
     1503                elif not iBeg-iFin:     #peak above high limit - done
     1504                    break
    14831505                pos = refl[5]
    14841506                tanth = tand(pos/2.0)
    14851507                costh = cosd(pos/2.0)
    1486                 dsdU = tanth**2
    1487                 dsdV = tanth
    1488                 dsdW = 1.0
    1489                 dgdX = 1.0/costh
    1490                 dgdY = tanth
    1491                 Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
    1492                 iBeg = np.searchsorted(x,refl[5]-fmin)
    1493                 iFin = np.searchsorted(x,refl[5]+fmax)
    14941508                dMdpk = np.zeros(shape=(6,len(x)))
    14951509                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
    14961510                for i in range(1,5):
    14971511                    dMdpk[i][iBeg:iFin] += 100.*dx*Icorr*refl[8]*dMdipk[i]
    1498                 dMdpk[0][iBeg:iFin] += 100.*dx*Icorr*dMdipk[0]
     1512                dMdpk[0][iBeg:iFin] += 100.*dx*Icorr*refl[8]*dMdipk[0]/refl[3]
    14991513                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
    15001514                if Ka2:
    1501                     pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
     1515                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
    15021516                    kdelt = int((pos2-refl[5])/dx)               
    1503                     iBeg = min(lenX,iBeg+kdelt)
    1504                     iFin = min(lenX,iFin+kdelt)
    1505                     if iBeg-iFin:
    1506                         dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])
     1517                    iBeg2 = min(lenX,iBeg+kdelt)
     1518                    iFin2 = min(lenX,iFin+kdelt)
     1519                    if iBeg2-iFin2:
     1520                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
    15071521                        for i in range(1,5):
    1508                             dMdpk[i][iBeg:iFin] += 100.*dx*Icorr*refl[8]*kRatio*dMdipk2[i]
    1509                         dMdpk[0][iBeg:iFin] += 100.*dx*kRatio*dMdipk2[0]
    1510                         dMdpk[5][iBeg:iFin] += 100.*dx*dMdipk2[0]
    1511                         dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*Icorr*refl[8]}
     1522                            dMdpk[i][iBeg2:iFin2] += 100.*dx*Icorr*refl[8]*kRatio*dMdipk2[i]
     1523                        dMdpk[0][iBeg2:iFin2] += 100.*dx*Icorr*refl[8]*kRatio*dMdipk2[0]/refl[3]
     1524                        dMdpk[5][iBeg2:iFin2] += 100.*dx*Icorr*dMdipk2[0]
     1525                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*refl[8]}
    15121526                try:
    15131527                    idx = varylist.index(pfx+'PWLref:'+str(iref))
     
    15171531                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
    15181532                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
    1519                     hfx+'U':[dsdU,'sig'],hfx+'V':[dsdV,'sig'],hfx+'W':[dsdW,'sig'],
    1520                     hfx+'X':[dgdX,'gam'],hfx+'Y':[dgdY,'gam'],hfx+'SH/L':[1.0,'shl'],
     1533                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
     1534                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
    15211535                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
    15221536                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
    1523                     hfx+'DisplaceY':[dpdY,'pos'],phfx+'MD':[dIdMD,'int'],}
     1537                    hfx+'DisplaceY':[dpdY,'pos'],}
    15241538                for name in names:
    15251539                    if name in varylist:
    15261540                        item = names[name]
    15271541                        dMdv[varylist.index(name)] += item[0]*dervDict[item[1]]
     1542                for iPO in dIdPO:
     1543                    if iPO in varylist:
     1544                        dMdv[varylist.index(iPO)] += dIdPO[iPO]*dervDict['int']
    15281545                cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
    15291546                for name,dpdA in cellDervNames:
     
    15381555                    if name in varylist:
    15391556                        dMdv[varylist.index(name)] += gamDict[name]*dervDict['gam']                       
    1540             else:
    1541                 raise ValueError
     1557            elif 'T' in calcControls[hfx+'histType']:
     1558                print 'TOF Undefined at present'
     1559                raise ValueError    #no TOF yet
    15421560           
    15431561    return dMdv   
     
    16861704#    print 'indParmList: ',G2mv.indParmList
    16871705#    print 'fixedDict: ',G2mv.fixedDict
    1688     covFile = ospath.splitext(GPXfile)[0]+'.gpxcov'
    1689     file = open(covFile,'wb')
    1690     cPickle.dump(varyList,file,1)
    1691     cPickle.dump(result[0],file,1)
    1692     cPickle.dump(covMatrix,file,1)
    1693     file.close()
    16941706    sigDict = dict(zip(varyList,sig))
    16951707    SetPhaseData(parmDict,sigDict,Phases)
    16961708    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
    16971709    SetHistogramData(parmDict,sigDict,Histograms)
    1698     SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,[varyList,cov])
     1710    covData = {'variables':result[0],'varyList':varyList,'covariance':cov}
     1711    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData)
    16991712#for testing purposes!!!
    17001713#    file = open('structTestdata.dat','wb')
Note: See TracChangeset for help on using the changeset viewer.