def test_scaler_ops(self):
efrho = Efrho(1000 * u.cm)
efrho = efrho / 2
assert efrho == 500 * u.cm
python类cm()的实例源码
def test_quantity_ops(self):
efrho = Efrho(1000 * u.cm)
v = efrho * 2 * u.cm
assert v == 2000 * u.cm**2
assert not isinstance(v, Efrho)
def test_from_flam(self):
fluxd = 3.824064455850402e-15 * u.W / u.m**2 / u.um
wave = 10 * u.um
aper = 1 * u.arcsec
eph = dict(rh=1.5 * u.au, delta=1.0 * u.au)
efrho = Efrho.from_fluxd(wave, fluxd, aper, eph)
assert np.isclose(efrho.cm, 1000)
def test_fluxd(self):
efrho = Efrho(260.76955995377915, 'cm')
wave = 5 * u.um
aper = 1 * u.arcsec
eph = dict(rh=1.5 * u.au, delta=1.0 * u.au)
fluxd = efrho.fluxd(wave, aper, eph)
assert np.isclose(fluxd.value, 1e-16)
def to_phase(self, to_phase, from_phase, Phi=None):
"""Scale Af? to another phase angle.
Parameters
----------
to_phase : `~astropy.units.Quantity`
New target phase angle.
from_phase : `~astropy.units.Quantity`
Current target phase angle.
Phi : callable or `None`
Phase function, a callable object that takes a single
parameter, phase angle as a `~astropy.units.Quantity`, and
returns a scale factor. If `None`, `phase_HalleyMarcus` is
used. The phase function is expected to be 1.0 at 0 deg.
Returns
-------
afrho : `~Afrho`
The scaled Af? quantity.
Examples
--------
>>> from sbpy.activity import Afrho
>>> afrho = Afrho(10 * u.cm).to_phase(15 * u.deg, 0 * u.deg)
>>> afrho.cm # doctest: +FLOAT_CMP
5.87201
"""
if Phi is None:
Phi = phase_HalleyMarcus
return self * Phi(to_phase) / Phi(from_phase)
def _transform_to_yt(self, map_3d, zrange, boundary_clipping=(0, 0, 0)):
"""
Reshape data structure to something yt can work with.
Parameters
----------
map_3d : `np.array`
3D+x,y,z array holding the x,y,z components of the extrapolated field
zrange : `astropy.Quantity`
Spatial range of the extrapolated field
boundary_clipping : `tuple`, optional
The extrapolated volume has a layer of ghost cells in each dimension. This tuple of
(nx,ny,nz) tells how many cells to contract the volume and map in each direction.
"""
# reshape the magnetic field data
_tmp = map_3d[boundary_clipping[0]:-boundary_clipping[0],
boundary_clipping[1]:-boundary_clipping[1],
boundary_clipping[2]:-boundary_clipping[2], :]
# some annoying and cryptic translation between yt and SunPy
data = dict(
Bx=(np.swapaxes(_tmp[:, :, :, 1], 0, 1), 'T'),
By=(np.swapaxes(_tmp[:, :, :, 0], 0, 1), 'T'),
Bz=(np.swapaxes(_tmp[:, :, :, 2], 0, 1), 'T'))
# trim the boundary hmi map appropriately
lcx, rcx = self.hmi_map.xrange + self.hmi_map.scale.axis1*u.Quantity([boundary_clipping[0], -boundary_clipping[0]], u.pixel)
lcy, rcy = self.hmi_map.yrange + self.hmi_map.scale.axis2*u.Quantity([boundary_clipping[1], -boundary_clipping[1]], u.pixel)
bottom_left = SkyCoord(lcx, lcy, frame=self.hmi_map.coordinate_frame)
top_right = SkyCoord(rcx, rcy, frame=self.hmi_map.coordinate_frame)
self.clipped_hmi_map = self.hmi_map.submap(bottom_left, top_right)
# create the bounding box
zscale = np.diff(zrange)[0]/(map_3d.shape[2]*u.pixel)
clipped_zrange = zrange + zscale*u.Quantity([boundary_clipping[2]*u.pixel, -boundary_clipping[2]*u.pixel])
bbox = np.array([self._convert_angle_to_length(self.clipped_hmi_map.xrange).value,
self._convert_angle_to_length(self.clipped_hmi_map.yrange).value,
self._convert_angle_to_length(clipped_zrange).value])
# assemble the dataset
return yt.load_uniform_grid(data, data['Bx'][0].shape, bbox=bbox, length_unit=yt.units.cm,
geometry=('cartesian', ('x', 'y', 'z')))
def calculate_counts_simple(channel, loop, *args):
"""
Calculate the AIA intensity using only the temperature response functions.
"""
response_function = (splev(np.ravel(loop.electron_temperature), channel['temperature_response_spline'])
* u.count*u.cm**5/u.s/u.pixel)
counts = np.reshape(np.ravel(loop.density**2)*response_function, loop.density.shape)
return counts
def __init__(self, name, coordinates, field_strength):
# set unique label for loop object
self.name = name
# Load in cartesian coordinates, assign units as centimeters
self.coordinates = coordinates*u.cm
# Load in field strength along the field line; convert from Tesla to Gauss
self.field_strength = (np.array(field_strength)*u.T).to(u.Gauss)
def __init__(self, density: u.cm**(-3), *args, **kwargs):
super().__init__(*args, **kwargs)
self.temperature = self[0].temperature
self.density = density
self.resolved_wavelengths = kwargs.get('resolved_wavelengths', {})
# Cannot have empty abundances so replace them as needed
default_abundance = kwargs.get('default_abundance_dataset', 'sun_photospheric_2009_asplund')
for ion in self._ion_list:
if ion.abundance is None:
ion._dset_names['abundance_filename'] = default_abundance
def __call__(self, T, metal, vsini=0.0, return_xypoint=True, **kwargs):
"""
Given parameters, return an interpolated spectrum
If return_xypoint is False, then it will only return
a numpy.ndarray with the spectrum
Before interpolating, we will do some error checking to make
sure the requested values fall within the grid
"""
# Scale the requested values
T = (T - self.T_scale[0]) / self.T_scale[1]
metal = (metal - self.metal_scale[0]) / self.metal_scale[1]
# Get the minimum and maximum values in the grid
T_min = min(self.grid[:, 0])
T_max = max(self.grid[:, 0])
metal_min = min(self.grid[:, 1])
metal_max = max(self.grid[:, 1])
input_list = (T, metal)
# Check to make sure the requested values fall within the grid
if (T_min <= T <= T_max and
metal_min <= metal <= metal_max):
y = self.interpolator(input_list)
else:
if self.debug:
warnings.warn("The requested parameters fall outside the model grid. Results may be unreliable!")
print(T, T_min, T_max)
print(metal, metal_min, metal_max)
y = self.NN_interpolator(input_list)
# Test to make sure the result is valid. If the requested point is
# outside the Delaunay triangulation, it will return NaN's
if np.any(np.isnan(y)):
if self.debug:
warnings.warn("Found NaNs in the interpolated spectrum! Falling back to Nearest Neighbor")
y = self.NN_interpolator(input_list)
model = DataStructures.xypoint(x=self.xaxis, y=y)
vsini *= units.km.to(units.cm)
model = Broaden.RotBroad(model, vsini, linear=self.rebin)
# Return the appropriate object
if return_xypoint:
return model
else:
return model.y
def __init__(self, filename_cube, filename_white=None, pixelsize=0.2 * u.arcsec, n_fig=1,
flux_units=1E-20 * u.erg / u.s / u.cm ** 2 / u.angstrom, vmin=None, vmax=None, wave_cal='air'):
"""
Parameters
----------
filename_cube: string
Name of the MUSE datacube .fits file
filename_white: string
Name of the MUSE white image .fits file
pixel_size : float or Quantity, optional
Pixel size of the datacube, if float it assumes arcsecs.
Default is 0.2 arcsec
n_fig : int, optional
XXXXXXXX
flux_units : Quantity
XXXXXXXXXX
"""
# init
self.color = False
self.cmap = ""
self.flux_units = flux_units
self.n = n_fig
plt.close(self.n)
self.wave_cal = wave_cal
self.filename = filename_cube
self.filename_white = filename_white
self.load_data()
self.white_data = fits.open(self.filename_white)[1].data
self.hdulist_white = fits.open(self.filename_white)
self.white_data = np.where(self.white_data < 0, 0, self.white_data)
if not vmin:
self.vmin=np.nanpercentile(self.white_data,0.25)
else:
self.vmin = vmin
if not vmax:
self.vmax=np.nanpercentile(self.white_data,98.)
else:
self.vmax = vmax
self.gc2 = aplpy.FITSFigure(self.filename_white, figure=plt.figure(self.n))
self.gc2.show_grayscale(vmin=self.vmin, vmax=self.vmax)
# self.gc = aplpy.FITSFigure(self.filename, slices=[1], figure=plt.figure(20))
self.pixelsize = pixelsize
gc.enable()
# plt.close(20)
print("MuseCube: Ready!")
def __init__(self, filename_cube, filename_white=None, pixelsize=0.2 * u.arcsec, n_fig=1,
flux_units=1E-20 * u.erg / u.s / u.cm ** 2 / u.angstrom, vmin=None, vmax=None, wave_cal='air'):
"""
Parameters
----------
filename_cube: string
Name of the MUSE datacube .fits file
filename_white: string
Name of the MUSE white image .fits file
pixel_size : float or Quantity, optional
Pixel size of the datacube, if float it assumes arcsecs.
Default is 0.2 arcsec
n_fig : int, optional
XXXXXXXX
flux_units : Quantity
XXXXXXXXXX
"""
# init
self.color = False
self.cmap = ""
self.flux_units = flux_units
self.n = n_fig
plt.close(self.n)
self.wave_cal = wave_cal
self.filename = filename_cube
self.filename_white = filename_white
self.load_data()
self.white_data = fits.open(self.filename_white)[1].data
self.hdulist_white = fits.open(self.filename_white)
self.white_data = np.where(self.white_data < 0, 0, self.white_data)
if not vmin:
self.vmin=np.nanpercentile(self.white_data,0.25)
else:
self.vmin = vmin
if not vmax:
self.vmax=np.nanpercentile(self.white_data,98.)
else:
self.vmax = vmax
self.gc2 = aplpy.FITSFigure(self.filename_white, figure=plt.figure(self.n))
self.gc2.show_grayscale(vmin=self.vmin, vmax=self.vmax)
# self.gc = aplpy.FITSFigure(self.filename, slices=[1], figure=plt.figure(20))
self.pixelsize = pixelsize
gc.enable()
# plt.close(20)
print("MuseCube: Ready!")
def from_fluxd(cls, wave_or_freq, fluxd, aper, eph, Tscale=1.1,
T=None):
"""Initialize from flux density.
Assumes the small angle approximation.
Parameters
----------
wave_or_freq : `~astropy.units.Quantity`
Wavelengths or frequencies of the observation.
fluxd : `~astropy.units.Quantity`
Flux density per unit wavelength or frequency.
aper : `~astropy.units.Quantity`, `~Aperture`
Aperture of the observation as a circular radius (length
or angular units), or as an sbpy `Aperture` class.
eph : dictionary-like or `~Ephem`
Ephemerides of the comet, describing heliocentric and
geocentric distances as `~astropy.units.Quantity` via
keywords `rh` and `delta`. `rh` is not required when `aper`
is in units of length.
Tscale : float, optional
Blackbody temperature scale factor. Ignored if `T` is
provided.
T : `~astropy.units.Quantity`, optional
Use this temperature for the Planck function.
Example
-------
>>> from sbpy.activity import Efrho
>>> import astropy.units as u
>>> wave = 15.8 * u.um
>>> fluxd = 6.52 * u.mJy
>>> aper = 11.1 * u.arcsec
>>> eph = {'rh': 4.42 * u.au, 'delta': 4.01 * u.au}
>>> efrho = Efrho.from_fluxd(wave, fluxd, aper, eph=eph)
>>> efrho.cm # doctest: +FLOAT_CMP
120.00836963059808
"""
fluxd1cm = Efrho(1 * u.cm).fluxd(wave_or_freq, aper, eph=eph,
Tscale=Tscale, T=T)
fluxd1cm = fluxd1cm.to(fluxd.unit, u.spectral_density(wave_or_freq))
return Efrho((fluxd / fluxd1cm).decompose() * u.cm)
def make_slope_map(self, temperature_bounds=None, em_threshold=None, rsquared_tolerance=0.5):
"""
Create map of emission measure slopes by fitting :math:`\mathrm{EM}\sim T^a` for a
given temperature range. Only those pixels for which the minimum :math:`\mathrm{EM}`
across all temperature bins is above some threshold value.
.. warning:: This method provides no measure of the goodness of the fit. Some slope values
may not provide an accurate fit to the data.
"""
if temperature_bounds is None:
temperature_bounds = u.Quantity((1e6, 4e6), u.K)
if em_threshold is None:
em_threshold = u.Quantity(1e25, u.cm**(-5))
# cut on temperature
temperature_bin_centers = (self.temperature_bin_edges[:-1] + self.temperature_bin_edges[1:])/2.
index_temperature_bounds = np.where(np.logical_and(temperature_bin_centers >= temperature_bounds[0],
temperature_bin_centers <= temperature_bounds[1]))
temperature_fit = temperature_bin_centers[index_temperature_bounds].value
# unwrap to 2D and threshold
data = self.as_array()*u.Unit(self[0].meta['bunit'])
flat_data = data.reshape(np.prod(data.shape[:2]), temperature_bin_centers.shape[0])
index_data_threshold = np.where(np.min(flat_data[:,index_temperature_bounds[0]], axis=1) >= em_threshold)
flat_data_threshold = flat_data.value[index_data_threshold[0],:][:,index_temperature_bounds[0]]
# very basic but vectorized fitting
_, rss_flat, _, _, _ = np.polyfit(np.log10(temperature_fit), np.log10(flat_data_threshold.T), 0, full=True)
coefficients, rss, _, _, _ = np.polyfit(np.log10(temperature_fit), np.log10(flat_data_threshold.T), 1, full=True)
rsquared = 1. - rss/rss_flat
slopes = np.where(rsquared >= rsquared_tolerance, coefficients[0], 0.)
# rebuild into a map
slopes_flat = np.zeros(flat_data.shape[0])
slopes_flat[index_data_threshold[0]] = slopes
slopes_2d = np.reshape(slopes_flat, data.shape[:2])
base_meta = self[0].meta.copy()
base_meta['temp_a'] = temperature_fit[0]
base_meta['temp_b'] = temperature_fit[-1]
base_meta['bunit'] = ''
base_meta['detector'] = r'$\mathrm{EM}(T)$ slope'
base_meta['comment'] = 'Linear fit to log-transformed LOS EM'
plot_settings = self[0].plot_settings.copy()
plot_settings['norm'] = None
return GenericMap(slopes_2d, base_meta, plot_settings=plot_settings)
