Changeset 5114


Ignore:
Timestamp:
Dec 21, 2021 8:53:00 AM (22 months ago)
Author:
vondreele
Message:

various fixes for EDX data fitting by LeBail?
various fixes for SAD analysis - plotting updates & importer
removed a couple of unused lambdas from G2sasd

Location:
trunk
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIImpsubs.py

    r4919 r5114  
    131131    return sInt,resList
    132132
     133def ComputeFobsSqEDbatch(profList):
     134    sInt = 0
     135    resList = []
     136    for refl,iref in profList:
     137        icod = ComputeFobsSqED(refl,iref)
     138        if type(icod) is tuple:
     139            resList.append((icod[0],iref))
     140            sInt += icod[1]
     141        elif icod == -1:
     142            resList.append((None,iref))
     143        elif icod == -2:
     144            break
     145    return sInt,resList
     146
    133147def ComputeFobsSqCW(refl,iref):
    134148    yp = np.zeros(len(x)) # not masked
    … …  
    197211    return refl8im,refl[11+im]*refl[9+im]
    198212
     213def ComputeFobsSqED(refl,iref):
     214    yp = np.zeros(len(x)) # not masked
     215    refl8im = 0
     216    Wd,fmin,fmax = G2pwd.getWidthsED(refl[5+im],refl[6+im])
     217    iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
     218    iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
     219    if not iBeg+iFin:       #peak below low limit - skip peak
     220        return 0
     221    if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
     222        return -1
     223    elif not iBeg-iFin:     #peak above high limit - done
     224        return -2
     225    if iBeg < iFin:
     226        fp,sumfp = G2pwd.getPsVoigt(refl[5+im],refl[6+im]*1.e4,0.001,x[iBeg:iFin])
     227        yp[iBeg:iFin] = 100.*refl[9+im]*fp*cw[iBeg:iFin]/sumfp
     228    refl8im = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin],0.0))
     229    return refl8im,refl[9+im]
    199230################################################################################
    200231# Powder Profile computation
    … …  
    232263        yc[iBeg:iFin] += refl[11+im]*refl[9+im]*fp
    233264    return yc
     265
     266def ComputePwdrProfED(profList):
     267    'Compute the peaks profile for a set of TOF peaks and add into the yc array'
     268    for pos,refl,iBeg,iFin in profList:
     269        fp = G2pwd.getPsVoigt(pos,refl[6+im]*1.e4,0.001,x[iBeg:iFin])[0]
     270        yc[iBeg:iFin] += refl[9+im]*fp
     271    return yc

  • trunk/GSASIIphsGUI.py

    r5112 r5114  
    39113911            if delList: delList += ', '
    39123912            delList += data['Atoms'][i][0]
    3913         dlg = wx.MessageDialog(G2frame,
    3914             'Do you want to delete atom(s): {}?'.format(delList),
    3915                     'Confirm delete',
    3916                     wx.YES|wx.NO)
     3913        dlg = wx.MessageDialog(G2frame,'Do you want to delete atom(s): {}?'.format(delList),
     3914            'Confirm delete',wx.YES|wx.NO)
    39173915        try:
    39183916            dlg.CenterOnParent()
    … …  
    62526250            print('build of fullrmc file {} completed'.format(rname))
    62536251        elif G2frame.RMCchoice == 'RMCProfile':
    6254             pName = generalData['Name'].replace(' ','_')
     6252            if ' ' in generalData['Name']:
     6253                wx.MessageDialog(G2frame,'ERROR: Phase name has space; change phase name','Bad phase name',wx.ICON_ERROR).ShowModal()
     6254                G2frame.dataWindow.FRMCDataEdit.Enable(G2G.wxID_RUNRMC,False)
     6255                return
     6256            pName = generalData['Name']
    62556257            G2frame.dataWindow.FRMCDataEdit.Enable(G2G.wxID_RUNRMC,True)
    62566258            RMCPdict = data['RMC']['RMCProfile']
    … …  
    62956297                G2frame.dataWindow.FRMCDataEdit.Enable(G2G.wxID_RUNRMC,False)
    62966298        elif G2frame.RMCchoice == 'PDFfit':
    6297             pName = generalData['Name'].replace(' ','_')
     6299            if ' ' in generalData['Name']:
     6300                wx.MessageDialog(G2frame,'ERROR: Phase name has space; change phase name','Bad phase name',wx.ICON_ERROR).ShowModal()
     6301                G2frame.dataWindow.FRMCDataEdit.Enable(G2G.wxID_RUNRMC,False)
     6302                return
    62986303            G2frame.dataWindow.FRMCDataEdit.Enable(G2G.wxID_RUNRMC,True)
    62996304            RMCPdict = data['RMC']['PDFfit']

  • trunk/GSASIIplot.py

    r5112 r5114  
    35793579                        Plot.set_yscale("log",nonpositive='mask') # >=3.3
    35803580                    except:
    3581                         Plot.set_yscale("log",nonposy='mask')
     3581                        Plot.set_yscale("log",nonpositive='mask')
    35823582                    if np.any(W>0.):
    35833583                        lims = [np.min(np.trim_zeros(W))/2.,np.max(Y)*2.]
    … …  
    36203620                            Plot.set_yscale("log",nonpositive='mask') # >=3.3
    36213621                        except:
    3622                             Plot.set_yscale("log",nonposy='mask')
     3622                            Plot.set_yscale("log",nonpositive='mask')
    36233623                        Plot.plot(X,Y,marker=pP,color=colors[0],linewidth=lW,
    36243624                            picker=True,pickradius=3.,clip_on=Clip_on,label=incCptn('obs'))
    … …  
    36323632                            Plot.set_yscale("log",nonpositive='mask')
    36333633                        except:
    3634                             Plot.set_xscale("log",nonposx='mask')
    3635                             Plot.set_yscale("log",nonposy='mask')
     3634                            Plot.set_xscale("log",nonpositive='mask')
     3635                            Plot.set_yscale("log",nonpositive='mask')
    36363636                        if G2frame.ErrorBars:
    36373637                            if Page.plotStyle['sqPlot']:
    … …  
    37213721                        except:
    37223722                            Plot.semilogy(X,Y,color=mcolors.cmap(icolor),
    3723                                 picker=False,nonposy='mask')
     3723                                picker=False,nonpositive='mask')
    37243724                    elif plottype in ['SASD','REFD']:
    37253725                        try:
    … …  
    37283728                        except:
    37293729                            Plot.semilogy(X,Y,color=mcolors.cmap(icolor),
    3730                                 picker=False,nonposy='mask')
     3730                                picker=False,nonpositive='mask')
    37313731                else:
    37323732                    if 'PWDR' in plottype:
    … …  
    37383738                        except:
    37393739                            Plot.loglog(X,Y,mcolors.cmap(icolor),
    3740                                 picker=False,nonposy='mask')
     3740                                picker=False,nonpositive='mask')
    37413741                        Plot.set_ylim(bottom=np.min(np.trim_zeros(Y))/2.,top=np.max(Y)*2.)
    37423742                           
    … …  
    38643864        if G2frame.Contour: # for contour plots expand y-axis to include all histograms
    38653865            G2frame.xylim = (G2frame.xylim[0], (0.,len(PlotList)-1.))
    3866         Plot.set_xlim(G2frame.xylim[0])
    3867         Plot.set_ylim(G2frame.xylim[1])
     3866        if 'PWDR' in plottype:
     3867            Plot.set_xlim(G2frame.xylim[0])
     3868            Plot.set_ylim(G2frame.xylim[1])
    38683869        Page.toolbar.push_current()
    38693870        Page.ToolBarDraw()
    … …  
    61586159    PatternId = G2frame.PatternId
    61596160    data = G2frame.GPXtree.GetItemPyData(G2gd.GetGPXtreeItemId(G2frame,PatternId, 'Models'))
    6160     Bins,Dbins,BinMag = data['Size']['Distribution']
     6161    try:
     6162        Bins,Dbins,BinMag = data['Size']['Distribution']
     6163    except ValueError:  #no data found; skip plot
     6164        return
    61616165    Plot.set_title('Size Distribution')
    61626166    Plot.set_xlabel(r'$D, \AA$',fontsize=14)
    61636167    Plot.set_ylabel(r'$Volume\ distribution,\ f(D)$',fontsize=14)
    61646168    if data['Size']['logBins']:
    6165         Plot.set_xscale("log",nonposx='mask')
     6169        Plot.set_xscale("log",nonpositive='mask')
    61666170        Plot.set_xlim([np.min(2.*Bins)/2.,np.max(2.*Bins)*2.])
    61676171    Plot.bar(2.*Bins-Dbins,BinMag,2.*Dbins,facecolor='white',edgecolor='green')       #plot diameters

  • trunk/GSASIIpwdGUI.py

    r5113 r5114  
    50075007                refs = np.vstack((refList.T[:17+Super],I100,MuStr,CrSize)).T
    50085008            elif 'E' in Inst['Type'][0]:
    5009                 refs = np.vstack((refList.T[:11+Super],I100,MuStr,CrSize)).T
     5009                refs = np.vstack((refList.T[:12+Super],I100,MuStr,CrSize)).T        #last two not shown for now
    50105010        rowLabels = [str(i) for i in range(len(refs))]
    50115011        Types = (4+Super)*[wg.GRID_VALUE_LONG,]+4*[wg.GRID_VALUE_FLOAT+':10,4',]+ \
    … …  
    77017701        if 'C' in inst['Type'][0]:
    77027702            wave = G2mth.getWave(inst)
    7703             keV = 12.397639/wave
     7703            keV = G2mth.wavekE(wave)
    77047704            qLimits = [tth2q(fullLimits[0],wave),tth2q(fullLimits[1],wave)]
    77057705            polariz = inst['Polariz.'][1]

  • trunk/GSASIIsasd.py

    r4919 r5114  
    4343npatand = lambda x: 180.*np.arctan(x)/np.pi
    4444npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
    45 npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
    46 npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
     45# npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
     46# npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
    4747   
    4848###############################################################################

  • trunk/GSASIIstrIO.py

    r5112 r5114  
    33253325            iFin = min(iBeg+9,len(insKeys))
    33263326            for item in insKeys[iBeg:iFin]:
    3327                 if item not in ['Type','Source']:
     3327                if item not in ['Type','Source','Bank']:
    33283328                    ptlbls += '%12s' % (item)
    33293329                    ptstr += '%12.6f' % (Inst[item][1])
    3330                     if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
     3330                    if item in ['Lam1','Lam2','Azimuth','fltPath']:
    33313331                        varstr += 12*' '
    33323332                    else:
    … …  
    33943394            if 'XC' in Type or 'XB' in Type:
    33953395                if pfx+'Lam1' in instDict:
    3396                     controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
     3396                    controlDict[pfx+'keV'] = G2mth.wavekE(instDict[pfx+'Lam1'])
    33973397                else:
    3398                     controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
     3398                    controlDict[pfx+'keV'] = G2mth.wavekE(instDict[pfx+'Lam'])           
    33993399            histDict.update(instDict)
    34003400            histVary += insVary
    … …  
    34313431            if 'X' in Inst['Type'][0]:
    34323432                histDict[pfx+'Lam'] = Inst['Lam'][1]
    3433                 controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
     3433                controlDict[pfx+'keV'] = G2mth.wavekE(histDict[pfx+'Lam'])
    34343434            elif 'NC' in Inst['Type'][0] or 'NB' in Inst['Type'][0]:                   
    34353435                histDict[pfx+'Lam'] = Inst['Lam'][1]
    … …  
    35733573            iFin = min(iBeg+9,len(insKeys))
    35743574            for name in insKeys[iBeg:iFin]:
    3575                 if name not in  ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
     3575                if name not in  ['Type','Lam1','Lam2','Azimuth','Source','fltPath','Bank']:
    35763576                    ptlbls += '%12s' % (name)
    35773577                    ptstr += '%12.6f' % (Inst[name][1])
    … …  
    36493649                pFile.write(' Instrument type: %s\n'%Sample['Type'])
    36503650                if calcControls != None:    #skipped in seqRefine
    3651                     if 'X' in Inst['Type'][0]:
     3651                    if 'X' in Inst['Type'][0] and 'E' not in Inst['Type'][0]:
    36523652                        PrintFprime(calcControls['FFtables'],pfx,pFile)
    36533653                    elif 'NC' in Inst['Type'][0]:

  • trunk/GSASIIstrMath.py

    r5112 r5114  
    26932693            Mgam = 0.018*refl[4+im]**2*tand(refl[5+im]/2.)*np.sqrt(Sum)/np.pi
    26942694    elif 'E' in calcControls[hfx+'histType']:
    2695         return 0.,0.
     2695        return 0.,0.        #wrong - fix later if wanted
    26962696    elif 'T' in calcControls[hfx+'histType']:       #All checked & OK
    26972697        #crystallite size
    … …  
    28052805        sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
    28062806    elif 'E' in calcControls[hfx+'histType']:
    2807         return {},{}
     2807        return {},{}            #wrong - fix later if wanted
    28082808    else:   #'T'OF - All checked & OK
    28092809        if calcControls[phfx+'SizeType'] == 'isotropic':    #OK
    … …  
    28942894        else:               #Debye-Scherrer - simple but maybe not right
    28952895            pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+(parmDict[hfx+'DisplaceY']+parmDict[phfx+'LayerDisp'])*sind(pos))
     2896    elif 'E' in calcControls[hfx+'histType']:
     2897        pos = G2mth.wavekE(2.0*d*sind(parmDict[hfx+'2-theta']/2.0))+parmDict[hfx+'ZE']+parmDict[hfx+'YE']*d+parmDict[hfx+'XE']*d**2
    28962898    elif 'T' in calcControls[hfx+'histType']:
    28972899        pos = parmDict[hfx+'difC']*d+parmDict[hfx+'difA']*d**2+parmDict[hfx+'difB']/d+parmDict[hfx+'Zero']
    … …  
    29302932            dpdYd = -shft*sind(pos)
    29312933            return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd,dpdV
     2934    elif 'E' in calcControls[hfx+'histType']:
     2935        tth = parmDict[hfx+'2-theta']/2.0
     2936        dpdZE = 1.0
     2937        dpdYE = dsp
     2938        dpdXE = dsp**2
     2939        dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])*refl[5+im]*dsp**2/2.
     2940        dpdTTh = -12.397639/(2.*tand(tth)*dsp*sind(tth))
     2941        return dpdA,dpdTTh,dpdXE,dpdYE,dpdZE
    29322942    elif 'T' in calcControls[hfx+'histType']:
    29332943        dpdA = -np.array([h**2,k**2,l**2,h*k,h*l,k*l])*parmDict[hfx+'difC']*dsp**3/2.
    … …  
    29462956    h,k,l = refl[:3]
    29472957    if laue in ['m3','m3m']:
    2948         Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
    2949             refl[4+im]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
     2958        hkl1 = (h**2+k**2+l**2)
     2959        hkl2 = ((h*k)**2+(h*l)**2+(k*l)**2)/hkl1**2
     2960        Dij = parmDict[phfx+'D11']*hkl1+parmDict[phfx+'eA']*hkl2
    29502961    elif laue in ['6/m','6/mmm','3m1','31m','3']:
    29512962        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
    … …  
    29692980    if 'C' in calcControls[hfx+'histType'] or 'B' in calcControls[hfx+'histType']:
    29702981        return -180.*Dij*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
     2982    elif 'E' in calcControls[hfx+'histType']:
     2983        return Dij*(12.397639/(2.*dsp*sind(tth)))
    29712984    else:
    2972         return -Dij*parmDict[hfx+'difC']*refl[4+im]**2/2.
     2985        return -Dij*parmDict[hfx+'difC']*0.5*refl[4+im]**2
    29732986           
    29742987def GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict):
    … …  
    29792992    if laue in ['m3','m3m']:
    29802993        dDijDict = {phfx+'D11':h**2+k**2+l**2,
    2981             phfx+'eA':refl[4+im]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
     2994            phfx+'eA':((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
    29822995    elif laue in ['6/m','6/mmm','3m1','31m','3']:
    29832996        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
    … …  
    30023015        for item in dDijDict:
    30033016            dDijDict[item] *= 180.0*refl[4+im]**2*tand(refl[5+im]/2.0)/np.pi
     3017    elif 'E' in calcControls[hfx+'histType']:
     3018        for item in dDijDict:
     3019            dDijDict[item] *= refl[5+im]*refl[4+im]**2/2.
    30043020    else:
    30053021        for item in dDijDict:
    … …  
    31253141                                    if parmDict[phfx+'LeBail']:
    31263142                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
     3143                elif 'E' in calcControls[hfx+'histType']:
     3144                    for iref,refl in enumerate(refDict['RefList']):
     3145                        if useMP:
     3146                            profArgs[iref%G2mp.ncores].append((refl,iref))
     3147                        else:
     3148                            icod = G2mp.ComputeFobsSqED(refl,iref)
     3149                            if type(icod) is tuple:
     3150                                refl[8+im] = icod[0]
     3151                                sumInt += icod[1]
     3152                                if parmDict[phfx+'LeBail']:
     3153                                    refl[9+im] = refl[8+im]
     3154                            elif icod == -1:
     3155                                refl[3+im] *= -1
     3156                                nExcl += 1
     3157                            elif icod == -2:
     3158                                break
     3159                    if useMP:
     3160                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqEDbatch,profArgs):
     3161                            sumInt += sInt
     3162                            for refl8im,irefl in resList:
     3163                                if refl8im is None:
     3164                                    refDict['RefList'][irefl][3+im] *= -1
     3165                                    nExcl += 1
     3166                                else:
     3167                                    refDict['RefList'][irefl][8+im] = refl8im
     3168                                    if parmDict[phfx+'LeBail']:
     3169                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
    31273170                elif 'B' in calcControls[hfx+'histType']:
    31283171                    for iref,refl in enumerate(refDict['RefList']):
    … …  
    32083251   
    32093252    def GetReflSigGamED(refl,im,G,GB,phfx,calcControls,parmDict):
    3210         sig = parmDict[hfx+'A']*refl[4+im]**2+parmDict[hfx+'B']*refl[4+im]+parmDict[hfx+'C']
     3253        sig = parmDict[hfx+'A']*refl[5+im]**2+parmDict[hfx+'B']*refl[5+im]+parmDict[hfx+'C']
    32113254        gam = 0.001
    3212         Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
    3213         sig += Ssig
    3214         gam += Sgam
     3255        # Ssig,Sgam = GetSampleSigGam(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
     3256        # sig += Ssig
     3257        # gam += Sgam
    32153258        return sig,gam
    32163259       
    … …  
    33083351                Uniq = np.inner(refl[:3],SGMT)
    33093352                refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,phfx,calcControls,parmDict)         #corrected reflection position
     3353                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
    33103354                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
    33113355                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
    … …  
    33543398                        fp2 = G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))[0]
    33553399                        yc[iBeg:iFin] += refl[11+im]*refl[9+im]*kRatio*fp2       #and here
     3400        elif 'E' in calcControls[hfx+'histType']:
     3401           
     3402            for iref,refl in enumerate(refDict['RefList']):
     3403                if im:
     3404                    h,k,l,m = refl[:4]
     3405                else:
     3406                    h,k,l = refl[:3]
     3407                Uniq = np.inner(refl[:3],SGMT)
     3408                refl[5+im] = GetReflPos(refl,im,0.0,A,pfx,hfx,phfx,calcControls,parmDict)         #corrected reflection position
     3409                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
     3410                refl[6+im:8+im] = GetReflSigGamED(refl,im,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
     3411                refl[11+im] = 1.0       #no intensity corrections; default = 1.0
     3412                 
     3413                if Phase['General'].get('doPawley'):
     3414                    try:
     3415                        if im:
     3416                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
     3417                        else:
     3418                            pInd = pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
     3419                        refl[9+im] = parmDict[pInd]
     3420                    except KeyError:
     3421                        continue
     3422                Wd,fmin,fmax = G2pwd.getWidthsED(refl[5+im],refl[6+im])
     3423                iBeg = np.searchsorted(x,refl[5+im]-fmin)
     3424                iFin = np.searchsorted(x,refl[5+im]+fmax)
     3425                if not iBeg+iFin:       #peak below low limit - skip peak
     3426                    continue
     3427                elif not iBeg-iFin:     #peak above high limit - done
     3428                    break
     3429                elif iBeg > iFin:   #bad peak coeff - skip
     3430                    badPeak = True
     3431                    continue
     3432                if useMP:
     3433                    profArgs[iref%ncores].append((refl[5+im],refl,iBeg,iFin,1.))
     3434                else:
     3435                    fp = G2pwd.getPsVoigt(refl[5+im],refl[6+im]*1.e4,.001,ma.getdata(x[iBeg:iFin]))[0]
     3436                    yc[iBeg:iFin] += refl[9+im]*fp
     3437           
    33563438        elif 'B' in calcControls[hfx+'histType']:
    33573439            for iref,refl in enumerate(refDict['RefList']):
    … …  
    34063488                refl[5+im] = GetReflPos(refl,im,0.0,A,pfx,hfx,phfx,calcControls,parmDict)         #corrected reflection position - #TODO - what about tabluated offset?
    34073489                Lorenz = sind(abs(parmDict[hfx+'2-theta'])/2)*refl[4+im]**4                                                #TOF Lorentz correction
    3408 #                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
     3490                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
    34093491                refl[6+im:8+im] = GetReflSigGamTOF(refl,im,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
    34103492                refl[12+im:14+im] = GetReflAlpBet(refl,im,hfx,parmDict)             #TODO - skip if alp, bet tabulated?
    … …  
    34473529        elif useMP and 'B' in calcControls[hfx+'histType']:
    34483530            for y in MPpool.imap_unordered(G2mp.ComputePwdrProfPink,profArgs):
     3531                yc += y
     3532            MPpool.terminate()
     3533        elif useMP and 'E' in calcControls[hfx+'histType']:
     3534            for y in MPpool.imap_unordered(G2mp.ComputePwdrProfED,profArgs):
    34493535                yc += y
    34503536            MPpool.terminate()
    … …  
    36013687                wave = refl[14+im]
    36023688            if 'E' in calcControls[hfx+'histType']:
    3603                 dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = [0.,0.,0.,0.,0.,0.,0.,0.]
     3689                dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = [[],[],[],[],[],[],[],[]]
    36043690            else:
    36053691                dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
    … …  
    36113697                Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im]/1.e4,refl[7+im]/100.)
    36123698            elif 'E' in calcControls[hfx+'histType']:
    3613                 Wd,fmin,fmax = G2pwd.getWidthsED(refl[5+im],refl[10+im])
     3699                Wd,fmin,fmax = G2pwd.getWidthsED(refl[5+im],refl[6+im]*1.e4)
    36143700            iBeg = np.searchsorted(x,refl[5+im]-fmin)
    36153701            iFin = np.searchsorted(x,refl[5+im]+fmax)
    … …  
    36653751                    break
    36663752                dMdpk = np.zeros(shape=(4,lenBF))
    3667                 dMdipk = G2pwd.getdPsVoigt(refl[5+im],refl[12+im]*10**4,0.001,ma.getdata(x[iBeg:iFin]))
     3753                dMdipk = G2pwd.getdPsVoigt(refl[5+im],refl[6+im]*10**4,0.001,ma.getdata(x[iBeg:iFin]))
    36683754                for i in range(4):
    3669                     dMdpk[i] = refl[11+im]*refl[9+im]*dMdipk[i]
    3670                 dervDict = {'int':dMdpk[0],'pos':-dMdpk[1],'sig':dMdpk[4]*1.e4}   
     3755                    dMdpk[i] = refl[9+im]*dMdipk[i]
     3756                dervDict = {'int':dMdpk[0],'pos':-dMdpk[1],'sig':dMdpk[2]*1.e4}   
    36713757            if Phase['General'].get('doPawley'):
    36723758                dMdpw = np.zeros(len(x))
    … …  
    37133799                    hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],hfx+'DisplaceX':[dpdX,'pos'],hfx+'DisplaceY':[dpdY,'pos']}
    37143800            elif 'E' in calcControls[hfx+'histType']:
    3715                 dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY,dpdV = GetReflPosDerv(refl,im,wave,A,pfx,hfx,phfx,calcControls,parmDict)
    3716                 names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],hfx+'Lam':[dpdw,'pos'],
    3717                     hfx+'Zero':[dpdZ,'pos'],hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'Z':[1.0,'gam'],
    3718                     hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],hfx+'Polariz.':[dIdpola,'int'],
    3719                     hfx+'alpha-0':[1.0,'alp'],hfx+'alpha-1':[tanth,'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[tanth,'bet'],
    3720                     hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],hfx+'DisplaceX':[dpdX,'pos'],hfx+'DisplaceY':[dpdY,'pos']}
    3721                
     3801                dpdA,dpdTTh,dpdXE,dpdYE,dpdZE = GetReflPosDerv(refl,im,0.0,A,pfx,hfx,phfx,calcControls,parmDict)
     3802                names = {hfx+'Scale':[dIdsh,'int'],
     3803                    hfx+'A':[refl[5+im]**2,'sig'],hfx+'B':[refl[5+im],'sig'],hfx+'C':[1.0,'sig'],
     3804                    hfx+'XE':[dpdXE,'pos'],hfx+'YE':[dpdYE,'pos'],hfx+'ZE':[dpdZE,'pos']   }
    37223805            for name in names:
    37233806                item = names[name]

  • trunk/imports/G2sad_xye.py

    r3136 r5114  
    7676                try:
    7777                    data = [float(val) for val in vals]
    78                     x.append(float(data[0]))
    7978                    f = float(data[1])
    8079                    if f <= 0.0:
    81                         del x[-1]
    82                         continue
     80                        break
     81                        # del x[-1]
     82                        # continue
    8383                    elif len(vals) > 2:
    8484                        y.append(float(data[1]))
    … …  
    8787                        y.append(float(data[1]))
    8888                        w.append(1.0/float(data[1]))
    89                 except ValueError:
    90                     msg = 'Error in line '+str(i+1)
     89                    x.append(float(data[0]))
     90                except ValueError:
     91                    msg = 'Error in line :%s'%S
    9192                    print (msg)
    9293                    continue
    … …  
    170171                try:
    171172                    data = [float(val) for val in vals]
    172                     x.append(float(data[0])/10.)        #convert nm-1 to A-1
    173173                    f = float(data[1])
    174174                    if f <= 0.0:
    175                         x.pop()
    176                         continue
     175                        break
     176                        # x.pop()
     177                        # continue
    177178                    elif len(vals) > 2:
    178179                        y.append(float(data[1]))
    … …  
    181182                        y.append(float(data[1]))
    182183                        w.append(1.0/float(data[1]))
    183                 except ValueError:
    184                     msg = 'Error in line '+str(i+1)
     184                    x.append(float(data[0])/10.)        #convert nm-1 to A-1
     185                except ValueError:
     186                    msg = 'Error in line :%s'%S
    185187                    print (msg)
    186188                    continue
    … …  
    212214        return True
    213215
    214 class txt_CWNeutronReaderClass(G2obj.ImportSmallAngleData):
    215     'Routines to import neutron CW q SAXD data from a .nsad or .ndat file'
     216class txt_NeutronReaderClass(G2obj.ImportSmallAngleData):
     217    'Routines to import neutron q SAXD data from a .nsad or .ndat file'
    216218    def __init__(self):
    217219        super(self.__class__,self).__init__( # fancy way to self-reference
    218220            extensionlist=('.nsad','.ndat'),
    219221            strictExtension=False,
    220             formatName = 'q (A-1) step neutron CW QIE data',
     222            formatName = 'q (A-1) step neutron QIE data',
    221223            longFormatName = 'q (A-1) stepped neutron CW text data file in Q,I,E order; E optional'
    222224            )
    … …  
    252254            if len(S) == 1:     #skip blank line
    253255                continue
    254             if '=' in S:
     256            if '=' in S or '#' in S:
    255257                self.comments.append(S[:-1])
    256258                if 'wave' in S.split('=')[0].lower():
    … …  
    264266                try:
    265267                    data = [float(val) for val in vals]
    266                     x.append(float(data[0]))
    267268                    f = float(data[1])
    268269                    if f <= 0.0:
    269                         y.append(0.0)
    270                         w.append(1.0)
     270                        break
     271                        # y.append(0.0)
     272                        # w.append(1.0)
    271273                    elif len(vals) > 2:
    272274                        y.append(float(data[1]))
    … …  
    275277                        y.append(float(data[1]))
    276278                        w.append(1.0/float(data[1]))
    277                 except ValueError:
    278                     msg = 'Error in line '+str(i+1)
     279                    x.append(float(data[0]))
     280                except ValueError:
     281                    msg = 'Error in line :%s'%S
    279282                    print (msg)
    280283                    continue
    … …  
    308311        return True
    309312
    310 class txt_nmCWNeutronReaderClass(G2obj.ImportSmallAngleData):
    311     'Routines to import neutron CW q in nm-1 SAXD data from a .nsad or .ndat file'
     313class txt_nmNeutronReaderClass(G2obj.ImportSmallAngleData):
     314    'Routines to import neutron q in nm-1 SAXD data from a .nsad or .ndat file'
    312315    def __init__(self):
    313316        super(self.__class__,self).__init__( # fancy way to self-reference
    314317            extensionlist=('.nsad','.ndat'),
    315318            strictExtension=False,
    316             formatName = 'q (nm-1) step neutron CW QIE data',
    317             longFormatName = 'q (nm-1) stepped neutron CW text data file in Q,I,E order; E optional'
     319            formatName = 'q (nm-1) step neutron QIE data',
     320            longFormatName = 'q (nm-1) stepped neutron text data file in Q,I,E order; E optional'
    318321            )
    319322
    … …  
    348351            if len(S) == 1:     #skip blank line
    349352                continue
    350             if '=' in S:
     353            if '=' in S or '#' in S:
    351354                self.comments.append(S[:-1])
    352355                if 'wave' in S.split('=')[0].lower():
    … …  
    360363                try:
    361364                    data = [float(val) for val in vals]
    362                     x.append(float(data[0])/10.)    #convert to A-1
    363365                    f = float(data[1])
    364366                    if f <= 0.0:
    365                         y.append(0.0)
    366                         w.append(1.0)
     367                        break
     368                        # y.append(0.0)
     369                        # w.append(1.0)
    367370                    elif len(vals) > 2:
    368371                        y.append(float(data[1]))
    … …  
    371374                        y.append(float(data[1]))
    372375                        w.append(1.0/float(data[1]))
    373                 except ValueError:
    374                     msg = 'Error in line '+str(i+1)
     376                    x.append(float(data[0])/10.)    #convert to A-1
     377                except ValueError:
     378                    msg = 'Error in line :%s'%S
    375379                    print (msg)
    376380                    continue

  • trunk/testDeriv.py

    r4521 r5114  
    242242            mmin = np.min(dMdV[names.index(name)])
    243243            mmax = np.max(dMdV[names.index(name)])
    244             print('parameter:',name,self.parmDict[name],delt,mmin,mmax)
    245244            if name in self.varylist:
    246245                self.values[self.varylist.index(name)] -= delt
    … …  
    263262                self.parmDict[name] -= delt   
    264263            Mn = (M1-M0)/(2.*abs(delt))
     264            print('parameter:',name,self.parmDict[name],delt,mmin,mmax,np.sum(M0),np.sum(M1),np.sum(Mn))
    265265            hplot.plot(Mn,'r',label='numeric deriv')
    266 #            hplot.plot(M2-Mn,'g',label='diff')
    267 #            GSASIIpath.IPyBreak()
    268266            hplot.legend(loc='best')           
    269267           
Note: See TracChangeset for help on using the changeset viewer.