Changeset 647


Ignore:
Timestamp:
Jun 5, 2012 11:03:48 AM (9 years ago)
Author:
vondreele
Message:

HKLF mult fix
finish findOffset - works & reasonable speed

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r646 r647  
    612612    DH = []
    613613    Dphi = []
    614     while float(F) > 0.5*Fmax:
     614    while i < 20:
    615615        F = Flist[i]
    616616        hkl = np.unravel_index(Fdict[F],hklShape)
     
    626626            dH = H-hkl
    627627            dang = ang-ang0
     628            if np.any(np.abs(dH)-8 > 0):    #keep low order DHs
     629                continue
    628630            DH.append(dH)
    629631            Dphi.append((dang+0.5) % 1.0)
    630632        i += 1
    631 #    for item in zip(DH,Dphi):
    632 #        print item
    633633    DH = np.array(DH)
    634634    Dphi = np.array(Dphi)
     635#    for item in zip(DH,Dphi):
     636#        print item[0],'%.4f'%(item[1])
    635637    DX = np.zeros(3)
    636     X,Y,Z = np.mgrid[0:1:steps[0]*(0+1j),0:1:steps[1]*(0+1j),0:1:steps[2]*(0+1j)]
     638    X,Y,Z = np.mgrid[0:1:10j,0:1:10j,0:1:10j]
    637639    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
    638640    Mmin = 1.e10
    639641   
    640     for xyz in XYZ:
     642    for xyz in XYZ:             #do a global search for best roll
    641643        M = np.sum(calcPhase(xyz,DH,Dphi)**2)
    642644        if M < Mmin:
     
    644646            Mmin = M
    645647   
    646     print DX,Mmin
    647648    result = so.leastsq(calcPhase,DX,full_output=True,args=(DH,Dphi))
     649#    for item in zip(DH,Dphi,result[2]['fvec']):
     650#        print item[0],'%.4f %.4f'%(item[1],item[2])
    648651    chisq = np.sum(result[2]['fvec']**2)
    649652    DX = np.array(np.rint(-result[0]*steps),dtype='i')
    650     print chisq,DX
    651     print result[0],result[0]*steps
     653    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
    652654    return DX
    653655
     
    725727            break
    726728    np.seterr(**old)
    727     print 'Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
    728     print 'No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')
     729    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
     730    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')
    729731    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
    730732    roll = findOffset(SGData,CErho,CEhkl)
     
    824826        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
    825827        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
    826             args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.0001,maxcyc=100)
     828            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
    827829        x1 = result[0]
    828830        if np.any(x1 < 0):
     
    839841                peaks.append(peak)
    840842                mags.append(x1[0])
    841             if len(peaks) > 100:
     843            if len(peaks) > 300:
    842844                break
    843         rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]] = peakFunc(result[0],rX,rY,rZ,rhoPeak,res,SGData['SGLaue'])
     845        rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]] = peakFunc(x1,rX,rY,rZ,rhoPeak,res,SGData['SGLaue'])
    844846        rho = np.roll(np.roll(np.roll(rho,-rMI[2],axis=2),-rMI[1],axis=1),-rMI[0],axis=0)
    845847    return np.array(peaks),np.array([mags,]).T
  • trunk/GSASIIphsGUI.py

    r646 r647  
    34973497            H = ref[:3]
    34983498            ref[4] = np.sqrt(1./G2lat.calc_rDsq2(H,G))
    3499             iabsnt,mulp,Uniq,phi = G2spc.GenHKLf(H,SGData)
    3500             ref[3] = mulp/2             #convert from powder mulp.
    3501             ref[11] = Uniq
    3502             ref[12] = phi
     3499            iabsnt,ref[3],ref[11],ref[12] = G2spc.GenHKLf(H,SGData)
    35033500        G2frame.PatternTree.SetItemPyData(Id,[histoName,reflData])
    35043501       
Note: See TracChangeset for help on using the changeset viewer.