Python interpolate.interp2d函数代码示例

def applyInterpolatedNoise(image, targets, noise_factor=10, stride=64):

    size = image.shape[0]

    x = np.arange(0, size, stride)
    y = np.arange(0, size, stride)

    interp_size = x.shape[0]

    delta_x = np.random.normal(scale=noise_factor, size=(interp_size, interp_size))
    delta_y = np.random.normal(scale=noise_factor, size=(interp_size, interp_size))

    f_x = interpolate.interp2d(x, y, delta_x, kind='cubic')
    f_y = interpolate.interp2d(x, y, delta_y, kind='cubic')

    noise_x = f_x(np.arange(0, size), np.arange(0, size))
    noise_y = f_y(np.arange(0, size), np.arange(0, size))

    grid_x = np.asarray([range(size)]*size)
    grid_y = np.asarray([size*[i] for i in range(size)])

    x_jitter = grid_x + noise_x
    y_jitter = grid_y + noise_y

    labels = [None for t in range(targets.shape[0])]
    for t in range(targets.shape[0]):
        labels[t] = lookupNearest(targets[t, :, :], x_jitter, y_jitter)

    return bilinear_interpolate(image, x_jitter, y_jitter), labels

def cdcl(v, nrpm):
    """Gives the value of the drag and lift coefficients as function of velocity and spin.
    Is only valid for velocities between 13.7 and 88.1 m/s, and spins between 2000 and 6000 rpm.
    The value is determined with linear interpolation of the data given by Bearman and Harvey, 
    Golf Ball Aerodynamics, volume 27, Aeronautival Quarterly, 1976."""
    
    from numpy import array, linspace
    from scipy.interpolate import interp2d
    import sys
    
    if (v<13.7 or v>88.1):
        sys.exit('v is out of bounds. Must be between 13.7 and 88.41 m/s.')
        return
    
    if (nrpm<2000 or nrpm>6000):
        sys.exit('nrpm is out of bounds. Must be between 2000 and 6000 rpm.')
        return
    
    v0 = array([13.7, 21.6, 29.9, 38.4, 46.9, 55.2, 63.1, 71.9, 80.2, 88.1])
    nrpm0 = linspace(2000,6000,21)
    
    CD = array([[0.3624,0.2885,0.2765,0.2529,0.2472,0.2481,0.2467,0.2470,0.2470,0.2470],[0.3806,0.3102,0.2853,0.2590,0.2507,0.2498,0.2485,0.2486,0.2484,0.2484],[0.3954,0.3288,0.2937,0.2649,0.2543,0.2516,0.2504,0.2502,0.2497,0.2497],[0.4070,0.3443,0.3018,0.2708,0.2580,0.2535,0.2522,0.2518,0.2511,0.2511],[0.4153,0.3566,0.3095,0.2765,0.2617,0.2556,0.2541,0.2534,0.2524,0.2524],[0.4203,0.3658,0.3169,0.2822,0.2655,0.2578,0.2560,0.2550,0.2538,0.2538],[0.4120,0.3719,0.3240,0.2876,0.2693,0.2602,0.2579,0.2566,0.2551,0.2551],[0.3960,0.3749,0.3308,0.2930,0.2732,0.2627,0.2599,0.2582,0.2565,0.2565],[0.3876,0.3766,0.3372,0.2983,0.2772,0.2653,0.2619,0.2598,0.2578,0.2578],[0.4100,0.3854,0.3433,0.3034,0.2811,0.2681,0.2639,0.2614,0.2592,0.2592],[0.4288,0.4003,0.3490,0.3084,0.2852,0.2710,0.2659,0.2630,0.2605,0.2605],[0.4445,0.4082,0.3544,0.3133,0.2893,0.2741,0.2680,0.2646,0.2619,0.2619],[0.4575,0.4153,0.3595,0.3180,0.2934,0.2772,0.2701,0.2662,0.2632,0.2632],[0.4682,0.4215,0.3643,0.3227,0.2976,0.2806,0.2722,0.2678,0.2646,0.2646],[0.4772,0.4269,0.3687,0.3272,0.3019,0.2840,0.2743,0.2694,0.2659,0.2659],[0.4848,0.4314,0.3728,0.3316,0.3062,0.2876,0.2765,0.2710,0.2673,0.2673],[0.4914,0.4350,0.3765,0.3358,0.3105,0.2913,0.2787,0.2726,0.2686,0.2686],[0.4976,0.4377,0.3799,0.3400,0.3149,0.2952,0.2809,0.2742,0.2700,0.2700],[0.5039,0.4395,0.3830,0.3440,0.3194,0.2992,0.2831,0.2758,0.2713,0.2713],[0.5105,0.4405,0.3858,0.3479,0.3239,0.3034,0.2854,0.2774,0.2727,0.2727],[0.5180,0.4406,0.3882,0.3517,0.3285,0.3076,0.2877,0.2790,0.2740,0.2740]])
    CL = array([[0.1040,0.1846,0.2460,0.1984,0.1762,0.1538,0.1418,0.1360,0.1280,0.1276],[0.1936,0.2318,0.2590,0.2089,0.1824,0.1603,0.1476,0.1405,0.1321,0.1324],[0.2608,0.2694,0.2715,0.2191,0.1886,0.1668,0.1533,0.1450,0.1362,0.1367],[0.3090,0.2986,0.2835,0.2289,0.1949,0.1731,0.1589,0.1494,0.1403,0.1404],[0.3418,0.3205,0.2947,0.2384,0.2011,0.1794,0.1646,0.1539,0.1444,0.1436],[0.3624,0.3362,0.3049,0.2475,0.2073,0.1856,0.1702,0.1584,0.1485,0.1464],[0.3743,0.3470,0.3140,0.2562,0.2135,0.1916,0.1757,0.1629,0.1526,0.1488],[0.3808,0.3541,0.3217,0.2644,0.2197,0.1976,0.1813,0.1673,0.1567,0.1508],[0.3854,0.3584,0.3280,0.2722,0.2259,0.2035,0.1868,0.1718,0.1608,0.1524],[0.3915,0.3614,0.3325,0.2795,0.2322,0.2092,0.1923,0.1763,0.1649,0.1539],[0.4005,0.3640,0.3347,0.2862,0.2384,0.2149,0.1977,0.1807,0.1690,0.1582],[0.4010,0.3696,0.3380,0.2925,0.2446,0.2205,0.2032,0.1852,0.1731,0.1624],[0.4026,0.3748,0.3412,0.2982,0.2508,0.2259,0.2086,0.1897,0.1772,0.1666],[0.4050,0.3797,0.3440,0.3033,0.2570,0.2313,0.2139,0.1941,0.1813,0.1707],[0.4084,0.3845,0.3468,0.3078,0.2632,0.2366,0.2193,0.1986,0.1854,0.1749],[0.4130,0.3894,0.3496,0.3119,0.2695,0.2418,0.2246,0.2031,0.1895,0.1791],[0.4192,0.3947,0.3527,0.3149,0.2757,0.2469,0.2298,0.2076,0.1936,0.1833],[0.4271,0.4004,0.3563,0.3184,0.2819,0.2518,0.2351,0.2120,0.1977,0.1875],[0.4371,0.4069,0.3607,0.3233,0.2881,0.2567,0.2403,0.2165,0.2018,0.1916],[0.4494,0.4143,0.3660,0.3300,0.2943,0.2615,0.2455,0.2210,0.2059,0.1958],[0.4644,0.4227,0.3724,0.3393,0.3005,0.2662,0.2507,0.2254,0.2100,0.2000]])
    
    cd = interp2d(v0, nrpm0, CD, kind='linear')
    cl = interp2d(v0, nrpm0, CL, kind='linear')
    return cd(v, nrpm), cl(v, nrpm)

 def __init__(self, filename):
     nc = Dataset(filename)
     self.debug_mode = 1.
     print self.debug_mode
     self.dt = nc.variables["time"][1]-nc.variables["time"][0]
     self.xlon = nc.variables["xlon"][:]
     self.xlat = nc.variables["xlat"][:]
     N_lon = self.xlon.shape[1]
     N_lat = self.xlat.shape[0]
     #self.idx_to_lon = interp1d(numpy.arange(N_lon), self.xlon[0,:])
     #self.idx_to_lat = interp1d(numpy.arange(N_lat), self.xlat[:,0])
     self.idx_to_lon = interp2d(numpy.arange(N_lon), numpy.arange(N_lat), self.xlon)
     self.idx_to_lat = interp2d(numpy.arange(N_lon), numpy.arange(N_lat), self.xlat)
     self.left = self.xlon.min()
     self.right = self.xlon.max()
     self.bottom = self.xlat.min()+18
     self.top = self.xlat.max()-8
     self.m = Basemap(llcrnrlon=self.left,
                      llcrnrlat=self.bottom,
                      urcrnrlon=self.right,
                      urcrnrlat=self.top,
                      projection='cyl', resolution='l')
     self.ocean_mask = maskoceans(self.xlon, self.xlat, self.xlon).mask
     self.time_min = 0.
     self.time_max = 0.
     self.no_ingested = 0
     self.masks = numpy.empty((0,)+self.ocean_mask.shape)
     self.tc_table = numpy.empty((0,5))
     self.time = numpy.empty((0,))

    def __init__(self, textureData):
        cols = len(textureData[0])
        rows = len(textureData)
        self.cols = cols
        self.rows = rows

        indicators = np.empty_like(textureData)

        for i in range(rows):
            for j in range(cols):
                if textureData[i][j] == -1:
                    indicators[i][j] = 1.0
                else:
                    indicators[i][j] = 0.0

        marginX = 1.0/(2.0 * cols)
        marginY = 1.0/(2.0 * rows)

        x = np.linspace(marginX, 1.0-marginX, cols)
        y = np.linspace(marginY, 1.0-marginY, rows)

        self.indicator = interpolate.interp2d(x, y, indicators, kind='linear')

        data = textureData
        data = np.vstack((data[0], data))
        data = np.vstack((data, data[-1]))
        data = np.column_stack((data[:,0],data))
        data = np.column_stack((data, data[:,-1]))
        self.textureData = data

        x = np.linspace(-marginX, 1.0+marginX, cols+2)
        y = np.linspace(-marginY, 1.0+marginY, rows+2)

        self.f = interpolate.interp2d(x, y, data, kind='linear')

