Getting Error "can Only Apply 'exp' Function To Dimensionless Quantities", Not Sure How To Fix This

I have been working on converting some MatLab code over to python for one of my professors (not an assignment just working on putting together some stuff) and I am stuck on this on

Solution 1:

The astropy constants are instances of classes. Try extracting the "value" for each before using them as an argument to np.exp():

import astropy.constantsasconstimport numpy as np

h = const.h.value
c = const.c.value
k = const.k_B.value
l = np.linspace(0, 1.5e-6, 1500);
T1 = 3750IM0 = ((2*h*c**2)/(l**5))/(np.exp(h*c/( h*c/(k*T1*l)))-1)

However please note there are numerical problems with IM0. The denominator is zero over all l.

Solution 2:

The accepted answer is fine and works so far as you're confident that the bare values you're using are in the correct dimensions. But in general using .value and throwing away the units information can be potentially dangerous--using Quantities with units ensures that all your calculations are being done with compatible units.

Let's look at just the exponent in your exponential, which normally should be a dimensionless quantity.

First note that all the constants you're using from Astropy have units:

>>> from astropy.constants import h, c, k_B                                     
>>> h                                                                           
<<class'astropy.constants.codata2018.CODATA2018'> name='Planck constant' value=6.62607015e-34 uncertainty=0.0 unit='J s' reference='CODATA 2018'>
>>> c                                                                           
<<class'astropy.constants.codata2018.CODATA2018'> name='Speed of light in vacuum' value=299792458.0 uncertainty=0.0 unit='m / s' reference='CODATA 2018'>
>>> k_B                                                                         
<<class'astropy.constants.codata2018.CODATA2018'> name='Boltzmann constant' value=1.380649e-23 uncertainty=0.0 unit='J / K' reference='CODATA 2018'>

You then declared some unitless values and mixed them in with those:

>>>T1 = 3750>>>l = np.linspace(0, 1.5e-6, 1500)>>>h*c/(h*c/(k_B*T1*l))                                                        
/home/embray/.virtualenvs/astropy/lib/python3.6/site-packages/astropy/units/quantity.py:481: RuntimeWarning: divide by zero encountered in true_divide
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
<Quantity [0.00000000e+00, 5.18088768e-29, 1.03617754e-28, ...,
           7.75578885e-26, 7.76096974e-26, 7.76615063e-26] J / K>

Giving a result in joules per kelvin, from k_B which needs to be canceled out with some values in the correct unit. I'm guessing T1 is supposed to be temperature in Kelvin (I'm not sure about l but let's say it's a thermodynamic beta in J though you should double-check whatever units this value is supposed to be in).

So what you probably want to do is declare these values with appropriate units (as an aside, you can avoid the annoying divide-by-zero by defining some epsilon and using that as the start of your range):

>>>from astropy import units as u>>>eps = np.finfo(float).eps>>>T1 = 3750 * u.K                                                             >>>l = np.linspace(eps, 1.5e-6, 1500) * (u.J**-1)

Now your exponent is properly a dimensionless quantity:

>>>h*c/(h*c/(k_B*T1*l))                                                        
<Quantity [1.14962123e-35, 5.18088882e-29, 1.03617765e-28, ...,
           7.75578885e-26, 7.76096974e-26, 7.76615063e-26]>
>>>np.exp(h*c/(h*c/(k_B*T1*l)))                                                
<Quantity [1., 1., 1., ..., 1., 1., 1.]>

(In this case the dimensionless values were all so close to zero that the exponents round to 1--if this is not correct then you'll want to check some of my assumptions about your units).

In any case, this is how the library is intended to be used, and the error you're getting is a deliberate safecheck against your assumptions.

Update: I saw in your other question that you gave some more context to your problem, in particular specifying that l are wavelengths in meters (this was my first guess but I wasn't sure based on the equation you gave).

Actually you can avoid directly using h and c in your Plank equation by taking advantage of equivalencies. Here you can define l as wavelengths in meters:

>>>l = np.linspace(eps, 1.5e-6, 1500) * u.m

and convert this directly to spectral energy:

>>>E = l.to(u.J, equivalencies=u.spectral())                                   >>>E                                                                           
<Quantity [8.94615682e-10, 1.98512112e-16, 9.92560670e-17, ...,
           1.32606651e-19, 1.32518128e-19, 1.32429724e-19] J>

Then write the exponent in your Plank's law equation like:

>>> np.exp(E / (k_B * T1))                                                      
/home/embray/.virtualenvs/astropy/lib/python3.6/site-packages/astropy/units/quantity.py:481: RuntimeWarning: overflow encountered in exp
  result= super().__array_ufunc__(function, method, *arrays, **kwargs)
<Quantity [        inf,         inf,         inf, ..., 12.95190431,
           12.92977839, 12.90771972]>

(here it gives some infinities near low wavelengths, but you can avoid this by clipping to a larger lower bound).