Changeset 357


Ignore:
Timestamp:
Aug 30, 2011 9:38:39 AM (11 years ago)
Author:
vondreele
Message:

GSASIIpwd.py - add polarization derivative
GSASIIstruct.py - almost complete derivatives for Pawley refinement checked; lacks derivatives for ellipsoidal size parameters

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIpwd.py

    r355 r357  
    5050npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
    5151   
    52        
    53 def factorize(num):
    54     ''' Provide prime number factors for integer num
    55     Returns dictionary of prime factors (keys) & power for each (data)
    56     '''
    57     factors = {}
    58     orig = num
    59 
    60     # we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1
    61     i, sqi = 2, 4
    62     while sqi <= num:
    63         while not num%i:
    64             num /= i
    65             factors[i] = factors.get(i, 0) + 1
    66 
    67         sqi += 2*i + 1
    68         i += 1
    69 
    70     if num != 1 and num != orig:
    71         factors[num] = factors.get(num, 0) + 1
    72 
    73     if factors:
    74         return factors
    75     else:
    76         return {num:1}          #a prime number!
    77            
    78 def makeFFTsizeList(nmin=1,nmax=1023,thresh=15):
    79     ''' Provide list of optimal data sizes for FFT calculations
    80     Input:
    81         nmin: minimum data size >= 1
    82         nmax: maximum data size > nmin
    83         thresh: maximum prime factor allowed
    84     Returns:
    85         list of data sizes where the maximum prime factor is < thresh
    86     '''
    87     plist = []
    88     nmin = max(1,nmin)
    89     nmax = max(nmin+1,nmax)
    90     for p in range(nmin,nmax):
    91         if max(factorize(p).keys()) < thresh:
    92             plist.append(p)
    93     return plist
    94 
     52#GSASII pdf calculation routines
     53       
    9554def Transmission(Geometry,Abs,Diam):
    9655#Calculate sample transmission
    … …  
    176135#    return (Pola*npcosd(Azm)**2+(1.-Pola)*npsind(Azm)**2)*npcosd(Tth)**2+ \
    177136#        Pola*npsind(Azm)**2+(1.-Pola)*npcosd(Azm)**2
    178     return ((1.0-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+   \
     137    pola = ((1.0-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+   \
    179138        (1.0-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
     139    dpdPola = npsind(Azm)**2*npsind(Tth)**2
     140    return pola,dpdPola
    180141   
    181142def Oblique(ObCoeff,Tth):
    … …  
    407368        Amu += el['FormulaNo']*el['mu']
    408369       
    409 #GSASII peak fitting routine: Finger, Cox & Jephcoat model       
     370def CalcPDF(data,inst,xydata):
     371    auxPlot = []
     372    import copy
     373    import scipy.fftpack as ft
     374    #subtract backgrounds - if any
     375    xydata['IofQ'] = copy.deepcopy(xydata['Sample'])
     376    if data['Sample Bkg.']['Name']:
     377        xydata['IofQ'][1][1] += (xydata['Sample Bkg.'][1][1]+
     378            data['Sample Bkg.']['Add'])*data['Sample Bkg.']['Mult']
     379    if data['Container']['Name']:
     380        xycontainer = (xydata['Container'][1][1]+data['Container']['Add'])*data['Container']['Mult']
     381        if data['Container Bkg.']['Name']:
     382            xycontainer += (xydata['Container Bkg.'][1][1]+
     383                data['Container Bkg.']['Add'])*data['Container Bkg.']['Mult']
     384        xydata['IofQ'][1][1] += xycontainer
     385    #get element data & absorption coeff.
     386    ElList = data['ElList']
     387    Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
     388    #Apply angle dependent corrections
     389    Tth = xydata['Sample'][1][0]
     390    dt = (Tth[1]-Tth[0])
     391    xydata['IofQ'][1][1] /= Absorb(data['Geometry'],Abs,data['Diam'],Tth)
     392    xydata['IofQ'][1][1] /= Polarization(inst['Polariz.'],Tth,Azm=inst['Azimuth'])[0]
     393    if data['DetType'] == 'Image plate':
     394        xydata['IofQ'][1][1] *= Oblique(data['ObliqCoeff'],Tth)
     395    XY = xydata['IofQ'][1]   
     396    #convert to Q
     397    hc = 12.397639
     398    if 'Lam' in inst:
     399        wave = inst['Lam']
     400    else:
     401        wave = inst['Lam1']
     402    keV = hc/wave
     403    minQ = npT2q(Tth[0],wave)
     404    maxQ = npT2q(Tth[-1],wave)   
     405    Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
     406    dq = Qpoints[1]-Qpoints[0]
     407    XY[0] = npT2q(XY[0],wave)   
     408#    Qdata = np.nan_to_num(si.griddata(XY[0],XY[1],Qpoints,method='linear')) #only OK for scipy 0.9!
     409    T = si.interp1d(XY[0],XY[1],bounds_error=False,fill_value=0.0)      #OK for scipy 0.8
     410    Qdata = T(Qpoints)
     411   
     412    qLimits = data['QScaleLim']
     413    minQ = np.searchsorted(Qpoints,qLimits[0])
     414    maxQ = np.searchsorted(Qpoints,qLimits[1])
     415    newdata = []
     416    xydata['IofQ'][1][0] = Qpoints
     417    xydata['IofQ'][1][1] = Qdata
     418    for item in xydata['IofQ'][1]:
     419        newdata.append(item[:maxQ])
     420    xydata['IofQ'][1] = newdata
     421   
     422
     423    xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
     424    FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2)  #these are <f^2>,<f>^2,Cf
     425    Q = xydata['SofQ'][1][0]
     426    ruland = Ruland(data['Ruland'],wave,Q,CF)
     427#    auxPlot.append([Q,ruland,'Ruland'])     
     428    CF *= ruland
     429#    auxPlot.append([Q,CF,'CF-Corr'])
     430    scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
     431    xydata['SofQ'][1][1] *= scale
     432    xydata['SofQ'][1][1] -= CF
     433    xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
     434    scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
     435    xydata['SofQ'][1][1] *= scale
     436   
     437    xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
     438    xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
     439    if data['Lorch']:
     440        xydata['FofQ'][1][1] *= LorchWeight(Q)
     441   
     442    xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
     443    nR = len(xydata['GofR'][1][1])
     444    xydata['GofR'][1][1] = -dq*np.imag(ft.fft(xydata['FofQ'][1][1],4*nR)[:nR])
     445    xydata['GofR'][1][0] = 0.5*np.pi*np.linspace(0,nR,nR)/qLimits[1]
     446   
     447       
     448    return auxPlot
     449       
     450#GSASII peak fitting routines: Finger, Cox & Jephcoat model       
     451
     452def factorize(num):
     453    ''' Provide prime number factors for integer num
     454    Returns dictionary of prime factors (keys) & power for each (data)
     455    '''
     456    factors = {}
     457    orig = num
     458
     459    # we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1
     460    i, sqi = 2, 4
     461    while sqi <= num:
     462        while not num%i:
     463            num /= i
     464            factors[i] = factors.get(i, 0) + 1
     465
     466        sqi += 2*i + 1
     467        i += 1
     468
     469    if num != 1 and num != orig:
     470        factors[num] = factors.get(num, 0) + 1
     471
     472    if factors:
     473        return factors
     474    else:
     475        return {num:1}          #a prime number!
     476           
     477def makeFFTsizeList(nmin=1,nmax=1023,thresh=15):
     478    ''' Provide list of optimal data sizes for FFT calculations
     479    Input:
     480        nmin: minimum data size >= 1
     481        nmax: maximum data size > nmin
     482        thresh: maximum prime factor allowed
     483    Returns:
     484        list of data sizes where the maximum prime factor is < thresh
     485    '''
     486    plist = []
     487    nmin = max(1,nmin)
     488    nmax = max(nmin+1,nmax)
     489    for p in range(nmin,nmax):
     490        if max(factorize(p).keys()) < thresh:
     491            plist.append(p)
     492    return plist
    410493
    411494np.seterr(divide='ignore')
    … …  
    9521035    PeaksPrint(parmDict,sigDict,varyList)
    9531036   
    954 def CalcPDF(data,inst,xydata):
    955     auxPlot = []
    956     import copy
    957     import scipy.fftpack as ft
    958     #subtract backgrounds - if any
    959     xydata['IofQ'] = copy.deepcopy(xydata['Sample'])
    960     if data['Sample Bkg.']['Name']:
    961         xydata['IofQ'][1][1] += (xydata['Sample Bkg.'][1][1]+
    962             data['Sample Bkg.']['Add'])*data['Sample Bkg.']['Mult']
    963     if data['Container']['Name']:
    964         xycontainer = (xydata['Container'][1][1]+data['Container']['Add'])*data['Container']['Mult']
    965         if data['Container Bkg.']['Name']:
    966             xycontainer += (xydata['Container Bkg.'][1][1]+
    967                 data['Container Bkg.']['Add'])*data['Container Bkg.']['Mult']
    968         xydata['IofQ'][1][1] += xycontainer
    969     #get element data & absorption coeff.
    970     ElList = data['ElList']
    971     Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
    972     #Apply angle dependent corrections
    973     Tth = xydata['Sample'][1][0]
    974     dt = (Tth[1]-Tth[0])
    975     xydata['IofQ'][1][1] /= Absorb(data['Geometry'],Abs,data['Diam'],Tth)
    976     xydata['IofQ'][1][1] /= Polarization(inst['Polariz.'],Tth,Azm=inst['Azimuth'])
    977     if data['DetType'] == 'Image plate':
    978         xydata['IofQ'][1][1] *= Oblique(data['ObliqCoeff'],Tth)
    979     XY = xydata['IofQ'][1]   
    980     #convert to Q
    981     hc = 12.397639
    982     if 'Lam' in inst:
    983         wave = inst['Lam']
    984     else:
    985         wave = inst['Lam1']
    986     keV = hc/wave
    987     minQ = npT2q(Tth[0],wave)
    988     maxQ = npT2q(Tth[-1],wave)   
    989     Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
    990     dq = Qpoints[1]-Qpoints[0]
    991     XY[0] = npT2q(XY[0],wave)   
    992 #    Qdata = np.nan_to_num(si.griddata(XY[0],XY[1],Qpoints,method='linear')) #only OK for scipy 0.9!
    993     T = si.interp1d(XY[0],XY[1],bounds_error=False,fill_value=0.0)      #OK for scipy 0.8
    994     Qdata = T(Qpoints)
    995    
    996     qLimits = data['QScaleLim']
    997     minQ = np.searchsorted(Qpoints,qLimits[0])
    998     maxQ = np.searchsorted(Qpoints,qLimits[1])
    999     newdata = []
    1000     xydata['IofQ'][1][0] = Qpoints
    1001     xydata['IofQ'][1][1] = Qdata
    1002     for item in xydata['IofQ'][1]:
    1003         newdata.append(item[:maxQ])
    1004     xydata['IofQ'][1] = newdata
    1005    
    1006 
    1007     xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
    1008     FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2)  #these are <f^2>,<f>^2,Cf
    1009     Q = xydata['SofQ'][1][0]
    1010     ruland = Ruland(data['Ruland'],wave,Q,CF)
    1011 #    auxPlot.append([Q,ruland,'Ruland'])     
    1012     CF *= ruland
    1013 #    auxPlot.append([Q,CF,'CF-Corr'])
    1014     scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
    1015     xydata['SofQ'][1][1] *= scale
    1016     xydata['SofQ'][1][1] -= CF
    1017     xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
    1018     scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
    1019     xydata['SofQ'][1][1] *= scale
    1020    
    1021     xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
    1022     xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
    1023     if data['Lorch']:
    1024         xydata['FofQ'][1][1] *= LorchWeight(Q)
    1025    
    1026     xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
    1027     nR = len(xydata['GofR'][1][1])
    1028     xydata['GofR'][1][1] = -dq*np.imag(ft.fft(xydata['FofQ'][1][1],4*nR)[:nR])
    1029     xydata['GofR'][1][0] = 0.5*np.pi*np.linspace(0,nR,nR)/qLimits[1]
    1030    
    1031        
    1032     return auxPlot
    1033        
    10341037#testing data
    10351038NeedTestData = True

  • trunk/GSASIIstruct.py

    r350 r357  
    5757   
    5858def ShowControls(Controls):
    59     print ' Controls:'
    60     if Controls['bandWidth']:
    61         print ' Band approximated least squares matrix refinement, width: ',Controls['bandWidth']
    62     else:
    63         print ' Full matrix least squares refinement'
    64     print ' Maximum number of refinement cycles: ',Controls['Ncycles']
    65     print ' Minimum sum(shift/esd)^2 for convergence: ','%.2f'%(Controls['minSumShftESD'])
    66     print ' The Marquardt damping factor: ','%.2f'%(Controls['Marquardt'])
    67     print ' Maximum allowed atom shift: ','%.2f'%(Controls['maxShift']),'A'
    68     if Controls['restraintWeight']:
    69         print ' The weights of the restraint histograms will be modified to normalize their contribution','\n'
    70     else:
    71         print ' The restraint histogram weights are fixed','\n'
     59    print ' Least squares controls:'
     60    print ' Derivative type: ',Controls['deriv type']
     61    print ' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
    7262   
    7363def GetPhaseNames(GPXfile):
    … …  
    328318                if isum == jsum:
    329319                    G2mv.StoreEquivalence(varyJ,(varyI,))
     320def cellVary(pfx,SGData):
     321    if SGData['SGLaue'] in ['-1',]:
     322        return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
     323    elif SGData['SGLaue'] in ['2/m',]:
     324        if SGData['SGUniq'] == 'a':
     325            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
     326        elif SGData['SGUniq'] == 'b':
     327            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
     328        else:
     329            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
     330    elif SGData['SGLaue'] in ['mmm',]:
     331        return [pfx+'A0',pfx+'A1',pfx+'A2']
     332    elif SGData['SGLaue'] in ['4/m','4/mmm']:
     333        return [pfx+'A0',pfx+'A2']
     334    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
     335        return [pfx+'A0',pfx+'A2']
     336    elif SGData['SGLaue'] in ['3R', '3mR']:
     337        return [pfx+'A0',pfx+'A3']                       
     338    elif SGData['SGLaue'] in ['m3m','m3']:
     339        return [pfx+'A0']
     340   
    330341       
    331342def GetPhaseData(PhaseData):
    332343   
    333     def cellVary(pfx,SGData):
    334         if SGData['SGLaue'] in ['-1',]:
    335             return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
    336         elif SGData['SGLaue'] in ['2/m',]:
    337             if SGData['SGUniq'] == 'a':
    338                 return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
    339             elif SGData['SGUniq'] == 'b':
    340                 return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
    341             else:
    342                 return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
    343         elif SGData['SGLaue'] in ['mmm',]:
    344             return [pfx+'A0',pfx+'A1',pfx+'A2']
    345         elif SGData['SGLaue'] in ['4/m','4/mmm']:
    346             return [pfx+'A0',pfx+'A2']
    347         elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
    348             return [pfx+'A0',pfx+'A2']
    349         elif SGData['SGLaue'] in ['3R', '3mR']:
    350             return [pfx+'A0',pfx+'A3']                       
    351         elif SGData['SGLaue'] in ['m3m','m3']:
    352             return [pfx+'A0']
    353        
    354344       
    355345    print ' Phases:'
    … …  
    498488            for i,refl in enumerate(pawleyRef):
    499489                key = pfx+'PWLref:'+str(i)
    500                 refl[6] = parmDict[key]
     490                refl[6] = abs(parmDict[key])        #suppress negative Fsq
    501491                if key in sigDict:
    502492                    refl[7] = sigDict[key]
    … …  
    692682                ptlbls += '%12s' % (name)
    693683                ptstr += '%12.6f' % (hapData[4][i])
    694                 sigstr += '%12.6f' % (mustrainSig[1][i])
     684                if mustrainSig[1][i]:
     685                    sigstr += '%12.6f' % (mustrainSig[1][i])
     686                else:
     687                    sigstr += 12*' '
    695688            print ptlbls
    696689            print ptstr
    … …  
    779772    def GetSampleParms(hId,Sample):
    780773        sampVary = []
    781         pfx = ':'+str(hId)+':'       
    782         sampDict = {pfx+'Gonio. radius':Sample['Gonio. radius']}
     774        hfx = ':'+str(hId)+':'       
     775        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius']}
    783776        Type = Sample['Type']
    784777        if 'Bragg' in Type:             #Bragg-Brentano
    785778            for item in ['Scale','Shift','Transparency']:       #surface roughness?, diffuse scattering?
    786                 sampDict[pfx+item] = Sample[item][0]
     779                sampDict[hfx+item] = Sample[item][0]
    787780                if Sample[item][1]:
    788                     sampVary.append(pfx+item)
     781                    sampVary.append(hfx+item)
    789782        elif 'Debye' in Type:        #Debye-Scherrer
    790783            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
    791                 sampDict[pfx+item] = Sample[item][0]
     784                sampDict[hfx+item] = Sample[item][0]
    792785                if Sample[item][1]:
    793                     sampVary.append(pfx+item)
     786                    sampVary.append(hfx+item)
    794787        return Type,sampDict,sampVary
    795788       
    … …  
    883876   
    884877    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
    885         backVals = Background[3:]
    886         backSig = [0 for i in range(len(Background[3:]))]
    887         for i in range(len(backVals)):
    888             backVals[i] = parmDict[pfx+'Back:'+str(i)]
     878        lenBack = len(Background[3:])
     879        backSig = [0 for i in range(lenBack)]
     880        for i in range(lenBack):
     881            Background[3+i] = parmDict[pfx+'Back:'+str(i)]
    889882            if Background[1]:
    890883                backSig[i] = sigDict[pfx+'Back:'+str(i)]
    … …  
    972965            hId = Histogram['hId']
    973966            pfx = ':'+str(hId)+':'
    974            
    975967            Background = Histogram['Background'][0]
    976968            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
    … …  
    10791071   
    10801072    def GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse):
     1073        costh = cosd(refl[5]/2.)
    10811074        if calcControls[phfx+'SizeType'] == 'isotropic':
    1082             gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:0']*cosd(refl[5]/2.))
     1075            gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:0']*costh)
    10831076        elif calcControls[phfx+'SizeType'] == 'uniaxial':
    10841077            H = np.array(refl[:3])
    10851078            P = np.array(calcControls[phfx+'SizeAxis'])
    10861079            cosP,sinP = G2lat.CosSinAngle(H,P,G)
    1087             gam = (1.8*wave/np.pi)*parmDict[phfx+'Size:0']*parmDict[phfx+'Size:1']
    1088             gam /= np.sqrt((cosP*parmDict[phfx+'Size:1'])**2+(sinP*parmDict[phfx+'Size:0'])**2)*cosd(refl[5]/2.)
     1080            gam = (1.8*wave/np.pi)/parmDict[phfx+'Size:0']*parmDict[phfx+'Size:1']
     1081            gam *= np.sqrt((cosP*parmDict[phfx+'Size:1'])**2+(sinP*parmDict[phfx+'Size:0'])**2)/costh
    10891082        else:           #ellipsoidal crystallites
    10901083            H = np.array(refl[:3])
    1091             gam += 1.8*wave/(np.pi*cosd(refl[5])*np.inner(H,np.inner(sizeEllipse,H)))           
     1084            gam += 1.8*wave/(np.pi*costh*np.inner(H,np.inner(sizeEllipse,H)))           
    10921085        if calcControls[phfx+'MustrainType'] == 'isotropic':
    10931086            gam += 0.018*parmDict[phfx+'Mustrain:0']*tand(refl[5]/2.)/np.pi
    … …  
    11091102    def GetIntensityCorr(refl,phfx,hfx,calcControls,parmDict):
    11101103        Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
    1111         Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
     1104        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
    11121105       
    11131106        return Icorr
    … …  
    11881181                elif not iBeg-iFin:     #peak above high limit - done
    11891182                    return yc,yb
    1190                 yc[iBeg:iFin] += G2pwd.getFCJVoigt3(refl[5],Icorr*refl[8],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
     1183                yc[iBeg:iFin] += Icorr*refl[8]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
    11911184                if Ka2:
    11921185                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
    … …  
    11941187                    iBeg = np.searchsorted(x,pos2-fmin)
    11951188                    iFin = np.searchsorted(x,pos2+fmax)
    1196                     yc[iBeg:iFin] += G2pwd.getFCJVoigt3(pos2,Icorr*refl[8]*kRatio,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
     1189                    yc[iBeg:iFin] += Icorr*refl[8]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
    11971190            else:
    11981191                raise ValueError
    11991192    return yc,yb   
    12001193           
    1201 
    1202    
     1194def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
     1195   
     1196    def GetSampleGamDerv(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse):
     1197        gamDict = {}
     1198        costh = cosd(refl[5]/2.)
     1199        tanth = tand(refl[5]/2.)
     1200        if calcControls[phfx+'SizeType'] == 'isotropic':
     1201            gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:0']*costh)
     1202            gamDict[phfx+'Size:0'] = -gam/parmDict[phfx+'Size:0']
     1203        elif calcControls[phfx+'SizeType'] == 'uniaxial':
     1204            H = np.array(refl[:3])
     1205            P = np.array(calcControls[phfx+'SizeAxis'])
     1206            cosP,sinP = G2lat.CosSinAngle(H,P,G)
     1207            Si = parmDict[phfx+'Size:0']
     1208            Sa = parmDict[phfx+'Size:1']
     1209            gami = (1.8*wave/np.pi)/(Si*Sa)
     1210            sqtrm = np.sqrt((cosP*Sa)**2+(sinP*Si)**2)
     1211            gam = gami*sqtrm/costh           
     1212            gamDict[phfx+'Size:0'] = gami*Si*sinP**2/(sqtrm*costh)-gam/Si
     1213            gamDict[phfx+'Size:1'] = gami*Sa*cosP**2/(sqtrm*costh)-gam/Sa         
     1214        else:           #ellipsoidal crystallites - do numerically?
     1215            H = np.array(refl[:3])
     1216            gam = 1.8*wave/(np.pi*costh*np.inner(H,np.inner(sizeEllipse,H)))
     1217                       
     1218        if calcControls[phfx+'MustrainType'] == 'isotropic':
     1219            gamDict[phfx+'Mustrain:0'] =  0.018*tanth/np.pi           
     1220        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
     1221            H = np.array(refl[:3])
     1222            P = np.array(calcControls[phfx+'MustrainAxis'])
     1223            cosP,sinP = G2lat.CosSinAngle(H,P,G)
     1224            Si = parmDict[phfx+'Mustrain:0']
     1225            Sa = parmDict[phfx+'Mustrain:1']
     1226            gami = 0.018*Si*Sa*tanth/np.pi
     1227            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
     1228            gam = gami/sqtrm
     1229            gamDict[phfx+'Mustrain:0'] = gam/Si-gami*Si*cosP**2/sqtrm**3
     1230            gamDict[phfx+'Mustrain:1'] = gam/Sa-gami*Sa*sinP**2/sqtrm**3
     1231        else:       #generalized - P.W. Stephens model
     1232            pwrs = calcControls[phfx+'MuPwrs']
     1233            const = 0.018*refl[4]**2*tanth
     1234            for i,pwr in enumerate(pwrs):
     1235                gamDict[phfx+'Mustrain:'+str(i)] = const*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
     1236        return gamDict
     1237       
     1238    def GetIntensityDerv(refl,phfx,hfx,calcControls,parmDict):
     1239        Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
     1240        pola,dpdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
     1241        dIdpola = Icorr*dpdPola
     1242        Icorr *= pola
     1243        dIdsh = Icorr/parmDict[hfx+'Scale']
     1244        dIdsp = Icorr/parmDict[phfx+'Scale']
     1245       
     1246        return Icorr,dIdsh,dIdsp,dIdpola
     1247       
     1248    def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
     1249        dpr = 180./np.pi
     1250        h,k,l = refl[:3]
     1251        dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
     1252        dst = np.sqrt(dstsq)
     1253        pos = refl[5]
     1254        const = dpr/np.sqrt(1.0-wave*dst/4.0)
     1255        dpdw = const*dst
     1256        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
     1257        dpdA *= const*wave/(2.0*dst)
     1258        dpdZ = 1.0
     1259        const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
     1260        if 'Bragg' in calcControls[hfx+'instType']:
     1261            dpdSh = -4.*const*cosd(pos/2.0)
     1262            dpdTr = -1.e-7*const*sind(pos)
     1263            return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
     1264        else:               #Debye-Scherrer - simple but maybe not right
     1265            dpdXd = -const*cosd(pos)
     1266            dpdYd = -const*sind(pos)
     1267            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
     1268           
     1269    def cellVaryDerv(pfx,SGData,dpdA):
     1270        if SGData['SGLaue'] in ['-1',]:
     1271            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
     1272                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
     1273        elif SGData['SGLaue'] in ['2/m',]:
     1274            if SGData['SGUniq'] == 'a':
     1275                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
     1276            elif SGData['SGUniq'] == 'b':
     1277                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
     1278            else:
     1279                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
     1280        elif SGData['SGLaue'] in ['mmm',]:
     1281            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
     1282        elif SGData['SGLaue'] in ['4/m','4/mmm']:
     1283            return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]]
     1284        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
     1285            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]]
     1286        elif SGData['SGLaue'] in ['3R', '3mR']:
     1287            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
     1288        elif SGData['SGLaue'] in ['m3m','m3']:
     1289            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]]
     1290   
     1291    lenX = len(x)               
     1292    hId = Histogram['hId']
     1293    hfx = ':%d:'%(hId)
     1294    bakType = calcControls[hfx+'bakType']
     1295    dMdv = np.zeros(shape=(len(varylist),len(x)))
     1296    if hfx+'Back:0' in varylist:
     1297        dMdb = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
     1298        bBpos =varylist.index(hfx+'Back:0')
     1299        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
     1300       
     1301    if 'C' in calcControls[hfx+'histType']:   
     1302        dx = x[1]-x[0]
     1303        shl = max(parmDict[hfx+'SH/L'],0.002)
     1304        Ka2 = False
     1305        if hfx+'Lam1' in parmDict.keys():
     1306            wave = parmDict[hfx+'Lam1']
     1307            Ka2 = True
     1308            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
     1309            kRatio = parmDict[hfx+'I(L2)/I(L1)']
     1310        else:
     1311            wave = parmDict[hfx+'Lam']
     1312    else:
     1313        print 'TOF Undefined at present'
     1314        raise ValueError
     1315    for phase in Histogram['Reflection Lists']:
     1316        refList = Histogram['Reflection Lists'][phase]
     1317        Phase = Phases[phase]
     1318        SGData = Phase['General']['SGData']
     1319        pId = Phase['pId']
     1320        pfx = '%d::'%(pId)
     1321        phfx = '%d:%d:'%(pId,hId)
     1322        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
     1323        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
     1324        sizeEllipse = []
     1325        if calcControls[phfx+'SizeType'] == 'ellipsoidal':
     1326            sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)])
     1327        for iref,refl in enumerate(refList):
     1328            if 'C' in calcControls[hfx+'histType']:
     1329                Icorr,dIdsh,dIdsp,dIdpola = GetIntensityDerv(refl,phfx,hfx,calcControls,parmDict)
     1330                hkl = refl[:3]
     1331                pos = refl[5]
     1332                tanth = tand(pos/2.0)
     1333                costh = cosd(pos/2.0)
     1334                dsdU = tanth**2
     1335                dsdV = tanth
     1336                dsdW = 1.0
     1337                dgdX = 1.0/costh
     1338                dgdY = tanth
     1339                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
     1340                iBeg = np.searchsorted(x,refl[5]-fmin)
     1341                iFin = np.searchsorted(x,refl[5]+fmax)
     1342                dMdpk = np.zeros(shape=(6,len(x)))
     1343                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
     1344                for i in range(1,5):
     1345                    dMdpk[i][iBeg:iFin] += 100.*dx*Icorr*refl[8]*dMdipk[i]
     1346                dMdpk[0][iBeg:iFin] += 100.*dx*Icorr*dMdipk[0]
     1347                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
     1348                if Ka2:
     1349                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
     1350                    kdelt = int((pos2-refl[5])/dx)               
     1351                    iBeg = min(lenX,iBeg+kdelt)
     1352                    iFin = min(lenX,iFin+kdelt)
     1353                    if iBeg-iFin:
     1354                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,xdata[iBeg:iFin])
     1355                        for i in range(1,5):
     1356                            dMdpk[i][iBeg:iFin] += 100.*dx*Icorr*refl[8]*kRatio*dMdipk2[i]
     1357                        dMdpk[0][iBeg:iFin] += 100.*dx*kRatio*dMdipk2[0]
     1358                        dMdpk[5][iBeg:iFin] += 100.*dx*dMdipk2[0]
     1359                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*Icorr*refl[8]}
     1360                try:
     1361                    idx = varylist.index(pfx+'PWLref:'+str(iref))
     1362                    dMdv[idx] = dervDict['int']
     1363                except ValueError:
     1364                    pass
     1365                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
     1366                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
     1367                    hfx+'U':[dsdU,'sig'],hfx+'V':[dsdV,'sig'],hfx+'W':[dsdW,'sig'],
     1368                    hfx+'X':[dgdX,'gam'],hfx+'Y':[dgdY,'gam'],hfx+'SH/L':[1.0,'shl'],
     1369                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
     1370                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
     1371                    hfx+'DisplaceY':[dpdY,'pos'],}
     1372                for name in names:
     1373                    if name in varylist:
     1374                        item = names[name]
     1375                        dMdv[varylist.index(name)] += item[0]*dervDict[item[1]]
     1376                cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
     1377                for name,dpdA in cellDervNames:
     1378                    if name in varylist:
     1379                        dMdv[varylist.index(name)] += dpdA*dervDict['pos']
     1380                gamDict = GetSampleGamDerv(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse)
     1381                for name in gamDict:
     1382                    if name in varylist:
     1383                        dMdv[varylist.index(name)] += gamDict[name]*dervDict['gam']                       
     1384            else:
     1385                raise ValueError
     1386           
     1387    return dMdv   
     1388                   
    12031389def Refine(GPXfile):
     1390   
     1391    def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
     1392        parmdict.update(zip(varylist,values))
     1393        G2mv.Dict2Map(parmDict)
     1394        Histograms,Phases = HistoPhases
     1395        dMdv = np.empty(0)
     1396        for histogram in Histograms:
     1397            if 'PWDR' in histogram[:4]:
     1398                Histogram = Histograms[histogram]
     1399                hId = Histogram['hId']
     1400                hfx = ':%d:'%(hId)
     1401                Limits = calcControls[hfx+'Limits']
     1402                x,y,w,yc,yb,yd = Histogram['Data']
     1403                xB = np.searchsorted(x,Limits[0])
     1404                xF = np.searchsorted(x,Limits[1])
     1405                dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
     1406                    varylist,Histogram,Phases,calcControls,pawleyLookup)
     1407                if len(dMdv):
     1408                    dMdv = np.concatenate((dMdv,dMdvh))
     1409                else:
     1410                    dMdv = dMdvh
     1411        return dMdv
    12041412   
    12051413    def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
    … …  
    12291437                yc[xB:xF] *= parmDict[hfx+'Scale']
    12301438                yc[xB:xF] += yb[xB:xF]
    1231                 yd[xB:xF] = y[xB:xF]-yc[xB:xF]
     1439                yd[xB:xF] = yc[xB:xF]-y[xB:xF]          #yc-yo then all dydv have no '-' needed
    12321440                Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
    12331441                M = np.concatenate((M,np.sqrt(w[xB:xF])*(yd[xB:xF])))
    … …  
    12451453    parmDict = {}
    12461454    calcControls = {}   
    1247 #    Controls = GetControls(GPXfile)
    1248 #    ShowControls(Controls)           
     1455    Controls = GetControls(GPXfile)
     1456    ShowControls(Controls)           
    12491457    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
    12501458    if not Phases:
    … …  
    12781486        Size = dlg.GetSize()
    12791487        dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
     1488        Ftol = Controls['min dM/M']
    12801489        try:
    1281             result = so.leastsq(errRefine,values,full_output=True,ftol=0.0001,epsfcn=1.e-8,
    1282                 args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
     1490            if Controls['deriv type'] == 'analytic':
     1491                result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,ftol=Ftol,col_deriv=True,
     1492                    args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
     1493                ncyc = int(result[2]['nfev']/2)               
     1494            else:           #'numeric'
     1495                result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,
     1496                    args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
     1497                ncyc = int(result[2]['nfev']/len(varyList))
    12831498        finally:
    12841499            dlg.Destroy()
    1285         runtime = time.time()-begin   
     1500        runtime = time.time()-begin
    12861501        chisq = np.sum(result[2]['fvec']**2)
    1287         ncyc = int(result[2]['nfev']/len(varyList))
    12881502        Values2Dict(parmDict, varyList, result[0])
    12891503        G2mv.Dict2Map(parmDict)
    … …  
    13161530    SetHistogramData(parmDict,sigDict,Histograms)
    13171531    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases)
     1532#for testing purposes!!!
     1533#    import cPickle
     1534#    file = open('structTestdata.dat','wb')
     1535#    cPickle.dump(parmDict,file,1)
     1536#    cPickle.dump(varyList,file,1)
     1537#    for histogram in Histograms:
     1538#        if 'PWDR' in histogram[:4]:
     1539#            Histogram = Histograms[histogram]
     1540#    cPickle.dump(Histogram,file,1)
     1541#    cPickle.dump(Phases,file,1)
     1542#    cPickle.dump(calcControls,file,1)
     1543#    cPickle.dump(pawleyLookup,file,1)
     1544#    file.close()
    13181545
    13191546def main():
Note: See TracChangeset for help on using the changeset viewer.