def driftcorr(A, ux, uy, interpolation='cubic'):
    ''' 
    drift correction on 2D image or 3D DOS map, use ux, uy calculated by driftmap. Crop edges of A_corr if needed.
    interpolation: 'linear', 'cubic' or 'quintic'. Default is 'cubic'
    '''

    A_corr = np.zeros_like(A)
    s = A.shape[-1]
    t = np.arange(s, dtype='float')
    x, y = np.meshgrid(t, t)
    xnew = (x - ux).ravel()
    ynew = (y - uy).ravel()
    tmp = np.zeros(s**2)

    if len(A.shape) is 2:
        tmp_f = interp2d(t, t, A, kind=interpolation)
        for ix in range(tmp.size):
            tmp[ix] = tmp_f(xnew[ix], ynew[ix])
        A_corr = tmp.reshape(s, s)
        return A_corr
    elif len(A.shape) is 3:
        for iz, layer in enumerate(A):
            tmp_f = interp2d(t, t, layer, kind=interpolation)
            for ix in range(tmp.size):
                tmp[ix] = tmp_f(xnew[ix], ynew[ix])
            A_corr[iz] = tmp.reshape(s, s)
            print('Processing slice %d/%d...'%(iz+1, A.shape[0]), end='\r')
        return A_corr
    else:
        print('ERR: Input must be 2D or 3D numpy array!')

    def rebin_data(self, grid, use_psf=True):
        """Calculates the center of mass of the grid and then
        rebins so that the center pixel really is the center of the array
        For this we do a 2-d interpolation on the grid
        """
        a = psf_fitter.psffit(abs(grid), circle=False, rotate=1)
        xcen = a[2]
        ycen = a[2]
        xlen, ylen = grid.shape
        xval = arange(xlen)
        yval = arange(ylen)

        xint = interp1d(xval, self.xpos_abs)
        yint = interp1d(yval, self.ypos_abs)

        xintcen = self.xmax_pos-xint(xcen)
        yintcen = self.ymax_pos-yint(ycen)

        print self.xmax_pos, xintcen, self.ymax_pos, yintcen
        f_real = interp2d(self.xpos_rel, self.ypos_rel, real(grid))
        f_imag = interp2d(self.xpos_rel, self.ypos_rel, imag(grid))

        xnew = self.xpos_rel - xintcen
        ynew = self.ypos_rel - yintcen

        recen_grid = f_real(xnew, ynew) + 1j*f_imag(xnew, ynew)

        print nd.center_of_mass(abs(recen_grid))

        return recen_grid

