Changeset 1595


Ignore:
Timestamp:
Dec 4, 2014 4:07:45 PM (9 years ago)
Author:
vondreele
Message:

change ValEsd? to round at 3.16228 =sqrt(10); same as old GSAS
add mod. Vector names to list
begin Pawley refinement for super structures - use an offset to allow shift in refl values due to extra index.
allow cell refinement from index peaks & entered cell parameters

Location:
trunk
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIImath.py

    r1581 r1595  
    13901390    '''
    13911391    # Note: this routine is Python 3 compatible -- I think
     1392    cutoff = 3.16228    #=(sqrt(10); same as old GSAS   was 1.95
    13921393    if math.isnan(value): # invalid value, bail out
    13931394        return '?'
    … …  
    13971398    elif esd != 0:
    13981399        # transform the esd to a one or two digit integer
    1399         l = math.log10(abs(esd)) % 1
    1400         # cut off of 19 1.9==>(19) but 1.95==>(2) N.B. log10(1.95) = 0.2900...
    1401         if l < 0.290034611362518: l+= 1.       
     1400        l = math.log10(abs(esd)) % 1.
     1401        if l < math.log10(cutoff): l+= 1.       
    14021402        intesd = int(round(10**l)) # esd as integer
    14031403        # determine the number of digits offset for the esd

  • trunk/GSASIIobj.py

    r1496 r1595  
    12531253        'A([0-5])' : 'Reciprocal metric tensor component \\1',
    12541254        'Vol' : 'Unit cell volume',
     1255        'mV([0-2])' : 'Modulation vector component \\1',
    12551256        # Atom vars (p::<var>:a)
    12561257        'dA([xyz])$' : 'change to atomic coordinate, \\1',

  • trunk/GSASIIphsGUI.py

    r1594 r1595  
    50565056                for h,k,l,m,d in HKLd:
    50575057                    ext,mul = G2spc.GenHKLf([h,k,l],SGData)[:2]
    5058                     if not ext:
     5058                    if m or not ext:
    50595059                        mul *= 2        #for powder multiplicity
    50605060                        PawleyPeaks.append([h,k,l,m,mul,d,False,100.0,1.0])
    … …  
    50835083            return
    50845084        Vst = 1.0/data['General']['Cell'][7]     #Get volume
     5085        generalData = data['General']
     5086        im = 0
     5087        if generalData['Type'] in ['modulated','magnetic',]:
     5088            im = 1
    50855089        HistoNames = filter(lambda a:Histograms[a]['Use']==True,Histograms.keys())
    50865090        PatternId = G2gd.GetPatternTreeItemId(G2frame,G2frame.root,HistoNames[0])
    … …  
    50965100        try:
    50975101            for ref in Refs:
    5098                 pos = 2.0*asind(wave/(2.0*ref[4]))
     5102                pos = 2.0*asind(wave/(2.0*ref[4+im]))
    50995103                if 'Bragg' in Sample['Type']:
    51005104                    pos -= const*(4.*Sample['Shift'][0]*cosd(pos/2.0)+ \
    … …  
    51095113                    # we multiply the observed peak height by sqrt(8 ln 2)/(FWHM*sqrt(pi)) to determine the value of Icorr*F^2
    51105114                    # then divide by Icorr to get F^2.
    5111                     ref[6] = (xdata[1][indx]-xdata[4][indx])*gconst/(FWHM*np.sqrt(np.pi))  #Area of Gaussian is height * FWHM * sqrt(pi)
     5115                    ref[6+im] = (xdata[1][indx]-xdata[4][indx])*gconst/(FWHM*np.sqrt(np.pi))  #Area of Gaussian is height * FWHM * sqrt(pi)
    51125116                    Lorenz = 1./(2.*sind(xdata[0][indx]/2.)**2*cosd(xdata[0][indx]/2.))           #Lorentz correction
    51135117                    pola = 1.0
    … …  
    51175121                        pola = 1.0
    51185122                    # Include histo scale and volume in calculation
    5119                     ref[6] /= (Sample['Scale'][0] * Vst * Lorenz * pola * ref[3])
     5123                    ref[6+im] /= (Sample['Scale'][0] * Vst * Lorenz * pola * ref[3+im])
    51205124                except IndexError:
    51215125                    pass

  • trunk/GSASIIpwdGUI.py

    r1594 r1595  
    23012301    def OnSpcSel(event):
    23022302        controls[13] = spcSel.GetString(spcSel.GetSelection())
     2303        G2frame.dataFrame.RefineCell.Enable(True)
    23032304       
    23042305    def SetCellValue(Obj,ObjId,value):
    … …  
    29402941        2*[wg.GRID_VALUE_FLOAT+':10,2',]+[wg.GRID_VALUE_FLOAT+':10,3',]+ \
    29412942        [wg.GRID_VALUE_FLOAT+':10,3',]
    2942     superLabels = ['M1','M2','M3']
    29432943    if HKLF:
    29442944        colLabels = ['H','K','L','mul','d','Fosq','sig','Fcsq','FoTsq','FcTsq','phase','ExtC',]
    … …  
    29472947            Types += 2*[wg.GRID_VALUE_FLOAT+':10,3',]
    29482948        if Super:
    2949             for i in range(Super):
    2950                 colLabels.insert(3+i,superLabels[i])
     2949            colLabels.insert(3,'M')
    29512950    else:
    29522951        if 'C' in Inst['Type'][0]:
    … …  
    29572956            Types += 7*[wg.GRID_VALUE_FLOAT+':10,3',]
    29582957        if Super:
    2959             for i in range(Super):
    2960                 colLabels.insert(3+i,superLabels[i])
     2958            colLabels.insert(3,'M')
    29612959           
    29622960    G2frame.PeakTable = G2gd.Table(refs,rowLabels=rowLabels,colLabels=colLabels,types=Types)

  • trunk/GSASIIstrIO.py

    r1594 r1595  
    973973        if cell[0]:
    974974            phaseVary += cellVary(pfx,SGData)
     975        SSGtext = []    #no superstructure
     976        im = 0
    975977        if General['Type'] in ['modulated','magnetic']:
     978            im = 1      #refl offset
    976979            Vec,vRef,maxH = General['SuperVec']
    977980            phaseDict.update({pfx+'mV0':Vec[0],pfx+'mV1':Vec[1],pfx+'mV2':Vec[2]})
    … …  
    10801083                PrintBLtable(BLtable)
    10811084                print >>pFile,''
    1082                 for line in SGtext: print >>pFile,line
    1083                 if len(SGtable):
    1084                     for item in SGtable:
    1085                         line = ' %s '%(item)
    1086                         print >>pFile,line   
     1085                if len(SSGtext):    #if superstructure
     1086                    for line in SSGtext: print >>pFile,line
     1087                    if len(SSGtable):
     1088                        for item in SSGtable:
     1089                            line = ' %s '%(item)
     1090                            print >>pFile,line   
     1091                    else:
     1092                        print >>pFile,' ( 1)    %s'%(SSGtable[0])
    10871093                else:
    1088                     print >>pFile,' ( 1)    %s'%(SGtable[0])
     1094                    for line in SGtext: print >>pFile,line
     1095                    if len(SGtable):
     1096                        for item in SGtable:
     1097                            line = ' %s '%(item)
     1098                            print >>pFile,line   
     1099                    else:
     1100                        print >>pFile,' ( 1)    %s'%(SGtable[0])
    10891101                PrintRBObjects(resRBData,vecRBData)
    10901102                PrintAtoms(General,Atoms)
    1091                 print >>pFile,'\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
    1092                     ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
    1093                     '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
     1103                print >>pFile,'\n Unit cell: a = %.5f'%(cell[1]),' b = %.5f'%(cell[2]),' c = %.5f'%(cell[3]), \
     1104                    ' alpha = %.3f'%(cell[4]),' beta = %.3f'%(cell[5]),' gamma = %.3f'%(cell[6]), \
     1105                    ' volume = %.3f'%(cell[7]),' Refine?',cell[0]
     1106                if len(SSGtext):    #if superstructure
     1107                    print >>pFile,'\n Modulation vector: mV0 = %.4f'%(Vec[0]),' mV1 = %.4f'%(Vec[1]),   \
     1108                        ' mV2 = %.4f'%(Vec[2]),' max mod. index = %d'%(maxH),' Refine?',vRef
    10941109                PrintTexture(textureData)
    10951110                if name in RestraintDict:
    … …  
    10981113                   
    10991114        elif PawleyRef:
     1115            if Print:
     1116                print >>pFile,'\n Phase name: ',General['Name']
     1117                print >>pFile,135*'-'
     1118                print >>pFile,''
     1119                if len(SSGtext):    #if superstructure
     1120                    for line in SSGtext: print >>pFile,line
     1121                    if len(SSGtable):
     1122                        for item in SSGtable:
     1123                            line = ' %s '%(item)
     1124                            print >>pFile,line   
     1125                    else:
     1126                        print >>pFile,' ( 1)    %s'%(SSGtable[0])
     1127                else:
     1128                    for line in SGtext: print >>pFile,line
     1129                    if len(SGtable):
     1130                        for item in SGtable:
     1131                            line = ' %s '%(item)
     1132                            print >>pFile,line   
     1133                    else:
     1134                        print >>pFile,' ( 1)    %s'%(SGtable[0])
     1135                print >>pFile,'\n Unit cell: a = %.5f'%(cell[1]),' b = %.5f'%(cell[2]),' c = %.5f'%(cell[3]), \
     1136                    ' alpha = %.3f'%(cell[4]),' beta = %.3f'%(cell[5]),' gamma = %.3f'%(cell[6]), \
     1137                    ' volume = %.3f'%(cell[7]),' Refine?',cell[0]
     1138                if len(SSGtext):    #if superstructure
     1139                    print >>pFile,'\n Modulation vector: mV0 = %.4f'%(Vec[0]),' mV1 = %.4f'%(Vec[1]),   \
     1140                        ' mV2 = %.4f'%(Vec[2]),' max mod. index = %d'%(maxH),' Refine?',vRef
    11001141            pawleyVary = []
    11011142            for i,refl in enumerate(PawleyRef):
    1102                 phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
    1103                 pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
    1104                 if refl[5]:
     1143                phaseDict[pfx+'PWLref:'+str(i)] = refl[6+im]
     1144                if im:
     1145                    pawleyLookup[pfx+'%d,%d,%d,%d'%(refl[0],refl[1],refl[2],refl[3])] = i
     1146                else:
     1147                    pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
     1148                if refl[5+im]:
    11051149                    pawleyVary.append(pfx+'PWLref:'+str(i))
    11061150            GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary)      #does G2mv.StoreEquivalence
    … …  
    15521596            print >>pFile,ptstr
    15531597            print >>pFile,sigstr
     1598        ik = 6  #for Pawley stuff below
     1599        if General['Type'] in ['modulated','magnetic']:
     1600            ik = 7
     1601            Vec,vRef,maxH = General['SuperVec']
     1602            if vRef:
     1603                print >>pFile,' New modulation vector:'
     1604                namstr = '  names :'
     1605                ptstr =  '  values:'
     1606                sigstr = '  esds  :'
     1607                for var in ['mV0','mV1','mV2']:
     1608                    namstr += '%12s'%(pfx+var)
     1609                    ptstr += '%12.4f'%(parmDict[pfx+var])
     1610                    if pfx+var in sigDict:
     1611                        sigstr += '%12.4f'%(sigDict[pfx+var])
     1612                    else:
     1613                        sigstr += 12*' '
     1614                print >>pFile,namstr
     1615                print >>pFile,ptstr
     1616                print >>pFile,sigstr
    15541617           
    15551618        General['Mass'] = 0.
    … …  
    15581621            for i,refl in enumerate(pawleyRef):
    15591622                key = pfx+'PWLref:'+str(i)
    1560                 refl[6] = parmDict[key]
     1623                refl[ik] = parmDict[key]
    15611624                if key in sigDict:
    1562                     refl[7] = sigDict[key]
     1625                    refl[ik+1] = sigDict[key]
    15631626                else:
    1564                     refl[7] = 0
     1627                    refl[ik+1] = 0
    15651628        else:
    15661629            VRBIds = RBIds['Vector']
    … …  
    17431806        cell = Phases[phase]['General']['Cell'][1:7]
    17441807        A = G2lat.cell2A(cell)
     1808        if Phases[phase]['General']['Type'] in ['modulated','magnetic']:
     1809            SSGData = Phases[phase]['General']['SSGData']
     1810            Vec,x,maxH = Phases[phase]['General']['SuperVec']
    17451811        pId = Phases[phase]['pId']
    17461812        histoList = HistoPhase.keys()
    … …  
    18401906                    PrintHStrain(hapData['HStrain'],SGData)
    18411907                    if hapData['Babinet']['BabA'][0]:
    1842                         PrintBabinet(hapData['Babinet'])
    1843                 HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))  #+DIJS
     1908                        PrintBabinet(hapData['Babinet'])                       
     1909                if Phases[phase]['General']['Type'] in ['modulated','magnetic']:
     1910                    HKLd = G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A)
     1911                else:
     1912                    HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
    18441913                if resetRefList:
    18451914                    refList = []
    18461915                    Uniq = []
    18471916                    Phi = []
    1848                     for h,k,l,d in HKLd:
    1849                         ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
    1850                         mul *= 2      # for powder overlap of Friedel pairs
    1851                         if ext:
    1852                             continue
    1853                         if 'C' in inst['Type'][0]:
    1854                             pos = 2.0*asind(wave/(2.0*d))+Zero
    1855                             if limits[0] < pos < limits[1]:
    1856                                 refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,0.0, 0.0,0.0,1.0,1.0,1.0])
    1857                                 #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
    1858                                 Uniq.append(uniq)
    1859                                 Phi.append(phi)
    1860                         elif 'T' in inst['Type'][0]:
    1861                             pos = inst['difC'][1]*d+inst['difA'][1]*d**2+inst['difB'][1]/d+Zero
    1862                             if limits[0] < pos < limits[1]:
    1863                                 wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
    1864                                 refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,wave, 1.0,1.0,1.0])
    1865                                 # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
    1866                                 Uniq.append(uniq)
    1867                                 Phi.append(phi)
     1917                    if Phases[phase]['General']['Type'] in ['modulated','magnetic']:
     1918                        for h,k,l,m,d in HKLd:
     1919                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)    #is this right for SS refl.??
     1920                            mul *= 2      # for powder overlap of Friedel pairs
     1921                            if m or not ext:
     1922                                if 'C' in inst['Type'][0]:
     1923                                    pos = G2lat.Dsp2pos(inst,d)
     1924                                    if limits[0] < pos < limits[1]:
     1925                                        refList.append([h,k,l,m,mul,d, pos,0.0,0.0,0.0,0.0, 0.0,0.0,1.0,1.0,1.0])
     1926                                        #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
     1927                                        Uniq.append(uniq)
     1928                                        Phi.append(phi)
     1929                                elif 'T' in inst['Type'][0]:
     1930                                    pos = G2lat.Dsp2pos(inst,d)
     1931                                    if limits[0] < pos < limits[1]:
     1932                                        wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
     1933                                        refList.append([h,k,l,m,mul,d, pos,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,wave, 1.0,1.0,1.0])
     1934                                        # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
     1935                                        Uniq.append(uniq)
     1936                                        Phi.append(phi)
     1937                    else:
     1938                        for h,k,l,d in HKLd:
     1939                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
     1940                            mul *= 2      # for powder overlap of Friedel pairs
     1941                            if ext:
     1942                                continue
     1943                            if 'C' in inst['Type'][0]:
     1944                                pos = G2lat.Dsp2pos(inst,d)
     1945                                if limits[0] < pos < limits[1]:
     1946                                    refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,0.0, 0.0,0.0,1.0,1.0,1.0])
     1947                                    #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
     1948                                    Uniq.append(uniq)
     1949                                    Phi.append(phi)
     1950                            elif 'T' in inst['Type'][0]:
     1951                                pos = G2lat.Dsp2pos(inst,d)
     1952                                if limits[0] < pos < limits[1]:
     1953                                    wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
     1954                                    refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,wave, 1.0,1.0,1.0])
     1955                                    # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
     1956                                    Uniq.append(uniq)
     1957                                    Phi.append(phi)
    18681958                    Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0]}
    18691959            elif 'HKLF' in histogram:

  • trunk/GSASIIstrMath.py

    r1591 r1595  
    798798    return dFdvDict
    799799   
    800 def SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varyList):
     800def SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varyList):
    801801    ''' Single crystal extinction function; returns extinction & derivative
    802802    '''
    … …  
    804804    dervDict = {}
    805805    if calcControls[phfx+'EType'] != 'None':
    806         SQ = 1/(4.*ref[4]**2)
     806        SQ = 1/(4.*ref[4+im]**2)
    807807        if 'C' in parmDict[hfx+'Type']:           
    808808            cos2T = 1.0-2.*SQ*parmDict[hfx+'Lam']**2           #cos(2theta)
    809809        else:   #'T'
    810             cos2T = 1.0-2.*SQ*ref[12]**2                       #cos(2theta)           
     810            cos2T = 1.0-2.*SQ*ref[12+im]**2                       #cos(2theta)           
    811811        if 'SXC' in parmDict[hfx+'Type']:
    812812            AV = 7.9406e5/parmDict[pfx+'Vol']**2
    813813            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
    814814            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
    815             PLZ = AV*P12*ref[7]*parmDict[hfx+'Lam']**2
     815            PLZ = AV*P12*ref[7+im]*parmDict[hfx+'Lam']**2
    816816        elif 'SNT' in parmDict[hfx+'Type']:
    817817            AV = 1.e7/parmDict[pfx+'Vol']**2
    818818            PL = SQ
    819             PLZ = AV*ref[7]*ref[12]**2
     819            PLZ = AV*ref[7+im]*ref[12+im]**2
    820820        elif 'SNC' in parmDict[hfx+'Type']:
    821821            AV = 1.e7/parmDict[pfx+'Vol']**2
    … …  
    829829                PLZ *= calcControls[phfx+'Tbar']
    830830            else: #'T'
    831                 PLZ *= ref[13]      #t-bar
     831                PLZ *= ref[13+im]      #t-bar
    832832        if 'Primary' in calcControls[phfx+'EType']:
    833833            PLZ *= 1.5
    … …  
    860860
    861861        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
    862             dervDict[phfx+'Ep'] = -ref[7]*PLZ*PF3
     862            dervDict[phfx+'Ep'] = -ref[7+im]*PLZ*PF3
    863863        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
    864             dervDict[phfx+'Es'] = -ref[7]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
     864            dervDict[phfx+'Es'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
    865865        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
    866             dervDict[phfx+'Eg'] = -ref[7]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
     866            dervDict[phfx+'Eg'] = -ref[7+im]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
    867867               
    868868    return 1./extCor,dervDict
    … …  
    908908    return newAtomDict
    909909   
    910 def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
     910def SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
    911911    'Spherical harmonics texture'
    912912    IFCoup = 'Bragg' in calcControls[hfx+'instType']
    … …  
    914914        tth = parmDict[hfx+'2-theta']
    915915    else:
    916         tth = refl[5]
     916        tth = refl[5+im]
    917917    odfCor = 1.0
    918918    H = refl[:3]
    … …  
    931931    return odfCor
    932932   
    933 def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
     933def SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict):
    934934    'Spherical harmonics texture derivatives'
    935935    if 'T' in calcControls[hfx+'histType']:
    936936        tth = parmDict[hfx+'2-theta']
    937937    else:
    938         tth = refl[5]
     938        tth = refl[5+im]
    939939    FORPI = 4.0*np.pi
    940940    IFCoup = 'Bragg' in calcControls[hfx+'instType']
    … …  
    960960    return odfCor,dFdODF,dFdSA
    961961   
    962 def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
     962def SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
    963963    'spherical harmonics preferred orientation (cylindrical symmetry only)'
    964964    if 'T' in calcControls[hfx+'histType']:
    965965        tth = parmDict[hfx+'2-theta']
    966966    else:
    967         tth = refl[5]
     967        tth = refl[5+im]
    968968    odfCor = 1.0
    969969    H = refl[:3]
    … …  
    985985    return np.squeeze(odfCor)
    986986   
    987 def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
     987def SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict):
    988988    'spherical harmonics preferred orientation derivatives (cylindrical symmetry only)'
    989989    if 'T' in calcControls[hfx+'histType']:
    990990        tth = parmDict[hfx+'2-theta']
    991991    else:
    992         tth = refl[5]
     992        tth = refl[5+im]
    993993    FORPI = 12.5663706143592
    994994    odfCor = 1.0
    … …  
    10131013    return odfCor,dFdODF
    10141014   
    1015 def GetPrefOri(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
     1015def GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
    10161016    'March-Dollase preferred orientation correction'
    10171017    POcorr = 1.0
    … …  
    10271027    return POcorr
    10281028   
    1029 def GetPrefOriDerv(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
     1029def GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict):
    10301030    'Needs a doc string'
    10311031    POcorr = 1.0
    … …  
    10451045    else:   #spherical harmonics
    10461046        if calcControls[phfx+'SHord']:
    1047             POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
     1047            POcorr,POderv = SHPOcalDerv(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
    10481048    return POcorr,POderv
    10491049   
    1050 def GetAbsorb(refl,hfx,calcControls,parmDict):
     1050def GetAbsorb(refl,im,hfx,calcControls,parmDict):
    10511051    'Needs a doc string'
    10521052    if 'Debye' in calcControls[hfx+'instType']:
    10531053        if 'T' in calcControls[hfx+'histType']:
    1054             return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14],parmDict[hfx+'2-theta'],0,0)
     1054            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],parmDict[hfx+'2-theta'],0,0)
    10551055        else:
    1056             return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
     1056            return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
    10571057    else:
    1058         return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5])
    1059    
    1060 def GetAbsorbDerv(refl,hfx,calcControls,parmDict):
     1058        return G2pwd.SurfaceRough(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im])
     1059   
     1060def GetAbsorbDerv(refl,im,hfx,calcControls,parmDict):
    10611061    'Needs a doc string'
    10621062    if 'Debye' in calcControls[hfx+'instType']:
    10631063        if 'T' in calcControls[hfx+'histType']:
    1064             return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14],parmDict[hfx+'2-theta'],0,0)
     1064            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption']*refl[14+im],parmDict[hfx+'2-theta'],0,0)
    10651065        else:
    1066             return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
     1066            return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5+im],0,0)
    10671067    else:
    1068         return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5]))
     1068        return np.array(G2pwd.SurfaceRoughDerv(parmDict[hfx+'SurfRoughA'],parmDict[hfx+'SurfRoughB'],refl[5+im]))
    10691069       
    1070 def GetPwdrExt(refl,pfx,phfx,hfx,calcControls,parmDict):
     1070def GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict):
    10711071    'Needs a doc string'
    10721072    coef = np.array([-0.5,0.25,-0.10416667,0.036458333,-0.0109375,2.8497409E-3])
    … …  
    10741074    if 'T' in calcControls[hfx+'histType']:
    10751075        sth2 = sind(parmDict[hfx+'2-theta']/2.)**2
    1076         wave = refl[14]
     1076        wave = refl[14+im]
    10771077    else:   #'C'W
    1078         sth2 = sind(refl[5]/2.)**2
     1078        sth2 = sind(refl[5+im]/2.)**2
    10791079        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
    10801080    c2th = 1.-2.0*sth2
    1081     flv2 = refl[9]*(wave/parmDict[pfx+'Vol'])**2
     1081    flv2 = refl[9+im]*(wave/parmDict[pfx+'Vol'])**2
    10821082    if 'X' in calcControls[hfx+'histType']:
    10831083        flv2 *= 0.079411*(1.0+c2th**2)/2.0
    … …  
    10951095    return exb*sth2+exl*(1.-sth2)
    10961096   
    1097 def GetPwdrExtDerv(refl,pfx,phfx,hfx,calcControls,parmDict):
     1097def GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict):
    10981098    'Needs a doc string'
    10991099    delt = 0.001
    11001100    parmDict[phfx+'Extinction'] += delt
    1101     plus = GetPwdrExt(refl,pfx,phfx,hfx,calcControls,parmDict)
     1101    plus = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
    11021102    parmDict[phfx+'Extinction'] -= 2.*delt
    1103     minus = GetPwdrExt(refl,pfx,phfx,hfx,calcControls,parmDict)
     1103    minus = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
    11041104    parmDict[phfx+'Extinction'] += delt
    11051105    return (plus-minus)/(2.*delt)   
    11061106   
    1107 def GetIntensityCorr(refl,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
     1107def GetIntensityCorr(refl,im,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
    11081108    'Needs a doc string'    #need powder extinction!
    1109     Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
     1109    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3+im]               #scale*multiplicity
    11101110    if 'X' in parmDict[hfx+'Type']:
    1111         Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
     1111        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])[0]
    11121112    POcorr = 1.0
    11131113    if pfx+'SHorder' in parmDict:                 #generalized spherical harmonics texture
    1114         POcorr = SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
     1114        POcorr = SHTXcal(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
    11151115    elif calcControls[phfx+'poType'] == 'MD':         #March-Dollase
    1116         POcorr = GetPrefOri(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
     1116        POcorr = GetPrefOri(uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)
    11171117    elif calcControls[phfx+'SHord']:                #cylindrical spherical harmonics
    1118         POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
     1118        POcorr = SHPOcal(refl,im,g,phfx,hfx,SGData,calcControls,parmDict)
    11191119    Icorr *= POcorr
    11201120    AbsCorr = 1.0
    1121     AbsCorr = GetAbsorb(refl,hfx,calcControls,parmDict)
     1121    AbsCorr = GetAbsorb(refl,im,hfx,calcControls,parmDict)
    11221122    Icorr *= AbsCorr
    1123     ExtCorr = GetPwdrExt(refl,pfx,phfx,hfx,calcControls,parmDict)
     1123    ExtCorr = GetPwdrExt(refl,im,pfx,phfx,hfx,calcControls,parmDict)
    11241124    Icorr *= ExtCorr
    11251125    return Icorr,POcorr,AbsCorr,ExtCorr
    11261126   
    1127 def GetIntensityDerv(refl,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
     1127def GetIntensityDerv(refl,im,wave,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
    11281128    'Needs a doc string'    #need powder extinction derivs!
    11291129    dIdsh = 1./parmDict[hfx+'Scale']
    11301130    dIdsp = 1./parmDict[phfx+'Scale']
    11311131    if 'X' in parmDict[hfx+'Type']:
    1132         pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
     1132        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5+im],parmDict[hfx+'Azimuth'])
    11331133        dIdPola /= pola
    11341134    else:       #'N'
    … …  
    11381138    dIdPO = {}
    11391139    if pfx+'SHorder' in parmDict:
    1140         odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
     1140        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,im,g,pfx,hfx,SGData,calcControls,parmDict)
    11411141        for iSH in dFdODF:
    11421142            dFdODF[iSH] /= odfCor
    … …  
    11441144            dFdSA[i] /= odfCor
    11451145    elif calcControls[phfx+'poType'] == 'MD' or calcControls[phfx+'SHord']:
    1146         POcorr,dIdPO = GetPrefOriDerv(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)       
     1146        POcorr,dIdPO = GetPrefOriDerv(refl,im,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict)       
    11471147        for iPO in dIdPO:
    11481148            dIdPO[iPO] /= POcorr
    11491149    if 'T' in parmDict[hfx+'Type']:
    1150         dFdAb = GetAbsorbDerv(refl,hfx,calcControls,parmDict)*wave/refl[16] #wave/abs corr
    1151         dFdEx = GetPwdrExtDerv(refl,pfx,phfx,hfx,calcControls,parmDict)/refl[17]    #/ext corr
     1150        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[16+im] #wave/abs corr
     1151        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[17+im]    #/ext corr
    11521152    else:
    1153         dFdAb = GetAbsorbDerv(refl,hfx,calcControls,parmDict)*wave/refl[13] #wave/abs corr
    1154         dFdEx = GetPwdrExtDerv(refl,pfx,phfx,hfx,calcControls,parmDict)/refl[14]    #/ext corr       
     1153        dFdAb = GetAbsorbDerv(refl,im,hfx,calcControls,parmDict)*wave/refl[13+im] #wave/abs corr
     1154        dFdEx = GetPwdrExtDerv(refl,im,pfx,phfx,hfx,calcControls,parmDict)/refl[14+im]    #/ext corr       
    11551155    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx
    11561156       
    1157 def GetSampleSigGam(refl,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
     1157def GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
    11581158    'Needs a doc string'
    11591159    if 'C' in calcControls[hfx+'histType']:     #All checked & OK
    1160         costh = cosd(refl[5]/2.)
     1160        costh = cosd(refl[5+im]/2.)
    11611161        #crystallite size
    11621162        if calcControls[phfx+'SizeType'] == 'isotropic':
    … …  
    11751175        #microstrain               
    11761176        if calcControls[phfx+'MustrainType'] == 'isotropic':
    1177             Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi
     1177            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
    11781178        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
    11791179            H = np.array(refl[:3])
    … …  
    11821182            Si = parmDict[phfx+'Mustrain;i']
    11831183            Sa = parmDict[phfx+'Mustrain;a']
    1184             Mgam = 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
     1184            Mgam = 0.018*Si*Sa*tand(refl[5+im]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
    11851185        else:       #generalized - P.W. Stephens model
    11861186            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
    … …  
    11881188            for i,strm in enumerate(Strms):
    11891189                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
    1190             Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*np.sqrt(Sum)/np.pi
     1190            Mgam = 0.018*refl[4+im]**2*tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
    11911191    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
    11921192        #crystallite size
    11931193        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
    1194             Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4]**2/parmDict[phfx+'Size;i']
     1194            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
    11951195        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
    11961196            H = np.array(refl[:3])
    11971197            P = np.array(calcControls[phfx+'SizeAxis'])
    11981198            cosP,sinP = G2lat.CosSinAngle(H,P,G)
    1199             Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
     1199            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a'])
    12001200            Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
    12011201        else:           #ellipsoidal crystallites   #OK
    … …  
    12031203            H = np.array(refl[:3])
    12041204            lenR = G2pwd.ellipseSize(H,Sij,GB)
    1205             Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4]**2/lenR
     1205            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/lenR
    12061206        #microstrain               
    12071207        if calcControls[phfx+'MustrainType'] == 'isotropic':    #OK
    1208             Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4]*parmDict[phfx+'Mustrain;i']
     1208            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
    12091209        elif calcControls[phfx+'MustrainType'] == 'uniaxial':   #OK
    12101210            H = np.array(refl[:3])
    … …  
    12131213            Si = parmDict[phfx+'Mustrain;i']
    12141214            Sa = parmDict[phfx+'Mustrain;a']
    1215             Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
     1215            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa/np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
    12161216        else:       #generalized - P.W. Stephens model  OK
    12171217            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
    … …  
    12191219            for i,strm in enumerate(Strms):
    12201220                Sum += parmDict[phfx+'Mustrain:'+str(i)]*strm
    1221             Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4]**3
     1221            Mgam = 1.e-6*parmDict[hfx+'difC']*np.sqrt(Sum)*refl[4+im]**3
    12221222           
    12231223    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
    … …  
    12261226    return sig,gam
    12271227       
    1228 def GetSampleSigGamDerv(refl,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
     1228def GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict):
    12291229    'Needs a doc string'
    12301230    gamDict = {}
    12311231    sigDict = {}
    12321232    if 'C' in calcControls[hfx+'histType']:         #All checked & OK
    1233         costh = cosd(refl[5]/2.)
    1234         tanth = tand(refl[5]/2.)
     1233        costh = cosd(refl[5+im]/2.)
     1234        tanth = tand(refl[5+im]/2.)
    12351235        #crystallite size derivatives
    12361236        if calcControls[phfx+'SizeType'] == 'isotropic':
    … …  
    12671267        #microstrain derivatives               
    12681268        if calcControls[phfx+'MustrainType'] == 'isotropic':
    1269             Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi
     1269            Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5+im]/2.)/np.pi
    12701270            gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
    12711271            sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
    … …  
    12861286            sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
    12871287        else:       #generalized - P.W. Stephens model
    1288             const = 0.018*refl[4]**2*tanth/np.pi
     1288            const = 0.018*refl[4+im]**2*tanth/np.pi
    12891289            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
    12901290            Sum = 0
    … …  
    13011301    else:   #'T'OF - All checked & OK
    13021302        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
    1303             Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4]**2/parmDict[phfx+'Size;i']
     1303            Sgam = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2/parmDict[phfx+'Size;i']
    13041304            gamDict[phfx+'Size;i'] = -Sgam*parmDict[phfx+'Size;mx']/parmDict[phfx+'Size;i']
    13051305            sigDict[phfx+'Size;i'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])**2/(ateln2*parmDict[phfx+'Size;i'])
    13061306        elif calcControls[phfx+'SizeType'] == 'uniaxial':   #OK
    1307             const = 1.e-4*parmDict[hfx+'difC']*refl[4]**2
     1307            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
    13081308            H = np.array(refl[:3])
    13091309            P = np.array(calcControls[phfx+'SizeAxis'])
    … …  
    13211321            sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
    13221322        else:           #OK  ellipsoidal crystallites
    1323             const = 1.e-4*parmDict[hfx+'difC']*refl[4]**2
     1323            const = 1.e-4*parmDict[hfx+'difC']*refl[4+im]**2
    13241324            Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
    13251325            H = np.array(refl[:3])
    … …  
    13341334        #microstrain derivatives               
    13351335        if calcControls[phfx+'MustrainType'] == 'isotropic':
    1336             Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4]*parmDict[phfx+'Mustrain;i']
    1337             gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
     1336            Mgam = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*parmDict[phfx+'Mustrain;i']
     1337            gamDict[phfx+'Mustrain;i'] =  1.e-6*refl[4+im]*parmDict[hfx+'difC']*parmDict[phfx+'Mustrain;mx']   #OK
    13381338            sigDict[phfx+'Mustrain;i'] =  2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])**2/(ateln2*parmDict[phfx+'Mustrain;i'])       
    13391339        elif calcControls[phfx+'MustrainType'] == 'uniaxial':
    … …  
    13431343            Si = parmDict[phfx+'Mustrain;i']
    13441344            Sa = parmDict[phfx+'Mustrain;a']
    1345             gami = 1.e-6*parmDict[hfx+'difC']*refl[4]*Si*Sa
     1345            gami = 1.e-6*parmDict[hfx+'difC']*refl[4+im]*Si*Sa
    13461346            sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
    13471347            Mgam = gami/sqtrm
    … …  
    13551355            pwrs = calcControls[phfx+'MuPwrs']
    13561356            Strms = G2spc.MustrainCoeff(refl[:3],SGData)
    1357             const = 1.e-6*parmDict[hfx+'difC']*refl[4]**3
     1357            const = 1.e-6*parmDict[hfx+'difC']*refl[4+im]**3
    13581358            Sum = 0
    13591359            for i,strm in enumerate(Strms):
    … …  
    13701370    return sigDict,gamDict
    13711371       
    1372 def GetReflPos(refl,wave,A,hfx,calcControls,parmDict):
     1372def GetReflPos(refl,im,wave,A,hfx,calcControls,parmDict):
    13731373    'Needs a doc string'
    13741374    h,k,l = refl[:3]
    13751375    d = 1./np.sqrt(G2lat.calc_rDsq(np.array([h,k,l]),A))
    13761376
    1377     refl[4] = d
     1377    refl[4+im] = d
    13781378    if 'C' in calcControls[hfx+'histType']:
    13791379        pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
    … …  
    13891389    return pos
    13901390
    1391 def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
     1391def GetReflPosDerv(refl,im,wave,A,hfx,calcControls,parmDict):
    13921392    'Needs a doc string'
    13931393    dpr = 180./np.pi
    … …  
    13971397    dsp = 1./dst
    13981398    if 'C' in calcControls[hfx+'histType']:
    1399         pos = refl[5]-parmDict[hfx+'Zero']
     1399        pos = refl[5+im]-parmDict[hfx+'Zero']
    14001400        const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
    14011401        dpdw = const*dst
    … …  
    14191419        return dpdA,dpdZ,dpdDC,dpdDA,dpdDB
    14201420           
    1421 def GetHStrainShift(refl,SGData,phfx,hfx,calcControls,parmDict):
     1421def GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict):
    14221422    'Needs a doc string'
    14231423    laue = SGData['SGLaue']
    … …  
    14261426    if laue in ['m3','m3m']:
    14271427        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
    1428             refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
     1428            refl[4+im]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
    14291429    elif laue in ['6/m','6/mmm','3m1','31m','3']:
    14301430        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
    … …  
    14471447            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
    14481448    if 'C' in calcControls[hfx+'histType']:
    1449         return -180.*Dij*refl[4]**2*tand(refl[5]/2.0)/np.pi
     1449        return -180.*Dij*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
    14501450    else:
    1451         return -Dij*parmDict[hfx+'difC']*refl[4]**2/2.
     1451        return -Dij*parmDict[hfx+'difC']*refl[4+im]**2/2.
    14521452           
    1453 def GetHStrainShiftDerv(refl,SGData,phfx,hfx,calcControls,parmDict):
     1453def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
    14541454    'Needs a doc string'
    14551455    laue = SGData['SGLaue']
    … …  
    14581458    if laue in ['m3','m3m']:
    14591459        dDijDict = {phfx+'D11':h**2+k**2+l**2,
    1460             phfx+'eA':refl[4]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
     1460            phfx+'eA':refl[4+im]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
    14611461    elif laue in ['6/m','6/mmm','3m1','31m','3']:
    14621462        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
    … …  
    14811481    if 'C' in calcControls[hfx+'histType']:
    14821482        for item in dDijDict:
    1483             dDijDict[item] *= 180.0*refl[4]**2*tand(refl[5]/2.0)/np.pi
     1483            dDijDict[item] *= 180.0*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
    14841484    else:
    14851485        for item in dDijDict:
    1486             dDijDict[item] *= -parmDict[hfx+'difC']*refl[4]**3/2.
     1486            dDijDict[item] *= -parmDict[hfx+'difC']*refl[4+im]**3/2.
    14871487    return dDijDict
    14881488   
    … …  
    15191519            for phase in refLists:
    15201520                Phase = Phases[phase]
     1521                im = 0
     1522                if Phase['General']['Type'] in ['modulated','magnetic']:
     1523                    im = 1
    15211524                pId = Phase['pId']
    15221525                phfx = '%d:%d:'%(pId,hId)
    … …  
    15291532                    if 'C' in calcControls[hfx+'histType']:
    15301533                        yp = np.zeros_like(yb)
    1531                         Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
    1532                         iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
    1533                         iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
     1534                        Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
     1535                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
     1536                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
    15341537                        iFin2 = iFin
    15351538                        if not iBeg+iFin:       #peak below low limit - skip peak
    … …  
    15381541                            break
    15391542                        elif iBeg < iFin:
    1540                             yp[iBeg:iFin] = refl[11]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))    #>90% of time spent here
     1543                            yp[iBeg:iFin] = refl[11+im]*refl[9+im]*G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))    #>90% of time spent here
    15411544                            if Ka2:
    1542                                 pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
    1543                                 Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
     1545                                pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
     1546                                Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
    15441547                                iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
    15451548                                iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
    1546                                 yp[iBeg2:iFin2] += refl[11]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg2:iFin2]))        #and here
    1547                             refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11]*(1.+kRatio)),0.0))
     1549                                yp[iBeg2:iFin2] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2]))        #and here
     1550                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11+im]*(1.+kRatio)),0.0))
    15481551                    elif 'T' in calcControls[hfx+'histType']:
    15491552                        yp = np.zeros_like(yb)
    1550                         Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7])
    1551                         iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
    1552                         iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
     1553                        Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
     1554                        iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
     1555                        iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
    15531556                        if iBeg < iFin:
    1554                             yp[iBeg:iFin] = refl[11]*refl[9]*G2pwd.getEpsVoigt(refl[5],refl[12],refl[13],refl[6],refl[7],ma.getdata(x[iBeg:iFin]))  #>90% of time spent here
    1555                             refl[8] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11],0.0))
    1556                     Fo = np.sqrt(np.abs(refl[8]))
    1557                     Fc = np.sqrt(np.abs(refl[9]))
     1557                            yp[iBeg:iFin] = refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))  #>90% of time spent here
     1558                            refl[8+im] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
     1559                    Fo = np.sqrt(np.abs(refl[8+im]))
     1560                    Fc = np.sqrt(np.abs(refl[9]+im))
    15581561                    sumFo += Fo
    1559                     sumFosq += refl[8]**2
     1562                    sumFosq += refl[8+im]**2
    15601563                    sumdF += np.abs(Fo-Fc)
    1561                     sumdFsq += (refl[8]-refl[9])**2
     1564                    sumdFsq += (refl[8+im]-refl[9+im])**2
    15621565                Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
    15631566                Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
    … …  
    15711574    'Needs a doc string'
    15721575   
    1573     def GetReflSigGamCW(refl,wave,G,GB,phfx,calcControls,parmDict):
     1576    def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
    15741577        U = parmDict[hfx+'U']
    15751578        V = parmDict[hfx+'V']
    … …  
    15771580        X = parmDict[hfx+'X']
    15781581        Y = parmDict[hfx+'Y']
    1579         tanPos = tand(refl[5]/2.0)
    1580         Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
     1582        tanPos = tand(refl[5+im]/2.0)
     1583        Ssig,Sgam = GetSampleSigGam(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
    15811584        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
    15821585        sig = max(0.001,sig)
    1583         gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam     #save peak gamma
     1586        gam = X/cosd(refl[5+im]/2.0)+Y*tanPos+Sgam     #save peak gamma
    15841587        gam = max(0.001,gam)
    15851588        return sig,gam
    15861589               
    1587     def GetReflSigGamTOF(refl,G,GB,phfx,calcControls,parmDict):
    1588         sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4]**2+   \
    1589             parmDict[hfx+'sig-2']*refl[4]**4+parmDict[hfx+'sig-q']/refl[4]**2
    1590         gam = parmDict[hfx+'X']*refl[4]+parmDict[hfx+'Y']*refl[4]**2
    1591         Ssig,Sgam = GetSampleSigGam(refl,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
     1590    def GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict):
     1591        sig = parmDict[hfx+'sig-0']+parmDict[hfx+'sig-1']*refl[4+im]**2+   \
     1592            parmDict[hfx+'sig-2']*refl[4+im]**4+parmDict[hfx+'sig-q']/refl[4+im]**2
     1593        gam = parmDict[hfx+'X']*refl[4+im]+parmDict[hfx+'Y']*refl[4+im]**2
     1594        Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
    15921595        sig += Ssig
    15931596        gam += Sgam
    15941597        return sig,gam
    15951598       
    1596     def GetReflAlpBet(refl,hfx,parmDict):
    1597         alp = parmDict[hfx+'alpha']/refl[4]
    1598         bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4]**4+parmDict[hfx+'beta-q']/refl[4]**2
     1599    def GetReflAlpBet(refl,im,hfx,parmDict):
     1600        alp = parmDict[hfx+'alpha']/refl[4+im]
     1601        bet = parmDict[hfx+'beta-0']+parmDict[hfx+'beta-1']/refl[4+im]**4+parmDict[hfx+'beta-q']/refl[4+im]**2
    15991602        return alp,bet
    16001603       
    … …  
    16261629        SGData = Phase['General']['SGData']
    16271630        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     1631        im = 0
     1632        if Phase['General']['Type'] in ['modulated','magnetic']:
     1633            SSGData = Phase['General']['SSGData']
     1634            SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
     1635            im = 1  #offset in SS reflection list
     1636            #??
    16281637        Dij = GetDij(phfx,SGData,parmDict)
    16291638        A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
    … …  
    16391648        for iref,refl in enumerate(refDict['RefList']):
    16401649            if 'C' in calcControls[hfx+'histType']:
    1641                 h,k,l = refl[:3]
     1650                if im:
     1651                    h,k,l,m = refl[:4]
     1652                else:
     1653                    h,k,l = refl[:3]
    16421654                Uniq = np.inner(refl[:3],SGMT)
    1643                 refl[5] = GetReflPos(refl,wave,A,hfx,calcControls,parmDict)         #corrected reflection position
    1644                 Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
    1645 #                refl[5] += GetHStrainShift(refl,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
    1646                 refl[6:8] = GetReflSigGamCW(refl,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
    1647                 refl[11:15] = GetIntensityCorr(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
    1648                 refl[11] *= Vst*Lorenz
     1655                refl[5+im] = GetReflPos(refl,im,wave,A,hfx,calcControls,parmDict)         #corrected reflection position
     1656                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
     1657#                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
     1658                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
     1659                refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
     1660                refl[11+im] *= Vst*Lorenz
    16491661                 
    16501662                if Phase['General'].get('doPawley'):
    16511663                    try:
    1652                         pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
    1653                         refl[9] = parmDict[pInd]
     1664                        if im:
     1665                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
     1666                        else:
     1667                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
     1668                        refl[9+im] = parmDict[pInd]
    16541669                    except KeyError:
    16551670#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
    16561671                        continue
    1657                 Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
    1658                 iBeg = np.searchsorted(x,refl[5]-fmin)
    1659                 iFin = np.searchsorted(x,refl[5]+fmax)
     1672                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
     1673                iBeg = np.searchsorted(x,refl[5+im]-fmin)
     1674                iFin = np.searchsorted(x,refl[5+im]+fmax)
    16601675                if not iBeg+iFin:       #peak below low limit - skip peak
    16611676                    continue
    … …  
    16651680                    badPeak = True
    16661681                    continue
    1667                 yc[iBeg:iFin] += refl[11]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))    #>90% of time spent here
     1682                yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))    #>90% of time spent here
    16681683                if Ka2:
    1669                     pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
    1670                     Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
     1684                    pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
     1685                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
    16711686                    iBeg = np.searchsorted(x,pos2-fmin)
    16721687                    iFin = np.searchsorted(x,pos2+fmax)
    … …  
    16771692                    elif iBeg > iFin:   #bad peak coeff - skip
    16781693                        continue
    1679                     yc[iBeg:iFin] += refl[11]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))        #and here
     1694                    yc[iBeg:iFin] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))        #and here
    16801695            elif 'T' in calcControls[hfx+'histType']:
    16811696                h,k,l = refl[:3]
    16821697                Uniq = np.inner(refl[:3],SGMT)
    1683                 refl[5] = GetReflPos(refl,0.0,A,hfx,calcControls,parmDict)         #corrected reflection position
    1684                 Lorenz = sind(parmDict[hfx+'2-theta']/2)*refl[4]**4                                                #TOF Lorentz correction
    1685 #                refl[5] += GetHStrainShift(refl,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
    1686                 refl[6:8] = GetReflSigGamTOF(refl,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
    1687                 refl[12:14] = GetReflAlpBet(refl,hfx,parmDict)
    1688                 refl[11],refl[15],refl[16],refl[17] = GetIntensityCorr(refl,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
    1689                 refl[11] *= Vst*Lorenz
     1698                refl[5+im] = GetReflPos(refl,im,0.0,A,hfx,calcControls,parmDict)         #corrected reflection position
     1699                Lorenz = sind(parmDict[hfx+'2-theta']/2)*refl[4+im]**4                                                #TOF Lorentz correction
     1700#                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
     1701                refl[6+im:8+im] = GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
     1702                refl[12+im:14+im] = GetReflAlpBet(refl,im,hfx,parmDict)
     1703                refl[11+im],refl[15+im],refl[16+im],refl[17+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
     1704                refl[11+im] *= Vst*Lorenz
    16901705                if Phase['General'].get('doPawley'):
    16911706                    try:
    1692                         pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
    1693                         refl[9] = parmDict[pInd]
     1707                        if im:
     1708                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
     1709                        else:
     1710                            pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
     1711                        refl[9+im] = parmDict[pInd]
    16941712                    except KeyError:
    16951713#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
    16961714                        continue
    1697                 Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7])
    1698                 iBeg = np.searchsorted(x,refl[5]-fmin)
    1699                 iFin = np.searchsorted(x,refl[5]+fmax)
     1715                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
     1716                iBeg = np.searchsorted(x,refl[5+im]-fmin)
     1717                iFin = np.searchsorted(x,refl[5+im]+fmax)
    17001718                if not iBeg+iFin:       #peak below low limit - skip peak
    17011719                    continue
    … …  
    17051723                    badPeak = True
    17061724                    continue
    1707                 yc[iBeg:iFin] += refl[11]*refl[9]*G2pwd.getEpsVoigt(refl[5],refl[12],refl[13],refl[6],refl[7],ma.getdata(x[iBeg:iFin]))/cw[iBeg:iFin]
     1725                yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))/cw[iBeg:iFin]
    17081726#        print 'profile calc time: %.3fs'%(time.time()-time0)
    17091727    if badPeak:
    … …  
    17831801        SGData = Phase['General']['SGData']
    17841802        SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     1803        im = 0
     1804        if Phase['General']['Type'] in ['modulated','magnetic']:
     1805            SSGData = Phase['General']['SSGData']
     1806            SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
     1807            im = 1  #offset in SS reflection list
     1808            #??
    17851809        pId = Phase['pId']
    17861810        pfx = '%d::'%(pId)
    … …  
    18001824            Uniq = np.inner(refl[:3],SGMT)
    18011825            if 'T' in calcControls[hfx+'histType']:
    1802                 wave = refl[14]
    1803             dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
     1826                wave = refl[14+im]
     1827            dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
    18041828            if 'C' in calcControls[hfx+'histType']:        #CW powder
    1805                 Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
     1829                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
    18061830            else: #'T'OF
    1807                 Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5],refl[12],refl[13],refl[6],refl[7])
    1808             iBeg = np.searchsorted(x,refl[5]-fmin)
    1809             iFin = np.searchsorted(x,refl[5]+fmax)
     1831                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
     1832            iBeg = np.searchsorted(x,refl[5+im]-fmin)
     1833            iFin = np.searchsorted(x,refl[5+im]+fmax)
    18101834            if not iBeg+iFin:       #peak below low limit - skip peak
    18111835                continue
    18121836            elif not iBeg-iFin:     #peak above high limit - done
    18131837                break
    1814             pos = refl[5]
     1838            pos = refl[5+im]
    18151839            if 'C' in calcControls[hfx+'histType']:
    18161840                tanth = tand(pos/2.0)
    … …  
    18181842                lenBF = iFin-iBeg
    18191843                dMdpk = np.zeros(shape=(6,lenBF))
    1820                 dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin]))
     1844                dMdipk = G2pwd.getdFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))
    18211845                for i in range(5):
    1822                     dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11]*refl[9]*dMdipk[i]
     1846                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11+im]*refl[9+im]*dMdipk[i]
    18231847                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
    18241848                if Ka2:
    1825                     pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
     1849                    pos2 = refl[5+im]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
    18261850                    iBeg2 = np.searchsorted(x,pos2-fmin)
    18271851                    iFin2 = np.searchsorted(x,pos2+fmax)
    … …  
    18291853                        lenBF2 = iFin2-iBeg2
    18301854                        dMdpk2 = np.zeros(shape=(6,lenBF2))
    1831                         dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg2:iFin2]))
     1855                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2]))
    18321856                        for i in range(5):
    1833                             dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[11]*refl[9]*kRatio*dMdipk2[i]
    1834                         dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[11]*dMdipk2[0]
     1857                            dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[11+im]*refl[9+im]*kRatio*dMdipk2[i]
     1858                        dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[11+im]*dMdipk2[0]
    18351859                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
    18361860            else:   #'T'OF
    … …  
    18391863                    break
    18401864                dMdpk = np.zeros(shape=(6,lenBF))
    1841                 dMdipk = G2pwd.getdEpsVoigt(refl[5],refl[12],refl[13],refl[6],refl[7],ma.getdata(x[iBeg:iFin]))
     1865                dMdipk = G2pwd.getdEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))
    18421866                for i in range(6):
    1843                     dMdpk[i] += refl[11]*refl[9]*dMdipk[i]      #cw[iBeg:iFin]*
     1867                    dMdpk[i] += refl[11+im]*refl[9+im]*dMdipk[i]      #cw[iBeg:iFin]*
    18441868                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}           
    18451869            if Phase['General'].get('doPawley'):
    18461870                dMdpw = np.zeros(len(x))
    18471871                try:
    1848                     pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
     1872                    if im:
     1873                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
     1874                    else:
     1875                        pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
    18491876                    idx = varylist.index(pIdx)
    1850                     dMdpw[iBeg:iFin] = dervDict['int']/refl[9]
     1877                    dMdpw[iBeg:iFin] = dervDict['int']/refl[9+im]
    18511878                    if Ka2: #not for TOF either
    1852                         dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9]
     1879                        dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9+im]
    18531880                    dMdv[idx] = dMdpw
    18541881                except: # ValueError:
    18551882                    pass
    18561883            if 'C' in calcControls[hfx+'histType']:
    1857                 dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
     1884                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,im,wave,A,hfx,calcControls,parmDict)
    18581885                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
    18591886                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
    … …  
    18711898                names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
    18721899                    hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
    1873                     hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4],'gam'],hfx+'Y':[refl[4]**2,'gam'],
    1874                     hfx+'alpha':[1./refl[4],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4]**4,'bet'],
    1875                     hfx+'beta-q':[1./refl[4]**2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4]**2,'sig'],
    1876                     hfx+'sig-2':[refl[4]**4,'sig'],hfx+'sig-q':[1./refl[4]**2,'sig'],
     1900                    hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4+im],'gam'],hfx+'Y':[refl[4+im]**2,'gam'],
     1901                    hfx+'alpha':[1./refl[4+im],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4+im]**4,'bet'],
     1902                    hfx+'beta-q':[1./refl[4+im]**2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4+im]**2,'sig'],
     1903                    hfx+'sig-2':[refl[4+im]**4,'sig'],hfx+'sig-q':[1./refl[4+im]**2,'sig'],
    18771904                    hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
    18781905            for name in names:
    … …  
    19241951                    if Ka2:
    19251952                        depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
    1926             dDijDict = GetHStrainShiftDerv(refl,SGData,phfx,hfx,calcControls,parmDict)
     1953            dDijDict = GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict)
    19271954            for name in dDijDict:
    19281955                if name in varylist:
    … …  
    19351962                        depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
    19361963            if 'C' in calcControls[hfx+'histType']:
    1937                 sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
     1964                sigDict,gamDict = GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
    19381965            else:   #'T'OF
    1939                 sigDict,gamDict = GetSampleSigGamDerv(refl,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
     1966                sigDict,gamDict = GetSampleSigGamDerv(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
    19401967            for name in gamDict:
    19411968                if name in varylist:
    … …  
    19571984                        depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
    19581985            for name in ['BabA','BabU']:
    1959                 if refl[9]:
     1986                if refl[9+im]:
    19601987                    if phfx+name in varylist:
    1961                         dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9]
     1988                        dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9+im]
    19621989                        if Ka2:
    1963                             dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9]
     1990                            dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9+im]
    19641991                    elif phfx+name in dependentVars:                   
    1965                         depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9]
     1992                        depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']/refl[9+im]
    19661993                        if Ka2:
    1967                             depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9]                 
     1994                            depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']/refl[9+im]                 
    19681995            if not Phase['General'].get('doPawley'):
    19691996                #do atom derivatives -  for RB,F,X & U so far             
    1970                 corr = dervDict['int']/refl[9]
     1997                corr = dervDict['int']/refl[9+im]
    19711998                if Ka2:
    1972                     corr2 = dervDict2['int']/refl[9]
     1999                    corr2 = dervDict2['int']/refl[9+im]
    19732000                for name in varylist+dependentVars:
    19742001                    if '::RBV;' in name:
    … …  
    20062033    phfx = '%d:%d:'%(Phase['pId'],hId)
    20072034    SGData = Phase['General']['SGData']
     2035    im = 0
     2036    if Phase['General']['Type'] in ['modulated','magnetic']:
     2037        SSGData = Phase['General']['SSGData']
     2038        SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
     2039        im = 1  #offset in SS reflection list
     2040        #??
    20082041    A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
    20092042    G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
    … …  
    20192052    if calcControls['F**2']:
    20202053        for iref,ref in enumerate(refDict['RefList']):
    2021             if ref[6] > 0:
    2022                 dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1]
    2023                 w = 1.0/ref[6]
    2024                 if w*ref[5] >= calcControls['minF/sig']:
    2025                     wdf[iref] = w*(ref[5]-ref[7])
     2054            if ref[6+im] > 0:
     2055                dervDict = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1]
     2056                w = 1.0/ref[6+im]
     2057                if w*ref[5+im] >= calcControls['minF/sig']:
     2058                    wdf[iref] = w*(ref[5+im]-ref[7+im])
    20262059                    for j,var in enumerate(varylist):
    20272060                        if var in dFdvDict:
    2028                             dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11]
     2061                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
    20292062                    for var in dependentVars:
    20302063                        if var in dFdvDict:
    2031                             depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11]
     2064                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
    20322065                    if phfx+'Scale' in varylist:
    2033                         dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*ref[11]
     2066                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9+im]*ref[11+im]
    20342067                    elif phfx+'Scale' in dependentVars:
    2035                         depDerivDict[phfx+'Scale'][iref] = w*ref[9]*ref[11]
     2068                        depDerivDict[phfx+'Scale'][iref] = w*ref[9+im]*ref[11+im]
    20362069                    for item in ['Ep','Es','Eg']:
    20372070                        if phfx+item in varylist and phfx+item in dervDict:
    2038                             dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11]  #OK
     2071                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im]  #OK
    20392072                        elif phfx+item in dependentVars and phfx+item in dervDict:
    2040                             depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11]  #OK
     2073                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]  #OK
    20412074                    for item in ['BabA','BabU']:
    20422075                        if phfx+item in varylist:
    2043                             dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11]
     2076                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
    20442077                        elif phfx+item in dependentVars:
    2045                             depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11]
     2078                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
    20462079    else:
    20472080        for iref,ref in enumerate(refDict['RefList']):
    2048             if ref[5] > 0.:
     2081            if ref[5+im] > 0.:
    20492082                dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist+dependentVars)[1]
    2050                 Fo = np.sqrt(ref[5])
    2051                 Fc = np.sqrt(ref[7])
    2052                 w = 1.0/ref[6]
     2083                Fo = np.sqrt(ref[5+im])
     2084                Fc = np.sqrt(ref[7+im])
     2085                w = 1.0/ref[6+im]
    20532086                if 2.0*Fo*w*Fo >= calcControls['minF/sig']:
    20542087                    wdf[iref] = 2.0*Fo*w*(Fo-Fc)
    20552088                    for j,var in enumerate(varylist):
    20562089                        if var in dFdvDict:
    2057                             dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11]
     2090                            dMdvh[j][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
    20582091                    for var in dependentVars:
    20592092                        if var in dFdvDict:
    2060                             depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11]
     2093                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
    20612094                    if phfx+'Scale' in varylist:
    2062                         dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9]*ref[11]
     2095                        dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9+im]*ref[11+im]
    20632096                    elif phfx+'Scale' in dependentVars:
    2064                         depDerivDict[phfx+'Scale'][iref] = w*ref[9]*ref[11]                           
     2097                        depDerivDict[phfx+'Scale'][iref] = w*ref[9+im]*ref[11+im]                           
    20652098                    for item in ['Ep','Es','Eg']:
    20662099                        if phfx+item in varylist and phfx+item in dervDict:
    2067                             dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11]  #correct
     2100                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im]  #correct
    20682101                        elif phfx+item in dependentVars and phfx+item in dervDict:
    2069                             depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11]
     2102                            depDerivDict[phfx+item][iref] = w*dervDict[phfx+item]/ref[11+im]
    20702103                    for item in ['BabA','BabU']:
    20712104                        if phfx+item in varylist:
    2072                             dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11]
     2105                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
    20732106                        elif phfx+item in dependentVars:
    2074                             depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11]
     2107                            depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
    20752108    return dMdvh,depDerivDict,wdf
    20762109   
    … …  
    22692302            phfx = '%d:%d:'%(Phase['pId'],hId)
    22702303            SGData = Phase['General']['SGData']
     2304            im = 0
     2305            if Phase['General']['Type'] in ['modulated','magnetic']:
     2306                SSGData = Phase['General']['SSGData']
     2307                SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
     2308                im = 1  #offset in SS reflection list
     2309                #??
    22712310            A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
    22722311            G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
    … …  
    22842323            if calcControls['F**2']:
    22852324                for i,ref in enumerate(refDict['RefList']):
    2286                     if ref[6] > 0:
    2287                         ref[11] = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
    2288                         w = 1.0/ref[6]
    2289                         ref[7] = parmDict[phfx+'Scale']*ref[9]*ref[11]  #correct Fc^2 for extinction
    2290                         ref[8] = ref[5]/(parmDict[phfx+'Scale']*ref[11])
    2291                         if w*ref[5] >= calcControls['minF/sig']:
    2292                             Fo = np.sqrt(ref[5])
     2325                    if ref[6+im] > 0:
     2326                        ref[11+im] = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
     2327                        w = 1.0/ref[6+im]
     2328                        ref[7+im] = parmDict[phfx+'Scale']*ref[9+im]*ref[11+im]  #correct Fc^2 for extinction
     2329                        ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
     2330                        if w*ref[5+im] >= calcControls['minF/sig']:
     2331                            Fo = np.sqrt(ref[5+im])
    22932332                            sumFo += Fo
    2294                             sumFo2 += ref[5]
    2295                             sumdF += abs(Fo-np.sqrt(ref[7]))
    2296                             sumdF2 += abs(ref[5]-ref[7])
     2333                            sumFo2 += ref[5+im]
     2334                            sumdF += abs(Fo-np.sqrt(ref[7+im]))
     2335                            sumdF2 += abs(ref[5+im]-ref[7+im])
    22972336                            nobs += 1
    2298                             df[i] = -w*(ref[5]-ref[7])
    2299                             sumwYo += (w*ref[5])**2
     2337                            df[i] = -w*(ref[5+im]-ref[7+im])
     2338                            sumwYo += (w*ref[5+im])**2
    23002339            else:
    23012340                for i,ref in enumerate(refDict['RefList']):
    2302                     if ref[5] > 0.:
    2303                         ref[11] = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
    2304                         ref[7] = parmDict[phfx+'Scale']*ref[9]*ref[11]    #correct Fc^2 for extinction
    2305                         ref[8] = ref[5]/(parmDict[phfx+'Scale']*ref[11])
    2306                         Fo = np.sqrt(ref[5])
    2307                         Fc = np.sqrt(ref[7])
    2308                         w = 2.0*Fo/ref[6]
     2341                    if ref[5+im] > 0.:
     2342                        ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
     2343                        ref[7+im] = parmDict[phfx+'Scale']*ref[9+im]*ref[11+im]    #correct Fc^2 for extinction
     2344                        ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
     2345                        Fo = np.sqrt(ref[5+im])
     2346                        Fc = np.sqrt(ref[7+im])
     2347                        w = 2.0*Fo/ref[6+im]
    23092348                        if w*Fo >= calcControls['minF/sig']:
    23102349                            sumFo += Fo
    2311                             sumFo2 += ref[5]
     2350                            sumFo2 += ref[5+im]
    23122351                            sumdF += abs(Fo-Fc)
    2313                             sumdF2 += abs(ref[5]-ref[7])
     2352                            sumdF2 += abs(ref[5+im]-ref[7+im])
    23142353                            nobs += 1
    23152354                            df[i] = -w*(Fo-Fc)
Note: See TracChangeset for help on using the changeset viewer.