Changeset 411


Ignore:
Timestamp:
Nov 11, 2011 2:14:35 PM (11 years ago)
Author:
vondreele
Message:

fix very old (!) bug in psvfcj.for
implement neutron scattering lengths in GSASII including the dozen anomalous scattering isotopes
isotope choice is in General
so GSASII will now do neutron CW Rietveld refinements
some cleanup of the constraints GUI
remove varyList from GSASIImapvars.py globals

Location:
trunk
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASII.py

    r406 r411  
    639639                name = wx.TextCtrl(panel,-1,item[1],size=wx.Size(200,20))
    640640                name.SetEditable(False)
    641                 scale = wx.TextCtrl(panel,id,str(item[0]),style=wx.TE_PROCESS_ENTER)
     641                scale = wx.TextCtrl(panel,id,'%.3f'%(item[0]),style=wx.TE_PROCESS_ENTER)
    642642                scale.Bind(wx.EVT_TEXT_ENTER,self.OnScaleChange)
    643643                scale.Bind(wx.EVT_KILL_FOCUS,self.OnScaleChange)
    … …  
    680680                if value and '-' not in value[0]:
    681681                    print 'bad input - numbers only'
    682                     self.FindWindowById(id).SetValue('0.0')
     682                    self.FindWindowById(id).SetValue('0.000')
    683683           
    684684        def OnOk(self,event):
    … …  
    14791479        parmDict = {}
    14801480        Histograms,Phases = self.GetUsedHistogramsAndPhasesfromTree()
    1481         Natoms,phaseVary,phaseDict,pawleyLookup,FFtable = G2str.GetPhaseData(Phases,Print=False)       
     1481        Natoms,phaseVary,phaseDict,pawleyLookup,FFtable,BLtable = G2str.GetPhaseData(Phases,Print=False)       
    14821482        hapVary,hapDict,controlDict = G2str.GetHistogramPhaseData(Phases,Histograms,Print=False)
    14831483        histVary,histDict,controlDict = G2str.GetHistogramData(Histograms,Print=False)

  • trunk/GSASIIElem.py

    r409 r411  
    252252def ScatFac(El, SQ):
    253253    """compute value of form factor
    254     @param El: element dictionary  defined in GetFormFactorCoeff
     254    @param El: element dictionary defined in GetFormFactorCoeff
    255255    @param SQ: (sin-theta/lambda)**2
    256256    @return: real part of form factor
    … …  
    261261    return np.sum(fa[:,np.newaxis]*np.exp(t)[:],axis=0)+El['fc']
    262262       
    263 def 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    
     263def BlenRes(BLdata,wave):
     264    FP = []
     265    FPP = []
     266    Emev = 81.80703/wave**2
     267    for BL in BLdata:
     268        if len(BL) >= 6:
     269            Emev = 81.80703/wave**2
     270            G2 = BL[5]**2
     271            T = [Emev-BL[4],0,0]
     272            D = [T**2+G2,0,0]
     273            fp = T/D
     274            fpp = 1.0/D
     275            if len(BL) == 8:
     276                T = Emev-BL[7]
     277                D = T**2+G2
     278                fp += BL[6]*T/D
     279                fpp += BL[6]/D
     280            if len(BL) == 10:
     281                T = Emev-BL[9]
     282                D = T**2+G2
     283                fp += BL[8]*T/D
     284                fpp += BL[8]/D
     285            FP.append(BL[2]*fp)
     286            FPP.append(-BL[3]*fpp)
     287        else:
     288            FP.append(0.0)
     289            FPP.append(0.0)
     290    return np.array(FP),np.array(FPP)
    283291   
    284292def ComptonFac(El,SQ):

  • trunk/GSASIIgrid.py

    r408 r411  
    783783    Histograms,Phases = self.GetUsedHistogramsAndPhasesfromTree()
    784784    AtomDict = dict([Phases[phase]['pId'],Phases[phase]['Atoms']] for phase in Phases)
    785     Natoms,phaseVary,phaseDict,pawleyLookup,FFtable = G2str.GetPhaseData(Phases,Print=False)
     785    Natoms,phaseVary,phaseDict,pawleyLookup,FFtable,BLtable = G2str.GetPhaseData(Phases,Print=False)
    786786    phaseList = []
    787787    for item in phaseDict:
    … …  
    852852                constr = []
    853853                for item in varbs[1:]:
    854 #                    constr += [[{FrstVarb:1.0},{item:-1.0},None,None],]
    855854                    constr += [[[1.0,FrstVarb],[-1.0,item],None,None],]
    856855                return constr           #multiple constraints
    857856            elif 'function' in constType:
    858 #                constr = map(dict,zip(varbs,[1.0 for i in range(len(varbs))]))
    859857                constr = map(list,zip([1.0 for i in range(len(varbs))],varbs))
    860858                return [constr+[0.0,False]]         #just one constraint
    861859            else:       #'constraint'
    862860                constr = map(list,zip([1.0 for i in range(len(varbs))],varbs))
    863 #                constr = map(dict,zip(varbs,[1.0 for i in range(len(varbs))]))
    864861                return [constr+[0.0,None]]          #just one constraint
    865862        return []
    … …  
    951948                constSizer.Add(constDel)
    952949                eqString = ' FIXED   '+item[0][1]+'   '
    953                 constSizer.Add((5,0),0)               
    954950            else:
    955951                constEdit = wx.Button(pageDisplay,-1,'Edit',style=wx.BU_EXACTFIT)
    … …  
    959955                constSizer.Add(constDel)
    960956                if isinstance(item[-1],bool):
    961                     constRef = wx.CheckBox(pageDisplay,-1,label=' Refine?')                   
    962                     constRef.SetValue(item[-1])
    963                     constRef.Bind(wx.EVT_CHECKBOX,OnConstRef)
    964                     Indx[constRef.GetId()] = item
    965                     constSizer.Add(constRef,0,wx.ALIGN_CENTER_VERTICAL)
    966957                    eqString = ' FUNCT   '
    967958                elif isinstance(item[-2],float):
    968                     constSizer.Add((5,5),0)
    969959                    eqString = ' CONSTR  '
    970960                else:
    971                     constSizer.Add((5,5),0)
    972961                    eqString = ' EQUIV   '
    973962                for term in item[:-2]:
    … …  
    978967                    eqString += ' = 0   '
    979968            constSizer.Add(wx.StaticText(pageDisplay,-1,eqString),0,wx.ALIGN_CENTER_VERTICAL)
     969            if isinstance(item[-1],bool):
     970                constRef = wx.CheckBox(pageDisplay,-1,label=' Refine?')                   
     971                constRef.SetValue(item[-1])
     972                constRef.Bind(wx.EVT_CHECKBOX,OnConstRef)
     973                Indx[constRef.GetId()] = item
     974                constSizer.Add(constRef,0,wx.ALIGN_CENTER_VERTICAL)
     975            else:
     976                constSizer.Add((5,5),0)
    980977        return constSizer
    981978               

  • trunk/GSASIImapvars.py

    r407 r411  
    9090     each group. This contains both parameters used in parameter
    9191     redefinitions as well as names of generated new parameters.
    92    varyList: a list of parameters that have been flagged to be varied.
    9392   
    9493   arrayList: a list by group of relationship matrices to relate
    … …  
    116115fixedDict = {} # a dictionary containing the fixed values corresponding to defined parameter equations
    117116               # key is original ascii string, value is float
    118 varyList = [] # a list of varied constraints
    119117
    120118# compile regular expressions used for parsing input
    … …  
    132130def InitVars():
    133131    '''Initializes all constraint information'''
    134     global dependentParmList,arrayList,invarrayList,indParmList,fixedDict,varyList,consNum
     132    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict,consNum
    135133    dependentParmList = [] # contains a list of parameters in each group
    136134    arrayList = [] # a list of of relationship matrices
    … …  
    138136    indParmList = [] # a list of names for the new parameters
    139137    fixedDict = {} # a dictionary containing the fixed values corresponding to defined parameter equations
    140     varyList = [] # a list of varied constraints
    141138    consNum = 0 # number of the next constraint to be created
    142139
    … …  
    301298    return groups,ParmList
    302299
    303 def GenerateConstraints(groups,parmlist,constrDict,constrFlag,fixedList):
     300def GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList):
    304301    '''Takes a list of relationship entries comprising a group of constraints and
    305302    builds the relationship lists and their inverse and stores them in global variables
    306303    '''
    307     global dependentParmList,arrayList,invarrayList,indParmList,varyList,consNum
     304    global dependentParmList,arrayList,invarrayList,indParmList,consNum
    308305    for group,varlist in zip(groups,parmlist):
    309306        VaryFree = False
    … …  
    337334    # key is original ascii string, value is float
    338335    for fixedval in fixedList:
    339         if fixedval is not None:
     336        if fixedval:
    340337            fixedDict[fixedval] = float(fixedval)
    341338
    … …  
    385382    return
    386383
    387 def VarRemapShow():
     384def VarRemapShow(varyList):
    388385    '''List out the saved relationships.
    389386    Returns a string containing the details.
    390387    '''
    391388    s = 'Mapping relations:\n'
    392     global dependentParmList,arrayList,invarrayList,indParmList,fixedDict,varyList
     389    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
    393390    for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList):
    394391        i = 0
    … …  
    443440                pass
    444441        if multarr is None: continue
    445         valuelist = [parmdict[var] for var in varlist]
     442        valuelist = [parmDict[var] for var in varlist]
    446443        parmDict.update(zip(mapvars,
    447444                            np.dot(multarr,np.array(valuelist)))
    … …  
    450447    parmDict.update(fixedDict)
    451448
    452 def Dict2Map(parmDict):
     449def Dict2Map(parmDict,varyList):
     450    #I think this needs fixing to update parmDict with new values
     451    #   from the constraints based on what is in varyList - RVD
    453452    '''Convert the remapped values back to the original parameters
    454453   
    … …  
    458457    '''
    459458    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
    460     # reset fixed values (should not be needed, but very quick)
     459    # reset fixed values (should not be needed, but very quick)
     460    # - this seems to update parmDict with {'0':0.0} & {'1.0':1.0} - probably not what was intended
    461461    parmDict.update(fixedDict)
    462462    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
    … …  
    675675
    676676    before = parmdict.copy()
    677     Dict2Map(parmdict)
     677    Dict2Map(parmdict,[])
    678678    print 'after Dict2Map'
    679679    print '  key / before / after'

  • trunk/GSASIIphsGUI.py

    r409 r411  
    173173                generalData['AngleRadii'].append(Info['Arad'])
    174174                generalData['vdWRadii'].append(Info['Vdrad'])
    175                 generalData['AtomMass'].append(Info['Mass'])
     175                if atom[ct] in generalData['Isotope']:
     176                    generalData['AtomMass'].append(Info['Isotopes'][generalData['Isotope'][atom[ct]]][0])
     177                else:
     178                    generalData['AtomMass'].append(Info['Mass'])
    176179                generalData['NoAtoms'][atom[ct]] = atom[cs-1]*float(atom[cs+1])
    177180                generalData['Color'].append(Info['Color'])
    … …  
    328331        def OnIsotope(event):
    329332            Obj = event.GetEventObject()
    330             generalData['Isotope'][Indx[Obj.GetId()]] = Obj.GetValue() #mass too
     333            item = Indx[Obj.GetId()]
     334            isotope = Obj.GetValue()
     335            generalData['Isotope'][item] = isotope
     336            indx = generalData['AtomTypes'].index(item)
     337            data['General']['AtomMass'][indx] = generalData['Isotopes'][item][isotope][0]
     338            UpdateGeneral()
    331339                                               
    332340        cellGUIlist = [[['m3','m3m'],4,zip([" Unit cell: a = "," Vol = "],["%.5f","%.3f"],[True,False],[0,0])],

  • trunk/GSASIIpwd.py

    r376 r411  
    698698   
    699699    Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
     700#    Df = pyd.pypsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
    700701    Df /= np.sum(Df)
    701702    return Df
    … …  
    704705   
    705706    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
     707#    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
    706708    sumDf = np.sum(Df)
    707709    return Df,dFdp,dFds,dFdg,dFdsh

  • trunk/GSASIIstruct.py

    r409 r411  
    6565   
    6666def GetConstraints(GPXfile):
    67     constrList = []
     67    constList = []
    6868    file = open(GPXfile,'rb')
    6969    while True:
    … …  
    7676            constDict = datum[1]
    7777            for item in constDict:
    78                 constrList += constDict[item]
     78                constList += constDict[item]
    7979    file.close()
    80     return constrList
     80    constDict = []
     81    constFlag = []
     82    fixedList = []
     83    for item in constList:
     84        if item[-2]:
     85            fixedList.append(str(item[-2]))
     86        else:
     87            fixedList.append('0')
     88        if item[-1]:
     89            constFlag.append(['VARY'])
     90        else:
     91            constFlag.append([])
     92        itemDict = {}
     93        for term in item[:-2]:
     94            itemDict[term[1]] = term[0]
     95        constDict.append(itemDict)
     96    return constDict,constFlag,fixedList
    8197   
    8298def GetPhaseNames(GPXfile):
    … …  
    392408    isotope = General['Isotope']
    393409    for El in atomTypes:
    394         BLtable[El] = isotopes[El][isotope[El]]
     410        BLtable[El] = [isotope[El],isotopes[El][isotope[El]]]
    395411    return BLtable
    396412       
    … …  
    452468           
    453469    def PrintFFtable(FFtable):
    454         print '\n Scattering factors:'
     470        print '\n X-ray scattering factors:'
    455471        print '   Symbol     fa                                      fb                                      fc'
    456472        print 99*'-'
    … …  
    461477            print ' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' %  \
    462478                (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
     479               
     480    def PrintBLtable(BLtable):
     481        print '\n Neutron scattering factors:'
     482        print '   Symbol   isotope       mass       b       resonant terms'
     483        print 99*'-'
     484        for Ename in BLtable:
     485            bldata = BLtable[Ename]
     486            isotope = bldata[0]
     487            mass = bldata[1][0]
     488            blen = bldata[1][1]
     489            bres = []
     490            if len(bldata[1]) > 2:
     491                bres = bldata[1][2:]
     492            line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
     493            for item in bres:
     494                line += '%10.5g'%(item)
     495            print line
    463496               
    464497    def PrintAtoms(General,Atoms):
    … …  
    485518       
    486519    def PrintTexture(textureData):
    487         print '\n Spherical harmonics texture: Order:' + \
    488             str(textureData['Order'])+' Refine? '+str(textureData['SH Coeff'][0])
     520        topstr = '\n Spherical harmonics texture: Order:' + \
     521            str(textureData['Order'])
     522        if textureData['Order']:
     523            print topstr+' Refine? '+str(textureData['SH Coeff'][0])
     524        else:
     525            print topstr
     526            return
    489527        names = ['omega','chi','phi']
    490528        line = '\n'
    … …  
    527565        except KeyError:
    528566            PawleyRef = []
    529         if Print: print '\n Phase name: ',General['Name']
    530567        SGData = General['SGData']
    531568        SGtext = G2spc.SGPrint(SGData)
    532         if Print:
    533             for line in SGtext: print line
    534569        cell = General['Cell']
    535570        A = G2lat.cell2A(cell[1:7])
    … …  
    537572        if cell[0]:
    538573            phaseVary += cellVary(pfx,SGData)
    539         if Print: print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
    540             ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
    541             '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
    542         if Print and FFtable:
    543             PrintFFtable(FFtable)
    544574        Natoms[pfx] = 0
    545575        if Atoms:
    … …  
    583613#            elif General['Type'] == 'macromolecular':
    584614
     615               
    585616            if 'SH Texture' in General:
    586617                textureData = General['SH Texture']
    … …  
    597628               
    598629            if Print:
     630                print '\n Phase name: ',General['Name']
     631                print 135*'-'
     632                PrintFFtable(FFtable)
     633                PrintBLtable(BLtable)
     634                print ''
     635                for line in SGtext: print line
    599636                PrintAtoms(General,Atoms)
     637                print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
     638                    ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
     639                    '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
    600640                if 'SH Texture' in General:
    601641                    PrintTexture(textureData)
    … …  
    15221562                keys[key][iatm] = parmDict[parm]
    15231563        FFdata.append(FFtables[Tdata[iatm]])
    1524         BLdata.append(BLtables[Tdata[iatm]])
     1564        BLdata.append(BLtables[Tdata[iatm]][1])
    15251565    return FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
    15261566   
    … …  
    15441584    FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict)
    15451585    if 'N' in parmDict[hfx+'Type']:
    1546         FP = 0.
    1547         FPP = 0.
     1586        FP,FPP = G2el.BlenRes(BLdata,parmDict[hfx+'Lam'])
    15481587    else:
    15491588        FP = np.array([El[hfx+'FP'] for El in FFdata])
    … …  
    16471686        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
    16481687        if not SGData['SGInv']:
    1649             dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)
     1688            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #problem here if occ=0 for some atom
    16501689            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)         
    16511690            dfbdui = np.sum(-SQfactor*fb,axis=2)
    … …  
    23302369def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
    23312370    parmdict.update(zip(varylist,values))
    2332     G2mv.Dict2Map(parmdict)
     2371    G2mv.Dict2Map(parmdict,varylist)
    23332372    Histograms,Phases = HistoPhases
    23342373    dMdv = np.empty(0)
    … …  
    23532392    parmdict.update(zip(varylist,values))
    23542393    Values2Dict(parmdict, varylist, values)
    2355     G2mv.Dict2Map(parmdict)
     2394    G2mv.Dict2Map(parmdict,varylist)
    23562395    Histograms,Phases = HistoPhases
    23572396    M = np.empty(0)
    … …  
    24052444    Controls = GetControls(GPXfile)
    24062445    ShowControls(Controls)           
    2407     constrList = GetConstraints(GPXfile)
     2446    constrDict,constrFlag,fixedList = GetConstraints(GPXfile)
    24082447    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
    24092448    if not Phases:
    … …  
    24282467    parmDict.update(histDict)
    24292468    GetFprime(calcControls,Histograms)
    2430     for item in constrList: print item
    2431     constrDict,constrFlag,fixedList = G2mv.InputParse([])        #constraints go here?
    24322469    groups,parmlist = G2mv.GroupConstraints(constrDict)
    2433     G2mv.GenerateConstraints(groups,parmlist,constrDict,constrFlag,fixedList)
     2470    G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList)
     2471    print G2mv.VarRemapShow(varyList)
    24342472    G2mv.Map2Dict(parmDict,varyList)
    24352473
    … …  
    24532491        chisq = np.sum(result[2]['fvec']**2)
    24542492        Values2Dict(parmDict, varyList, result[0])
    2455         G2mv.Dict2Map(parmDict)
     2493        G2mv.Dict2Map(parmDict,varyList)
    24562494        newAtomDict = ApplyXYZshifts(parmDict,varyList)
    24572495       
    … …  
    24802518                    break
    24812519
    2482 #    print 'dependentParmList: ',G2mv.dependentParmList
    2483 #    print 'arrayList: ',G2mv.arrayList
    2484 #    print 'invarrayList: ',G2mv.invarrayList
    2485 #    print 'indParmList: ',G2mv.indParmList
    2486 #    print 'fixedDict: ',G2mv.fixedDict
     2520    print 'dependentParmList: ',G2mv.dependentParmList
     2521    print 'arrayList: ',G2mv.arrayList
     2522    print 'invarrayList: ',G2mv.invarrayList
     2523    print 'indParmList: ',G2mv.indParmList
     2524    print 'fixedDict: ',G2mv.fixedDict
     2525    print 'test1'
    24872526    GetFobsSq(Histograms,Phases,parmDict,calcControls)
     2527    print 'test2'
    24882528    sigDict = dict(zip(varyList,sig))
    24892529    covData = {'variables':result[0],'varyList':varyList,'sig':sig,
    … …  
    24942534    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData)
    24952535#for testing purposes!!!
    2496 #    file = open('structTestdata.dat','wb')
    2497 #    cPickle.dump(parmDict,file,1)
    2498 #    cPickle.dump(varyList,file,1)
    2499 #    for histogram in Histograms:
    2500 #        if 'PWDR' in histogram[:4]:
    2501 #            Histogram = Histograms[histogram]
    2502 #    cPickle.dump(Histogram,file,1)
    2503 #    cPickle.dump(Phases,file,1)
    2504 #    cPickle.dump(calcControls,file,1)
    2505 #    cPickle.dump(pawleyLookup,file,1)
    2506 #    file.close()
     2536    file = open('structTestdata.dat','wb')
     2537    cPickle.dump(parmDict,file,1)
     2538    cPickle.dump(varyList,file,1)
     2539    for histogram in Histograms:
     2540        if 'PWDR' in histogram[:4]:
     2541            Histogram = Histograms[histogram]
     2542    cPickle.dump(Histogram,file,1)
     2543    cPickle.dump(Phases,file,1)
     2544    cPickle.dump(calcControls,file,1)
     2545    cPickle.dump(pawleyLookup,file,1)
     2546    file.close()
    25072547
    25082548def SeqRefine(GPXfile,dlg):
    … …  
    25162556    Controls = GetControls(GPXfile)
    25172557    ShowControls(Controls)           
    2518     constrList = GetConstraints(GPXfile)
     2558    constrDict,constrFlag,fixedList = GetConstraints(GPXfile)
    25192559    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
    25202560    if not Phases:
    … …  
    25802620        constrDict,constrFlag,fixedList = G2mv.InputParse([])        #constraints go here?
    25812621        groups,parmlist = G2mv.GroupConstraints(constrDict)
    2582         G2mv.GenerateConstraints(groups,parmlist,constrDict,constrFlag,fixedList)
     2622        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList)
    25832623        G2mv.Map2Dict(parmDict,varyList)
    25842624   
    … …  
    26002640            chisq = np.sum(result[2]['fvec']**2)
    26012641            Values2Dict(parmDict, varyList, result[0])
    2602             G2mv.Dict2Map(parmDict)
     2642            G2mv.Dict2Map(parmDict,varyList)
    26032643            newAtomDict = ApplyXYZshifts(parmDict,varyList)
    26042644           

  • trunk/fsource/powsubs/psvfcj.for

    r353 r411  
    126126c
    127127      IF (A .ne. 0.0) then
    128         Einfl = Acosd(cos2THETA) ! 2phi(infl) FCJ eq 5 (degrees)
     128        Einfl = Acosd(cos2THETA)
    129129        tmp2 = 1.0 + ApB2
    130130        tmp = SQRT(tmp2)*cos2THETA
    131 c
    132 C Treat case where A is zero - Set Einfl = 2theta
    133 c
    134         if (A.eq.0.0) Einfl = Acosd(cos2THETA)
    135131        if (abs(tmp) .le. 1.0) then
    136           Emin = Acosd(tmp)      ! 2phi(min) FCJ eq 4 (degrees)
     132          Emin = Acosd(tmp)
    137133          tmp1 = tmp2*(1.0 - tmp2*(1.0-sin2THETA2))
    138134        else
     135            print *,'tmp > 1.0'
    139136          tmp1 = 0.0
    140137          if (tmp .gt. 0.0) then
    … …  
    213210          tanDELTA = tand(Delta)
    214211          cosDELTA2 = cosDELTA*cosDELTA     
    215           tmp1 = cosDELTA2 - cos2THETA2
    216           tmp2 = sin2THETA2 - sinDelta * sinDELTA
    217           tmp = tmp2
    218           if ( ttheta.gt.4500.0 ) tmp = tmp1
     212          if ( ttheta.le.4500.0 .or. ttheta.gt.13500.0 ) then
     213            tmp = sin2THETA2-sinDelta**2
     214          else
     215            tmp = cosDELTA2-cos2THETA2
     216          end if
    219217          if (tmp .gt. 0) then
    220218            tmp1 = 1.0/SQRT(tmp)
    221219            F = abs(cos2THETA) * tmp1
    222             dFdA = cosDELTA*cos2THETA*sinDELTA*dDELTAdA *
    223      1           (tmp1*tmp1*tmp1)
     220            dFdA = cosDELTA*cos2THETA*sinDELTA*dDELTAdA*tmp1**3
    224221          else
    225222            F = 1.0
    … …  
    228225c       
    229226c Calculate G(Delta,2theta) [G = W /(h cos(delta) ] [ FCJ eq. 7(a) and 7(b) ]
    230 c       
    231           if(abs(delta-emin) .gt. abs(einfl-emin))then
    232 c
    233 C N.B. this is the only place where d()/dA <> d()/dB
    234 c
     227c
     228          if ( ttheta.gt.9000.0 ) then                  !are you kidding!! this worked
     229!          if(abs(delta-emin) .gt. abs(einfl-emin))then
    235230              G = 2.0*A*F*RcosDELTA
    236231              dGdA = 2.0*A*RcosDELTA*(dFdA + F*tanDELTA*dDELTAdA)
    237           else                                            ! delta .le. einfl .or. min(A,B) .eq. 0
     232          else
    238233            G = (-1.0 + ApB*F) * RcosDELTA
    239234            dGdA = RcosDELTA*(F - tanDELTA*dDELTAdA
Note: See TracChangeset for help on using the changeset viewer.