calculated hyperbola rings correctly plotted, but no observed points. No fitting of hyperbola rings possible - wrong math for locating points & for FitDetector - ellipse only.

Author: vondreele

GSASII/GSASIIimage.py

@@ -97,7 +97,7 @@ def ellipse_axis_length( p ):
         down2=(b*b-a*c)*( (a-c)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
         res1=np.sqrt(up/down1)
         res2=np.sqrt(up/down2)
-        return np.array([ res2,res1])
+        return np.array([ res2,res1,0.0])
 
     xy = np.array(xy)
     x = np.asarray(xy.T[0])[:,np.newaxis]
@@ -371,22 +371,23 @@ def ellipseC():
             if debug: 
                 theta = npacosd(1./np.sqrt(1.+(radii[0]/radii[1])**2))
                 print ('hyperbola:',theta)
-            return 200.*np.pi
+            return 720.
         apb = radii[1]+radii[0]
         amb = radii[1]-radii[0]
         return np.pi*apb*(1+3*(amb/apb)**2/(10+np.sqrt(4-3*(amb/apb)**2)))
-
     cent,phi,radii = ellipse
+    if radii[1] < 0.:
+        return None
     cphi = cosd(phi-90.)        #convert to major axis rotation
     sphi = sind(phi-90.)
     ring = []
     C = int(ellipseC())*mul         #ring circumference in mm
     azm = []
     for i in range(0,C,1):      #step around ring in 1mm increments
         a = 360.*i/C
-        if radii[1] <= 0:        #parabola or hyperbola
-            x = radii[1]*coshd(a-phi+90.)        #major axis
-            y = radii[0]*sinhd(a-phi+90.)
+        if radii[0] <= 0:        #parabola or hyperbola
+            x = -radii[0]/cosd(a-phi)        #major axis
+            y = radii[1]*tand(a-phi)
         else:
             x = radii[1]*cosd(a-phi+90.)        #major axis
             y = radii[0]*sind(a-phi+90.)
@@ -412,7 +413,7 @@ def GetEllipse2(tth,dxy,dist,cent,tilt,phi):
     radii[0] (b-minor axis) set < 0. for hyperbola
 
     '''
-    radii = [0,0]
+    radii = [0,0,0]
     stth = sind(tth)
     cosb = cosd(tilt)
     tanb = tand(tilt)
@@ -427,6 +428,7 @@ def GetEllipse2(tth,dxy,dist,cent,tilt,phi):
         vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
         radii[0] = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2.      #+minor axis
         radii[1] = (vplus+vminus)/2.                                    #major axis
+        radii[2] = tth                                                  #save for ellipse; might be useful
         zdis = (fplus-fminus)/2.
     else:   #hyperbola!
         f = d*abs(tanb)*stth/(cosb+stth)
@@ -435,10 +437,11 @@ def GetEllipse2(tth,dxy,dist,cent,tilt,phi):
         eps = (v-f)/(delt-v)
         radii[0] = -eps*(delt-f)/np.sqrt(eps**2-1.)                     #-minor axis
         radii[1] = eps*(delt-f)/(eps**2-1.)                             #major axis
+        radii[2] = tth                                                  #save 2-theta for hyperbola
         if tilt > 0:
             zdis = f+radii[1]*eps
         else:
-            zdis = -f
+            zdis = -f-radii[1]*eps
 #NB: zdis is || to major axis & phi is rotation of minor axis
 #thus shift from beam to ellipse center is [Z*sin(phi),-Z*cos(phi)]
     elcent = [cent[0]+zdis*sind(phi),cent[1]-zdis*cosd(phi)]
@@ -457,21 +460,40 @@ def GetEllipse(dsp,data):
     dxy = peneCorr(tth,dep,dist)
     return GetEllipse2(tth,dxy,dist,cent,tilt,phi)
 
+def LinePlaneCollision(planeNormal, planePoint, rayDirection, rayPoint, epsilon=1e-6):
+
+	ndotu = planeNormal.dot(rayDirection)
+	if ndotu < epsilon:
+		return None
+	w = rayPoint - planePoint
+	si = -planeNormal.dot(w) / ndotu
+	Psi = w + si * rayDirection + planePoint
+	return Psi
+
+def GetDetSamAngle(x,y,data):
+    def costth(xyz,d0):
+        ''' compute cos of angle between vectors; xyz not normalized, d0 normalized'''
+        u = xyz/nl.norm(xyz,axis=-1)[:,:,nxs]
+        return np.dot(u,d0)
+#zero detector 2-theta: tested with tilted images - perfect integrations
+    dx = x-data['center'][0]
+    dy = y-data['center'][1]
+    tilt = data['tilt']
+    dist = data['distance']/npcosd(tilt)    #sample-beam intersection point on detector plane
+    T = makeMat(tilt,0)         #detector tilt matrix
+    R = makeMat(data['rotation'],2)     #rotation of tilt axis matrix
+    MN = np.inner(R,np.inner(R,T))      #should be detector transformation matrix; why not np.inner(R,T)
+    d001 = np.array([0.,0.,1.])         #vector along z (beam direction); normal to untilted detector plane
+    r001 = np.inner(d001,MN)            #should rotate vector same as detector
+    dxyz0 = np.inner(np.dstack([dx,dy,np.zeros_like(dx)]),MN)    #transform detector pixel x,y by tilt/rotate
+    dxyz0 += np.array([0.,0.,dist])         #shift away from sample
+    ctth0 = costth(dxyz0,r001)              #cos of angle between detector normal & sample-pixel vector
+    return ctth0
+
 def GetDetectorXY(dsp,azm,data):
     '''Get detector x,y position from d-spacing (dsp), azimuth (azm,deg)
     & image controls dictionary (data) - new version
     '''
-    
-    def LinePlaneCollision(planeNormal, planePoint, rayDirection, rayPoint, epsilon=1e-6):
-
-    	ndotu = planeNormal.dot(rayDirection)
-    	if ndotu < epsilon:
-    		return None
-    	w = rayPoint - planePoint
-    	si = -planeNormal.dot(w) / ndotu
-    	Psi = w + si * rayDirection + planePoint
-    	return Psi
-
     dist = data['distance']/npcosd(data['tilt'])    #sample-beam intersection point on detector plane
     cent = data['center']
     phi = data['rotation']
@@ -699,7 +721,7 @@ def EdgeFinder(image,data):
     return zip(tax,tay)
 
 def CalcRings(G2frame,ImageZ,data,masks):
-    ''' Not used anywhere?'''
+    ''' '''
     pixelSize = data['pixelSize']
     scalex = 1000./pixelSize[0]
     scaley = 1000./pixelSize[1]
@@ -730,7 +752,6 @@ def CalcRings(G2frame,ImageZ,data,masks):
     else:
         absent = ()
     HKL = G2lat.sortHKLd(HKL,True,False)
-    wave = data['wavelength']
     frame = masks['Frames']
     tam = ma.make_mask_none(ImageZ.shape)
     if frame:
@@ -808,6 +829,7 @@ def ImageRecalibrate(G2frame,ImageZ,data,masks,getRingsOnly=False):
     tam = ma.make_mask_none(ImageZ.shape)
     if frame:
         tam = ma.mask_or(tam,ma.make_mask(np.abs(polymask(data,[frame,])-255)))
+    hyperbola = False
     for iH,H in enumerate(HKL):
         if debug:   print (H)
         dsp = H[3]
@@ -826,11 +848,14 @@ def ImageRecalibrate(G2frame,ImageZ,data,masks,getRingsOnly=False):
             continue
         else:                   #no more rings beyond edge of detector
             data['ellipses'].append([])
+            if Ring is None:
+                hyperbola = True
             continue
     if not data['rings']:
         G2fil.G2Print ('no rings found; try lower Min ring I/Ib',mode='warn')
         return []
-
+    if hyperbola:
+        print('Hyperbola found, outer rings not fitted')
     rings = np.concatenate((data['rings']),axis=0)
     if getRingsOnly:
         return rings,HKL
@@ -1516,8 +1541,6 @@ def ImageIntegrate(image,data,masks,blkSize=128,returnN=False,useTA=None,useMask
         H1 = LRazm
     if 'SASD' not in data['type']:
         H0 *= np.array(G2pwd.Polarization(data['PolaVal'][0],H2[:-1],0.)[0])
-    # if np.abs(data['det2theta']) < 1.0:         #small angle approx only; not appropriate for detectors at large 2-theta
-        # H0 /= np.abs(npcosd(H2[:-1]+data['tilt']-np.abs(data['det2theta'])))**4           #parallax correction (why **4?)
     if 'SASD' in data['type']:
         H0 /= npcosd(H2[:-1])           #one more for small angle scattering data?
     if data['Oblique'][1]:

GSASII/GSASIIplot.py

@@ -5511,13 +5511,33 @@ def OnImRelease(event):
                     Plot.plot(xring,yring,'.',color=colors[N%NC])
                     N += 1
             for ellipse in Data['ellipses']:      #what about hyperbola?
-                cent,phi,[width,height],col = ellipse
+                try:
+                    cent,phi,[width,height,tth],col = ellipse
+                except ValueError:
+                    cent,phi,[width,height],col = ellipse
+                    tth = 0.0
                 if width > 0:       #ellipses
                     try:  # angle was changed to a keyword at some point, needed in mpl 3.8
                         Plot.add_artist(Ellipse([cent[0],cent[1]],2*width,2*height,angle=phi,ec=col,fc='none'))
                     except: # but keep the old version as a patch (5/20/24) in case old call needed for old MPL
                         Plot.add_artist(Ellipse([cent[0],cent[1]],2*width,2*height,phi,ec=col,fc='none'))
                     Plot.text(cent[0],cent[1],'+',color=col,ha='center',va='center')
+                elif tth:       #future hyperbola plot
+                    dsp =0.5*Data['wavelength']/npsind(tth/2.0)
+                    darc = Data['rotation']
+#                    Azm = np.arange(-10.-darc,190.5-darc,.5)
+                    Azm = np.arange(0.,360.5,.5)
+                    xyH = []
+                    for azm in Azm:
+                        xy = G2img.GetDetectorXY(dsp,azm,Data)
+                        if np.any(xy):
+                            xyH.append(xy)
+                    if len(xyH):
+                        xyH = np.array(xyH)
+                        xH,yH = xyH.T
+                        Plot.plot(xH,yH,color=col)
+                    Plot.text(cent[0],cent[1],'+',color=col,ha='center',va='center')
+                    continue
         if G2frame.PickId and G2frame.GPXtree.GetItemText(G2frame.PickId) in ['Stress/Strain',]:
             for N,ring in enumerate(StrSta['d-zero']):
                 if 'ImxyCalc' in ring: