Code covered by the BSD License

### Highlights from Abel Inversion Algorithm

5.0
5.0 | 9 ratings Rate this file 28 Downloads (last 30 days) File Size: 4.17 KB File ID: #43639 Version: 1.5

# Abel Inversion Algorithm

### Carsten Killer (view profile)

26 Sep 2013 (Updated )

Fourier-based reconstruction of an unknown radial distribution assuming cylindrical symmetry.

File Information
Description

The reconstruction of the radial density distribution of a cylindrically symmetric object is a common task in different area of physics (e.g. plasma physics). Typically, optical measurements of objects like plasma columns or flames are integrated along the line of sight.
To obtain the underlying distribution from a measured projection, the inverse Abel transform has to be calculated.
This submission provides a Fourier-based algorithm which extracts the radial (2D) distribution from a one-dimensional projection measurement. The algorithm has been proposed and published by G. Pretzler (Z. Naturforsch. 46a, 639 (1991)). Compared to earlier approaches towards Abel inversion, this algorithm is relatively insensitive to noise and errors in the determination of the object's center (see G. Pretzler et al. , Z. Naturforsch. 47a, 955 (1994)).
The fundamental idea is to fit the whole measured profile to a set of cos-expansion-based integrals. (In the conventional approach, in contrast, the radial distribution is obtained by starting at the edges and incrementally iterating towards the center - making it more prone to errors.)

Required Products MATLAB
MATLAB release MATLAB 8.0 (R2012b)
MATLAB Search Path
`/`
05 Apr 2016 Wenting He

### Wenting He (view profile)

Thanks a lot! Because my data at the edge is fluctuant, so I cannot so I cannot subtract the background by "h2 = h - h(end)". But I'll try to reduce the effect caused by my data at the edge.
It's really nice to discuss with you!

Comment only
05 Apr 2016 Carsten Killer

### Carsten Killer (view profile)

20 data points is not that much, so the negative result in the center might be an unavoidable result of your limited spatial resolution. Improving the number of data points to 950 by interpolation/fitting does not really add new information, so I am not surprised that this does not change much.

Regarding substracting the background: What I had in my mind does not change the number of data points. If your data is in a 20-element vector "h", just create a new data vector

"h2 = h - h(end)"

which will still contain 20 data points but goes to zero at the edge.

Besides this suggestion, I am afraid that you have to accept the negative results in the center which are probably caused by the spatial resolution of your dataset

Comment only
04 Apr 2016 Wenting He

### Wenting He (view profile)

Thank you very very much Carsten for the patient and professional respond.
I have 20 data points for one Abel inversion. And If I set Nu as 6,then I got negative reconstructed value in the center but if I Nu is 4, then I won't get negative value.
I also tried to subtract the background, then I only have 16 or 13 data points. but it didn't help too much. I still got negative data even I set Nu as 4.Only that I set Nu smaller than 4.
what's more, I also tried to fit the data first and then subtract the background. then I have about 950 data points. however, even in this case, I can obtain positive value only that I set Nu smaller than 6. Therefore, I think Nu have a significant influence on the value in the center of plasma plume but not the number of my data points. Anyway, I don't need very accurate reconstructed value at the edge because the emission at the edge is naturally low.
in my case, do you still have better suggestions?

Comment only
04 Apr 2016 Carsten Killer

### Carsten Killer (view profile)

How many data points do you have? In any case, Nu should be considerably smaller than the number of data points.

If the data does not fall to zero at the edge you will get wrong reconstruction results at the edge of your plasma since the basic geometry of Abel inversion requires this boundary condition. Whether this error is restricted to the edge region of your data set or whether it might also affect the plasma center depends, again, on the size of your data set and Nu.

As a test, you could manually substract the background light from your dataset, so that the data at the edge becomes zero. Then you can check whether you still get negative emissivities in the center

Comment only
01 Apr 2016 Wenting He

### Wenting He (view profile)

