Changeset 4104 for trunk/GSASIIimage.py


Ignore:
Timestamp:
Aug 20, 2019 9:06:16 PM (2 years ago)
Author:
toby
Message:

new mode for combined fit of wavelength to a set of images

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIimage.py

    r4102 r4104  
    138138
    139139def FitDetector(rings,varyList,parmDict,Print=True,covar=False):
    140     '''Fit detector calibration paremeters
     140    '''Fit detector calibration parameters
    141141
    142142    :param np.array rings: vector of ring positions
     
    225225        return [chisq,vals,sigList]
    226226
     227def FitMultiDist(rings,varyList,parmDict,Print=True,covar=False):
     228    '''Fit detector calibration parameters with multi-distance data
     229
     230    :param np.array rings: vector of ring positions (x,y,dist,d-space)
     231    :param list varyList: calibration parameters to be refined
     232    :param dict parmDict: calibration parameters
     233    :param bool Print: set to True (default) to print the results
     234    :param bool covar: set to True to return the covariance matrix (default is False)
     235    :returns: [chisq,vals,sigDict] unless covar is True, then
     236        [chisq,vals,sigDict,coVarMatrix] is returned
     237    '''
     238       
     239    def CalibPrint(parmDict,sigDict,chisq,Npts):
     240        print ('Image Parameters: chi**2: %12.3g, Np: %d'%(chisq,Npts))
     241        ptlbls = 'names :'
     242        ptstr =  'values:'
     243        sigstr = 'esds  :'
     244        names = ['wavelength', 'dep', 'phi', 'tilt']
     245        if 'deltaDist' in parmDict:
     246            names += ['deltaDist']
     247        for name in names:
     248            if name == 'wavelength':
     249                fmt = '%12.6f'
     250            elif name == 'dep':
     251                fmt = '%12.4f'
     252            else:
     253                fmt = '%12.3f'
     254               
     255            ptlbls += "%s" % (name.rjust(12))
     256            if name == 'phi':
     257                ptstr += fmt % (parmDict[name]%360.)
     258            else:
     259                ptstr += fmt % (parmDict[name])
     260            if name in sigDict:
     261                sigstr += fmt % (sigDict[name])
     262            else:
     263                sigstr += 12*' '
     264        print (ptlbls)
     265        print (ptstr)
     266        print (sigstr)
     267        ptlbls = 'names :'
     268        ptstr =  'values:'
     269        sigstr = 'esds  :'
     270        for d in sorted(set([i[5:] for i in parmDict.keys() if 'det-X' in i]),key=lambda x:int(x)):
     271            fmt = '%12.3f'
     272            for key in 'det-X','det-Y','delta':
     273                name = key+d
     274                if name not in parmDict: continue
     275                ptlbls += "%12s" % name
     276                ptstr += fmt % (parmDict[name])
     277                if name in sigDict:
     278                    sigstr += fmt % (sigDict[name])
     279                else:
     280                    sigstr += 12*' '
     281                if len(ptlbls) > 68:
     282                    print()
     283                    print (ptlbls)
     284                    print (ptstr)
     285                    print (sigstr)
     286                    ptlbls = 'names :'
     287                    ptstr =  'values:'
     288                    sigstr = 'esds  :'
     289        if len(ptlbls) > 0:
     290            print()
     291            print (ptlbls)
     292            print (ptstr)
     293            print (sigstr)
     294
     295    def ellipseCalcD(B,xyd,varyList,parmDict):
     296        x,y,dist,dsp = xyd
     297        varyDict = dict(zip(varyList,B))
     298        parms = {}
     299        for parm in parmDict:
     300            if parm in varyList:
     301                parms[parm] = varyDict[parm]
     302            else:
     303                parms[parm] = parmDict[parm]
     304        # create arrays with detector center values
     305        detX = np.array([parms['det-X'+str(int(d))] for d in dist])
     306        detY = np.array([parms['det-Y'+str(int(d))] for d in dist])
     307        if 'deltaDist' in parms:
     308            deltaDist = parms['deltaDist']
     309        else:
     310            deltaDist = np.array([parms['delta'+str(int(d))] for d in dist])
     311               
     312        phi = parms['phi']-90.               #get rotation of major axis from tilt axis
     313        tth = 2.0*npasind(parms['wavelength']/(2.*dsp))
     314        phi0 = npatan2d(y-detY,x-detX)
     315        dxy = peneCorr(tth,parms['dep'],dist-deltaDist,parms['tilt'],phi0)
     316        stth = npsind(tth)
     317        cosb = npcosd(parms['tilt'])
     318        tanb = nptand(parms['tilt'])       
     319        tbm = nptand((tth-parms['tilt'])/2.)
     320        tbp = nptand((tth+parms['tilt'])/2.)
     321        d = (dist-deltaDist)+dxy
     322        fplus = d*tanb*stth/(cosb+stth)
     323        fminus = d*tanb*stth/(cosb-stth)
     324        vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
     325        vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
     326        R0 = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2.      #+minor axis
     327        R1 = (vplus+vminus)/2.                                    #major axis
     328        zdis = (fplus-fminus)/2.
     329        Robs = np.sqrt((x-detX)**2+(y-detY)**2)
     330        rsqplus = R0**2+R1**2
     331        rsqminus = R0**2-R1**2
     332        R = rsqminus*npcosd(2.*phi0-2.*phi)+rsqplus
     333        Q = np.sqrt(2.)*R0*R1*np.sqrt(R-2.*zdis**2*npsind(phi0-phi)**2)
     334        P = 2.*R0**2*zdis*npcosd(phi0-phi)
     335        Rcalc = (P+Q)/R
     336        return (Robs-Rcalc)*10.        #why 10? does make "chi**2" more reasonable
     337       
     338    p0 = [parmDict[key] for key in varyList]
     339    result = leastsq(ellipseCalcD,p0,args=(rings.T,varyList,parmDict),full_output=True,ftol=1.e-8)
     340    chisq = np.sum(result[2]['fvec']**2)/(rings.shape[0]-len(p0))   #reduced chi^2 = M/(Nobs-Nvar)
     341    parmDict.update(zip(varyList,result[0]))
     342    vals = list(result[0])
     343    if chisq > 1:
     344        sig = list(np.sqrt(chisq*np.diag(result[1])))
     345    else:
     346        sig = list(np.sqrt(np.diag(result[1])))
     347    sigDict = {name:s for name,s in zip(varyList,sig)}
     348    if Print:
     349        CalibPrint(parmDict,sigDict,chisq,rings.shape[0])
     350    if covar:
     351        return [chisq,vals,sigDict,result[1]]
     352    else:
     353        return [chisq,vals,sigDict]
     354   
    227355def ImageLocalMax(image,w,Xpix,Ypix):
    228356    'Needs a doc string'
     
    556684    return tam.T
    557685   
    558 def ImageRecalibrate(G2frame,ImageZ,data,masks):
     686def ImageRecalibrate(G2frame,ImageZ,data,masks,getRingsOnly=False):
    559687    '''Called to repeat the calibration on an image, usually called after
    560688    calibration is done initially to improve the fit.
     
    564692    :param dict data: the Controls dict for the image
    565693    :param dict masks: a dict with masks
    566     :returns: a list containing vals,varyList,sigList,parmDict,covar
     694    :returns: a list containing vals,varyList,sigList,parmDict,covar or rings
     695      (with an array of x, y, and d-space values) if getRingsOnly is True
     696      or an empty list, in case of an error
    567697    '''
    568698    import ImageCalibrants as calFile
    569     G2fil.G2Print ('Image recalibration:')
     699    if not getRingsOnly:
     700        G2fil.G2Print ('Image recalibration:')
    570701    time0 = time.time()
    571702    pixelSize = data['pixelSize']
     
    640771       
    641772    rings = np.concatenate((data['rings']),axis=0)
     773    if getRingsOnly:
     774        return rings
    642775    [chisq,vals,sigList,covar] = FitDetector(rings,varyList,parmDict,True,True)
    643776    data['wavelength'] = parmDict['wave']
Note: See TracChangeset for help on using the changeset viewer.