Python MCUtils.angularSeparation方法代码示例

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import angularSeparation [as 别名]
def bg_mask_sources(band,ra0,dec0,ras,decs,responses,sources,maskradius=1.5):
    # At present, masks to 1.5 sigma where FWHM = 2.3548*sigma
    for i in range(len(sources['ra'])):
        ix = np.where(mc.angularSeparation(sources['ra'][i], sources['dec'][i],
                ras,decs)>=(maskradius/2.3548)*np.median(sources['fwhm'][i,:]))
        ras, decs, responses = ras[ix], decs[ix], responses[ix]
    return ras,decs,responses

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import angularSeparation [as 别名]
def read_photons(photonfile,ra0,dec0,tranges,radius,verbose=0,
        colnames=['t','x','y','xa','ya','q','xi','eta','ra','dec','flags']):
    """Read a photon list file and return a python dict() with the expected
    format.
    """
    if verbose:
        print 'Reading photon list file: {f}'.format(f=photonfile)
    data = pd.io.parsers.read_csv(photonfile,names=colnames)
    ra,dec = np.array(data['ra']),np.array(data['dec'])
    angsep = mc.angularSeparation(ra0,dec0,ra,dec)
    ix = np.array([])
    for trange in tranges:
        if verbose:
            print trange
        cut = np.where((angsep<=radius) & (np.isfinite(angsep)))[0]
        ix = np.concatenate((ix,cut),axis=0)
    events = {'t':np.array(data['t'][ix])/tscale,
              'ra':np.array(data['ra'][ix]),
              'dec':np.array(data['dec'][ix]),
              'xi':np.array(data['xi'][ix]),
              'eta':np.array(data['eta'][ix]),
              'x':np.array(data['x'][ix]),
              'y':np.array(data['y'][ix])}
    return events

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import angularSeparation [as 别名]
def quickmag(band, ra0, dec0, tranges, radius, annulus=None, data={},
             stepsz=None, verbose=0, maskdepth=20.0,
             maskradius=1.5,detsize=1.25,coadd=False, photonfile=None):
    if verbose:
        mc.print_inline("Retrieving all of the target events.")
    trange = [np.array(tranges).min(),np.array(tranges).max()]
    try:
        searchradius = annulus[1]
    except TypeError:
        searchradius = radius
    data = pullphotons(band, ra0, dec0, tranges, searchradius,
                       verbose=verbose, photonfile=photonfile)
    if verbose:
        mc.print_inline("Isolating source from background.")
    angSep = mc.angularSeparation(ra0, dec0, data['ra'], data['dec'])
    if verbose:
        mc.print_inline("Binning data according to requested depth.")
    # Multiple ways of defining bins
    if coadd:
        bins = np.array(trange)
    elif stepsz:
        bins = np.append(np.arange(min(trange), max(trange), stepsz),
                                                                max(trange))
    else:
        bins = np.unique(np.array(tranges).flatten())
    # This is equivalent in function to np.digitize(data['t'],bins) except
    # that it's much, much faster. See numpy issue #2656.
    ix = np.searchsorted(bins,data['t'],"right")
    # Initialize histogrammed arrays
    # FIXME: allocate these from a dict of constructors
    lcurve_cols = ['counts', 'sources', 'bg_counts','responses',
                   'detxs', 'detys', 't0_data', 't1_data', 't_mean', 'racent',
                   'deccent']
    lcurve = {'params':gphot_params(band,[ra0,dec0],radius,annulus=annulus,
                                    verbose=verbose,
                                    detsize=detsize,stepsz=stepsz,
                                    trange=trange,maskdepth=maskdepth,
                                    maskradius=maskradius)}
    for col in lcurve_cols:
        lcurve[col] = np.zeros(len(bins)-1)
    # FIXME: Bottleneck. There's probably a way to do this without looping.
    # Don't bother looping through anything with no data.
    lcurve['bg'] = {'simple':np.zeros(len(bins)-1),
                    'cheese':np.zeros(len(bins)-1)}
    if annulus is not None:
        lcurve['bg']['sources'] = bg_sources(band,ra0,dec0,annulus[1],
                                             maskdepth=maskdepth)
        lcurve['bg']['eff_area'] = cheese_bg_area(band,ra0,dec0,annulus,
                                                  lcurve['bg']['sources'])
    else:
        lcurve['bg']['sources'] = None
        lcurve['bg']['eff_area'] = 0.
    if verbose:
        mc.print_inline("Populating histograms.")
    for cnt,i in enumerate(np.unique(ix)):
        # Exclude data outside of the bins in searchsorted.
        if i-1<0 or i==len(bins):
            continue
        if verbose:
            mc.print_inline('Binning {i} of {l}.'.format(
                                                i=cnt,l=len(np.unique(ix))))
        t_ix = np.where(ix==i)
        # TODO: Optionally limit data to specific parts of detector.
        rad_ix = np.where((angSep <= radius) & (ix == i))
        # NOTE: This checks for the dim edge case where you have photons in
        #  the annulus but not in the aperture.
        if not len(rad_ix[0]):
            continue
        lcurve['t0_data'][i-1] = data['t'][rad_ix].min()
        lcurve['t1_data'][i-1] = data['t'][rad_ix].max()
        lcurve['t_mean'][i-1] = data['t'][rad_ix].mean()
        lcurve['counts'][i-1] = len(rad_ix[0])
        lcurve['sources'][i-1] = (1./data['response'][rad_ix]).sum()
        lcurve['responses'][i-1] = data['response'][rad_ix].mean()
        lcurve['detxs'][i-1] = data['col'][rad_ix].mean()
        lcurve['detys'][i-1] = data['row'][rad_ix].mean()
        lcurve['racent'][i-1] = data['ra'][rad_ix].mean()
        lcurve['deccent'][i-1] = data['dec'][rad_ix].mean()
        if annulus is not None:
            ann_ix = np.where((angSep > annulus[0]) &
                              (angSep <= annulus[1]) & (ix == i))
            lcurve['bg_counts'][i-1] = len(ann_ix[0])
            # Background is reported as counts within the aperture
            lcurve['bg']['simple'][i-1] = (mc.area(radius) *
                (1./data['response'][ann_ix]).sum() /
                (mc.area(annulus[1])-mc.area(annulus[0])))
            lcurve['bg']['cheese'][i-1] = cheese_bg(band, ra0, dec0, radius,
                annulus, data['ra'][t_ix], data['dec'][t_ix],
                data['response'][t_ix], maskdepth=maskdepth,
                eff_area=lcurve['bg']['eff_area'],
                sources=lcurve['bg']['sources'])
        else:
            lcurve['bg_counts'][i-1]=0.
            lcurve['bg']['simple'][i-1]=0.
            lcurve['bg']['cheese'][i-1]=0.
    # Only return bins that contain data.
    ix = np.where((np.isfinite(lcurve['sources'])) &
                  (np.array(lcurve['sources']) > 0))
    lcurve['t0'] = bins[ix]
    lcurve['t1'] = bins[ix[0]+1]