def get_translate(workdir=None):

    filename = os.path.join(workdir, "Vicalloy/Fe-Co-V_140922a_META_DATA.csv")
    compdata_f = pd.read_csv(filename, sep='\t').dropna()
    print compdata_f.head()
    x = compdata_f["Xnom (mm)"].values
    y = compdata_f["Ynom (mm)"].values
    Co_concentration = compdata_f["Co (at%)"].values
    Fe_concentration = compdata_f["Fe (at%)"].values
    V_concentration  = compdata_f["V (at%)"].values
    method = 'linear'
    # method = 'nearest'

    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        Co_concI = interp2d(x,y,Co_concentration, kind = method)
        Fe_concI = interp2d(x,y,Fe_concentration, kind = method)
        V_concI = interp2d(x,y,V_concentration , kind = method)

    def translate(key):
        manip_z, manip_y = key
        sample_y = manip_z - 69.5
        sample_x = (manip_y +8) *2
        Co = Co_concI(sample_x,sample_y)[0]/100.
        Fe = Fe_concI(sample_x,sample_y)[0]/100.
        V  = V_concI(sample_x,sample_y)[0]/100.
        return (
            "Fe{:.2f}Co{:.2f}V{:.2f}".format(Fe,Co,V),
            sample_x, sample_y
        )


    return translate

