Changeset 1163
- Timestamp:
- Dec 7, 2013 3:00:51 PM (9 years ago)
- Location:
- trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIimage.py
r1159 r1163 202 202 R1 = (vplus+vminus)/2. #major axis 203 203 zdis = (fplus-fminus)/2. 204 X = x-x0-zdis*npsind(phi) 204 X = x-x0-zdis*npsind(phi) #shift obsd. ellipses to ellipse center= 0,0 205 205 Y = y-y0-zdis*npcosd(phi) 206 XR = X*npcosd(phi)+Y*npsind(phi) 207 YR = X*npsind(phi)+Y*npcosd(phi) 208 return (XR/R0)**2+(YR/R1)**2-1 206 XR = X*npcosd(phi)+Y*npsind(phi) #rotate by phi to put major axis along x 207 YR = -X*npsind(phi)+Y*npcosd(phi) 208 M = (XR/R0)**2+(YR/R1)**2-1 #=0 for all points exactly on ellipse 209 return M 209 210 210 211 names = ['dist','det-X','det-Y','tilt','phi','dep','wave'] … … 299 300 sinb = sind(tilt) 300 301 d = dist+dxy 301 if tth+ tilt< 90.: #ellipse302 if tth+abs(tilt) < 90.: #ellipse 302 303 fplus = d*tanb*stth/(cosb+stth) 303 304 fminus = d*tanb*stth/(cosb-stth) … … 334 335 'Needs a doc string' 335 336 336 from scipy.optimize import fsolve337 def ellip(xy,*args):338 azm,phi,R0,R1,A,B = args339 cp = cosd(phi)340 sp = sind(phi)341 x,y = xy342 out = []343 out.append(y-x*tand(azm))344 out.append(R0**2*((x+A)*sp-(y+B)*cp)**2+R1**2*((x+A)*cp+(y+B)*sp)**2-(R0*R1)**2)345 return out346 def hyperb(xy,*args):347 azm,phi,R0,R1,A,B = args348 cp = cosd(phi)349 sp = sind(phi)350 x,y = xy351 out = []352 out.append(y+x*tand(azm))353 out.append(R0**2*((x+A)*sp-(y+B)*cp)**2-R1**2*((x+A)*cp+(y+B)*sp)**2-(R0*R1)**2)354 return out355 337 elcent,phi,radii = GetEllipse(dsp,data) 338 phi = data['rotation']-90. #to give rotation of major axis 339 tilt = data['tilt'] 340 dist = data['distance'] 356 341 cent = data['center'] 357 phi = data['rotation'] 358 wave = data['wavelength'] 359 tilt = data['tilt'] 360 dist = data['distance']/cosd(tilt) 361 dep = data['DetDepth'] 362 tth = 2.0*asind(wave/(2.*dsp)) 363 dxy = peneCorr(tth,dep) 342 tth = 2.0*asind(data['wavelength']/(2.*dsp)) 364 343 ttth = tand(tth) 365 radius = (dist+dxy)*ttth366 344 stth = sind(tth) 345 ctth = cosd(tth) 346 sinb = sind(tilt) 367 347 cosb = cosd(tilt) 368 if tth+abs(tilt) < 90.: 369 R1 = (dist+dxy)*stth*cosd(tth)*cosb/(cosb**2-stth**2) 370 R0 = math.sqrt(R1*radius*cosb) 371 zdis = R1*ttth*tand(tilt) 372 A = zdis*sind(phi) 373 B = -zdis*cosd(phi) 374 xy0 = [radius*cosd(azm),radius*sind(azm)] 375 xy = fsolve(ellip,xy0,args=(azm,phi,R0,R1,A,B))+cent 376 else: 377 R1 = (dist+dxy)*stth*cosd(tth)*cosb/(cosb**2+stth**2) 378 R0 = math.sqrt(R1*radius*cosb) 379 zdis = R1*ttth*tand(tilt) 380 A = zdis*sind(phi) 381 B = -zdis*cosd(phi) 382 xy0 = [radius*cosd(azm),radius*sind(azm)] 383 xy = fsolve(hyperb,xy0,args=(azm,phi,R0,R1,A,B))+cent 348 tanb = tand(tilt) 349 if radii[0] > 0.: 350 fplus = dist*tanb*stth/(cosb+stth) 351 fminus = dist*tanb*stth/(cosb-stth) 352 zdis = 0. 353 if tilt: 354 zdis = tilt*(fplus-fminus)/(2.*abs(tilt)) 355 rsqplus = radii[0]**2+radii[1]**2 356 rsqminus = radii[0]**2-radii[1]**2 357 R = rsqminus*cosd(2.*azm-2.*phi)+rsqplus 358 Q = np.sqrt(2.)*radii[0]*radii[1]*np.sqrt(R-2.*zdis**2*sind(azm)**2) 359 P = zdis*(rsqminus*cosd(azm-2.*phi)+rsqplus*cosd(azm)) 360 radius = (P+Q)/R 361 xy = np.array([radius*cosd(azm),radius*sind(azm)]) 362 xy += cent 363 else: #hyperbola - both branches (one is way off screen!) 364 f = dist*tanb*stth/(cosb+stth) 365 v = dist*(tanb+tand(tth-tilt)) 366 delt = dist*stth*(1+stth*cosb)/(sinb*cosb*(stth+cosb)) 367 ecc = (v-f)/(delt-v) 368 R = radii[1]*(ecc**2-1)/(1-ecc*cosd(azm)) 369 xy = [R*cosd(azm)-f,R*sind(azm)] 370 xy = -np.array([xy[0]*cosd(phi)-xy[1]*sind(phi),xy[0]*sind(phi)+xy[1]*cosd(phi)]) 371 xy += cent 384 372 return xy 385 373 … … 409 397 DY = np.sqrt(dx**2+dy**2-Z**2) 410 398 D = (DX**2+DY**2)/dist**2 #for geometric correction = 1/cos(2theta)^2 if tilt=0. 411 tth = npatan d(DY/DX) #depth corr not correct for tilted detector399 tth = npatan2d(DY,DX) #depth corr not correct for tilted detector 412 400 dsp = wave/(2.*npsind(tth/2.)) 413 401 azm = (npatan2d(dy,dx)+azmthoff+720.)%360. … … 637 625 print 'WARNING - selected ring was fitted as a circle' 638 626 print ' - if detector was tilted we suggest you skip this ring - WARNING' 627 tilt = 0.01 639 628 #1st estimate of dist: sample to detector normal to plane 640 629 data['distance'] = dist = radii[0]**2/(ttth*radii[1]) … … 642 631 zdisp = radii[1]*ttth*tand(tilt) 643 632 zdism = radii[1]*ttth*tand(-tilt) 633 print 'zdis',zdisp,zdism 644 634 #cone axis position; 2 choices. Which is right? 645 635 centp = [elcent[0]+zdisp*sind(phi),elcent[1]+zdisp*cosd(phi)] 646 636 centm = [elcent[0]+zdism*sind(phi),elcent[1]+zdism*cosd(phi)] 637 print 'cent',centp,centm 647 638 #check get same ellipse parms either way 648 639 #now do next ring; estimate either way & check sum of Imax/Imin in 3x3 block around each point … … 650 641 tth = 2.0*asind(wave/(2.*dsp)) 651 642 ellipsep = GetEllipse2(tth,0.,dist,centp,tilt,phi) 652 Ringp,delt,sumIp = makeRing(dsp,ellipsep,2,cutoff,scalex,scaley,self.ImageZ) 643 # data['ellipses'].append(copy.deepcopy(ellipsep+('g',))) 644 Ringp,delt,sumIp = makeRing(dsp,ellipsep,3,cutoff,scalex,scaley,self.ImageZ) 653 645 if len(Ringp) > 10: 654 646 outEp,errp = FitRing(Ringp,True) 655 647 else: 656 648 errp = 1e6 657 # print '+',ellipsep,errp649 print '+',sumIp,errp 658 650 ellipsem = GetEllipse2(tth,0.,dist,centm,-tilt,phi) 659 Ringm,delt,sumIm = makeRing(dsp,ellipsem,2,cutoff,scalex,scaley,self.ImageZ) 651 # data['ellipses'].append(copy.deepcopy(ellipsem+('r',))) 652 Ringm,delt,sumIm = makeRing(dsp,ellipsem,3,cutoff,scalex,scaley,self.ImageZ) 660 653 if len(Ringm) > 10: 661 654 outEm,errm = FitRing(Ringm,True) 662 655 else: 663 656 errm = 1e6 664 # print '-',ellipsem,errm 665 if errp < errm: 657 print '-',sumIm,errm 658 # if errp < errm: 659 if sumIp > sumIm: 666 660 data['tilt'] = tilt 667 661 data['center'] = centp -
trunk/GSASIIimgGUI.py
r1158 r1163 464 464 465 465 def OnLRazim(event): 466 Lazm =int(G2frame.Lazim.GetValue()) 466 Lazm = int(G2frame.Lazim.GetValue())%360 467 Razm = int(G2frame.Razim.GetValue())%360 468 if Lazm > Razm: 469 Razm += 360 467 470 if data['fullIntegrate']: 468 G2frame.Razim.SetValue("%6d" % (Lazm+360)) 469 Razm = int(G2frame.Razim.GetValue()) 470 if Lazm > Razm: 471 Lazm -= 360 471 Razm = Lazm+360 472 G2frame.Lazim.SetValue("%6d" % (Lazm)) 473 G2frame.Razim.SetValue("%6d" % (Razm)) 472 474 data['LRazimuth'] = [Lazm,Razm] 473 475 G2plt.PlotExposedImage(G2frame,event=event) -
trunk/GSASIIplot.py
r1153 r1163 2323 2323 elif 'line3' in itemPicked: 2324 2324 Data['LRazimuth'][0] = int(azm) 2325 if Data['fullIntegrate']:2326 Data['LRazimuth'][1] = Data['LRazimuth'][0]+3602327 2325 elif 'line4' in itemPicked and not Data['fullIntegrate']: 2328 2326 Data['LRazimuth'][1] = int(azm) 2329 2327 2328 Data['LRazimuth'][0] %= 360 2329 Data['LRazimuth'][1] %= 360 2330 2330 if Data['LRazimuth'][0] > Data['LRazimuth'][1]: 2331 Data['LRazimuth'][0] -= 360 2332 2333 azLim = np.array(Data['LRazimuth']) 2334 if np.any(azLim>360): 2335 azLim -= 360 2336 Data['LRazimuth'] = list(azLim) 2331 Data['LRazimuth'][1] += 360 2332 if Data['fullIntegrate']: 2333 Data['LRazimuth'][1] = Data['LRazimuth'][0]+360 2337 2334 2338 2335 if Data['IOtth'][0] > Data['IOtth'][1]: … … 2494 2491 xyI = [] 2495 2492 for azm in Azm: 2496 xyI.append(G2img.GetDetectorXY(dspI,azm,Data)) #what about hyperbola? 2497 xyI = np.array(xyI) 2498 arcxI,arcyI = xyI.T 2499 Plot.plot(arcxI,arcyI,picker=3) 2493 xy = G2img.GetDetectorXY(dspI,azm,Data) 2494 if np.any(xy): 2495 xyI.append(xy) 2496 if len(xyI): 2497 xyI = np.array(xyI) 2498 arcxI,arcyI = xyI.T 2499 Plot.plot(arcxI,arcyI,picker=3) 2500 2500 if ellO: 2501 2501 xyO = [] 2502 2502 for azm in Azm: 2503 xyO.append(G2img.GetDetectorXY(dspO,azm,Data)) #what about hyperbola? 2504 xyO = np.array(xyO) 2505 arcxO,arcyO = xyO.T 2506 Plot.plot(arcxO,arcyO,picker=3) 2503 xy = G2img.GetDetectorXY(dspO,azm,Data) 2504 if np.any(xy): 2505 xyO.append(xy) 2506 if len(xyO): 2507 xyO = np.array(xyO) 2508 arcxO,arcyO = xyO.T 2509 Plot.plot(arcxO,arcyO,picker=3) 2507 2510 if ellO and ellI: 2508 2511 Plot.plot([arcxI[0],arcxO[0]],[arcyI[0],arcyO[0]],picker=3)
Note: See TracChangeset
for help on using the changeset viewer.