Changeset 357 for trunk/GSASIIpwd.py


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

File:
1 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
Note: See TracChangeset for help on using the changeset viewer.