def ijcoast(coast, grd):

    if type(grd).__name__ == "ROMS_Grid":
        lon = grd.hgrid.lon_vert
        lat = grd.hgrid.lat_vert
    if type(grd).__name__ == "CGrid_geo":
        lon = grd.lon_vert
        lat = grd.lat_vert

    ijcoast = []

    for k in range(coast.shape[0]):
        if np.isnan(coast[k, 0]):
            ijcoast.append([np.nan, np.nan])
        else:
            iindex, jindex = find_nearestgridpoints(coast[k, 0], coast[k, 1], grd, Cpos="vert")
            if iindex:
                i, j = np.meshgrid(iindex, jindex)
                x = lon[j, i]
                y = lat[j, i]
                funct_i = interpolate.interp2d(x.flatten(), y.flatten(), i.flatten())
                funct_j = interpolate.interp2d(x.flatten(), y.flatten(), j.flatten())
                i_coast = funct_i(coast[k, 0], coast[k, 1])[0]
                j_coast = funct_j(coast[k, 0], coast[k, 1])[0]
                ijcoast.append([i_coast, j_coast])

    return np.array(ijcoast)

def interpolate_isentropic_T_rho_Vp_Vs(T_ref, filename):
    pressures, temperatures, property_array = table_from_tab(filename)
    s, rho, vp, vs = property_array
    entropy = interpolate.interp2d(pressures, temperatures, s, kind='linear')
    density = interpolate.interp2d(pressures, temperatures, rho, kind='linear')
    v_p = interpolate.interp2d(pressures, temperatures, vp, kind='linear')
    v_s = interpolate.interp2d(pressures, temperatures, vs, kind='linear')
    S_ref = entropy(pressures.min(), T_ref)
    if isnan(S_ref):
        print('Oh no! For some reason the entropy interpolation on your perplex file returns NaN at the minimum pressure.')
        print('Conditions: P='+str(pressures.min()/1.e5)+' bar, T='+str(T_ref)+' K.')
        exit()
    adiabatic_temperatures = np.empty_like(pressures)
    adiabatic_densities = np.empty_like(pressures)
    adiabatic_Vps = np.empty_like(pressures)
    adiabatic_Vss = np.empty_like(pressures)
    for i, s_p in enumerate(zip(*s)):
        entropy_slice = interpolate.interp1d(s_p, temperatures, kind='linear')
        adiabatic_temperatures[i] = entropy_slice(S_ref)
        adiabatic_densities[i] = density(pressures[i], adiabatic_temperatures[i])
        adiabatic_Vps[i] = v_p(pressures[i], adiabatic_temperatures[i])
        adiabatic_Vss[i] = v_s(pressures[i], adiabatic_temperatures[i])
    
    adiabat_density = interpolate.interp1d(pressures, adiabatic_densities, kind='linear')
    adiabat_Vp = interpolate.interp1d(pressures, adiabatic_Vps, kind='linear')
    adiabat_Vs = interpolate.interp1d(pressures, adiabatic_Vss, kind='linear')
    adiabat_temperature = interpolate.interp1d(pressures, adiabatic_temperatures, kind='linear')

    return adiabat_temperature, adiabat_density, adiabat_Vp, adiabat_Vs

