femwind/findq.m

   1 function q=findq(a,n)
   2 % q=findq(a,n,tol)
   3 % find q such that 1+q+q^2 +...+q^(n-1)=a if a>1 else return q=1
   4 % in:
   5 %   a,n size 1
   6 % out:
   7 %   q size 1
   8 if a>n,
   9     q=1.01;
  10     for k=1:10000,
  11         qnew = (1+a*(q-1))^(1/n);
  12         diff = qnew - q;
  13         q = qnew;
  14         % fprintf('iteration %i q=%g diff=%g\n',k,q,diff)
  15         if abs(diff)<eps
  16             break
  17         end
  18     end
  19 else
  20     q=1;
  21 end