def __init__(self, A, h, alpha, CL, CD, rho, smoothing=0, k_spline=3):
self.A = A
self.h = h
self.Asp = 2*self.h**2/self.A
self.rho = rho
# Input lift and drag data
self.n = len(alpha)
self.alpha = alpha
self.CL = CL
self.CD = CD
self.k_spline = k_spline
# -------- Spline interpolation -----------------------------
if len(self.alpha.shape) == 1:
self.interpolationMethod = 'spline'
self.nrControlParameters = 1
self.CL_spline = interpolate.splrep(self.alpha, self.CL, s=smoothing, k=self.k_spline)
self.CD_spline = interpolate.splrep(self.alpha, self.CD, s=smoothing, k=self.k_spline)
elif len(self.alpha.shape) == 2:
self.interpolationMethod = 'griddata'
self.nrControlParameters = 2
self.CL_RbfModel = interpolate.Rbf(self.alpha[:, 0],self.alpha[:, 1], self.CL, smooth=smoothing)
self.CD_RbfModel = interpolate.Rbf(self.alpha[:, 0],self.alpha[:, 1], self.CD, smooth=smoothing)
python类Rbf()的实例源码
def CLCD(self, alpha):
if self.interpolationMethod == 'spline':
CL = interpolate.splev(alpha, self.CL_spline)
CD = interpolate.splev(alpha, self.CD_spline)
elif self.interpolationMethod == 'Rbf':
CL = self.CL_RbfModel(alpha[0], alpha[1])
CD = self.CD_RbfModel(alpha[0], alpha[1])
elif self.interpolationMethod == 'griddata':
CL = interpolate.griddata(self.alpha, self.CL, alpha)
CD = interpolate.griddata(self.alpha, self.CD, alpha)
if self.nrControlParameters == 2:
return float(CL), float(CD)
else:
return CL, CD
def populateMissingData(self,approach="Smooth",ilog=None):
'''
This function is used to interpolate missing data in the image.
'''
if approach == 'Smooth':
# first run a median filter over the array, then smooth the result.
#xmin,xmax,ymin,ymax,zmin,zmax = self.getRange()
mask = np.array(self.flags,dtype=np.bool)
z = self.getZImage().asMatrix2D()
median = nd.median_filter(z,size=(15,15))
mask = mask.flatten()
z = z.flatten()
median = median.flatten()
z[ mask==False ] = median[ mask==False ]
if ilog != None:
ilog.log(pv.Image(median.reshape(self.width,self.height)),label="Median")
ilog.log(pv.Image(z.reshape(self.width,self.height)),label="ZMedian")
mask = mask.flatten()
z = z.flatten()
median = median.flatten()
for i in range(5):
tmp = z.copy()
smooth = nd.gaussian_filter(z.reshape(self.width,self.height),2.0).flatten()
z[ mask==False ] = smooth[ mask==False ]
print "Iteration:",i,(z-tmp).max(),(z-tmp).min()
ilog.log(pv.Image(z.reshape(self.width,self.height)),label="ZSmooth%02d"%i)
ilog.log(pv.Image((z-tmp).reshape(self.width,self.height)),label="ZSmooth%02d"%i)
if approach == 'RBF':
mask = np.array(self.flags,dtype=np.bool)
mask = mask.flatten()
x = np.arange(self.width).reshape(self.width,1)
x = x*np.ones((1,self.height))
x = x.flatten()
y = np.arange(self.height).reshape(1,self.height)
y = y*np.ones((self.width,1))
y = y.flatten()
z = self.z.copy()
z = z.flatten()
print "Coords:"
print len(mask)
print len(x[mask])
print len(y[mask])
print len(z[mask])
# this produces an error. Probably has too much data
it.Rbf(x[mask],y[mask],z[mask])
pass
def compute_ck2004_response(self, path, verbose=False):
"""
Computes Castelli & Kurucz (2004) intensities across the entire
range of model atmospheres.
@path: path to the directory containing ck2004 SEDs
@verbose: switch to determine whether computing progress should
be printed on screen
Returns: n/a
"""
models = glob.glob(path+'/*M1.000*')
Nmodels = len(models)
# Store the length of the filename extensions for parsing:
offset = len(models[0])-models[0].rfind('.')
Teff, logg, abun = np.empty(Nmodels), np.empty(Nmodels), np.empty(Nmodels)
InormE, InormP = np.empty(Nmodels), np.empty(Nmodels)
if verbose:
print('Computing Castelli & Kurucz (2004) passband intensities for %s:%s. This will take a while.' % (self.pbset, self.pbname))
for i, model in enumerate(models):
#~ spc = np.loadtxt(model).T -- waaay slower
spc = np.fromfile(model, sep=' ').reshape(-1,2).T
Teff[i] = float(model[-17-offset:-12-offset])
logg[i] = float(model[-11-offset:-9-offset])/10
sign = 1. if model[-9-offset]=='P' else -1.
abun[i] = sign*float(model[-8-offset:-6-offset])/10
spc[0] /= 1e10 # AA -> m
spc[1] *= 1e7 # erg/s/cm^2/A -> W/m^3
wl = spc[0][(spc[0] >= self.ptf_table['wl'][0]) & (spc[0] <= self.ptf_table['wl'][-1])]
fl = spc[1][(spc[0] >= self.ptf_table['wl'][0]) & (spc[0] <= self.ptf_table['wl'][-1])]
fl *= self.ptf(wl)
flP = fl*wl
InormE[i] = np.log10(fl.sum()/self.ptf_area*(wl[1]-wl[0])) # energy-weighted intensity
InormP[i] = np.log10(flP.sum()/self.ptf_photon_area*(wl[1]-wl[0])) # photon-weighted intensity
if verbose:
if 100*i % (len(models)) == 0:
print('%d%% done.' % (100*i/(len(models)-1)))
# Store axes (Teff, logg, abun) and the full grid of Inorm, with
# nans where the grid isn't complete.
self._ck2004_axes = (np.unique(Teff), np.unique(logg), np.unique(abun))
self._ck2004_energy_grid = np.nan*np.ones((len(self._ck2004_axes[0]), len(self._ck2004_axes[1]), len(self._ck2004_axes[2]), 1))
self._ck2004_photon_grid = np.nan*np.ones((len(self._ck2004_axes[0]), len(self._ck2004_axes[1]), len(self._ck2004_axes[2]), 1))
for i, I0 in enumerate(InormE):
self._ck2004_energy_grid[Teff[i] == self._ck2004_axes[0], logg[i] == self._ck2004_axes[1], abun[i] == self._ck2004_axes[2], 0] = I0
for i, I0 in enumerate(InormP):
self._ck2004_photon_grid[Teff[i] == self._ck2004_axes[0], logg[i] == self._ck2004_axes[1], abun[i] == self._ck2004_axes[2], 0] = I0
# Tried radial basis functions but they were just terrible.
#~ self._log10_Inorm_ck2004 = interpolate.Rbf(self._ck2004_Teff, self._ck2004_logg, self._ck2004_met, self._ck2004_Inorm, function='linear')
self.content.append('ck2004')
self.atmlist.append('ck2004')