Morning Carsten, thank you for your reply. I do get a negative value in the center of plasma column once I set the Nu larger than 5. And if I set Nu as 20, then I will get this warning "Matrix is close to singular or badly scaled. Results may be inaccurate." Maybe I just have a few data. Do you think it is a problem that my data doesn't fall to zero at the edge of plasma column?

Comment only
31 Mar 2016 Carsten Killer

### Carsten Killer (view profile)

Dear Wenting He,
the question for the best choice for Nu is indeed not easy to answer. It depends on your problem, e.g. on the number of data points available. While in general a larger Nu allows for more details to be reconstructed you also emphasize noise artifacts and other uncertainties (e.g. a slightly incorrect determination of the radial center of your plasma)

Of course you are right that a negative plasma emissivity is unphysical. Do you still get negative results for larger Nu (e.g. 10,15,20)? Can you rule out that you have other processes influencing your measurement such as re-absorption of light? (such things might also lead to the appearance of negative emissivities).

Comment only
30 Mar 2016 Wenting He

### Wenting He (view profile)

Hello Carsten, Thanks for sharing your code. It helps a lot. But I still have a question that if I set Nu as 6, I will obtain a negative reconstructed distribution in the center of plasma column. Because my measured profile have a low value in the center of plasma column. I don't think this is correct because the smallest emissivity should be equal or larger than zero. what's the reason for this? PS: if I set Nu as 4 or smaller than 4, the reconstructed distribution won't be negative. But in your explanation, 4 is too low which causes low=pass filtering effect. It looks that the result is highly dependent on the Nu. So how do I know which is the best Nu for my data?

Comment only
16 Mar 2016 Carsten Killer

### Carsten Killer (view profile)

Hi,
"integral" was introduced in Matlab2012. If you use an earlier version of Matlab, just replace "integral" with "quadgk", which should also work fine.

Comment only
11 Mar 2016 sc sc

### sc sc (view profile)

hi,
some functions used in your algorism seems cannot be found, such as 'integral' in compute_expansion.m . Did I should download some other toolbox?

Comment only
29 Oct 2015 Carsten Killer

### Carsten Killer (view profile)

@ot tl: B and L are defined according to equation (8) in the original article by Pretzler, which is in the form A*L=B.
Hence, in my code L is the matrix containing the product of each integral h_n with each other integral h_n. Further, B is defined as the product of the measured data with the integral h_n. On both sides of the equation, we sum up over all positions (index k).

Comment only
28 Oct 2015 ot tl

### ot tl (view profile)

Hello Carsten, Thanks for sharing your code. I'm trying to go over your code in order to learn how to implement a theoretical mathematical algorithm in MATLAB. However I did not manage to understand the code regarding the solution of the equation system: A*L=B. Can you please explain why did you define the L matrix and B vector the way you did ? ( in particular: what's the meaning of the sums in these definitions ? ) . Thanks a lot !!!

Comment only
18 Jul 2015 Carsten Killer

### Carsten Killer (view profile)

@F. Salehi: Actually the only reason for using this particular basis function is that it was used in the original article by Pretzler (see description), which inspired this code.
You are right, it should principally also work for simpler basis functions. I have not investigated whether the choice of the basis function has any practical impact on the algorithm.

Comment only
16 Jul 2015 F.Salehi

### F.Salehi (view profile)

Hi Carsten, Thanks for sharing your scripts here.
I was wondering if there is any specific reason for using (1 - (-1)^n*cos(n*pi*x/R)) as your basis function for expansion instead of simply using cos(n*pi*x/R) ?
This way your basis will be orthogonal.

Comment only
06 Apr 2015 dj327

### dj327 (view profile)

Hi Carsten, the problem of isqcuivefit not running is solved. A classic 1*N matrix or N*1 matrix problem. Although I made sure the input of h is 1*N matrix as requested, dimension and h and myfun inside lsqcurvefit(myfun,x0,hn,h) don't match...so I just simply transposed h

