a ߙfbe@sdZddlZddlZddlmZddlmZddl m Z ddl m Z ddl mZmZgd ZGd d d e ZGd d d e ZGddde ZGddde ZdS)z$ Models that have physical origins. N) constants)units)AstropyUserWarning)Fittable1DModel) ParameterInputParameterError) BlackBodyDrude1D Plummer1DNFWcseZdZdZeddejddZeddddZd Z d e iZ ej ej d ejejejZej ej d ejejejej ej d ejejejd Zfd dZddZeddZddZeddZeddZeddZZS)r aZ Blackbody model using the Planck function. Parameters ---------- temperature : `~astropy.units.Quantity` ['temperature'] Blackbody temperature. scale : float or `~astropy.units.Quantity` ['dimensionless'] Scale factor. If dimensionless, input units will assumed to be in Hz and output units in (erg / (cm ** 2 * s * Hz * sr). If not dimensionless, must be equivalent to either (erg / (cm ** 2 * s * Hz * sr) or erg / (cm ** 2 * s * AA * sr), in which case the result will be returned in the requested units and the scale will be stripped of units (with the float value applied). Notes ----- Model formula: .. math:: B_{\nu}(T) = A \frac{2 h \nu^{3} / c^{2}}{exp(h \nu / k T) - 1} Examples -------- >>> from astropy.modeling import models >>> from astropy import units as u >>> bb = models.BlackBody(temperature=5000*u.K) >>> bb(6000 * u.AA) # doctest: +FLOAT_CMP .. plot:: :include-source: import numpy as np import matplotlib.pyplot as plt from astropy.modeling.models import BlackBody from astropy import units as u from astropy.visualization import quantity_support bb = BlackBody(temperature=5778*u.K) wav = np.arange(1000, 110000) * u.AA flux = bb(wav) with quantity_support(): plt.figure() plt.semilogx(wav, flux) plt.axvline(bb.nu_max.to(u.AA, equivalencies=u.spectral()).value, ls='--') plt.show() g@rzBlackbody temperaturedefaultminunit description?z Scale factorrrrTx)SNUZSLAMcs~|dd}t|drd|jtjsd|j}||jtdtjsRt d||j |d<||_ n|j|_ t j |i|S)Nscalerrz8scale units not dimensionless or in surface brightness: )gethasattrr is_equivalentudimensionless_unscaled _native_unitsspectral_densityAA ValueErrorvalue _output_unitssuper__init__)selfargskwargsrZ output_units __class__?lib/python3.9/site-packages/astropy/modeling/physical_models.pyr$_s  zBlackBody.__init__c Cst|tjst|tj}n|}t|tjs@t||jd}n|}ttt2tj|tjt j d}t|tj}Wdn1s0Yt |dkrt d|t t |rt |dkrtdttj|tj|}t |} dtj|dtjd | tj} |jjdurZt|d sL||jj}|tjj}|| |jt|} t|d r| S| jS) aEvaluate the model. Parameters ---------- x : float, `~numpy.ndarray`, or `~astropy.units.Quantity` ['frequency'] Frequency at which to compute the blackbody. If no units are given, this defaults to Hz (or AA if `scale` was initialized with units equivalent to erg / (cm ** 2 * s * AA * sr)). temperature : float, `~numpy.ndarray`, or `~astropy.units.Quantity` Temperature of the blackbody. If no units are given, this defaults to Kelvin. scale : float, `~numpy.ndarray`, or `~astropy.units.Quantity` ['dimensionless'] Desired scale for the blackbody. Returns ------- y : number or ndarray Blackbody spectrum. The units are determined from the units of ``scale``. .. note:: Use `numpy.errstate` to suppress Numpy warnings, if desired. .. warning:: Output values might contain ``nan`` and ``inf``. Raises ------ ValueError Invalid temperature. ZeroDivisionError Wavelength is zero (when converting to frequency). r)ZdtypeNrz Temperature should be positive: z4Input contains invalid wavelength/frequency value(s)g@rr) isinstancerQuantityK input_unitsZadd_enabled_equivalenciesspectral temperatureHznpZfloat64anyr allZisfinitewarningswarnrconsthk_BZexpm1csrrrrtorr!r"r) r%rr2rZin_tempZin_xZfreqZtempZ log_boltzZboltzm1Zbb_nuyr*r*r+evaluateps6'  , &   zBlackBody.evaluatecCs6|j|jdr"|jdtjiS|jdtjiSdS)Nrr)r"r_native_output_unitsinputsrr3rr%r*r*r+r0szBlackBody.input_unitscCs dtjiS)Nr2)rr/r%Z inputs_unitZ outputs_unitr*r*r+_parameter_units_for_data_unitssz)BlackBody._parameter_units_for_data_unitscCs\|jjdur|jjtj}n|jj}|tj|j dt j }|tj tj dtjS)zBolometric flux.Nr)rrZquantityr>rrr!r9Zsigma_sbr2r4piergcms)r%rZnative_bolfluxr*r*r+bolometric_fluxs  zBlackBody.bolometric_fluxcCs tj|jS)z=Peak wavelength when the curve is expressed as power density.)r9Zb_wienr2rCr*r*r+ lambda_maxszBlackBody.lambda_maxcCsdtj|jtjS)zN@)r9r;r2r:rCr*r*r+nu_maxszBlackBody.nu_max)__name__ __module__ __qualname____doc__rrr/r2r _input_units_allow_dimensionlessr1Zinput_units_equivalenciesrHrIrJr3r=rrrAr$r@propertyr0rErKrLrM __classcell__r*r*r(r+r s(5 "   Z   r c@seZdZdZedddZedddZedddZeddZ ed d Z e d d Z d dZ e ddZejddZdddZdS)r a@ Drude model based one the behavior of electons in materials (esp. metals). Parameters ---------- amplitude : float Peak value x_0 : float Position of the peak fwhm : float Full width at half maximum Model formula: .. math:: f(x) = A \frac{(fwhm/x_0)^2}{((x/x_0 - x_0/x)^2 + (fwhm/x_0)^2} Examples -------- .. plot:: :include-source: import numpy as np import matplotlib.pyplot as plt from astropy.modeling.models import Drude1D fig, ax = plt.subplots() # generate the curves and plot them x = np.arange(7.5 , 12.5 , 0.1) dmodel = Drude1D(amplitude=1.0, fwhm=1.0, x_0=10.0) ax.plot(x, dmodel(x)) ax.set_xlabel('x') ax.set_ylabel('F(x)') plt.show() rz Peak ValuerrzPosition of the peakzFull width at half maximumcCs0|||d||||d||dS)z6 One dimensional Drude model function rr*)r amplitudex_0fwhmr*r*r+r@ s  zDrude1D.evaluatecCs||d||||d||d}d||d|||d|d| |d|||||d|d|d}d|||d|}|||gS)z5 Drude1D model function derivatives. rrr,r*)rrVrWrXZ d_amplitudeZd_x_0Zd_fwhmr*r*r+ fit_deriv+s&, zDrude1D.fit_derivcCs"|jjdurdS|jd|jjiSNr)rWrrBrCr*r*r+r0Bs zDrude1D.input_unitscCs*||jd||jd||jddS)Nr)rWrXrV)rBoutputsrDr*r*r+rEHs   z'Drude1D._parameter_units_for_data_unitscCs"|jjdurdS|jd|jjiSr[)rVrr\rCr*r*r+ return_unitsOs zDrude1D.return_unitscCst|dkrtddS)z Ensure `x_0` is not 0.rz!0 is not an allowed value for x_0N)r4r5r)r%valr*r*r+rWUsz Drude1D.x_02cCs |j}||j}||||fS)zTuple defining the default ``bounding_box`` limits, ``(x_low, x_high)``. Parameters ---------- factor : float The multiple of FWHM used to define the limits. )rWrX)r%ZfactorZx0Zdxr*r*r+ bounding_box[s  zDrude1D.bounding_boxN)r_)rNrOrPrQrrVrWrX staticmethodr@rZrSr0rEr]Z validatorr`r*r*r*r+r s )        r c@sTeZdZdZedddZedddZeddZedd Z e d d Z d d Z dS)r aOne dimensional Plummer density profile model. Parameters ---------- mass : float Total mass of cluster. r_plum : float Scale parameter which sets the size of the cluster core. Notes ----- Model formula: .. math:: \rho(r)=\frac{3M}{4\pi a^3}(1+\frac{r^2}{a^2})^{-5/2} References ---------- .. [1] https://ui.adsabs.harvard.edu/abs/1911MNRAS..71..460P rzTotal mass of clusterrUz7Scale parameter which sets the size of the cluster corecCs.d|dtj|dd||ddS)z9 Evaluate plummer density profile model. r,rFrrgr4rG)rmassr_plumr*r*r+r@szPlummer1D.evaluatecCsxddtj|d||ddd}d||dd||ddtj|dd||dd}||gS) z. Plummer1D model derivatives. r,rFrrg@ g @rb)rrcrdZd_massZd_r_plumr*r*r+rZs *.zPlummer1D.fit_derivcCs2|jjdur|jjdurdS|jd|jjiSdSr[)rcrrdrBrCr*r*r+r0szPlummer1D.input_unitscCs0||jd||jdd||jddS)Nrr,)rcrd)r\rBrDr*r*r+rEs z)Plummer1D._parameter_units_for_data_unitsN) rNrOrPrQrrcrdrar@rZrSr0rEr*r*r*r+r js     r cseZdZdZeddejddZeddddZeddddZ d Z e ej ej ej e j d d ffd d ZddZddZeddZddZeddZddZeddZeddZeddZed d!Zd"d#Zed$d%Zed&d'Zd(d)ZZS)*r u Navarro–Frenk–White (NFW) profile - model for radial distribution of dark matter. Parameters ---------- mass : float or `~astropy.units.Quantity` ['mass'] Mass of NFW peak within specified overdensity radius. concentration : float Concentration of the NFW profile. redshift : float Redshift of the NFW profile. massfactor : tuple or str Mass overdensity factor and type for provided profiles: Tuple version: ("virial",) : virial radius ("critical", N) : radius where density is N times that of the critical density ("mean", N) : radius where density is N times that of the mean density String version: "virial" : virial radius "Nc" : radius where density is N times that of the critical density (e.g. "200c") "Nm" : radius where density is N times that of the mean density (e.g. "500m") cosmo : :class:`~astropy.cosmology.Cosmology` Background cosmology for density calculation. If None, the default cosmology will be used. Notes ----- Model formula: .. math:: \rho(r)=\frac{\delta_c\rho_{c}}{r/r_s(1+r/r_s)^2} References ---------- .. [1] https://arxiv.org/pdf/astro-ph/9508025 .. [2] https://en.wikipedia.org/wiki/Navarro%E2%80%93Frenk%E2%80%93White_profile .. [3] https://en.wikipedia.org/wiki/Virial_mass rz-Peak mass within specified overdensity radiusr Z ConcentrationrgZRedshiftT)criticalNc  s|durddlm}|}||||t|tjsFt|tj}n|}|||| ||t j f|||d|dS)Nr)default_cosmologyrc concentrationredshift) Zastropy.cosmologyrir_density_deltar-rr.M_sun _radius_s _density_sr#r$) r%rcrkrl massfactorcosmor'riin_massr(r*r+r$s    z NFW.__init__cCspt|dr|}nt|tj}|||||j}||||td|d}t|drf|S|jSdS)a2 One dimensional NFW profile function Parameters ---------- r : float or `~astropy.units.Quantity` ['length'] Radial position of density to be calculated for the NFW profile. mass : float or `~astropy.units.Quantity` ['mass'] Mass of NFW peak within specified overdensity radius. concentration : float Concentration of the NFW profile. redshift : float Redshift of the NFW profile. Returns ------- density : float or `~astropy.units.Quantity` ['density'] NFW profile mass density at location ``r``. The density units are: [``mass`` / ``r`` ^3] Notes ----- .. warning:: Output values might contain ``nan`` and ``inf``. rrrN) rrr.kpcror>rrpr!)r%rrcrkrlin_rZradius_reducedZdensityr*r*r+r@s   z NFW.evaluatec Cst|tr|ddkr,d}|d}n\|ddkrNt|d}d}n:|ddkrpt|d}d}ntd t|dd nzp|dkrd}|}nR|d dks|d dkrt|dd }|d }ntd t|d Wn,ttfy&td t|dYn0|dkrt||d}dt j dd|d|d}|| ||_ n@|dkr|| ||_ n$|dkr|| ||||_ |j S)z* Calculate density delta. rZvirialNrgrr<Zmeanmz Massfactor 'z,' not one of 'critical', 'mean', or 'virial'z Massfactor z/ string not of the form '#m', '#c', or 'virial'z not a tuple or stringrg2@rgT@gC@) r-tuplelowerfloatr strAttributeError TypeErrorZOmr4rGZcritical_density density_delta)r%rqrrrlZdeltaZmasstypeZOm_cZd_cr*r*r+rm(sD        "  zNFW._density_deltacCstd||d|S)z Dimensionless volume integral of the NFW profile, used as an intermediate step in some calculations for this model. Notes ----- Model formula: .. math:: A_{NFW} = [\ln(1+y) - \frac{y}{1+y}] r)r4log)r?r*r*r+A_NFWbs z NFW.A_NFWcCsPt|tjst|tj}n|}|dtj|||d|||_|jS)z= Calculate scale density of the NFW profile. @r,) r-rr.rnr4rGror density_sr%rcrkrsr*r*r+rpqs zNFW._density_scCs|jS)ze Scale density of the NFW profile. Often written in the literature as :math:`\rho_s` )rrCr*r*r+ rho_scalesz NFW.rho_scalecCsPt|tjst|tj}n|}d|dtj|jd||_|jtj S)z< Calculate scale radius of the NFW profile. g@rgUUUUUU?) r-rr.rnr4rGrradius_sr>rtrr*r*r+ros z NFW._radius_scCs|jS)z2 Scale radius of the NFW profile. )rrCr*r*r+r_sszNFW.r_scCs |j|jS)zn Mass factor defined virial radius of the NFW profile (R200c for M200c, Rvir for Mvir, etc.). )rrkrCr*r*r+r_virialsz NFW.r_virialcCs |jdS)z6 Radius of maximum circular velocity. g7L@)rrCr*r*r+r_maxsz NFW.r_maxcCs ||jS)z, Maximum circular velocity. )circular_velocityrrCr*r*r+v_maxsz NFW.v_maxc Cst|dr|}nt|tj}t|jtj |j d|jj tj d|j }||j |j }t|d| |j||| |j}| tjtj S)a Circular velocities of the NFW profile. Parameters ---------- r : float or `~astropy.units.Quantity` ['length'] Radial position of velocity to be calculated for the NFW profile. Returns ------- velocity : float or `~astropy.units.Quantity` ['speed'] NFW profile circular velocity at location ``r``. The velocity units are: [km / s] Notes ----- Model formula: .. math:: v_{circ}(r)^2 = \frac{1}{x}\frac{\ln(1+cx)-(cx)/(1+cx)}{\ln(1+c)-c/(1+c)} .. math:: x = r/r_s .. warning:: Output values might contain ``nan`` and ``inf``. rr,r)rrr.rtr4Zsqrtrcr9Gr>rrJrrrkZkm)r%rurvZ v_profileZreduced_radiusZvelocityr*r*r+rs ,zNFW.circular_velocitycCs|jdtjiSr[)rBrrtrCr*r*r+r0szNFW.input_unitscCsZ|jjdur0|jdtj|j|jddiS|jd|jj|j|jddiSdS)Nrr,)rcrr\rrnr0rBrCr*r*r+r]s $zNFW.return_unitscCstjdddS)Nrj)rrnrDr*r*r+rEsz#NFW._parameter_units_for_data_units) rNrOrPrQrrrnrcrkrlrRr.rrr$r@rmrarrprSrrorrrrrr0r]rErTr*r*r(r+r sB. 0:      1  r )rQr7Znumpyr4Zastropyrr9rrZastropy.utils.exceptionsrcorer parametersrr__all__r r r r r*r*r*r+s    _x7