def refine_data(lon, lat, f, refine):
    lon = lon[0, :]
    lat = lat[:, 0]

    dlon = lon[1] - lon[0]
    dlat = lat[1] - lat[0]

    lat_hi = np.arange(lat[0],lat[-1],dlat/refine)
    lon_hi = np.arange(lon[0],lon[-1],dlon/refine)

    nx = len(lon_hi)
    ny = len(lat_hi)

    a = np.array(f.mask).astype(int)

    f[np.isnan(f)] = 100000

    ipol = interp2d(lon, lat, f)
    apol = interp2d(lon, lat, a)
    f = ipol(lon_hi, lat_hi)
    a = apol(lon_hi, lat_hi)
    f = np.ma.masked_where(a>.2, f)


    lon_hi, lat_hi = np.meshgrid(lon_hi, lat_hi)
    return lon_hi, lat_hi, f

 def __init__(self,folder,T_unit='C',S_unit = 'kg/kg'):
     self.name = 'Liu'
     self.folder = folder
     self.T_unit = T_unit
     self.S_unit = S_unit
     #Simple dictionary for converting units of temperature and salinity
     #Pressure will always be in MPa
     self.S_multiply = {'kg/kg': 1.0, 'ppt': 0.001, 'wt. %': 0.01}#Liu calculation is in kg/kg
     self.T_add = {'K': -273.15, 'C': 0.0, 'F': 32}#Liu calculation is in C
     self.T_multiply = {'K': 1.0, 'C': 1.0, 'F': 5.0/9.0}#Liu calculation is in C
     #Load calculated dissociation pressure from file
     P_file = self.folder + 'Liu_method_Pmat.csv'
     A = np.loadtxt(P_file,delimiter=',')
     Temp = A[1:,0]#Celsius
     Sal = A[0,1:]#kg/kg
     Peq = A[1:,1:]#MPa
     self.P_func = interpolate.interp2d(Sal,Temp,Peq,kind='cubic')
     #Load calculated dissociation salinity from file
     S_file = self.folder + 'Liu_method_Smat.csv'
     B = np.loadtxt(S_file,delimiter=',')
     Pres = B[1:,0]#Celsius
     Temp2 = B[0,1:]#MPa
     Seq = B[1:,1:]*self.S_multiply[self.S_unit]#kg/kg
     Seq[Seq<1e-3]=0
     self.S_func = interpolate.interp2d(Temp2,Pres,Seq,kind='cubic')

