High accuracy needed

Share icon

Hi everybody,


I need to calculate some numbers based on Planck's Law. The formula can be found e.g. under http://en.wikipedia.org/wiki/Planck%27s_law. When I wirte the formula as

function PlancksLaw(lambda, T)
variable lambda // in m
variable T // in K

variable h = 6.62606896E-34 // J*s
variable c = 299792548 // m/s
variable k = 1.3806504 // J/k

return 2 * pi * h * c^2 / lambda^5 * 1/(exp(h*c/(lambda*k*T)) -1)
end

and I call it with

print PlancksLaw(550E-9, 5000)

, i.e. for green visible light and a D50 photographic standard lamp, it yields always inf. And that is obviously not true, because blackbodies do not radiate an infinite amount of energy. The problem can be brought down to the fact, that the formula calculates the exponential function of 5.23191e-23, which returns 1. That is obviously wrong, because it should be something around 1.00000000000000000000005. In many cases this minor difference does not play any important role, but here it does.

So the question is, how can I turn the Igor precision high enough to make this piece of phyical calculation work.


Many thanks and best regards,
Joerg Kunze.




Although I never used it, the operation APMath could help you here. It allows "arbitrary precision calculation of basic mathematical expressions."

Hope that helps,
Thomas

February 14, 2011 at 06:58 am - Permalink

joerg.kunze wrote:
I need to calculate some numbers based on Planck's Law. The formula can be found e.g. under http://en.wikipedia.org/wiki/Planck%27s_law. ... The problem can be brought down to the fact, that the formula calculates the exponential function of 5.23191e-23, which returns 1. That is obviously wrong, because it should be something around 1.00000000000000000000005. In many cases this minor difference does not play any important role, but here it does.

So the question is, how can I turn the Igor precision high enough to make this piece of phyical calculation work.


Independently of the method chosen to improve precision, the answer 1.00000000000000000000005 is fundamentally outside the bounds of precision that you can reliably determine based on the input precision of your constants. The best you can do is 1.00000000.

Also, if I am not mistaken, k has a value closer to 1.38E-28 J/K.

BTW, you might be interested in the Planck Distribution Plot Package http://www.igorexchange.com/project/PlanckDistributionDemo.

--
J. J. Weimer
Chemistry / Chemical & Materials Engineering, UAHuntsville

February 14, 2011 at 09:24 am - Permalink

Igor variables are stored in double-precision IEEE floating point format which can represent about 15 decimal digits of precision.

Therefore, 1.00000000000000000000005 can not be represented in double-precision floating point and evaluates to 1.0.

As Thomas said, you might be able to use APMath. Execute:
DisplayHelpTopic "APMath"

February 14, 2011 at 09:51 am - Permalink

Hello Joerg,

A fundamental step in numerical computation is scaling your equations. In as much as it is nice to work in SI and keep your I/O in original units this makes no sense when your computer works in double precision. So, start by re-writing hc/k. BTW: your K above is missing a crucial e-23...

Variable myConst=hc/k= 0.013594376153594717388

As you can see, despite the large initial exponents, this value if manageable. Now you can compute exp(myConst/(lambda*temp)

Furthermore, since your temperatures are on the order of 1000 you should use that to scale your wavelength input to the micron range.

I hope this helps,

A.G.
WaveMetrics, Inc.

February 14, 2011 at 11:17 am - Permalink

Many thanks. The E-28 makes indeed a lot of difference!

Sometimes it is really hard to recognize tho pretty obvious things...

Best regards,
Joerg.


February 16, 2011 at 08:06 am - Permalink