Changeset 1443 for trunk/GSASIIpwd.py


Ignore:
Timestamp:
Jul 28, 2014 4:53:56 PM (7 years ago)
Author:
vondreele
Message:

add calibration of lam, difC, etc. from index peak positions
new plotCalib routine

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIpwd.py

    r1439 r1443  
    11231123    backVary += peaksVary   
    11241124    return bakType,backDict,backVary
     1125   
     1126def DoCalibInst(IndexPeaks,Inst):
     1127   
     1128    def SetInstParms():
     1129        dataType = Inst['Type'][0]
     1130        insVary = []
     1131        insNames = []
     1132        insVals = []
     1133        for parm in Inst:
     1134            insNames.append(parm)
     1135            insVals.append(Inst[parm][1])
     1136            if parm in ['Lam','difC','difA','difB','Zero',]:
     1137                if Inst[parm][2]:
     1138                    insVary.append(parm)
     1139        instDict = dict(zip(insNames,insVals))
     1140        return dataType,instDict,insVary
     1141       
     1142    def GetInstParms(parmDict,Inst,varyList):
     1143        for name in Inst:
     1144            Inst[name][1] = parmDict[name]
     1145       
     1146    def InstPrint(Inst,sigDict):
     1147        print 'Instrument Parameters:'
     1148        if 'C' in Inst['Type'][0]:
     1149            ptfmt = "%12.6f"
     1150        else:
     1151            ptfmt = "%12.3f"
     1152        ptlbls = 'names :'
     1153        ptstr =  'values:'
     1154        sigstr = 'esds  :'
     1155        for parm in Inst:
     1156            if parm in  ['Lam','difC','difA','difB','Zero',]:
     1157                ptlbls += "%s" % (parm.center(12))
     1158                ptstr += ptfmt % (Inst[parm][1])
     1159                if parm in sigDict:
     1160                    sigstr += ptfmt % (sigDict[parm])
     1161                else:
     1162                    sigstr += 12*' '
     1163        print ptlbls
     1164        print ptstr
     1165        print sigstr
     1166       
     1167    def errPeakPos(values,peakDsp,peakPos,peakWt,dataType,parmDict,varyList):
     1168        parmDict.update(zip(varyList,values))
     1169        return np.sqrt(peakWt)*(G2lat.getPeakPos(dataType,parmDict,peakDsp)-peakPos)
     1170
     1171    peakPos = []
     1172    peakDsp = []
     1173    peakWt = []
     1174    for peak in IndexPeaks:
     1175        if peak[2] and peak[3]:
     1176            peakPos.append(peak[0])
     1177            peakDsp.append(peak[8])
     1178            peakWt.append(1/peak[1])
     1179    peakPos = np.array(peakPos)
     1180    peakDsp = np.array(peakDsp)
     1181    peakWt = np.array(peakWt)
     1182    dataType,insDict,insVary = SetInstParms()
     1183    parmDict = {}
     1184    parmDict.update(insDict)
     1185    varyList = insVary
     1186    while True:
     1187        begin = time.time()
     1188        values =  np.array(Dict2Values(parmDict, varyList))
     1189        result = so.leastsq(errPeakPos,values,full_output=True,ftol=0.0001,
     1190            args=(peakDsp,peakPos,peakWt,dataType,parmDict,varyList))
     1191        ncyc = int(result[2]['nfev']/2)
     1192        runtime = time.time()-begin   
     1193        chisq = np.sum(result[2]['fvec']**2)
     1194        Values2Dict(parmDict, varyList, result[0])
     1195        GOF = chisq/(len(peakPos)-len(varyList))       #reduced chi^2
     1196        print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',len(peakPos),' Number of parameters: ',len(varyList)
     1197        print 'calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
     1198        print 'chi**2 = %12.6g, reduced chi**2 = %6.2f'%(chisq,GOF)
     1199        try:
     1200            sig = np.sqrt(np.diag(result[1])*GOF)
     1201            if np.any(np.isnan(sig)):
     1202                print '*** Least squares aborted - some invalid esds possible ***'
     1203            break                   #refinement succeeded - finish up!
     1204        except ValueError:          #result[1] is None on singular matrix
     1205            print '**** Refinement failed - singular matrix ****'
     1206        return
     1207       
     1208    sigDict = dict(zip(varyList,sig))
     1209    GetInstParms(parmDict,Inst,varyList)
     1210    InstPrint(Inst,sigDict)
    11251211           
    11261212def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,oneCycle=False,controls=None,dlg=None):
     
    14111497        elif FitPgm == 'BFGS':
    14121498            print 'Other program here'
    1413             return
     1499            return {}
    14141500       
    14151501    sigDict = dict(zip(varyList,sig))
     
    14231509    GetPeaksParms(Inst,parmDict,Peaks,varyList)   
    14241510    PeaksPrint(dataType,parmDict,sigDict,varyList)
     1511    return sigDict
    14251512
    14261513def calcIncident(Iparm,xdata):
Note: See TracChangeset for help on using the changeset viewer.