2
$\begingroup$

I have a data set of intensities and corresponding wavelengths. I want to fit it to the Planck black body radiation curve so that I can estimate the temperature for the given data . I have used the filter while getting these data, So, it will be good to use data from 400 nm to 800 nm to block the effect of filter.

Here is small sample of the data.

360.84  7.3516E-01
361.06  7.6002E-01
361.27  7.7199E-01
361.49  7.7312E-01
361.70  7.5713E-01
361.92  7.5512E-01 

where 1st column gives wavelength in nm and the 2nd columns give the corresponding intensities calculated from the ocean optics usb4000 spectrometer. The full data can be get from (https://pastebin.ca/3956256). I need to plot with Planck black body radiation to estimate the temperature related to the given spectra. Please only take the data in the range from 400nm to 800nm.

data = Import["link of my data in the text file of my desktop", "Table"];
pts = data[[#]] & /@ Table[i, {i, 1, 1458, 50}];
h = 6.62607004*10^(-34);
c = 299792458;
kb = 1.38064852*10^(-23);

Planks[l_, T_, A_] := (2*h*c^2 10^45)/(A l^5 (Exp[(10^9 h*c)/(l*kb*T)] - 1));
fittesana2 = FindFit[data, Planks[l, T, A], {{A, 1*10^8}, {T, what?}}, l];
Plot[Planks[l, T, A] /. fittesana2, {l, 400, 900}, 
  PlotStyle -> Red, 
  PlotRange -> All, 
  Epilog :> {Blue, PointSize[0.015], Point[pts]}]

Pfit = NonlinearModelFit[data, Planks[l, T, A], {{A, 1*10^8}, {T, what?}}, l]
Plot[Pfit[l], {l, 400, 900}, 
  PlotStyle -> Red, 
  PlotRange -> All, 
  Epilog :> {Blue, PointSize[0.015], Point[pts]}]

Normal[Pfit]

Pfit["ANOVATable"]

Pfit["ParameterTable"]

Pfit["FitCurvatureTable"]

The above code does not work. Please give me some suggestions

$\endgroup$
  • $\begingroup$ It is unfortunately not possible to say anything about this, when you don't provide a data-set that can be fitted. The 6 sample points are not enough. What exactly does not work? $\endgroup$ – halirutan Jan 7 '18 at 0:40
  • $\begingroup$ We really need the full data set. You might use Pastebin to share it. Alternatively, it appears that you are extracting 30 points from your full data set. Give us those points in the question and modify your code (i.e., pts = {list of 30, ...}). $\endgroup$ – Jack LaVigne Jan 7 '18 at 0:54
  • $\begingroup$ You are right, i need to give full data. The data is very much big, So i edited the question with the inclusion of the full data shared with the Pastebin. Thank you so much. $\endgroup$ – rabink Jan 7 '18 at 2:09
5
$\begingroup$

I have played with your data and you were not far from having this running. First I plotted the data you offered and noticed there were some (automatic?) change of scale at high wavelengths. Also, I had to edit the file by hand to make it readable by Mathematica (1. I changed the Enn notation to *10^nn 2. I added a space to separate the x and y coordinates at the end of the file to make them readable by the Import 3. the Drop[ ,-1] is to remove the last entry which was empty: End-of-file?). I used:

data1 = Drop[Import[NotebookDirectory[] <> "Spectral file - Black 
Body.txt", "Data"], -1]

data$raw = {ToExpression[#[[1]]], ToExpression[#[[2]]]} & /@ data1

enter image description here

So I have trimmed the data$raw file with:

data = Select[data$raw, 500 < #[[1]] < 850 &]

I changed your pts definition using Range (It's more Mathematica like!) but I don't use it.

pts = data[[#]] & /@ Range[1, Length[data], 50];
h = 6.62607004*10^(-34);
c = 299792458;
kb = 1.38064852*10^(-23);

Planks[l_, T_, A_] := (2*h*c^2 A)/(l^5 (Exp[(10^9 h*c)/(l*kb*T)] - 1));
Plot[Planks[l, T, A] /. {A -> 10^36, T -> 1700}, {l, 400, 900}]

fittesana2 = FindFit[data, Planks[l, T, A] + bck, {{A, 10^36}, 
{T, 1700}, {bck, 0.1}}, l]

Plot[(Planks[l, T, A] + bck) /. fittesana2, {l, 400, 900}, PlotStyle -> Red, 
PlotRange -> All, Epilog -> {Blue, PointSize[0.015], Point[data]}]


Pfit = NonlinearModelFit[data, Planks[l, T, A] + bck, {{A, 10^36}, 
{T, 1700}, {bck, 0.1}}, l]

Plot[Pfit[l], {l, 400, 900}, PlotStyle -> Red, PlotRange -> All, 
Epilog -> {Blue, PointSize[0.015], Point[data]}]

Note that I rewrote your equation to have the parameter A at the numerator (to control scale fitting from the denominator is an added challenge). I got the following curve with Planks[]

enter image description here

I find that A set at 10^36 gives you a curve of the same order of magnitude than the experimental one.

For the fitting, I have added a small background because your experimental curve does not go to zero and I got this (I got the same with Pfit). The fitting complains about accuracy not reached but at least it is a start for you! I got the following fit parameters: {A -> 1.13599*10^37, T -> 1174.69, bck -> 0.217298}

enter image description here

$\endgroup$
  • $\begingroup$ Hi, Gwanguy thank you so much for your idea that i could start off. I had some questions regarding it. (1) while trimming the real data, i used {data = Select[pastebin.ca/3956256, 500 < #[[1]] < 850 &], but it doesn't trimmed. My raw data is in the txt format. (2) While plotting, {Plot[Planks[l, T, A] /. {A -> 10^36, T -> 1700}, {l, 400, 900}]}, I assumed the value T=1700. But actually i need to calculate the value of T with the given data set fitting to Planck. Is there a way to calculate the value of T with the given data set provided.Thanks $\endgroup$ – rabink Jan 7 '18 at 14:11
  • $\begingroup$ When I copied and paste your data, I got a RTF file, so I had to convert it to PlainText and converted it with ToExpression to extract the numbers. Now I don't know why the Select would not work if it is what you did. About the temperature, I think once A and T give you a reasonable amplitude, you can start fitting. I don't see another way of getting the temperature. $\endgroup$ – Gwanguy Jan 7 '18 at 18:29
3
$\begingroup$

I downloaded the data and used your constants and equations

h = 6.62607004*10^(-34);
c = 299792458;
kb = 1.38064852*10^(-23);

plancks[l_, t_, a_] := (2*h*c^2 10^45)/(a l^5 (Exp[(10^9 h*c)/(l*kb*t)] - 1))

(I used small case letters, a good rule to avoid conflicts with Mathematica's symbols).

Here is your equation in TraditionalForm. You are using 10^9 to convert the wavelength to meters and the parameter a to adjust the amplitude of the measured spectrum. You are trying to infer the temperature and the amplitude adjustment parameter.

Mathematica graphics

Here is a plot of black body radiation at your starting temperature of 1700 degrees kelvin.

Plot[
 plancks[λ, 1700, 1 ],
 {λ, 0 , 5000},
 PlotStyle -> Black
 ]

Mathematica graphics

and here is your data

ListLinePlot[data, PlotRange -> All, PlotStyle -> Black]

Mathematica graphics

In order to more readily compare them, let's limit the x-axis scale to the maximum wavelength of your data.

{lmin, lmax} = MinMax[data[[All, 1]]]
(* {345.55, 1040.56} *)

{ampMin, ampMax} = MinMax[data[[All, 2]]]
(* {0., 137.43} *)

Plot[
 plancks[λ, 1700, 1 ],
 {λ, lmin, lmax},
 PlotStyle -> Black
 ]

Mathematica graphics

At this point I would make the statement that Mathematica is working fine in the sense that it finds the best fit parameters for your data (demonstrated below).

The real problem is that your data and model don't match well.

Is it intended that you only fit out to the first upturn or maybe only the last upturn?

At any rate, with regard to estimating the parameter a and the temperature, I assumed that you were only seeing a portion of the upturn and estimated the amplitude adjustment parameter to be approximately 2*10^8.

You need to give the software some help and constrain a to a reasonable value about that estimate (this is important using the raw data).

nlm = NonlinearModelFit[
  data,
  {
   plancks[l, t, a],
   1*10^8 < a < 1*10^9
   },
  {
   {t, 1700},
   {a, 2*10^8}
   },
  l
  ]

nlm["BestFitParameters"]
(* {t -> 1352.47, a -> 1.97905*10^8} *)

It does the best job possible. Here are two figures.

Plot[
 nlm[l],
 {l, 340, 1050},
 PlotStyle -> Red,
 PlotRange -> {{lmin, lmax}, {ampMin, ampMax}},
 Epilog :> {Blue, PointSize[0.015], Point[data]}
 ]

Mathematica graphics

Plot[
 nlm[l],
 {l, 340, 1050},
 PlotStyle -> Red,
 PlotRange -> All,
 Epilog :> {Blue, PointSize[0.01], Point[data]}
 ]

Mathematica graphics

Again, the problem is with your data. Looks like there is some junk at the start (which you might discard) and something that might be gain changes in the measured spectrum at longer wavelengths.

Update

To limit the data from say 500 to 880 use

data2 = Select[data, 500 < #[[1]] < 880 &];

Then use the same procedure replacing data with data2.

nlm = NonlinearModelFit[
  data2,
  {
   plancks[l, t, a],
   1*10^8 < a < 1*10^9
   },
  {
   {t, 1700},
   {a, 2*10^8}
   },
  l
  ]

nlm["BestFitParameters"]
(* {t -> 1284.47, a -> 1.98349*10^8} *)

Plot[
 nlm[l],
 {l, 500, 880},
 PlotStyle -> Red,
 PlotRange -> All,
 Epilog :> {Blue, PointSize[0.01], Point[data2]}
 ]

Mathematica graphics

This looks slightly better but there still appear to be problems.

$\endgroup$
  • $\begingroup$ thanks Jack, I put the temperature T=1700 to plot the data. Actually i used filter to get the data. So, the data from 400 to (800 or 880) nm is good to fit it to planck radiation. My problem is , i need to estimate the temperature T with the given set of data when it fitted to Planck. How can i do it? $\endgroup$ – rabink Jan 7 '18 at 16:59
  • 1
    $\begingroup$ This is a much trickier problem than you think. What do you mean by "intensity"? Power per wavelength interval? Power per frequency interval? Photocurrent in some bandwidth? If you're using a filter, the filter's detailed transmission has to go into your model (no filter has a perfectly clean bandpass). Does your photosensor have dark current, and can you correct for it or model it? $\endgroup$ – John Doty Jan 7 '18 at 17:07
  • $\begingroup$ The temperature is being estimated (1284 K). Using the trimmed data (500 to 880) appears ,more reasonable but there are still problems. At any rate, this is a site about Mathematica. You are responsible for making decisions about models and data. We help with how to use the software. $\endgroup$ – Jack LaVigne Jan 7 '18 at 17:15
  • 2
    $\begingroup$ I think the strange peaks may be Wood's anomalies. This cheap and terrible spectrometer has a very nonuniform responsivity, and in fact Ocean Optics sells calibration sources for radiometric applications. Attempting to make anything of this data without a radiometric calibration is a fool's errand, I'm afraid. $\endgroup$ – Oleksandr R. Jan 7 '18 at 18:03
  • 1
    $\begingroup$ @rabink Mathematical solvers that fit parameters to data can be classified as global or local. FindFit and nlmModelFit are local solvers. Typically they need a reasonable starting value or they can get stuck in what is known as a local minimum. In the second example above, the starting value for t and a are 1700 and 2*10^8 respectively. nlmModelFit is simultaneously determining the best fit parameters for t and a which it found to be 1284 and 1.98*10^8 respectively. $\endgroup$ – Jack LaVigne Jan 7 '18 at 20:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.