why does this trig iteration generates complex numbers?

im trying to write a program that analyses angles in a four bar mechanism given by three members d1&d2 (a kinked bar- kink at angle beta), d3 and d4. d1 and d4 are hinged at the origin and a point (h,v) respectively, where v is below the origin. The driving angle between the horiz and d1 is alpha, the angles between the horiz and d3, and the horiz and d4 are phi and theta respectively. I want to output theta, but i need to iterate to find theta and phi for a given set of d's. currently, my program is outputting complex values of theta and i cant figure out why... my suspicion is that the iteration doesnt work. THis is the code:

d1 = 2;
d2 = 0.5;
v = 2;
h = 5;
beta = pi()/4;
%set minimum length with flap and hinge to hinge lengths
h2h_len = sqrt(v*v + h*h)
flap_len = sqrt(d1*d1 + d2*d2)
min_len = h2h_len + flap_len
d3 = 1;
%loop for later calculations...
% while (d3 < min_len + 3)
  d4 = min_len + 2 - d3;
  alpha = pi()/4;
  %while (alpha < pi()/2)
      Ah = -d1*cos(alpha) + d2*cos(alpha - beta);
      Av = d1*sin(alpha) - d2*sin(alpha - beta);
      theta = pi()/4;
%ITERATION SECTION STARTS HERE
      count = 1;
      while (count <= 50)
          theta = asin((sqrt(d3*d3 - ((h-Ah)-d4*cos(theta))^2) + v + Av)/d4);
          count = count + 1;
      end
      theta
      phi = acos(((h - Ah) - d4*cos(theta))/d3)
      %alpha = alpha + 0.0873;
  %end
%     d3 = d3 + 1;
%end
end

[José-Luis]

For real elements of x outside the range [-1,1], asin(X) is complex.

Therefore your function might evaluation to one or more such values.

Answer 1

[Roger Stafford]

You are apparently trying to solve the trigonometric equation

   sin(theta) = (sqrt(d3*d3 - ((h-Ah)-d4*cos(theta))^2) + v + Av)/d4

but the iterative method you are using is definitely not the way to accomplish this. You have no reason to believe that your iteration would ever converge to a solution even if you started with a close approximation. Other techniques are certainly called for here. By the appropriate algebraic manipulation that equation can be put in this form:

   d4*sin(theta)-v-Av = sqrt(d3^2-(h-Ah-d4*cos(theta))^2)

By squaring both sides and further manipulation it becomes

   2*(v+Av)*d4*sin(theta)+2*((h-Ah)*d4*cos(theta) 
     = d4^2-d3^2+(v+Av)^2+(h-Ah)^2

This equation is of the form

   A*sin(theta)+B*cos(theta) = C

where A = 2*(v+Av)*d4, B = 2*((h-Ah)*d4, and C = d4^2-d3^2+(v+Av)^2+(h-Ah)^2. Its solution can be found using matlab's 'atan2' function. If we evaluate

   eta = atan2(A,B)

it will satisfy sin(eta) = A/sqrt(A^2+B^2) and cos(eta) = B/sqrt(A^2+B^2) so that

   cos(theta-eta) = cos(theta)*cos(eta)+sin(theta)*sin(eta)
     = (B*cos(theta)+A*sin(theta))/sqrt(A^2+B^2)
     = C/sqrt(A^2+B^2)

This finally allows a solution

   theta = eta + acos(C/sqrt(A^2+B^2)) = atan2(A,B) + acos(C/sqrt(A^2+B^2))

Of course you can also use iteration with the 'fzero' function, but that is iteration of a very different kind than you have attempted here.

[Roger Stafford]

I should have added that

   theta = atan2(A,B) - acos(C/sqrt(A^2+B^2))

is also a solution, as well as with the addition or subtraction of any multiple of 2*pi.

Also it should be noted that while all of these are certainly solutions to the equation

   2*(v+Av)*d4*sin(theta)+2*((h-Ah)*d4*cos(theta) 
     = d4^2-d3^2+(v+Av)^2+(h-Ah)^2

they may not all be solutions to

   sin(theta) = (sqrt(d3*d3-((h-Ah)-d4*cos(theta))^2)+v+Av)/d4

but some would instead satisfy

   sin(theta) = (-sqrt(d3*d3-((h-Ah)-d4*cos(theta))^2)+v+Av)/d4

Answer 2

[the cyclist]

You could set a condition breakpoint that checks for the first time theta goes complex, and trace back through the value that go into the calculation, and what went awry. Maybe you have a tiny bit of roundoff error that takes you out of the expected input range.

See this page for details on debugging, setting breakpoints, etc, if you don't know how to do that.