Comment only
02 Apr 2015 dj327

### dj327 (view profile)

02 Apr 2015 Carsten Killer

### Carsten Killer (view profile)

@dj327: That depends on whether the line of sight is a full secant line (as in most experimental situations) or just a half secant line, which would generally be sufficient for the problem. The code is written for the latter situation (with hindsight maybe not a good choice, but I do not want to change the submitted code for this minor issue for reasons of consistency).

Please note that this issue is already addressed in the last paragraph in the code description above ("UPDATE: ...").

Comment only
02 Apr 2015 dj327

### dj327 (view profile)

I believe you have missed a factor of 2. I need to divide f_rec by 2 to get the expected result.

Comment only
06 Nov 2014 Carsten Killer

### Carsten Killer (view profile)

@aneps: Just process every row (or column, depending on the symmetry axis in your system) separatey. For example, you can call the abel_inversion function 100 times for the 100 rows/colmuns, using a simple loop.

Comment only
06 Nov 2014 aneps

### aneps (view profile)

My data input (H) is a colormap; a matrix of size 100 x 100. When I tried to use this program to perform inverse Abel transform it shows error : "Error using .*
Matrix dimensions must agree" .

I understand the input data set should be 1 X N matrix. How can I use this program in my case where I have an intensity map of mxn matrix (mine is 100x100)? Thank you.

Comment only
20 Oct 2014 surya Narayan

### surya Narayan (view profile)

Whenever i mask the data i have does not start from the axis symmetric center but at a position of some pixels off from the center. How to use abel inversion in this case? Please suggest some references.

Comment only
20 Oct 2014 Carsten Killer

### Carsten Killer (view profile)

Could you please specify what exactly you mean with "masked" ? (In any case, your data has to be smooth enough to comply with the general constraints of the Abel problem.)

Comment only
18 Oct 2014 surya Narayan

### surya Narayan (view profile)

I am using the same code for abel inversion of the axis symmetric data. I have a masked region at the center. Abel inversion is giving me a wrong values if i have a mask. What should i do to make the code work? please send me some references if u have any.

18 Oct 2014 surya Narayan

### surya Narayan (view profile)

12 Oct 2014 Carsten Killer

### Carsten Killer (view profile)

Hi Laurel,
if the symmetry axis is in the center of your image, the image has to be "cut in half" for Abel processing, i.e. the left and right half are calculated independently. Please remember to use something like fliplr for the left half, so that the first element of the input vector corresponds to the center of the original image.
In your case, an input vector therefore seems to contain only 64 data points. That is not much - I am not sure how well the algorithm will perform in this case.
If are still not sure whether you use the function correctly , you can send a sample image to killer@physik.uni-greifswald.de and I will take a look.
Cheers, Carsten

Comment only
12 Oct 2014 Laurel

### Laurel (view profile)

Hello,
I'm trying to perform an Abel inversion on a set photographs (128x800) of a flame issuing from a jet (to analyze for emissions). The cylindrical symmetry is oriented along the longer (800) axis. Per previous comments, I've been trying run the Abel function one row at a time (running from 0 to 128). However, the produced results do not appear to be accurate. I also tried cutting the image in half and started my iteration from the centerline of the jet. I am unsure if I am using the function in the correct manner.
email: lpaxton@usc.edu
Thank you

Comment only
14 Aug 2014 Carsten Killer

### Carsten Killer (view profile)

Recently, the journal (Zeitschrift fÃ¼r Naturforschung) has been digitalized. The paper on Abel inversion can be found on page 639 at

http://www.znaturforsch.com/aa/v46a/c46a.htm

