1 /*
2  * $Id: series.i,v 1.1 2005-09-18 22:06:07 dhmunro Exp $
3  * Routines for handling geometric series (e.g.- tapered zoning).
4  */
5 /* Copyright (c) 2005, The Regents of the University of California.
6  * All rights reserved.
7  * This file is part of yorick (http://yorick.sourceforge.net).
8  * Read the accompanying LICENSE file for details.
9  */
10 
series_s(r,n)11 func series_s(r, n)
12 /* DOCUMENT series_s(r, n)
13      returns the sum s of the finite geometric series
14         1 + r + r^2 + r^3 + ... + r^n
15      Using n<0 is equivalent to using the reciprocal of r, that is,
16         series_s(r, -n) == series_s(1./r, n)
17      R or N or both may be arrays, as long as they are conformable.
18    SEE ALSO: series_r, series_n
19  */
20 {
21   /* force conformability immediately */
22   dims= dimsof(r, n);
23   rr= r + array(0.0, dims);
24   nn= n + array(0, dims);
25 
26   /* form result array, initialized to n==0 result */
27   s= array(1.0, dims);
28 
29   /* subdivide into n==0, n>0, n<0, and r==1.0 cases */
30   rnot= rr!=1.0;
31 
32   mask= (nn>0 & rnot);   /* n>0 case */
33   list= where(mask);
34   if (numberof(list)) {
35     rrx= rr(list);
36     s= merge((rrx^(1+nn(list))-1.0)/(rrx-1.0), s(where(!mask)), mask);
37   }
38 
39   mask= (nn<0 & rnot);   /* n<0 case */
40   list= where(mask);
41   if (numberof(list)) {
42     rrx= 1.0/rr(list);
43     s= merge((rrx^(1-nn(list))-1.0)/(rrx-1.0), s(where(!mask)), mask);
44   }
45 
46   list= where(!rnot);   /* r==1.0 case */
47   if (numberof(list)) s= merge(s(where(rnot)), abs(nn(list)), rnot);
48 
49   return s;
50 }
51 
series_r(s,n)52 func series_r(s, n)
53 /* DOCUMENT series_r(s, n)
54      returns the ratio r of the finite geometric series, given the sum s:
55         1 + r + r^2 + r^3 + ... + r^n = s
56      Using n<0 will return the the reciprocal of n>0 result, that is,
57         series_r(s, -n) == 1.0/series_r(s, n)
58      If n==0, returns s-1 (the n==1 result).
59      S or N or both may be arrays, as long as they are conformable.
60    SEE ALSO: series_s, series_n
61  */
62 {
63   /* force conformability immediately */
64   dims= dimsof(s, n);
65   ss= s + array(0.0, dims);
66   nn= n + array(0, dims);
67 
68   /* form result array, initialized to abs(n)<2 result */
69   r= ss-1.0;
70 
71   /* work only with n>0 -- take reciprocals at end if necessary */
72   nneg= nn<0;
73   nn= abs(nn)+1;
74   nbig= nn>2;
75 
76   /* compute an approximate result which has exact values and
77      derivatives at s==1, s==n, and s->infinity --
78      different approximations apply for s>n and s<n */
79   mask= nbig & (ss>nn);
80   list= where(mask);
81   if (numberof(list)) {
82     sx= ss(list);
83     nx= nn(list);
84     pow= 1.0/(nx-1.0);
85     npow= nx^pow - 1.0;
86     n2r= 1.0/(nx-2.0);
87     A= (2.0-nx*npow)*n2r;
88     B= (2.0*npow-nx*pow)*nx*n2r;
89     r= merge(sx^pow - pow + A*(nx/sx)^pow + B/sx, r(where(!mask)), mask);
90   }
91   mask= nbig & (ss<=nn);
92   list= where(mask);
93   if (numberof(list)) {
94     sx= ss(list);
95     nx= nn(list);
96     sn= (sx-1.0)/(nx-1.0);
97     n2r= 1.0/(nx*nx);
98     r= merge(1.0 - 1.0/sx + n2r*sn*sn*(nx+1.0 - sn), r(where(!mask)), mask);
99   }
100 
101   /* Polish the approximation using Newton-Raphson iterations.
102      There are never many of these, so do the entire vector together.  */
103   mask= nbig & (ss!=nn);
104   list= where(mask);
105   if (numberof(list)) {
106     rx= r(list);
107     ss= ss(list);
108     nn= nn(list);
109     for (;;) {
110       rr= rx-1.0;
111       rn= rx^(nn-1);
112       rrss= rr*ss;
113       delta= rrss - (rx*rn-1.0);
114       if (allof(abs(delta)<=1.e-9*abs(rrss))) break;
115       rx+= delta/(nn*rn-ss);
116     }
117     /* try to get it to machine precision */
118     if (anyof(delta)) rx+= delta/(nn*rn-ss);
119     r= merge(rx, r(where(!mask)), mask);
120   }
121 
122   list= where(nneg);
123   if (numberof(list)) r= merge(1.0/r(list), r(where(!nneg)), nneg);
124 
125   return r;
126 }
127 
series_n(r,s)128 func series_n(r, s)
129 /* DOCUMENT series_n(r, s)
130      returns the minimum number n of terms required for the geometric
131      series
132         1 + r + r^2 + r^3 + ... + r^n = s
133      to reach at least the given value s.  An alternate viewpoint is
134      that n is the minimum number of terms required to achieve the
135      sum s, with a ratio no larger than r.
136      Returns 0 if r<1 and s>1/(1-r), or if s<1.
137      The routine makes the most sense for r>1 and s substantially
138      greater than 1.  The intended use is to determine the minimum
139      number of zones required to span a given thickness t with a given
140      minimum zone size z, and maximum taper ratio r (assumed >1 here):
141         n= series_n(r, t/z);
142      With this n, you have the option of adjusting r or z downwards
143      (using series_r or series_s, respectively) to achieve the final
144      desired zoning.
145      R or S or both may be arrays, as long as they are conformable.
146    SEE ALSO: series_s, series_r
147  */
148 {
149   n= 1.0 + (r-1.0)*s;
150   bad= (n<=0.0 | s<1.0);
151   list= where(bad);
152   if (numberof(list))
153     n= merge(array(1.0, numberof(list)), n(where(!bad)), bad);
154   mask= r==1.0 & !bad;
155   list= where(mask);
156   if (numberof(list))
157     n= merge((s+0.0*n)(list), n(where(!mask)), mask);
158   mask= !mask & !bad;
159   list= where(mask);
160   if (numberof(list))
161     n= merge(log(n(list))/log(r(list)), n(where(!mask)), mask);
162   return long(ceil(n))-1;
163 }
164