#.........这里部分代码省略.........

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import angularSeparation [as 别名]
def bg_mask_annulus(band,ra0,dec0,annulus,ras,decs,responses):
    ix = np.where((mc.angularSeparation(ra0,dec0,ras,decs)>=annulus[0]) &
                  (mc.angularSeparation(ra0,dec0,ras,decs)<=annulus[1]))
    return ras[ix],decs[ix],responses[ix]

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import angularSeparation [as 别名]
                                                                 markersize=5)
    # dRA
    plt.subplot(2,2,2,yticks=[],xticks=[],ylim=[-0.004*a,0.004*a])
    plt.hist(dRA*np.cos(data[band]['ra'])*a,bins=500,orientation='horizontal',
                                                                    color='k')
    # dDec
    plt.subplot(2,2,3,yticks=[],xlim=[-0.004*a,0.004*a])
    plt.xlabel('{d}Dec (arcsec)'.format(d=r'$\Delta$'))
    plt.gca().invert_yaxis()
    plt.hist(dDec*a,bins=500,color='k')
    fig.savefig('../calpaper/src/dRA_v_dDec({band})'.format(band=band))

fig = plt.figure(figsize=(8*scl,4*scl))
fig.subplots_adjust(left=0.12,right=0.95,wspace=0.02,bottom=0.15,top=0.9)
for i,band in enumerate(bands):
    delta = mc.angularSeparation(data[band]['ra'],data[band]['dec'],
                                 data[band]['racent'],data[band]['deccent'])
    plt.subplot(1,2,i+1,yticks=[],xlim=[0.*a,0.002*a])
    plt.title('{band} Angular Separation (arcsec)'.format(
                                                    band=band,d=r'$\Delta$'))
    plt.hist(delta*a,bins=500,range=[0.*a,0.002*a],color='k')
    fig.savefig('../calpaper/src/angSep({band}).png'.format(band=band))

###############################################################################
"""Deadtime Sanity Checks
According to the calibration paper, the FUV deadtime correction should be
small (~ a few percent), but it is actually bigger than the NUV correction.
"""
fig = plt.figure(figsize=(8,4))
fig.subplots_adjust(left=0.12,right=0.95,wspace=0.02,bottom=0.15,top=0.9)
for i,band in enumerate(bands):
    plt.subplot(1,2,i+1,yticks=[])

# 需要导入模块: import MCUtils [as 别名]
# 或者: from MCUtils import angularSeparation [as 别名]
 def test_angularSeparation(self):
         self.assertAlmostEqual(
             mc.angularSeparation(10,20,10.1,20.1),0.13720278279273748)