source: trunk/src/seismic_processing/rayloc_ew/robust_util.f @ 3168

Revision 3168, 54.9 KB checked in by withers, 11 years ago (diff)

Synced with hydra version of rayloc

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1      subroutine adder(navsum,ihold)
2c
3c $$$$$ calls drdwr $$$$$
4c
5c   Adder constructs minus the gradient of the dispersion function and
6c   also the mean of the travel-time derivitives. Navsum is the number
7c   of data to be used and ihold is the number of constraints.  Ihold
8c   will be 0 for a hypocentral step and 1 for an epicentral step.
9c   Programmed on 10 December 1983 by R. Buland.
10c
11      include 'locprm.inc'
12      double precision a,b,c,x,s,avdr
13c   Physical constants.
14      common /circle/ pi,rad,deg,radius
15c   Storage for the normal equations.
16      common /matrix/ a(3,4),b(3),c(4),x(3),s(3),avdr(3)
17c   Hypocentral parameters.
18      include 'hypocm.inc'
19c   Station-reading parameters.
20      include 'pickcm.inc'
21c   Scores.
22      common /sccom/ scgen(maxph)
23      data eps/1e-5/
24c
25      m1=3-ihold
26c   Zero out array storage.
27      do 1 k=1,m1
28      b(k)=0d0
29 1    avdr(k)=0d0
30c   Loop over the stations accumulating the sums.
31      do 2 i=1,navsum
32c   Pop one of the sorted residuals.
33      ierr=lpop(i,j,res)
34c   Access station parameter storage.
35      call drdwr(1,j,1,5)
36c   C is the contribution to the normal equations.
37      c(1)=deg*cos(azim)*dtdd
38      c(2)=-deg*sin(azim)*dtdd
39      c(3)=dtdh
40c   Update the sums.
41      do 3 k=1,m1
42 3    b(k)=b(k)+scgen(i)*c(k)
43      c(2)=c(2)*olats
44c     print *,'c',j,(c(k),k=1,3)
45      do 6 k=1,m1
46 6    avdr(k)=avdr(k)+c(k)
47 2    continue
48      b(1)=b(1)/radius
49      b(2)=b(2)/radius
50c   Compute the mean.
51      fac=0.
52      do 4 k=1,m1
53 4    fac=fac+b(k)*b(k)
54      cn=1./navsum
55      fac=1./sqrt(amax1(fac,eps))
56      do 5 k=1,m1
57      b(k)=fac*b(k)
58 5    avdr(k)=cn*avdr(k)
59      return
60      end
61
62      subroutine cholin(ndim,n,a,jstage,ifbad)
63c
64c $$$$$ calls only library routines $$$$$
65c
66c   Routine for the inversion and factorization of a positive definite,
67c   symmetric matrix based on an Algol routine given by Wilkinson and
68c   Reinsch, Linear Algebra, p16-17.  Cholin may be used to recover the
69c   factor L in the Cholesky factorization A = L*L-trans, or the inverse
70c   of L, or the inverse of A found by Cholesky.  Storage has been
71c   transposed from the Algol original for the convenience of Fortran
72c   users.  As a result, it is convienent to supply L-trans or inverse
73c   L-trans when L or inverse L respectively have been called for.  A
74c   must be dimension at least A(ndim,n+1).  On entry, the input matrix
75c   must occupy the first n columns, or if desired, just the lower left
76c   triangle.  On exit, the upper triangle of the array comprised of
77c   columns 2,3...n+1 of A contains the desired result.  Jstage
78c   indicates the desired calculation:  1 to determine L-trans, 2 to
79c   determine inverse L-trans, 0 or 3 to determine inverse A.  If jstage
80c   is -2 or -3, then inverse L-trans or inverse A respectively is
81c   determined assuming that L-trans has already been found in a
82c   previous call to cholsl.  Note that the lower triangle of the
83c   orignal A matrix is not destroyed.  Ifbad is an error flag which
84c   will be set to 0 if the calculation has been sucessful or to 1 if
85c   the matrix is not positive definite.  Translated from Algol by
86c   R. L. Parker.
87c
88      implicit double precision (a-h,o-z)
89      dimension a(ndim,1)
90c
91c   Stage 1:  form factors (L-trans is stored above the diagonal).
92c
93      ifbad=0
94      istage=iabs(jstage)
95      if (jstage.lt.0) go to 1550
96      do 1500 i=1,n
97      i1=i + 1
98      do 1400 j=i,n
99      j1=j + 1
100      x=a(j,i)
101      if (i.eq.1) go to 1200
102      do 1100 kback=2,i
103      k=i-kback+1
104 1100 x=x - a(k,j1)*a(k,i1)
105 1200 if (i.ne.j) go to 1300
106      if (x.le.0d0) go to 4000
107      y=1d0/dsqrt(x)
108      a(i,i1)=y
109      go to 1400
110 1300 a(i,j1)=x*y
111 1400 continue
112 1500 continue
113      if (istage.ne.1) go to 2000
114 1550 do 1600 i=1,n
115 1600 a(i,i+1)=1d0/a(i,i+1)
116      if (istage.eq.1) return
117c
118c   Stage 2:  construct inverse L-trans.
119c
120 2000 do 2600 i=1,n
121      i1=i+1
122      if (i1.gt.n) go to 2600
123      do 2500 j=i1,n
124      z=0d0
125      j1=j + 1
126      j2=j - 1
127      do 2200 kback=i,j2
128      k=j2-kback+i
129 2200 z=z - a(k,j1)*a(i,k+1)
130      a(i,j1)=z*a(j,j1)
131 2500 continue
132 2600 continue
133      if (istage.eq.2) return
134c
135c   Stage 3:  construct the inverse of A above the diagonal.
136c
137      do 3500 i=1,n
138      do 3500 j=i,n
139      z=0d0
140      n1=n + 1
141      j1=j + 1
142      do 3200 k=j1,n1
143 3200 z=z + a(j,k)*a(i,k)
144      a(i,j+1)=z
145 3500 continue
146      return
147c   Error exit.
148 4000 ifbad=1
149      return
150      end
151
152      subroutine cholsl(ndim,n,a,x,b,istage)
153c
154c $$$$$ calls no other routine $$$$$
155c
156c   Routine for solving the n x n linear system A*x = b by back-
157c   substitution assuming that A(ndim,n+1) has already been factored by
158c   calling cholin with jstage = 1.  See cholin for details on the
159c   organization of A.  The right-hand-side vector, b(n), must be given
160c   and the unknown vector, x(n), will be determined.  Istage is the
161c   mode of the solution:  1 to not complete the back-substitution,
162c   returning x = inverse(L-trans)*b; otherwiser complete the solution
163c   as described above.  Translated from Algol by R. L. Parker.
164c
165      implicit double precision (a-h,o-z)
166      dimension a(ndim,1),x(n),b(n)
167c
168      do 2500 i=1,n
169      z=b(i)
170      if (i.eq.1) go to 2500
171      do 2200 kk=2,i
172      k=i-kk+1
173 2200 z=z - a(k,i+1)*x(k)
174 2500 x(i)=z/a(i,i+1)
175      if (istage.eq.1) return
176c   Complete the back-substitution.
177      do 2700 ii=1,n
178      i=n-ii+1
179      z=x(i)
180      i1=i+1
181      if (i.eq.n) go to 2700
182      do 2600 k=i1,n
183 2600 z= z - a(i,k+1)*x(k)
184 2700 x(i)=z/a(i,i+1)
185      return
186      end
187
188      subroutine drdwr(irw,irec,i1,i2)
189c
190c $$$$$ calls no other routine $$$$$
191c
192c   Drdwr manages station parameter storage for reloc.  Up to max
193c   stations will be stored in main memory.  Additional stations will
194c   be stored on disc.  The storage media is transparant to reloc and
195c   the parameter max need be set only in drdwr.  If irw = 1, a read
196c   into commons /stadat/ and /stadtc/ will be done from station record
197c   number irec.  Only elements i1,i1+1,...,i2 of common /stadat/ will
198c   be copied.  If irw = 2, a write from the commons will be done to
199c   record irec.  The partial transfers are simply to save time on
200c   transfers to and from main memory.  Transfers to and from disc are
201c   necessarily always complete.  Note that stations should have record
202c   numbers:  irec = 1,2,3,...,nsta, where nsta is the number of
203c   stations.  Records 1 - max will be in memory and records (max+1) -
204c   nsta will be disc records 1 - (nsta-max).  Rewritten by R. Buland on
205c   19 July 1983 from a similar routine by E. R. Engdahl.
206c
207      include 'locprm.inc'
208        character*22 flgs,cstor
209c   Logical unit numbers.
210      common /unit/ min,mout,merr
211c   Communication with outside world (station parameter data).
212      common /stadat/ sdat(13)
213      common /stadtc/ flgs
214c   Main memory storage.
215      common /bnstor/ rstor(22,maxph)
216      common /chstor/ cstor(maxph)
217c       integer stdout, stdin
218c
219c       stdout = 6
220c       stdin  = 5
221c
222c   Branch on operation type.
223      if (irw.ne.1) go to 1
224c
225c   Do a read operation.
226c
227      if(irec.gt.maxph) go to 2
228c   Transfer from memory.
229      do 3 i=i1,i2
230      sdat(i)=rstor(i,irec)
231 3      continue
232      flgs=cstor(irec)
233      return
234c
235c   Do a write operation.
236c
237 1    if(irec.gt.maxph) go to 2
238c   Transfer to memory.
239      do 5 i=i1,i2
240      rstor(i,irec)=sdat(i)
2415       continue
242      cstor(irec)=flgs
243      return
244c   Write to error report.
2452       write(merr,*)'Maximun stations exceeded!! '
246        write(merr,*)'irec= ',irec,' > maxph = ',maxph
247        call exit(101)
248      end
249
250      function dsprsn(n,a)
251c
252c $$$$$ calls scgen $$$$$
253c
254c   Given n residuals, a(1) <= a(2) <= ... <= a(n), dsprsn returns the
255c   dispersion of the residuals as defined by Jaeckel,L.A.,1972,
256c   "Estimating regression coefficients by minimizing the dispersion of
257c   the residuals":AMS,v43,p1449-58.  Scgen is the score generating
258c   function which must be nondecreasing and the sum of the scores must
259c   be zero.  Programmed on 10 December 1983 by R. Buland.
260c
261      include 'locprm.inc'
262      dimension a(1)
263      common/sccom/scgen(maxph)
264c
265      d=0.
266      do 1 i=1,n
267 1    d=d+scgen(i)*a(i)
268      dsprsn=d
269      return
270      end
271
272      subroutine errelp(n,m,se)
273c
274c $$$$$ calls f10 and jacobi $$$$$
275c
276c   Errelp computes the azimuth, plunge, and semi-length of the
277c   principal axes of the three dimensional (two if depth is held)
278c   spatial error ellipsoid at a 90% confidence level.  Asymptotically, 
279c   the rank-sum estimator produces parameters (i.e., the hypocenter) 
280c   that are chi-squared distributed.  An empirical model of 
281c   travel-time residuals has been used to compute the constants.  N is 
282c   the number of data used, m is the number of degrees of freedom 
283c   (including origin time), and se is the standard deviation of the 
284c   residuals.  Programmed on 28 October 1983 by R. Buland.
285c
286      include 'locprm.inc'
287      double precision a,b,c,x,s,ev(3,3),emax,ce,tol,tol2,avdr
288c   Physical constants.
289      common /circle/ pi,rad,deg,radius
290c   Storage for the normal equations.
291      common /matrix/ a(3,4),b(3),c(4),x(3),s(3),avdr(3)
292c   Principal axies.
293      include 'statcm.inc'
294      data tol,tol2/1d-15,2d-15/
295c
296c   Initialize the principal axies.
297      do 1 j=1,3
298      do 1 i=1,3
299 1    prax(i,j)=0.
300      avh=0.
301c   If all parameters are held, we are done.
302      if(m.lt.3) return
303c   M1 is the number of degrees of freedom of the spatial projection of
304c   the correlation matrix.
305      m1=m-1
306c   Decompose the correlation matrix.
307      call jacobi(m1,3,3,a,b,ev,c,x)
308      cn=rsc3
309      if(m1.lt.3) cn=rsc2
310c
311c   Loop over the eigenvalues.
312c
313      do 2 i=1,m1
314c   Sort the eigenvalues into order of decreasing size.
315      emax=-1d0
316      do 3 j=1,m1
317      if(b(j).le.emax) go to 3
318      k=j
319      emax=b(j)
320 3    continue
321c   Find the azimuth, plunge, and standard deviation for each axis.
322      ce=1d0
323      if(m1.gt.2) ce=dsign(1d0,ev(3,k))
324      if(dabs(ev(1,k))+dabs(ev(2,k)).gt.tol2) prax(i,1)=deg*datan2(ce*
325     1 ev(2,k),-ce*ev(1,k))
326      if(prax(i,1).lt.0.) prax(i,1)=prax(i,1)+360.
327      if((dabs(ev(3,k)).le.tol).and.(prax(i,1).gt.180.)) prax(i,1)=
328     1 prax(i,1)-180.
329      if(m1.gt.2) prax(i,2)=deg*dasin(dmin1(ce*ev(3,k),1d0))
330      prax(i,3)=cn*dsqrt(dmax1(b(k),0d0))
331c   Guarantee that this axis is not found again.
332 2    b(k)=-1d0
333c
334c   Compute avh.
335      if(m.le.3) then
336c   If depth was held, we already have all the information we need.
337        avh=rsc1*sqrt(prax(1,3)*prax(2,3))/cn
338      else
339c   If depth was free, we need to re-decompose the latitude-longitude
340c   correlation sub-matrix.
341        call jacobi(2,3,3,a,b,ev,c,x)
342        avh=rsc1*dsqrt(dsqrt(dmax1(b(1)*b(2),0d0)))
343      endif
344      return
345      end
346
347      function estlin(navsum,ihold,xx)
348c
349c $$$$$ calls drdwr $$$$$
350c
351c   Estlin returns a linear estimate of the dispersion function at a
352c   position z(k)=y(k)+xx*x(k), where y(k) is the current position.
353c   Programmed on 12 December 1983 by R. Buland.
354c
355      double precision a,b,c,x,s,avdr
356      dimension dr(3)
357c
358c   Logical units numbers
359c
360      common /unit/ min,mout,merr
361c
362c   Storage for the normal equations.
363c
364      common /matrix/ a(3,4),b(3),c(4),x(3),s(3),avdr(3)
365c   Hypocentral parameters.
366      include 'hypocm.inc'
367c   Station-reading parameters.
368      include 'pickcm.inc'
369c
370      m1=3-ihold
371      n=0
372c   Zero out array storage.
373      do 1 k=1,m1
374 1    dr(k)=xx*x(k)
375c     print *,'dr ',dr
376c   Loop over the stations accumulating the sums.
377      do 2 i=1,navsum
378c   Pop one of the sorted residuals.
379      ierr=lpop(i,j,res)
380c   Access station parameter storage.
381      call drdwr(1,j,1,5)
382c     print *,'f2 ',j,-cos(azim)*dtdd,sin(azim)*olats*dtdd,-dtdh
383      r0=res-(dr(1)*cos(azim)-dr(2)*sin(azim)*olats)*dtdd-dr(3)*dtdh
384c     print *,'res ',j,res,r0
385      res=r0
386 2    ierr=lpush(n,j,res)
387      ierr=lsrt(navsum,av,sd,chisq,'TT')
388      estlin=chisq
389      return
390      end
391
392      subroutine jacobi(n,l1,l2,a,d,v,b,z)
393c
394c $$$$$ calls only library routines $$$$$
395c
396c   Jacobi finds all eigenvalues and eigenvectors of the n by n real
397c   symmetric matrix a by Jacobi's method.  The unordered eigenvalues
398c   are returned in real n vector d and the corresponding eigenvectors
399c   are in the respective columns of n by n real matrix v.  B and z are
400c   scratch arrays.  The arrays must be dimensioned at least a(l,n),
401c   d(n), v(n,n), b(n), and z(n).  Note that the inner dimension of
402c   array a is l.  Taken from Linear Algebra, Volume II, by Wilkinson
403c   and Reinsch, 1971, Springer-Verlag.  The algorithm is published
404c   (in algol) by H. Rutishauser on pages 202-211.  Converted to Fortran
405c   by R. Buland, 30 October 1978.
406c
407      implicit double precision (a-h,o-z)
408      integer p,q,p1,p2,q1,q2
409      double precision a(l1,n),d(n),v(l2,n),b(n),z(n)
410      n1=n-1
411      nn=n*n
412c   Initialize array storage.
413      do 1 p=1,n
414      do 2 q=1,n
415 2    v(p,q)=0d0
416      v(p,p)=1d0
417      b(p)=a(p,p)
418      d(p)=b(p)
419 1    z(p)=0d0
420c
421c   Make up to 50 passes rotating each off diagonal element.
422c
423      do 3 i=1,50
424      sm=0d0
425      do 4 p=1,n1
426      p1=p+1
427      do 4 q=p1,n
428 4    sm=sm+dabs(a(p,q))
429c   Exit if all off diagonal elements have underflowed.
430      if(sm.eq.0d0) go to 13
431      tresh=0d0
432      if(i.lt.4) tresh=.2d0*sm/nn
433c
434c   Loop over each off diagonal element.
435c
436      do 5 p=1,n1
437      p1=p+1
438      p2=p-1
439      do 5 q=p1,n
440      q1=q+1
441      q2=q-1
442      g=1d2*dabs(a(p,q))
443c   Skip this element if it has already underflowed.
444      if((i.le.4).or.(dabs(d(p))+g.ne.dabs(d(p))).or.
445     1 (dabs(d(q))+g.ne.dabs(d(q)))) go to 6
446      a(p,q)=0d0
447      go to 5
448 6    if(dabs(a(p,q)).le.tresh) go to 5
449c   Compute the rotation.
450      h=d(q)-d(p)
451      if(dabs(h)+g.eq.dabs(h)) go to 7
452      theta=.5d0*h/a(p,q)
453      t=1d0/(dabs(theta)+dsqrt(1d0+theta*theta))
454      if(theta.lt.0d0) t=-t
455      go to 14
456 7    t=a(p,q)/h
457 14   c=1d0/dsqrt(1d0+t*t)
458      s=t*c
459      tau=s/(1d0+c)
460c   Rotate the diagonal.
461      h=t*a(p,q)
462      z(p)=z(p)-h
463      z(q)=z(q)+h
464      d(p)=d(p)-h
465      d(q)=d(q)+h
466      a(p,q)=0d0
467c   Rotate the off diagonal.  Note that only the upper diagonal elements
468c   are touched.  This allows the recovery of matrix a later.
469      if(p2.lt.1) go to 15
470      do 8 j=1,p2
471      g=a(j,p)
472      h=a(j,q)
473      a(j,p)=g-s*(h+g*tau)
474 8    a(j,q)=h+s*(g-h*tau)
475 15   if(q2.lt.p1) go to 16
476      do 9 j=p1,q2
477      g=a(p,j)
478      h=a(j,q)
479      a(p,j)=g-s*(h+g*tau)
480 9    a(j,q)=h+s*(g-h*tau)
481 16   if(n.lt.q1) go to 17
482      do 10 j=q1,n
483      g=a(p,j)
484      h=a(q,j)
485      a(p,j)=g-s*(h+g*tau)
486 10   a(q,j)=h+s*(g-h*tau)
487c   Rotate the eigenvector matrix.
488 17   do 11 j=1,n
489      g=v(j,p)
490      h=v(j,q)
491      v(j,p)=g-s*(h+g*tau)
492 11   v(j,q)=h+s*(g-h*tau)
493 5    continue
494c   Reset the temporary storage for the next rotation pass.
495      do 3 p=1,n
496      d(p)=b(p)+z(p)
497      b(p)=d(p)
498 3    z(p)=0d0
499c
500c   All finished.  Prepare for exiting by reconstructing the upper
501c   triangle of a from the untouched lower triangle.
502c
503 13   do 12 p=1,n1
504      p1=p+1
505      do 12 q=p1,n
506 12   a(p,q)=a(q,p)
507      return
508      end
509
510      function lsrt(l,av,sd,chisq,tag)
511c
512c   Function Lsrt sorts the l data in array del2 (reordering ind2
513c   accordingly) and computes the series median, av, spread, sd, and
514c   dispersion, chisq.  Tag is a two character ID printed with errors
515c   to help identify the source of the problem.
516c
517c   At most, min2 data may be considered.
518      include 'locprm.inc'
519      parameter (min2=maxph)
520      character*2 tag
521      character*1 pflg
522      character*20 alpha
523        common /unit/ min, mout, merr
524c   Elimination scratch storage.
525      common /elimsc/ ind2(min2),del2(min2)
526c   Station flags.
527      common /stadtc/ alpha,pflg
528c
529c   Assume success.
530      lsrt=1
531c   Check the number of data.
532      if(l-1)16,17,18
533c   No data.
534 16   av=0.
535      lsrt=-1
536      return
537c   One data.
538 17   av=del2(1)
539      return
540c   More than one datum.  Issue an error message if necessary.
541 18   if(mfl.gt.0) write(merr,100)mfl,tag
542 100  format(/1x,'Reloc:  scratch storage full -',i4,1x,a,'''s ',
543     1 'residualed.'/)
544c   Sort all data in the scratch storage.
545      call magsrt(l,del2,ind2)
546c   Compute the median, spread, and dispersion.
547      av=xmed(l,del2,sd)
548      chisq=dsprsn(l,del2)
549      return
550c
551c   Entry lpush increments storage pointer n and saves index, ndex, and
552c   datum, dat, in arrays ind2 and del2.  If successful, +1 is returned. 
553c   If no more storage is available, -1 is returned, pflg is set to
554c   residual, and a the error count is bumped.
555c
556      entry lpush(n,ndex,dat)
557c   Assume success.
558      lpush=1
559c   On the first datum, initialize the error count.
560      if(n.le.0) mfl=0
561c   See if elimination scratch storage is full.
562      if(n.ge.min2) go to 13
563c   Push the datum onto the stack.
564      n=n+1
565      ind2(n)=ndex
566      del2(n)=dat
567      return
568c   Scratch storage is full.  Set the error flag.
569 13   lpush=-1
570c   Change the bad datum flag from 'u' to ' '.
571      pflg=' '
572c   Up the error count.
573      mfl=mfl+1
574      return
575c
576c   Entry lpop returns the n th saved values of index, ndex, and datum,
577c   dat, from arrays ind2 and del2.  If successful, +1 is returned.  If
578c   n is beyond array storage, -1 is returned.
579c
580      entry lpop(n,ndex,dat)
581c   Assume success.
582      lpop=1
583c   See if n is beyond array storage.
584      if(n.ge.min2) go to 15
585c   Pop the datum from the stack.
586      ndex=ind2(n)
587      dat=del2(n)
588      return
589c   N is too big.  Set the error flag.
590 15   lpop=-1
591      return
592      end
593
594      subroutine lintry(navsum,ihold,chisq,deld,dmin,dmax)
595c
596c $$$$$ calls drdwr $$$$$
597c
598c   Lintry uses a Taylor series linearization of travel time to find a
599c   minimum in the dispersion function in the direction of the travel-
600c   time gradient.  Programmed on 12 December 1983 by R. Buland.
601c
602      double precision a,b,c,x,s,avdr
603c   Storage for the normal equations.
604      common /matrix/ a(3,4),b(3),c(4),x(3),s(3),avdr(3)
605c
606      m1=3-ihold
607      del2=deld
608      x0=0.
609      d0=chisq
610      x1=deld
611      d1=estlin(navsum,ihold,x1)
612      if(d0.lt.d1) go to 1
613      x2=2.*deld
614      d2=estlin(navsum,ihold,x2)
615 3    if(d1.lt.d2) go to 2
616      x0=x1
617      d0=d1
618      x1=x2
619      d1=d2
620      if(x1.ge.dmax) go to 13
621      del2=2.*del2
622      x2=amax1(x1+del2,dmax)
623      d2=estlin(navsum,ihold,x2)
624      go to 3
625 1    x2=x1
626      d2=d1
627      x1=.5*(x0+x2)
628      d1=estlin(navsum,ihold,x1)
629      if(x1.le.dmin) go to 14
630      if(d0.lt.d1) go to 1
631 2    xx=.5*(x0+x1)
632      if((x2-x0)/x1.le..15.or.x2-x0.le.dmin) go to 13
633      dd=estlin(navsum,ihold,xx)
634      if(dd.lt.d1) go to 4
635      x0=xx
636      d0=dd
637      xx=.5*(x1+x2)
638      dd=estlin(navsum,ihold,xx)
639      if(dd.lt.d1) go to 5
640      x2=xx
641      d2=dd
642      go to 2
643 4    x2=x1
644      d2=d1
645      x1=xx
646      d1=dd
647      go to 2
648 5    x0=x1
649      d0=d1
650      x1=xx
651      d1=dd
652      go to 2
653 14   if(d1.lt.d0) go to 13
654      x1=0.
655      d1=chisq
656 13   do 6 i=1,m1
657 6    x(i)=x1*x(i)
658      deld=x1
659      return
660      end
661
662      subroutine magsrt(n,rkey,ndex)
663c
664c $$$$$ calls no other routine $$$$$
665c
666c   Magsrt sorts the n elements of array rkey so that rkey(i),
667c   i = 1, 2, 3, ..., n are in asending order.  Auxillary integer*4
668c   array ndex is sorted in parallel.  Magsrt is a trivial modification
669c   of ACM algorithm 347:  "An efficient algorithm for sorting with
670c   minimal storage" by R. C. Singleton.  Array rkey (and ndex) is
671c   sorted in place in order n*alog2(n) operations.  Coded on
672c   8 March 1979 by R. Buland.  Modified to handle real*4 data on
673c   27 September 1983 by R. Buland.
674c
675      dimension rkey(1),ndex(1),il(10),iu(10)
676c   Note:  il and iu implement a stack containing the upper and
677c   lower limits of subsequences to be sorted independently.  A
678c   depth of k allows for n<=2**(k+1)-1.
679      if(n.le.1) return
680      r=.375
681      m=1
682      i=1
683      j=n
684c
685c   The first section interchanges low element i, middle element ij,
686c   and high element j so they are in order.
687c
688 5    if(i.ge.j) go to 70
689 10   k=i
690c   Use a floating point modification, r, of Singleton's bisection
691c   strategy (suggested by R. Peto in his verification of the
692c   algorithm for the ACM).
693      if(r.gt..58984375) go to 11
694      r=r+.0390625
695      go to 12
696 11   r=r-.21875
697 12   ij=i+(j-i)*r
698      if(rkey(i).le.rkey(ij)) go to 20
699      rit=rkey(ij)
700      rkey(ij)=rkey(i)
701      rkey(i)=rit
702      it=ndex(ij)
703      ndex(ij)=ndex(i)
704      ndex(i)=it
705 20   l=j
706      if(rkey(j).ge.rkey(ij)) go to 39
707      rit=rkey(ij)
708      rkey(ij)=rkey(j)
709      rkey(j)=rit
710      it=ndex(ij)
711      ndex(ij)=ndex(j)
712      ndex(j)=it
713      if(rkey(i).le.rkey(ij)) go to 39
714      rit=rkey(ij)
715      rkey(ij)=rkey(i)
716      rkey(i)=rit
717      it=ndex(ij)
718      ndex(ij)=ndex(i)
719      ndex(i)=it
720 39   tmpkey=rkey(ij)
721      go to 40
722c
723c   The second section continues this process.  K counts up from i and
724c   l down from j.  Each time the k element is bigger than the ij
725c   and the l element is less than the ij, then interchange the
726c   k and l elements.  This continues until k and l meet.
727c
728 30   rit=rkey(l)
729      rkey(l)=rkey(k)
730      rkey(k)=rit
731      it=ndex(l)
732      ndex(l)=ndex(k)
733      ndex(k)=it
734 40   l=l-1
735      if(rkey(l).gt.tmpkey) go to 40
736 50   k=k+1
737      if(rkey(k).lt.tmpkey) go to 50
738      if(k.le.l) go to 30
739c
740c   The third section considers the intervals i to l and k to j.  The
741c   larger interval is saved on the stack (il and iu) and the smaller
742c   is remapped into i and j for another shot at section one.
743c
744      if(l-i.le.j-k) go to 60
745      il(m)=i
746      iu(m)=l
747      i=k
748      m=m+1
749      go to 80
750 60   il(m)=k
751      iu(m)=j
752      j=l
753      m=m+1
754      go to 80
755c
756c   The fourth section pops elements off the stack (into i and j).  If
757c   necessary control is transfered back to section one for more
758c   interchange sorting.  If not we fall through to section five.  Note
759c   that the algorighm exits when the stack is empty.
760c
761 70   m=m-1
762      if(m.eq.0) return
763      i=il(m)
764      j=iu(m)
765 80   if(j-i.ge.11) go to 10
766      if(i.eq.1) go to 5
767      i=i-1
768c
769c   The fifth section is the end game.  Final sorting is accomplished
770c   (within each subsequence popped off the stack) by rippling out
771c   of order elements down to their proper positions.
772c
773 90   i=i+1
774      if(i.eq.j) go to 70
775      if(rkey(i).le.rkey(i+1)) go to 90
776      k=i
777      kk=k+1
778      rib=rkey(kk)
779      ib=ndex(kk)
780 100  rkey(kk)=rkey(k)
781      ndex(kk)=ndex(k)
782      kk=k
783      k=k-1
784      if(rib.lt.rkey(k)) go to 100
785      rkey(kk)=rib
786      ndex(kk)=ib
787      go to 90
788      end
789
790      subroutine nclose(nphase,ihold,navsum,nstuse,ifl,ierr)
791c
792c $$$$$ calls ellip, delaz, drdwr, lest, lpush,
793c             normeq $$$$$
794c
795c   Nclose performs all calculations for output which are not needed for
796c   the Gauss-Newton iteration itself.  This includes setting remaining
797c   free flags to use, computing distance-azimuth-residuals for
798c   residualed first arrivals, re-computing pP-P depths, and computing
799c   hypocentral parameter standard deviations.
800c   Nphase is the number of associated data, ihold is the number of
801c   constraint equations, navsum is the number of usable data, and ifl
802c   is a matrix decomposition flag.  If ifl is 1, the solution will be
803c   presumed to have changed since the last matrix decomposition.  If
804c   ifl is 2, the previous decomposition will be used.  If all goes
805c   well, nclose will return +1.
806c
807      include 'locprm.inc'
808        logical phasid,log
809c       integer stdin, stdout
810      double precision a,b,c,x,s,avdr
811c   Physical constants.
812      common /circle/ pi,rad,deg,radius
813c   Logical units.
814      common /unit/ min,mout,merr
815c   Storage for the normal equations.
816      common /matrix/ a(3,4),b(3),c(4),x(3),s(3),avdr(3)
817c   Hypocentral parameters.
818      include 'hypocm.inc'
819c   Station-reading parameters.
820      include 'pickcm.inc'
821c
822c       stdin = 5
823c       stdout= 6
824c
825c   Assume success.
826c
827      ierr=1
828c
829c   If all hypocentral parameters are held, there is no need to
830c   decompose the matrix as all standard deviations must be zero anyway.
831c   If the matrix was singular, there is nothing to decompose.
832c
833      if(ihold.eq.4.or.ng.eq.'b') go to 2
834      m=4-ihold
835c
836c   If there are insufficient equations for a solution, set a flag and
837c   skip the matrix decomposition.
838c
839      if(m-1.le.navsum.and.nstuse.ge.3) go to 1
840      ng='s'
841      go to 2
842c
843c   Compute the correlation matrix (A**-1) from the matrix
844c   decomposition
845c
8461     ierr=normeq(navsum,ihold)
847c
848c   Loop over all station-readings.
849c
8502     do 3 j = 1, nphase
851c
852c   Fetch data from station parameter storage.
853c
854      call drdwr(1,j,1,13)
855c
856c   Skip stations with bad stations or which have been noassed.
857c
858      if( pflg .eq. 'u' ) go to 3
859c
860c   Compute distance and azimuth.
861c
862      call delaz(olats,olatc,olons,olonc,slats,slatc,slons,slonc,
863        1                       delta,azim)
864c
865      log=phasid()
866c
867c   All done.  Push data back out to station parameter storage.
868c
869      call drdwr(2,j,1,13)
8703     continue
871c
872c   Clean up event level data.
873c
874c   Compute hypocentral parameter standard deviations.  If depth is
875c   held, all parameters are held, or the solution is bad, zero
876c   appropriate standard deviations.
877c
878      sedep=0.
879      if(ihold.eq.4.or.ng.eq.'s') se=0.
880      if(ihold.ne.4.and.ng.ne.'s'.and.ng.ne.'b') go to 10
881      selat=0.
882      selon=0.
883      setim=0.
884      call errelp(navsum,0,se)
885      return
886c   Convert the units of the correlation matrix to km**2.
887 10   rs=radius-odep
888      rl=rs*olats
889      m1=3-ihold
890      do 11 k=1,m1
891      a(1,k)=rs*a(1,k)
892      a(k,1)=rs*a(k,1)
893      a(2,k)=rl*a(2,k)
894 11   a(k,2)=rl*a(k,2)
895c   Otherwise, compute standard deviations from the diagonal elements of
896c   the correlation matrix.
897      setim=rsc1*se
898      selat=rsc1*dsqrt(dmax1(a(1,1),0d0))
899      selon=rsc1*dsqrt(dmax1(a(2,2),0d0))
900      if(ihold.eq.0) sedep=rsc1*dsqrt(dmax1(a(3,3),0d0))
901c     print *,'setim selat selon sedep',setim,selat,selon,sedep
902      call errelp(navsum,m,se)
903      return
904c   Flag a problem with the pP-P depth.
905 13   ierr=-1
906      return
907      end
908
909      function normeq(navsum,ihold)
910c
911c $$$$$ calls cholin, cholsl, and drdwr $$$$$
912c
913c   Normeq constructs and solves the normal equations to determine the
914c   next hypocentral step.  Nphases are processed.  All usable
915c   (not residualed or noassed) data are used.  M is the number of
916c   equations solved.  It will be 4 for a hypocentral step and 3 for an
917c   epicentral step.  If all goes well normeq returns 1.  Normeq returns
918c   -1 if the matrix is singular.  If ifl is 1, the system of equations
919c   is decomposed and  solved.  If ifl is 2, the correlation matrix is
920c   computed assuming that the matrix has already been decomposed.  If
921c   ifl is 3, the system of equations is decomposed, but not solved.
922c   Rewritten on 21 July 1983 by R. Buland from a fragment of a routine
923c   by E. R. Engdahl.
924c
925      double precision a,b,c,x,s,avdr,eps,fac
926c   Physical constants.
927      common /circle/ pi,rad,deg,radius
928c   Storage for the normal equations.
929      common /matrix/ a(3,4),b(3),c(4),x(3),s(3),avdr(3)
930c   Hypocentral parameters.
931      include 'hypocm.inc'
932c   Station-reading parameters.
933      include 'pickcm.inc'
934      data eps,fac0/1d-5,31.06341/
935c
936c   Assume success.
937      normeq=1
938      jfl=1
939      m1=3-ihold
940c   Zero out array storage.
941      do 1 k=1,m1
942      do 1 l=k,m1
943 1    a(l,k)=0d0
944c   Loop over the stations accumulating the normal matrix and data
945c   vector.
946      do 2 i=1,navsum
947c   Pop one of the sorted residuals.
948      ierr=lpop(i,j,res)
949c   Access station parameter storage.
950      call drdwr(1,j,1,5)
951      if(pflg.ne.'u') go to 2
952c   C is the contribution to the normal equations.
953      c(1)=deg*cos(azim)*dtdd-avdr(1)
954      c(2)=-deg*sin(azim)*dtdd*olats-avdr(2)
955      c(3)=dtdh-avdr(3)
956c     print *,'normeq',j,avdr,c
957c   Update the normal equations.
958      do 3 k=1,m1
959      do 3 l=k,m1
960 3    a(l,k)=a(l,k)+c(k)*c(l)
961 2    continue
962c     print 500,((a(l,k),k=1,m1),l=1,m1)
963c500  format(/(1x,1p3e15.4))
964c
965c   Solve the system.
966c
967c   Weight the system of equations by diagonal matrix S to improve
968c   conditioning.  Note that here as above, only the lower triangle of
969c   the normal matrix is needed by the Cholesky routine used below.
970      do 5 k=1,m1
971 5    s(k)=1d0/dmax1(dsqrt(a(k,k)),eps)
972      do 6 k=1,m1
973      do 6 l=k,m1
974 6    a(l,k)=s(l)*s(k)*a(l,k)
975c   Compute the Cholesky decomposition of the normal matrix.  The lower
976c   triangle of array A is unchanged.  The upper triangle of A is
977c   replaced by the Cholesky factor (offset by one column so everything
978c   will fit).
979 7    call cholin(3,m1,a,1,ierr)
980c   Branch if there is no error.
981      if(ierr.eq.0) go to 8
982c   Bail out if the fixup has already failed.
983      if(jfl.le.0) go to 13
984      jfl=-1
985c   Add eps times the identity matrix to the normal matrix to improve
986c   conditioning and try again.
987      do 9 k=1,m1
988 9    a(k,k)=a(k,k)+eps
989      go to 7
990c
991c   Transform the Cholesky factor in the upper triangle of A into the
992c   inverse of the normal matrix (the covariance matrix) in place.
993 8    call cholin(3,m1,a,-3,ierr)
994c   Remove the weighting from the covariance matrix, shift it over to
995c   replace the original a matrix, and reconstruct the lower triangle.
996      do 11 k=1,m1
997      do 11 l=k,m1
998      a(k,l)=s(k)*s(l)*a(k,l+1)
999 11   a(l,k)=a(k,l)
1000      return
1001c   Flag a singular solution.
1002 13   normeq=-1
1003      return
1004      end
1005
1006      function nstep(nsta,ihold,loop,iter,navsum,av,chisq,deld)
1007c
1008c $$$$$ calls normeq and stats $$$$$
1009c
1010c   Assuming that travel-time residuals and derivitives with respect to
1011c   distance and depth have already been computed, nstep finds the next
1012c   Gauss-Newton step, updates the hypocenter, and assesses the goodness
1013c   of fit of the new solution.  In the assessment phase, travel-time
1014c   residuals and derivitives for the new solution are computed.  Nsta
1015c   is the number of associated stations, ihold is the number of
1016c   constraint equations, loop is a character flag indicating the
1017c   iteration stage, iter is the current Gauss-Newton iteration number,
1018c   navsum is the number of usable arrivals, av is the average travel-
1019c   time residual, chisq is the sum of the squares of the residuals, and
1020c   deld is the length of the spatial step made in kilometers.  Nsta,
1021c   loop and iter0 are not changed.  Navsum is assumed to be valid on
1022c   entry, but is recomputed internally.  Av, chisq, and deld are
1023c   computed each time by nstep.  If nstep returns -1, depth was
1024c   constrained by nstep and ihold has been modified accordingly.
1025c   Otherwise, nstep returns +1.  Rewritten on 25 July 1983 by R. Buland
1026c   from a fragment of a routine by E. R. Engdahl.
1027c
1028      character*1 loop
1029      double precision a,b,c,x,s,avdr
1030      dimension hyprm(8),hypsav(4)
1031c   Logical units on input and output files
1032        common /unit/ min, mout, merr
1033c   Physical constants.
1034      common /circle/ pi,rad,deg,radius
1035c   Storage for the normal equations.
1036      common /matrix/ a(3,4),b(3),c(4),x(3),s(3),avdr(3)
1037c   Hypocentral parameters.
1038      include 'hypocm.inc'
1039c   Commands.
1040      include 'cmndcm.inc'
1041c   Back up hypocenter.
1042      common /hypfix/ hypbck(8)
1043c   Iteration control parameters.
1044      common /itcntl/ miter1,miter2,miter3,tol1,tol2,tol3,hmax1,hmax2,
1045     1 hmax3,tol,hmax,zmin,zmax,zshal,zdeep,rchi,roff,utol
1046      equivalence (osec,hyprm)
1047      data almost,eps,savz,addz,nexp/1.1,1e-4,-1.,0.,1/
1048c
1049c   Assume success.
1050      nstep=1
1051c   Find number of equations.
1052      m=4-ihold
1053      m1=m-1
1054c   If there aren't enough equations, bail out.
1055      if(m1.gt.navsum) return
1056c   Initialize output parameters.
1057      delh=0.
1058      delz=0.
1059      if(iter.le.1) cn=50.
1060c   Addz accounts for depth changes imposed outside nstep in the tag
1061c   line printed below.
1062      if(savz.ge.0.) addz=odep-savz
1063c   The damping factor, roff, is fiddled each time nstep is called.
1064c   This "randomization" avoids excessive iteration in certain
1065c   pathlogical cases.
1066      if(rchi.gt.roff) go to 33
1067      rchi=rchi+.0390625
1068      go to 34
1069 33   rchi=rchi-.21875
1070c
1071c   Use steepest descents to find the minimum of the piecewise linear
1072c   rank-sum penelty function.
1073c
1074c   Set the base hypocenter.
1075 34   do 3 i=1,m
1076 3    hypsav(i)=hyprm(i)
1077      cn=amax1(cn,2.*tol)
1078      x(1)=b(1)/(rad*radius)
1079      x(2)=b(2)/(rad*radius*olats)
1080      x(3)=b(3)
1081      call lintry(navsum,ihold,chisq,cn,tol,hmax)
1082c     print *,'x',(x(i),i=1,m1)
1083c   Set the horizontal and vertical steps (in km).
1084      delh=sqrt(amax1((sngl(x(1))*rad*radius)**2+
1085     1 (sngl(x(2))*rad*radius*olats)**2,0.))
1086      if(ihold.eq.0) delz=x(3)
1087      savdel=sqrt(amax1(delh*delh+delz*delz,0.))
1088c     print *,'del',delh,delz,savdel
1089c
1090c   Damp the step on maximum horizontal distance.
1091c
1092      if(delh.le.hmax) go to 8
1093      rat=hmax/delh
1094      do 9 i=1,m1
1095 9    x(i)=rat*x(i)
1096      delh=hmax
1097c
1098c   Update the depth.
1099c
1100 8    if(ihold.ne.0) go to 4
1101c
1102c   Trap airquakes.
1103      if(odep+x(3).ge.zmin) go to 10
1104      if(abs(odep-zmin).le.eps) go to 27
1105c   If we aren't yet at zmin, damp the step to just reach zmin.
1106      rat=(zmin-odep)/x(3)
1107      do 18 i=1,m1
1108 18   x(i)=rat*x(i)
1109      delh=rat*delh
1110c   If the damped step would be interpreted as convergence, constrain
1111c   depth.
1112      if(sqrt(amax1(delh*delh+sngl(x(3))*sngl(x(3)),0.))-tol) 20,20,11
1113c   If we are already at zmin, move along the z = zmin surface.
1114 27   if(delh.le.tol) go to 20
1115      x(3)=0d0
1116      go to 11
1117c
1118c   Trap lower mantle quakes.
1119 10   if(odep+x(3).le.zmax) go to 11
1120      if(abs(zmax-odep).le.eps) go to 28
1121c   If we aren't yet at zmax, damp step to just get there.
1122      rat=(zmax-odep)/x(3)
1123      do 19 i=1,m1
1124 19   x(i)=rat*x(i)
1125      delh=rat*delh
1126c   If the damped step would be interpreted as convergence, constrain
1127c   depth.
1128      if(sqrt(amax1(delh*delh+sngl(x(3))*sngl(x(3)),0.))-tol) 21,21,11
1129c   If we are already at zmax, move along the z = zmax surface.
1130 28   if(delh.le.tol) go to 21
1131      x(3)=0d0
1132 11   delz=x(3)
1133      go to 4
1134c
1135c   On a bad depth, reset it to be held at either zshal or zdeep as
1136c   apropriate and set flags.
1137c
1138 20   rat=zshal
1139      go to 22
1140 21   rat=zdeep
1141 22   if(rat.ne.hypbck(4)) go to 29
1142c   If we have already converged at the held depth, reset the hypocenter
1143c   to the known solution.
1144      delh=sqrt(amax1(((hypbck(2)-hypsav(2))*rad*radius)**2+((hypbck(3)-
1145     1 hypsav(3))*rad*radius*olats)**2,0.))
1146      delz=hypbck(4)-hypsav(4)
1147      do 30 i=1,8
1148 30   hyprm(i)=hypbck(i)
1149      nstep=0
1150      go to 31
1151c   Otherwise, use the current epicenter, reset depth to the desired
1152c   value, and find an approximately compatible origin time.
1153 29   do 24 i=1,m
1154 24   hyprm(i)=hypsav(i)
1155      delh=0.
1156      delz=rat-odep
1157      osec=osec+.1*delz
1158      odep=rat
1159      nstep=-1
1160c   Set the depth held flags.
1161 31   ihold=1
1162      m=3
1163      m1=2
1164      hmax=hmax1
1165      dep ='c'
1166c   Go recompute residuals and derivitives for the next iteration.
1167      go to 1
1168c
1169c   Update the hypocentral parameters.
1170c
1171 4    do 12 i=1,m1
1172 12   hyprm(i+1)=hypsav(i+1)+x(i)
1173      odep=amin1(amax1(odep,zmin),zmax)
1174c   Guarantee legal co-latitude and longitude.
1175      if(olat.ge.0..and.olat.le.180.) go to 15
1176      olon=olon+180.
1177      if(olat.lt.0.) olat=abs(olat)
1178      if(olat.gt.180.) olat=360.-olat
1179 15   if(olon.lt.0.) olon=olon+360.
1180 17   if(olon.le.360.) go to 16
1181      olon=olon-360.
1182      go to 17
1183c   Calculate epicentral sines and cosines.
1184 16   olats=sin(olat*rad)
1185      olatc=cos(olat*rad)
1186      olons=sin(olon*rad)
1187      olonc=cos(olon*rad)
1188c
1189c   Recompute travel-time residuals, derivitives, and chi**2.
1190c
1191 1    call stats(nsta,ihold,navsum,av,chi2)
1192      hyprm(1)=hyprm(1)+av
1193c   Deld is the convergence parameter.
1194      deld=sqrt(delh*delh+delz*delz)
1195c
1196c   If chi**2 has decreased or damping has been unsuccessful, get out.
1197c
1198      if(chi2.lt.chisq.or.nstep.le.0.or.ng.ne.' ') go to 6
1199c   If we can't damp the step any more, bail out.
1200      if(rchi*deld.le.tol) go to 25
1201c   If chi**2 has increased, damp the step length by factor rchi.
1202      do 7 i=1,m1
1203 7    x(i)=rchi*x(i)
1204      delh=rchi*delh
1205      delz=rchi*delz
1206      go to 4
1207c   Trap an unstable solution on a clipped depth.
1208 25   if(abs(odep-zmin).le..01.and.ihold.eq.0) go to 20
1209      if(abs(odep-zmax).le..01.and.ihold.eq.0) go to 21
1210c   Flag an unstable solution.  Go around one more time to recover the
1211c   starting point which can't be improved.
1212      ng='u'
1213      if(savdel.le.utol) ng='q'
1214      if(chi2.le.almost*chisq.and.deld.le.tol) ng='o'
1215c   Zero the step and go around again to recompute the original
1216c   residuals.
1217      do 26 i=1,m1
1218 26   x(i)=0d0
1219      delh=0.
1220      delz=0.
1221      go to 4
1222c
1223c   Compute the event standard deviation.
1224 6    chisq=chi2
1225      rms=0.0
1226      if(navsum.gt.0) rms=chisq/navsum
1227c
1228c   Write a summary line to the standard output.
1229      rlat=90.-olat
1230      rlon=olon
1231      if(rlon.gt.180.) rlon=rlon-360.
1232      delz=delz+addz
1233      savz=odep
1234      write(merr,100,iostat=ios)loop,iter,navsum,rlat,rlon,odep,
1235     1 delh,delz,rms,ng,dep
1236 100  format('nstep: ',a1,i3,i5,f9.3,f10.3,f7.1,' delh=',f6.1,
1237     1       ' delz=',f7.1,' rms=',f10.4,1x,'ng= ',a1,'dep= ',a1)
1238 13   return
1239      end
1240
1241      subroutine scopt(n,sc)
1242c
1243c   Generate scores for n data.  Note that the scores don't depend on 
1244c   the data value, only on it's position.
1245c
1246      dimension sc(1),p(29),s(29)
1247      data p/0.,.1375,.1625,.1875,.2125,.2375,.2625,.2875,.3125,.3375,
1248     1 .3625,.3875,.4125,.4375,.4625,.4875,.5125,.5375,.5625,.5875,
1249     2 .6125,.6375,.6625,.6875,.7125,.7375,.7625,.7875,1./
1250      data s/0.0775,0.0775,0.1546,0.5328,0.8679,1.1714,1.4542,1.7266,
1251     1 1.9987,2.2802,2.5803,2.9068,3.2657,3.6603,4.0912,4.5554,5.0470,
1252     2 5.5572,6.0754,6.5906,7.0919,7.5702,8.0194,8.4365,8.8223,9.1812,
1253     3 9.5207,9.5974,9.5974/
1254c
1255      cn=1./(n+1)
1256      k1=1
1257      k=2
1258      av=0.
1259      do 1 i=1,n
1260      x=cn*i
1261 3    if(x.le.p(k)) go to 2
1262      k1=k
1263      k=k+1
1264      go to 3
1265 2    sc(i)=(x-p(k1))*(s(k)-s(k1))/(p(k)-p(k1))+s(k1)
1266 1    av=av+sc(i)
1267      av=av/n
1268      do 4 i=1,n
1269 4    sc(i)=sc(i)-av
1270      return
1271      end
1272
1273      subroutine setup(nphase,ihold,navsum,av,chisq)
1274c
1275c $$$$$ calls ellip, delaz, drdwr $$$$$
1276c
1277c   Setup initializes the location process.  This involves calculating
1278c   distance and azimuth, travel-time and derivitives, setting the phase
1279c   flag, and implementing the pkp commands.  nphase is the number of
1280c   associated phases, ihold is the
1281c   number of constraints on the solution, navsum is the number of
1282c   usable first arrivals, av is the average travel-time residual, and
1283c   chisq is the sum of the squared residuals.  Navsum, av,and chisq are
1284c   computed by setup.  Rewritten on 22 July 1983 by R. Buland from a
1285c   fragment of a routine by E. R. Engdahl.
1286c
1287      include 'locprm.inc'
1288        logical phasid
1289c       integer stdin, stdout
1290        common /unit/ min, mout, merr
1291c   Physical constants.
1292      common /circle/ pi,rad,deg,radius
1293c   Hypocentral parameters.
1294      include 'hypocm.inc'
1295c   Event level analyst commands.
1296      include 'cmndcm.inc'
1297c   Station-reading parameters.
1298      include 'pickcm.inc'
1299c   Scores.
1300      common /sccom/ scgen(maxph)
1301c
1302c       stdout = 6
1303c       stdin = 5
1304c
1305c   Initialize counter
1306c
1307      navsum=0
1308c
1309c   Loop over all station-readings.
1310c
1311        jj = 0
1312      do 2 j=1,nphase
1313c
1314c   Fetch data from station parameter storage.
1315c
1316      call drdwr(1,j,6,11)
1317c
1318c   Skip stations which can't be used in the location.
1319c
1320        if( pflg .ne. 'u' ) go to 2
1321c
1322c   Compute distance, azimuth and traveltime of phase.
1323c
1324      call delaz(olats,olatc,olons,olonc,slats,slatc,slons,slonc,
1325        1                       delta,azim)
1326c
1327c   Apply the DRES commands.
1328      do k=1,5
1329        if(dres(k).eq.'r'.and.delta.ge.dellim(1,k).and.delta.le.
1330     1   dellim(2,k)) then
1331          pflg=' '
1332          go to 8
1333        endif
1334      enddo
1335c
1336c   ID the phases and compute residuals.
1337      if(phasid()) then
1338c
1339c   Apply the RRES command.
1340                if(rres.eq.'r'.and.(res.lt.reslim(1).or.res.gt.reslim(2))) then
1341                  pflg=' '
1342                  go to 8
1343                endif
1344c   Push the residual for statistical purposes.
1345                ierr=lpush(navsum,j,res)
1346        endif
1347c
1348c   Restore data to station parameter storage.
1349c
1350 8    call drdwr(2,j,1,5)
1351c       call drdwr(2,j,22,22)
1352 2    continue
1353c
1354      call scopt(navsum,scgen)
1355c
1356c   Compute all statistical parameters.
1357c
1358      ierr=lsrt(navsum,av,se,chisq,'TT')
1359      call adder(navsum,ihold)
1360c       write(merr,*)'setup: ',' nphase= ',nphase,' navsum= ',navsum
1361      return
1362      end
1363
1364      subroutine stats(nphase,ihold,navsum,av,chisq)
1365c
1366c $$$$$ calls ellip, delaz, drdwr $$$$$
1367c
1368c   Stats computes travel-time residuals and derivitives for the current
1369c   solution.  Nphase is the number of associated stations, ihold is the
1370c   number of constraint equations, navsum is the number of usable
1371c   arrivals, av is the average travel-time residual, and chisq is the
1372c   sum of the squared residuals.  Nphase and ihold are not changed.
1373c   Navsum, av, and chisq are computed by stats.  Rewritten on
1374c   5 August 1983 by R. Buland from a fragment of a routine by
1375c   E. R. Engdahl.
1376c
1377      include 'locprm.inc'
1378        logical phasid
1379c       integer stdout,stdin
1380c   Physical constants.
1381      common /circle/ pi,rad,deg,radius
1382c   Logical units.
1383      common /unit/ min,mout,merr
1384c   Hypocentral parameters.
1385      include 'hypocm.inc'
1386c   Station-reading parameters.
1387      include 'pickcm.inc'
1388c
1389c       stdout = 6
1390c       stdin = 5
1391c
1392c   Initialize statistical counters.
1393c
1394      navsum=0
1395c
1396c   Loop over all station-readings.
1397      do 1 j=1,nphase
1398c   Fetch data from station parameter storage.
1399      call drdwr(1,j,6,11)
1400c   Don't bother with bad stations or noassed or residualed readings.
1401        if( pflg .ne. 'u' ) go to 1
1402c   Get the new distance and azimuth.
1403      call delaz(olats,olatc,olons,olonc,slats,slatc,slons,slonc,
1404        1                       delta,azim)
1405c
1406c   ID the phase and compute the residual.
1407        if(phasid()) then
1408c
1409c   Update the statistical sums.
1410                ierr=lpush(navsum,j,res)
1411        endif
1412c
1413c   Save computed data in station parameter storage.
1414c
1415      call drdwr(2,j,1,5)
14161       continue
1417c
1418c   Compute the average travel-time residual and chi**2.
1419c
1420      ierr=lsrt(navsum,av,se,chisq,'TT')
1421      call adder(navsum,ihold)
1422      return
1423      end
1424
1425      function stirf(u,t0,t1,t2)
1426c
1427c   Compute an asymptotic approximation to the Stirling function.
1428c
1429      stirf=t1+0.5*u*(t2-t0)+0.5*(u**2)*(t2-2.0*t1+t0)
1430      return
1431      end
1432 
1433      function xmed(n,a,sd)
1434c
1435c $$$$$ calls only library routines $$$$$
1436c
1437c   Given the n numbers, a(1) <= a(2) <= ... <= a(n), xmed returns
1438c   robust estimates of center and spread of the data about the center.
1439c   The array a is not altered.  Programmed on 9 December 1983 by
1440c   R. Buland.
1441c
1442      dimension a(1)
1443      data xno,xne/1.482580,.7412898/
1444c   See if there is enough data to work with.
1445      if(n-1)13,12,14
1446c   No, return.
1447 13   xmed=0.
1448      sd=0.
1449      return
1450c   One datum is a special case.
1451 12   xmed=a(1)
1452      sd=0.
1453      return
1454c   N > 1, the median is a robust estimate of center.
1455 14   m=n/2
1456      md=mod(n,2)
1457      if(md.le.0) xmed=.5*(a(m)+a(m+1))
1458      if(md.gt.0) xmed=a(m+1)
1459c   Find the smallest positive value.
1460      do 1 i=1,n
1461      if(a(i)-xmed.ge.0.) go to 2
1462 1    continue
1463      i=n
1464c   Find the smallest absolute value.
1465 2    if(i.gt.1.and.abs(a(i-1)-xmed).lt.abs(a(i)-xmed)) i=i-1
1466      k=i
1467      j=1
1468      if(j.ge.m) go to 15
1469c   Implicitly order the series in order of increasing absolute value by
1470c   sucessively moving i up or k down until the m th term is found.
1471 6    if(k.le.1) go to 3
1472      if(i.ge.n) go to 4
1473      if(abs(a(i+1)-xmed).gt.abs(a(k-1)-xmed)) go to 5
1474c   Move i up.
1475      i=i+1
1476      j=j+1
1477      if(j.lt.m) go to 6
1478c   The m th term is found, mark it.
1479 15   l1=i
1480      go to 8
1481c   Move k down.
1482 5    k=k-1
1483      j=j+1
1484      if(j.lt.m) go to 6
1485c   The m th term is found, mark it.
1486      l1=k
1487c   Find and mark the m+1 th term.
1488 8    if(k.le.1) go to 9
1489      if(i.ge.n) go to 10
1490      if(abs(a(i+1)-xmed)-abs(a(k-1)-xmed))9,9,10
1491 9    l2=i+1
1492      go to 11
1493 10   l2=k-1
1494      go to 11
1495c   We have run out of negative terms, the m th term must be up.
1496 3    l1=i+m-j
1497      l2=l1+1
1498      go to 11
1499c   We have run out of positive terms, the m th term must be down.
1500 4    l1=k-(m-j)
1501      l2=l1-1
1502c   The spread is the normalized median of the absolute series.
1503 11   if(md.le.0) sd=xne*(abs(a(l1)-xmed)+abs(a(l2)-xmed))
1504      if(md.gt.0) sd=xno*abs(a(l2)-xmed)
1505      return
1506      end
1507
1508c   Withers modified declaration of phasid function to prefent f90
1509c   from thinking its main.  (gotta love .for)
1510c     logical function phasid
1511      logical function phasid()
1512
1513c
1514c   Re-identify the current phase and return its travel-time residual.
1515c
1516      implicit none
1517        real depmin,cn
1518        parameter(depmin=1.,cn=5.)
1519c   Miscellaneous constants.
1520      include 'locprm.inc'
1521      character*8 phcd1(maxtt), tmpcod
1522        character*20 prmmod
1523        logical first
1524        integer*4 n,i,j,ntrbl,mprm,ln
1525c       integer*4 min,mout,merr
1526      real*4 zs,tt1(maxtt),dtdd1(maxtt),dtdh1(maxtt),dddp1(maxtt),
1527     1 ellip,elcor,pi,rad,deg,radius,tcor,ecor,zs0,usrc(2),
1528     2 travt(maxtt),win(maxtt),spd,amp
1529c   Physical constants.
1530      common /circle/ pi,rad,deg,radius
1531c   Logical units.
1532c     common /unit/ min,mout,merr
1533c   Hypocentral parameters.
1534      include 'hypocm.inc'
1535c   Station-reading parameters.
1536      include 'pickcm.inc'
1537        data first,zs0/.true.,-1./
1538c
1539c   The first time through initialize statistical parameters and the
1540c   ellipticity routine.
1541        if(first) then
1542                first=.false.
1543                call freeunit(mprm)
1544                call getprm(mprm,prmmod)
1545                close(mprm)
1546                tcor=ellip('P       ',100.0,zs,glat,45.0)
1547c               print *,'Phasid: ellip set'
1548        endif
1549c
1550c   Initialize depth for the calculation of the traveltimes and partial
1551c      derivatives
1552        zs = amax1(odep,depmin)
1553        if(zs.ne.zs0) then
1554                zs0=zs
1555                call depset(zs,usrc)
1556c               print *,'Phasid: depset - zs glat =',zs,glat
1557        endif
1558c
1559c   Compute travel times of candidate phases.
1560c       print *,'Phasid: scode depth delta azim elev = ',sta,zs,delta,
1561c       1 azim*deg,elev
1562      call trtm(delta,maxtt,n,tt1,dtdd1,dtdh1,dddp1,phcd1)
1563c
1564c   Compute and apply travel-time corrections.
1565        do i = 1, n
1566c               print *,'Trtm: phcd tt dtdd dtdh = ',phcd1(i),tt1(i),dtdd1(i),
1567c       1        dtdh1(i)
1568                ecor=elcor(phcd1(i),elev,dtdd1(i))
1569                if((phcd1(i).eq.'Pg'.or.phcd1(i).eq.'Pb'.or.phcd1(i).eq.'Pn'.or.
1570        1        phcd1(i).eq.'P').and.dtdh1(i).gt.0.) then
1571                        tcor=ellip('Pup',delta,zs,glat,azim*deg)
1572c                       print *,'Pup'
1573                else if((phcd1(i).eq.'Sg'.or.phcd1(i).eq.'Sb'.or.
1574        1        phcd1(i).eq.'Sn'.or.phcd1(i).eq.'S').and.dtdh1(i).gt.0.) then
1575                        tcor=ellip('Sup',delta,zs,glat,azim*deg)
1576c                       print *,'Sup'
1577                else
1578                        tcor=ellip(phcd1(i),delta,zs,glat,azim*deg)
1579                endif
1580                travt(i)=tt1(i)+(tcor+ecor)+osec
1581c   Base the association window on half the spread.
1582          call getstt(phcd1(i),delta,dtdd1(i),spd,amp)
1583          win(i)=cn*spd
1584c               print *,'Phasid: phase tcor ecor travt win = ',
1585c       1        phcd1(i),tcor,ecor,travt(i),win(i)
1586        enddo
1587c       pause
1588c
1589c   Identify the phase using the Chicxulub algorithm.
1590c
1591c   Handle a generic P (phase class 'p').
1592      if(phase.eq.'P') then
1593        i=0
1594        res=1e30
1595        do j=1,n
1596          if(phcd1(j).eq.'Pg'.or.phcd1(j).eq.'Pb'.or.phcd1(j).eq.
1597        1        'Pn'.or.phcd1(j).eq.'P'.or.phcd1(j).eq.'Pdif'.or.
1598     2     phcd1(j)(1:3).eq.'PKP'.or.phcd1(j).eq.'PKiKP') then
1599            if(abs(trvtim-travt(j)).lt.res) then
1600              i=j
1601              res=abs(trvtim-travt(j))
1602            endif
1603          endif
1604        enddo
1605c   Handle a generic S (phase class 's').
1606      else if(phase.eq.'S') then
1607        i=0
1608        res=1e30
1609        do j=1,n
1610          if(phcd1(j).eq.'Sg'.or.phcd1(j).eq.'Sb'.or.phcd1(j).eq.
1611        1        'Sn'.or.phcd1(j).eq.'S') then
1612            if(abs(trvtim-travt(j)).lt.res) then
1613              i=j
1614              res=abs(trvtim-travt(j))
1615            endif
1616          endif
1617        enddo
1618c   Handle a generic PKP (phase class 'k').
1619      else if(phase.eq.'PKP') then
1620        i=0
1621        res=1e30
1622        do j=1,n
1623          if(phcd1(j)(1:3).eq.'PKP'.or.phcd1(j).eq.'PKiKP') then
1624            if(abs(trvtim-travt(j)).lt.res) then
1625              i=j
1626              res=abs(trvtim-travt(j))
1627            endif
1628          endif
1629        enddo
1630c   Handle a generic PP (phase class '2').
1631      else if(phase.eq.'PP') then
1632        i=0
1633        res=1e30
1634        do j=1,n
1635          if(phcd1(j).eq.'PgPg'.or.phcd1(j).eq.'PbPb'.or.phcd1(j).eq.
1636     1     'PnPn'.or.phcd1(j).eq.'PP') then
1637            if(abs(trvtim-travt(j)).lt.res) then
1638              i=j
1639              res=abs(trvtim-travt(j))
1640            endif
1641          endif
1642        enddo
1643c   Handle a generic core phase (phase class 'c').
1644      else if(phase.eq.'pPKP'.or.phase.eq.'sPKP'.or.phase.eq.'SKP'.or.
1645     1 phase.eq.'PKKP'.or.phase.eq.'P''P'''.or.phase.eq.'SKS') then
1646        ln=ntrbl(phase)
1647        i=0
1648        res=1e30
1649        do j=1,n
1650          if(phcd1(j)(1:ln).eq.phase) then
1651            if(abs(trvtim-travt(j)).lt.res) then
1652              i=j
1653              res=abs(trvtim-travt(j))
1654            endif
1655          endif
1656        enddo
1657c   Handle an unidentified phase (phase class ' ').
1658      else if(phase.eq.' ') then
1659        i=0
1660        res=1e30
1661        do j=1,n
1662c         if(abs(trvtim-travt(j)).lt.10.) print *,'phase res win = ',
1663c    1     phcd1(j),trvtim-travt(j),win(j),delta
1664          if(abs(trvtim-travt(j)).le.win(j)) then
1665            if(abs(trvtim-travt(j)).lt.res) then
1666              i=j
1667              res=abs(trvtim-travt(j))
1668            endif
1669          endif
1670        enddo
1671c       if(i.gt.0) print *,'Select ',phcd1(i)
1672c   Otherwise, look for an exact match (phase class 'e').
1673      else
1674        do i=1,n
1675          if(phase.eq.phcd1(i)) go to 1
1676        enddo
1677        i=0
1678      endif
1679c
1680c   Take care of special cases.
1681      if(i.eq.0.and.phase.ne.' ') then
1682c   If an explicit P crustal phase doesn't exist, let it float.
1683        if(phase.eq.'Pg'.or.phase.eq.'Pb'.or.phase.eq.'Pn') then
1684          i=0
1685          res=1e30
1686          do j=1,n
1687            if(phcd1(j).eq.'Pg'.or.phcd1(j).eq.'Pb'.or.phcd1(j).eq.
1688     1       'Pn') then
1689              if(abs(trvtim-travt(j)).lt.res.and.
1690        1                abs(trvtim-travt(j)).lt.win(j)) then
1691                i=j
1692                res=abs(trvtim-travt(j))
1693              endif
1694            endif
1695          enddo
1696c   If an explicit S crustal phase doesn't exist, let it float.
1697        else if(phase.eq.'Sg'.or.phase.eq.'Sb'.or.phase.eq.'Sn') then
1698          i=0
1699          res=1e30
1700          do j=1,n
1701            if(phcd1(j).eq.'Sg'.or.phcd1(j).eq.'Sb'.or.phcd1(j).eq.
1702     1       'Sn') then
1703              if(abs(trvtim-travt(j)).lt.res.and.
1704        1                abs(trvtim-travt(j)).lt.win(j)) then
1705                i=j
1706                res=abs(trvtim-travt(j))
1707              endif
1708            endif
1709          enddo
1710c   Allow natural extensions of P and PKP phases.
1711        else
1712          tmpcod=' '
1713          if(phase.eq.'Pg') then
1714            tmpcod='Pb'
1715          else if(phase.eq.'Pn') then
1716            tmpcod='P'
1717          else if(phase.eq.'Sg') then
1718            tmpcod='Sb'
1719          else if(phase.eq.'Sn') then
1720            tmpcod='S'
1721          else if(phase.eq.'PKPdf') then
1722            tmpcod='PKiKP'
1723          else if(phase.eq.'PKiKP') then
1724            tmpcod='PKPdf'
1725          else if(phase.eq.'PKPbc') then
1726            tmpcod='PKPdif'
1727          else if(phase.eq.'PKPdif') then
1728            tmpcod='PKPbc'
1729          else if(phase.eq.'Pdif') then
1730            tmpcod='P'
1731          endif
1732          if(tmpcod.ne.' ') then
1733            do i=1,n
1734              if(tmpcod.eq.phcd1(i)) go to 1
1735            enddo
1736          endif
1737        endif
1738      endif
1739c
1740c   Got it.  Compute the residual.
1741c
17421       if(i.gt.0) then
1743c               print *,'Phasid: phase phcd res = ',phase,phcd1(i),
1744c       1        trvtim-travt(i)
1745                phase=phcd1(i)
1746                res=trvtim-travt(i)
1747                dtdd=dtdd1(i)
1748                dtdh=dtdh1(i)
1749                phasid=.true.
1750        else
1751c
1752c Phase not found for that distance and depth
1753c
1754                pflg = ' '
1755                res=0.
1756                phasid=.false.
1757        endif
1758        return
1759        end
1760
1761        real*4 function elcor(phase,elev,dtdd)
1762c
1763c   Return the elevation correction.
1764c
1765        implicit none
1766      real vp,vs,kmpd
1767      parameter(vp=5.80,vs=3.46,kmpd=111.19493)
1768        character*(*) phase
1769        integer*4 ntrbl,j
1770        real*4 elev,dtdd
1771c
1772c   Figure out the arriving wave type.
1773        do j=ntrbl(phase),1,-1
1774                if(phase(j:j).eq.'P') then
1775c   This phase arrives as a P.  Do the correction using Vp.
1776                        elcor=(elev/vp)*sqrt(abs(1.-
1777     1           amin1((vp*dtdd/kmpd)**2,1.)))
1778c                       print *,'Elcor: P corr = ',phase,elcor
1779                        return
1780                endif
1781                if(phase(j:j).eq.'S') then
1782c   This phase arrives as an S.  Do the correction using Vs.
1783                        elcor=(elev/vs)*sqrt(abs(1.-
1784     1           amin1((vs*dtdd/kmpd)**2,1.)))
1785c                       print *,'Elcor: S corr = ',phase,elcor
1786                        return
1787                endif
1788        enddo
1789      return
1790        end
1791
1792      subroutine getprm(nprm,modnam)
1793c
1794c   Read the parameter file.
1795c
1796      implicit none
1797      character*(*) modnam
1798      character*20 flnm
1799      integer*4 nprm,i,j,k,ntrbl
1800      include 'resprm.inc'
1801      data flnm/'res_ac.prm'/
1802      character*1001 dirname
1803      common /my_dir/dirname
1804c
1805c    Withers modified to be as other open statements in rayloc_ew
1806c     open(nprm,file=flnm,access='sequential',form='formatted',
1807c    1 status='old',readonly,shared)
1808c
1809c    Note that this param file is new to this version of rayloc_ew
1810c    (as of 20071205).  It should live in the same place as the model
1811c    (e.g. ak135.*).  Withers
1812c
1813c       Make sure we use dirname if one is provided
1814      if(dirname(1:2).eq.'. ') then
1815        open(nprm,file=flnm,access='sequential',form='formatted',
1816     >       status='old')
1817      else
1818        open(nprm,file=dirname(1:ntrbl(dirname))//'/'//flnm,
1819     >       access='sequential',form='formatted',status='old')
1820      endif
1821
1822      read(nprm,101)modnam
1823 101  format(1x,a)
1824c     print *,'Getprm:  modnam = ',modnam
1825c
1826      do i=1,mxbr
1827         read(nprm,102,end=1)phnm(i),ms(i),me(i),xm(i),cut(i)
1828 102     format(1x,a8,2i8,f8.1,f8.0)
1829c        write(*,102)phnm(i),ms(i),me(i),xm(i),cut(i)
1830         do k=1,3
1831            read(nprm,103)mo(k,i),(prm(j,k,i),j=1,mo(k,i))
1832 103        format(6x,i5,1p4e15.7:/(11x,1p4e15.7))
1833c           write(*,103)mo(k,i),(prm(j,k,i),j=1,mo(k,i))
1834c             pause
1835         enddo
1836      enddo
1837      i=mxbr+1
1838c
1839 1    nbr=i-1
1840      if(nbr.lt.mxbr) phnm(nbr+1)=' '
1841c     print *,'Getprm:  modnam nbr = ',modnam(1:ntrbl(modnam)),nbr
1842      close(nprm)
1843      return
1844      end
1845
1846      subroutine getstt(phcd,delta,dtdd,spd,amp)
1847c
1848c   Given a phase code, phcd, source-receiver distance, delta (in deg),
1849c   and dt/d(delta), dtdd (in s/deg), Getstt returns the observed
1850c   spread, spd (sec), and relative amplitude, amp of the phase at that
1851c   distance.  Note that the spread is the robust equivalent of one
1852c   standard deviation.  If statistics aren't available, spd and amp are 
1853c   set to zero.
1854c
1855      implicit none
1856      character*(*) phcd
1857      integer*4 j,iupcor,margin
1858      real*4 delta,dtdd,spd,amp,xcor,tcor,x0,evlply
1859      include 'resprm.inc'
1860c
1861c   Find the desired phase in the statistics common.
1862      do j=1,nbr
1863        if(phcd.eq.phnm(j)) then
1864c   Figure the correction to surface focus.
1865          if(iupcor(phcd(1:1),dtdd,xcor,tcor).le.0) then
1866c           print *,'Iupcor failed on phase ',phcd
1867            xcor=0.
1868          endif
1869c         print *,'Getstt:  phcd delta xcor = ',phcd,delta,xcor
1870c   Correct the phase to surface focus.
1871          if(phcd(1:1).eq.'p'.or.phcd(1:1).eq.'s') then
1872c   For depth phases, the deeper, the closer.
1873            x0=delta-xcor
1874          else
1875c   For other phases, the deeper, the farther.
1876            x0=delta+xcor
1877          endif
1878c   If statistics are defined, do the interpolation.
1879          if(phnm(j).ne.phnm(min0(j+1,mxbr)).or.j.ge.nbr) then
1880            margin=15
1881          else
1882            margin=0
1883          endif
1884          if(x0.ge.ms(j)-16.and.x0.le.me(j)+margin) then
1885            x0=x0-xm(j)
1886            spd=amin1(amax1(evlply(x0,mo(3,j),prm(1,3,j)),.7),10.)
1887            amp=amax1(evlply(x0,mo(1,j),prm(1,1,j)),0.)
1888c           print *,'Getstt:  x0 spd amp =',x0,spd,amp,ms(j),me(j)
1889c           pause
1890            return
1891          endif
1892        endif
1893      enddo
1894c   If we don't find anything, return zero.
1895      spd=0.
1896      amp=0.
1897      return
1898      end
1899
1900      function evlply(x,m,b)
1901c
1902c   Evlply evaluates the polynomial.
1903c
1904      dimension b(1)
1905c
1906      xx=1.
1907      evlply=b(1)
1908      do j=2,m
1909         xx=xx*x
1910         evlply=evlply+xx*b(j)
1911      enddo
1912      return
1913      end
Note: See TracBrowser for help on using the repository browser.