1. Actually, noise is accounted for by using a finite number of cos-expansions. The less expansions you use, the more noise is filtered out (real physical information is of course neglected on the same scale as the noise, so you shouldn't use a too little number of expansion.)

2. You are right, center shift is not included in this code.

Comment only
14 Aug 2014 nudt

### nudt (view profile)

The programm is helpful, but several factors such as noise and center shift are not taken into account.
Could you please share reference 1 to me so that I can understand the theoretical background? my email is jianguotan@hotmail.com
Thank you

21 Jul 2014 Carsten Killer

### Carsten Killer (view profile)

It is not necessary to figure out a function. Just use a row of your image (with the symmetry axis at the first element of the vector) as the first input of the main function ("h").

(the code works numerically (not analytically), e.g. the polynomial sample function in GENERATE_TEST_DATA is also discretized onto a grid and is not treated as an analytical function)

Comment only
18 Jul 2014 Dimas Avila

### Dimas Avila (view profile)

Thanks for the speedy reply Carsten! So for each row of my image should I figure out a function that characterizes that row and input that?

Comment only
18 Jul 2014 Carsten Killer

### Carsten Killer (view profile)

@Dimas Avila: The polynomial function is just an example given in the function GENERATE_TEST_DATA. (the result of this function is shown in the sample image at the top of this page). This function is only called if you do not provide input data.

The data vector that you want to perform the Abel inversion on is the first input variable in the main function. Please perform the Abel inversion for each row/column (depending on whether your symmertry axis is vertical or horizontal) of your image separately.

Comment only
17 Jul 2014 Dimas Avila

### Dimas Avila (view profile)

I was trying to perform an Abel inversion on a set of flame photographs I took to look for emission radicals. I was wondering how I can use this program to perform this. I see there is a function which generates a polynomial distribution function and then takes the abel inversion of that function. I just don't know how to do this for an image.

Comment only
08 Jul 2014 Carsten Killer

### Carsten Killer (view profile)

@Adnan: I don't see any way in which Abel inversion would be related to your problem. You might want to consider some edge detection algorithm (Matlab offers plenty of them) to find the shape of the person in each frame.

Comment only

@Carsten. In my case i have video of a person who is walking in front of camera (we can see the side view of person) and there is nothing at the background. so all the frames consist of each step he took for walking. so i want to apply Abel inversion projection on each frames (similar like making features vector)

Comment only
07 Jul 2014 Carsten Killer

### Carsten Killer (view profile)

@Adnan: I am not sure what you mean with "finding changes". Of course, in different measurement images, the Abel inversion results of these images will probably differ as well. But please consider that the Abel inversion is only suited for a special situation, i.e. the investigation of symmetrical objects. Of course you can apply the Abel inversion code individually for each row or column, if your data is some kind of line-of-sight integrated measurement of a symmetrical object/system.
If the axis of symmetry is vertical, you have to consider the rows of your images. Please keep in mind that the axis of symmetry is at the start of your input vector (row) , i.e. the first element of the row vector corresponds to the center of your system r=0 and the last element corresponds the the edge, r=R.

Comment only

@ Carsten and Borek, Thanks for your answers. Actually i want to use Abel transform on sequence of images i have to find the changes in that ... so i want to in that case is it possible to apply Abel inversion directly on 2D row wise and column wise..??

Comment only
07 Jul 2014 Borek

### Borek (view profile)

As Carsten wrote, result is a distribution
f(r), where r is the diameter of a symetric circle. It means, the distribution is same for angle phi from 0 to 2pi. To plot f(r,phi) in 2D, the same parametrization can be used. If you need something more sofisticated, see assymetric abel inversion, which can be used for assymetric signals for example for simple plot of irradiation profile of plasma in tokamaks.

Comment only
07 Jul 2014 Carsten Killer

### Carsten Killer (view profile)

@Adnan Farooq: Could you maybe pose a more specific question? To summarize, this code takes your input F(y) and calculates the inverse Abel transform, resulting in the corresponding f(r). (using the standard notation, e.g. used in http://en.wikipedia.org/wiki/Abel_transform ).

Please note that this code only performs the calculation for one vector, i.e. one row (or column, depending on your systems's geometry) of your image. If you want to transform a complete 2D image, every row (or column...) has to be processed individually.

Comment only

Hello Carsten

Thank you for submitting for code. I am studying about Abel Transform and want to use for 2D image.. i want to know that does this code supports?

10 Jun 2014 Carsten Killer

### Carsten Killer (view profile)

@surya Narayan:
1) In your case, only half of the data has to be used, i.e. the symmentry axis is at the beginning of your data vector. From the geometrical point of view, both halves (left and right) are independent. For perfect symmetry, they will of course yield identical results. Therefore, the difference between the left and right half isa good indicator for the quality of symmetry in your system.

2) The radius is important for the scaling of the results. For any arbitrary value of R you will get a qualitatively correct result from the Abel inversion (except for very low values of R << 1 , which should be avoided). But, if you need a quantitatively correct result, you have to use the (outer) radius of your object. (e.g. if you investigated some glowing or light absorbing stuff in a glass cylinder which has a radius of 50 mm , just use R=50 and your result will be scaled to mm units)

Comment only
10 Jun 2014 surya Narayan

### surya Narayan (view profile)

I have 2-D unwrapped phase values (axis is at the middle of the data)

1. I am entering the axis symmetric data row wise(vertical symmetry is present) as full row. Is this right or should i give only one half of the data?

2. What is the radius that has to be given as input?

can anyone help me out?

23 May 2014 Sina Tebianian

### Sina Tebianian (view profile)

22 May 2014 Sina Tebianian

### Sina Tebianian (view profile)

Thanks Carsten, Sorry I wrote by mistake "arrow" instead of array. Anyway my data has a physical range that goes from 0.48 to 1. So I understand that you suggest that I subtract the minimum value of the vector from all its components. Once I run the function, I will add the minimum value to the output data again. I figured the normalization once I went through the mfile. Many thanks for your help.
Sina

Comment only
22 May 2014 Carsten Killer

### Carsten Killer (view profile)

Hi Sina,
about the data range: If you use the plot option of my code, both your original data and the calculated radial profile will be plotted normalized to 1 (to make it easer to visually compare the profiles). If you are interested in the correctly scaled data, please use the output (f_rec=abel_inversion(...) ) and plot it manually. (of course, the code uses your original data for the calculations and not the normalized data shown in the figure).

About your data: I am not sure what you mean with an "arrow". However, the boundary value of your data vector (i.e. the last value) should ideally be zero (or something close). This condition is crucial for the algorithm and is imposed by the basic geometry of the situation. If there is some kind of constant background information behind your data, you might want to consider substracting 0.5 in order to create a data vector from 0 to 0.1 (it depends on your specific situation whether this is legitimate).
I hope this helps a bit. Feel free to ask further questions if anything remains unclear.
Carsten

Comment only
22 May 2014 Sina Tebianian

### Sina Tebianian (view profile)

Hi Carsten,
Thanks for your useful code. I tried to run the program with an arrow of my data but the plots produce data that are completely unexpected. The measured profile that must reflect my data seems to be shifted up. My data range goes from 0.51 to 0.6 but then on the profile I see that the measured profile goes from 0.8 to 1. Also the reconstructed profile is completely different from the measured one. I am wondering if I have not understood part of the code. My vector is basically 1x301 that goes from 0.5 to 0.6 therefore the radial profile should not be very different from these values.
Sina

Comment only
22 May 2014 Mohammed Ibrahim

### Mohammed Ibrahim (view profile)

@Carsten Killer

I am not able to access reference 1, based on which the program is created. It would be great if you can share it with me.
ibrahim@aero.iisc.ernet.in
thanks

Comment only
15 May 2014 Carsten Killer

### Carsten Killer (view profile)

@Mohammed Ibrahim: The "measured distribution" is F(y) and the "reconstructed distribution" is f(r).

(I am sorry for any confusion caused. I thought these denotation for F(y) and f(r) to be practical description for applications in data analysis of physics experiments)

Comment only
15 May 2014 Mohammed Ibrahim

### Mohammed Ibrahim (view profile)

@carsten kiler

when i run the abel inversion for a test case i get measured and reconstructed distribution F(y). How do i get the radial distribution f(r)?

Comment only
15 Feb 2014 Erzberger

### Erzberger (view profile)

24 Jan 2014 Carsten Killer

### Carsten Killer (view profile)

@Akash:
By default, lsqcurvefit is not called unless you set the fifth input parameter to 1. If you try the same dataset without using lsqcurvefit (e.g. setting the fifth input to 0), do you also get an error?
Furthermore, do you also get this error when using the lsqcurvefit routine with the built-in test function (i.e. calling "abel_inversion([],[],10,1,1) ")?

Comment only
23 Jan 2014 Akash

### Akash (view profile)

Dear Sir
Thanks, Now i am suffering from another kind of problems like

??? Error using ==> lsqcurvefit at 253
Function value and YDATA sizes are
incommensurate.

Error in ==> solve_lsq at 17
A=lsqcurvefit(myfun,x0,hn,h);

Error in ==> abel_inversion at 79
A=solve_lsq(h,hn,x0); % solve for
amplitudes A

Comment only
22 Jan 2014 Carsten Killer

### Carsten Killer (view profile)

Hello Akash,
I recently found out that the integral function was introduced in Matlab2012 and is not available in older versions.
However, simply replacing "integral" by "quadgk" should do the job in older Matlab versions.

Comment only
22 Jan 2014 Akash

### Akash (view profile)

Hello Sir
I am not able to run this program. Their is an error
" Undefined function or method 'integral' for input arguments of type 'function_handle"
I think integral function is missing. Since i am new to matlab so not 100% sure what is going on. Please help me.

Comment only
29 Nov 2013 Carsten Killer

### Carsten Killer (view profile)

@R: I will try it with Octave as soon as possible. With Matlab, your example works for me: If I do this
x=0:0.01:6;
h=2*sqrt(6^2-x.^2);
and then start the function:
abel_inversion(h,6,10)
The result is in my case a more or less constant line at y=2 for 0<x<6 (not minding the numerical insecurities now), which is exactly what one would expect.

From the list in Wolfram Mathworld, you might expect the result to be a constant line at y=1 instead of y=2. This difference of the factor 2 is just the consequence of considering the complete half circle or just one quarter segment of the circle (see e.g. in the sketch in https://en.wikipedia.org/wiki/Abel_transform )

So, your line at \$y approx 2\$ should be the right result (or at least the approximation to the analytical result y=2 constant). What do you mean with "a few other messy lines"?

Comment only
28 Nov 2013 R

### R (view profile)

Does it work with Octave? I tried a couple of analytic function/inversion pairs and it's not even close. For instance: f(x)=2*sqrt(a^2-x^2); Wolfram Mathworld says the inversion sould be PI_a(x), where PI_a(x)=0<x<a? 1:0. With a=6 and number of cos-expansions=10, on a reconstructed distribution graph, I get a line \$x \approx 2\$ for 0<=x<=6 and a few other messy lines.

Comment only
28 Nov 2013 Carsten Killer

### Carsten Killer (view profile)

@aneps: OK, I think now I see what your problem is about. The program will not work with scattered data, i.e. a list of x-y-pairs. You will have to bin your data (map it onto a 2D-grid), e.g. with bilinear interpolation.
After doing this, you will have a 2D-matrix which is some kind of a density plot of your measurement (like taking a picture with a camera).
If y is the axis of symmetry, each row of the matrix can be used as input H to obtain the radial distribution at that particular heigth (y-coordinate).

(note that if the axis of symmetry is located in the center of your matrix and not at x=0, you just have to use "half" of the rows. In an ideal measurement, with a perfect symmetry of your observable, the left and right half of the matrix would of course be identical. But in reality this will probably not be the case.)

The radius R is only important if one wants the result to be scaled to the units suiting your experimental setup. If you just want to check wether the algorith works for you, just set R=10 or any other arbitrary value.

The upper frequency defines how many cos-expansions will be used. (or, more precisely, up to which order the expansion is calculated). A low value for UPF will lead to a relatively rough approximation of your data, but it is on the other hend insensitive to (high-frequency) noise. Choosing higher values for UPF will result in a more and more exact calculation, which accounts for ever smaller nuances in your data.

Comment only
27 Nov 2013 aneps

### aneps (view profile)

@Carsten Killer: Thank you for the reply. I obtain the data by projecting the electron wave function on a plane detector (electron velocity map imaging). The axis of symmetry is Y axis (the polarization axis of the laser used). This is the 2D projection of a 3D wave function. The recorded data (points on which the electron collapses for each measurement) is a 2D distribution of points. My program records the coordinates (x,y) of each points. So, I have a big data of x,y values. SO, here how should I give the input arguments? I mean, how should I choose H, Radius R(which is unpredictable in this case), upper frequency etc?

Comment only
27 Nov 2013 Carsten Killer

### Carsten Killer (view profile)

@aneps: How is your data distributed along the x- and y-direction? If y is the axis of symmetry (e.g. the vertical coordinate in an upright standing cylindrical object), then I assume you want to perform the Abel inversion in x-direction.
Just take each row of your image and compute the Abel inversion for each row (i.e. for each y-value, you can of course do this in a loop) seperately, because the algorithm accepts only 1D-vectors as input.

Concerning the polar coordinates: You can put your measured data (which is in cartesian coordinates, I guess) directly into the algorithm.
However, if the axis of symmetry in your dataset is not along the x- or y-axis, you would have the perform some kind of transformation to produce a dataset that has an "easy" symmetry (i.e. in x- or y-direction).

Comment only
27 Nov 2013 aneps

### aneps (view profile)

I have a data points with X and Y coordinates. Can I use this algorithm to find the abel inversion of the data? Should I convert the data into polar coordinates? Could you please advise me how to use this algorithm to get the abel inversion for my data?

Comment only
23 Oct 2013 James

### James (view profile)

I tested this inversion against 3 known analytical abel pairs and found good agreement with this numerical solution. I found that with smooth data with good signal/noise, increasing upf provided a better answer.

23 Oct 2013 James

### James (view profile)

I was able to get the code working by replacing lsqcurvefit with a standard non-linear equations solver and simply supplying the residuals. Everything seems to be ok but I still need to through a few theoretical test curves at it. Thanks!

Comment only
22 Oct 2013 Carsten Killer

### Carsten Killer (view profile)

@James: Thank you for your suggestion, which is quite helpful. I changed that part of the algorithm - now the optimization toolbox is no longer necessary.

Comment only
21 Oct 2013 James

### James (view profile)

Would it be possible to use this code without having the optimization toolbox (lsqcurvefit function)? Does it simply find a scaled coefficient for each of the vectors?

Comment only
23 Oct 2013 1.1

Instead of using lsqcurvefit (which requires the Optimization Toolbox), the final equation system can now be solved algebraically

18 Feb 2014 1.3

Description has been updated. A paragraph regarding the interpretation of the result with respect to the real (physical) units in experimental situations has been added.

05 Jan 2016 1.4

I have been made aware of an error in the original article of G. Pretzler which results in missing factor 2. This factor 2 is now included in the code and I have deleted the paragraph in the description above which dealt with this issue.

07 Mar 2016 1.5

Improvement of the integration scheme (thanks to Benjamin Tadsen). Using the coordinate transformation r=sqrt(t^2+x^2), the full integration limits can be evaluated (in previous versions, an epsilon was substracted to prevent diverging integrals).