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