Changeset 795


Ignore:
Timestamp:
Nov 3, 2012 3:22:36 PM (10 years ago)
Author:
vondreele
Message:

implement TOF input, peak search & fitting
auto peak search
convert instrument parms to dictionary
add charge flip on max rho

Location:
trunk
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASII.py

    r794 r795  
    6161import GSASIIgrid as G2gd
    6262import GSASIIplot as G2plt
     63import GSASIIpwd as G2pwd
    6364import GSASIIpwdGUI as G2pdG
    6465import GSASIIspc as G2spc
    … …  
    169170            id=wx.ID_ANY, kind=wx.ITEM_NORMAL,text='Make new PDFs')
    170171        self.MakePDF.append(item)
    171         item.Enable(False)
     172#        item.Enable(False)
    172173        self.Bind(wx.EVT_MENU, self.OnMakePDFs, id=item.GetId())
    173174       
    … …  
    557558            S = File.readline()               
    558559        File.close()
    559         inst = dict(zip(newItems,zip(newVals,newVals,len(newVals)*[False,])))
    560         for item in inst:
    561             inst[item] = list(inst[item])
    562         return inst
     560        return G2IO.makeInstDict(newItems,newVals,len(newVals)*[False,])
    563561       
    564562    def ReadPowderIparm(self,instfile,bank,databanks,rd):
    … …  
    658656                    data.extend([0.0,0.0,0.002,azm])                                      #OK defaults if fxn #3 not 1st in iprm file
    659657                codes.extend([0,0,0,0,0,0,0])
    660                 inst = dict(zip(names,zip(data,data,codes)))
    661                 for item in inst:
    662                     inst[item] = list(inst[item])
    663                 return inst
     658                return G2IO.makeInstDict(names,data,codes)
    664659            elif 'T' in DataType:
    665660                names = ['Type','2-theta','difC','difA','Zero','alpha','beta-0','beta-1','var-inst','X','Y','Azimuth']
    666661                codes = [0,0,0,0,0,0,0,0,0,0,0,0]
    667                 azm = Iparm.get('INS  1DETAZM')
    668                 if azm is None: #not in this Iparm file
    669                     azm = 0.0
    670                 else:
    671                     azm = float(azm)
     662                azm = 0.
     663                if 'INS  1DETAZM' in Iparm:
     664                    azm = float(Iparm['INS  1DETAZM'])
    672665                s = Iparm['INS  1BNKPAR'].split()
    673666                data.extend([G2IO.sfloat(s[1]),])               #2-theta for bank
    … …  
    677670                pfType = int(s[0])
    678671                s = Iparm['INS  1PRCF11'].split()
    679                 if pfType == 1:
     672                if abs(pfType) == 1:
    680673                    data.extend([G2IO.sfloat(s[1]),G2IO.sfloat(s[2]),G2IO.sfloat(s[3])])
    681674                    s = Iparm['INS  1PRCF12'].split()
    682675                    data.extend([G2IO.sfloat(s[1]),0.0,0.0,azm])
    683                 elif pfType in [3,4,5]:
     676                elif abs(pfType) in [3,4,5]:
    684677                    data.extend([G2IO.sfloat(s[0]),G2IO.sfloat(s[1]),G2IO.sfloat(s[2])])
    685                     if pfType == 4:
     678                    if abs(pfType) == 4:
    686679                        data.extend([G2IO.sfloat(s[3]),0.0,0.0,azm])
    687680                    else:
    688681                        s = Iparm['INS  1PRCF12'].split()
    689                         data.extend([G2IO.sfloat(s[0]),0.0,0.0,azm])
    690                 inst = dict(zip(names,zip(data,data,codes)))
    691                 for item in inst:
    692                     inst[item] = list(inst[item])
    693                 return inst
     682                        data.extend([G2IO.sfloat(s[0]),0.0,0.0,azm])                       
     683                Inst = G2IO.makeInstDict(names,data,codes)
     684                if pfType < 0:
     685                    Ipab = 'INS  1PAB'+str(-pfType)
     686                    Npab = int(Iparm[Ipab+'  '].strip())
     687                    Inst['Pdabc'] = []
     688                    for i in range(Npab):
     689                        k = Ipab+str(i+1).rjust(2)
     690                        s = Iparm[k].split()
     691                        Inst['Pdabc'].append([float(t) for t in s])
     692                    Inst['Pdabc'] = np.array(Inst['Pdabc'])
     693                if 'INS  1I ITYP' in Iparm:
     694                    Ityp = int(Iparm['INS  1I ITYP'].split()[0])
     695                    if Ityp in [1,2,3,4,5]:
     696                        Inst['Itype'] = Ityp
     697                        Icoeff = []
     698                        Iesd = []
     699                        Icovar = []                   
     700                        for i in range(3):
     701                            s = Iparm['INS  1ICOFF'+str(i+1)].split()
     702                            Icoeff += [float(S) for S in s]
     703                            s = Iparm['INS  1IECOF'+str(i+1)].split()
     704                            Iesd += [float(S) for S in s]
     705                        for i in range(8):
     706                            s = Iparm['INS  1IECOR'+str(i+1)].split()
     707                            Icovar += [float(S) for S in s]
     708                        Inst['Icoeff'] = Icoeff
     709                        Inst['Iesd'] = Iesd
     710                        Inst['Icovar'] = Icovar
     711                return Inst
    694712
    695713        # stuff we might need from the reader
    … …  
    846864                self.PatternTree.AppendItem(Id,text='Comments'),
    847865                rd.comments)
     866            if 'T' in instParmList['Type'][0]:
     867                if not rd.clockWd and rd.GSAS:
     868                    rd.powderdata[0] *= 100.
     869                cw = np.diff(rd.powderdata[0])
     870                rd.powderdata[0] = rd.powderdata[0][:-1]+cw/2.
     871                rd.powderdata[1] = rd.powderdata[1][:-1]/cw
     872                rd.powderdata[2] = rd.powderdata[2][:-1]
     873                rd.powderdata[3] = rd.powderdata[3][:-1]
     874                rd.powderdata[4] = rd.powderdata[4][:-1]
     875                rd.powderdata[5] = rd.powderdata[5][:-1]
     876                if 'Itype' in instParmList:
     877                    incident = G2pwd.calcIncident(instParmList,rd.powderdata[0])
    848878            Tmin = min(rd.powderdata[0])
    849879            Tmax = max(rd.powderdata[0])
    … …  
    11291159                names = ['Type','Lam','Zero']
    11301160                codes = [0,0,0]
    1131                 inst = dict(zip(names,zip(data,data,codes)))
    1132                 for item in inst:
    1133                     inst[item] = list(inst[item])
     1161                inst = G2IO.makeInstDict(names,data,codes)
    11341162                self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Instrument Parameters'),inst)
    11351163                self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Comments'),comments)

  • trunk/GSASIIIO.py

    r794 r795  
    3737    else:
    3838        return 0
     39
     40def makeInstDict(names,data,codes):
     41    inst = dict(zip(names,zip(data,data,codes)))
     42    for item in inst:
     43        inst[item] = list(inst[item])
     44    return inst
     45
    3946
    4047def FileDlgFixExt(dlg,file):
    … …  
    443450#patch
    444451                if datus[0] == 'Instrument Parameters' and not isinstance(datus[1],dict):
    445                     datus[1] = dict(zip(datus[1][3],zip(datus[1][0],datus[1][1],datus[1][2])))
     452                    if 'PWDR' in datum[0]:
     453                        datus[1] = dict(zip(datus[1][3],zip(datus[1][0],datus[1][1],datus[1][2])))
     454                    else:
     455                        datus[1] = dict(zip(datus[1][2],zip(datus[1][0],datus[1][1])))
    446456                    for item in datus[1]:               #zip makes tuples - now make lists!
    447457                        datus[1][item] = list(datus[1][item])

  • trunk/GSASIIgrid.py

    r792 r795  
    8686] = [wx.NewId() for item in range(6)]
    8787
    88 [ wxID_UNDO,wxID_LSQPEAKFIT,wxID_LSQONECYCLE,wxID_RESETSIGGAM,wxID_CLEARPEAKS,
    89 ] = [wx.NewId() for item in range(5)]
     88[ wxID_UNDO,wxID_LSQPEAKFIT,wxID_LSQONECYCLE,wxID_RESETSIGGAM,wxID_CLEARPEAKS,wxID_AUTOSEARCH,
     89] = [wx.NewId() for item in range(6)]
    9090
    9191[  wxID_INDXRELOAD, wxID_INDEXPEAKS, wxID_REFINECELL, wxID_COPYCELL, wxID_MAKENEWPHASE,
    … …  
    537537        self.PeakEdit = wx.Menu(title='')
    538538        self.PeakMenu.Append(menu=self.PeakEdit, title='Peak Fitting')
     539        self.AutoSearch = self.PeakEdit.Append(help='Automatic peak search',
     540            id=wxID_AUTOSEARCH, kind=wx.ITEM_NORMAL,text='Auto search')
    539541        self.UnDo = self.PeakEdit.Append(help='Undo last least squares refinement',
    540542            id=wxID_UNDO, kind=wx.ITEM_NORMAL,text='UnDo')

  • trunk/GSASIImath.py

    r779 r795  
    780780        CEsig = np.std(CErho)
    781781        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
     782        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem! make 20. adjustible
    782783        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
    783784        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
    … …  
    980981    return Ind
    981982   
     983def setPeakparms(Parms,pos,mag,ifQ=False):
     984    ins = {}
     985    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
     986        for x in ['U','V','W','X','Y']:
     987            ins[x] = Parms[x][1]
     988        if ifQ:                              #qplot - convert back to 2-theta
     989            pos = 2.0*asind(pos*wave/(4*math.pi))
     990        sig = ins['U']*tand(pos/2.0)**2+ins['V']*tand(pos/2.0)+ins['W']
     991        gam = ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)           
     992        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
     993    else:
     994        dsp = pos/Parms['difC'][1]
     995        if 'Pdabc' in Parms:
     996            for x in ['var-inst','X','Y']:
     997                ins[x] = Parms[x][1]
     998            Pdabc = Parms['Pdabc'].T
     999            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
     1000            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
     1001        else:
     1002            for x in ['alpha','beta-0','beta-0','var-inst','X','Y']:
     1003                ins[x] = Parms[x][1]
     1004            alp = ins['alpha']*dsp
     1005            bet = ins['beta-0']+ins['beta-0']/dsp**4
     1006        sig = ins['var-inst']*dsp**2
     1007        gam = ins['X']*dsp+ins['Y']*dsp**2
     1008        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
     1009    return XY
     1010       
     1011def getWave(Parms):
     1012    try:
     1013        return Parms['Lam'][1]
     1014    except KeyError:
     1015        return Parms['Lam1'][1]
     1016   
    9821017def prodQQ(QA,QB):
    9831018    ''' Grassman quaternion product

  • trunk/GSASIIphsGUI.py

    r789 r795  
    403403        if 'Flip' not in generalData:
    404404            generalData['Flip'] = {'RefList':'','Resolution':0.5,'Norm element':'None',
    405                 'k-factor':0.1}
     405                'k-factor':0.1,'k-Max':20.}
    406406        if 'doPawley' not in generalData:
    407407            generalData['doPawley'] = False
    … …  
    856856               
    857857        def FlipSizer():
     858            if 'k-Max' not in Flip: Flip['k-Max'] = 20.
    858859           
    859860            def OnRefList(event):
    … …  
    885886                kFactor.SetValue("%.3f"%(Flip['k-factor']))          #reset in case of error
    886887           
     888            def OnkMax(event):
     889                try:
     890                    res = float(kMax.GetValue())
     891                    if res >= 10.:
     892                        Flip['k-Max'] = res
     893                except ValueError:
     894                    pass
     895                kMax.SetValue("%.0f"%(Flip['k-Max']))          #reset in case of error
     896
    887897            refList = data['Histograms'].keys()
    888898            flipSizer = wx.BoxSizer(wx.VERTICAL)
    … …  
    893903            refList.Bind(wx.EVT_COMBOBOX,OnRefList)
    894904            lineSizer.Add(refList,0,wx.ALIGN_CENTER_VERTICAL)
     905            lineSizer.Add(wx.StaticText(dataDisplay,label=' Normalizing element: '),0,wx.ALIGN_CENTER_VERTICAL)
     906            normElem = wx.Button(dataDisplay,label=Flip['Norm element'],style=wx.TE_READONLY)
     907            normElem.Bind(wx.EVT_BUTTON,OnNormElem)
     908            lineSizer.Add(normElem,0,wx.ALIGN_CENTER_VERTICAL)
    895909            flipSizer.Add(lineSizer,0,wx.ALIGN_CENTER_VERTICAL)
    896910            line2Sizer = wx.BoxSizer(wx.HORIZONTAL)
    897             line2Sizer.Add(wx.StaticText(dataDisplay,label=' Normalizing element: '),0,wx.ALIGN_CENTER_VERTICAL)
    898             normElem = wx.Button(dataDisplay,label=Flip['Norm element'],style=wx.TE_READONLY)
    899             normElem.Bind(wx.EVT_BUTTON,OnNormElem)
    900             line2Sizer.Add(normElem,0,wx.ALIGN_CENTER_VERTICAL)
    901911            line2Sizer.Add(wx.StaticText(dataDisplay,label=' Resolution: '),0,wx.ALIGN_CENTER_VERTICAL)
    902912            flipRes =  wx.TextCtrl(dataDisplay,value='%.2f'%(Flip['Resolution']),style=wx.TE_PROCESS_ENTER)
    … …  
    909919            kFactor.Bind(wx.EVT_KILL_FOCUS,OnkFactor)
    910920            line2Sizer.Add(kFactor,0,wx.ALIGN_CENTER_VERTICAL)
     921            line2Sizer.Add(wx.StaticText(dataDisplay,label=' k-Max (<10.0): '),0,wx.ALIGN_CENTER_VERTICAL)
     922            kMax = wx.TextCtrl(dataDisplay,value='%.0f'%(Flip['k-Max']),style=wx.TE_PROCESS_ENTER)
     923            kMax.Bind(wx.EVT_TEXT_ENTER,OnkMax)       
     924            kMax.Bind(wx.EVT_KILL_FOCUS,OnkMax)
     925            line2Sizer.Add(kMax,0,wx.ALIGN_CENTER_VERTICAL)
    911926            flipSizer.Add(line2Sizer,0,wx.ALIGN_CENTER_VERTICAL)
    912927            return flipSizer
    … …  
    964979                G2frame.dataFrame.AtomsMenu.Enable(item,False)
    965980        Items = [G2gd.wxID_ATOMVIEWINSERT, G2gd.wxID_ATOMSVIEWADD]
    966         if data['Drawing']['showABC']:
     981        if 'showABC' in data['Drawing']:
    967982            for item in Items:
    968983                G2frame.dataFrame.AtomsMenu.Enable(item,True)
    … …  
    16501665                else:
    16511666                    X = G2spc.ApplyStringOps(opr,SGData,atom[3:6])
    1652                     atomInfo = [atom[:2]+list(X)+oldatom[5:9]+atom[9:]+[oldatom[-1]]][0]
     1667                    atomInfo = [atom[:2]+list(X)+oldatom[5:9]+atom[9:]+oldatom[17:]][0]
    16531668            else:
    16541669                atomInfo = [atom[:2]+atom[3:6]+['1',]+['vdW balls',]+
    … …  
    40194034        xdata = G2frame.PatternTree.GetItemPyData(PatternId)[1]
    40204035        Inst = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Instrument Parameters'))
    4021         Inst = dict(zip(Inst[3],Inst[1]))
    40224036        Sample = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Sample Parameters'))
    4023         if 'Lam' in Inst:
    4024             wave = Inst['Lam']
    4025         else:
    4026             wave = Inst['Lam1']
    4027         posCorr = Inst['Zero']
     4037        wave = G2mth.getWave(Inst)
     4038        posCorr = Inst['Zero'][1]
    40284039        const = 9.e-2/(np.pi*Sample['Gonio. radius'])                  #shifts in microns
    40294040       
    … …  
    40434054                    ref[6] = FWHM*(xdata[1][indx]-xdata[4][indx])*cosd(pos/2.)**3/dx
    40444055                    Lorenz = 1./(2.*sind(xdata[0][indx]/2.)**2*cosd(xdata[0][indx]/2.))           #Lorentz correction
    4045                     pola,dIdPola = G2pwd.Polarization(Inst['Polariz.'],xdata[0][indx],Inst['Azimuth'])
     4056                    pola,dIdPola = G2pwd.Polarization(Inst['Polariz.'][1],xdata[0][indx],Inst['Azimuth'][1])
    40464057                    ref[6] /= (Lorenz*pola*ref[3])
    40474058                except IndexError:

  • trunk/GSASIIplot.py

    r792 r795  
    323323    global HKL
    324324
    325     def getWave(Parms):
    326         try:
    327             return Parms['Lam'][1]
    328         except KeyError:
    329             return Parms['Lam1'][1]
    330    
    331325    def OnKeyBox(event):
    332326        if G2frame.G2plotNB.nb.GetSelection() == G2frame.G2plotNB.plotList.index('Powder Patterns'):
    … …  
    436430                Parms = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Instrument Parameters'))
    437431                if 'C' in Parms['Type'][0]:
    438                     wave = getWave(Parms)
     432                    wave = G2mth.getWave(Parms)
    439433                    if G2frame.qPlot:
    440434                        try:
    … …  
    484478            return
    485479        if 'C' in Parms['Type'][0]:
    486             wave = getWave(Parms)
     480            wave = G2mth.getWave(Parms)
    487481        else:
    488482            difC = Parms['difC'][1]
    … …  
    497491            if ind.all() != [0]:                                    #picked a data point
    498492                data = G2frame.PatternTree.GetItemPyData(G2frame.PickId)
    499                 if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
    500                     ins = [Parms[x][1] for x in ['U','V','W','X','Y']]
    501                     if G2frame.qPlot:                              #qplot - convert back to 2-theta
    502                         xy[0] = 2.0*asind(xy[0]*wave/(4*math.pi))
    503                     sig = ins[0]*tand(xy[0]/2.0)**2+ins[1]*tand(xy[0]/2.0)+ins[2]
    504                     gam = ins[3]/cosd(xy[0]/2.0)+ins[4]*tand(xy[0]/2.0)           
    505                     XY = [xy[0],0, xy[1],1, sig,0, gam,0]       #default refine intensity 1st
    506                 else:
    507                     XY = [xy[0],0,xy[1],1,1,0,1,0]
     493                XY = G2mth.setPeakparms(Parms,xy[0],xy[1])
    508494                data.append(XY)
    509495                G2pdG.UpdatePeakGrid(G2frame,data)
    … …  
    536522        Parms = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Instrument Parameters'))
    537523        if 'C' in Parms['Type'][0]:
    538             wave = getWave(Parms)
     524            wave = G2mth.getWave(Parms)
    539525        else:
    540526            difC = Parms['difC'][1]
    … …  
    690676        Parms = ParmList[N]
    691677        if 'C' in Parms['Type'][0]:
    692             wave = getWave(Parms)
     678            wave = G2mth.getWave(Parms)
    693679        else:
    694680            difC = Parms['difC'][1]
    … …  
    781767        Parms = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Instrument Parameters'))
    782768        if 'C' in Parms['Type'][0]:
    783             wave = getWave(Parms)
     769            wave = G2mth.getWave(Parms)
    784770        else:
    785771            difC = Parms['difC'][1]

  • trunk/GSASIIpwd.py

    r794 r795  
    2828import GSASIIgrid as G2gd
    2929import GSASIIIO as G2IO
     30import GSASIImath as G2mth
    3031import pypowder as pyd
    3132
    … …  
    195196    return sumNoAtoms/Vol
    196197           
    197 def MultGetQ(Tth,MuT,Geometry,b=88.0,a=0.01):
    198     NS = 500
    199     Gama = np.linspace(0.,np.pi/2.,NS,False)[1:]
    200     Cgama = np.cos(Gama)[:,np.newaxis]
    201     Sgama = np.sin(Gama)[:,np.newaxis]
    202     CSgama = 1.0/Sgama
    203     Delt = Gama[1]-Gama[0]
    204     emc = 7.94e-26
    205     Navo = 6.023e23
    206     Cth = npcosd(Tth/2.0)
    207     CTth = npcosd(Tth)
    208     Sth = npcosd(Tth/2.0)
    209     STth = npsind(Tth)
    210     CSth = 1.0/Sth
    211     CSTth = 1.0/STth
    212     SCth = 1.0/Cth
    213     SCTth = 1.0/CTth
    214     if 'Bragg' in Geometry:
    215         G = 8.0*Delt*Navo*emc*Sth/((1.0-CTth**2)*(1.0-np.exp(-2.0*MuT*CSth)))
    216         Y1 = np.pi
    217         Y2 = np.pi/2.0
    218         Y3 = 3.*np.pi/8. #3pi/4?
    219         W = 1.0/(Sth+np.fabs(Sgama))
    220         W += np.exp(-MuT*CSth)*(2.0*np.fabs(Sgama)*np.exp(-MuT*np.fabs(CSgama))-
    221             (Sth+np.fabs(Sgama))*np.exp(-MuT*CSth))/(Sth**2-Sgama**2)
    222         Fac0 = Sth**2*Sgama**2
    223         X = Fac0+(Fac0+CTth)**2/2
    224         Y = Cgama**2*Cth**2*(1.0-Fac0-CTth)
    225         Z = Cgama**4*Cth**4/2.0
    226         E = 2.0*(1.0-a)/(b*Cgama/Cth)
    227         F1 = (2.0+b*(1.0+Sth*Sgama))/(b*Cth*Cgama) #trouble if < 1
    228         F2 = (2.0+b*(1.0-Sth*Sgama))/(b*Cth*Cgama)
    229         T1 = np.pi/np.sqrt(F1**2-1.0)
    230         T2 = np.pi/np.sqrt(F2**2-1.0)
    231         Y4 = T1+T2
    232         Y5 = F1**2*T1+F2**2*T2-np.pi*(F1+F2)
    233         Y6 = F1**4*T1+F2**4*T2-np.pi*(F1+F2)/2.0-np.pi*(F1**3+F2**3)
    234         Y7 = (T2-T1)/(F1-F2)
    235         YT = F2**2*T2-F1**2*T1
    236         Y8 = Y1+YT/(F1-F2)
    237         Y9 = Y2+(F2**4*T2-F1**4*T1)/(F1-F2)+Y1*((F1+F2)**2-F1*F2)
    238         M = (a**2*(X*Y1+Y*Y2+Z*Y3)+a*E*(X*Y4+Y*Y5+Z*Y6)+E**2*(X*Y7+Y*Y8+Z*Y9))*Cgama
    239        
    240         Q = np.sum(W*M,axis=0)
    241         return Q*G*NS/(NS-1.)
    242 #
    243 #      cos2th=2.0d*costh^2 - 1.0d
    244 #      G= delta * 8.0d * Na * emc * sinth/(1.0d + cos2th^2)/(1.0d - exp(-2.0d*mut*cscth))
    245 #      Y1=3.1415926d
    246 #      Y2=Y1*0.5d
    247 #      Y3=Y2*0.75d
    248 #      for i=1,num_steps-1 do begin
    249 #         cosgama=double(cos(gama[i]))
    250 #         singama=double(sin(gama[i]))
    251 #         cscgama=1.0d / singama
    252 #
    253 #         W=1.0d/(sinth+abs(singama))
    254 #         W=W+exp(-1.0*mut*cscth)*(2.0d*abs(singama)*exp(-1.0d*mut*abs(cscgama))- $
    255 #                 (sinth+abs(singama))*exp(-1.0d*mut*cscth))/(sinth^2-singama^2)
    256 #
    257 #         factor0=sinth^2*singama^2
    258 #         X=factor0+(factor0+cos2th)^2/2.0d
    259 #         Y=cosgama^2*(1.0d - factor0-cos2th)*costh^2
    260 #         Z=cosgama^4/2.0d*costh^4
    261 #         E=2.0d*(1.0-a)/b/cosgama/costh
    262 #
    263 #         F1=1.0d/b/cosgama*(2.0d + b*(1.0+sinth*singama))/costh
    264 #         F2=1.0d/b/cosgama*(2.0d + b*(1.0-sinth*singama))/costh
    265 #
    266 #         T1=3.14159/sqrt(F1^2-1.0d)
    267 #         T2=3.14159/sqrt(F2^2-1.0d)
    268 #         Y4=T1+T2
    269 #         Y5=F1^2*T1+F2^2*T2-3.14159*(F1+F2)
    270 #         Y6=F1^4*T1+F2^4*T2-3.14159*(F1+F2)/2.0-3.14159*(F1^3+F2^3)
    271 #         Y7=(T2-T1)/(F1-F2)
    272 #         Y8=Y1+(F2^2*T2-F1^2*T1)/(F1-F2)
    273 #         Y9=Y2+(F2^4*T2-F1^4*T1)/(F1-F2)+Y1*((F1+F2)^2-F1*F2)
    274 #         M=(a^2*(X*Y1+Y*Y2+Z*Y3)+a*E*(X*Y4+Y*Y5+Z*Y6)+E^2* $
    275 #                      (X*Y7+Y*Y8+Z*Y9))*cosgama
    276 #
    277 #         Q=Q+W*M
    278 #
    279 #      endfor
    280 #      Q=double(num_steps)/Double(num_steps-1)*Q*G
    281 #      end
    282     elif 'Cylinder' in Geometry:
    283         Q = 0.
    284     elif 'Fixed' in Geometry:   #Dwiggens & Park, Acta Cryst. A27, 264 (1971) with corrections
    285         EMA = np.exp(-MuT*(1.0-SCTth))
    286         Fac1 = (1.-EMA)/(1.0-SCTth)
    287         G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*Fac1)
    288         Fac0 = Cgama/(1-Sgama*SCTth)
    289         Wp = Fac0*(Fac1-(EMA-np.exp(-MuT*(CSgama-SCTth)))/(CSgama-1.0))
    290         Fac0 = Cgama/(1.0+Sgama*SCTth)
    291         Wm = Fac0*(Fac1+(np.exp(-MuT*(1.0+CSgama))-1.0)/(CSgama+1.0))
    292         X = (Sgama**2+CTth**2*(1.0-Sgama**2+Sgama**4))/2.0
    293         Y = Sgama**3*Cgama*STth*CTth
    294         Z = Cgama**2*(1.0+Sgama**2)*STth**2/2.0
    295         Fac2 = 1.0/(b*Cgama*STth)
    296         U = 2.0*(1.0-a)*Fac2
    297         V = (2.0+b*(1.0-CTth*Sgama))*Fac2
    298         Mp = 2.0*np.pi*(a+2.0*(1.0-a)/(2.0+b*(1.0-Sgama)))*(a*X+a*Z/2.0-U*Y+U*(X+Y*V+Z*V**2)/np.sqrt(V**2-1.0)-U*Z*V)
    299         V = (2.0+b*(1.0+CTth*Sgama))*Fac2
    300         Y = -Y
    301         Mm = 2.0*np.pi*(a+2.0*(1.0-a)/(2.0+b*(1.0+Sgama)))*(a*X+a*Z/2.0-U*Y+U*(X+Y*V+Z*V**2)/np.sqrt(V**2-1.0)-U*Z*V)
    302         Q = np.sum(Wp*Mp+Wm*Mm,axis=0)
    303         return Q*G*NS/(NS-1.)
    304     elif 'Tilting' in Geometry:
    305         EMA = np.exp(-MuT*SCth)
    306         G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*EMA)
    307 #          Wplus = expmutsecth/(1.0d - singama*secth) + singama/mut/(1.0 -singama*secth)/(1.0-singama*secth)* $
    308 #                                                       (Exp(-1.0*mut*cscgama) - expmutsecth)
    309 #          Wminus = expmutsecth/(1.0d + singama*secth) - singama/mut/(1.0 +singama*secth)/(1.0+singama*secth)* $
    310 #                                                        expmutsecth*(1.0d - expmutsecth*Exp(-1.0*mut*cscgama))
    311         Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)/(1.0-Sgama*SCth)*(np.exp(-MuT*CSgama)-EMA)
    312 #        Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)**2*(np.exp(-MuT*CSgama)-EMA)
    313         Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)/(1.0+Sgama*SCth)*EMA*(1.0-EMA*np.exp(-MuT*CSgama))
    314 #        Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)**2*EMA*(1.0-EMA*np.exp(-MuT*CSgama))
    315         X = 0.5*(Cth**2*(Cth**2*Sgama**4-4.0*Sth**2*Cgama**2)+1.0)
    316         Y = Cgama**2*(1.0+Cgama**2)*Cth**2*Sth**2
    317         Z = 0.5*Cgama**4*Sth**4
    318 #          X = 0.5*(costh*costh*(costh*costh*singama*singama*singama*singama - $
    319 #                           4.0*sinth*sinth*cosgama*cosgama) +1.0d)
    320 #
    321 #          Y = cosgama*cosgama*(1.0 + cosgama*cosgama)*costh*costh*sinth*sinth
    322 #
    323 #          Z= 0.5 *cosgama*cosgama*cosgama*cosgama* (sinth^4)
    324 #
    325         AlP = 2.0+b*(1.0-Cth*Sgama)
    326         AlM = 2.0+b*(1.0+Cth*Sgama)
    327 #          alphaplus = 2.0 + b*(1.0 - costh*singama)
    328 #          alphaminus = 2.0 + b*(1.0 + costh*singama)
    329         BeP = np.sqrt(np.fabs(AlP**2-(b*Cgama*Sth)**2))
    330         BeM = np.sqrt(np.fabs(AlM**2-(b*Cgama*Sth)**2))
    331 #          betaplus = Sqrt(Abs(alphaplus*alphaplus - b*b*cosgama*cosgama*sinth*sinth))
    332 #          betaminus = Sqrt(Abs(alphaminus*alphaminus - b*b*cosgama*cosgama*sinth*sinth))
    333         Mp = Cgama*(np.pi*a**2*(2.0*X+Y+0.75*Z)+(2.0*np.pi*(1.0-a))*(1.0-a+a*AlP)* \
    334             (4.0*X/AlP/BeP+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlP/BeP-1.0)+
    335             2.0/b**4*AlP/BeP*AlP**2-2.0/b**4*AlP**2-Cgama**2/b/b*Sth*2))
    336 #          Mplus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + $
    337 #                   (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaplus)* $
    338 #                   (4.0*x/alphaplus/betaplus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaplus/betaplus -1.0) + $
    339 #                   2.0/(b^4)*alphaplus/betaplus*alphaplus*alphaplus - 2.0/(b^4)*alphaplus*alphaplus - $
    340 #                   cosgama*cosgama/b/b*sinth*sinth))
    341         Mm =Cgama*(np.pi*a**2*(2.0*X+Y+0.75*Z)+(2.0*np.pi*(1.0-a))*(1.0-a+a*AlM)* \
    342             (4.0*X/AlM/BeM+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlM/BeM-1.0)+
    343             2.0/b**4*AlM/BeM*AlM**2-2.0/b**4*AlM**2-Cgama**2/b/b*Sth*2))
    344 #          Mminus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + $
    345 #                   (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaminus)* $
    346 #                   (4.0*x/alphaminus/betaminus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaminus/betaminus -1.0) + $
    347 #                   2.0/(b^4)*alphaminus/betaminus*alphaminus*alphaminus - 2.0/(b^4)*alphaminus*alphaminus - $
    348 #                   cosgama*cosgama/b/b*sinth*sinth))
    349         Q = np.sum(Wp*Mp+Wm*Mm,axis=0)
    350         return Q*G*NS/(NS-1.)
    351 #       expmutsecth = Exp(-1.0*mut*secth)
    352 #       G= delta * 2.0 * Na * emc /(1.0+costth^2)/expmutsecth
    353 #       for i=1, num_steps-1 do begin
    354 #          cosgama=double(cos(gama[i]))
    355 #          singama=double(sin(gama[i]))
    356 #          cscgama=1.0d/singama
    357 #
    358 #
    359 #     ; print, "W", min(wplus), max(wplus), min(wminus), max(wminus)
    360 #
    361 #
    362 #
    363 #
    364 #    ;               print, a,b
    365 #  ; print, "M", min(mplus), max(mplus), min(mminus), max(mminus)
    366 #          Q=Q+ Wplus*Mplus + Wminus*Mminus
    367 #      endfor
    368 #      Q=double(num_steps)/double(num_steps-1)*Q*G
    369 #   ;   print, min(q), max(q)
    370 #      end
    371 
    372 def MultiScattering(Geometry,ElList,Tth):
    373     BN = BD = 0.0
    374     Amu = 0.0
    375     for El in ElList:
    376         el = ElList[El]
    377         BN += el['Z']*el['FormulaNo']
    378         BD += el['FormulaNo']
    379         Amu += el['FormulaNo']*el['mu']
     198#def MultGetQ(Tth,MuT,Geometry,b=88.0,a=0.01):
     199#    NS = 500
     200#    Gama = np.linspace(0.,np.pi/2.,NS,False)[1:]
     201#    Cgama = np.cos(Gama)[:,np.newaxis]
     202#    Sgama = np.sin(Gama)[:,np.newaxis]
     203#    CSgama = 1.0/Sgama
     204#    Delt = Gama[1]-Gama[0]
     205#    emc = 7.94e-26
     206#    Navo = 6.023e23
     207#    Cth = npcosd(Tth/2.0)
     208#    CTth = npcosd(Tth)
     209#    Sth = npcosd(Tth/2.0)
     210#    STth = npsind(Tth)
     211#    CSth = 1.0/Sth
     212#    CSTth = 1.0/STth
     213#    SCth = 1.0/Cth
     214#    SCTth = 1.0/CTth
     215#    if 'Bragg' in Geometry:
     216#        G = 8.0*Delt*Navo*emc*Sth/((1.0-CTth**2)*(1.0-np.exp(-2.0*MuT*CSth)))
     217#        Y1 = np.pi
     218#        Y2 = np.pi/2.0
     219#        Y3 = 3.*np.pi/8. #3pi/4?
     220#        W = 1.0/(Sth+np.fabs(Sgama))
     221#        W += np.exp(-MuT*CSth)*(2.0*np.fabs(Sgama)*np.exp(-MuT*np.fabs(CSgama))-
     222#            (Sth+np.fabs(Sgama))*np.exp(-MuT*CSth))/(Sth**2-Sgama**2)
     223#        Fac0 = Sth**2*Sgama**2
     224#        X = Fac0+(Fac0+CTth)**2/2
     225#        Y = Cgama**2*Cth**2*(1.0-Fac0-CTth)
     226#        Z = Cgama**4*Cth**4/2.0
     227#        E = 2.0*(1.0-a)/(b*Cgama/Cth)
     228#        F1 = (2.0+b*(1.0+Sth*Sgama))/(b*Cth*Cgama) #trouble if < 1
     229#        F2 = (2.0+b*(1.0-Sth*Sgama))/(b*Cth*Cgama)
     230#        T1 = np.pi/np.sqrt(F1**2-1.0)
     231#        T2 = np.pi/np.sqrt(F2**2-1.0)
     232#        Y4 = T1+T2
     233#        Y5 = F1**2*T1+F2**2*T2-np.pi*(F1+F2)
     234#        Y6 = F1**4*T1+F2**4*T2-np.pi*(F1+F2)/2.0-np.pi*(F1**3+F2**3)
     235#        Y7 = (T2-T1)/(F1-F2)
     236#        YT = F2**2*T2-F1**2*T1
     237#        Y8 = Y1+YT/(F1-F2)
     238#        Y9 = Y2+(F2**4*T2-F1**4*T1)/(F1-F2)+Y1*((F1+F2)**2-F1*F2)
     239#        M = (a**2*(X*Y1+Y*Y2+Z*Y3)+a*E*(X*Y4+Y*Y5+Z*Y6)+E**2*(X*Y7+Y*Y8+Z*Y9))*Cgama
     240#       
     241#        Q = np.sum(W*M,axis=0)
     242#        return Q*G*NS/(NS-1.)
     243##
     244##      cos2th=2.0d*costh^2 - 1.0d
     245##      G= delta * 8.0d * Na * emc * sinth/(1.0d + cos2th^2)/(1.0d - exp(-2.0d*mut*cscth))
     246##      Y1=3.1415926d
     247##      Y2=Y1*0.5d
     248##      Y3=Y2*0.75d
     249##      for i=1,num_steps-1 do begin
     250##         cosgama=double(cos(gama[i]))
     251##         singama=double(sin(gama[i]))
     252##         cscgama=1.0d / singama
     253##
     254##         W=1.0d/(sinth+abs(singama))
     255##         W=W+exp(-1.0*mut*cscth)*(2.0d*abs(singama)*exp(-1.0d*mut*abs(cscgama))- $
     256##                 (sinth+abs(singama))*exp(-1.0d*mut*cscth))/(sinth^2-singama^2)
     257##
     258##         factor0=sinth^2*singama^2
     259##         X=factor0+(factor0+cos2th)^2/2.0d
     260##         Y=cosgama^2*(1.0d - factor0-cos2th)*costh^2
     261##         Z=cosgama^4/2.0d*costh^4
     262##         E=2.0d*(1.0-a)/b/cosgama/costh
     263##
     264##         F1=1.0d/b/cosgama*(2.0d + b*(1.0+sinth*singama))/costh
     265##         F2=1.0d/b/cosgama*(2.0d + b*(1.0-sinth*singama))/costh
     266##
     267##         T1=3.14159/sqrt(F1^2-1.0d)
     268##         T2=3.14159/sqrt(F2^2-1.0d)
     269##         Y4=T1+T2
     270##         Y5=F1^2*T1+F2^2*T2-3.14159*(F1+F2)
     271##         Y6=F1^4*T1+F2^4*T2-3.14159*(F1+F2)/2.0-3.14159*(F1^3+F2^3)
     272##         Y7=(T2-T1)/(F1-F2)
     273##         Y8=Y1+(F2^2*T2-F1^2*T1)/(F1-F2)
     274##         Y9=Y2+(F2^4*T2-F1^4*T1)/(F1-F2)+Y1*((F1+F2)^2-F1*F2)
     275##         M=(a^2*(X*Y1+Y*Y2+Z*Y3)+a*E*(X*Y4+Y*Y5+Z*Y6)+E^2* $
     276##                      (X*Y7+Y*Y8+Z*Y9))*cosgama
     277##
     278##         Q=Q+W*M
     279##
     280##      endfor
     281##      Q=double(num_steps)/Double(num_steps-1)*Q*G
     282##      end
     283#    elif 'Cylinder' in Geometry:
     284#        Q = 0.
     285#    elif 'Fixed' in Geometry:   #Dwiggens & Park, Acta Cryst. A27, 264 (1971) with corrections
     286#        EMA = np.exp(-MuT*(1.0-SCTth))
     287#        Fac1 = (1.-EMA)/(1.0-SCTth)
     288#        G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*Fac1)
     289#        Fac0 = Cgama/(1-Sgama*SCTth)
     290#        Wp = Fac0*(Fac1-(EMA-np.exp(-MuT*(CSgama-SCTth)))/(CSgama-1.0))
     291#        Fac0 = Cgama/(1.0+Sgama*SCTth)
     292#        Wm = Fac0*(Fac1+(np.exp(-MuT*(1.0+CSgama))-1.0)/(CSgama+1.0))
     293#        X = (Sgama**2+CTth**2*(1.0-Sgama**2+Sgama**4))/2.0
     294#        Y = Sgama**3*Cgama*STth*CTth
     295#        Z = Cgama**2*(1.0+Sgama**2)*STth**2/2.0
     296#        Fac2 = 1.0/(b*Cgama*STth)
     297#        U = 2.0*(1.0-a)*Fac2
     298#        V = (2.0+b*(1.0-CTth*Sgama))*Fac2
     299#        Mp = 2.0*np.pi*(a+2.0*(1.0-a)/(2.0+b*(1.0-Sgama)))*(a*X+a*Z/2.0-U*Y+U*(X+Y*V+Z*V**2)/np.sqrt(V**2-1.0)-U*Z*V)
     300#        V = (2.0+b*(1.0+CTth*Sgama))*Fac2
     301#        Y = -Y
     302#        Mm = 2.0*np.pi*(a+2.0*(1.0-a)/(2.0+b*(1.0+Sgama)))*(a*X+a*Z/2.0-U*Y+U*(X+Y*V+Z*V**2)/np.sqrt(V**2-1.0)-U*Z*V)
     303#        Q = np.sum(Wp*Mp+Wm*Mm,axis=0)
     304#        return Q*G*NS/(NS-1.)
     305#    elif 'Tilting' in Geometry:
     306#        EMA = np.exp(-MuT*SCth)
     307#        G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*EMA)
     308##          Wplus = expmutsecth/(1.0d - singama*secth) + singama/mut/(1.0 -singama*secth)/(1.0-singama*secth)* $
     309##                                                       (Exp(-1.0*mut*cscgama) - expmutsecth)
     310##          Wminus = expmutsecth/(1.0d + singama*secth) - singama/mut/(1.0 +singama*secth)/(1.0+singama*secth)* $
     311##                                                        expmutsecth*(1.0d - expmutsecth*Exp(-1.0*mut*cscgama))
     312#        Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)/(1.0-Sgama*SCth)*(np.exp(-MuT*CSgama)-EMA)
     313##        Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)**2*(np.exp(-MuT*CSgama)-EMA)
     314#        Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)/(1.0+Sgama*SCth)*EMA*(1.0-EMA*np.exp(-MuT*CSgama))
     315##        Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)**2*EMA*(1.0-EMA*np.exp(-MuT*CSgama))
     316#        X = 0.5*(Cth**2*(Cth**2*Sgama**4-4.0*Sth**2*Cgama**2)+1.0)
     317#        Y = Cgama**2*(1.0+Cgama**2)*Cth**2*Sth**2
     318#        Z = 0.5*Cgama**4*Sth**4
     319##          X = 0.5*(costh*costh*(costh*costh*singama*singama*singama*singama - $
     320##                           4.0*sinth*sinth*cosgama*cosgama) +1.0d)
     321##
     322##          Y = cosgama*cosgama*(1.0 + cosgama*cosgama)*costh*costh*sinth*sinth
     323##
     324##          Z= 0.5 *cosgama*cosgama*cosgama*cosgama* (sinth^4)
     325##
     326#        AlP = 2.0+b*(1.0-Cth*Sgama)
     327#        AlM = 2.0+b*(1.0+Cth*Sgama)
     328##          alphaplus = 2.0 + b*(1.0 - costh*singama)
     329##          alphaminus = 2.0 + b*(1.0 + costh*singama)
     330#        BeP = np.sqrt(np.fabs(AlP**2-(b*Cgama*Sth)**2))
     331#        BeM = np.sqrt(np.fabs(AlM**2-(b*Cgama*Sth)**2))
     332##          betaplus = Sqrt(Abs(alphaplus*alphaplus - b*b*cosgama*cosgama*sinth*sinth))
     333##          betaminus = Sqrt(Abs(alphaminus*alphaminus - b*b*cosgama*cosgama*sinth*sinth))
     334#        Mp = Cgama*(np.pi*a**2*(2.0*X+Y+0.75*Z)+(2.0*np.pi*(1.0-a))*(1.0-a+a*AlP)* \
     335#            (4.0*X/AlP/BeP+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlP/BeP-1.0)+
     336#            2.0/b**4*AlP/BeP*AlP**2-2.0/b**4*AlP**2-Cgama**2/b/b*Sth*2))
     337##          Mplus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + $
     338##                   (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaplus)* $
     339##                   (4.0*x/alphaplus/betaplus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaplus/betaplus -1.0) + $
     340##                   2.0/(b^4)*alphaplus/betaplus*alphaplus*alphaplus - 2.0/(b^4)*alphaplus*alphaplus - $
     341##                   cosgama*cosgama/b/b*sinth*sinth))
     342#        Mm =Cgama*(np.pi*a**2*(2.0*X+Y+0.75*Z)+(2.0*np.pi*(1.0-a))*(1.0-a+a*AlM)* \
     343#            (4.0*X/AlM/BeM+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlM/BeM-1.0)+
     344#            2.0/b**4*AlM/BeM*AlM**2-2.0/b**4*AlM**2-Cgama**2/b/b*Sth*2))
     345##          Mminus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + $
     346##                   (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaminus)* $
     347##                   (4.0*x/alphaminus/betaminus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaminus/betaminus -1.0) + $
     348##                   2.0/(b^4)*alphaminus/betaminus*alphaminus*alphaminus - 2.0/(b^4)*alphaminus*alphaminus - $
     349##                   cosgama*cosgama/b/b*sinth*sinth))
     350#        Q = np.sum(Wp*Mp+Wm*Mm,axis=0)
     351#        return Q*G*NS/(NS-1.)
     352##       expmutsecth = Exp(-1.0*mut*secth)
     353##       G= delta * 2.0 * Na * emc /(1.0+costth^2)/expmutsecth
     354##       for i=1, num_steps-1 do begin
     355##          cosgama=double(cos(gama[i]))
     356##          singama=double(sin(gama[i]))
     357##          cscgama=1.0d/singama
     358##
     359##
     360##     ; print, "W", min(wplus), max(wplus), min(wminus), max(wminus)
     361##
     362##
     363##
     364##
     365##    ;               print, a,b
     366##  ; print, "M", min(mplus), max(mplus), min(mminus), max(mminus)
     367##          Q=Q+ Wplus*Mplus + Wminus*Mminus
     368##      endfor
     369##      Q=double(num_steps)/double(num_steps-1)*Q*G
     370##   ;   print, min(q), max(q)
     371##      end
     372
     373#def MultiScattering(Geometry,ElList,Tth):
     374#    BN = BD = 0.0
     375#    Amu = 0.0
     376#    for El in ElList:
     377#        el = ElList[El]
     378#        BN += el['Z']*el['FormulaNo']
     379#        BD += el['FormulaNo']
     380#        Amu += el['FormulaNo']*el['mu']
    380381       
    381382def CalcPDF(data,inst,xydata):
    … …  
    402403    MuR = Abs*data['Diam']/20.0
    403404    xydata['IofQ'][1][1] /= Absorb(data['Geometry'],MuR,Tth)
    404     xydata['IofQ'][1][1] /= Polarization(inst['Polariz.'],Tth,Azm=inst['Azimuth'])[0]
     405    xydata['IofQ'][1][1] /= Polarization(inst['Polariz.'][1],Tth,Azm=inst['Azimuth'][1])[0]
    405406    if data['DetType'] == 'Image plate':
    406407        xydata['IofQ'][1][1] *= Oblique(data['ObliqCoeff'],Tth)
    … …  
    408409    #convert to Q
    409410    hc = 12.397639
    410     if 'Lam' in inst:
    411         wave = inst['Lam']
    412     else:
    413         wave = inst['Lam1']
     411    wave = G2mth.getWave(inst)
    414412    keV = hc/wave
    415413    minQ = npT2q(Tth[0],wave)
    … …  
    587585fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
    588586               
    589 def getWidths(pos,sig,gam,shl):
     587def getWidthsCW(pos,sig,gam,shl):
    590588    widths = [np.sqrt(sig)/100.,gam/200.]
    591589    fwhm = 2.355*widths[0]+2.*widths[1]
    … …  
    596594    return widths,fmin,fmax
    597595   
     596def getWidthsTOF(pos,alp,bet,sig,gam):
     597    lnf = 3.3      # =log(0.001/2)
     598    widths = [np.sqrt(sig),gam]
     599    fwhm = 2.355*widths[0]+2.*widths[1]
     600    fmin = 10.*fwhm*(1.+1./alp)   
     601    fmax = 10.*fwhm*(1.+1./bet)
     602    return widths,fmin,fmax
     603   
    598604def getFWHM(TTh,Inst):
    599605    sig = lambda Th,U,V,W: 1.17741*math.sqrt(max(0.001,U*tand(Th)**2+V*tand(Th)+W))*math.pi/180.
    600606    gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))*math.pi/180.
    601607    gamFW = lambda s,g: math.exp(math.log(s**5+2.69269*s**4*g+2.42843*s**3*g**2+4.47163*s**2*g**3+0.07842*s*g**4+g**5)/5.)
    602     s = sig(TTh/2.,Inst['U'],Inst['V'],Inst['W'])*100.
    603     g = gam(TTh/2.,Inst['X'],Inst['Y'])*100.
    604     return gamFW(g,s)
    605    
     608    s = sig(TTh/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])*100.
     609    g = gam(TTh/2.,Inst['X'][1],Inst['Y'][1])*100.
     610    return gamFW(g,s)   
    606611               
    607612def getFCJVoigt(pos,intens,sig,gam,shl,xdata):   
    608613    DX = xdata[1]-xdata[0]
    609     widths,fmin,fmax = getWidths(pos,sig,gam,shl)
     614    widths,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
    610615    x = np.linspace(pos-fmin,pos+fmin,256)
    611616    dx = x[1]-x[0]
    … …  
    664669            bakInt = si.interp1d(bakPos,bakVals,'linear')
    665670            yb = bakInt(xdata)
     671    if 'difC' in parmDict:
     672        ff = 1.
     673    else:       
     674        try:
     675            wave = parmDict[pfx+'Lam']
     676        except KeyError:
     677            wave = parmDict[pfx+'Lam1']
     678        q = 4.0*np.pi*npsind(xdata/2.0)/wave
     679        SQ = (q/(4.*np.pi))**2
     680        FF = G2elem.GetFormFactorCoeff('Si')[0]
     681        ff = np.array(G2elem.ScatFac(FF,SQ)[0])**2
    666682    iD = 0       
    667     try:
    668         wave = parmDict[pfx+'Lam']
    669     except KeyError:
    670         wave = parmDict[pfx+'Lam1']
    671     q = 4.0*np.pi*npsind(xdata/2.0)/wave
    672     SQ = (q/(4.*np.pi))**2
    673     FF = G2elem.GetFormFactorCoeff('Si')[0]
    674     ff = np.array(G2elem.ScatFac(FF,SQ)[0])**2
    675683    while True:
    676684        try:
    … …  
    690698            pkG = parmDict[pfx+'BkPkgam:'+str(iD)]
    691699            shl = 0.002
    692             Wd,fmin,fmax = getWidths(pkP,pkS,pkG,shl)
     700            Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,shl)
    693701            iBeg = np.searchsorted(xdata,pkP-fmin)
    694702            iFin = np.searchsorted(xdata,pkP+fmax)
    … …  
    745753                        np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
    746754                        np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
     755    if 'difC' in parmDict:
     756        ff = 1.
     757    else:
     758        try:
     759            wave = parmDict[pfx+'Lam']
     760        except KeyError:
     761            wave = parmDict[pfx+'Lam1']
     762        q = 4.0*np.pi*npsind(xdata/2.0)/wave
     763        SQ = (q/(4*np.pi))**2
     764        FF = G2elem.GetFormFactorCoeff('Si')[0]
     765        ff = np.array(G2elem.ScatFac(FF,SQ)[0])
    747766    iD = 0       
    748     try:
    749         wave = parmDict[pfx+'Lam']
    750     except KeyError:
    751         wave = parmDict[pfx+'Lam1']
    752     q = 4.0*np.pi*npsind(xdata/2.0)/wave
    753     SQ = (q/(4*np.pi))**2
    754     FF = G2elem.GetFormFactorCoeff('Si')[0]
    755     ff = np.array(G2elem.ScatFac(FF,SQ)[0])
    756767    while True:
    757768        try:
    … …  
    776787            pkG = parmDict[pfx+'BkPkgam:'+str(iD)]
    777788            shl = 0.002
    778             Wd,fmin,fmax = getWidths(pkP,pkS,pkG,shl)
     789            Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,shl)
    779790            iBeg = np.searchsorted(xdata,pkP-fmin)
    780791            iFin = np.searchsorted(xdata,pkP+fmax)
    … …  
    804815    return Df,dFdp,dFds,dFdg,dFdsh
    805816
     817def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
     818    Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
     819    Df /= np.sum(Df)
     820    return Df 
     821   
     822def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
     823    Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
     824    sumDf = np.sum(Df)
     825    return Df,dFdp,dFda,dFdb,dFds,dFdg   
     826
    806827def ellipseSize(H,Sij,GB):
    807828    HX = np.inner(H.T,GB)
    … …  
    831852    return HKLs
    832853
    833 def getPeakProfile(parmDict,xdata,varyList,bakType):
     854def getPeakProfile(dataType,parmDict,xdata,varyList,bakType):
    834855   
    835856    yb = getBackground('',parmDict,bakType,xdata)
    836857    yc = np.zeros_like(yb)
    837     dx = xdata[1]-xdata[0]
    838     U = parmDict['U']
    839     V = parmDict['V']
    840     W = parmDict['W']
    841     X = parmDict['X']
    842     Y = parmDict['Y']
    843     shl = max(parmDict['SH/L'],0.002)
    844     Ka2 = False
    845     if 'Lam1' in parmDict.keys():
    846         Ka2 = True
    847         lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
    848         kRatio = parmDict['I(L2)/I(L1)']
    849     iPeak = 0
    850     while True:
    851         try:
    852             pos = parmDict['pos'+str(iPeak)]
    853             theta = (pos-parmDict['Zero'])/2.0
    854             intens = parmDict['int'+str(iPeak)]
    855             sigName = 'sig'+str(iPeak)
    856             if sigName in varyList:
    857                 sig = parmDict[sigName]
    858             else:
    859                 sig = U*tand(theta)**2+V*tand(theta)+W
    860             sig = max(sig,0.001)          #avoid neg sigma
    861             gamName = 'gam'+str(iPeak)
    862             if gamName in varyList:
    863                 gam = parmDict[gamName]
    864             else:
    865                 gam = X/cosd(theta)+Y*tand(theta)
    866             gam = max(gam,0.001)             #avoid neg gamma
    867             Wd,fmin,fmax = getWidths(pos,sig,gam,shl)
    868             iBeg = np.searchsorted(xdata,pos-fmin)
    869             lenX = len(xdata)
    870             if not iBeg:
    871                 iFin = np.searchsorted(xdata,pos+fmin)
    872             elif iBeg == lenX:
    873                 iFin = iBeg
    874             else:
    875                 iFin = min(lenX,iBeg+int((fmin+fmax)/dx))
    876             if not iBeg+iFin:       #peak below low limit
     858    if 'C' in dataType:
     859        dx = xdata[1]-xdata[0]
     860        U = parmDict['U']
     861        V = parmDict['V']
     862        W = parmDict['W']
     863        X = parmDict['X']
     864        Y = parmDict['Y']
     865        shl = max(parmDict['SH/L'],0.002)
     866        Ka2 = False
     867        if 'Lam1' in parmDict.keys():
     868            Ka2 = True
     869            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
     870            kRatio = parmDict['I(L2)/I(L1)']
     871        iPeak = 0
     872        while True:
     873            try:
     874                pos = parmDict['pos'+str(iPeak)]
     875                theta = (pos-parmDict['Zero'])/2.0
     876                intens = parmDict['int'+str(iPeak)]
     877                sigName = 'sig'+str(iPeak)
     878                if sigName in varyList:
     879                    sig = parmDict[sigName]
     880                else:
     881                    sig = U*tand(theta)**2+V*tand(theta)+W
     882                sig = max(sig,0.001)          #avoid neg sigma
     883                gamName = 'gam'+str(iPeak)
     884                if gamName in varyList:
     885                    gam = parmDict[gamName]
     886                else:
     887                    gam = X/cosd(theta)+Y*tand(theta)
     888                gam = max(gam,0.001)             #avoid neg gamma
     889                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
     890                iBeg = np.searchsorted(xdata,pos-fmin)
     891                lenX = len(xdata)
     892                if not iBeg:
     893                    iFin = np.searchsorted(xdata,pos+fmin)
     894                elif iBeg == lenX:
     895                    iFin = iBeg
     896                else:
     897                    iFin = min(lenX,iBeg+int((fmin+fmax)/dx))
     898                if not iBeg+iFin:       #peak below low limit
     899                    iPeak += 1
     900                    continue
     901                elif not iBeg-iFin:     #peak above high limit
     902                    return yb+yc
     903                yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
     904                if Ka2:
     905                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
     906                    kdelt = int((pos2-pos)/dx)               
     907                    iBeg = min(lenX,iBeg+kdelt)
     908                    iFin = min(lenX,iFin+kdelt)
     909                    if iBeg-iFin:
     910                        yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
    877911                iPeak += 1
    878                 continue
    879             elif not iBeg-iFin:     #peak above high limit
     912            except KeyError:        #no more peaks to process
    880913                return yb+yc
    881             yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
    882             if Ka2:
    883                 pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
    884                 kdelt = int((pos2-pos)/dx)               
    885                 iBeg = min(lenX,iBeg+kdelt)
    886                 iFin = min(lenX,iFin+kdelt)
    887                 if iBeg-iFin:
    888                     yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
    889             iPeak += 1
    890         except KeyError:        #no more peaks to process
    891             return yb+yc
     914    else:
     915        difC = parmDict['difC']
     916        alp0 = parmDict['alpha']
     917        bet0 = parmDict['beta-0']
     918        bet1 = parmDict['beta-1']
     919        sig0 = parmDict['var-inst']
     920        X = parmDict['X']
     921        Y = parmDict['Y']
     922        iPeak = 0
     923        while True:
     924            try:
     925                pos = parmDict['pos'+str(iPeak)]               
     926                tof = pos-parmDict['Zero']
     927                dsp = tof/difC
     928                intens = parmDict['int'+str(iPeak)]
     929                alpName = 'alp'+str(iPeak)
     930                if alpName in varyList:
     931                    alp = parmDict[alpName]
     932                else:
     933                    alp = alp0*dsp
     934                betName = 'bet'+str(iPeak)
     935                if betName in varyList:
     936                    bet = parmDict[betName]
     937                else:
     938                    bet = bet0+bet1/dsp**4
     939                sigName = 'sig'+str(iPeak)
     940                if sigName in varyList:
     941                    sig = parmDict[sigName]
     942                else:
     943                    sig = sig0*dsp**2
     944                gamName = 'gam'+str(iPeak)
     945                if gamName in varyList:
     946                    gam = parmDict[gamName]
     947                else:
     948                    gam = X*dsp**2+Y*dsp
     949                gam = max(gam,0.001)             #avoid neg gamma
     950                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
     951                iBeg = np.searchsorted(xdata,pos-fmin)
     952                iFin = np.searchsorted(xdata,pos+fmax)
     953                lenX = len(xdata)
     954                if not iBeg:
     955                    iFin = np.searchsorted(xdata,pos+fmax)
     956                elif iBeg == lenX:
     957                    iFin = iBeg
     958                else:
     959                    iFin = np.searchsorted(xdata,pos+fmax)
     960                if not iBeg+iFin:       #peak below low limit
     961                    iPeak += 1
     962                    continue
     963                elif not iBeg-iFin:     #peak above high limit
     964                    return yb+yc
     965                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
     966                iPeak += 1
     967            except KeyError:        #no more peaks to process
     968                return yb+yc
    892969           
    893 def getPeakProfileDerv(parmDict,xdata,varyList,bakType):
     970def getPeakProfileDerv(dataType,parmDict,xdata,varyList,bakType):
    894971# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
    895972    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
    … …  
    909986            ip = names.index(parm)
    910987            dMdv[varyList.index(name)] = dMdpk[4*int(id)+ip]
    911     dx = xdata[1]-xdata[0]
    912     U = parmDict['U']
    913     V = parmDict['V']
    914     W = parmDict['W']
    915     X = parmDict['X']
    916     Y = parmDict['Y']
    917     shl = max(parmDict['SH/L'],0.002)
    918     Ka2 = False
    919     if 'Lam1' in parmDict.keys():
    920         Ka2 = True
    921         lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
    922         kRatio = parmDict['I(L2)/I(L1)']
    923     iPeak = 0
    924     while True:
    925         try:
    926             pos = parmDict['pos'+str(iPeak)]
    927             theta = (pos-parmDict['Zero'])/2.0
    928             intens = parmDict['int'+str(iPeak)]
    929             sigName = 'sig'+str(iPeak)
    930             tanth = tand(theta)
    931             costh = cosd(theta)
    932             if sigName in varyList:
    933                 sig = parmDict[sigName]
    934             else:
    935                 sig = U*tanth**2+V*tanth+W
    936                 dsdU = tanth**2
    937                 dsdV = tanth
    938                 dsdW = 1.0
    939             sig = max(sig,0.001)          #avoid neg sigma
    940             gamName = 'gam'+str(iPeak)
    941             if gamName in varyList:
    942                 gam = parmDict[gamName]
    943             else:
    944                 gam = X/costh+Y*tanth
    945                 dgdX = 1.0/costh
    946                 dgdY = tanth
    947             gam = max(gam,0.001)             #avoid neg gamma
    948             Wd,fmin,fmax = getWidths(pos,sig,gam,shl)
    949             iBeg = np.searchsorted(xdata,pos-fmin)
    950             lenX = len(xdata)
    951             if not iBeg:
    952                 iFin = np.searchsorted(xdata,pos+fmin)
    953             elif iBeg == lenX:
    954                 iFin = iBeg
    955             else:
    956                 iFin = min(lenX,iBeg+int((fmin+fmax)/dx))
    957             if not iBeg+iFin:       #peak below low limit
     988    if 'C' in dataType:
     989        dx = xdata[1]-xdata[0]
     990        U = parmDict['U']
     991        V = parmDict['V']
     992        W = parmDict['W']
     993        X = parmDict['X']
     994        Y = parmDict['Y']
     995        shl = max(parmDict['SH/L'],0.002)
     996        Ka2 = False
     997        if 'Lam1' in parmDict.keys():
     998            Ka2 = True
     999            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
     1000            kRatio = parmDict['I(L2)/I(L1)']
     1001        iPeak = 0
     1002        while True:
     1003            try:
     1004                pos = parmDict['pos'+str(iPeak)]
     1005                theta = (pos-parmDict['Zero'])/2.0
     1006                intens = parmDict['int'+str(iPeak)]
     1007                sigName = 'sig'+str(iPeak)
     1008                tanth = tand(theta)
     1009                costh = cosd(theta)
     1010                if sigName in varyList:
     1011                    sig = parmDict[sigName]
     1012                else:
     1013                    sig = U*tanth**2+V*tanth+W
     1014                    dsdU = tanth**2
     1015                    dsdV = tanth
     1016                    dsdW = 1.0
     1017                sig = max(sig,0.001)          #avoid neg sigma
     1018                gamName = 'gam'+str(iPeak)
     1019                if gamName in varyList:
     1020                    gam = parmDict[gamName]
     1021                else:
     1022                    gam = X/costh+Y*tanth
     1023                    dgdX = 1.0/costh
     1024                    dgdY = tanth
     1025                gam = max(gam,0.001)             #avoid neg gamma
     1026                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
     1027                iBeg = np.searchsorted(xdata,pos-fmin)
     1028                lenX = len(xdata)
     1029                if not iBeg:
     1030                    iFin = np.searchsorted(xdata,pos+fmin)
     1031                elif iBeg == lenX:
     1032                    iFin = iBeg
     1033                else:
     1034                    iFin = min(lenX,iBeg+int((fmin+fmax)/dx))
     1035                if not iBeg+iFin:       #peak below low limit
     1036                    iPeak += 1
     1037                    continue
     1038                elif not iBeg-iFin:     #peak above high limit
     1039                    break
     1040                dMdpk = np.zeros(shape=(6,len(xdata)))
     1041                dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
     1042                for i in range(1,5):
     1043                    dMdpk[i][iBeg:iFin] += 100.*dx*intens*dMdipk[i]
     1044                dMdpk[0][iBeg:iFin] += 100.*dx*dMdipk[0]
     1045                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
     1046                if Ka2:
     1047                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
     1048                    kdelt = int((pos2-pos)/dx)               
     1049                    iBeg = min(lenX,iBeg+kdelt)
     1050                    iFin = min(lenX,iFin+kdelt)
     1051                    if iBeg-iFin:
     1052                        dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
     1053                        for i in range(1,5):
     1054                            dMdpk[i][iBeg:iFin] += 100.*dx*intens*kRatio*dMdipk2[i]
     1055                        dMdpk[0][iBeg:iFin] += 100.*dx*kRatio*dMdipk2[0]
     1056                        dMdpk[5][iBeg:iFin] += 100.*dx*dMdipk2[0]
     1057                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
     1058                for parmName in ['pos','int','sig','gam']:
     1059                    try:
     1060                        idx = varyList.index(parmName+str(iPeak))
     1061                        dMdv[idx] = dervDict[parmName]
     1062                    except ValueError:
     1063                        pass
     1064                if 'U' in varyList:
     1065                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
     1066                if 'V' in varyList:
     1067                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
     1068                if 'W' in varyList:
     1069                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
     1070                if 'X' in varyList:
     1071                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
     1072                if 'Y' in varyList:
     1073                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
     1074                if 'SH/L' in varyList:
     1075                    dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
     1076                if 'I(L2)/I(L1)' in varyList:
     1077                    dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
    9581078                iPeak += 1
    959                 continue
    960             elif not iBeg-iFin:     #peak above high limit
     1079            except KeyError:        #no more peaks to process
    9611080                break
    962             dMdpk = np.zeros(shape=(6,len(xdata)))
    963             dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
    964             for i in range(1,5):
    965                 dMdpk[i][iBeg:iFin] += 100.*dx*intens*dMdipk[i]
    966             dMdpk[0][iBeg:iFin] += 100.*dx*dMdipk[0]
    967             dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
    968             if Ka2:
    969                 pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
    970                 kdelt = int((pos2-pos)/dx)               
    971                 iBeg = min(lenX,iBeg+kdelt)
    972                 iFin = min(lenX,iFin+kdelt)
    973                 if iBeg-iFin:
    974                     dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
    975                     for i in range(1,5):
    976                         dMdpk[i][iBeg:iFin] += 100.*dx*intens*kRatio*dMdipk2[i]
    977                     dMdpk[0][iBeg:iFin] += 100.*dx*kRatio*dMdipk2[0]
    978                     dMdpk[5][iBeg:iFin] += 100.*dx*dMdipk2[0]
    979                     dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
    980             for parmName in ['pos','int','sig','gam']:
    981                 try:
    982                     idx = varyList.index(parmName+str(iPeak))
    983                     dMdv[idx] = dervDict[parmName]
    984                 except ValueError:
    985                     pass
    986             if 'U' in varyList:
    987                 dMdv[varyList.index('U')] += dsdU*dervDict['sig']
    988             if 'V' in varyList:
    989                 dMdv[varyList.index('V')] += dsdV*dervDict['sig']
    990             if 'W' in varyList:
    991                 dMdv[varyList.index('W')] += dsdW*dervDict['sig']
    992             if 'X' in varyList:
    993                 dMdv[varyList.index('X')] += dgdX*dervDict['gam']
    994             if 'Y' in varyList:
    995                 dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
    996             if 'SH/L' in varyList:
    997                 dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
    998             if 'I(L2)/I(L1)' in varyList:
    999                 dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
    1000             iPeak += 1
    1001         except KeyError:        #no more peaks to process
    1002             break
     1081    else:
     1082        difC = parmDict['difC']
     1083        alp0 = parmDict['alpha']
     1084        bet0 = parmDict['beta-0']
     1085        bet1 = parmDict['beta-1']
     1086        sig0 = parmDict['var-inst']
     1087        X = parmDict['X']
     1088        Y = parmDict['Y']
     1089        iPeak = 0
     1090        while True:
     1091            try:
     1092                pos = parmDict['pos'+str(iPeak)]               
     1093                tof = pos-parmDict['Zero']
     1094                dsp = tof/difC
     1095                intens = parmDict['int'+str(iPeak)]
     1096                alpName = 'alp'+str(iPeak)
     1097                if alpName in varyList:
     1098                    alp = parmDict[alpName]
     1099                else:
     1100                    alp = alp0*dsp
     1101                    dada0 = dsp
     1102                betName = 'bet'+str(iPeak)
     1103                if betName in varyList:
     1104                    bet = parmDict[betName]
     1105                else:
     1106                    bet = bet0+bet1/dsp**4
     1107                    dbdb0 = 1.
     1108                    dbdb1 = 1/dsp**4
     1109                sigName = 'sig'+str(iPeak)
     1110                if sigName in varyList:
     1111                    sig = parmDict[sigName]
     1112                else:
     1113                    sig = sig0*dsp**2
     1114                    dsds0 = dsp**2
     1115                gamName = 'gam'+str(iPeak)
     1116                if gamName in varyList:
     1117                    gam = parmDict[gamName]
     1118                else:
     1119                    gam = X*dsp**2+Y*dsp
     1120                    dsdX = dsp**2
     1121                    dsdY = dsp
     1122                gam = max(gam,0.001)             #avoid neg gamma
     1123                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
     1124                iBeg = np.searchsorted(xdata,pos-fmin)
     1125                lenX = len(xdata)
     1126                if not iBeg:
     1127                    iFin = np.searchsorted(xdata,pos+fmax)
     1128                elif iBeg == lenX:
     1129                    iFin = iBeg
     1130                else:
     1131                    iFin = np.searchsorted(xdata,pos+fmax)
     1132                if not iBeg+iFin:       #peak below low limit
     1133                    iPeak += 1
     1134                    continue
     1135                elif not iBeg-iFin:     #peak above high limit
     1136                    break
     1137                dMdpk = np.zeros(shape=(7,len(xdata)))
     1138                dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
     1139                for i in range(1,6):
     1140                    dMdpk[i][iBeg:iFin] += intens*dMdipk[i]
     1141                dMdpk[0][iBeg:iFin] += dMdipk[0]
     1142                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
     1143                for parmName in ['pos','int','alp','bet','sig','gam']:
     1144                    try:
     1145                        idx = varyList.index(parmName+str(iPeak))
     1146                        dMdv[idx] = dervDict[parmName]
     1147                    except ValueError:
     1148                        pass
     1149                if 'alpha' in varyList:
     1150                    dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
     1151                if 'beta-0' in varyList:
     1152                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
     1153                if 'beta-1' in varyList:
     1154                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
     1155                if 'var-inst' in varyList:
     1156                    dMdv[varyList.index('var-inst')] += dsds0*dervDict['sig']
     1157                if 'X' in varyList:
     1158                    dMdv[varyList.index('X')] += dsdX*dervDict['gam']
     1159                if 'Y' in varyList:
     1160                    dMdv[varyList.index('Y')] += dsdY*dervDict['gam']         #problem here
     1161                iPeak += 1
     1162            except KeyError:        #no more peaks to process
     1163                break
    10031164    return dMdv
    10041165       
    … …  
    11431304            insNames.append(parm)
    11441305            insVals.append(Inst[parm][1])
    1145             if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)'] and Inst[parm][2]:
    1146                 insVary.append(parm)
     1306            if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
     1307                'beta-0','beta-1','var-inst',] and Inst[parm][2]:
     1308                    insVary.append(parm)
    11471309        instDict = dict(zip(insNames,insVals))
    11481310        instDict['X'] = max(instDict['X'],0.01)
    11491311        instDict['Y'] = max(instDict['Y'],0.01)
    1150         instDict['SH/L'] = max(instDict['SH/L'],0.002)
     1312        if 'SH/L' in instDict:
     1313            instDict['SH/L'] = max(instDict['SH/L'],0.002)
    11511314        return dataType,instDict,insVary
    11521315       
    … …  
    11601323                pos = parmDict['pos'+str(iPeak)]
    11611324                if sigName not in varyList:
    1162                     parmDict[sigName] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
     1325                    if 'C' in Inst['Type'][0]:
     1326                        parmDict[sigName] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
     1327                    else:
     1328                        dsp = pos/Inst['difC'][1]
     1329                        parmDict[sigName] = parmDict['var-inst']*dsp**2
    11631330                gamName = 'gam'+str(iPeak)
    11641331                if gamName not in varyList:
    1165                     parmDict[gamName] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
     1332                    if 'C' in Inst['Type'][0]:
     1333                        parmDict[gamName] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
     1334                    else:
     1335                        dsp = pos/Inst['difC'][1]
     1336                        parmDict[sigName] = parmDict['X']*dsp**2+parmDict['Y']*dsp
    11661337                iPeak += 1
    11671338            except KeyError:
    … …  
    11751346        sigstr = 'esds  :'
    11761347        for parm in Inst:
    1177             if parm in  ['U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
     1348            if parm in  ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
     1349                'beta-0','beta-1','var-inst',]:
    11781350                ptlbls += "%s" % (parm.center(12))
    11791351                ptstr += ptfmt % (Inst[parm][1])
    … …  
    11861358        print sigstr
    11871359
    1188     def SetPeaksParms(Peaks):
     1360    def SetPeaksParms(dataType,Peaks):
    11891361        peakNames = []
    11901362        peakVary = []
    11911363        peakVals = []
    1192         names = ['pos','int','sig','gam']
     1364        if 'C' in dataType:
     1365            names = ['pos','int','sig','gam']
     1366        else:
     1367            names = ['pos','int','alpha','beta','sig','gam']
    11931368        for i,peak in enumerate(Peaks):
    1194             for j in range(4):
     1369            for j,name in enumerate(names):
    11951370                peakVals.append(peak[2*j])
    1196                 parName = names[j]+str(i)
     1371                parName = name+str(i)
    11971372                peakNames.append(parName)
    11981373                if peak[2*j+1]:
    … …  
    12001375        return dict(zip(peakNames,peakVals)),peakVary
    12011376               
    1202     def GetPeaksParms(parmDict,Peaks,varyList):
    1203         names = ['pos','int','sig','gam']
     1377    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
     1378        if 'C' in Inst['Type'][0]:
     1379            names = ['pos','int','sig','gam']
     1380        else:
     1381            names = ['pos','int','alpha','beta','sig','gam']
    12041382        for i,peak in enumerate(Peaks):
    1205             for j in range(4):
    1206                 pos = parmDict['pos'+str(i)]
     1383            pos = parmDict['pos'+str(i)]
     1384            if 'difC' in Inst:
     1385                dsp = pos/Inst['difC'][1]
     1386            for j in range(len(names)):
    12071387                parName = names[j]+str(i)
    12081388                if parName in varyList:
    12091389                    peak[2*j] = parmDict[parName]
     1390                elif 'alpha' in parName:
     1391                    peak[2*j] = parmDict['alpha']*dsp
     1392                elif 'beta' in parName:
     1393                    peak[2*j] = parmDict['beta-0']+parmDict['beta-1']/dsp**4
    12101394                elif 'sig' in parName:
    1211                     peak[2*j] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
     1395                    if 'C' in Inst['Type'][0]:
     1396                        peak[2*j] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
     1397                    else:
     1398                        peak[2*j] = parmDict['var-inst']*dsp**2
    12121399                elif 'gam' in parName:
    1213                     peak[2*j] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
     1400                    if 'C' in Inst['Type'][0]:
     1401                        peak[2*j] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
     1402                    else:
     1403                        peak[2*j] = parmDict['X']*dsp**2+parmDict['Y']*dsp
    12141404                       
    1215     def PeaksPrint(parmDict,sigDict,varyList):
     1405    def PeaksPrint(dataType,parmDict,sigDict,varyList):
    12161406        print 'Peak coefficients:'
    1217         names = ['pos','int','sig','gam']
    1218         head = 15*' '
     1407        if 'C' in dataType:
     1408            names = ['pos','int','sig','gam']
     1409        else:
     1410            names = ['pos','int','alpha','beta','sig','gam']           
     1411        head = 13*' '
    12191412        for name in names:
    1220             head += name.center(12)+'esd'.center(12)
     1413            head += name.center(10)+'esd'.center(10)
    12211414        print head
    1222         ptfmt = {'pos':"%12.5f",'int':"%12.1f",'sig':"%12.3f",'gam':"%12.3f"}
     1415        if 'C' in dataType:
     1416            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
     1417        else:
     1418            ptfmt = {'pos':"%10.2f",'int':"%10.1f",'alpha':"%10.3f",'beta':"%10.5f",'sig':"%10.3f",'gam':"%10.3f"}
    12231419        for i,peak in enumerate(Peaks):
    12241420            ptstr =  ':'
    1225             for j in range(4):
     1421            for j in range(len(names)):
    12261422                name = names[j]
    12271423                parName = name+str(i)
    … …  
    12321428                else:
    12331429#                    ptstr += G2IO.ValEsd(parmDict[parName],0.0)
    1234                     ptstr += 12*' '
     1430                    ptstr += 10*' '
    12351431            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
    12361432               
    1237     def devPeakProfile(values, xdata, ydata, weights, parmdict, varylist,bakType,dlg):
     1433    def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
    12381434        parmdict.update(zip(varylist,values))
    1239         return np.sqrt(weights)*getPeakProfileDerv(parmdict,xdata,varylist,bakType)
     1435        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
    12401436           
    1241     def errPeakProfile(values, xdata, ydata, weights, parmdict, varylist,bakType,dlg):       
     1437    def errPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):       
    12421438        parmdict.update(zip(varylist,values))
    1243         M = np.sqrt(weights)*(getPeakProfile(parmdict,xdata,varylist,bakType)-ydata)
     1439        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
    12441440        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
    12451441        if dlg:
    … …  
    12621458    bakType,bakDict,bakVary = SetBackgroundParms(Background)
    12631459    dataType,insDict,insVary = SetInstParms(Inst)
    1264     peakDict,peakVary = SetPeaksParms(Peaks)
     1460    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
    12651461    parmDict = {}
    12661462    parmDict.update(bakDict)
    … …  
    12741470            try:
    12751471                result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
    1276                     args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg))
     1472                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
    12771473                ncyc = int(result[2]['nfev']/2)
    12781474            finally:
    … …  
    13051501    sigDict = dict(zip(varyList,sig))
    13061502    yb[xBeg:xFin] = getBackground('',parmDict,bakType,x[xBeg:xFin])
    1307     yc[xBeg:xFin] = getPeakProfile(parmDict,x[xBeg:xFin],varyList,bakType)
     1503    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)
    13081504    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
    13091505    GetBackgroundParms(parmDict,Background)
    … …  
    13111507    GetInstParms(parmDict,Inst,varyList,Peaks)
    13121508    InstPrint(Inst,sigDict)
    1313     GetPeaksParms(parmDict,Peaks,varyList)   
    1314     PeaksPrint(parmDict,sigDict,varyList)
     1509    GetPeaksParms(Inst,parmDict,Peaks,varyList)   
     1510    PeaksPrint(dataType,parmDict,sigDict,varyList)
     1511
     1512def calcIncident(Iparm,xdata):
     1513    Itype = Iparm['Itype']
     1514    Icoef = Iparm['Icoeff']
     1515    Iesd = Iparm['Iesd']
     1516    Icovar = Iparm['Icovar']
     1517    intens = np.zeros_like(xdata)
     1518    x = xdata/1000.
     1519    if Itype == 1:
     1520        intens = Icoef[0]
     1521        for i in range(1,10,2):
     1522            intens += Icoef[i]*np.exp(-Icoef[i+1]*x**((i+1)/2))
     1523       
     1524    elif Itype == 4:
     1525        intens = Icoef[0]
     1526        intens += (Icoef[1]/x**5)*np.exp(-Icoef[2]/(x**2))
     1527       
     1528    return intens
     1529
    13151530   
    13161531#testing data

  • trunk/GSASIIpwdGUI.py

    r794 r795  
    1313import wx.grid as wg
    1414import numpy as np
     15import numpy.ma as ma
    1516import math
    1617import time
    … …  
    2021import GSASIIpath
    2122GSASIIpath.SetVersionNumber("$Revision$")
     23import GSASIImath as G2mth
    2224import GSASIIpwd as G2pwd
    2325import GSASIIIO as G2IO
    … …  
    2931import GSASIIElemGUI as G2elemGUI
    3032import GSASIIElem as G2elem
    31 
    3233VERY_LIGHT_GREY = wx.Colour(235,235,235)
    33 
    3434# trig functions in degrees
    3535sind = lambda x: math.sin(x*math.pi/180.)
    … …  
    5858        'Temperature':300.,'Pressure':1.0,'Humidity':0.0,
    5959        'Voltage':0.0,'Force':0.0,'Gonio. radius':200.0,
    60         'Omega':0.0,'Chi':0.0,'Phi':0.0}   
    61            
     60        'Omega':0.0,'Chi':0.0,'Phi':0.0}
     61                         
    6262################################################################################
    6363#####  Powder Peaks
    … …  
    6767    if G2frame.dataDisplay:
    6868        G2frame.dataFrame.Clear()
     69       
     70    def OnAutoSearch(event):
     71        PatternId = G2frame.PatternId
     72        PickId = G2frame.PickId
     73        limits = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Limits'))[1]
     74        inst = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Instrument Parameters'))
     75        profile = G2frame.PatternTree.GetItemPyData(PatternId)[1]
     76        x0 = profile[0]
     77        iBeg = np.searchsorted(x0,limits[0])
     78        iFin = np.searchsorted(x0,limits[1])
     79        x = x0[iBeg:iFin]
     80        y0 = profile[1][iBeg:iFin]
     81        y1 = copy.copy(y0)
     82        ysig = np.std(y1)
     83        offset = [-1,1]
     84        ymask = ma.array(y0,mask=(y0<ysig))
     85        for off in offset:
     86            ymask = ma.array(ymask,mask=(ymask-np.roll(y0,off)<=0.))
     87        indx = ymask.nonzero()
     88        mags = ymask[indx]
     89        poss = x[indx]
     90        for pos,mag in zip(poss,mags):
     91            data.append(G2mth.setPeakparms(inst,pos,mag))
     92        UpdatePeakGrid(G2frame,data)
     93        G2plt.PlotPatterns(G2frame)       
    6994   
    7095    def OnUnDo(event):
    … …  
    149174        PatternId = G2frame.PatternId
    150175        PickId = G2frame.PickId
     176        Inst = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Instrument Parameters'))
    151177        peaks = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Peak List'))
    152178        if not peaks:
    153179            G2frame.ErrorDialog('No peaks!','Nothing to do!')
    154180            return
    155         Inst = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Instrument Parameters'))
    156181        for peak in peaks:
    157182            if 'C' in Inst['Type'][0]:
    … …  
    159184                peak[6] = Inst['X'][1]/cosd(peak[0]/2.0)+Inst['Y'][1]*tand(peak[0]/2.0)
    160185            else:
    161                 pass
     186                dsp = peak[0]/Inst['difC'][0]
     187                peak[4] = Inst['alpha'][0]*dsp
     188                peak[6] = Inst['beta-0'][0]+Inst['beta-1'][0]/dsp**4
     189                peak[8] = Inst['var-inst'][0]*dsp**2
     190                peak[10] = Inst['X'][0]*dsp**2+Inst['Y'][0]*dsp
    162191        UpdatePeakGrid(G2frame,peaks)
    163192               
    … …  
    178207       for r in range(G2frame.dataDisplay.GetNumberRows()):
    179208           for c in range(G2frame.dataDisplay.GetNumberCols()):
    180                if G2frame.dataDisplay.GetColLabelValue(c) in ['position','intensity','sigma','gamma']:
     209               if G2frame.dataDisplay.GetColLabelValue(c) in ['position','intensity','alpha','beta','sigma','gamma']:
    181210                   if float(G2frame.dataDisplay.GetCellValue(r,c)) < 0.:
    182211                       G2frame.dataDisplay.SetCellBackgroundColour(r,c,wx.RED)
    … …  
    239268    if not G2frame.dataFrame.GetStatusBar():
    240269        Status = G2frame.dataFrame.CreateStatusBar()
     270    G2frame.Bind(wx.EVT_MENU, OnAutoSearch, id=G2gd.wxID_AUTOSEARCH)
    241271    G2frame.Bind(wx.EVT_MENU, OnUnDo, id=G2gd.wxID_UNDO)
    242272    G2frame.Bind(wx.EVT_MENU, OnLSQPeakFit, id=G2gd.wxID_LSQPEAKFIT)
    … …  
    250280    G2frame.PickTable = []
    251281    rowLabels = []
     282    PatternId = G2frame.PatternId
     283    Inst = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Instrument Parameters'))
    252284    for i in range(len(data)): rowLabels.append(str(i+1))
    253     colLabels = ['position','refine','intensity','refine','sigma','refine','gamma','refine']
    254     Types = [wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_BOOL,
    255         wg.GRID_VALUE_FLOAT+':10,1',wg.GRID_VALUE_BOOL,
    256         wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_BOOL,
    257         wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_BOOL]
     285    if 'C' in Inst['Type'][0]:
     286        colLabels = ['position','refine','intensity','refine','sigma','refine','gamma','refine']
     287        Types = [wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_BOOL,
     288            wg.GRID_VALUE_FLOAT+':10,1',wg.GRID_VALUE_BOOL,
     289            wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_BOOL,
     290            wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_BOOL]
     291    else:
     292        colLabels = ['position','refine','intensity','refine','alpha','refine',
     293            'beta','refine','sigma','refine','gamma','refine']
     294        Types = [wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_BOOL,
     295            wg.GRID_VALUE_FLOAT+':10,1',wg.GRID_VALUE_BOOL,
     296            wg.GRID_VALUE_FLOAT+':10,4',wg.GRID_VALUE_BOOL,
     297            wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_BOOL,
     298            wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_BOOL,
     299            wg.GRID_VALUE_FLOAT+':10,5',wg.GRID_VALUE_BOOL]
    258300    T = []
    259301    for peak in data:
    … …  
    642684                    peak[6] = insVal['X']/cosd(peak[0]/2.0)+insVal['Y']*tand(peak[0]/2.0)
    643685            else:
    644                 pass
    645 #                for peak in peaks:
    646                    
     686                for peak in peaks:
     687                    dsp = peak[0]/insVal['difC']
     688                    peak[4] = insVal['alpha']*dsp
     689                    peak[6] = insVal['beta-0']+insVal['beta-1']/dsp**4
     690                    peak[8] = insVal['var-inst']*dsp**2
     691                    peak[10] = insVal['X']*dsp**2+insVal['Y']*dsp
    647692                   
    648693    def OnLoad(event):
    … …  
    673718                    S = File.readline()               
    674719                File.close()
    675                 data = dict(zip(newItems,zip(newVals,newVals,len(newVals)*[False,])))
    676                 for item in data:
    677                     data[item] = list(data[item])
     720                inst = G2IO.makeInstDict(newItems,newVals,len(newVals)*[False,])
    678721                G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId,'Instrument Parameters'),data)
    679722                RefreshInstrumentGrid(event,doAnyway=True)          #to get peaks updated
    … …  
    940983        else:                                   #time of flight (neutrons)
    941984            instSizer.Add(wx.StaticText(G2frame.dataDisplay,-1,' Azimuth: %7.2f'%(insVal['Azimuth'])),0,wx.ALIGN_CENTER_VERTICAL)
     985            instSizer.Add(wx.StaticText(G2frame.dataDisplay,-1,' 2-theta: %7.2f'%(insVal['2-theta'])),0,wx.ALIGN_CENTER_VERTICAL)
     986            instSizer.Add((5,5),0)
     987            instSizer.Add((5,5),0)
    942988            for item in ['difC','difA','Zero','alpha','beta-0','beta-1','var-inst','X','Y']:
    943989                fmt = '%10.3f'
    … …  
    952998                itemVal.Bind(wx.EVT_KILL_FOCUS,OnItemValue)
    953999                instSizer.Add(itemVal,0,wx.ALIGN_CENTER_VERTICAL)
    954                 itemRef = wx.CheckBox(G2frame.dataDisplay,wx.ID_ANY,label=' Refine?')
    955                 itemRef.SetValue(bool(insRef[item]))
    956                 RefObj[itemRef.GetId()] = item
    957                 itemRef.Bind(wx.EVT_CHECKBOX, OnItemRef)
    958                 instSizer.Add(itemRef,0,wx.ALIGN_CENTER_VERTICAL)
     1000                if not ifHisto and item in ['difC','difA','Zero',]:
     1001                    instSizer.Add((5,5),0)
     1002                else:
     1003                    itemRef = wx.CheckBox(G2frame.dataDisplay,wx.ID_ANY,label=' Refine?')
     1004                    itemRef.SetValue(bool(insRef[item]))
     1005                    RefObj[itemRef.GetId()] = item
     1006                    itemRef.Bind(wx.EVT_CHECKBOX, OnItemRef)
     1007                    instSizer.Add(itemRef,0,wx.ALIGN_CENTER_VERTICAL)
    9591008       
    9601009    else:                       #single crystal data
    … …  
    12641313    IndexId = G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Index Peak List')
    12651314    Inst = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Instrument Parameters'))
    1266     try:
    1267         wave = Inst['Lam'][1]
    1268     except KeyError:
    1269         wave = Inst['Lam1'][0]
     1315    wave = G2mth.getWave(Inst)
    12701316   
    12711317    def RefreshIndexPeaksGrid(event):
    … …  
    13691415        'P 4/m m m','F m m m','I m m m','C m m m','P m m m','C 2/m','P 2/m','P -1']
    13701416    Inst = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Instrument Parameters'))
    1371     if 'Lam' in Inst:
    1372         wave = Inst['Lam'][1]
    1373     else:
    1374         wave = Inst['Lam1'][0]
    1375        
     1417    wave = G2mth.getWave(Inst)
     1418   
    13761419    def SetLattice(controls):
    13771420        ibrav = bravaisSymb.index(controls[5])

  • trunk/GSASIIstruct.py

    r792 r795  
    16971697            controlDict[pfx+'wtFactor'] =Histogram['wtFactor']
    16981698            Inst = Histogram['Instrument Parameters']
    1699             controlDict[pfx+'histType'] = Inst[1][0]
    1700             histDict[pfx+'Lam'] = Inst[1][1]
     1699            controlDict[pfx+'histType'] = Inst['Type'][0]
     1700            histDict[pfx+'Lam'] = Inst['Lam'][1]
    17011701            controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']                   
    17021702    return histVary,histDict,controlDict
    … …  
    25702570                    if 'C' in calcControls[hfx+'histType']:
    25712571                        yp = np.zeros_like(yb)
    2572                         Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
     2572                        Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
    25732573                        iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
    25742574                        iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
    … …  
    25772577                        if Ka2:
    25782578                            pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
    2579                             Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
     2579                            Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
    25802580                            iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
    25812581                            iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
    … …  
    26652665#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
    26662666                        continue
    2667                 Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
     2667                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
    26682668                iBeg = np.searchsorted(x,refl[5]-fmin)
    26692669                iFin = np.searchsorted(x,refl[5]+fmax)
    … …  
    26752675                if Ka2:
    26762676                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
    2677                     Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
     2677                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
    26782678                    iBeg = np.searchsorted(x,pos2-fmin)
    26792679                    iFin = np.searchsorted(x,pos2+fmax)
    … …  
    27732773                h,k,l = refl[:3]
    27742774                dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
    2775                 Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
     2775                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
    27762776                iBeg = np.searchsorted(x,refl[5]-fmin)
    27772777                iFin = np.searchsorted(x,refl[5]+fmax)

  • trunk/fsource/pypowder.for

    r768 r795  
    1717      REAL*4 TTHETA,SIG,GAM,SPH
    1818      INTEGER*4 NPTS,I
    19       FW = (2.355*SQRT(SIG)+GAM)/100.0
    20       FMIN = 10.0*(-FW-SPH*COSD(TTHETA))
    21       FMAX = 15.0*FW
    2219      DO I=0,NPTS-1
    2320        CALL PSVFCJ(DTT(I)*100.,TTHETA*100.,SIG,GAM,SPH,
    … …  
    5552      REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),SIGPART(0:NPTS-1),
    5653     1  GAMPART(0:NPTS-1),SLPART(0:NPTS-1),PRFUNC(0:NPTS-1)
    57       FW = (2.355*SQRT(SIG)+GAM)/100.0
    58       FMIN = 10.0*(-FW-SPH*COSD(TTHETA))
    59       FMAX = 15.0*FW
    6054      DO I=0,NPTS-1
    6155        CALL PSVFCJ(DTT(I)*100.,TTHETA*100.,SIG,GAM,SPH,
    … …  
    8478      REAL*4 TTHETA,SIG,GAM,SPH
    8579      REAL*4 DTT(0:NPTS-1),PRFUNC(0:NPTS-1)
    86       FW = (2.355*SQRT(SIG)+GAM)/100.0
    87       FMIN = 10.0*(-FW-SPH*COSD(TTHETA))
    88       FMAX = 15.0*FW
    8980      DO I=0,NPTS-1
    9081        CALL PSVFCJO(DTT(I)*100.,TTHETA*100.,SIG,GAM,SPH/2.0,SPH/2.0,
    … …  
    122113      REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),SIGPART(0:NPTS-1),
    123114     1  GAMPART(0:NPTS-1),SLPART(0:NPTS-1),PRFUNC(0:NPTS-1)
    124       FW = (2.355*SQRT(SIG)+GAM)/100.0
    125       FMIN = 10.0*(-FW-SPH*COSD(TTHETA))
    126       FMAX = 15.0*FW
    127115      DO I=0,NPTS-1
    128116        CALL PSVFCJO(DTT(I)*100.,TTHETA*100.,SIG,GAM,SHL/2.,SHL/2.,
    … …  
    130118          SLPART(I) = SPART
    131119        DPRDT(I) = DPRDT(I)*100.
     120      END DO
     121      RETURN
     122      END
     123
     124      SUBROUTINE PYEPSVOIGT(NPTS,DTT,ALP,BET,SIG,GAM,PRFUNC)
     125C DTT in microsec
     126C RETURNS FUNCTION
     127Cf2py intent(in) NPTS
     128Cf2py intent(in) DTT
     129cf2py depend(NPTS) DTT
     130Cf2py intent(in) ALP
     131Cf2py intent(in) BET
     132Cf2py intent(in) SIG
     133Cf2py intent(in) GAM
     134Cf2py intent(out) PRFUNC
     135Cf2py depend(NPTS) PRFUNC
     136
     137      INTEGER*4 NPTS
     138      REAL*4 ALP,BET,SIG,GAM,SHL
     139      REAL*4 DTT(0:NPTS-1),PRFUNC(0:NPTS-1),DPRDT,ALPPART,
     140     1  BETPART,SIGPART,GAMPART
     141      DO I=0,NPTS-1
     142        CALL EPSVOIGT(DTT(I),ALP,BET,SIG,GAM,PRFUNC(I),DPRDT,
     143     1    ALPPART,BETPART,SIGPART,GAMPART)
     144      END DO
     145      RETURN
     146      END
     147
     148      SUBROUTINE PYDEPSVOIGT(NPTS,DTT,ALP,BET,SIG,GAM,PRFUNC,
     149     1  DPRDT,ALPPART,BETPART,SIGPART,GAMPART)
     150C DTT in microsec
     151C RETURNS FUNCTION & DERIVATIVES
     152Cf2py intent(in) NPTS
     153Cf2py intent(in) DTT
     154cf2py depend(NPTS) DTT
     155Cf2py intent(in) ALP
     156Cf2py intent(in) BET
     157Cf2py intent(in) SIG
     158Cf2py intent(in) GAM
     159Cf2py intent(out) PRFUNC
     160Cf2py depend(NPTS) PRFUNC
     161Cf2py intent(out) DPRDT
     162Cf2py depend(NPTS) DPRDT
     163Cf2py intent(out) ALPPART
     164Cf2py depend(NPTS) ALPPART
     165Cf2py intent(out) BETPART
     166Cf2py depend(NPTS) BETPART
     167Cf2py intent(out) SIGPART
     168Cf2py depend(NPTS) SIGPART
     169Cf2py intent(out) GAMPART
     170Cf2py depend(NPTS) GAMPART
     171
     172      INTEGER*4 NPTS
     173      REAL*4 ALP,BET,SIG,GAM,SHL
     174      REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),ALPPART(0:NPTS-1),
     175     1  BETPART(0:NPTS-1),SIGPART(0:NPTS-1),
     176     1  GAMPART(0:NPTS-1),PRFUNC(0:NPTS-1)
     177      DO I=0,NPTS-1
     178        CALL EPSVOIGT(DTT(I),ALP,BET,SIG,GAM,PRFUNC(I),DPRDT(I),
     179     1    ALPPART(I),BETPART(I),SIGPART(I),GAMPART(I))
    132180      END DO
    133181      RETURN

  • trunk/imports/G2pwd_fxye.py

    r792 r795  
    2323            longFormatName = 'GSAS powder data files'
    2424            )
     25        self.clockWd = None
    2526
    2627    # Validate the contents -- look for a bank line
    … …  
    2829        #print 'ContentsValidator: '+self.formatName
    2930        for i,line in enumerate(filepointer):
    30             self.TimeMap = False
     31            self.GSAS = True
    3132            if i==0: # first line is always a comment
    3233                continue
    … …  
    3738            if line[:4] == 'BANK':
    3839                return True
    39             if line [:8] == 'TIME_MAP':          #LANSCE TOF data
    40                 self.TimeMap = True
     40            elif line[:7] == 'Monitor': continue
     41            elif line [:8] == 'TIME_MAP':          #LANSCE TOF data
    4142                return True
    4243            else:
    … …  
    4849
    4950    def Reader(self,filename,filepointer, ParentFrame=None, **kwarg):
    50         clockWd = None
    5151       
    5252        def GetFXYEdata(File,Pos,Bank):
    … …  
    9393            File.seek(Pos)
    9494            cons = Bank.split()
    95             if clockWd:
     95            if self.clockWd:
    9696                start = 0
    9797                step = 1
    … …  
    119119                S = File.readline()
    120120            N = len(x)
    121             if clockWd:
    122                 x,cw = Tmap2TOF(self.TimeMap,clockWd)
    123                 x = x[:-2]+cw[:-1]/2.
    124                 return [x,np.array(y)/cw[:-1],np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
    125             else:
    126                 return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
     121            if self.clockWd:
     122                x = Tmap2TOF(self.TimeMap,clockWd)
     123            return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
    127124       
    128125        def GetSTDdata(File,Pos,Bank):
    … …  
    130127            cons = Bank.split()
    131128            Nch = int(cons[2])
    132             if clockWd:
     129            if self.clockWd:
    133130                start = 0
    134131                step = 1
    135132            else:
    136133                start = float(cons[5])/100.0               #CW: from centidegrees to degrees
    137                 step = float(cons[6])/100.0
     134                step = float(cons[6])/100.0                 #NB TOF 0.1*ms!
    138135            x = []
    139136            y = []
    … …  
    158155                S = File.readline()
    159156            N = len(x)
    160             if clockWd:
    161                 x,cw = Tmap2TOF(self.TimeMap,clockWd)
    162                 x = x[:-2]+cw[:-1]/2.
    163                 return [x,np.array(y)/cw[:-1],np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
    164             else:
    165                 return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
     157            if self.clockWd:
     158                x = Tmap2TOF(self.TimeMap,self.clockWd)[:-2]
     159            return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
    166160           
    167161        def GetTimeMap(File,Pos,TimeMap):
    … …  
    196190                Tch,T,Step = tmap
    197191            TOF = np.array(TOF)*clockWd
    198             chWdt = np.diff(TOF)
    199             return TOF,chWdt
     192            return TOF
    200193
    201194        x = []
    … …  
    244237                        Pos.append(filepointer.tell())
    245238                    if S[:8] == 'TIME_MAP':
    246                         self.TimeMap,clockWd = GetTimeMap(filepointer,filepointer.tell(),S)
     239                        self.TimeMap,self.clockWd = GetTimeMap(filepointer,filepointer.tell(),S)
    247240            except Exception as detail:
    248241                print self.formatName+' scan error:'+str(detail) # for testing
    … …  
    286279        try:
    287280            if 'FXYE' in Bank:
    288                 self.powderdata = GetFXYEdata(filepointer,
    289                                               Pos[selblk],Banks[selblk])
     281                self.powderdata = GetFXYEdata(filepointer,Pos[selblk],Banks[selblk])
    290282            elif 'FXY' in Bank:
    291                 self.powderdata = GetFXYdata(filepointer,
    292                                              Pos[selblk],Banks[selblk])
     283                self.powderdata = GetFXYdata(filepointer,Pos[selblk],Banks[selblk])
    293284            elif 'ESD' in Bank:
    294                 self.powderdata = GetESDdata(filepointer,
    295                                              Pos[selblk],Banks[selblk])
     285                self.powderdata = GetESDdata(filepointer,Pos[selblk],Banks[selblk])
    296286            elif 'STD' in Bank:
    297                 self.powderdata = GetSTDdata(filepointer,
    298                                              Pos[selblk],Banks[selblk])
     287                self.powderdata = GetSTDdata(filepointer,Pos[selblk],Banks[selblk])
    299288            else:
    300                 self.powderdata = GetSTDdata(filepointer,
    301                                              Pos[selblk],Banks[selblk])
     289                self.powderdata = GetSTDdata(filepointer,Pos[selblk],Banks[selblk])
    302290        except Exception as detail:
    303291            print self.formatName+' read error:'+str(detail) # for testing
Note: See TracChangeset for help on using the changeset viewer.