Changeset 1148


Ignore:
Timestamp:
Nov 23, 2013 5:29:28 PM (9 years ago)
Author:
vondreele
Message:

major changes to 2-D detector calibration; now works for strongly tilted detectors & short sample to detector distances.
Distance is now defined as sample to detector plane. Previously it was sample to intercept of detector plane with incident beam (Bragg cone axis).
The "penetration" parameter is still suspect.

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIIO.py

    r1147 r1148  
    410410    File = open(filename,'rb')
    411411    dataType = 5
     412    center = [None,None]
     413    wavelength = None
     414    distance = None
    412415    try:
    413416        Meta = open(filename+'.metadata','Ur')
    … …  
    464467        pixy = (marFrame.pixelsizeX/1000.0,marFrame.pixelsizeY/1000.0)
    465468        head = marFrame.outputHead()
     469# extract resonable wavelength from header
     470        wavelength = marFrame.sourceWavelength*1e-5
     471        wavelength = (marFrame.opticsWavelength > 0) and marFrame.opticsWavelength*1e-5 or wavelength
     472        wavelength = (wavelength <= 0) and None or wavelength
     473# extract resonable distance from header
     474        distance = (marFrame.startXtalToDetector+marFrame.endXtalToDetector)*5e-4
     475        distance = (distance <= 0) and marFrame.xtalToDetector*1e-3 or distance
     476        distance = (distance <= 0) and None or distance
     477# extract resonable center from header
     478        center = [marFrame.beamX*marFrame.pixelsizeX*1e-9,marFrame.beamY*marFrame.pixelsizeY*1e-9]
     479        center = (center[0] != 0 and center[1] != 0) and center or [None,None]
     480#print head,tifType,pixy
    466481    elif 272 in IFD:
    467482        ifd = IFD[272]
    … …  
    561576       
    562577    image = np.reshape(image,(sizexy[1],sizexy[0]))
    563     center = [pixy[0]*sizexy[0]/2000,pixy[1]*sizexy[1]/2000]
    564     data = {'pixelSize':pixy,'wavelength':0.10,'distance':100.0,'center':center,'size':sizexy}
     578    center = (not center[0]) and [pixy[0]*sizexy[0]/2000,pixy[1]*sizexy[1]/2000] or center
     579    wavelength = (not wavelength) and 0.10 or wavelength
     580    distance = (not distance) and 100.0 or distance
     581    data = {'pixelSize':pixy,'wavelength':wavelength,'distance':distance,'center':center,'size':sizexy}
    565582    File.close()   
    566583    if imageOnly:

  • trunk/GSASIIimage.py

    r1133 r1148  
    1515
    1616needs some minor refactoring to remove wx code
     17actually not easy in this case wx.ProgressDialog needs #of blocks to process when started
     18not known in G2imageGUI.
    1719'''
    1820
    … …  
    4446npasind = lambda x: 180.*np.arcsin(x)/np.pi
    4547npcosd = lambda x: np.cos(x*np.pi/180.)
     48npacosd = lambda x: 180.*np.arccos(x)/np.pi
    4649nptand = lambda x: np.tan(x*np.pi/180.)
    4750npatand = lambda x: 180.*np.arctan(x)/np.pi
    … …  
    9194    else:
    9295        err,parms = FitCircle(ring)
    93     return parms
     96    return parms,err
    9497       
    9598def FitCircle(ring):
    … …  
    118121        if det < 0.:
    119122            print 'hyperbola!'
     123            print B
    120124            return 0
    121125        elif det == 0.:
    122126            print 'parabola!'
     127            print B
    123128            return 0
    124129        cent = [(B[1]*B[3]-2.*(1.-B[0])*B[2])/det, \
    … …  
    150155    return err,makeParmsEllipse(bb)
    151156   
    152 def FitDetector(rings,varyList,parmDict):
     157def FitDetector(rings,varyList,parmDict,Print=True):
    153158    'Needs a doc string'
    154159       
    … …  
    173178       
    174179    def ellipseCalcD(B,xyd,varyList,parmDict):
    175         x = xyd[0]
    176         y = xyd[1]
    177         dsp = xyd[2]
     180        x,y,dsp = xyd
    178181        wave = parmDict['wave']
    179182        if 'dep' in varyList:
    … …  
    187190        dxy = peneCorr(tth,dep)
    188191        ttth = nptand(tth)
    189         radius = (dist+dxy)*ttth
    190192        stth = npsind(tth)
    191193        cosb = npcosd(tilt)
    192         R1 = (dist+dxy)*stth*npcosd(tth)*cosb/(cosb**2-stth**2)
    193         R0 = np.sqrt(R1*radius*cosb)
    194         zdis = R1*ttth*nptand(tilt)
    195         X = x-x0+zdis*npsind(phi)
     194        tanb = nptand(tilt)       
     195        tbm = nptand((tth-tilt)/2.)
     196        tbp = nptand((tth+tilt)/2.)
     197        sinb = npsind(tilt)
     198        d = dist+dxy
     199        fplus = d*tanb*stth/(cosb+stth)
     200        fminus = d*tanb*stth/(cosb-stth)
     201        vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
     202        vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
     203        R0 = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2.      #+minor axis
     204        R1 = (vplus+vminus)/2.                                    #major axis
     205        zdis = (fplus-fminus)/2.
     206        X = x-x0-zdis*npsind(phi)
    196207        Y = y-y0-zdis*npcosd(phi)
    197         XR = X*npcosd(phi)-Y*npsind(phi)
     208        XR = X*npcosd(phi)+Y*npsind(phi)
    198209        YR = X*npsind(phi)+Y*npcosd(phi)
    199210        return (XR/R0)**2+(YR/R1)**2-1
    … …  
    202213    fmt = ['%12.2f','%12.2f','%12.2f','%12.2f','%12.2f','%12.3f','%12.5f']
    203214    p0 = [parmDict[key] for key in varyList]
    204     result = leastsq(ellipseCalcD,p0,args=(rings.T,varyList,parmDict),full_output=True)
     215    result = leastsq(ellipseCalcD,p0,args=(rings.T,varyList,parmDict),full_output=True,ftol=1.e-8)
     216    chisq = np.sum(result[2]['fvec']**2)
    205217    parmDict.update(zip(varyList,result[0]))
    206218    vals = list(result[0])
    … …  
    212224            sigList[i] = sig[varyList.index(name)]
    213225    ValSig = zip(names,fmt,vals,sig)
    214     CalibPrint(ValSig)
    215 #    try:
    216 #        print 'Trial refinement of wavelength - not used for calibration'
    217 #        p0 = result[0]
    218 #        print p0
    219 #        print parms
    220 #        p0 = np.append(p0,parms[0])
    221 #        resultW = leastsq(ellipseCalcW,p0,args=(rings.T,parms[1:]),full_output=True)
    222 #        resultW[0][3] = np.mod(result[0][3],360.0)          #remove random full circles
    223 #        sig = np.sqrt(np.sum(ellipseCalcW(resultW[0],rings.T)**2))
    224 #        ValSig = zip(names,fmt,resultW[0],sig*np.sqrt(np.diag(resultW[1])))
    225 #        CalibPrint(ValSig)
    226 #        return result[0],resultW[0][-1]       
    227 #    except ValueError:
    228 #        print 'Bad refinement - no result'
    229 #        return result[0],wave
     226    if Print:
     227        CalibPrint(ValSig)
     228    return chi,chisq
    230229                   
    231230def ImageLocalMax(image,w,Xpix,Ypix):
    … …  
    247246def makeRing(dsp,ellipse,pix,reject,scalex,scaley,image):
    248247    'Needs a doc string'
     248    def ellipseC():
     249        'compute estimate of ellipse circumference'
     250        if radii[0] < 0:        #hyperbola
     251            theta = npacosd(1./np.sqrt(1.+(radii[0]/radii[1])**2))
     252            print theta
     253            return 0
     254        apb = radii[1]+radii[0]
     255        amb = radii[1]-radii[0]
     256        return np.pi*apb*(1+3*(amb/apb)**2/(10+np.sqrt(4-3*(amb/apb)**2)))
    249257    cent,phi,radii = ellipse
    250258    cphi = cosd(phi)
    251259    sphi = sind(phi)
    252260    ring = []
     261    sumI = 0
    253262    amin = 0
    254263    amax = 360
    255     for a in range(0,360,1):
     264    for a in range(0,int(ellipseC()),1):
    256265        x = radii[0]*cosd(a)
    257266        y = radii[1]*sind(a)
    … …  
    260269        X,Y,I,J = ImageLocalMax(image,pix,X,Y)
    261270        if I and J and I/J > reject:
     271            sumI += I/J
    262272            X += .5                             #set to center of pixel
    263273            Y += .5
    … …  
    269279                ring.append([X,Y,dsp])
    270280    delt = amax-amin
    271     if len(ring) < 10:             #want more than 20 deg
    272         return [],delt > 90
    273     return ring,delt > 90
     281    if len(ring) < 10:             #want more than 10 deg
     282        return [],(delt > 90),0
     283    return ring,(delt > 90),sumI/len(ring)
    274284   
    275285def makeIdealRing(ellipse,azm=None):
    … …  
    286296       
    287297    a = np.linspace(aR[0],aR[1],aR[2])
    288     x = radii[0]*npcosd(a-phi)
    289     y = radii[1]*npsind(a-phi)
    290     X = (cphi*x-sphi*y+cent[0])
    291     Y = (sphi*x+cphi*y+cent[1])
    292     return zip(X,Y)
    293                
     298    slr = radii[0]**2/radii[1]
     299    rat2 = (radii[0]/radii[1])**2
     300    if radii[0] > 0.:       #ellipse
     301        ecc = np.sqrt(max(1.-rat2,0.))
     302    else:                   #hyperbola
     303        ecc = np.sqrt(1.+rat2)
     304    rad = slr/(1.+ecc*npcosd(a))
     305    xy = np.array([rad*npsind(a)+cent[0],rad*npcosd(a)+cent[1]])
     306    return xy
     307                   
    294308def calcDist(radii,tth):
    295309    'Needs a doc string'
    … …  
    309323        zdis = radii[1]*ttth*cosb/sinb
    310324        return zdis,cosB
     325       
     326def GetEllipse2(tth,dxy,dist,cent,tilt,phi):
     327    '''uses Dandelin spheres to find ellipse or hyperbola parameters from detector geometry
     328    on output
     329    radii[0] (b-minor axis) set < 0. for hyperbola
     330   
     331    '''
     332    radii = [0,0]
     333    ttth = tand(tth)
     334    stth = sind(tth)
     335    ctth = cosd(tth)
     336    cosb = cosd(tilt)
     337    tanb = tand(tilt)
     338    tbm = tand((tth-tilt)/2.)
     339    tbp = tand((tth+tilt)/2.)
     340    sinb = sind(tilt)
     341    d = dist+dxy
     342    if tth+tilt < 90.:      #ellipse
     343        fplus = d*tanb*stth/(cosb+stth)
     344        fminus = d*tanb*stth/(cosb-stth)
     345        vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
     346        vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
     347        radii[0] = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2.      #+minor axis
     348        radii[1] = (vplus+vminus)/2.                                    #major axis
     349        zdis = (fplus-fminus)/2.
     350    else:   #hyperbola!
     351        f = d*tanb*stth/(cosb+stth)
     352        v = d*(tanb+tand(tth-tilt))
     353        delt = d*stth*(1+stth*cosb)/(sinb*cosb*(stth+cosb))
     354        eps = (v-f)/(delt-v)
     355        radii[0] = -eps*(delt-f)/np.sqrt(eps**2-1.)                     #-minor axis
     356        radii[1] = eps*(delt-f)/(eps**2-1.)                             #major axis
     357        zdis = f+radii[1]*eps
     358    elcent = [cent[0]+zdis*sind(phi),cent[1]+zdis*cosd(phi)]
     359    return elcent,phi,radii
    311360   
    312361def GetEllipse(dsp,data):
    313     'Needs a doc string'
    314     dist = data['distance']
     362    '''uses Dandelin spheres to find ellipse or hyperbola parameters from detector geometry
     363    as given in image controls dictionary (data) and a d-spacing (dsp)
     364    '''
    315365    cent = data['center']
    316366    tilt = data['tilt']
    317367    phi = data['rotation']
    318368    dep = data['DetDepth']
    319     radii = [0,0]
    320369    tth = 2.0*asind(data['wavelength']/(2.*dsp))
    321     ttth = tand(tth)
    322     stth = sind(tth)
    323     ctth = cosd(tth)
    324     cosb = cosd(tilt)
    325370    dxy = peneCorr(tth,dep)
    326     radius = ttth*(dist+dxy)
    327     radii[1] = (dist+dxy)*stth*ctth*cosb/(cosb**2-stth**2)
    328     if radii[1] > 0:
    329         radii[0] = math.sqrt(radii[1]*radius*cosb)
    330         zdis = radii[1]*ttth*tand(tilt)
    331         elcent = [cent[0]-zdis*sind(phi),cent[1]+zdis*cosd(phi)]
    332         return elcent,phi,radii
    333     else:
    334         print 'bad ellipse - radii:',radii
    335         return False
     371    dist = data['distance']
     372    return GetEllipse2(tth,dxy,dist,cent,tilt,phi)
    336373       
    337374def GetDetectorXY(dsp,azm,data):
    338375    'Needs a doc string'
     376   
    339377    from scipy.optimize import fsolve
    340378    def func(xy,*args):
    … …  
    349387    elcent,phi,radii = GetEllipse(dsp,data)
    350388    cent = data['center']
    351     tilt = data['tilt']
    352389    phi = data['rotation']
    353390    wave = data['wavelength']
    354     dist = data['distance']
     391    tilt = data['tilt']
     392    dist = data['distance']/cosd(tilt)
    355393    dep = data['DetDepth']
    356394    tth = 2.0*asind(wave/(2.*dsp))
    … …  
    371409def GetDetXYfromThAzm(Th,Azm,data):
    372410    'Needs a doc string'
    373     dsp = data['wavelength']/(2.0*npsind(Th))
    374    
     411    dsp = data['wavelength']/(2.0*npsind(Th))   
    375412    return GetDetectorXY(dsp,azm,data)
    376413                   
    377414def GetTthAzmDsp(x,y,data):
    378     'Needs a doc string'
     415    'Needs a doc string - checked OK for ellipses dont know about hyperbola'
    379416    wave = data['wavelength']
    380     dist = data['distance']
    381417    cent = data['center']
    382418    tilt = data['tilt']
     419    dist = data['distance']/cosd(tilt)
    383420    phi = data['rotation']
    384421    dep = data['DetDepth']
    … …  
    392429    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
    393430    dxy = peneCorr(tth,dep)
    394     tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z+dxy))
     431    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z+dxy)) #depth corr not correct for tilted detector
    395432    dsp = wave/(2.*npsind(tth/2.))
    396     azm = (npatan2d(dx,-dy)+azmthoff+720.)%360.
     433    azm = (npatan2d(dy,dx)+azmthoff+720.)%360.
    397434    return tth,azm,dsp
    398435   
    … …  
    428465def EdgeFinder(image,data):
    429466    '''this makes list of all x,y where I>edgeMin suitable for an ellipse search?
     467    Not currently used but might be useful in future?
    430468    '''
    431469    import numpy.ma as ma
    … …  
    465503    for bravais,cell in zip(Bravais,Cells):
    466504        A = G2lat.cell2A(cell)
    467         hkl = G2lat.GenHBravais(dmin,bravais,A)[skip:]
     505        hkl = G2lat.GenHBravais(dmin,bravais,A)
    468506        HKL += hkl
    469507    HKL = G2lat.sortHKLd(HKL,True,False)
    … …  
    471509    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
    472510        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
     511    Found = False
    473512    for H in HKL:
    474513        dsp = H[3]
    475514        ellipse = GetEllipse(dsp,data)
    476         Ring,delt = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
     515        Ring,delt,sumI = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
    477516        if Ring:
    478 #            numZ = len(Ring)
    479517            data['rings'].append(np.array(Ring))
    480518            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
    481         else:
     519            Found = True
     520        elif not Found:         #skipping inner rings, keep looking until ring found
    482521            continue
     522        else:                   #no more rings beyond edge of detector
     523            break
    483524    rings = np.concatenate((data['rings']),axis=0)
    484525    if data['DetDepthRef']:
    … …  
    518559    #fit start points on inner ring
    519560    data['ellipses'] = []
    520     outE = FitRing(ring,True)
     561    outE,err = FitRing(ring,True)
     562    fmt  = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f'
     563    fmt2 = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f, chi^2: %.3f, N: %d'
    521564    if outE:
    522         print 'start ellipse:',outE
     565        print fmt%('start ellipse: ',outE[0][0],outE[0][1],outE[1],outE[2][0],outE[2][1])
    523566        ellipse = outE
    524567    else:
    … …  
    526569       
    527570    #setup 360 points on that ring for "good" fit
    528     Ring,delt = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
     571    Ring,delt,sumI = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
    529572    if Ring:
    530         ellipse = FitRing(Ring,delt)
    531         Ring,delt = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)    #do again
    532         ellipse = FitRing(Ring,delt)
     573        ellipse,err = FitRing(Ring,delt)
     574        Ring,delt,sumI = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)    #do again
     575        ellipse,err = FitRing(Ring,delt)
    533576    else:
    534577        print '1st ring not sufficiently complete to proceed'
    535578        return False
    536     print 'inner ring:',ellipse     #cent,phi,radii
    537     data['center'] = copy.copy(ellipse[0])           #not right!! (but useful for now)
     579    print fmt2%('inner ring:    ',ellipse[0][0],ellipse[0][1],ellipse[1],ellipse[2][0],ellipse[2][1],err,len(Ring))     #cent,phi,radii
    538580    data['ellipses'].append(ellipse[:]+('r',))
     581    data['rings'].append(np.array(Ring))
    539582    G2plt.PlotImage(self,newImage=True)
    540583   
    … …  
    542585    data['rings'] = []
    543586    data['ellipses'] = []
    544     if not data['calibrant']:
     587    if not data['calibrant']  or 'None' in data['calibrant']:
    545588        print 'no calibration material selected'
    546589        return True
    … …  
    548591    skip = data['calibskip']
    549592    dmin = data['calibdmin']
     593#generate reflection set
    550594    Bravais,Cells = calFile.Calibrants[data['calibrant']][:2]
    551595    HKL = []
    … …  
    556600    HKL = G2lat.sortHKLd(HKL,True,False)
    557601    wave = data['wavelength']
    558     cent = data['center']
    559     elcent,phi,radii = ellipse
     602#set up 1st ring
     603    elcent,phi,radii = ellipse              #from fit of 1st ring
    560604    dsp = HKL[0][3]
    561605    tth = 2.0*asind(wave/(2.*dsp))
    562     ttth = tand(tth)
    563     data['distance'] = dist = calcDist(radii,tth)
    564     radius = dist*tand(tth)
    565     zdis,cosB = calcZdisCosB(radius,tth,radii)
    566     cent1 = []
    567     cent2 = []
    568     xSum = 0
    569     ySum = 0
    570     zxSum = 0
    571     zySum = 0
    572     phiSum = 0
    573     tiltSum = 0
    574     distSum = 0
    575     Zsum = 0
     606    ttth = nptand(tth)
     607    stth = npsind(tth)
     608    ctth = npcosd(tth)
     609#1st estimate of tilt; assume ellipse - don't know sign though
     610    tilt = npasind(np.sqrt(max(0.,1.-(radii[0]/radii[1])**2))*ctth)
     611#1st estimate of dist: sample to detector normal to plane
     612    data['distance'] = dist = radii[0]**2/(ttth*radii[1])
     613#ellipse to cone axis (x-ray beam); 2 choices depending on sign of tilt
     614    zdisp = radii[1]*ttth*tand(tilt)
     615    zdism = radii[1]*ttth*tand(-tilt)
     616#cone axis position; 2 choices. Which is right?
     617    centp = [elcent[0]+zdisp*sind(phi),elcent[1]+zdisp*cosd(phi)]
     618    centm = [elcent[0]+zdism*sind(phi),elcent[1]+zdism*cosd(phi)]
     619#check get same ellipse parms either way
     620#now do next ring; estimate either way & check sum of Imax/Imin in 3x3 block around each point
     621    dsp = HKL[1][3]
     622    tth = 2.0*asind(wave/(2.*dsp))
     623    ellipsep = GetEllipse2(tth,0.,dist,centp,tilt,phi)
     624    Ringp,delt,sumIp = makeRing(dsp,ellipsep,2,cutoff,scalex,scaley,self.ImageZ)
     625    outEp,errp = FitRing(Ringp,True)
     626    print '+',ellipsep,errp
     627    ellipsem = GetEllipse2(tth,0.,dist,centm,-tilt,phi)
     628    Ringm,delt,sumIm = makeRing(dsp,ellipsem,2,cutoff,scalex,scaley,self.ImageZ)
     629    outEm,errm = FitRing(Ringm,True)
     630    print '-',ellipsem,errm
     631    if errp < errm:
     632        data['tilt'] = tilt
     633        data['center'] = centp
     634        negTilt = 1
     635    else:
     636        data['tilt'] = -tilt
     637        data['center'] = centm
     638        negTilt = -1
     639    data['rotation'] = phi
     640    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
     641        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
     642    varyList = ['dist','det-X','det-Y','tilt','phi']
     643    if data['DetDepthRef']:
     644        varyList.append('dep')
    576645    for i,H in enumerate(HKL):
    577646        dsp = H[3]
    578         tth = 2.0*asind(0.5*wave/dsp)
    579         stth = sind(tth)
    580         ctth = cosd(tth)
    581         ttth = tand(tth)
    582         radius = dist*ttth
    583         elcent,phi,radii = ellipse
    584         if i:
    585             radii[1] = dist*stth*ctth*cosB/(cosB**2-stth**2)
    586             radii[0] = math.sqrt(radii[1]*radius*cosB)
    587             zdis,cosB = calcZdisCosB(radius,tth,radii)
    588             zsinp = zdis*sind(phi)
    589             zcosp = zdis*cosd(phi)
    590             cent = data['center']
    591             elcent = [cent[0]-zsinp,cent[1]+zcosp]
    592             ellipse = (elcent,phi,radii)
    593         ratio = radii[1]/radii[0]
    594         Ring,delt = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
     647        tth = 2.0*asind(wave/(2.*dsp))
     648        print 'HKLD:',H[:4],'2-theta: %.4f'%(tth)
     649        elcent,phi,radii = ellipse = GetEllipse(dsp,data)
     650        data['ellipses'].append(copy.deepcopy(ellipse+('g',)))
     651        print fmt%('predicted ellipse:',elcent[0],elcent[1],phi,radii[0],radii[1])
     652        Ring,delt,sumI = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
    595653        if Ring:
    596             numZ = len(Ring)
    597654            data['rings'].append(np.array(Ring))
    598             newellipse = FitRing(Ring,delt)
    599             elcent,phi,radii = newellipse               
    600             if abs(phi) > 45. and phi < 0.:
    601                 phi += 180.
    602             dist = calcDist(radii,tth)
    603             distR = 1.-dist/data['distance']
    604             if abs(distR) > 0.1:
    605 #                print dsp,dist,data['distance'],distR,len(Ring),delt
    606                 break
    607             if distR > 0.001:
    608                 print 'Wavelength too large?'
    609             elif distR < -0.001:
    610                 print 'Wavelength too small?'
    611             else:
    612                 ellipse = newellipse
    613             zdis,cosB = calcZdisCosB(radius,tth,radii)
    614             Tilt = acosd(cosB)          # 0 <= tilt <= 90
    615             zsinp = zdis*sind(ellipse[1])
    616             zcosp = zdis*cosd(ellipse[1])
    617             cent1.append(np.array([elcent[0]+zsinp,elcent[1]-zcosp]))
    618             cent2.append(np.array([elcent[0]-zsinp,elcent[1]+zcosp]))
     655            rings = np.concatenate((data['rings']),axis=0)
    619656            if i:
    620                 d1 = cent1[-1]-cent1[-2]        #get shift of 2 possible center solutions
    621                 d2 = cent2[-1]-cent2[-2]
    622                 if np.dot(d2,d2) > np.dot(d1,d1):  #right solution is the larger shift
    623                     data['center'] = cent1[-1]
    624                 else:
    625                     data['center'] = cent2[-1]
    626                 Zsum += numZ
    627                 phiSum += numZ*phi
    628                 distSum += numZ*dist
    629                 xSum += numZ*data['center'][0]
    630                 ySum += numZ*data['center'][1]
    631                 tiltSum += numZ*abs(Tilt)
    632                 if not np.all(checkEllipse(Zsum,distSum,xSum,ySum,dist,data['center'][0],data['center'][1])):
    633                     print 'Bad fit for ring # %i. Try reducing Pixel search range'%(i)
    634             cent = data['center']
    635 #            print ('for ring # %2i @ d-space %.4f: dist %.3f rotate %6.2f tilt %6.2f Xcent %.3f Ycent %.3f Npts %d'
    636 #                %(i,dsp,dist,phi-90.,Tilt,cent[0],cent[1],numZ))
     657                chi,chisq = FitDetector(rings,varyList,parmDict,False)
     658                data['distance'] = parmDict['dist']
     659                data['center'] = [parmDict['det-X'],parmDict['det-Y']]
     660                data['rotation'] = np.mod(parmDict['phi'],360.0)
     661                data['tilt'] = parmDict['tilt']
     662                data['DetDepth'] = parmDict['dep']
     663                elcent,phi,radii = ellipse = GetEllipse(dsp,data)
     664                print fmt2%('fitted ellipse:   ',elcent[0],elcent[1],phi,radii[0],radii[1],chisq,len(rings))
    637665            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
    638666        else:
     667            print 'insufficient number of points in this ellipse to fit'
    639668            break
    640669    G2plt.PlotImage(self,newImage=True)
    … …  
    642671    if 2*radii[1] < .9*fullSize:
    643672        print 'Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib'
    644     if not Zsum:
    645         print 'Only one ring fitted. Check your wavelength.'
    646         return False
    647     data['center'] = [xSum/Zsum,ySum/Zsum]
    648     data['distance'] = distSum/Zsum
    649    
    650     #possible error if no. of rings < 3! Might need trap here
    651     d1 = cent1[-1]-cent1[1]             #compare last ring to 2nd ring
    652     d2 = cent2[-1]-cent2[1]
    653     Zsign = 1
    654     len1 = math.sqrt(np.dot(d1,d1))
    655     len2 = math.sqrt(np.dot(d2,d2))
    656     t1 = d1/len1
    657     t2 = d2/len2
    658     if len2 > len1:
    659         if -135. < atan2d(t2[1],t2[0]) < 45.:
    660             Zsign = -1
    661     else:
    662         if -135. < atan2d(t1[1],t1[0]) < 45.:
    663             Zsign = -1
    664    
    665     data['tilt'] = Zsign*tiltSum/Zsum
    666     data['rotation'] = phiSum/Zsum
    667     varyList = ['dist','det-X','det-Y','tilt','phi']
    668     parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
    669         'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
    670     rings = np.concatenate((data['rings']),axis=0)
    671     if data['DetDepthRef']:
    672         varyList.append('dep')
     673    N = len(data['ellipses'])
     674    if N > 2:
    673675        FitDetector(rings,varyList,parmDict)
    674     data['distance'] = parmDict['dist']
    675     data['center'] = [parmDict['det-X'],parmDict['det-Y']]
    676     data['rotation'] = np.mod(parmDict['phi'],360.0)
    677     data['tilt'] = parmDict['tilt']
    678     data['DetDepth'] = parmDict['dep']
    679     N = len(data['ellipses'])
    680     data['ellipses'] = []           #clear away individual ellipse fits
     676        data['distance'] = parmDict['dist']
     677        data['center'] = [parmDict['det-X'],parmDict['det-Y']]
     678        data['rotation'] = np.mod(parmDict['phi'],360.0)
     679        data['tilt'] = parmDict['tilt']
     680        data['DetDepth'] = parmDict['dep']
    681681    for H in HKL[:N]:
    682682        ellipse = GetEllipse(H[3],data)
    … …  
    837837    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
    838838        ellipse = GetEllipse(ring['Dset'],StaControls)
    839         Ring,delt = makeRing(ring['Dset'],ellipse,ring['pixLimit'],ring['cutoff'],scalex,scaley,Image)
     839        Ring,delt,sumI = makeRing(ring['Dset'],ellipse,ring['pixLimit'],ring['cutoff'],scalex,scaley,Image)
    840840        Ring = np.array(Ring).T
    841841        ring['ImxyObs'] = np.array(Ring[:2])      #need to apply masks to this to eliminate bad points

  • trunk/GSASIIimgGUI.py

    r1135 r1148  
    8787        else:
    8888            G2frame.Integrate = G2img.ImageIntegrate(G2frame.ImageZ,data,masks)
    89         G2plt.PlotIntegration(G2frame,newPlot=True)
     89#        G2plt.PlotIntegration(G2frame,newPlot=True)
    9090        G2IO.SaveIntegration(G2frame,G2frame.PickId,data)
    9191        for item in G2frame.MakePDF: item.Enable(True)
    … …  
    203203                for key in keys:
    204204                    if key in ['rotation']:
    205                         File.write(key+':'+str(data[key]-90.)+'\n')                       
     205                        File.write(key+':'+str(data[key])+'\n')                       
    206206                    else:
    207207                        File.write(key+':'+str(data[key])+'\n')
    … …  
    229229                        save[key] = val
    230230                    elif key in ['rotation']:
    231                         save[key] = float(val)+90.
     231                        save[key] = float(val)
    232232                    elif key in ['center',]:
    233233                        if ',' in val:
    … …  
    416416        calibSizer.Add(wx.StaticText(parent=G2frame.dataDisplay,label=' Tilt rotation'),0,
    417417            wx.ALIGN_CENTER_VERTICAL)
    418         rotSel = wx.TextCtrl(parent=G2frame.dataDisplay,value=("%9.3f"%(data['rotation']-90.)),style=wx.TE_READONLY) #kluge to get rotation from vertical - see GSASIIimage
     418        rotSel = wx.TextCtrl(parent=G2frame.dataDisplay,value=("%9.3f"%(data['rotation'])),style=wx.TE_READONLY)
    419419        rotSel.SetBackgroundColour(VERY_LIGHT_GREY)
    420420        calibSizer.Add(rotSel,0,wx.ALIGN_CENTER_VERTICAL)

  • trunk/GSASIIplot.py

    r1134 r1148  
    24942494                xyI = []
    24952495                for azm in Azm:
    2496                     xyI.append(G2img.GetDetectorXY(dspI,azm-90.,Data))
     2496                    xyI.append(G2img.GetDetectorXY(dspI,azm,Data))
    24972497                xyI = np.array(xyI)
    24982498                arcxI,arcyI = xyI.T
    … …  
    25012501                xyO = []
    25022502                for azm in Azm:
    2503                     xyO.append(G2img.GetDetectorXY(dspO,azm-90.,Data))
     2503                    xyO.append(G2img.GetDetectorXY(dspO,azm,Data))
    25042504                xyO = np.array(xyO)
    25052505                arcxO,arcyO = xyO.T
    … …  
    25262526                    N += 1           
    25272527            for ellipse in Data['ellipses']:
     2528#                print 'plot:',ellipse
    25282529                cent,phi,[width,height],col = ellipse
    25292530                Plot.add_artist(Ellipse([cent[0],cent[1]],2*width,2*height,phi,ec=col,fc='none'))
    … …  
    26502651        extent=[ysc[0],ysc[-1],xsc[-1],xsc[0]],aspect='auto')
    26512652    colorBar = Page.figure.colorbar(Img)
    2652     if Data['ellipses']:           
    2653         for ellipse in Data['ellipses']:
    2654             ring = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
    2655             x,y = np.hsplit(ring,2)
    2656             tth,azm = G2img.GetTthAzm(x,y,Data)
    2657 #            azm = np.where(azm < 0.,azm+360,azm)
    2658             Plot.plot(tth,azm,'b,')
     2653#    if Data['ellipses']:           
     2654#        for ellipse in Data['ellipses']:
     2655#            x,y = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
     2656#            tth,azm = G2img.GetTthAzm(x,y,Data)
     2657##            azm = np.where(azm < 0.,azm+360,azm)
     2658#            Plot.plot(tth,azm,'b,')
    26592659    if not newPlot:
    26602660        Page.toolbar.push_current()

  • trunk/fsource/SConstruct

    r771 r1148  
    6262    elif GFORTpath != "":
    6363      FCompiler='gfortran'
     64#      LDFLAGS = '-static-libgfortran -static-libgcc'
    6465    else:
    6566      print 'No Fortran compiler in path'
Note: See TracChangeset for help on using the changeset viewer.