1 INTEGER FUNCTION ignpoi(mu)
84 REAL a0,a1,a2,a3,a4,a5,a6,a7,b1,b2,c,c0,c1,c2,c3,d,del,difmuk,e,
85 + fk,fx,fy,g,muold,muprev,omega,p,p0,px,py,q,s,t,u,v,x,xx
87 INTEGER j,k,kflag,l,ll,m
94 EXTERNAL ranf,sexpo,snorm
97 INTRINSIC abs,alog,
exp,float,ifix,max0,min0,sign,
sqrt
101 SAVE s, d, l, ll, omega, c3, c2, c1, c0, c, m, p, q, p0,
102 + a0, a1, a2, a3, a4, a5, a6, a7, fact, pp, muprev, muold
109 DATA muprev,muold/-1.0e37,-1.0e37/
110 DATA a0,a1,a2,a3,a4,a5,a6,a7/-.5,.3333333,-.2500068,.2000118,
111 + -.1661269,.1421878,-.1384794,.1250060/
112 DATA fact/1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880./
117 IF (mu.EQ.muprev) go to 10
118 IF (mu.LT.10.0) go to 120
137 10 g = mu + s*snorm()
138 IF (g.LT.0.0) go to 20
143 IF (ignpoi.GE.ll)
RETURN
150 IF (d*u.GE.difmuk*difmuk*difmuk)
RETURN
159 20
IF (mu.EQ.muold) go to 30
166 c1 = b1 - 6.*b2 + 45.*c3
167 c0 = 1. - b1 + 3.*b2 - 15.*c3
169 30
IF (g.LT.0.0) go to 50
178 40
IF (fy-u*fy.LE.py*
exp(px-fx))
RETURN
188 IF (t.LE. (-.6744)) go to 50
189 ignpoi = ifix(mu+s*t)
200 60
IF (c*
abs(u).GT.py*
exp(px+e)-fy*
exp(fx+e)) go to 50
206 70
IF (ignpoi.GE.10) go to 80
208 py = mu**ignpoi/fact(ignpoi+1)
215 80 del = .8333333e-1/fk
216 del = del - 4.8*del*del*del
218 IF (
abs(v).LE.0.25) go to 90
219 px = fk*alog(1.0+v) - difmuk - del
222 90 px = fk*v*v* (((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0) -
224 100 py = .3989423/
sqrt(fk)
225 110 x = (0.5-difmuk)/s
228 fy = omega* (((c3*xx+c2)*xx+c1)*xx+c0)
229 IF (kflag.LE.0) go to 40
236 IF (mu.EQ.muold) go to 130
238 IF (mu.GE.0.0) go to 125
239 WRITE (*,*)
'MU < 0 in IGNPOI - ABORT'
240 WRITE (*,*)
'Value of MU: ',mu
241 CALL xstopx(
'MU < 0 in IGNPOI - ABORT')
260 IF (l.EQ.0) go to 150
262 IF (u.GT.0.458) j = min0(l,m)
264 IF (u.LE.pp(k)) go to 180
266 IF (l.EQ.35) go to 130
276 IF (u.LE.q) go to 170
octave_value sqrt(void) const