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

Revision 1669, 47.4 KB checked in by friberg, 15 years ago (diff)

First commit of rayloc_ew in EW-CENTRAL CVS

  • 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*21 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      character*8 phcd1(maxtt),ophase, ph
809c       integer stdin, stdout
810      real*4 tt1(maxtt),dtdd1(maxtt),dtdh1(maxtt),dddp1(maxtt),usrc(2),
811     1 ellip, elcor
812      double precision a,b,c,x,s,avdr
813c   Physical constants.
814      common /circle/ pi,rad,deg,radius
815c   Logical units.
816      common /unit/ min,mout,merr
817c   Storage for the normal equations.
818      common /matrix/ a(3,4),b(3),c(4),x(3),s(3),avdr(3)
819c   Hypocentral parameters.
820      include 'hypocm.inc'
821c   Station-reading parameters.
822      include 'pickcm.inc'
823        data dtds/19.96/
824c
825c       stdin = 5
826c       stdout= 6
827c
828c   Assume success.
829c
830      ierr=1
831c
832c   If all hypocentral parameters are held, there is no need to
833c   decompose the matrix as all standard deviations must be zero anyway.
834c   If the matrix was singular, there is nothing to decompose.
835c
836      if(ihold.eq.4.or.ng.eq.'b') go to 2
837      m=4-ihold
838c
839c   If there are insufficient equations for a solution, set a flag and
840c   skip the matrix decomposition.
841c
842      if(m-1.le.navsum.and.nstuse.ge.3) go to 1
843      ng='s'
844      go to 2
845c
846c   Compute the correlation matrix (A**-1) from the matrix
847c   decomposition
848c
8491     ierr=normeq(navsum,ihold)
850c
851c Initialize the depth
852c
8532       if (odep .le. 0.0) odep = 1.0
854        zs = odep
855      call depset(zs,usrc)
856c
857c   Loop over all station-readings.
858c
859      do 3 j = 1, nphase
860c
861c   Fetch data from station parameter storage.
862c
863      call drdwr(1,j,1,13)
864c
865c   Skip stations with bad stations or which have been noassed.
866c
867      if( pflg .eq. 'u' ) go to 3
868c
869c   Compute distance and azimuth.
870c
871      call delaz(olats,olatc,olons,olonc,slats,slatc,slons,slonc,
872        1                       delta,azim)
873c
874        n = 0
875      call trtm(delta,maxtt,n,tt1,dtdd1,dtdh1,dddp1,phcd1)
876c
877c Find the computed traveltimes and partials for phases other than P and S
878c
879        do 111 i = 1, n
880                ph = phcd1(i)
881                ilen = ntrbl(phase)
882                if(phase(1:ilen) .eq. ph(1:ilen)) then
883                        travt = tt1(i)
884                        dtdd = dtdd1(i)
885                        dtdh = dtdh1(i)
886                        ophase = phase
887                        phase = ph
888                        go to 100
889                endif
890111     continue
891c
892c Phase not found for that distance and depth
893c
894        pflg = ' '
895        res=0.
896        go to 5
897c
898100     continue
899c
900c This was the original call to get the ellipticity and elevation correction,
901c   which is replaced by subroutine used by Engdahl.  Correction for the
902c   elevation was taken out of subroutine datum
903c     call datum(elev,dtdd,olat,delta,azim*deg,tcor,slatc)
904c
905        azimuth = azim*deg
906        tcor=ellip(phase,delta,zs,glat,azimuth)
907        if(dtdd.gt.dtds) dtdd=dtds
908      elcor=(elev/5.57)*sqrt (1.0-(5.57*dtdd/(radius*rad))**2)
909      res=trvtim-osec-travt-(tcor+elcor)
910c
911c
912c       write(merr,333)j,sta,ophase,ph,pflg,delta,azimuth,
913c       2                               trvtim,res,dtdd,dtdh
914c333    format('nclose: ',i5,a8,1x,a8,1x,a8,1x,a1,1x,f6.2,1x,f5.1,1x,
915c    1          f8.2,1x,f6.2,1x,f8.4,1x,f8.4)
916c
917c   All done.  Push data back out to station parameter storage.
918c
9195     call drdwr(2,j,1,13)
9203     continue
921c
922c   Clean up event level data.
923c
924c   Compute hypocentral parameter standard deviations.  If depth is
925c   held, all parameters are held, or the solution is bad, zero
926c   appropriate standard deviations.
927c
928      sedep=0.
929      if(ihold.eq.4.or.ng.eq.'s') se=0.
930      if(ihold.ne.4.and.ng.ne.'s'.and.ng.ne.'b') go to 10
931      selat=0.
932      selon=0.
933      setim=0.
934      call errelp(navsum,0,se)
935      return
936c   Convert the units of the correlation matrix to km**2.
937 10   rs=radius-odep
938      rl=rs*olats
939      m1=3-ihold
940      do 11 k=1,m1
941      a(1,k)=rs*a(1,k)
942      a(k,1)=rs*a(k,1)
943      a(2,k)=rl*a(2,k)
944 11   a(k,2)=rl*a(k,2)
945c   Otherwise, compute standard deviations from the diagonal elements of
946c   the correlation matrix.
947      setim=rsc1*se
948      selat=rsc1*dsqrt(dmax1(a(1,1),0d0))
949      selon=rsc1*dsqrt(dmax1(a(2,2),0d0))
950      if(ihold.eq.0) sedep=rsc1*dsqrt(dmax1(a(3,3),0d0))
951c     print *,'setim selat selon sedep',setim,selat,selon,sedep
952      call errelp(navsum,m,se)
953      return
954c   Flag a problem with the pP-P depth.
955 13   ierr=-1
956      return
957      end
958
959      function normeq(navsum,ihold)
960c
961c $$$$$ calls cholin, cholsl, and drdwr $$$$$
962c
963c   Normeq constructs and solves the normal equations to determine the
964c   next hypocentral step.  Nphases are processed.  All usable
965c   (not residualed or noassed) data are used.  M is the number of
966c   equations solved.  It will be 4 for a hypocentral step and 3 for an
967c   epicentral step.  If all goes well normeq returns 1.  Normeq returns
968c   -1 if the matrix is singular.  If ifl is 1, the system of equations
969c   is decomposed and  solved.  If ifl is 2, the correlation matrix is
970c   computed assuming that the matrix has already been decomposed.  If
971c   ifl is 3, the system of equations is decomposed, but not solved.
972c   Rewritten on 21 July 1983 by R. Buland from a fragment of a routine
973c   by E. R. Engdahl.
974c
975      double precision a,b,c,x,s,avdr,eps
976c     double fac
977c   Physical constants.
978      common /circle/ pi,rad,deg,radius
979c   Storage for the normal equations.
980      common /matrix/ a(3,4),b(3),c(4),x(3),s(3),avdr(3)
981c   Hypocentral parameters.
982      include 'hypocm.inc'
983c   Station-reading parameters.
984      include 'pickcm.inc'
985      data eps,fac0/1d-5,31.06341/
986c
987c   Assume success.
988      normeq=1
989      jfl=1
990      m1=3-ihold
991c   Zero out array storage.
992      do 1 k=1,m1
993      do 1 l=k,m1
994 1    a(l,k)=0d0
995c   Loop over the stations accumulating the normal matrix and data
996c   vector.
997      do 2 i=1,navsum
998c   Pop one of the sorted residuals.
999      ierr=lpop(i,j,res)
1000c   Access station parameter storage.
1001      call drdwr(1,j,1,5)
1002      if(pflg.ne.'u') go to 2
1003c   C is the contribution to the normal equations.
1004      c(1)=deg*cos(azim)*dtdd-avdr(1)
1005      c(2)=-deg*sin(azim)*dtdd*olats-avdr(2)
1006      c(3)=dtdh-avdr(3)
1007c     print *,'normeq',j,avdr,c
1008c   Update the normal equations.
1009      do 3 k=1,m1
1010      do 3 l=k,m1
1011 3    a(l,k)=a(l,k)+c(k)*c(l)
1012 2    continue
1013c     print 500,((a(l,k),k=1,m1),l=1,m1)
1014c500  format(/(1x,1p3e15.4))
1015c
1016c   Solve the system.
1017c
1018c   Weight the system of equations by diagonal matrix S to improve
1019c   conditioning.  Note that here as above, only the lower triangle of
1020c   the normal matrix is needed by the Cholesky routine used below.
1021      do 5 k=1,m1
1022 5    s(k)=1d0/dmax1(dsqrt(a(k,k)),eps)
1023      do 6 k=1,m1
1024      do 6 l=k,m1
1025 6    a(l,k)=s(l)*s(k)*a(l,k)
1026c   Compute the Cholesky decomposition of the normal matrix.  The lower
1027c   triangle of array A is unchanged.  The upper triangle of A is
1028c   replaced by the Cholesky factor (offset by one column so everything
1029c   will fit).
1030 7    call cholin(3,m1,a,1,ierr)
1031c   Branch if there is no error.
1032      if(ierr.eq.0) go to 8
1033c   Bail out if the fixup has already failed.
1034      if(jfl.le.0) go to 13
1035      jfl=-1
1036c   Add eps times the identity matrix to the normal matrix to improve
1037c   conditioning and try again.
1038      do 9 k=1,m1
1039 9    a(k,k)=a(k,k)+eps
1040      go to 7
1041c
1042c   Transform the Cholesky factor in the upper triangle of A into the
1043c   inverse of the normal matrix (the covariance matrix) in place.
1044 8    call cholin(3,m1,a,-3,ierr)
1045c   Remove the weighting from the covariance matrix, shift it over to
1046c   replace the original a matrix, and reconstruct the lower triangle.
1047      do 11 k=1,m1
1048      do 11 l=k,m1
1049      a(k,l)=s(k)*s(l)*a(k,l+1)
1050 11   a(l,k)=a(k,l)
1051      return
1052c   Flag a singular solution.
1053 13   normeq=-1
1054      return
1055      end
1056
1057      function nstep(nsta,ihold,loop,iter,navsum,av,chisq,deld)
1058c
1059c $$$$$ calls normeq and stats $$$$$
1060c
1061c   Assuming that travel-time residuals and derivitives with respect to
1062c   distance and depth have already been computed, nstep finds the next
1063c   Gauss-Newton step, updates the hypocenter, and assesses the goodness
1064c   of fit of the new solution.  In the assessment phase, travel-time
1065c   residuals and derivitives for the new solution are computed.  Nsta
1066c   is the number of associated stations, ihold is the number of
1067c   constraint equations, loop is a character flag indicating the
1068c   iteration stage, iter is the current Gauss-Newton iteration number,
1069c   navsum is the number of usable arrivals, av is the average travel-
1070c   time residual, chisq is the sum of the squares of the residuals, and
1071c   deld is the length of the spatial step made in kilometers.  Nsta,
1072c   loop and iter0 are not changed.  Navsum is assumed to be valid on
1073c   entry, but is recomputed internally.  Av, chisq, and deld are
1074c   computed each time by nstep.  If nstep returns -1, depth was
1075c   constrained by nstep and ihold has been modified accordingly.
1076c   Otherwise, nstep returns +1.  Rewritten on 25 July 1983 by R. Buland
1077c   from a fragment of a routine by E. R. Engdahl.
1078c
1079      character*1 loop
1080      double precision a,b,c,x,s,avdr
1081      dimension hyprm(8),hypsav(4)
1082c   Logical units on input and output files
1083        common /unit/ min, mout, merr
1084c   Physical constants.
1085      common /circle/ pi,rad,deg,radius
1086c   Storage for the normal equations.
1087      common /matrix/ a(3,4),b(3),c(4),x(3),s(3),avdr(3)
1088c   Hypocentral parameters.
1089      include 'hypocm.inc'
1090c   Commands.
1091      include 'cmndcm.inc'
1092c   Back up hypocenter.
1093      common /hypfix/ hypbck(8)
1094c   Iteration control parameters.
1095      common /itcntl/ miter1,miter2,miter3,tol1,tol2,tol3,hmax1,hmax2,
1096     1 hmax3,tol,hmax,zmin,zmax,zshal,zdeep,rchi,roff,utol
1097      equivalence (osec,hyprm)
1098      data almost,eps,savz,addz,nexp/1.1,1e-4,-1.,0.,1/
1099c
1100c   Assume success.
1101      nstep=1
1102c   Find number of equations.
1103      m=4-ihold
1104      m1=m-1
1105c   If there aren't enough equations, bail out.
1106      if(m1.gt.navsum) return
1107c   Initialize output parameters.
1108      delh=0.
1109      delz=0.
1110      if(iter.le.1) cn=50.
1111c   Addz accounts for depth changes imposed outside nstep in the tag
1112c   line printed below.
1113      if(savz.ge.0.) addz=odep-savz
1114c   The damping factor, roff, is fiddled each time nstep is called.
1115c   This "randomization" avoids excessive iteration in certain
1116c   pathlogical cases.
1117      if(rchi.gt.roff) go to 33
1118      rchi=rchi+.0390625
1119      go to 34
1120 33   rchi=rchi-.21875
1121c
1122c   Use steepest descents to find the minimum of the piecewise linear
1123c   rank-sum penelty function.
1124c
1125c   Set the base hypocenter.
1126 34   do 3 i=1,m
1127 3    hypsav(i)=hyprm(i)
1128      cn=amax1(cn,2.*tol)
1129      x(1)=b(1)/(rad*radius)
1130      x(2)=b(2)/(rad*radius*olats)
1131      x(3)=b(3)
1132      call lintry(navsum,ihold,chisq,cn,tol,hmax)
1133c     print *,'x',(x(i),i=1,m1)
1134c   Set the horizontal and vertical steps (in km).
1135      delh=sqrt(amax1((sngl(x(1))*rad*radius)**2+
1136     1 (sngl(x(2))*rad*radius*olats)**2,0.))
1137      if(ihold.eq.0) delz=x(3)
1138      savdel=sqrt(amax1(delh*delh+delz*delz,0.))
1139c     print *,'del',delh,delz,savdel
1140c
1141c   Damp the step on maximum horizontal distance.
1142c
1143      if(delh.le.hmax) go to 8
1144      rat=hmax/delh
1145      do 9 i=1,m1
1146 9    x(i)=rat*x(i)
1147      delh=hmax
1148c
1149c   Update the depth.
1150c
1151 8    if(ihold.ne.0) go to 4
1152c
1153c   Trap airquakes.
1154      if(odep+x(3).ge.zmin) go to 10
1155      if(abs(odep-zmin).le.eps) go to 27
1156c   If we aren't yet at zmin, damp the step to just reach zmin.
1157      rat=(zmin-odep)/x(3)
1158      do 18 i=1,m1
1159 18   x(i)=rat*x(i)
1160      delh=rat*delh
1161c   If the damped step would be interpreted as convergence, constrain
1162c   depth.
1163      if(sqrt(amax1(delh*delh+sngl(x(3))*sngl(x(3)),0.))-tol) 20,20,11
1164c   If we are already at zmin, move along the z = zmin surface.
1165 27   if(delh.le.tol) go to 20
1166      x(3)=0d0
1167      go to 11
1168c
1169c   Trap lower mantle quakes.
1170 10   if(odep+x(3).le.zmax) go to 11
1171      if(abs(zmax-odep).le.eps) go to 28
1172c   If we aren't yet at zmax, damp step to just get there.
1173      rat=(zmax-odep)/x(3)
1174      do 19 i=1,m1
1175 19   x(i)=rat*x(i)
1176      delh=rat*delh
1177c   If the damped step would be interpreted as convergence, constrain
1178c   depth.
1179      if(sqrt(amax1(delh*delh+sngl(x(3))*sngl(x(3)),0.))-tol) 21,21,11
1180c   If we are already at zmax, move along the z = zmax surface.
1181 28   if(delh.le.tol) go to 21
1182      x(3)=0d0
1183 11   delz=x(3)
1184      go to 4
1185c
1186c   On a bad depth, reset it to be held at either zshal or zdeep as
1187c   apropriate and set flags.
1188c
1189 20   rat=zshal
1190      go to 22
1191 21   rat=zdeep
1192 22   if(rat.ne.hypbck(4)) go to 29
1193c   If we have already converged at the held depth, reset the hypocenter
1194c   to the known solution.
1195      delh=sqrt(amax1(((hypbck(2)-hypsav(2))*rad*radius)**2+((hypbck(3)-
1196     1 hypsav(3))*rad*radius*olats)**2,0.))
1197      delz=hypbck(4)-hypsav(4)
1198      do 30 i=1,8
1199 30   hyprm(i)=hypbck(i)
1200      nstep=0
1201      go to 31
1202c   Otherwise, use the current epicenter, reset depth to the desired
1203c   value, and find an approximately compatible origin time.
1204 29   do 24 i=1,m
1205 24   hyprm(i)=hypsav(i)
1206      delh=0.
1207      delz=rat-odep
1208      osec=osec+.1*delz
1209      odep=rat
1210      nstep=-1
1211c   Set the depth held flags.
1212 31   ihold=1
1213      m=3
1214      m1=2
1215      hmax=hmax1
1216      dep ='c'
1217c   Go recompute residuals and derivitives for the next iteration.
1218      go to 1
1219c
1220c   Update the hypocentral parameters.
1221c
1222 4    do 12 i=1,m1
1223 12   hyprm(i+1)=hypsav(i+1)+x(i)
1224      odep=amin1(amax1(odep,zmin),zmax)
1225c   Guarantee legal co-latitude and longitude.
1226      if(olat.ge.0..and.olat.le.180.) go to 15
1227      olon=olon+180.
1228      if(olat.lt.0.) olat=abs(olat)
1229      if(olat.gt.180.) olat=360.-olat
1230 15   if(olon.lt.0.) olon=olon+360.
1231 17   if(olon.le.360.) go to 16
1232      olon=olon-360.
1233      go to 17
1234c   Calculate epicentral sines and cosines.
1235 16   olats=sin(olat*rad)
1236      olatc=cos(olat*rad)
1237      olons=sin(olon*rad)
1238      olonc=cos(olon*rad)
1239c
1240c   Recompute travel-time residuals, derivitives, and chi**2.
1241c
1242 1    call stats(nsta,ihold,navsum,av,chi2)
1243      hyprm(1)=hyprm(1)+av
1244c   Deld is the convergence parameter.
1245      deld=sqrt(delh*delh+delz*delz)
1246c
1247c   If chi**2 has decreased or damping has been unsuccessful, get out.
1248c
1249      if(chi2.lt.chisq.or.nstep.le.0.or.ng.ne.' ') go to 6
1250c   If we can't damp the step any more, bail out.
1251      if(rchi*deld.le.tol) go to 25
1252c   If chi**2 has increased, damp the step length by factor rchi.
1253      do 7 i=1,m1
1254 7    x(i)=rchi*x(i)
1255      delh=rchi*delh
1256      delz=rchi*delz
1257      go to 4
1258c   Trap an unstable solution on a clipped depth.
1259 25   if(abs(odep-zmin).le..01.and.ihold.eq.0) go to 20
1260      if(abs(odep-zmax).le..01.and.ihold.eq.0) go to 21
1261c   Flag an unstable solution.  Go around one more time to recover the
1262c   starting point which can't be improved.
1263      ng='u'
1264      if(savdel.le.utol) ng='q'
1265      if(chi2.le.almost*chisq.and.deld.le.tol) ng='o'
1266c   Zero the step and go around again to recompute the original
1267c   residuals.
1268      do 26 i=1,m1
1269 26   x(i)=0d0
1270      delh=0.
1271      delz=0.
1272      go to 4
1273c
1274c   Compute the event standard deviation.
1275 6    chisq=chi2
1276      rms=0.0
1277      if(navsum.gt.0) rms=chisq/navsum
1278c
1279c   Write a summary line to the standard output.
1280      rlat=90.-olat
1281      rlon=olon
1282      if(rlon.gt.180.) rlon=rlon-360.
1283      delz=delz+addz
1284      savz=odep
1285      write(merr,100,iostat=ios)loop,iter,navsum,rlat,rlon,odep,
1286     1 delh,delz,rms,ng,dep
1287 100  format('nstep: ',a1,i3,i5,f9.3,f10.3,f7.1,' delh=',f6.1,
1288     1       ' delz=',f7.1,' rms=',f10.4,1x,'ng= ',a1,'dep= ',a1)
1289 13   return
1290      end
1291
1292      subroutine scopt(n,sc)
1293c
1294c   Generate scores for n data.  Note that the scores don't depend on 
1295c   the data value, only on it's position.
1296c
1297      dimension sc(1),p(29),s(29)
1298      data p/0.,.1375,.1625,.1875,.2125,.2375,.2625,.2875,.3125,.3375,
1299     1 .3625,.3875,.4125,.4375,.4625,.4875,.5125,.5375,.5625,.5875,
1300     2 .6125,.6375,.6625,.6875,.7125,.7375,.7625,.7875,1./
1301      data s/0.0775,0.0775,0.1546,0.5328,0.8679,1.1714,1.4542,1.7266,
1302     1 1.9987,2.2802,2.5803,2.9068,3.2657,3.6603,4.0912,4.5554,5.0470,
1303     2 5.5572,6.0754,6.5906,7.0919,7.5702,8.0194,8.4365,8.8223,9.1812,
1304     3 9.5207,9.5974,9.5974/
1305c
1306      cn=1./(n+1)
1307      k1=1
1308      k=2
1309      av=0.
1310      do 1 i=1,n
1311      x=cn*i
1312 3    if(x.le.p(k)) go to 2
1313      k1=k
1314      k=k+1
1315      go to 3
1316 2    sc(i)=(x-p(k1))*(s(k)-s(k1))/(p(k)-p(k1))+s(k1)
1317 1    av=av+sc(i)
1318      av=av/n
1319      do 4 i=1,n
1320 4    sc(i)=sc(i)-av
1321      return
1322      end
1323
1324      subroutine setup(nphase,ihold,navsum,av,chisq)
1325c
1326c $$$$$ calls ellip, delaz, drdwr $$$$$
1327c
1328c   Setup initializes the location process.  This involves calculating
1329c   distance and azimuth, travel-time and derivitives, setting the phase
1330c   flag, and implementing the pkp commands.  nphase is the number of
1331c   associated phases, ihold is the
1332c   number of constraints on the solution, navsum is the number of
1333c   usable first arrivals, av is the average travel-time residual, and
1334c   chisq is the sum of the squared residuals.  Navsum, av,and chisq are
1335c   computed by setup.  Rewritten on 22 July 1983 by R. Buland from a
1336c   fragment of a routine by E. R. Engdahl.
1337c
1338      include 'locprm.inc'
1339      character*8 ophase, ph,phcd1(maxtt)
1340c       integer stdin, stdout
1341      real*4 tt1(maxtt),dtdd1(maxtt),dtdh1(maxtt),dddp1(maxtt),usrc(2),
1342     1 ellip, elcor
1343        common /unit/ min, mout, merr
1344c   Physical constants.
1345      common /circle/ pi,rad,deg,radius
1346c   Hypocentral parameters.
1347      include 'hypocm.inc'
1348c   Event level analyst commands.
1349      include 'cmndcm.inc'
1350c   Station-reading parameters.
1351      include 'pickcm.inc'
1352c   Scores.
1353      common /sccom/ scgen(maxph)
1354        data dtds/19.96/
1355c
1356c       stdout = 6
1357c       stdin = 5
1358c
1359c   Initialize counter
1360c
1361      navsum=0
1362c
1363c   Initialize depth for the calculation of the traveltimes and partial
1364c      derivatives
1365c
1366      if (odep .le. 0.0) odep = 1.0
1367        zs = odep
1368      call depset(zs,usrc)
1369c
1370c The first time that ellip is called, it initializes the routine and then
1371c   tcor is calculated on subsequent calls.  See further down in the
1372c   subroutine.
1373c
1374        tcor=ellip( 'P       ',100.0,zs,glat,45.0)
1375c
1376c   Loop over all station-readings.
1377c
1378        jj = 0
1379      do 2 j=1,nphase
1380c
1381c   Fetch data from station parameter storage.
1382c
1383      call drdwr(1,j,6,11)
1384c
1385c   Skip stations which can't be used in the location.
1386c
1387        if( pflg .ne. 'u' ) go to 2
1388c
1389c   Compute distance, azimuth and traveltime of phase.
1390c
1391      call delaz(olats,olatc,olons,olonc,slats,slatc,slons,slonc,
1392        1                       delta,azim)
1393c
1394c   Apply the DRES commands.
1395      do k=1,5
1396        if(dres(k).eq.'r'.and.delta.ge.dellim(1,k).and.delta.le.
1397     1   dellim(2,k)) then
1398          pflg=' '
1399          go to 8
1400        endif
1401      enddo
1402c
14036     call trtm(delta,maxtt,n,tt1,dtdd1,dtdh1,dddp1,phcd1)
1404c
1405c This part of the code is stupid.  Since some of the phases read in might not
1406c  be identified by branch (e.g PKPdf), we find the first phase that is closest
1407c  to read in phase.  For example, if the phase read in is "P", but the first
1408c  phase in time is "Pn" than the phase is rename to "Pn".
1409c
1410        do 111 i = 1, n
1411                ph = phcd1(i)
1412                ilen = ntrbl(phase)
1413                if(phase(1:ilen) .eq. ph(1:ilen)) then
1414                        travt = tt1(i)
1415                        dtdd = dtdd1(i)
1416                        dtdh = dtdh1(i)
1417                        ophase = phase
1418                        phase = ph
1419                        go to 100
1420                endif
1421111     continue
1422c
1423c Phase not found for that distance and depth
1424c
1425        pflg = ' '
1426c       write(merr,334)sta,phase
1427c334    format('setup: ','Sta= ',a5,' no phase computed for ',a8)
1428        go to 8
1429c
1430100     continue
1431c
1432c Calculate the ellipticity correction for different phases.
1433c
1434        azimuth = azim*deg
1435        tcor=ellip(phase,delta,zs,glat,azimuth)
1436      if(dtdd.gt.dtds) dtdd=dtds
1437      elcor=(elev/5.57)*sqrt (1.0-(5.57*dtdd/(radius*rad))**2)
1438      res=trvtim-osec-travt-(tcor+elcor)
1439c       res0 = res
1440c
1441c   Apply the RRES command.
1442      if(rres.eq.'r'.and.(res.lt.reslim(1).or.res.gt.reslim(2))) then
1443        pflg=' '
1444        go to 8
1445      endif
1446c
1447      ierr=lpush(navsum,j,res)
1448c       write(merr,333)j,navsum,sta,ophase,ph,pflg,delta,azimuth,
1449c       2                               trvtim,res,dtdd,dtdh
1450c333    format('setup: ',i5,i5,1x,a8,1x,a8,1x,a8,1x,a1,1x,f6.2,1x,f5.1,1x,
1451c    1          f8.2,1x,f6.2,1x,f8.4,1x,f8.4)
1452c
1453c   Restore data to station parameter storage.
1454c
1455 8    call drdwr(2,j,1,5)
1456c       call drdwr(2,j,22,22)
1457 2    continue
1458c
1459      call scopt(navsum,scgen)
1460c
1461c   Compute all statistical parameters.
1462c
1463      ierr=lsrt(navsum,av,se,chisq,'TT')
1464      call adder(navsum,ihold)
1465c       write(merr,*)'setup: ',' nphase= ',nphase,' navsum= ',navsum
1466      return
1467      end
1468
1469      subroutine stats(nphase,ihold,navsum,av,chisq)
1470c
1471c $$$$$ calls ellip, delaz, drdwr $$$$$
1472c
1473c   Stats computes travel-time residuals and derivitives for the current
1474c   solution.  Nphase is the number of associated stations, ihold is the
1475c   number of constraint equations, navsum is the number of usable
1476c   arrivals, av is the average travel-time residual, and chisq is the
1477c   sum of the squared residuals.  Nphase and ihold are not changed.
1478c   Navsum, av, and chisq are computed by stats.  Rewritten on
1479c   5 August 1983 by R. Buland from a fragment of a routine by
1480c   E. R. Engdahl.
1481c
1482      include 'locprm.inc'
1483      character*8 phcd1(maxtt),ph, ophase
1484c       integer stdout,stdin
1485      real*4 tt1(maxtt),dtdd1(maxtt),dtdh1(maxtt),dddp1(maxtt),usrc(2),
1486     1 ellip, elcor
1487c   Physical constants.
1488      common /circle/ pi,rad,deg,radius
1489c   Logical units.
1490      common /unit/ min,mout,merr
1491c   Hypocentral parameters.
1492      include 'hypocm.inc'
1493c   Station-reading parameters.
1494      include 'pickcm.inc'
1495        data dtds/19.96/
1496c
1497c       stdout = 6
1498c       stdin = 5
1499c
1500c   Initialize statistical counters.
1501c
1502      navsum=0
1503c
1504c Initialize the depth
1505c
1506        if (odep .le. 0.0) odep = 1.0
1507        zs = odep
1508      call depset(zs,usrc)
1509c   Loop over all station-readings.
1510      do 1 j=1,nphase
1511c   Fetch data from station parameter storage.
1512      call drdwr(1,j,6,11)
1513c   Don't bother with bad stations or noassed or residualed readings.
1514        if( pflg .ne. 'u' ) go to 1
1515c   Get the new distance and azimuth.
1516      call delaz(olats,olatc,olons,olonc,slats,slatc,slons,slonc,
1517        1                       delta,azim)
1518c
1519      call trtm(delta,maxtt,n,tt1,dtdd1,dtdh1,dddp1,phcd1)
1520c
1521c This part of the code is tricky.  Since some of the phases read in might not
1522c  be identified by branch (e.g PKPdf), we find the first phase that is closest
1523c  to read in phase.  For example, if the phase read in is "P", but the first
1524c  phase in time is "Pn" than the phase is rename to "Pn".
1525c
1526        do 111 i = 1, n
1527                ph = phcd1(i)
1528                ilen = ntrbl(phase)
1529                if(phase(1:ilen) .eq. ph(1:ilen)) then
1530                        travt = tt1(i)
1531                        dtdd = dtdd1(i)
1532                        dtdh = dtdh1(i)
1533                        ophase = phase
1534                        phase = ph
1535                        go to 100
1536                endif
1537111     continue
1538c
1539c Phase not found for that distance and depth
1540c
1541        pflg = ' '
1542        go to 8
1543c
1544100     continue
1545c
1546c
1547        azimuth = azim*deg
1548      if(dtdd.gt.dtds) dtdd=dtds
1549      elcor=(elev/5.57)*sqrt (1.0-(5.57*dtdd/(radius*rad))**2)
1550        tcor=ellip(phase,delta,zs,glat,azim)
1551      res=trvtim-osec-travt-(tcor+elcor)
1552c
1553c       write(merr,333)j,sta,ophase,ph,pflg,delta,azimuth,
1554c       2                               trvtim,res,dtdd,dtdh
1555c333    format('stats: ',i5,a8,1x,a8,1x,a8,1x,a1,1x,f6.2,1x,f5.1,1x,
1556c    1          f8.2,1x,f6.2,1x,f8.4,1x,f8.4)
1557c
1558c   Update the statistical sums.
1559c
1560      ierr=lpush(navsum,j,res)
1561c
1562c   Save computed data in station parameter storage.
1563c
15648     call drdwr(2,j,1,5)
15651       continue
1566c
1567c   Compute the average travel-time residual and chi**2.
1568c
1569      ierr=lsrt(navsum,av,se,chisq,'TT')
1570      call adder(navsum,ihold)
1571      return
1572      end
1573
1574      function stirf(u,t0,t1,t2)
1575c
1576c   Compute an asymptotic approximation to the Stirling function.
1577c
1578      stirf=t1+0.5*u*(t2-t0)+0.5*(u**2)*(t2-2.0*t1+t0)
1579      return
1580      end
1581 
1582      function xmed(n,a,sd)
1583c
1584c $$$$$ calls only library routines $$$$$
1585c
1586c   Given the n numbers, a(1) <= a(2) <= ... <= a(n), xmed returns
1587c   robust estimates of center and spread of the data about the center.
1588c   The array a is not altered.  Programmed on 9 December 1983 by
1589c   R. Buland.
1590c
1591      dimension a(1)
1592      data xno,xne/1.482580,.7412898/
1593c   See if there is enough data to work with.
1594      if(n-1)13,12,14
1595c   No, return.
1596 13   xmed=0.
1597      sd=0.
1598      return
1599c   One datum is a special case.
1600 12   xmed=a(1)
1601      sd=0.
1602      return
1603c   N > 1, the median is a robust estimate of center.
1604 14   m=n/2
1605      md=mod(n,2)
1606      if(md.le.0) xmed=.5*(a(m)+a(m+1))
1607      if(md.gt.0) xmed=a(m+1)
1608c   Find the smallest positive value.
1609      do 1 i=1,n
1610      if(a(i)-xmed.ge.0.) go to 2
1611 1    continue
1612      i=n
1613c   Find the smallest absolute value.
1614 2    if(i.gt.1.and.abs(a(i-1)-xmed).lt.abs(a(i)-xmed)) i=i-1
1615      k=i
1616      j=1
1617      if(j.ge.m) go to 15
1618c   Implicitly order the series in order of increasing absolute value by
1619c   sucessively moving i up or k down until the m th term is found.
1620 6    if(k.le.1) go to 3
1621      if(i.ge.n) go to 4
1622      if(abs(a(i+1)-xmed).gt.abs(a(k-1)-xmed)) go to 5
1623c   Move i up.
1624      i=i+1
1625      j=j+1
1626      if(j.lt.m) go to 6
1627c   The m th term is found, mark it.
1628 15   l1=i
1629      go to 8
1630c   Move k down.
1631 5    k=k-1
1632      j=j+1
1633      if(j.lt.m) go to 6
1634c   The m th term is found, mark it.
1635      l1=k
1636c   Find and mark the m+1 th term.
1637 8    if(k.le.1) go to 9
1638      if(i.ge.n) go to 10
1639      if(abs(a(i+1)-xmed)-abs(a(k-1)-xmed))9,9,10
1640 9    l2=i+1
1641      go to 11
1642 10   l2=k-1
1643      go to 11
1644c   We have run out of negative terms, the m th term must be up.
1645 3    l1=i+m-j
1646      l2=l1+1
1647      go to 11
1648c   We have run out of positive terms, the m th term must be down.
1649 4    l1=k-(m-j)
1650      l2=l1-1
1651c   The spread is the normalized median of the absolute series.
1652 11   if(md.le.0) sd=xne*(abs(a(l1)-xmed)+abs(a(l2)-xmed))
1653      if(md.gt.0) sd=xno*abs(a(l2)-xmed)
1654      return
1655      end
Note: See TracBrowser for help on using the repository browser.