3

I am a python newbie. I am trying to evaluate Planck equation using python. I have written a simple program. But when I give input to it is giving me an error. can anyone hep me where I am going wrong? Here is the program and the error:

Program:

from __future__ import division
from sympy.physics.units import *
from math import *
import numpy
from scipy.interpolate import interp1d

#Planck's Law evaluation at a single wavelength and temperature
def planck_law(wavelength,temperature):
    T=temperature
    f=c/wavelength
    h=planck
    k=boltzmann
    U=2*h/(c**3)*(f**3)/(exp(h*f/(k*T))-1)
    return U.evalf()

Input: I have imported the function as 'cp' and the input is as follows

value = (cp.planck_law(400,2000))

Error:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>`enter code here`
  File "Camera_performance.py", line 14, in planck_law
  U=2*h/(c**3)*(f**3)/(exp(h*f/(k*T))-1)
  File "/usr/lib/python2.7/dist-packages/sympy/core/expr.py", line 221, in __float__
  raise TypeError("can't convert expression to float")

TypeError: can't convert expression to float
Pascal Bugnion
  • 4,878
  • 1
  • 24
  • 29
Niroop
  • 53
  • 1
  • 1
  • 4
  • 1
    for a start: do not use sympy, since it is (very) slow, and uses its own datatypes. from scipy.constants you can get your constants as well, and they are numpy floats. Another thing is to use numpy.exp instead of math.exp. – usethedeathstar Mar 21 '14 at 09:50

3 Answers3

10

It seem like you are mixing namespaces, since your are using from ... import *. You wanted to use sympy.exp() but your code uses math.exp(). It is good practice to keep the namespaces separated, i.e. never use from ... import * - it might seem like more typing at first, but will produce much cleaner easier comprehensible code in the end. Try:

import sympy as sy
import sympy.physics.units as units

def planck_law(wavelength,temperature):
     """Planck's Law evaluation at a single wavelength and temperature """   
     T=temperature
     f=units.c/wavelength
     h=units.planck
     k=units.boltzmann
     U=2*h/(units.c**3)*(f**3)/(sy.exp(h*f/(k*T))-1)
     return U.evalf()

# Test:
print(planck_law(640e-9*units.m, 500*units.K))
# Result: 1.503553603007e-34*kg/(m*s)
Dietrich
  • 5,241
  • 3
  • 24
  • 36
1

Your term h*f/(k*T) isn't unitless, so you cannot pass it to exp(). Doesn't make sense in a physical context, either ;-)

You can get a result if you divide it by K and m:

exp(h*f/(k*T)/K/m)

but that surely makes no sense. It's just to make the program run (into nonsense).

I guess you will have to check your formula and find out how you really want to compute the value you want to pass to exp().

EDIT:

As Pascal pointed out, you are just missing the units for the arguments you pass. Try it with:

planck_law(400e-9*m,2000*K)

Returns:

3.20224927538564e-22*kg/(m*s)
Alfe
  • 56,346
  • 20
  • 107
  • 159
1

When you call your function, you need to pass units to the temperature and wavelength argument. Calling cp.planck_law(400*meters,2000*K) gives the expected 1.15133857387385e-33*kg/(m*s).

The problem arose with the exponential function, which, rightfully, expects to receive a dimensionless argument. Since the arguments you were passing did not include units, Sympy treated them as dimensionless floats, rather than a temperature and a length, as required.

Pascal Bugnion
  • 4,878
  • 1
  • 24
  • 29