Re: Help with mathematical function

[Robert A. Macy]

Stefan,

Thank you for your reply.  

The problem is the data only "looks" like a spiral.  There
is actual information in the way the spiral decays,
including "lumpiness".  

My problem was that the angle() function resets itself as
the spiral decays.  It is very easy to get erroneously huge
delta angle changes at those junctions.  

Luckily the angle steps are known to be small enough that
one should NEVER get too large a change in angle.  

So my solution was to add some simple "if" tests.
v(i,k) is an array of complex values...

for i=1:500
  for k=2:100
...processing v(i,k)...
delta_angle(i,k)=angle(v(i,k))-angle(v(i,k-1));
if { delta_angle(i,k) > pi() )
  delta_angle(i,k) = delta_angle(i,k) - 2*pi();
endif
if ( delta_angle(i,k) < -pi() )
  delta_angle(i,k) = delta_angle(i,k) + 2*pi();
endif
...more processing of results(i,k)...
  endfor
endfor


The tests still fail when the spiral is really close to
zero, but by then the contribution to error is quite small
since the absolute vector length is so small.  

The only problem now is that the loops to do this, the
for..for.. loops, take several minutes.  

Is there someway to use vector math?  And still include the
checks on the angles?  Or, at least disable the constant
testing of the sizes of the variables which goes on each
time the loop is exercised?

        - Robert -

On Sat, 9 Apr 2005 18:34:49 +0200
 Stefan van der Walt <address@hidden> wrote:
> Hi Robert
> 
> I assumed that the radius of the spiral decreases
> linearly towards
> zero, as one moves through 360 degrees.
> 
> The following script draws the spiral and calculates the
> area for
> angles between 0 and 360 (0 and 2*pi radian).
> 
> Will this do, or do you need to take the measured data
> into account?
> In that case, your method is nearly correct.  You should,
> however,
> approximate the area as that of a rectangular triangle:
> 
> 1/2 * radius^2 * delta-angle
> 
> (radius * delta-angle is the length of one side of the
> trangle, radius
> the length of the other)
> 
> This is the same method used to solve the line integral
> of the
> theoretical function,
> 
> integral 1/2 (2 pi - theta)^2 R^2 dtheta
> 
> Regards
> Stefan
> 
> ----
> R = 1;
> theta = linspace(0, 2*pi, 1000);
> 
> area = R^2 * theta / 6 .* (theta.^2 - 6*pi*theta +
> 12*pi^2);
> 
> automatic_replot = 0;
> subplot(121);
> axis("square");
> grid on;
> polar(theta, (2*pi - theta)*R, ";Spiral;");
> 
> subplot(122);
> plot(theta, area, ";Area under spiral;")
> ----
> 
> On Fri, Apr 08, 2005 at 04:08:15PM -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 -
> > 



-------------------------------------------------------------
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
-------------------------------------------------------------