def default_absorbers(Tatm, ozone_file = 'apeozone_cam3_5_54.nc'):
    '''Initialize a dictionary of well-mixed radiatively active gases
    All values are volumetric mixing ratios.

    Ozone is set to a climatology.

    All other gases are assumed well-mixed:

        - CO2
        - CH4
        - N2O
        - O2
        - CFC11
        - CFC12
        - CFC22
        - CCL4

    Specific values are based on the AquaPlanet Experiment protocols,
    except for O2 which is set the realistic value 0.21
    (affects the RRTMG scheme).
    '''
    absorber_vmr = {}
    absorber_vmr['CO2']   = 348. / 1E6
    absorber_vmr['CH4']   = 1650. / 1E9
    absorber_vmr['N2O']   = 306. / 1E9
    absorber_vmr['O2']    = 0.21
    absorber_vmr['CFC11'] = 0.
    absorber_vmr['CFC12'] = 0.
    absorber_vmr['CFC22'] = 0.
    absorber_vmr['CCL4']  = 0.

    datadir = os.path.join(os.path.dirname(__file__), 'data', 'ozone')
    ozonefilepath = os.path.join(datadir, ozone_file)
    #  Open the ozone data file
    print 'Getting ozone data from', ozonefilepath
    ozonedata = nc.Dataset(ozonefilepath)
    ozone_lev = ozonedata.variables['lev'][:]
    ozone_lat = ozonedata.variables['lat'][:]
    #  zonal and time average
    ozone_zon = np.mean(ozonedata.variables['OZONE'], axis=(0,3))
    ozone_global = np.average(ozone_zon, weights=np.cos(np.deg2rad(ozone_lat)), axis=1)
    lev = Tatm.domain.axes['lev'].points
    if Tatm.shape == lev.shape:
        # 1D interpolation on pressure levels using global average data
        f = interp1d(ozone_lev, ozone_global)
        #  interpolate data to model levels
        absorber_vmr['O3'] = f(lev)
    else:
        #  Attempt 2D interpolation in pressure and latitude
        f2d = interp2d(ozone_lat, ozone_lev, ozone_zon)
        try:
            lat = Tatm.domain.axes['lat'].points
            f2d = interp2d(ozone_lat, ozone_lev, ozone_zon)
            absorber_vmr['O3'] = f2d(lat, lev).transpose()
        except:
            print 'Interpolation of ozone data failed.'
            print 'Reverting to default O3.'
            absorber_vmr['O3'] = np.zeros_like(Tatm)
    return absorber_vmr

def streamline(u1,u2,d,x1_0,x2_0,dl=1.,fig=None,nmax=600):

    if fig==None:
        ax=plt.gca()
    else:
        ax=fig.figure.gca()

    xrange=[ax.get_xlim()[0],ax.get_xlim()[1]]
    yrange=[ax.get_ylim()[0],ax.get_ylim()[1]]

    if nmax == 600 and fig != None:
        nmax = fig.dpi * max([fig.fig_w,fig.fig_h])


    # Aspect ratio:
    r=(xrange[1]-xrange[0])/(yrange[1]-yrange[0])
    if r<1:
        ny=nmax
        nx=int(r*nmax)
    else:
        nx=nmax
        ny=int(nmax/r)
    nregrid = [nx,ny]

    CC=d.getCenterPoints()
    tmp0=np.complex(0,nregrid[0])
    tmp1=np.complex(0,nregrid[1])
    x=np.linspace(xrange[0],xrange[1],nregrid[0])
    y=np.linspace(yrange[0],yrange[1],nregrid[1])
    grid_x, grid_y = np.mgrid[xrange[0]:xrange[1]:tmp0, yrange[0]:yrange[1]:tmp1]

    u = griddata(CC, u1, (grid_x, grid_y), method='linear')
    v = griddata(CC, u2, (grid_x, grid_y), method='linear')
    uisnan=np.isnan(u)
    visnan=np.isnan(v)
    un = np.empty(np.shape(u))
    vn = np.empty(np.shape(v))
    un[uisnan] = griddata(CC, u1, (grid_x[uisnan], grid_y[uisnan]), method='nearest')
    vn[visnan] = griddata(CC, u2, (grid_x[visnan], grid_y[visnan]), method='nearest')
    u[uisnan]=un[uisnan]
    v[visnan]=vn[visnan]
# Normalize:
    mag=np.sqrt(u**2+v**2)
    u=u/mag
    v=v/mag

    fu = interpolate.interp2d(grid_x, grid_y, u, kind='cubic')
    fv = interpolate.interp2d(grid_x, grid_y, v, kind='cubic')
    x1=[x1_0]
    x2=[x2_0]
    while(x1[-1] >= xrange[0] and x1[-1] <= xrange[1] and
    x2[-1] >= yrange[0] and x2[-1] <= yrange[1]):
        dt=dl
        x1.append(x1[-1]+fu(x1[-1],x2[-1])*dt)
        x2.append(x2[-1]+fv(x1[-1],x2[-1])*dt)
    return [np.array(x1),np.array(x2)]

