151
152
153
154
155
156
157 INTEGER I0, N0, N0IN, PP, TTYPE
158 REAL DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU
159
160
161 REAL Z( * )
162
163
164
165
166
167 REAL CNST1, CNST2, CNST3
168 parameter( cnst1 = 0.5630e0, cnst2 = 1.010e0,
169 $ cnst3 = 1.050e0 )
170 REAL QURTR, THIRD, HALF, ZERO, ONE, TWO, HUNDRD
171 parameter( qurtr = 0.250e0, third = 0.3330e0,
172 $ half = 0.50e0, zero = 0.0e0, one = 1.0e0,
173 $ two = 2.0e0, hundrd = 100.0e0 )
174
175
176 INTEGER I4, NN, NP
177 REAL A2, B1, B2, GAM, GAP1, GAP2, S
178
179
180 INTRINSIC max, min, sqrt
181
182
183
184
185
186
187 IF( dmin.LE.zero ) THEN
188 tau = -dmin
189 ttype = -1
190 RETURN
191 END IF
192
193 nn = 4*n0 + pp
194 IF( n0in.EQ.n0 ) THEN
195
196
197
198 IF( dmin.EQ.dn .OR. dmin.EQ.dn1 ) THEN
199
200 b1 = sqrt( z( nn-3 ) )*sqrt( z( nn-5 ) )
201 b2 = sqrt( z( nn-7 ) )*sqrt( z( nn-9 ) )
202 a2 = z( nn-7 ) + z( nn-5 )
203
204
205
206 IF( dmin.EQ.dn .AND. dmin1.EQ.dn1 ) THEN
207 gap2 = dmin2 - a2 - dmin2*qurtr
208 IF( gap2.GT.zero .AND. gap2.GT.b2 ) THEN
209 gap1 = a2 - dn - ( b2 / gap2 )*b2
210 ELSE
211 gap1 = a2 - dn - ( b1+b2 )
212 END IF
213 IF( gap1.GT.zero .AND. gap1.GT.b1 ) THEN
214 s = max( dn-( b1 / gap1 )*b1, half*dmin )
215 ttype = -2
216 ELSE
217 s = zero
218 IF( dn.GT.b1 )
219 $ s = dn - b1
220 IF( a2.GT.( b1+b2 ) )
221 $ s = min( s, a2-( b1+b2 ) )
222 s = max( s, third*dmin )
223 ttype = -3
224 END IF
225 ELSE
226
227
228
229 ttype = -4
230 s = qurtr*dmin
231 IF( dmin.EQ.dn ) THEN
232 gam = dn
233 a2 = zero
234 IF( z( nn-5 ) .GT. z( nn-7 ) )
235 $ RETURN
236 b2 = z( nn-5 ) / z( nn-7 )
237 np = nn - 9
238 ELSE
239 np = nn - 2*pp
240 gam = dn1
241 IF( z( np-4 ) .GT. z( np-2 ) )
242 $ RETURN
243 a2 = z( np-4 ) / z( np-2 )
244 IF( z( nn-9 ) .GT. z( nn-11 ) )
245 $ RETURN
246 b2 = z( nn-9 ) / z( nn-11 )
247 np = nn - 13
248 END IF
249
250
251
252 a2 = a2 + b2
253 DO 10 i4 = np, 4*i0 - 1 + pp, -4
254 IF( b2.EQ.zero )
255 $ GO TO 20
256 b1 = b2
257 IF( z( i4 ) .GT. z( i4-2 ) )
258 $ RETURN
259 b2 = b2*( z( i4 ) / z( i4-2 ) )
260 a2 = a2 + b2
261 IF( hundrd*max( b2, b1 ).LT.a2 .OR. cnst1.LT.a2 )
262 $ GO TO 20
263 10 CONTINUE
264 20 CONTINUE
265 a2 = cnst3*a2
266
267
268
269 IF( a2.LT.cnst1 )
270 $ s = gam*( one-sqrt( a2 ) ) / ( one+a2 )
271 END IF
272 ELSE IF( dmin.EQ.dn2 ) THEN
273
274
275
276 ttype = -5
277 s = qurtr*dmin
278
279
280
281 np = nn - 2*pp
282 b1 = z( np-2 )
283 b2 = z( np-6 )
284 gam = dn2
285 IF( z( np-8 ).GT.b2 .OR. z( np-4 ).GT.b1 )
286 $ RETURN
287 a2 = ( z( np-8 ) / b2 )*( one+z( np-4 ) / b1 )
288
289
290
291 IF( n0-i0.GT.2 ) THEN
292 b2 = z( nn-13 ) / z( nn-15 )
293 a2 = a2 + b2
294 DO 30 i4 = nn - 17, 4*i0 - 1 + pp, -4
295 IF( b2.EQ.zero )
296 $ GO TO 40
297 b1 = b2
298 IF( z( i4 ) .GT. z( i4-2 ) )
299 $ RETURN
300 b2 = b2*( z( i4 ) / z( i4-2 ) )
301 a2 = a2 + b2
302 IF( hundrd*max( b2, b1 ).LT.a2 .OR. cnst1.LT.a2 )
303 $ GO TO 40
304 30 CONTINUE
305 40 CONTINUE
306 a2 = cnst3*a2
307 END IF
308
309 IF( a2.LT.cnst1 )
310 $ s = gam*( one-sqrt( a2 ) ) / ( one+a2 )
311 ELSE
312
313
314
315 IF( ttype.EQ.-6 ) THEN
316 g = g + third*( one-g )
317 ELSE IF( ttype.EQ.-18 ) THEN
318 g = qurtr*third
319 ELSE
320 g = qurtr
321 END IF
322 s = g*dmin
323 ttype = -6
324 END IF
325
326 ELSE IF( n0in.EQ.( n0+1 ) ) THEN
327
328
329
330 IF( dmin1.EQ.dn1 .AND. dmin2.EQ.dn2 ) THEN
331
332
333
334 ttype = -7
335 s = third*dmin1
336 IF( z( nn-5 ).GT.z( nn-7 ) )
337 $ RETURN
338 b1 = z( nn-5 ) / z( nn-7 )
339 b2 = b1
340 IF( b2.EQ.zero )
341 $ GO TO 60
342 DO 50 i4 = 4*n0 - 9 + pp, 4*i0 - 1 + pp, -4
343 a2 = b1
344 IF( z( i4 ).GT.z( i4-2 ) )
345 $ RETURN
346 b1 = b1*( z( i4 ) / z( i4-2 ) )
347 b2 = b2 + b1
348 IF( hundrd*max( b1, a2 ).LT.b2 )
349 $ GO TO 60
350 50 CONTINUE
351 60 CONTINUE
352 b2 = sqrt( cnst3*b2 )
353 a2 = dmin1 / ( one+b2**2 )
354 gap2 = half*dmin2 - a2
355 IF( gap2.GT.zero .AND. gap2.GT.b2*a2 ) THEN
356 s = max( s, a2*( one-cnst2*a2*( b2 / gap2 )*b2 ) )
357 ELSE
358 s = max( s, a2*( one-cnst2*b2 ) )
359 ttype = -8
360 END IF
361 ELSE
362
363
364
365 s = qurtr*dmin1
366 IF( dmin1.EQ.dn1 )
367 $ s = half*dmin1
368 ttype = -9
369 END IF
370
371 ELSE IF( n0in.EQ.( n0+2 ) ) THEN
372
373
374
375
376
377 IF( dmin2.EQ.dn2 .AND. two*z( nn-5 ).LT.z( nn-7 ) ) THEN
378 ttype = -10
379 s = third*dmin2
380 IF( z( nn-5 ).GT.z( nn-7 ) )
381 $ RETURN
382 b1 = z( nn-5 ) / z( nn-7 )
383 b2 = b1
384 IF( b2.EQ.zero )
385 $ GO TO 80
386 DO 70 i4 = 4*n0 - 9 + pp, 4*i0 - 1 + pp, -4
387 IF( z( i4 ).GT.z( i4-2 ) )
388 $ RETURN
389 b1 = b1*( z( i4 ) / z( i4-2 ) )
390 b2 = b2 + b1
391 IF( hundrd*b1.LT.b2 )
392 $ GO TO 80
393 70 CONTINUE
394 80 CONTINUE
395 b2 = sqrt( cnst3*b2 )
396 a2 = dmin2 / ( one+b2**2 )
397 gap2 = z( nn-7 ) + z( nn-9 ) -
398 $ sqrt( z( nn-11 ) )*sqrt( z( nn-9 ) ) - a2
399 IF( gap2.GT.zero .AND. gap2.GT.b2*a2 ) THEN
400 s = max( s, a2*( one-cnst2*a2*( b2 / gap2 )*b2 ) )
401 ELSE
402 s = max( s, a2*( one-cnst2*b2 ) )
403 END IF
404 ELSE
405 s = qurtr*dmin2
406 ttype = -11
407 END IF
408 ELSE IF( n0in.GT.( n0+2 ) ) THEN
409
410
411
412 s = zero
413 ttype = -12
414 END IF
415
416 tau = s
417 RETURN
418
419
420