| 1 | # -*- coding: utf-8 -*- |
|---|
| 2 | #GSASII cell indexing program: variation on that of A. Coehlo |
|---|
| 3 | # includes cell refinement from peak positions (not zero as yet) |
|---|
| 4 | ########### SVN repository information ################### |
|---|
| 5 | # $Date: 2017-01-20 16:22:57 +0000 (Fri, 20 Jan 2017) $ |
|---|
| 6 | # $Author: toby $ |
|---|
| 7 | # $Revision: 2650 $ |
|---|
| 8 | # $URL: trunk/GSASIIindex.py $ |
|---|
| 9 | # $Id: GSASIIindex.py 2650 2017-01-20 16:22:57Z toby $ |
|---|
| 10 | ########### SVN repository information ################### |
|---|
| 11 | ''' |
|---|
| 12 | *GSASIIindex: Cell Indexing Module* |
|---|
| 13 | =================================== |
|---|
| 14 | |
|---|
| 15 | Cell indexing program: variation on that of A. Coehlo |
|---|
| 16 | includes cell refinement from peak positions |
|---|
| 17 | ''' |
|---|
| 18 | |
|---|
| 19 | import math |
|---|
| 20 | import time |
|---|
| 21 | import numpy as np |
|---|
| 22 | import GSASIIpath |
|---|
| 23 | GSASIIpath.SetVersionNumber("$Revision: 2650 $") |
|---|
| 24 | import GSASIIlattice as G2lat |
|---|
| 25 | import GSASIIpwd as G2pwd |
|---|
| 26 | import GSASIIspc as G2spc |
|---|
| 27 | import GSASIImath as G2mth |
|---|
| 28 | import scipy.optimize as so |
|---|
| 29 | |
|---|
| 30 | # trig functions in degrees |
|---|
| 31 | sind = lambda x: math.sin(x*math.pi/180.) |
|---|
| 32 | asind = lambda x: 180.*math.asin(x)/math.pi |
|---|
| 33 | tand = lambda x: math.tan(x*math.pi/180.) |
|---|
| 34 | atand = lambda x: 180.*math.atan(x)/math.pi |
|---|
| 35 | atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi |
|---|
| 36 | cosd = lambda x: math.cos(x*math.pi/180.) |
|---|
| 37 | acosd = lambda x: 180.*math.acos(x)/math.pi |
|---|
| 38 | rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p) |
|---|
| 39 | #numpy versions |
|---|
| 40 | npsind = lambda x: np.sin(x*np.pi/180.) |
|---|
| 41 | npasind = lambda x: 180.*np.arcsin(x)/math.pi |
|---|
| 42 | npcosd = lambda x: np.cos(x*math.pi/180.) |
|---|
| 43 | nptand = lambda x: np.tan(x*math.pi/180.) |
|---|
| 44 | npatand = lambda x: 180.*np.arctan(x)/np.pi |
|---|
| 45 | npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi |
|---|
| 46 | rpd = np.pi/180. |
|---|
| 47 | |
|---|
| 48 | def scaleAbyV(A,V): |
|---|
| 49 | 'needs a doc string' |
|---|
| 50 | v = G2lat.calc_V(A) |
|---|
| 51 | scale = math.exp(math.log(v/V)/3.)**2 |
|---|
| 52 | for i in range(6): |
|---|
| 53 | A[i] *= scale |
|---|
| 54 | |
|---|
| 55 | def ranaxis(dmin,dmax): |
|---|
| 56 | 'needs a doc string' |
|---|
| 57 | import random as rand |
|---|
| 58 | return rand.random()*(dmax-dmin)+dmin |
|---|
| 59 | |
|---|
| 60 | def ran2axis(k,N): |
|---|
| 61 | 'needs a doc string' |
|---|
| 62 | import random as rand |
|---|
| 63 | T = 1.5+0.49*k/N |
|---|
| 64 | # B = 0.99-0.49*k/N |
|---|
| 65 | # B = 0.99-0.049*k/N |
|---|
| 66 | B = 0.99-0.149*k/N |
|---|
| 67 | R = (T-B)*rand.random()+B |
|---|
| 68 | return R |
|---|
| 69 | |
|---|
| 70 | #def ranNaxis(k,N): |
|---|
| 71 | # import random as rand |
|---|
| 72 | # T = 1.0+1.0*k/N |
|---|
| 73 | # B = 1.0-1.0*k/N |
|---|
| 74 | # R = (T-B)*rand.random()+B |
|---|
| 75 | # return R |
|---|
| 76 | |
|---|
| 77 | def ranAbyV(Bravais,dmin,dmax,V): |
|---|
| 78 | 'needs a doc string' |
|---|
| 79 | cell = [0,0,0,0,0,0] |
|---|
| 80 | bad = True |
|---|
| 81 | while bad: |
|---|
| 82 | bad = False |
|---|
| 83 | cell = rancell(Bravais,dmin,dmax) |
|---|
| 84 | G,g = G2lat.cell2Gmat(cell) |
|---|
| 85 | A = G2lat.Gmat2A(G) |
|---|
| 86 | if G2lat.calc_rVsq(A) < 1: |
|---|
| 87 | scaleAbyV(A,V) |
|---|
| 88 | cell = G2lat.A2cell(A) |
|---|
| 89 | for i in range(3): |
|---|
| 90 | bad |= cell[i] < dmin |
|---|
| 91 | return A |
|---|
| 92 | |
|---|
| 93 | def ranAbyR(Bravais,A,k,N,ranFunc): |
|---|
| 94 | 'needs a doc string' |
|---|
| 95 | R = ranFunc(k,N) |
|---|
| 96 | if Bravais in [0,1,2]: #cubic - not used |
|---|
| 97 | A[0] = A[1] = A[2] = A[0]*R |
|---|
| 98 | A[3] = A[4] = A[5] = 0. |
|---|
| 99 | elif Bravais in [3,4]: #hexagonal/trigonal |
|---|
| 100 | A[0] = A[1] = A[3] = A[0]*R |
|---|
| 101 | A[2] *= R |
|---|
| 102 | A[4] = A[5] = 0. |
|---|
| 103 | elif Bravais in [5,6]: #tetragonal |
|---|
| 104 | A[0] = A[1] = A[0]*R |
|---|
| 105 | A[2] *= R |
|---|
| 106 | A[3] = A[4] = A[5] = 0. |
|---|
| 107 | elif Bravais in [7,8,9,10]: #orthorhombic |
|---|
| 108 | A[0] *= R |
|---|
| 109 | A[1] *= R |
|---|
| 110 | A[2] *= R |
|---|
| 111 | A[3] = A[4] = A[5] = 0. |
|---|
| 112 | elif Bravais in [11,12]: #monoclinic |
|---|
| 113 | A[0] *= R |
|---|
| 114 | A[1] *= R |
|---|
| 115 | A[2] *= R |
|---|
| 116 | A[4] *= R |
|---|
| 117 | A[3] = A[5] = 0. |
|---|
| 118 | else: #triclinic |
|---|
| 119 | A[0] *= R |
|---|
| 120 | A[1] *= R |
|---|
| 121 | A[2] *= R |
|---|
| 122 | A[3] *= R |
|---|
| 123 | A[4] *= R |
|---|
| 124 | A[5] *= R |
|---|
| 125 | return A |
|---|
| 126 | |
|---|
| 127 | def rancell(Bravais,dmin,dmax): |
|---|
| 128 | 'needs a doc string' |
|---|
| 129 | if Bravais in [0,1,2]: #cubic |
|---|
| 130 | a = b = c = ranaxis(dmin,dmax) |
|---|
| 131 | alp = bet = gam = 90 |
|---|
| 132 | elif Bravais in [3,4]: #hexagonal/trigonal |
|---|
| 133 | a = b = ranaxis(dmin,dmax) |
|---|
| 134 | c = ranaxis(dmin,dmax) |
|---|
| 135 | alp = bet = 90 |
|---|
| 136 | gam = 120 |
|---|
| 137 | elif Bravais in [5,6]: #tetragonal |
|---|
| 138 | a = b = ranaxis(dmin,dmax) |
|---|
| 139 | c = ranaxis(dmin,dmax) |
|---|
| 140 | alp = bet = gam = 90 |
|---|
| 141 | elif Bravais in [7,8,9,10]: #orthorhombic - F,I,C,P - a<b<c convention |
|---|
| 142 | abc = [ranaxis(dmin,dmax),ranaxis(dmin,dmax),ranaxis(dmin,dmax)] |
|---|
| 143 | if Bravais in [7,8,10]: |
|---|
| 144 | abc.sort() |
|---|
| 145 | a = abc[0] |
|---|
| 146 | b = abc[1] |
|---|
| 147 | c = abc[2] |
|---|
| 148 | alp = bet = gam = 90 |
|---|
| 149 | elif Bravais in [11,12]: #monoclinic - C,P - a<c convention |
|---|
| 150 | ac = [ranaxis(dmin,dmax),ranaxis(dmin,dmax)] |
|---|
| 151 | if Bravais == 12: |
|---|
| 152 | ac.sort() |
|---|
| 153 | a = ac[0] |
|---|
| 154 | b = ranaxis(dmin,dmax) |
|---|
| 155 | c = ac[1] |
|---|
| 156 | alp = gam = 90 |
|---|
| 157 | bet = ranaxis(90.,140.) |
|---|
| 158 | else: #triclinic - a<b<c convention |
|---|
| 159 | abc = [ranaxis(dmin,dmax),ranaxis(dmin,dmax),ranaxis(dmin,dmax)] |
|---|
| 160 | abc.sort() |
|---|
| 161 | a = abc[0] |
|---|
| 162 | b = abc[1] |
|---|
| 163 | c = abc[2] |
|---|
| 164 | r = 0.5*b/c |
|---|
| 165 | alp = ranaxis(acosd(r),acosd(-r)) |
|---|
| 166 | r = 0.5*a/c |
|---|
| 167 | bet = ranaxis(acosd(r),acosd(-r)) |
|---|
| 168 | r = 0.5*a/b |
|---|
| 169 | gam = ranaxis(acosd(r),acosd(-r)) |
|---|
| 170 | return [a,b,c,alp,bet,gam] |
|---|
| 171 | |
|---|
| 172 | def calc_M20(peaks,HKL,ifX20=True): |
|---|
| 173 | 'needs a doc string' |
|---|
| 174 | diff = 0 |
|---|
| 175 | X20 = 0 |
|---|
| 176 | for Nobs20,peak in enumerate(peaks): |
|---|
| 177 | if peak[3]: |
|---|
| 178 | Qobs = 1.0/peak[7]**2 |
|---|
| 179 | Qcalc = 1.0/peak[8]**2 |
|---|
| 180 | diff += abs(Qobs-Qcalc) |
|---|
| 181 | elif peak[2]: |
|---|
| 182 | X20 += 1 |
|---|
| 183 | if Nobs20 == 19: |
|---|
| 184 | d20 = peak[7] |
|---|
| 185 | break |
|---|
| 186 | else: |
|---|
| 187 | d20 = peak[7] |
|---|
| 188 | Nobs20 = len(peaks) |
|---|
| 189 | for N20,hkl in enumerate(HKL): |
|---|
| 190 | if hkl[3] < d20: |
|---|
| 191 | break |
|---|
| 192 | Q20 = 1.0/d20**2 |
|---|
| 193 | if diff: |
|---|
| 194 | M20 = Q20/(2.0*diff) |
|---|
| 195 | else: |
|---|
| 196 | M20 = 0 |
|---|
| 197 | if ifX20: |
|---|
| 198 | M20 /= (1.+X20) |
|---|
| 199 | return M20,X20 |
|---|
| 200 | |
|---|
| 201 | def calc_M20SS(peaks,HKL): |
|---|
| 202 | 'needs a doc string' |
|---|
| 203 | diff = 0 |
|---|
| 204 | X20 = 0 |
|---|
| 205 | for Nobs20,peak in enumerate(peaks): |
|---|
| 206 | if peak[3]: |
|---|
| 207 | Qobs = 1.0/peak[8]**2 |
|---|
| 208 | Qcalc = 1.0/peak[9]**2 |
|---|
| 209 | diff += abs(Qobs-Qcalc) |
|---|
| 210 | elif peak[2]: |
|---|
| 211 | X20 += 1 |
|---|
| 212 | if Nobs20 == 19: |
|---|
| 213 | d20 = peak[8] |
|---|
| 214 | break |
|---|
| 215 | else: |
|---|
| 216 | d20 = peak[8] |
|---|
| 217 | Nobs20 = len(peaks) |
|---|
| 218 | for N20,hkl in enumerate(HKL): |
|---|
| 219 | if hkl[4] < d20: |
|---|
| 220 | break |
|---|
| 221 | Q20 = 1.0/d20**2 |
|---|
| 222 | if diff: |
|---|
| 223 | M20 = Q20/(2.0*diff) |
|---|
| 224 | else: |
|---|
| 225 | M20 = 0 |
|---|
| 226 | M20 /= (1.+X20) |
|---|
| 227 | return M20,X20 |
|---|
| 228 | |
|---|
| 229 | def sortM20(cells): |
|---|
| 230 | 'needs a doc string' |
|---|
| 231 | #cells is M20,X20,Bravais,a,b,c,alp,bet,gam |
|---|
| 232 | #sort highest M20 1st |
|---|
| 233 | T = [] |
|---|
| 234 | for i,M in enumerate(cells): |
|---|
| 235 | T.append((M[0],i)) |
|---|
| 236 | D = dict(zip(T,cells)) |
|---|
| 237 | T.sort() |
|---|
| 238 | T.reverse() |
|---|
| 239 | X = [] |
|---|
| 240 | for key in T: |
|---|
| 241 | X.append(D[key]) |
|---|
| 242 | return X |
|---|
| 243 | |
|---|
| 244 | def sortCells(cells,col): |
|---|
| 245 | #cells is M20,X20,Bravais,a,b,c,alp,bet,gam,volume |
|---|
| 246 | #sort smallest a,b,c,alpha,beta,gamma or volume 1st |
|---|
| 247 | T = [] |
|---|
| 248 | for i,M in enumerate(cells): |
|---|
| 249 | T.append((M[col],i)) |
|---|
| 250 | D = dict(zip(T,cells)) |
|---|
| 251 | T.sort() |
|---|
| 252 | X = [] |
|---|
| 253 | for key in T: |
|---|
| 254 | X.append(D[key]) |
|---|
| 255 | return X |
|---|
| 256 | |
|---|
| 257 | def findMV(peaks,controls,ssopt,Inst,dlg): |
|---|
| 258 | |
|---|
| 259 | def Val2Vec(vec,Vref,values): |
|---|
| 260 | Vec = [] |
|---|
| 261 | i = 0 |
|---|
| 262 | for j,r in enumerate(Vref): |
|---|
| 263 | if r: |
|---|
| 264 | if values.size > 1: |
|---|
| 265 | Vec.append(max(-1.,min(1.0,values[i]))) |
|---|
| 266 | else: |
|---|
| 267 | Vec.append(max(0.0,min(1.0,values))) |
|---|
| 268 | i += 1 |
|---|
| 269 | else: |
|---|
| 270 | Vec.append(vec[j]) |
|---|
| 271 | return np.array(Vec) |
|---|
| 272 | |
|---|
| 273 | def ZSSfunc(values,peaks,dmin,Inst,SGData,SSGData,vec,Vref,maxH,A,wave,Z,dlg=None): |
|---|
| 274 | Vec = Val2Vec(vec,Vref,values) |
|---|
| 275 | HKL = G2pwd.getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A) |
|---|
| 276 | Peaks = np.array(IndexSSPeaks(peaks,HKL)[1]).T |
|---|
| 277 | Qo = 1./Peaks[-2]**2 |
|---|
| 278 | Qc = G2lat.calc_rDsqZSS(Peaks[4:8],A,Vec,Z,Peaks[0],wave) |
|---|
| 279 | chi = np.sum((Qo-Qc)**2) |
|---|
| 280 | if dlg: |
|---|
| 281 | dlg.Pulse() |
|---|
| 282 | return chi |
|---|
| 283 | |
|---|
| 284 | def TSSfunc(values,peaks,dmin,Inst,SGData,SSGData,vec,Vref,maxH,A,difC,Z,dlg=None): |
|---|
| 285 | Vec = Val2Vec(vec,Vref,values) |
|---|
| 286 | HKL = G2pwd.getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A) |
|---|
| 287 | Peaks = np.array(IndexSSPeaks(peaks,HKL)[1]).T |
|---|
| 288 | Qo = 1./Peaks[-2]**2 |
|---|
| 289 | Qc = G2lat.calc_rDsqTSS(Peaks[4:8],A,vec,Z,Peaks[0],difC) |
|---|
| 290 | chi = np.sum((Qo-Qc)**2) |
|---|
| 291 | if dlg: |
|---|
| 292 | dlg.Pulse() |
|---|
| 293 | return chi |
|---|
| 294 | |
|---|
| 295 | if 'T' in Inst['Type'][0]: |
|---|
| 296 | difC = Inst['difC'][1] |
|---|
| 297 | else: |
|---|
| 298 | wave = G2mth.getWave(Inst) |
|---|
| 299 | SGData = G2spc.SpcGroup(controls[13])[1] |
|---|
| 300 | SSGData = G2spc.SSpcGroup(SGData,ssopt['ssSymb'])[1] |
|---|
| 301 | A = G2lat.cell2A(controls[6:12]) |
|---|
| 302 | Z = controls[1] |
|---|
| 303 | Vref = [True if x in ssopt['ssSymb'] else False for x in ['a','b','g']] |
|---|
| 304 | values = [] |
|---|
| 305 | ranges = [] |
|---|
| 306 | dT = 0.02 #seems to be a good choice |
|---|
| 307 | for v,r in zip(ssopt['ModVec'],Vref): |
|---|
| 308 | if r: |
|---|
| 309 | ranges += [slice(dT,1.-dT,dT),] #NB: unique part for (00g) & (a0g); (abg)? |
|---|
| 310 | values += [v,] |
|---|
| 311 | dmin = getDmin(peaks)-0.005 |
|---|
| 312 | if 'T' in Inst['Type'][0]: |
|---|
| 313 | result = so.brute(TSSfunc,ranges,finish=so.fmin_cg,full_output=True, |
|---|
| 314 | args=(peaks,dmin,Inst,SGData,SSGData,ssopt['ModVec'],Vref,1,A,difC,Z,dlg)) |
|---|
| 315 | else: |
|---|
| 316 | result = so.brute(ZSSfunc,ranges,finish=so.fmin_cg,full_output=True, |
|---|
| 317 | args=(peaks,dmin,Inst,SGData,SSGData,ssopt['ModVec'],Vref,1,A,wave,Z,dlg)) |
|---|
| 318 | return Val2Vec(ssopt['ModVec'],Vref,result[0]),result |
|---|
| 319 | |
|---|
| 320 | def IndexPeaks(peaks,HKL): |
|---|
| 321 | 'needs a doc string' |
|---|
| 322 | import bisect |
|---|
| 323 | N = len(HKL) |
|---|
| 324 | if N == 0: return False,peaks |
|---|
| 325 | hklds = list(np.array(HKL).T[3])+[1000.0,0.0,] |
|---|
| 326 | hklds.sort() # ascending sort - upper bound at end |
|---|
| 327 | hklmax = [0,0,0] |
|---|
| 328 | for ipk,peak in enumerate(peaks): |
|---|
| 329 | peak[4:7] = [0,0,0] #clear old indexing |
|---|
| 330 | peak[8] = 0. |
|---|
| 331 | if peak[2]: |
|---|
| 332 | i = bisect.bisect_right(hklds,peak[7]) # find peak position in hkl list |
|---|
| 333 | dm = peak[-2]-hklds[i-1] # peak to neighbor hkls in list |
|---|
| 334 | dp = hklds[i]-peak[-2] |
|---|
| 335 | pos = N-i # reverse the order |
|---|
| 336 | if dp > dm: pos += 1 # closer to upper than lower |
|---|
| 337 | if pos >= N: |
|---|
| 338 | break |
|---|
| 339 | hkl = HKL[pos] # put in hkl |
|---|
| 340 | if hkl[-1] >= 0: # peak already assigned - test if this one better |
|---|
| 341 | opeak = peaks[int(hkl[-1])] #hkl[-1] needs to be int here |
|---|
| 342 | dold = abs(opeak[-2]-hkl[3]) |
|---|
| 343 | dnew = min(dm,dp) |
|---|
| 344 | if dold > dnew: # new better - zero out old |
|---|
| 345 | opeak[4:7] = [0,0,0] |
|---|
| 346 | opeak[8] = 0. |
|---|
| 347 | else: # old better - do nothing |
|---|
| 348 | continue |
|---|
| 349 | hkl[-1] = ipk |
|---|
| 350 | peak[4:7] = hkl[:3] |
|---|
| 351 | peak[8] = hkl[3] # fill in d-calc |
|---|
| 352 | for peak in peaks: |
|---|
| 353 | peak[3] = False |
|---|
| 354 | if peak[2]: |
|---|
| 355 | if peak[-1] > 0.: |
|---|
| 356 | for j in range(3): |
|---|
| 357 | if abs(peak[j+4]) > hklmax[j]: hklmax[j] = abs(peak[j+4]) |
|---|
| 358 | peak[3] = True |
|---|
| 359 | if hklmax[0]*hklmax[1]*hklmax[2] > 0: |
|---|
| 360 | return True,peaks |
|---|
| 361 | else: |
|---|
| 362 | return False,peaks #nothing indexed! |
|---|
| 363 | |
|---|
| 364 | def IndexSSPeaks(peaks,HKL): |
|---|
| 365 | 'needs a doc string' |
|---|
| 366 | import bisect |
|---|
| 367 | N = len(HKL) |
|---|
| 368 | Peaks = np.copy(peaks) |
|---|
| 369 | if N == 0: return False,Peaks |
|---|
| 370 | if len(peaks[0]) == 9: #add m column if missing |
|---|
| 371 | Peaks = np.insert(Peaks,7,np.zeros_like(Peaks.T[0]),axis=1) |
|---|
| 372 | # for peak in Peaks: |
|---|
| 373 | # peak.insert(7,0) |
|---|
| 374 | hklds = list(np.array(HKL).T[4])+[1000.0,0.0,] |
|---|
| 375 | hklds.sort() # ascending sort - upper bound at end |
|---|
| 376 | hklmax = [0,0,0,0] |
|---|
| 377 | for ipk,peak in enumerate(Peaks): |
|---|
| 378 | peak[4:8] = [0,0,0,0] #clear old indexing |
|---|
| 379 | peak[9] = 0. |
|---|
| 380 | if peak[2]: #Use |
|---|
| 381 | i = bisect.bisect_right(hklds,peak[8]) # find peak position in hkl list |
|---|
| 382 | dm = peak[8]-hklds[i-1] # peak to neighbor hkls in list |
|---|
| 383 | dp = hklds[i]-peak[8] |
|---|
| 384 | pos = N-i # reverse the order |
|---|
| 385 | if dp > dm: pos += 1 # closer to upper than lower |
|---|
| 386 | if pos >= N: |
|---|
| 387 | print pos,N |
|---|
| 388 | break |
|---|
| 389 | hkl = HKL[pos] # put in hkl |
|---|
| 390 | if hkl[-1] >= 0: # peak already assigned - test if this one better |
|---|
| 391 | opeak = Peaks[hkl[-1]] |
|---|
| 392 | dold = abs(opeak[-2]-hkl[4]) |
|---|
| 393 | dnew = min(dm,dp) |
|---|
| 394 | if dold > dnew: # new better - zero out old |
|---|
| 395 | opeak[4:8] = [0,0,0,0] |
|---|
| 396 | opeak[9] = 0. |
|---|
| 397 | else: # old better - do nothing |
|---|
| 398 | continue |
|---|
| 399 | hkl[-1] = ipk |
|---|
| 400 | peak[4:8] = hkl[:4] |
|---|
| 401 | peak[9] = hkl[4] # fill in d-calc |
|---|
| 402 | for peak in Peaks: |
|---|
| 403 | peak[3] = False |
|---|
| 404 | if peak[2]: |
|---|
| 405 | if peak[-1] > 0.: |
|---|
| 406 | for j in range(4): |
|---|
| 407 | if abs(peak[j+4]) > hklmax[j]: hklmax[j] = abs(peak[j+4]) |
|---|
| 408 | peak[3] = True |
|---|
| 409 | if hklmax[0]*hklmax[1]*hklmax[2]*hklmax[3] > 0: |
|---|
| 410 | return True,Peaks |
|---|
| 411 | else: |
|---|
| 412 | return False,Peaks #nothing indexed! |
|---|
| 413 | |
|---|
| 414 | def Values2A(ibrav,values): |
|---|
| 415 | 'needs a doc string' |
|---|
| 416 | if ibrav in [0,1,2]: |
|---|
| 417 | return [values[0],values[0],values[0],0,0,0] |
|---|
| 418 | elif ibrav in [3,4]: |
|---|
| 419 | return [values[0],values[0],values[1],values[0],0,0] |
|---|
| 420 | elif ibrav in [5,6]: |
|---|
| 421 | return [values[0],values[0],values[1],0,0,0] |
|---|
| 422 | elif ibrav in [7,8,9,10]: |
|---|
| 423 | return [values[0],values[1],values[2],0,0,0] |
|---|
| 424 | elif ibrav in [11,12]: |
|---|
| 425 | return [values[0],values[1],values[2],0,values[3],0] |
|---|
| 426 | else: |
|---|
| 427 | return list(values[:6]) |
|---|
| 428 | |
|---|
| 429 | def A2values(ibrav,A): |
|---|
| 430 | 'needs a doc string' |
|---|
| 431 | if ibrav in [0,1,2]: |
|---|
| 432 | return [A[0],] |
|---|
| 433 | elif ibrav in [3,4,5,6]: |
|---|
| 434 | return [A[0],A[2]] |
|---|
| 435 | elif ibrav in [7,8,9,10]: |
|---|
| 436 | return [A[0],A[1],A[2]] |
|---|
| 437 | elif ibrav in [11,12]: |
|---|
| 438 | return [A[0],A[1],A[2],A[4]] |
|---|
| 439 | else: |
|---|
| 440 | return A |
|---|
| 441 | |
|---|
| 442 | def Values2Vec(ibrav,vec,Vref,val): |
|---|
| 443 | if ibrav in [3,4,5,6]: |
|---|
| 444 | Nskip = 2 |
|---|
| 445 | elif ibrav in [7,8,9,10]: |
|---|
| 446 | Nskip = 3 |
|---|
| 447 | elif ibrav in [11,12]: |
|---|
| 448 | Nskip = 4 |
|---|
| 449 | else: |
|---|
| 450 | Nskip = 6 |
|---|
| 451 | Vec = [] |
|---|
| 452 | i = 0 |
|---|
| 453 | for j,r in enumerate(Vref): |
|---|
| 454 | if r: |
|---|
| 455 | Vec.append(val[i+Nskip]) |
|---|
| 456 | i += 1 |
|---|
| 457 | else: |
|---|
| 458 | Vec.append(vec[j]) |
|---|
| 459 | return np.array(Vec) |
|---|
| 460 | |
|---|
| 461 | def FitHKL(ibrav,peaks,A,Pwr): |
|---|
| 462 | 'needs a doc string' |
|---|
| 463 | |
|---|
| 464 | def errFit(values,ibrav,d,H,Pwr): |
|---|
| 465 | A = Values2A(ibrav,values) |
|---|
| 466 | Qo = 1./d**2 |
|---|
| 467 | Qc = G2lat.calc_rDsq(H,A) |
|---|
| 468 | return (Qo-Qc)*d**Pwr |
|---|
| 469 | |
|---|
| 470 | def dervFit(values,ibrav,d,H,Pwr): |
|---|
| 471 | if ibrav in [0,1,2]: |
|---|
| 472 | derv = [H[0]*H[0]+H[1]*H[1]+H[2]*H[2],] |
|---|
| 473 | elif ibrav in [3,4,]: |
|---|
| 474 | derv = [H[0]*H[0]+H[1]*H[1]+H[0]*H[1],H[2]*H[2]] |
|---|
| 475 | elif ibrav in [5,6]: |
|---|
| 476 | derv = [H[0]*H[0]+H[1]*H[1],H[2]*H[2]] |
|---|
| 477 | elif ibrav in [7,8,9,10]: |
|---|
| 478 | derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2]] |
|---|
| 479 | elif ibrav in [11,12]: |
|---|
| 480 | derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[2]] |
|---|
| 481 | else: |
|---|
| 482 | derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[1],H[0]*H[2],H[1]*H[2]] |
|---|
| 483 | derv = -np.array(derv) |
|---|
| 484 | return (derv*d**Pwr).T |
|---|
| 485 | |
|---|
| 486 | Peaks = np.array(peaks).T |
|---|
| 487 | values = A2values(ibrav,A) |
|---|
| 488 | result = so.leastsq(errFit,values,Dfun=dervFit,full_output=True,ftol=0.000001, |
|---|
| 489 | args=(ibrav,Peaks[7],Peaks[4:7],Pwr)) |
|---|
| 490 | A = Values2A(ibrav,result[0]) |
|---|
| 491 | return True,np.sum(errFit(result[0],ibrav,Peaks[7],Peaks[4:7],Pwr)**2),A,result |
|---|
| 492 | |
|---|
| 493 | def errFitZ(values,ibrav,d,H,tth,wave,Z,Zref): |
|---|
| 494 | Zero = Z |
|---|
| 495 | if Zref: |
|---|
| 496 | Zero = values[-1] |
|---|
| 497 | A = Values2A(ibrav,values) |
|---|
| 498 | Qo = 1./d**2 |
|---|
| 499 | Qc = G2lat.calc_rDsqZ(H,A,Zero,tth,wave) |
|---|
| 500 | return (Qo-Qc) |
|---|
| 501 | |
|---|
| 502 | def dervFitZ(values,ibrav,d,H,tth,wave,Z,Zref): |
|---|
| 503 | if ibrav in [0,1,2]: |
|---|
| 504 | derv = [H[0]*H[0]+H[1]*H[1]+H[2]*H[2],] |
|---|
| 505 | elif ibrav in [3,4,]: |
|---|
| 506 | derv = [H[0]*H[0]+H[1]*H[1]+H[0]*H[1],H[2]*H[2]] |
|---|
| 507 | elif ibrav in [5,6]: |
|---|
| 508 | derv = [H[0]*H[0]+H[1]*H[1],H[2]*H[2]] |
|---|
| 509 | elif ibrav in [7,8,9,10]: |
|---|
| 510 | derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2]] |
|---|
| 511 | elif ibrav in [11,12]: |
|---|
| 512 | derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[2]] |
|---|
| 513 | else: |
|---|
| 514 | derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[1],H[0]*H[2],H[1]*H[2]] |
|---|
| 515 | if Zref: |
|---|
| 516 | derv.append(npsind(tth)*2.0*rpd/wave**2) |
|---|
| 517 | derv = -np.array(derv) |
|---|
| 518 | return derv.T |
|---|
| 519 | |
|---|
| 520 | def FitHKLZ(wave,ibrav,peaks,A,Z,Zref): |
|---|
| 521 | 'needs a doc string' |
|---|
| 522 | |
|---|
| 523 | Peaks = np.array(peaks).T |
|---|
| 524 | values = A2values(ibrav,A) |
|---|
| 525 | if Zref: |
|---|
| 526 | values.append(Z) |
|---|
| 527 | result = so.leastsq(errFitZ,values,Dfun=dervFitZ,full_output=True,ftol=0.0001, |
|---|
| 528 | args=(ibrav,Peaks[7],Peaks[4:7],Peaks[0],wave,Z,Zref)) |
|---|
| 529 | A = Values2A(ibrav,result[0][:6]) |
|---|
| 530 | if Zref: |
|---|
| 531 | Z = result[0][-1] |
|---|
| 532 | chisq = np.sum(errFitZ(result[0],ibrav,Peaks[7],Peaks[4:7],Peaks[0],wave,Z,Zref)**2) |
|---|
| 533 | return True,chisq,A,Z,result |
|---|
| 534 | |
|---|
| 535 | def errFitZSS(values,ibrav,d,H,tth,wave,vec,Vref,Z,Zref): |
|---|
| 536 | Zero = Z |
|---|
| 537 | if Zref: |
|---|
| 538 | Zero = values[-1] |
|---|
| 539 | A = Values2A(ibrav,values) |
|---|
| 540 | Vec = Values2Vec(ibrav,vec,Vref,values) |
|---|
| 541 | Qo = 1./d**2 |
|---|
| 542 | Qc = G2lat.calc_rDsqZSS(H,A,Vec,Zero,tth,wave) |
|---|
| 543 | return (Qo-Qc) |
|---|
| 544 | |
|---|
| 545 | def dervFitZSS(values,ibrav,d,H,tth,wave,vec,Vref,Z,Zref): |
|---|
| 546 | A = Values2A(ibrav,values) |
|---|
| 547 | Vec = Values2Vec(ibrav,vec,Vref,values) |
|---|
| 548 | HM = H[:3]+(H[3][:,np.newaxis]*Vec).T |
|---|
| 549 | if ibrav in [3,4,]: |
|---|
| 550 | derv = [HM[0]*HM[0]+HM[1]*HM[1]+HM[0]*HM[1],HM[2]*HM[2]] |
|---|
| 551 | elif ibrav in [5,6]: |
|---|
| 552 | derv = [HM[0]*HM[0]+HM[1]*HM[1],HM[2]*HM[2]] |
|---|
| 553 | elif ibrav in [7,8,9,10]: |
|---|
| 554 | derv = [HM[0]*HM[0],HM[1]*HM[1],HM[2]*HM[2]] |
|---|
| 555 | elif ibrav in [11,12]: |
|---|
| 556 | derv = [HM[0]*HM[0],HM[1]*HM[1],HM[2]*HM[2],HM[0]*HM[2]] |
|---|
| 557 | else: |
|---|
| 558 | derv = [HM[0]*HM[0],HM[1]*HM[1],HM[2]*HM[2],HM[0]*HM[1],HM[0]*HM[2],HM[1]*HM[2]] |
|---|
| 559 | if Vref[0]: |
|---|
| 560 | derv.append(2.*A[0]*HM[0]*H[3]+A[3]*HM[1]*H[3]+A[4]*HM[2]*H[3]) |
|---|
| 561 | if Vref[1]: |
|---|
| 562 | derv.append(2.*A[1]*HM[1]*H[3]+A[3]*HM[0]*H[3]+A[5]*HM[2]*H[3]) |
|---|
| 563 | if Vref[2]: |
|---|
| 564 | derv.append(2.*A[2]*HM[2]*H[3]+A[4]*HM[1]*H[3]+A[5]*HM[0]*H[3]) |
|---|
| 565 | if Zref: |
|---|
| 566 | derv.append(npsind(tth)*2.0*rpd/wave**2) |
|---|
| 567 | derv = -np.array(derv) |
|---|
| 568 | return derv.T |
|---|
| 569 | |
|---|
| 570 | def FitHKLZSS(wave,ibrav,peaks,A,V,Vref,Z,Zref): |
|---|
| 571 | 'needs a doc string' |
|---|
| 572 | |
|---|
| 573 | Peaks = np.array(peaks).T |
|---|
| 574 | values = A2values(ibrav,A) |
|---|
| 575 | for v,r in zip(V,Vref): |
|---|
| 576 | if r: |
|---|
| 577 | values.append(v) |
|---|
| 578 | if Zref: |
|---|
| 579 | values.append(Z) |
|---|
| 580 | result = so.leastsq(errFitZSS,values,Dfun=dervFitZSS,full_output=True,ftol=1.e-6, |
|---|
| 581 | args=(ibrav,Peaks[8],Peaks[4:8],Peaks[0],wave,V,Vref,Z,Zref)) |
|---|
| 582 | A = Values2A(ibrav,result[0]) |
|---|
| 583 | Vec = Values2Vec(ibrav,V,Vref,result[0]) |
|---|
| 584 | if Zref: |
|---|
| 585 | Z = result[0][-1] |
|---|
| 586 | chisq = np.sum(errFitZSS(result[0],ibrav,Peaks[8],Peaks[4:8],Peaks[0],wave,Vec,Vref,Z,Zref)**2) |
|---|
| 587 | return True,chisq,A,Vec,Z,result |
|---|
| 588 | |
|---|
| 589 | def errFitT(values,ibrav,d,H,tof,difC,Z,Zref): |
|---|
| 590 | Zero = Z |
|---|
| 591 | if Zref: |
|---|
| 592 | Zero = values[-1] |
|---|
| 593 | A = Values2A(ibrav,values) |
|---|
| 594 | Qo = 1./d**2 |
|---|
| 595 | Qc = G2lat.calc_rDsqT(H,A,Zero,tof,difC) |
|---|
| 596 | return (Qo-Qc) |
|---|
| 597 | |
|---|
| 598 | def dervFitT(values,ibrav,d,H,tof,difC,Z,Zref): |
|---|
| 599 | if ibrav in [0,1,2]: |
|---|
| 600 | derv = [H[0]*H[0]+H[1]*H[1]+H[2]*H[2],] |
|---|
| 601 | elif ibrav in [3,4,]: |
|---|
| 602 | derv = [H[0]*H[0]+H[1]*H[1]+H[0]*H[1],H[2]*H[2]] |
|---|
| 603 | elif ibrav in [5,6]: |
|---|
| 604 | derv = [H[0]*H[0]+H[1]*H[1],H[2]*H[2]] |
|---|
| 605 | elif ibrav in [7,8,9,10]: |
|---|
| 606 | derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2]] |
|---|
| 607 | elif ibrav in [11,12]: |
|---|
| 608 | derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[2]] |
|---|
| 609 | else: |
|---|
| 610 | derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[1],H[0]*H[2],H[1]*H[2]] |
|---|
| 611 | if Zref: |
|---|
| 612 | derv.append(np.ones_like(d)/difC) |
|---|
| 613 | derv = np.array(derv) |
|---|
| 614 | return derv.T |
|---|
| 615 | |
|---|
| 616 | def FitHKLT(difC,ibrav,peaks,A,Z,Zref): |
|---|
| 617 | 'needs a doc string' |
|---|
| 618 | |
|---|
| 619 | Peaks = np.array(peaks).T |
|---|
| 620 | values = A2values(ibrav,A) |
|---|
| 621 | if Zref: |
|---|
| 622 | values.append(Z) |
|---|
| 623 | result = so.leastsq(errFitT,values,Dfun=dervFitT,full_output=True,ftol=0.0001, |
|---|
| 624 | args=(ibrav,Peaks[7],Peaks[4:7],Peaks[0],difC,Z,Zref)) |
|---|
| 625 | A = Values2A(ibrav,result[0]) |
|---|
| 626 | if Zref: |
|---|
| 627 | Z = result[0][-1] |
|---|
| 628 | chisq = np.sum(errFitT(result[0],ibrav,Peaks[7],Peaks[4:7],Peaks[0],difC,Z,Zref)**2) |
|---|
| 629 | return True,chisq,A,Z,result |
|---|
| 630 | |
|---|
| 631 | def errFitTSS(values,ibrav,d,H,tof,difC,vec,Vref,Z,Zref): |
|---|
| 632 | Zero = Z |
|---|
| 633 | if Zref: |
|---|
| 634 | Zero = values[-1] |
|---|
| 635 | A = Values2A(ibrav,values) |
|---|
| 636 | Vec = Values2Vec(ibrav,vec,Vref,values) |
|---|
| 637 | Qo = 1./d**2 |
|---|
| 638 | Qc = G2lat.calc_rDsqTSS(H,A,Vec,Zero,tof,difC) |
|---|
| 639 | return (Qo-Qc) |
|---|
| 640 | |
|---|
| 641 | def dervFitTSS(values,ibrav,d,H,tof,difC,vec,Vref,Z,Zref): |
|---|
| 642 | A = Values2A(ibrav,values) |
|---|
| 643 | Vec = Values2Vec(ibrav,vec,Vref,values) |
|---|
| 644 | HM = H[:3]+(H[3][:,np.newaxis]*Vec).T |
|---|
| 645 | if ibrav in [3,4,]: |
|---|
| 646 | derv = [HM[0]*HM[0]+HM[1]*HM[1]+HM[0]*HM[1],HM[2]*HM[2]] |
|---|
| 647 | elif ibrav in [5,6]: |
|---|
| 648 | derv = [HM[0]*HM[0]+HM[1]*HM[1],HM[2]*HM[2]] |
|---|
| 649 | elif ibrav in [7,8,9,10]: |
|---|
| 650 | derv = [HM[0]*HM[0],HM[1]*HM[1],HM[2]*HM[2]] |
|---|
| 651 | elif ibrav in [11,12]: |
|---|
| 652 | derv = [HM[0]*HM[0],HM[1]*HM[1],HM[2]*HM[2],HM[0]*HM[2]] |
|---|
| 653 | else: |
|---|
| 654 | derv = [HM[0]*HM[0],HM[1]*HM[1],HM[2]*HM[2],HM[0]*HM[1],HM[0]*HM[2],HM[1]*HM[2]] |
|---|
| 655 | if Vref[0]: |
|---|
| 656 | derv.append(2.*A[0]*HM[0]*H[3]+A[3]*HM[1]*H[3]+A[4]*HM[2]*H[3]) |
|---|
| 657 | if Vref[1]: |
|---|
| 658 | derv.append(2.*A[1]*HM[1]*H[3]+A[3]*HM[0]*H[3]+A[5]*HM[2]*H[3]) |
|---|
| 659 | if Vref[2]: |
|---|
| 660 | derv.append(2.*A[2]*HM[2]*H[3]+A[4]*HM[1]*H[3]+A[5]*HM[0]*H[3]) |
|---|
| 661 | if Zref: |
|---|
| 662 | derv.append(np.ones_like(d)/difC) |
|---|
| 663 | derv = np.array(derv) |
|---|
| 664 | return derv.T |
|---|
| 665 | |
|---|
| 666 | def FitHKLTSS(difC,ibrav,peaks,A,V,Vref,Z,Zref): |
|---|
| 667 | 'needs a doc string' |
|---|
| 668 | |
|---|
| 669 | Peaks = np.array(peaks).T |
|---|
| 670 | values = A2values(ibrav,A) |
|---|
| 671 | for v,r in zip(V,Vref): |
|---|
| 672 | if r: |
|---|
| 673 | values.append(v) |
|---|
| 674 | if Zref: |
|---|
| 675 | values.append(Z) |
|---|
| 676 | result = so.leastsq(errFitTSS,values,Dfun=dervFitTSS,full_output=True,ftol=0.0001, |
|---|
| 677 | args=(ibrav,Peaks[8],Peaks[4:8],Peaks[0],difC,V,Vref,Z,Zref)) |
|---|
| 678 | A = Values2A(ibrav,result[0]) |
|---|
| 679 | Vec = Values2Vec(ibrav,V,Vref,result[0]) |
|---|
| 680 | if Zref: |
|---|
| 681 | Z = result[0][-1] |
|---|
| 682 | chisq = np.sum(errFitTSS(result[0],ibrav,Peaks[8],Peaks[4:8],Peaks[0],difC,V,Vref,Z,Zref)**2) |
|---|
| 683 | return True,chisq,A,Vec,Z,result |
|---|
| 684 | |
|---|
| 685 | def rotOrthoA(A): |
|---|
| 686 | 'needs a doc string' |
|---|
| 687 | return [A[1],A[2],A[0],0,0,0] |
|---|
| 688 | |
|---|
| 689 | def swapMonoA(A): |
|---|
| 690 | 'needs a doc string' |
|---|
| 691 | return [A[2],A[1],A[0],0,A[4],0] |
|---|
| 692 | |
|---|
| 693 | def oddPeak(indx,peaks): |
|---|
| 694 | 'needs a doc string' |
|---|
| 695 | noOdd = True |
|---|
| 696 | for peak in peaks: |
|---|
| 697 | H = peak[4:7] |
|---|
| 698 | if H[indx] % 2: |
|---|
| 699 | noOdd = False |
|---|
| 700 | return noOdd |
|---|
| 701 | |
|---|
| 702 | def halfCell(ibrav,A,peaks): |
|---|
| 703 | 'needs a doc string' |
|---|
| 704 | if ibrav in [0,1,2]: |
|---|
| 705 | if oddPeak(0,peaks): |
|---|
| 706 | A[0] *= 2 |
|---|
| 707 | A[1] = A[2] = A[0] |
|---|
| 708 | elif ibrav in [3,4,5,6]: |
|---|
| 709 | if oddPeak(0,peaks): |
|---|
| 710 | A[0] *= 2 |
|---|
| 711 | A[1] = A[0] |
|---|
| 712 | if oddPeak(2,peaks): |
|---|
| 713 | A[2] *=2 |
|---|
| 714 | else: |
|---|
| 715 | if oddPeak(0,peaks): |
|---|
| 716 | A[0] *=2 |
|---|
| 717 | if oddPeak(1,peaks): |
|---|
| 718 | A[1] *=2 |
|---|
| 719 | if oddPeak(2,peaks): |
|---|
| 720 | A[2] *=2 |
|---|
| 721 | return A |
|---|
| 722 | |
|---|
| 723 | def getDmin(peaks): |
|---|
| 724 | 'needs a doc string' |
|---|
| 725 | return peaks[-1][-2] |
|---|
| 726 | |
|---|
| 727 | def getDmax(peaks): |
|---|
| 728 | 'needs a doc string' |
|---|
| 729 | return peaks[0][-2] |
|---|
| 730 | |
|---|
| 731 | def refinePeaksZ(peaks,wave,ibrav,A,Zero,ZeroRef): |
|---|
| 732 | 'needs a doc string' |
|---|
| 733 | dmin = getDmin(peaks) |
|---|
| 734 | OK,smin,Aref,Z,result = FitHKLZ(wave,ibrav,peaks,A,Zero,ZeroRef) |
|---|
| 735 | Peaks = np.array(peaks).T |
|---|
| 736 | H = Peaks[4:7] |
|---|
| 737 | Peaks[8] = 1./np.sqrt(G2lat.calc_rDsqZ(H,Aref,Z,Peaks[0],wave)) |
|---|
| 738 | peaks = Peaks.T |
|---|
| 739 | HKL = G2lat.GenHBravais(dmin,ibrav,A) |
|---|
| 740 | M20,X20 = calc_M20(peaks,HKL) |
|---|
| 741 | return len(HKL),M20,X20,Aref,Z |
|---|
| 742 | |
|---|
| 743 | def refinePeaksT(peaks,difC,ibrav,A,Zero,ZeroRef): |
|---|
| 744 | 'needs a doc string' |
|---|
| 745 | dmin = getDmin(peaks) |
|---|
| 746 | OK,smin,Aref,Z,result = FitHKLT(difC,ibrav,peaks,A,Zero,ZeroRef) |
|---|
| 747 | Peaks = np.array(peaks).T |
|---|
| 748 | H = Peaks[4:7] |
|---|
| 749 | Peaks[8] = 1./np.sqrt(G2lat.calc_rDsqT(H,Aref,Z,Peaks[0],difC)) |
|---|
| 750 | peaks = Peaks.T |
|---|
| 751 | HKL = G2lat.GenHBravais(dmin,ibrav,A) |
|---|
| 752 | M20,X20 = calc_M20(peaks,HKL) |
|---|
| 753 | return len(HKL),M20,X20,Aref,Z |
|---|
| 754 | |
|---|
| 755 | def refinePeaksZSS(peaks,wave,Inst,SGData,SSGData,maxH,ibrav,A,vec,vecRef,Zero,ZeroRef): |
|---|
| 756 | 'needs a doc string' |
|---|
| 757 | dmin = getDmin(peaks) |
|---|
| 758 | OK,smin,Aref,Vref,Z,result = FitHKLZSS(wave,ibrav,peaks,A,vec,vecRef,Zero,ZeroRef) |
|---|
| 759 | Peaks = np.array(peaks).T |
|---|
| 760 | H = Peaks[4:8] |
|---|
| 761 | Peaks[9] = 1./np.sqrt(G2lat.calc_rDsqZSS(H,Aref,Vref,Z,Peaks[0],wave)) #H,A,vec,Z,tth,lam |
|---|
| 762 | peaks = Peaks.T |
|---|
| 763 | HKL = G2pwd.getHKLMpeak(dmin,Inst,SGData,SSGData,Vref,maxH,Aref) |
|---|
| 764 | M20,X20 = calc_M20SS(peaks,HKL) |
|---|
| 765 | return len(HKL),M20,X20,Aref,Vref,Z |
|---|
| 766 | |
|---|
| 767 | def refinePeaksTSS(peaks,difC,Inst,SGData,SSGData,maxH,ibrav,A,vec,vecRef,Zero,ZeroRef): |
|---|
| 768 | 'needs a doc string' |
|---|
| 769 | dmin = getDmin(peaks) |
|---|
| 770 | OK,smin,Aref,Vref,Z,result = FitHKLTSS(difC,ibrav,peaks,A,vec,vecRef,Zero,ZeroRef) |
|---|
| 771 | Peaks = np.array(peaks).T |
|---|
| 772 | H = Peaks[4:8] |
|---|
| 773 | Peaks[9] = 1./np.sqrt(G2lat.calc_rDsqTSS(H,Aref,Vref,Z,Peaks[0],difC)) |
|---|
| 774 | peaks = Peaks.T |
|---|
| 775 | HKL = G2pwd.getHKLMpeak(dmin,Inst,SGData,SSGData,Vref,maxH,Aref) |
|---|
| 776 | HKL = G2lat.GenHBravais(dmin,ibrav,A) |
|---|
| 777 | M20,X20 = calc_M20SS(peaks,HKL) |
|---|
| 778 | return len(HKL),M20,X20,Aref,Vref,Z |
|---|
| 779 | |
|---|
| 780 | def refinePeaks(peaks,ibrav,A,ifX20=True): |
|---|
| 781 | 'needs a doc string' |
|---|
| 782 | dmin = getDmin(peaks) |
|---|
| 783 | smin = 1.0e10 |
|---|
| 784 | pwr = 8 |
|---|
| 785 | maxTries = 10 |
|---|
| 786 | OK = False |
|---|
| 787 | tries = 0 |
|---|
| 788 | HKL = G2lat.GenHBravais(dmin,ibrav,A) |
|---|
| 789 | while len(HKL) > 2 and IndexPeaks(peaks,HKL)[0]: |
|---|
| 790 | Pwr = pwr - (tries % 2) |
|---|
| 791 | HKL = [] |
|---|
| 792 | tries += 1 |
|---|
| 793 | osmin = smin |
|---|
| 794 | oldA = A[:] |
|---|
| 795 | Vold = G2lat.calc_V(oldA) |
|---|
| 796 | OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr) |
|---|
| 797 | Vnew = G2lat.calc_V(A) |
|---|
| 798 | if Vnew > 2.0*Vold or Vnew < 2.: |
|---|
| 799 | A = ranAbyR(ibrav,oldA,tries+1,maxTries,ran2axis) |
|---|
| 800 | OK = False |
|---|
| 801 | continue |
|---|
| 802 | try: |
|---|
| 803 | HKL = G2lat.GenHBravais(dmin,ibrav,A) |
|---|
| 804 | except FloatingPointError: |
|---|
| 805 | A = oldA |
|---|
| 806 | OK = False |
|---|
| 807 | break |
|---|
| 808 | if len(HKL) == 0: break #absurd cell obtained! |
|---|
| 809 | rat = (osmin-smin)/smin |
|---|
| 810 | if abs(rat) < 1.0e-5 or not OK: break |
|---|
| 811 | if tries > maxTries: break |
|---|
| 812 | if OK: |
|---|
| 813 | OK,smin,A,result = FitHKL(ibrav,peaks,A,2) |
|---|
| 814 | Peaks = np.array(peaks).T |
|---|
| 815 | H = Peaks[4:7] |
|---|
| 816 | try: |
|---|
| 817 | Peaks[8] = 1./np.sqrt(G2lat.calc_rDsq(H,A)) |
|---|
| 818 | peaks = Peaks.T |
|---|
| 819 | except FloatingPointError: |
|---|
| 820 | A = oldA |
|---|
| 821 | |
|---|
| 822 | M20,X20 = calc_M20(peaks,HKL,ifX20) |
|---|
| 823 | return len(HKL),M20,X20,A |
|---|
| 824 | |
|---|
| 825 | def findBestCell(dlg,ncMax,A,Ntries,ibrav,peaks,V1,ifX20=True): |
|---|
| 826 | 'needs a doc string' |
|---|
| 827 | # dlg & ncMax are used for wx progress bar |
|---|
| 828 | # A != 0 find the best A near input A, |
|---|
| 829 | # A = 0 for random cell, volume normalized to V1; |
|---|
| 830 | # returns number of generated hkls, M20, X20 & A for best found |
|---|
| 831 | mHKL = [3,3,3, 5,5, 5,5, 7,7,7,7, 9,9, 10] |
|---|
| 832 | dmin = getDmin(peaks)-0.05 |
|---|
| 833 | amin = 2.5 |
|---|
| 834 | amax = 5.*getDmax(peaks) |
|---|
| 835 | Asave = [] |
|---|
| 836 | GoOn = True |
|---|
| 837 | Skip = False |
|---|
| 838 | if A: |
|---|
| 839 | HKL = G2lat.GenHBravais(dmin,ibrav,A[:]) |
|---|
| 840 | if len(HKL) > mHKL[ibrav]: |
|---|
| 841 | peaks = IndexPeaks(peaks,HKL)[1] |
|---|
| 842 | Asave.append([calc_M20(peaks,HKL,ifX20),A[:]]) |
|---|
| 843 | tries = 0 |
|---|
| 844 | while tries < Ntries and GoOn: |
|---|
| 845 | if A: |
|---|
| 846 | Abeg = ranAbyR(ibrav,A,tries+1,Ntries,ran2axis) |
|---|
| 847 | if ibrav in [11,12,13]: #monoclinic & triclinic |
|---|
| 848 | Abeg = ranAbyR(ibrav,A,tries/10+1,Ntries,ran2axis) |
|---|
| 849 | else: |
|---|
| 850 | Abeg = ranAbyV(ibrav,amin,amax,V1) |
|---|
| 851 | HKL = G2lat.GenHBravais(dmin,ibrav,Abeg) |
|---|
| 852 | Nc = len(HKL) |
|---|
| 853 | if Nc >= ncMax: |
|---|
| 854 | GoOn = False |
|---|
| 855 | else: |
|---|
| 856 | if dlg: |
|---|
| 857 | # GoOn,Skip = dlg.Update(100*Nc/ncMax) #wx error doesn't work in 32 bit versions! |
|---|
| 858 | GoOn = dlg.Update(100*Nc/ncMax)[0] |
|---|
| 859 | if Skip or not GoOn: |
|---|
| 860 | GoOn = False |
|---|
| 861 | break |
|---|
| 862 | |
|---|
| 863 | if IndexPeaks(peaks,HKL)[0] and len(HKL) > mHKL[ibrav]: |
|---|
| 864 | Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg,ifX20) |
|---|
| 865 | Asave.append([calc_M20(peaks,HKL,ifX20),Aref[:]]) |
|---|
| 866 | if ibrav == 9: #C-centered orthorhombic |
|---|
| 867 | for i in range(2): |
|---|
| 868 | Abeg = rotOrthoA(Abeg[:]) |
|---|
| 869 | Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg,ifX20) |
|---|
| 870 | HKL = G2lat.GenHBravais(dmin,ibrav,Aref) |
|---|
| 871 | peaks = IndexPeaks(peaks,HKL)[1] |
|---|
| 872 | Asave.append([calc_M20(peaks,HKL,ifX20),Aref[:]]) |
|---|
| 873 | elif ibrav == 11: #C-centered monoclinic |
|---|
| 874 | Abeg = swapMonoA(Abeg[:]) |
|---|
| 875 | Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg,ifX20) |
|---|
| 876 | HKL = G2lat.GenHBravais(dmin,ibrav,Aref) |
|---|
| 877 | peaks = IndexPeaks(peaks,HKL)[1] |
|---|
| 878 | Asave.append([calc_M20(peaks,HKL,ifX20),Aref[:]]) |
|---|
| 879 | else: |
|---|
| 880 | break |
|---|
| 881 | Nc = len(HKL) |
|---|
| 882 | tries += 1 |
|---|
| 883 | X = sortM20(Asave) |
|---|
| 884 | if X: |
|---|
| 885 | Lhkl,M20,X20,A = refinePeaks(peaks,ibrav,X[0][1],ifX20) |
|---|
| 886 | return GoOn,Skip,Lhkl,M20,X20,A |
|---|
| 887 | else: |
|---|
| 888 | return GoOn,Skip,0,0,0,0 |
|---|
| 889 | |
|---|
| 890 | def monoCellReduce(ibrav,A): |
|---|
| 891 | 'needs a doc string' |
|---|
| 892 | a,b,c,alp,bet,gam = G2lat.A2cell(A) |
|---|
| 893 | G,g = G2lat.A2Gmat(A) |
|---|
| 894 | if ibrav in [11]: |
|---|
| 895 | u = [0,0,-1] |
|---|
| 896 | v = [1,0,2] |
|---|
| 897 | anew = math.sqrt(np.dot(np.dot(v,g),v)) |
|---|
| 898 | if anew < a: |
|---|
| 899 | cang = np.dot(np.dot(u,g),v)/(anew*c) |
|---|
| 900 | beta = acosd(-abs(cang)) |
|---|
| 901 | A = G2lat.cell2A([anew,b,c,90,beta,90]) |
|---|
| 902 | else: |
|---|
| 903 | u = [-1,0,0] |
|---|
| 904 | v = [1,0,1] |
|---|
| 905 | cnew = math.sqrt(np.dot(np.dot(v,g),v)) |
|---|
| 906 | if cnew < c: |
|---|
| 907 | cang = np.dot(np.dot(u,g),v)/(a*cnew) |
|---|
| 908 | beta = acosd(-abs(cang)) |
|---|
| 909 | A = G2lat.cell2A([a,b,cnew,90,beta,90]) |
|---|
| 910 | return A |
|---|
| 911 | |
|---|
| 912 | def DoIndexPeaks(peaks,controls,bravais,dlg,ifX20=True): |
|---|
| 913 | 'needs a doc string' |
|---|
| 914 | |
|---|
| 915 | delt = 0.005 #lowest d-spacing cushion - can be fixed? |
|---|
| 916 | amin = 2.5 |
|---|
| 917 | amax = 5.0*getDmax(peaks) |
|---|
| 918 | dmin = getDmin(peaks)-delt |
|---|
| 919 | bravaisNames = ['Cubic-F','Cubic-I','Cubic-P','Trigonal-R','Trigonal/Hexagonal-P', |
|---|
| 920 | 'Tetragonal-I','Tetragonal-P','Orthorhombic-F','Orthorhombic-I','Orthorhombic-C', |
|---|
| 921 | 'Orthorhombic-P','Monoclinic-C','Monoclinic-P','Triclinic'] |
|---|
| 922 | tries = ['1st','2nd','3rd','4th','5th','6th','7th','8th','9th','10th'] |
|---|
| 923 | N1s = [1,1,1, 5,5, 5,5, 50,50,50,50, 50,50, 200] |
|---|
| 924 | N2s = [1,1,1, 2,2, 2,2, 2,2,2,2, 2,2, 4] |
|---|
| 925 | Nm = [1,1,1, 1,1, 1,1, 1,1,1,1, 2,2, 4] |
|---|
| 926 | notUse = 0 |
|---|
| 927 | for peak in peaks: |
|---|
| 928 | if not peak[2]: |
|---|
| 929 | notUse += 1 |
|---|
| 930 | Nobs = len(peaks)-notUse |
|---|
| 931 | zero,ncno = controls[1:3] |
|---|
| 932 | ncMax = Nobs*ncno |
|---|
| 933 | print "%s %8.3f %8.3f" % ('lattice parameter range = ',amin,amax) |
|---|
| 934 | print "%s %.4f %s %d %s %d" % ('Zero =',zero,'Nc/No max =',ncno,' Max Nc =',ncno*Nobs) |
|---|
| 935 | cells = [] |
|---|
| 936 | lastcell = np.zeros(7) |
|---|
| 937 | for ibrav in range(14): |
|---|
| 938 | begin = time.time() |
|---|
| 939 | if bravais[ibrav]: |
|---|
| 940 | print 'cell search for ',bravaisNames[ibrav] |
|---|
| 941 | print ' M20 X20 Nc a b c alpha beta gamma volume V-test' |
|---|
| 942 | V1 = controls[3] |
|---|
| 943 | bestM20 = 0 |
|---|
| 944 | topM20 = 0 |
|---|
| 945 | cycle = 0 |
|---|
| 946 | while cycle < 5: |
|---|
| 947 | if dlg: |
|---|
| 948 | dlg.Update(0,newmsg=tries[cycle]+" cell search for "+bravaisNames[ibrav]) |
|---|
| 949 | try: |
|---|
| 950 | GoOn = True |
|---|
| 951 | while GoOn: #Loop over increment of volume |
|---|
| 952 | N2 = 0 |
|---|
| 953 | while N2 < N2s[ibrav]: #Table 2 step (iii) |
|---|
| 954 | if ibrav > 2: |
|---|
| 955 | if not N2: |
|---|
| 956 | A = [] |
|---|
| 957 | GoOn,Skip,Nc,M20,X20,A = findBestCell(dlg,ncMax,A,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1,ifX20) |
|---|
| 958 | if Skip: |
|---|
| 959 | break |
|---|
| 960 | if A: |
|---|
| 961 | GoOn,Skip,Nc,M20,X20,A = findBestCell(dlg,ncMax,A[:],N1s[ibrav],ibrav,peaks,0,ifX20) |
|---|
| 962 | else: |
|---|
| 963 | GoOn,Skip,Nc,M20,X20,A = findBestCell(dlg,ncMax,0,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1,ifX20) |
|---|
| 964 | if Skip: |
|---|
| 965 | break |
|---|
| 966 | elif Nc >= ncMax: |
|---|
| 967 | GoOn = False |
|---|
| 968 | break |
|---|
| 969 | elif 3*Nc < Nobs: |
|---|
| 970 | N2 = 10 |
|---|
| 971 | break |
|---|
| 972 | else: |
|---|
| 973 | if not GoOn: |
|---|
| 974 | break |
|---|
| 975 | if 1.e6 > M20 > 1.0: #exclude nonsense |
|---|
| 976 | bestM20 = max(bestM20,M20) |
|---|
| 977 | A = halfCell(ibrav,A[:],peaks) |
|---|
| 978 | if ibrav in [12]: |
|---|
| 979 | A = monoCellReduce(ibrav,A[:]) |
|---|
| 980 | HKL = G2lat.GenHBravais(dmin,ibrav,A) |
|---|
| 981 | peaks = IndexPeaks(peaks,HKL)[1] |
|---|
| 982 | a,b,c,alp,bet,gam = G2lat.A2cell(A) |
|---|
| 983 | V = G2lat.calc_V(A) |
|---|
| 984 | if M20 >= 2.0: |
|---|
| 985 | cell = [M20,X20,ibrav,a,b,c,alp,bet,gam,V,False,False] |
|---|
| 986 | newcell = np.array(cell[3:10]) |
|---|
| 987 | if not np.allclose(newcell,lastcell): |
|---|
| 988 | print "%10.3f %3d %3d %10.5f %10.5f %10.5f %10.3f %10.3f %10.3f %10.2f %10.2f" \ |
|---|
| 989 | %(M20,X20,Nc,a,b,c,alp,bet,gam,V,V1) |
|---|
| 990 | cells.append(cell) |
|---|
| 991 | lastcell = np.array(cell[3:10]) |
|---|
| 992 | if not GoOn: |
|---|
| 993 | break |
|---|
| 994 | N2 += 1 |
|---|
| 995 | if Skip: |
|---|
| 996 | cycle = 10 |
|---|
| 997 | GoOn = False |
|---|
| 998 | break |
|---|
| 999 | if ibrav < 11: |
|---|
| 1000 | V1 *= 1.1 |
|---|
| 1001 | elif ibrav in range(11,14): |
|---|
| 1002 | V1 *= 1.05 |
|---|
| 1003 | if not GoOn: |
|---|
| 1004 | if bestM20 > topM20: |
|---|
| 1005 | topM20 = bestM20 |
|---|
| 1006 | if cells: |
|---|
| 1007 | V1 = cells[0][9] |
|---|
| 1008 | else: |
|---|
| 1009 | V1 = 25 |
|---|
| 1010 | ncMax += Nobs |
|---|
| 1011 | cycle += 1 |
|---|
| 1012 | print 'Restart search, new Max Nc = ',ncMax |
|---|
| 1013 | else: |
|---|
| 1014 | cycle = 10 |
|---|
| 1015 | finally: |
|---|
| 1016 | pass |
|---|
| 1017 | # dlg.Destroy() |
|---|
| 1018 | print '%s%s%s%s'%('finished cell search for ',bravaisNames[ibrav], \ |
|---|
| 1019 | ', elapsed time = ',G2lat.sec2HMS(time.time()-begin)) |
|---|
| 1020 | |
|---|
| 1021 | if cells: |
|---|
| 1022 | return True,dmin,cells |
|---|
| 1023 | else: |
|---|
| 1024 | return False,0,[] |
|---|
| 1025 | |
|---|
| 1026 | |
|---|
| 1027 | NeedTestData = True |
|---|
| 1028 | def TestData(): |
|---|
| 1029 | 'needs a doc string' |
|---|
| 1030 | global NeedTestData |
|---|
| 1031 | NeedTestData = False |
|---|
| 1032 | global TestData |
|---|
| 1033 | TestData = [12, [7.,8.70,10.86,90.,102.95,90.], [7.76006,8.706215,10.865679,90.,102.947,90.],3, |
|---|
| 1034 | [[2.176562137832974, 761.60902227696033, True, True, 0, 0, 1, 10.591300714328161, 10.589436], |
|---|
| 1035 | [3.0477561489789498, 4087.2956049071572, True, True, 1, 0, 0, 7.564238997554908, 7.562777], |
|---|
| 1036 | [3.3254921120068524, 1707.0253890991009, True, True, 1, 0, -1, 6.932650301411212, 6.932718], |
|---|
| 1037 | [3.428121546163426, 2777.5082170150563, True, True, 0, 1, 1, 6.725163158013632, 6.725106], |
|---|
| 1038 | [4.0379791325512118, 1598.4321673135987, True, True, 1, 1, 0, 5.709789097440156, 5.70946], |
|---|
| 1039 | [4.2511182350743937, 473.10955149057577, True, True, 1, 1, -1, 5.423637972781876, 5.42333], |
|---|
| 1040 | [4.354684330373451, 569.88528280256071, True, True, 0, 0, 2, 5.2947091882172534, 5.294718], |
|---|
| 1041 | [4.723324574319177, 342.73882372499997, True, True, 1, 0, -2, 4.881681587039431, 4.881592], |
|---|
| 1042 | [4.9014773581253994, 5886.3516356615492, True, True, 1, 1, 1, 4.704350709093183, 4.70413], |
|---|
| 1043 | [5.0970774474587275, 3459.7541692903033, True, True, 0, 1, 2, 4.523933797797693, 4.523829], |
|---|
| 1044 | [5.2971997607389518, 1290.0229964239879, True, True, 0, 2, 0, 4.353139557169142, 4.353108], |
|---|
| 1045 | [5.4161306205553847, 1472.5726977257755, True, True, 1, 1, -2, 4.257619398422479, 4.257944], |
|---|
| 1046 | [5.7277364698554161, 1577.8791668322888, True, True, 0, 2, 1, 4.026169751907777, 4.026193], |
|---|
| 1047 | [5.8500213058834163, 420.74210142657131, True, True, 1, 0, 2, 3.9420803081518443, 3.942219], |
|---|
| 1048 | [6.0986764166731708, 163.02160537058708, True, True, 2, 0, 0, 3.7814965150452537, 3.781389], |
|---|
| 1049 | [6.1126665157702753, 943.25461245706833, True, True, 1, 2, 0, 3.772849962062199, 3.772764], |
|---|
| 1050 | [6.2559260555056957, 250.55355015505376, True, True, 1, 2, -1, 3.6865353266375283, 3.686602], |
|---|
| 1051 | [6.4226243128279892, 5388.5560141098349, True, True, 1, 1, 2, 3.5909481979190283, 3.591214], |
|---|
| 1052 | [6.5346132446561134, 1951.6070344509026, True, True, 0, 0, 3, 3.5294722429440584, 3.529812], |
|---|
| 1053 | [6.5586952135236443, 259.65938178131034, True, True, 2, 1, -1, 3.516526936765838, 3.516784], |
|---|
| 1054 | [6.6509216222783722, 93.265376597376573, True, True, 2, 1, 0, 3.4678179073694952, 3.468369], |
|---|
| 1055 | [6.7152737044107722, 289.39386813803162, True, True, 1, 2, 1, 3.4346235125812807, 3.434648], |
|---|
| 1056 | [6.8594130457361899, 603.54959764648322, True, True, 0, 2, 2, 3.362534044860622, 3.362553], |
|---|
| 1057 | [7.0511627728884454, 126.43246447656593, True, True, 0, 1, 3, 3.2712038721790675, 3.271181], |
|---|
| 1058 | [7.077700845503319, 125.49742760019636, True, True, 1, 1, -3, 3.2589538988480626, 3.259037], |
|---|
| 1059 | [7.099393757363675, 416.55444885434633, True, True, 1, 2, -2, 3.2490085228959193, 3.248951], |
|---|
| 1060 | [7.1623933278642742, 369.27397110921817, True, True, 2, 1, -2, 3.2204673608202383, 3.220487], |
|---|
| 1061 | [7.4121734953058924, 482.84120827021826, True, True, 2, 1, 1, 3.1120858221599876, 3.112308]] |
|---|
| 1062 | ] |
|---|
| 1063 | global TestData2 |
|---|
| 1064 | TestData2 = [12, [0.15336547830008007, 0.017345499139401827, 0.008122368657493792, 0, 0.02893538955687591, 0], 3, |
|---|
| 1065 | [[2.176562137832974, 761.6090222769603, True, True, 0, 0, 1, 10.591300714328161, 11.095801], |
|---|
| 1066 | [3.0477561489789498, 4087.295604907157, True, True, 0, 1, 0, 7.564238997554908, 7.592881], |
|---|
| 1067 | [3.3254921120068524, 1707.025389099101, True, False, 0, 0, 0, 6.932650301411212, 0.0], |
|---|
| 1068 | [3.428121546163426, 2777.5082170150563, True, True, 0, 1, 1, 6.725163158013632, 6.266192], |
|---|
| 1069 | [4.037979132551212, 1598.4321673135987, True, False, 0, 0, 0, 5.709789097440156, 0.0], |
|---|
| 1070 | [4.251118235074394, 473.10955149057577, True, True, 0, 0, 2, 5.423637972781876, 5.5479], |
|---|
| 1071 | [4.354684330373451, 569.8852828025607, True, True, 0, 0, 2, 5.2947091882172534, 5.199754], |
|---|
| 1072 | [4.723324574319177, 342.738823725, True, False, 0, 0, 0, 4.881681587039431, 0.0], |
|---|
| 1073 | [4.901477358125399, 5886.351635661549, True, False, 0, 0, 0, 4.704350709093183, 0.0], |
|---|
| 1074 | [5.0970774474587275, 3459.7541692903033, True, True, 0, 1, 2, 4.523933797797693, 4.479534], |
|---|
| 1075 | [5.297199760738952, 1290.022996423988, True, True, 0, 1, 0, 4.353139557169142, 4.345087], |
|---|
| 1076 | [5.416130620555385, 1472.5726977257755, True, False, 0, 0, 0, 4.257619398422479, 0.0], |
|---|
| 1077 | [5.727736469855416, 1577.8791668322888, True, False, 0, 0, 0, 4.026169751907777, 0.0], |
|---|
| 1078 | [5.850021305883416, 420.7421014265713, True, False, 0, 0, 0, 3.9420803081518443, 0.0], |
|---|
| 1079 | [6.098676416673171, 163.02160537058708, True, True, 0, 2, 0, 3.7814965150452537, 3.796441], |
|---|
| 1080 | [6.112666515770275, 943.2546124570683, True, False, 0, 0, 0, 3.772849962062199, 0.0], |
|---|
| 1081 | [6.255926055505696, 250.55355015505376, True, True, 0, 0, 3, 3.6865353266375283, 3.6986], |
|---|
| 1082 | [6.422624312827989, 5388.556014109835, True, True, 0, 2, 1, 3.5909481979190283, 3.592005], |
|---|
| 1083 | [6.534613244656113, 191.6070344509026, True, True, 1, 0, -1, 3.5294722429440584, 3.546166], |
|---|
| 1084 | [6.558695213523644, 259.65938178131034, True, True, 0, 0, 3, 3.516526936765838, 3.495428], |
|---|
| 1085 | [6.650921622278372, 93.26537659737657, True, True, 0, 0, 3, 3.4678179073694952, 3.466503], |
|---|
| 1086 | [6.715273704410772, 289.3938681380316, True, False, 0, 0, 0, 3.4346235125812807, 0.0], |
|---|
| 1087 | [6.85941304573619, 603.5495976464832, True, True, 0, 1, 3, 3.362534044860622, 3.32509], |
|---|
| 1088 | [7.051162772888445, 126.43246447656593, True, True, 0, 1, 2, 3.2712038721790675, 3.352121], |
|---|
| 1089 | [7.077700845503319, 125.49742760019636, True, False, 0, 0, 0, 3.2589538988480626, 0.0], |
|---|
| 1090 | [7.099393757363675, 416.55444885434633, True, False, 0, 0, 0, 3.2490085228959193, 0.0], |
|---|
| 1091 | [7.162393327864274, 369.27397110921817, True, False, 0, 0, 0, 3.2204673608202383, 0.0], |
|---|
| 1092 | [7.412173495305892, 482.84120827021826, True, True, 0, 2, 2, 3.112085822159976, 3.133096]] |
|---|
| 1093 | ] |
|---|
| 1094 | # |
|---|
| 1095 | def test0(): |
|---|
| 1096 | if NeedTestData: TestData() |
|---|
| 1097 | ibrav,cell,bestcell,Pwr,peaks = TestData |
|---|
| 1098 | print 'best cell:',bestcell |
|---|
| 1099 | print 'old cell:',cell |
|---|
| 1100 | Peaks = np.array(peaks) |
|---|
| 1101 | HKL = Peaks[4:7] |
|---|
| 1102 | print calc_M20(peaks,HKL) |
|---|
| 1103 | A = G2lat.cell2A(cell) |
|---|
| 1104 | OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr) |
|---|
| 1105 | print 'new cell:',G2lat.A2cell(A) |
|---|
| 1106 | print 'x:',result[0] |
|---|
| 1107 | print 'cov_x:',result[1] |
|---|
| 1108 | print 'infodict:' |
|---|
| 1109 | for item in result[2]: |
|---|
| 1110 | print item,result[2][item] |
|---|
| 1111 | print 'msg:',result[3] |
|---|
| 1112 | print 'ier:',result[4] |
|---|
| 1113 | result = refinePeaks(peaks,ibrav,A) |
|---|
| 1114 | N,M20,X20,A = result |
|---|
| 1115 | print 'refinePeaks:',N,M20,X20,G2lat.A2cell(A) |
|---|
| 1116 | print 'compare bestcell:',bestcell |
|---|
| 1117 | # |
|---|
| 1118 | def test1(): |
|---|
| 1119 | if NeedTestData: TestData() |
|---|
| 1120 | ibrav,A,Pwr,peaks = TestData2 |
|---|
| 1121 | print 'bad cell:',G2lat.A2cell(A) |
|---|
| 1122 | print 'FitHKL' |
|---|
| 1123 | OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr) |
|---|
| 1124 | result = refinePeaks(peaks,ibrav,A) |
|---|
| 1125 | N,M20,X20,A = result |
|---|
| 1126 | print 'refinePeaks:',N,M20,X20,A |
|---|
| 1127 | # Peaks = np.array(peaks) |
|---|
| 1128 | # HKL = Peaks[4:7] |
|---|
| 1129 | # print calc_M20(peaks,HKL) |
|---|
| 1130 | # OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr) |
|---|
| 1131 | # print 'new cell:',G2lat.A2cell(A) |
|---|
| 1132 | # print 'x:',result[0] |
|---|
| 1133 | # print 'cov_x:',result[1] |
|---|
| 1134 | # print 'infodict:' |
|---|
| 1135 | # for item in result[2]: |
|---|
| 1136 | # print item,result[2][item] |
|---|
| 1137 | # print 'msg:',result[3] |
|---|
| 1138 | # print 'ier:',result[4] |
|---|
| 1139 | |
|---|
| 1140 | # |
|---|
| 1141 | if __name__ == '__main__': |
|---|
| 1142 | test0() |
|---|
| 1143 | test1() |
|---|
| 1144 | # test2() |
|---|
| 1145 | # test3() |
|---|
| 1146 | # test4() |
|---|
| 1147 | # test5() |
|---|
| 1148 | # test6() |
|---|
| 1149 | # test7() |
|---|
| 1150 | # test8() |
|---|
| 1151 | print "OK" |
|---|