def bad_make_star(mass, age, model='Burrows97'):
    
    '''
    Calculates stellar properties such as Teff, radii, logg and logL using evolutionary models
    '''

    if np.min(mass) < 0.0005:
        raise NameError('Mass below minimum mass of 0.0005Msun')

    if np.max(mass) > 0.1 and model=='Baraffe03':
       warnings.warn('Mass above maximum mass of 0.1Msun for Baraffe 2003. Using Burrows 1997 instead.')  
       model = 'Burrows97'

    if np.min(mass) > 0.2:
        raise NameError('Mass above maximum mass of 0.2Msun for Burrows 1997')    

    if model == 'Burrows97':
        #0.0005 - 0.2 Msun
        burrows = pd.read_pickle("burrows97.pickle")
        allages = burrows["Age (Gyr)"]
        allmasses = burrows["M/Ms"]
        teff = burrows["Teff"]
        radius = burrows["R/Rs"]
        logg = burrows["logg(cgs)"]
        logL = burrows["logL/Ls"]
        
    if model == 'Baraffe03':
        #0.0005 - 0.1 Msun
        baraffe = pd.read_pickle("baraffe03.pickle")
        allages = baraffe["Age (Gyr)"]
        allmasses = baraffe["M/Ms"]
        teff = baraffe["Teff"]
        radius = baraffe["R/Rs"]
        logg = baraffe["logg(cgs)"]
        logL = baraffe["logL/Ls"]

    interpteff = interpolate.interp2d(allages,allmasses,teff,kind='linear')
    interprad = interpolate.interp2d(allages,allmasses,radius,kind='linear')
    interplogg = interpolate.interp2d(allages,allmasses,logg,kind='linear')
    interplogL = interpolate.interp2d(allages,allmasses,logL,kind='linear')  
        
    mass = np.array(mass).flatten()
    age = np.array(age).flatten()
    
    newteff = np.array([interpteff(i,j) for i,j in zip(age,mass)])
    newrad = np.array([interprad(i,j) for i,j in zip(age,mass)])
    newlogg = np.array([interplogg(i,j) for i,j in zip(age,mass)])
    newlogL = np.array([interplogL(i,j) for i,j in zip(age,mass)])
    
    stardict = {'Teff (K)':newteff,'Radius (Rs)':newrad, 'log g':newlogg, 'log L':newlogL}    
    
    starTable = Table(stardict)
    
    return starTable

def interp(iqehist, newE):
    """compute a new IQE histogram using the new energy array
    
    * iqehist: input IQE histogram
    * newE: new energy centers array
    """
    from scipy import interpolate

    mask = iqehist.I != iqehist.I
    # find energy boundaries of dynamic range for each Q
    def get_boundary_indexes(a):
        nz = np.nonzero(a)[0]
        if nz.size:
            return nz[0], nz[-1]
        else:
            return 0, 0

    boundary_indexes = [get_boundary_indexes(row) for row in np.logical_not(mask)]
    try:
        E = iqehist.energy
    except:
        E = iqehist.E
    Eranges = [(E[ind1], E[ind2]) for ind1, ind2 in boundary_indexes]
    #
    iqehist.I[mask] = 0
    iqehist.E2[mask] = 0
    Q = iqehist.Q
    f = interpolate.interp2d(E, Q, iqehist.I, kind="linear")
    E2f = interpolate.interp2d(E, Q, iqehist.E2, kind="linear")
    dE = E[1] - E[0]
    Emin = E[0] // dE * dE
    Emax = E[-1] // dE * dE
    # Enew = np.arange(Emin, Emax+dE/2, dE)
    newS = f(newE, Q)
    newS_E2 = E2f(newE, Q)
    # create new histogram
    Eaxis = H.axis("E", newE, unit="meV")
    Qaxis = H.axis("Q", Q, unit="1./angstrom")
    newHist = H.histogram("IQE", [Qaxis, Eaxis], data=newS, errors=newS_E2)
    #
    for Erange, q in zip(Eranges, Q):
        Emin, Emax = Erange
        if Emin > newE[0]:
            Emin = min(Emin, newE[-1])
            newHist[q, (None, Emin)].I[:] = np.nan
            newHist[q, (None, Emin)].E2[:] = np.nan
        if Emax < newE[-1]:
            Emax = max(Emax, newE[0])
            newHist[q, (Emax, None)].I[:] = np.nan
            newHist[q, (Emax, None)].E2[:] = np.nan
        continue
    return newHist