我们从Python开源项目中，提取了以下44个代码示例，用于说明如何使用astropy.units.cm()。
def norm_apec(self, z): """ The normalization factor of the XSPEC APEC model assuming EM = 1 (i.e., n_e = n_H = 1 cm^-3, and V = 1 cm^3) norm = 1e-14 / (4*pi* (D_A * (1+z))^2) * int(n_e * n_H) dV unit: [cm^-5] This value will be used to calculate the cooling function values. References ---------- * XSPEC: APEC model https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/manual/XSmodelApec.html """ da = self.angular_diameter_distance(z, unit="cm") norm = 1e-14 / (4*math.pi * (da * (1+z))**2) return norm
def ionization_rate(self): """ Total ionization rate. Includes contributions from both direct ionization and excitation-autoionization See Also -------- direct_ionization_rate excitation_autoionization_rate """ di_rate = self.direct_ionization_rate() di_rate = np.zeros(self.temperature.shape)*u.cm**3/u.s if di_rate is None else di_rate ea_rate = self.excitation_autoionization_rate() ea_rate = np.zeros(self.temperature.shape)*u.cm**3/u.s if ea_rate is None else ea_rate return di_rate + ea_rate
def recombination_rate(self): """ Total recombination rate. Includes contributions from dielectronic recombination and radiative recombination. See Also -------- radiative_recombination_rate dielectronic_recombination_rate """ rr_rate = self.radiative_recombination_rate() rr_rate = np.zeros(self.temperature.shape)*u.cm**3/u.s if rr_rate is None else rr_rate dr_rate = self.dielectronic_recombination_rate() dr_rate = np.zeros(self.temperature.shape)*u.cm**3/u.s if dr_rate is None else dr_rate return rr_rate + dr_rate
def GetFluxRatio(sptlist, Tsec, xgrid): """ Returns the flux ratio between the secondary star of temperature Tsec and the (possibly multiple) primary star(s) given in the 'sptlist' list (given as spectral types) xgrid is a np.ndarray containing the x-coordinates to find the flux ratio at (in nm) """ prim_flux = np.zeros(xgrid.size) # Determine the flux from the primary star(s) for spt in sptlist: end = search("[0-9]", spt).end() T = MS.Interpolate(MS.Temperature, spt[:end]) R = MS.Interpolate(MS.Radius, spt[:end]) prim_flux += Planck(xgrid * units.nm.to(units.cm), T) * R ** 2 # Determine the secondary star flux s_spt = MS.GetSpectralType(MS.Temperature, Tsec) R = MS.Interpolate(MS.Radius, s_spt) sec_flux = Planck(xgrid * units.nm.to(units.cm), Tsec) * R ** 2 return sec_flux / prim_flux
def _filter_streamlines(self, streamline, close_threshold=0.05, loop_length_range: u.cm =[2.e+9, 5.e+10]*u.cm, **kwargs): """ Check extracted loop to make sure it fits given criteria. Return True if it passes. Parameters ---------- streamline : yt streamline object close_threshold : `float` percentage of domain width allowed between loop endpoints loop_length_range : `~astropy.Quantity` minimum and maximum allowed loop lengths (in centimeters) """ streamline = streamline[np.all(streamline != 0.0, axis=1)] loop_length = np.sum(np.linalg.norm(np.diff(streamline, axis=0), axis=1)) if np.fabs(streamline[0, 2] - streamline[-1, 2]) > close_threshold*self.extrapolated_3d_field.domain_width[2]: return False elif loop_length > loop_length_range[1].to(u.cm).value or loop_length < loop_length_range[0].to(u.cm).value: return False else: return True
def peek(self, figsize=(10, 10), color='b', alpha=0.75, print_to_file=None, **kwargs): """ Show extracted fieldlines overlaid on HMI image. """ fig = plt.figure(figsize=figsize) ax = fig.gca(projection=self.hmi_map) self.hmi_map.plot() ax.set_autoscale_on(False) for stream, _ in self.streamlines: ax.plot(self._convert_angle_to_length(stream[:, 0]*u.cm, working_units=u.arcsec).to(u.deg), self._convert_angle_to_length(stream[:, 1]*u.cm, working_units=u.arcsec).to(u.deg), alpha=alpha, color=color, transform=ax.get_transform('world')) if print_to_file is not None: plt.savefig(print_to_file, **kwargs) plt.show()
def flatten_emissivities(channel, emission_model): """ Compute product between wavelength response and emissivity for all ions """ flattened_emissivities = [] for ion in emission_model: wavelength, emissivity = emission_model.get_emissivity(ion) if wavelength is None or emissivity is None: flattened_emissivities.append(None) continue interpolated_response = splev(wavelength.value, channel['wavelength_response_spline'], ext=1) em_summed = np.dot(emissivity.value, interpolated_response) unit = emissivity.unit*u.count/u.photon*u.steradian/u.pixel*u.cm**2 flattened_emissivities.append(u.Quantity(em_summed, unit)) return flattened_emissivities
def __init__(self, *args, **kwargs): if len(args) == 1 and os.path.exists(args[0]): data, header, wavelength = self._restore_from_file(args[0], **kwargs) elif all([k in kwargs for k in ['data', 'header', 'wavelength']]): data = kwargs.get('data') header = kwargs.get('header') wavelength = kwargs.get('wavelength') else: raise ValueError('''EISCube can only be initialized with a valid FITS file or NumPy array with an associated wavelength and header.''') # check dimensions if data.shape[-1] != wavelength.shape[0]: raise ValueError('''Third dimension of data cube must have the same length as wavelength.''') self.meta = header.copy() self.wavelength = wavelength self.data = data self.cmap = kwargs.get('cmap', sunpy.cm.get_cmap('hinodexrt')) self._fix_header()
def load_results(self, loop): """ Load EBTEL output for a given loop object. Parameters ---------- loop : `synthesizAR.Loop` object """ # load text N_s = len(loop.field_aligned_coordinate) _tmp = np.loadtxt(loop.hydro_configuration['output_filename']) # reshape into a 1D loop structure with units time = _tmp[:, 0]*u.s electron_temperature = np.outer(_tmp[:, 1], np.ones(N_s))*u.K ion_temperature = np.outer(_tmp[:, 2], np.ones(N_s))*u.K density = np.outer(_tmp[:, 3], np.ones(N_s))*(u.cm**(-3)) velocity = np.outer(_tmp[:, -2], np.ones(N_s))*u.cm/u.s # flip sign of velocity at apex i_mirror = np.where(np.diff(loop.coordinates.value[:, 2]) > 0)[0][-1] + 2 velocity[:, i_mirror:] = -velocity[:, i_mirror:] return time, electron_temperature, ion_temperature, density, velocity
def non_equilibrium_ionization(self, time: u.s, temperature: u.K, density: u.cm**(-3), rate_matrix=None, initial_condition=None): """ Compute the ionization fraction in non-equilibrium for a given temperature and density timeseries. """ if rate_matrix is None: rate_matrix = self._rate_matrix() if initial_condition is None: initial_condition = self.equilibrium_ionization(rate_matrix=rate_matrix) interpolate_indices = [np.abs(self.temperature - t).argmin() for t in temperature] y = np.zeros(time.shape + (self.atomic_number+1,)) y[0, :] = initial_condition[interpolate_indices[0], :] identity = u.Quantity(np.eye(self.atomic_number + 1)) for i in range(1, time.shape[0]): dt = time[i] - time[i-1] term1 = identity - density[i]*dt/2.*rate_matrix[interpolate_indices[i], :, :] term2 = identity + density[i-1]*dt/2.*rate_matrix[interpolate_indices[i-1], :, :] y[i, :] = np.linalg.inv(term1) @ term2 @ y[i-1, :] y[i, :] = np.fabs(y[i, :]) y[i, :] /= y[i, :].sum() return u.Quantity(y)
def proton_collision_rate(self): """ Calculates the collision rate for de-exciting and exciting collisions for protons """ # Create scaled temperature--these are not stored in the file bt_t = np.vectorize(np.linspace, excluded=[0, 1], otypes='O')(0, 1, [ups.shape[0] for ups in self._psplups['bt_rate']]) # Get excitation rates directly from scaled data energy_ratio = np.outer(const.k_B.cgs*self.temperature, 1.0/self._psplups['delta_energy'].to(u.erg)) ex_rate = np.array(list(map(self.burgess_tully_descale, bt_t, self._psplups['bt_rate'], energy_ratio.T, self._psplups['bt_c'], self._psplups['bt_type']))) ex_rate = u.Quantity(np.where(ex_rate > 0., ex_rate, 0.), u.cm**3/u.s).T # Calculation de-excitation rates from excitation rate omega_upper = 2.*self._elvlc['J'][self._psplups['upper_level'] - 1] + 1. omega_lower = 2.*self._elvlc['J'][self._psplups['lower_level'] - 1] + 1. dex_rate = (omega_lower/omega_upper)*ex_rate*np.exp(1./energy_ratio) return dex_rate, ex_rate
def emissivity(self, density: u.cm**(-3), include_energy=False, **kwargs): """ Calculate emissivity for all lines as a function of temperature and density """ populations = self.level_populations(density, include_protons=kwargs.get('include_protons', True)) if populations is None: return (None, None) wavelengths = np.fabs(self._wgfa['wavelength']) # Exclude 0 wavelengths which correspond to two-photon decays upper_levels = self._wgfa['upper_level'][wavelengths != 0*u.angstrom] a_values = self._wgfa['A'][wavelengths != 0*u.angstrom] wavelengths = wavelengths[wavelengths != 0*u.angstrom] if include_energy: energy = const.h.cgs*const.c.cgs/wavelengths.to(u.cm) else: energy = 1.*u.photon emissivity = populations[:, :, upper_levels - 1]*(a_values*energy) return wavelengths, emissivity
def cm_per_pix(self, z): """ Calculate the transversal length (unit: cm) corresponding to 1 ACIS pixel (i.e., 0.492 arcsec) at the *angular diameter distance* of z. """ return self.kpc_per_pix(z) * au.kpc.to(au.cm)
def preprocessor(self, table, line, index): line = line.strip().split() if index == 0: filetype = int(line[0]) table.append([filetype]) if filetype == 1: self.dtypes = [int,float,float,float,float] self.units = [None,(u.cm**3)/u.s,None,u.K,u.K] self.headings = ['fit_type', 'A_fit', 'B_fit', 'T0_fit', 'T1_fit'] self.descriptions = ['fit type','A fit parameter','B fit parameter', 'T0 fit parameter','T1 fit parameter'] elif filetype == 2: self.dtypes = [int,float,float,float,float,float,float] self.units = [None,(u.cm**3)/u.s,None,u.K,u.K,None,u.K] self.headings = ['fit_type', 'A_fit', 'B_fit', 'T0_fit', 'T1_fit', 'C_fit', 'T2_fit'] self.descriptions = ['fit type','A fit parameter','B fit parameter', 'T0 fit parameter','T1 fit parameter','C fit parameter', 'T2 fit parameter'] elif filetype == 3: self.dtypes = [int,float,float] self.units = [None,(u.cm**3)/u.s,None] self.headings = ['fit_type', 'A_fit', 'eta_fit'] self.descriptions = ['fit type','A rad fit parameter','eta fit parameter'] else: raise ValueError('Unrecognized .rrparams filetype {}'.format(filetype)) else: if table[0][0] == 1 or table[0][0] == 2: table[0] += line[3:] else: table[0] += line[2:]
def preprocessor(self, table, line, index): line = line.strip().split() if index == 0: self._drparams_filetype = int(line[0]) if self._drparams_filetype == 1: # Badnell type self.dtypes = [int,float,float] self.units = [None,u.K,(u.cm**3)/u.s*(u.K**(3/2))] self.headings = ['fit_type', 'E_fit', 'c_fit'] self.descriptions = ['fit type', 'E fit parameter','c fit parameter'] elif self._drparams_filetype == 2: # Shull type self.dtypes = [int,float,float,float,float] self.units = [None,(u.cm**3)/u.s*(u.K**(3/2)),u.dimensionless_unscaled,u.K,u.K] self.headings = ['fit_type', 'A_fit', 'B_fit', 'T0_fit', 'T1_fit'] self.descriptions = ['fit type','A fit coefficient','B fit coefficient', 'T0 fit coefficient','T1 fit coefficient'] else: raise ValueError('Unrecognized drparams filetype {}'.format(self._drparams_filetype)) else: if self._drparams_filetype == 1: tmp = np.array(line[2:],dtype=float) if index%2 == 0: tmp_col = table[-1] for i in range(tmp.shape[0]): table.append([self._drparams_filetype,tmp_col[i],tmp[i]]) del table[0] else: table.append(tmp) else: table.append([self._drparams_filetype]+line[2:])
def get_stellar_mass(self, z=0, rnorm=1., cosmo=Planck15): """ Get M/L, L, stellar mass, as measured in R-band Returns: M/Lr, log10(Lr), log10(stellar_mass) """ from sedpy.observate import Filter import astropy.units as u import astropy.constants as const MLr, MLi, MLk = self.get_M2L() #plt.plot(self.wave, self.spec); print(MLr) rfilt = Filter('sdss_r0') norm_spec = self.get_model(in_place=False)*rnorm #rband_Jy = rfilt.obj_counts(self.wave, norm_spec)/rfilt.ab_zero_counts*3631 #rband_flam = rband_Jy/(3.34e4*rfilt.wave_pivot**2)#*u.erg/u.second/u.cm**2/u.AA #dL = Planck15.luminosity_distance(z) Mr = rfilt.ab_mag(self.wave, norm_spec) - cosmo.distmod(z).value Lr = 10**(-0.4*(Mr-rfilt.solar_ab_mag)) #Lr = (rband_flam*4*np.pi*dL.to(u.cm).value**2)/3.828e+33*rfilt.wave_pivot*(1+z) stellar_mass = Lr*MLr return MLr, np.log10(Lr), np.log10(stellar_mass)
def __init__(self, name, z, x, y): self.name = name self.z = z self.x = x self.y = y self.lDcm = cosmo.luminosity_distance(self.z)*u.Mpc.to(u.cm) / u.Mpc self.radToKpc = conv.arcsec_per_kpc_proper(self.z)*0.05/u.arcsec*u.kpc # setting up final xl file
def __init__(self, name, z, x, y): self.name = name self.z = z self.x = x self.y = y self.lDcm = cosmo.luminosity_distance(self.z)*u.Mpc.to(u.cm) / u.Mpc self.radToKpc = conv.arcsec_per_kpc_proper(self.z)*0.05/u.arcsec*u.kpc # We grab our galaxy data and make a list of galaxy objects
def __init__(self, name, z, x, y): self.name = name self.z = z self.x = x self.y = y self.lDcm = cosmo.luminosity_distance(self.z)*u.Mpc.to(u.cm) / u.Mpc self.radToKpc = conv.arcsec_per_kpc_proper(self.z)*0.05/u.arcsec*u.kpc # Grabbing and filling galaxy data
def test_init(self): afrho = Afrho(1000 * u.cm) assert afrho.value == 1000 assert afrho.unit == u.cm
def test_scaler_ops(self): afrho = Afrho(1000 * u.cm) afrho = afrho / 2 assert afrho == 500 * u.cm
def test_quantity_ops(self): afrho = Afrho(1000 * u.cm) v = afrho * 2 * u.cm assert v == 2000 * u.cm**2 assert not isinstance(v, Afrho)
def test_from_flam(self): fluxd = 6.730018324465526e-14 * u.W / u.m**2 / u.um aper = 1 * u.arcsec eph = dict(rh=1.5 * u.au, delta=1.0 * u.au) S = 1869 * u.W / u.m**2 / u.um afrho = Afrho.from_fluxd(None, fluxd, aper, eph, S=S) assert np.isclose(afrho.cm, 1000)
def test_from_fnu(self): fluxd = 6.161081515869728 * u.mJy nu = 2.998e14 / 11.7 * u.Hz aper = 1 * u.arcsec eph = dict(rh=1.5 * u.au, delta=1.0 * u.au) S = 1.711e14 * u.Jy afrho = Afrho.from_fluxd(nu, fluxd, aper, eph, S=S) assert np.isclose(afrho.cm, 1000.0)
def test_to_phase(self): afrho = Afrho(10 * u.cm).to_phase(15 * u.deg, 0 * u.deg) assert np.isclose(afrho.cm, 5.8720)
def test_init(self): efrho = Efrho(1000 * u.cm) assert efrho.value == 1000 assert efrho.unit == u.cm
def test_scaler_ops(self): efrho = Efrho(1000 * u.cm) efrho = efrho / 2 assert efrho == 500 * u.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 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)
