Re: Help with mathematical function

[Alan Linton]

On Fri, 2005-04-08 at 16:08 -0700, Robert A. Macy wrote:
> I am not a mathematician, so not sure I can ask this
> correctly, but...
> I need to calculate a "spiral" integral - calculate the
> area a vector, or "ray", sweeps out as a spiral collapses.
>  
> 
> The tip of the vector monotonically decreases down to near
> zero value.  with 500 of these spirals.  
> 
> The data array is composed of 500 by 100 complex values.
>  The data along each row (500 rows) spirals down to near
> zero values at the highest index value.  
> 
> I would like to calculate the area the ray sweeps out as
> the spiral collapses.  
> 
> I tried a simple approach like using the area defined by
> radius times the change in angle (delta angle)
> 
> average( abs(arrayvalues) ) * delta( angle(arrayvalues) )
> 
> and then integrated these values along each row.  The mesh
> plots looked promising, BUT...
> 
> The angle function has a limited range so their plus minus
> transitions played havoc with the results.  The spirals can
> be over 700 to 1000 degrees.  
> 
> Plus, there is noise in the data, so the spiral tends to
> "lump" around 0,0 not necessarily converge to it, which
> again causes wild angle transitions.  However, along a row
> between any two adjacent array values the data seems to be
> monotonically decreasing (spiralling down) to zero values.
>  
> 
> Does anyone know of a simple algorithm that will take each
> vectors' end points, calculate the area, and sum that up
> over the 100 values along each row?  
> 
>              - Robert -

If t is the unwrapped angle and r is the radius
let z=0.5*r.*r;
then area=sum(diff(t).*(z(1:end-1)+z(2:end)))/2
(using trapezoidal integration of z with respect to t)

Since there is noise in the data there is probably no point in using a
higher order integration method.

The following code demonstrates.
Hope this helps

Alan Linton

clear
grid('on')
t=(0:.1:8)*pi;
r=0.6*exp(0.15*t);
%r=0.6*t;

err=0.5;
x=r.*cos(t)+err*(2*rand(size(t))-1);
y=r.*sin(t)+err*(2*rand(size(t))-1);
xlabel('x')
ylabel('y')
title('raw xy data')
plot(x,y)
disp('press any key to continue')
fflush(stdout);
kbhit();

r=sqrt(x.*x+y.*y);
t=atan2(y,x);
t=t';
r=r';
xlabel('t')
ylabel('r')
title('raw tr data')
plot(t,r)
disp('press any key to continue')
fflush(stdout);
kbhit();

t=unwrap(t);
xlabel('x')
ylabel('y')
title('polynomial curve fit')
plot(x,y)
hold on
a=[ones(size(t)) t t.^2 t.^3 t.^4 t.^5 t.^6];
c=a\r;
rfit=a*c;
xfit=rfit.*cos(t);
yfit=rfit.*sin(t);
plot(xfit,yfit,'b')
hold off
disp('press any key to continue')
fflush(stdout);
kbhit();

xlabel('t')
ylabel('r')
title('unwrapped t')
plot(t,r)
hold on
plot(t,rfit,'b')
hold off
disp('press any key to continue')
fflush(stdout);
kbhit();

z=0.5*r.*r;
xlabel('t')
ylabel('0.5*r.*r')
title('curve to be integrated')
plot(t,z);
hold on
a=[t.^2 t.^3 t.^4 t.^5 t.^6];
c=a\z;
zfit=a*c;
plot(t,zfit,'b')
hold off
%trapezoidal integration
area=sum(diff(t).*(z(1:end-1)+z(2:end)))/2




-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------