Changeset 530 for trunk/GSASIIphsGUI.py
 Timestamp:
 Apr 6, 2012 1:30:45 PM (11 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIphsGUI.py
r529 r530 25 25 import GSASIIIO as G2IO 26 26 import GSASIIstruct as G2str 27 import GSASIImath as G2mth 27 28 import numpy as np 28 29 import numpy.linalg as nl … … 3445 3446 data['Drawing']['mapSize'] = 10. 3446 3447 print mapData['MapType']+' computed: rhomax = %.3f rhomin = %.3f'%(np.max(mapData['rho']),np.min(mapData['rho'])) 3447 G2plt.PlotStructure(G2frame,data) 3448 ## map printing for testing purposes 3449 # ix,jy,kz = mapData['rho'].shape 3450 # for k in range(kz): 3451 # print 'k = ',k 3452 # for j in range(jy): 3453 # line = '' 3454 # if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']: 3455 # line += (jyj)*' ' 3456 # for i in range(ix): 3457 # r = int(100*mapData['rho'][i,j,k]/mapData['rhoMax']) 3458 # line += '%4d'%(r) 3459 # print line+'\n' 3460 ### keep this 3448 G2plt.PlotStructure(G2frame,data) 3449 3450 def printRho(SGData,rho,rhoMax): 3451 # map printing for testing purposes 3452 ix,jy,kz = rho.shape 3453 for k in range(kz): 3454 print 'k = ',k 3455 for j in range(jy): 3456 line = '' 3457 if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']: 3458 line += (jyj)*' ' 3459 for i in range(ix): 3460 r = int(100*rho[i,j,k]/rhoMax) 3461 line += '%4d'%(r) 3462 print line+'\n' 3463 ## keep this 3461 3464 3462 3465 def OnSearchMaps(event): 3463 3466 3467 norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3) 3468 3469 def noDuplicate(xyz,peaks,SGData): #be careful where this is used  it's slow 3470 equivs = G2spc.GenAtom(xyz,SGData) 3471 xyzs = [equiv[0] for equiv in equivs] 3472 for x in xyzs: 3473 if True in [np.allclose(x,peak,atol=0.002) for peak in peaks]: 3474 return False 3475 return True 3476 3464 3477 def findRoll(rhoMask,mapHalf): 3465 3478 indx = np.array(ma.nonzero(rhoMask)).T … … 3467 3480 rhoMax = np.max(rhoList) 3468 3481 return mapHalfindx[np.argmax(rhoList)] 3482 3483 def rhoCalc(parms,rX,rY,rZ,res,SGLaue): 3484 Mag,x0,y0,z0,sig = parms 3485 if SGLaue in ['3','3m1','31m','6/m','6/mmm']: 3486 return norm*Mag*np.exp(((x0rX)**2+(y0rY)**2+(x0rX)*(y0rY)+(z0rZ)**2)/(2.*sig**2))/(sig*res**3) 3487 else: 3488 return norm*Mag*np.exp(((x0rX)**2+(y0rY)**2+(z0rZ)**2)/(2.*sig**2))/(sig*res**3) 3489 3490 def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue): 3491 Mag,x0,y0,z0,sig = parms 3492 M = rhorhoCalc(parms,rX,rY,rZ,res,SGLaue) 3493 return M 3494 3495 def peakHess(parms,rX,rY,rZ,rho,res,SGLaue): 3496 Mag,x0,y0,z0,sig = parms 3497 dMdv = np.zeros(([5,]+list(rX.shape))) 3498 delt = .01 3499 for i in range(5): 3500 parms[i] = delt 3501 rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue) 3502 parms[i] += 2.*delt 3503 rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue) 3504 parms[i] = delt 3505 dMdv[i] = (rhoCprhoCm)/(2.*delt) 3506 rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue) 3507 Vec = np.sum(np.sum(np.sum(dMdv*(rhorhoC),axis=3),axis=2),axis=1) 3508 dMdv = np.reshape(dMdv,(5,rX.size)) 3509 Hess = np.inner(dMdv,dMdv) 3510 3511 return Vec,Hess 3469 3512 3470 3513 generalData = data['General'] … … 3479 3522 rho = copy.copy(mapData['rho']) #don't mess up original 3480 3523 mapHalf = np.array(rho.shape)/2 3481 step = max(1.0,1./mapData['Resolution']) 3482 steps = 3*np.array(3*[step,]) 3524 res = mapData['Resolution'] 3525 incre = 1./np.array(rho.shape) 3526 step = max(1.0,1./res)+1 3527 steps = np.array(3*[step,]) 3483 3528 except KeyError: 3484 3529 print '**** ERROR  Fourier map not defined' 3485 3530 return 3486 3531 peaks = [] 3487 while True: 3488 rhoMask = ma.array(rho,mask=(rho<contLevel)) 3489 rMI = findRoll(rhoMask,mapHalf) 3490 rho = np.roll(np.roll(np.roll(rho,rMI[0],axis=0),rMI[1],axis=1),rMI[2],axis=2) 3491 rMM = mapHalfsteps 3492 rMP = mapHalf+steps 3493 3494 3495 rho = np.roll(np.roll(np.roll(rho,rMI[2],axis=2),rMI[1],axis=1),rMI[0],axis=0) 3496 break 3497 print 'search Fourier map' 3532 mags = [] 3533 print ' Begin fourier map search  can take some time' 3534 time0 = time.time() 3535 wx.BeginBusyCursor() 3536 try: 3537 while True: 3538 rhoMask = ma.array(rho,mask=(rho<contLevel)) 3539 if not ma.count(rhoMask): 3540 break 3541 rMI = findRoll(rhoMask,mapHalf) 3542 rho = np.roll(np.roll(np.roll(rho,rMI[0],axis=0),rMI[1],axis=1),rMI[2],axis=2) 3543 rMM = mapHalfsteps 3544 rMP = mapHalf+steps+1 3545 rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]] 3546 peakInt = np.sum(rhoPeak)*res**3 3547 rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]] 3548 x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0] #magnitude, position & width(sig) 3549 result = G2mth.HessianLSQ(peakFunc,x0,Hess=peakHess, 3550 args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.0001,maxcyc=100) 3551 x1 = result[0] 3552 if np.any(x1 < 0): 3553 break 3554 peak = (np.array(x1[1:4])rMI)*incre 3555 if not len(peaks): 3556 peaks.append(peak) 3557 mags.append(x1[0]) 3558 else: 3559 if noDuplicate(peak,peaks,SGData) and x1[0] > 0.: 3560 peaks.append(peak) 3561 mags.append(x1[0]) 3562 rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]] = peakFunc(result[0],rX,rY,rZ,rhoPeak,res,SGData['SGLaue']) 3563 rho = np.roll(np.roll(np.roll(rho,rMI[2],axis=2),rMI[1],axis=1),rMI[0],axis=0) 3564 finally: 3565 wx.EndBusyCursor() 3566 sortIdx = np.argsort(mags) 3567 print ' Map search peaks found:' 3568 print ' No. Mag. x y z' 3569 for j,i in enumerate(sortIdx[::1]): 3570 print ' %3d %8.3f %8.5f %8.5f %8.5f'%(j,mags[i],peaks[i][0],peaks[i][1],peaks[i][2]) 3571 print ' Map search finished, time = %.2fs'%(time.time()time0) 3498 3572 3499 3573 def OnTextureRefine(event):
Note: See TracChangeset
for help on using the changeset viewer.