 Timestamp:
 May 18, 2021 1:22:50 PM (4 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIconstrGUI.py
r4906 r4907 2279 2279 2280 2280 def onSetPlane(event): 2281 '''Compute leastsquares plane for selected atoms; 2282 move atoms so that LS plane aligned with xy plane, 2283 with minimum required change to x 2281 '''Finds eigen vector/matrix for best "ellipsoid" about atoms; 2282 rotate atoms so that smallest axis is along choice. 2284 2283 ''' 2285 2284 grid.completeEdits() 2286 2285 selList = [i==1 for i in rd.Phase['RBselection']] 2287 2286 XYZ = rd.Phase['RBcoords'][selList] 2288 if len(XYZ) < 3:2289 G2G.G2MessageBox(G2frame,'A plane requires three or more atoms',2290 2287 Natoms = len(XYZ) 2288 if Natoms < 3: 2289 G2G.G2MessageBox(G2frame,'A plane requires three or more atoms','Need more atoms') 2291 2290 return 2292 # fit 3 ways (in case of singularity) and take result with lowest residual 2293 X,Y,Z = [XYZ[:,i] for i in (0,1,2)] 2294 XZ = copy.copy(XYZ) 2295 XZ[:,1] = 1 2296 (a,d,b), resd, rank, sing = nl.lstsq(XZ, Y,rcond=1) 2297 resid_min = resd 2298 normal = a,1,b 2299 YZ = copy.copy(XYZ) 2300 YZ[:,0] = 1 2301 (d,a,b), resd, rank, sing = nl.lstsq(YZ, X,rcond=1) 2302 if resid_min > resd: 2303 resid_min = resd 2304 normal = 1,a,b 2305 XY = copy.copy(XYZ) 2306 XY[:,2] = 1 2307 (a,b,d), resd, rank, sing = nl.lstsq(XY, Z,rcond=1) 2308 if resid_min > resd: 2309 resid_min = resd 2310 normal = a,b,1 2311 # solve for ax + bx + z + c = 0 or equivalently ax + bx + c = z 2312 # try: 2313 # except: 2314 # G2G.G2MessageBox(G2frame, 2315 # 'Error computing plane; are atoms in a line?', 2316 # 'Computation error') 2317 # return 2318 if bntOpts['plane'] == 'xy': 2319 # new coordinate system is 2320 # ZP, z' normal to plane 2321 # YP, y' = z' cross x (= [0,1,b]) 2322 # XP, x' = (z' cross x) cross z' 2323 # this puts XP as close as possible to X with XP & YP in plane 2324 ZP = np.array(normal) 2325 ZP /= np.sqrt(np.sum(ZP**2)) 2326 YP = np.cross(ZP,[1,0,0]) 2327 if sum(YP*YP) < .1: # pathological condition: z' along x 2328 YP = np.cross(ZP,[0,1,0]) 2329 YP /= np.sqrt(np.sum(YP**2)) 2330 XP = np.cross(YP,ZP) 2331 elif bntOpts['plane'] == 'yz': 2332 # new coordinate system is 2333 # XP, x' normal to plane 2334 # ZP, z' = x' cross y 2335 # YP, y' = (x' cross y) cross x' 2336 # this puts y' as close as possible to y with z' & y' in plane 2337 XP = np.array(normal) 2338 XP /= np.sqrt(np.sum(XP**2)) 2339 ZP = np.cross(XP,[0,1,0]) 2340 if sum(ZP*ZP) < .1: # pathological condition: x' along y 2341 ZP = np.cross(XP,(0,0,1)) 2342 ZP /= np.sqrt(np.sum(ZP**2)) 2343 YP = np.cross(ZP,XP) 2344 elif bntOpts['plane'] == 'xz': 2345 # new coordinate system is 2346 # YP, y' normal to plane 2347 # ZP, z' = x cross y' 2348 # XP, y' = (x cross y') cross z' 2349 # this puts XP as close as possible to X with XP & YP in plane 2350 YP = np.array(normal) 2351 YP /= np.sqrt(np.sum(YP**2)) 2352 ZP = np.cross([1,0,0],YP) 2353 if sum(ZP*ZP) < .1: # pathological condition: y' along x 2354 ZP = np.cross([0,1,0],YP) 2355 ZP /= np.sqrt(np.sum(ZP**2)) 2356 XP = np.cross(YP,ZP) 2291 Zmat = np.zeros((3,3)) 2292 for xyz in XYZ: 2293 Zmat += np.outer(xyz.T,xyz) 2294 Evec,Emat = nl.eig(Zmat) 2295 Order = np.argsort(np.nan_to_num(Evec)) #shortlong order 2296 if bntOpts['plane'] == 'xy': #short along z 2297 trans = np.array([Emat[Order[2]],Emat[Order[1]],Emat[Order[0]]]) 2298 elif bntOpts['plane'] == 'yz': #short along x 2299 trans = np.array([Emat[Order[0]],Emat[Order[2]],Emat[Order[1]]]) 2300 elif bntOpts['plane'] == 'xz': #short along y 2301 trans = np.array([Emat[Order[1]],Emat[Order[0]],Emat[Order[2]]]) 2357 2302 else: 2358 2303 print('unexpected plane',bntOpts['plane']) 2359 2304 return 2360 trans = np.array((XP,YP,ZP))2361 2305 # update atoms in place 2362 2306 rd.Phase['RBcoords'][:] = np.inner(trans,rd.Phase['RBcoords']).T
Note: See TracChangeset
for help on using the changeset viewer.