Changeset 409


Ignore:
Timestamp:
Nov 8, 2011 4:27:01 PM (10 years ago)
Author:
vondreele
Message:

make it do CW neutrons - 1st pass at it

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIElem.py

    r380 r409  
    1919
    2020def GetFormFactorCoeff(El):
    21     """Read form factor coefficients from `atomdata.asc` file
     21    """Read X-ray form factor coefficients from `atomdata.asc` file
    2222
    2323    :param El: element 1-2 character symbol case irrevelant
    2424    :return: `FormFactors`: list of form factor dictionaries
    2525   
    26     Each form factor dictionary is:
     26    Each X-ray form factor dictionary is:
    2727   
    2828    * `Symbol`: 4 character element symbol with valence (e.g. 'NI+2')
     
    3030    * `fa`: 4 A coefficients
    3131    * `fb`: 4 B coefficients
    32     * `fc`: C coefficient
     32    * `fc`: C coefficient
     33   
    3334    """
    3435    ElS = El.upper()
     
    9394    S = '1'
    9495    AtomInfo = {}
     96    Isotopes = {}
    9597    Mass = []
    9698    while S:
     
    100102                if not Mass:                                 #picks 1st one; natural abundance or 1st isotope
    101103                    Mass = float(S[10:19])
    102                 if S[5:9] == '_SIZ':
     104                if S[6] in [' ','1','2','3','4','5','6','7','8','9']:                       
     105                    isoName = S[6:9]
     106                    if isoName == '   ':
     107                        isoName = 'Nat. Abund.'              #natural abundance
     108                    if S[76:78] in ['LS','BW']:     #special anomalous scattering length info
     109                        St = [S[10:19],S[19:25],S[25:31],S[31:38],S[38:44],S[44:50],
     110                            S[50:56],S[56:62],S[62:68],S[68:74],]
     111                        Vals = []
     112                        for item in St:
     113                            if item.strip():
     114                                Vals.append(float(item.strip()))
     115                        Isotopes[isoName.rstrip()] = Vals                       
     116                    else:
     117                        Isotopes[isoName.rstrip()] = [float(S[10:19]),float(S[19:25])]
     118                elif S[5:9] == '_SIZ':
    103119                    Z=int(S[:2])
    104120                    Symbol = S[3:5].strip().lower().capitalize()
     
    108124                    Color = ET.ElTable[Elements.index(Symbol)][6]
    109125    FFdata.close()
    110     AtomInfo={'Symbol':Symbol,'Mass':Mass,'Z':Z,'Drad':Drad,'Arad':Arad,'Vdrad':Vdrad,'Color':Color}   
     126    AtomInfo={'Symbol':Symbol,'Isotopes':Isotopes,'Mass':Mass,'Z':Z,'Drad':Drad,'Arad':Arad,'Vdrad':Vdrad,'Color':Color}   
    111127    return AtomInfo
    112128     
     
    244260    t = -fb[:,np.newaxis]*SQ
    245261    return np.sum(fa[:,np.newaxis]*np.exp(t)[:],axis=0)+El['fc']
     262       
     263def BlenFac(El,wave):
     264    pass
     265   
     266#        F(I) = BLEN(I)
     267#        IF ( BFAN(1,I).NE.0.0 ) THEN
     268#          EMEV = 81.80703/XRAY**2
     269#          GAM2 = BFAN(4,I)**2
     270#          T1 = EMEV-BFAN(3,I)
     271#          D1 = T1**2+GAM2
     272#          T2 = EMEV-BFAN(6,I)
     273#          D2 = T2**2+GAM2
     274#          T3 = EMEV-BFAN(8,I)
     275#          D3 = T3**2+GAM2
     276#          FP(I) = BFAN(1,I)*(T1/D1+BFAN(5,I)*T2/D2+BFAN(7,I)*T3/D3)
     277#          FPP(I) = -BFAN(2,I)*(1.0/D1+BFAN(5,I)/D2+BFAN(7,I)/D3)
     278#        ELSE
     279#          FP(I) = 0.0
     280#          FPP(I) = 0.0
     281#        END IF
     282   
    246283   
    247284def ComptonFac(El,SQ):
  • trunk/GSASIIphsGUI.py

    r407 r409  
    147147        atomData = data['Atoms']
    148148        generalData['AtomTypes'] = []
     149        generalData['Isotopes'] = {}
     150        if 'Isotope' not in generalData:
     151            generalData['Isotope'] = {}
    149152        generalData['NoAtoms'] = {}
    150153        generalData['BondRadii'] = []
     
    166169                Info = G2elem.GetAtomInfo(atom[ct])
    167170                generalData['AtomTypes'].append(atom[ct])
     171                generalData['Isotopes'][atom[ct]] = Info['Isotopes']
    168172                generalData['BondRadii'].append(Info['Drad'])
    169173                generalData['AngleRadii'].append(Info['Arad'])
     
    320324            except ValueError:
    321325                pass
    322             pawlVal.SetValue("%.3f"%(generalData['Pawley dmin']))          #reset in case of error           
    323                                    
     326            pawlVal.SetValue("%.3f"%(generalData['Pawley dmin']))          #reset in case of error       
     327           
     328        def OnIsotope(event):
     329            Obj = event.GetEventObject()
     330            generalData['Isotope'][Indx[Obj.GetId()]] = Obj.GetValue() #mass too
     331                                               
    324332        cellGUIlist = [[['m3','m3m'],4,zip([" Unit cell: a = "," Vol = "],["%.5f","%.3f"],[True,False],[0,0])],
    325333        [['3R','3mR'],6,zip([" a = "," alpha = "," Vol = "],["%.5f","%.3f","%.3f"],[True,True,False],[0,2,0])],
     
    393401        mainSizer.Add((5,5),0)
    394402       
     403        Indx = {}
    395404        if len(generalData['AtomTypes']):
    396405            mass = 0.
     
    413422            mainSizer.Add((5,5),0)
    414423           
    415             elemSizer = wx.FlexGridSizer(7,len(generalData['AtomTypes'])+1,1,1)
     424            elemSizer = wx.FlexGridSizer(8,len(generalData['AtomTypes'])+1,1,1)
    416425            elemSizer.Add(wx.StaticText(dataDisplay,label='Elements'),0,wx.ALIGN_CENTER_VERTICAL)
    417426            for elem in generalData['AtomTypes']:
     
    419428                typTxt.SetBackgroundColour(VERY_LIGHT_GREY)
    420429                elemSizer.Add(typTxt,0,wx.ALIGN_CENTER_VERTICAL)
     430            elemSizer.Add(wx.StaticText(dataDisplay,label='Isotope'),0,wx.ALIGN_CENTER_VERTICAL)
     431            for elem in generalData['AtomTypes']:
     432                choices = generalData['Isotopes'][elem].keys()
     433                if elem not in generalData['Isotope']:
     434                    generalData['Isotope'][elem] = 'Nat. Abund.'
     435                isoSel = wx.ComboBox(dataDisplay,-1,value=generalData['Isotope'][elem],choices=choices,
     436                    style=wx.CB_READONLY|wx.CB_DROPDOWN)
     437                isoSel.Bind(wx.EVT_COMBOBOX,OnIsotope)
     438                Indx[isoSel.GetId()] = elem
     439                elemSizer.Add(isoSel,1,wx.ALIGN_CENTER_VERTICAL)
    421440            elemSizer.Add(wx.StaticText(dataDisplay,label='No. per cell'),0,wx.ALIGN_CENTER_VERTICAL)
    422441            for elem in generalData['AtomTypes']:
  • trunk/GSASIIstruct.py

    r408 r409  
    380380    return FFtable
    381381   
     382def GetBLtable(General):
     383    ''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General
     384    input:
     385        General = dictionary of phase info.; includes AtomTypes & Isotopes
     386    return:
     387        BLtable = dictionary of scattering length data; key is atom type
     388    '''
     389    atomTypes = General['AtomTypes']
     390    BLtable = {}
     391    isotopes = General['Isotopes']
     392    isotope = General['Isotope']
     393    for El in atomTypes:
     394        BLtable[El] = isotopes[El][isotope[El]]
     395    return BLtable
     396       
    382397def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
    383398    if SGLaue in ['-1','2/m','mmm']:
     
    492507    phaseConstr = {}
    493508    pawleyLookup = {}
    494     FFtables = {}
     509    FFtables = {}                   #scattering factors - xrays
     510    BLtables = {}                   # neutrons
    495511    Natoms = {}
    496512    AtMults = {}
     
    503519        pfx = str(pId)+'::'
    504520        FFtable = GetFFtable(General)
     521        BLtable = GetBLtable(General)
    505522        FFtables.update(FFtable)
     523        BLtables.update(BLtable)
    506524        Atoms = PhaseData[name]['Atoms']
    507525        try:
     
    593611            phaseVary += pawleyVary
    594612               
    595     return Natoms,phaseVary,phaseDict,pawleyLookup,FFtables
     613    return Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
    596614   
    597615def getVCov(varyNames,varyList,covMatrix):
     
    14811499                PrintBackgroundSig(Background,backSig)
    14821500
    1483 def GetAtomFXU(pfx,FFtables,calcControls,parmDict):
     1501def GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict):
    14841502    Natoms = calcControls['Natoms'][pfx]
    14851503    Tdata = Natoms*[' ',]
     
    14881506    Fdata = np.zeros(Natoms)
    14891507    FFdata = []
     1508    BLdata = []
    14901509    Xdata = np.zeros((3,Natoms))
    14911510    dXdata = np.zeros((3,Natoms))
     
    15031522                keys[key][iatm] = parmDict[parm]
    15041523        FFdata.append(FFtables[Tdata[iatm]])
    1505     return FFdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
     1524        BLdata.append(BLtables[Tdata[iatm]])
     1525    return FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
    15061526   
    15071527def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
     
    15211541    Mast = twopisq*np.multiply.outer(ast,ast)
    15221542    FFtables = calcControls['FFtables']
    1523     FFdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,calcControls,parmDict)
    1524     FP = np.array([El[hfx+'FP'] for El in FFdata])
    1525     FPP = np.array([El[hfx+'FPP'] for El in FFdata])
     1543    BLtables = calcControls['BLtables']
     1544    FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict)
     1545    if 'N' in parmDict[hfx+'Type']:
     1546        FP = 0.
     1547        FPP = 0.
     1548    else:
     1549        FP = np.array([El[hfx+'FP'] for El in FFdata])
     1550        FPP = np.array([El[hfx+'FPP'] for El in FFdata])
    15261551    maxPos = len(SGData['SGOps'])
    15271552    Uij = np.array(G2lat.U6toUij(Uijdata))
     
    15311556        H = refl[:3]
    15321557        SQ = 1./(2.*refl[4])**2
    1533         FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata])
     1558        if 'N' in parmDict[hfx+'Type']:
     1559            FF = np.array([El[1] for El in BLdata])
     1560        else:       #'X'
     1561            FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata])
    15341562        SQfactor = 4.0*SQ*twopisq
    15351563        Uniq = refl[11]
     
    15611589    Mast = twopisq*np.multiply.outer(ast,ast)
    15621590    FFtables = calcControls['FFtables']
    1563     FFdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,calcControls,parmDict)
    1564     FP = np.array([El[hfx+'FP'] for El in FFdata])
    1565     FPP = np.array([El[hfx+'FPP'] for El in FFdata])
     1591    BLtables = calcControls['BLtables']
     1592    FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict)
     1593    if 'N' in parmDict[hfx+'Type']:
     1594        FP = 0.
     1595        FPP = 0.
     1596    else:
     1597        FP = np.array([El[hfx+'FP'] for El in FFdata])
     1598        FPP = np.array([El[hfx+'FPP'] for El in FFdata])
    15661599    maxPos = len(SGData['SGOps'])       
    15671600    Uij = np.array(G2lat.U6toUij(Uijdata))
     
    15751608        H = np.array(refl[:3])
    15761609        SQ = 1./(2.*refl[4])**2          # or (sin(theta)/lambda)**2
    1577         FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata])
     1610        if 'N' in parmDict[hfx+'Type']:
     1611            FF = np.array([El[1] for El in BLdata])
     1612        else:       #'X'
     1613            FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata])
    15781614        SQfactor = 8.0*SQ*np.pi**2
    15791615        Uniq = refl[11]
     
    17801816def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
    17811817    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
    1782     Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
     1818    if 'X' in parmDict[hfx+'Type']:
     1819        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
    17831820    Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
    17841821    if pfx+'SHorder' in parmDict:
     
    17891826    dIdsh = 1./parmDict[hfx+'Scale']
    17901827    dIdsp = 1./parmDict[phfx+'Scale']
    1791     pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
    1792     dIdPola /= pola
     1828    if 'X' in parmDict[hfx+'Type']:
     1829        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
     1830        dIdPola /= pola
     1831    else:       #'N'
     1832        dIdPola = 0.0
    17931833    POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
    17941834    for iPO in dIdPO:
     
    23752415        print ' *** Refine aborted ***'
    23762416        raise Exception
    2377     Natoms,phaseVary,phaseDict,pawleyLookup,FFtables = GetPhaseData(Phases)
     2417    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases)
    23782418    calcControls['Natoms'] = Natoms
    23792419    calcControls['FFtables'] = FFtables
     2420    calcControls['BLtables'] = BLtables
    23802421    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
    23812422    calcControls.update(controlDict)
     
    24852526        print ' *** Refine aborted ***'
    24862527        raise Exception
    2487     Natoms,phaseVary,phaseDict,pawleyLookup,FFtables = GetPhaseData(Phases,False)
     2528    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,False)
    24882529    if 'Seq Data' in Controls:
    24892530        histNames = Controls['Seq Data']
     
    25022543        calcControls['Natoms'] = Natoms
    25032544        calcControls['FFtables'] = FFtables
     2545        calcControls['BLtables'] = BLtables
    25042546        varyList = []
    25052547        parmDict = {}
Note: See TracChangeset for help on using the changeset viewer.