subroutine homest (storxy,nfixes,iolp,t,ind,davex,davey,dvarx
     1 ,dvary,dcov,darea)
      integer nfixes, iolp
      real storxy(3,nfixes)
      double precision t
      double precision davex, davey, dvarx, dvary, dcov
      real darea
c
c **********************************************************************
c     code name = homest    written by james e. dunn, 1972-1977.
c                          adapted for spss by gary c. white, may, 1983.
c                t    sampling interval within bursts (decimal).
c                nburst   number of bursts of data.
c                na   number of animals whose coordinates were observed.
c                nc   number of coordinates/animal (2 is assumed).
c                ind  nature of output required:
c                        0 for minimal output.
c                        1 for observed, fitted, and residual values.
c                        2 for likelihood functions.
c                        3 for both residuals and likelihood functions.
c                        4 to force 10 iterative steps (diagnostic only)
c
c     remarks:
c     (1) separate bursts are assumed to be statistically independent.
c     (2) an equal sampling interval is assumed within bursts.
c     (3) the present program is limited to a maximum of two coordinates
c         expansion may be accomplished by increasing the values in the
c         parameter statements.
c     (4) up to 300 bursts are allowed with present dimensions
c **********************************************************************
      integer ish1, idate(5)
      logical filopn
      integer maxcor, maxcr2, maxcr4, maxcr5, maxcr6, maxcr7, maxcr8
c     parameter (maxcor=2, maxcr2=maxcor*maxcor, maxcr4=maxcr2*maxcr2,
c    1 maxcr5=maxcr2*(maxcr2+1)/2, maxcr6=maxcor*2, maxcr7=maxcr6+1,
c    2 maxcr8=maxcr6*maxcr7+1)
      parameter (maxcor=2, maxcr2=4, maxcr4=28,
     1 maxcr5=10, maxcr6=4, maxcr7=5,
     2 maxcr8=21)
      parameter (maxbur=300)
      double precision a, o, s, sl, r, z, g, phi, b, c, d, tr, sini, e,
     1 u, sinn, crifn, crifn1, f, bb
      dimension z(20), g(maxcr2), phi(maxcr2), b(maxcr4), c(maxcr4), d
     1 (maxcr4), tr(maxcr4), sini(maxcr4), e(maxcr4), u(20), sinn(maxcr4
     2 ), f(maxcr2), bb(maxcr5), o(maxcr4), a(maxcr4), s(maxcr4), sl
     3 (maxcr4), r(maxcr4), n(maxbur)
      double precision ack(maxcr7,maxcr6), tck(maxcr6,maxcr6), xck
     1 (maxcr6,maxcr7)
c     equivalence (c(1),ack(1)),(c(maxcr8),xck(1)),(tck(1),f(1))
      real crvals
      double precision crval
      equivalence (crvals,crval)
      character errstr*8
      integer ier
c
      ish1=14
   10 ish1=ish1+1
      if (ish1.lt.99) go to 15
 1150 write (unit=iolp,fmt='(11x,a/11x,a,l4,i4,i4)')
     1 'error from '//errstr//'with opening a scratch file.',
     1 'values for filopn, ier and ish1 are :',filopn,ier,ish1
      call linect(2)
      return
 15   errstr='inquire'
      inquire (unit=ish1,opened=filopn,iostat=ier,err=1150)
      if (filopn) then
         go to 10
      else
         errstr='open'
         open (unit=ish1,status='scratch',form='unformatted',recl=1024
     1   ,iostat=ier,err=1150)
      endif
c
      if (ind.ge.0) then
            call dphead
            write (unit=iolp,fmt='(//10x,9hdata list)')
      call linect (3)
      endif
      ncoor=2
      na=1
      nc=na*ncoor
      if (nc.gt.maxcor) call error (902)
      call mdchi(0.95,real(nc),crvals,ier)
      crval=crvals
      nc2=nc*nc
      do 20 i=1,nc
      d(i)=0
      s(i)=0
   20 sl(i)=0
      do 30 i=1,nc2
      sini(i)=0
      a(i)=0
   30 r(i)=0
      rewind (unit=ish1)
      iobs=1
      nburst=0
   40 nburst=nburst+1
      if (ind.ge.0) then
         write (unit=iolp,fmt='(//)')
         call linect (2)
      endif
      if (nburst.gt.maxbur) then
         write (unit=iolp,fmt='(/1x,a,i4)')
     1   'number of bursts in data exceeds maximum of ',maxbur
         call linect(2)
         return
      endif
      n(nburst)=0
      call readxy (storxy(1,iobs),o,nc)
      if (ind.ge.0) then
         call homer9 (idate,storxy(3,iobs))
         write (unit=iolp,fmt=1160) idate,(o(i),i=1,nc)
         call linect (1+nc/6)
      endif
      do 50 i=1,nc
   50 d(i)=d(i)+o(i)
      k=0
      do 60 j=1,nc
      do 60 i=1,nc
      k=k+1
   60 sini(k)=sini(k)+o(i)*o(j)
c     read and print o'th obs.
   70 iobs=iobs+1
      if (iobs.gt.nfixes) go to 100
      call readxy (storxy(1,iobs),o(nc+1),nc)
      if (storxy(3,iobs)-storxy(3,iobs-1).gt.t) go to 40
      if (ind.ge.0) then
         call homer9 (idate,storxy(3,iobs))
         write (unit=iolp,fmt=1160) idate,(o(nc+i),i=1,nc)
         call linect (1+nc/6)
      endif
      n(nburst)=n(nburst)+1
      do 80 i=1,nc
      sl(i)=sl(i)+o(i)
      s(i)=s(i)+o(nc+i)
      k=(i-1)*nc
      do 80 j=1,nc
      a(k+j)=a(k+j)+o(i)*o(j)
   80 r(k+j)=r(k+j)+o(j)*o(nc+i)
      do 90 i=1,nc
   90 o(i)=o(nc+i)
      go to 70
  100 k=0
      do 120 j=1,nc
      do 120 i=1,nc
      k=k+1
      if (nburst.ne.1) go to 110
      sini(k)=0.
      go to 120
  110 sini(k)=(sini(k)-d(i)*d(j)/nburst)/(nburst-1)
  120 continue
      nobs=0
      k=0
      do 140 i=1,nburst
      if (n(i).eq.0) go to 130
      nobs=nobs+n(i)
      go to 140
  130 k=k+1
  140 continue
      if (ind.ge.0) then
         write (unit=iolp,fmt=1170) nburst,nobs,t
         call linect (8)
      endif
      fn=float(nobs)
      do 150 i=1,nc
      k=(i-1)*nc
      do 150 j=1,nc
      a(k+j)=a(k+j)-sl(i)*sl(j)/fn
  150 r(k+j)=r(k+j)-sl(j)*s(i)/fn
      ier=0
      eps=1.0e-10
      call dgelg (r,a,nc,nc,eps,ier)
c     gamma transpose resides columnwise in r.
c     columns of gamma transpose are rows of gamma.
      if (ind.ge.0) then
         write (unit=iolp,fmt=1180)
         call linect (1)
      endif
      write (unit=ish1) (r(i),i=1,nc2)
      do 170 i=1,nc
      l=(i-1)*nc
      do 160 j=1,nc
  160 ack(i,j)=r(l+j)
      if (ind.ge.0) then
         write (unit=iolp,fmt=1190) (r(l+j),j=1,nc)
         call linect (nc/6+1)
      endif
  170 continue
      if (na.gt.6) go to 190
c     calculate tridiagonal matrix tck from gamma.
      call lanczo (ack,tck,nc,z,o,xck,iolp,ind)
c     calculate the coefficients of a polynomial of order nc.
      call coeff (tck,nc,a,z,o)
c     calculate the eigenvalues.
      call dpolrt (a,b,nc,z,o,ier)
      if (ind.ge.0) then
         write (unit=iolp,fmt=1200)
         call linect (2)
      endif
      rabs=0.
      do 180 i=1,nc
         abei=dsqrt(z(i)**2+o(i)**2)
         rabs=amax1(rabs,abei)
         reb=alog(abei)/t
         cob=datan(o(i)/z(i))/t
         if (ind.ge.0) then
            write (unit=iolp,fmt=1210) z(i),o(i),abei,reb,cob
            call linect (2)
         endif
  180    continue
      if (rabs.lt.1.0) go to 190
      if (ind.ge.0) then
         write (unit=iolp,fmt=1220)
         call linect (2)
      endif
      go to 1140
  190 continue
c     now form vector 'nu'.
      do 200 i=1,nc
         d(i)=d(i)/float(nburst)
  200    z(i)=s(i)
      do 210 i=1,nc
         k=(i-1)*nc
         do 210 j=1,nc
  210       z(i)=z(i)-r(k+j)*sl(j)
      do 220 i=1,nc
  220    z(i)=z(i)/fn
c     'nu' is resident in z.
      if (ind.ge.0) then
         write (unit=iolp,fmt=1230) (z(i),i=1,nc)
         call linect (6+nc/7)
      endif
c     now form vector 'mu'.
      do 230 i=1,nc
         o(i)=z(i)
         k=(i-1)*nc
         do 230 j=1,nc
            l=(j-1)*nc+i
  230       a(k+j)=-r(l)
      do 240 i=1,nc2
  240    g(i)=-a(i)
      do 250 i=1,nc
         k=(i-1)*nc+i
  250    a(k)=a(k)+1.
      write (unit=ish1) (a(i),i=1,nc2)
      call dgelg (o,a,nc,1,eps,ier)
      if (ind.ge.0) then
         write (unit=iolp,fmt=1240) (o(i),i=1,nc)
         call linect (6+nc/7)
      endif
      do 260 i=1,nc
  260 u(i)=o(i)
c     now form matrix 'phi'.
      do 270 i=1,nc2
  270 a(i)=0
      iobs=0
      do 340 i=1,nburst
      iobs=iobs+1
      call readxy (storxy(1,iobs),o,nc)
      if (n(i).eq.0) go to 340
      l=n(i)
      do 330 j=1,l
      do 280 k=1,nc
  280 o(nc+k)=0
      do 290 ii=1,nc
      k=(ii-1)*nc
      do 290 jj=1,nc
  290 o(nc+ii)=o(nc+ii)+r(k+jj)*o(jj)
      iobs=iobs+1
      call readxy (storxy(1,iobs),o,nc)
      do 300 ii=1,nc
  300 o(nc+ii)=o(nc+ii)+z(ii)
      do 310 ii=1,nc
  310 o(nc+ii)=-o(nc+ii)+o(ii)
      do 320 ii=1,nc
      k=(ii-1)*nc
      do 320 jj=1,nc
  320 a(k+jj)=a(k+jj)+o(nc+ii)*o(nc+jj)
  330 continue
  340 continue
      write (unit=ish1) (d(i),i=1,nc)
      do 350 i=1,nc2
  350 a(i)=a(i)/fn
      if (ind.ge.0) then
         write (unit=iolp,fmt=1250)
         call linect (5)
      endif
      do 360 i=1,nc
      k=(i-1)*nc
      if (ind.ge.0) then
         write (unit=iolp,fmt=1190) (a(k+j),j=1,nc)
         call linect (2+nc/6)
      endif
  360 continue
      if (ncoor.ne.2) go to 380
      do 370 ll=1,na
      i1=2*(ll-1)*nc+(ll-1)*2+1
      i2=i1+1
      i4=((ll-1)*2+1)*nc+2*ll
      area=dsqrt(a(i1)*a(i4)-a(i2)**2)*18.82272672
      if (ind.ge.0) then
         write (unit=iolp,fmt=1260) ll,area
         call linect (2)
      endif
  370 continue
  380 continue
      do 390 i=1,nc2
  390 e(i)=a(i)
      l=0
      do 400 i=1,nc
      k=(i-1)*nc
      do 400 j=1,i
      l=l+1
  400 phi(l)=a(k+j)
c     now form matrix 'lambda'.
      call krprod (nc,nc,g,nc,nc,g,n1,n2,r)
      nn=n1*n2
      do 410 i=1,nn
  410 r(i)=-r(i)
      do 420 i=1,n1
      k=(i-1)*n1+i
  420 r(k)=r(k)+1.
      call dgelg (a,r,n1,1,eps,ier)
      if (ind.ge.0) then
         write (unit=iolp,fmt=1270)
         call linect (5)
      endif
      do 430 i=1,nc
      k=(i-1)*nc
      if (ind.ge.0) then
         write (unit=iolp,fmt=1190) (a(k+j),j=1,nc)
         call linect (2+nc/6)
      endif
  430 continue
      l=0
      do 440 i=1,nc
      k=(i-1)*nc
      do 440 j=1,i
      l=l+1
  440 r(l)=a(k+j)
      call deigen (r,o,nc,1)
      i2=nc*(nc+1)/2
      noitco=0
      if (r(i2).gt.0.0) go to 470
      if (ind.ge.0) then
         write (unit=iolp,fmt=1280)
         call linect (5)
      endif
      noitco=1
      call dgmprd (g,sini,o,nc,nc,nc)
      call dgmtra (g,r,nc,nc)
      call dgmprd (o,r,a,nc,nc,nc)
      do 450 i=1,nc2
  450 a(i)=a(i)+e(i)
      do 460 i=1,nc
      k=(i-1)*nc
      if (ind.ge.0) then
         write (unit=iolp,fmt=1190) (a(k+j),j=1,nc)
         call linect (2+nc/6)
      endif
  460 continue
  470 continue
      if (ncoor.ne.2) go to 490
      do 480 ll=1,na
      i1=2*(ll-1)*nc+(ll-1)*2+1
      i2=i1+1
      i4=((ll-1)*2+1)*nc+2*ll
      area=dsqrt(a(i1)*a(i4)-a(i2)**2)*18.82272672
      if (ind.ge.0) then
         write (unit=iolp,fmt=1290) ll,area
         call linect (3)
      endif
  480 continue
  490 continue
      write (unit=ish1) (a(i),i=1,nc2)
      write (unit=ish1) (e(i),i=1,nc2)
      call dsinv (phi,nc,eps,ier)
      if (ind.ge.0) then
         write (unit=iolp,fmt=1300) (d(i),i=1,nc)
         call linect (6+nc/7)
      endif
      if (ind.ge.0) then
         write (unit=iolp,fmt=1310)
         call linect (5)
      endif
      do 500 i=1,nc
      k=(i-1)*nc
      if (ind.ge.0) then
         write (unit=iolp,fmt=1190) (sini(k+j),j=1,nc)
         call linect (6+nc/7)
      endif
  500 continue
      l=0
      do 510 i=1,nc
      k=(i-1)*nc
      do 510 j=1,i
      l=l+1
  510 r(l)=a(k+j)
      call dsinv (r,nc,eps,ier)
      l=0
      do 520 i=1,nc
      do 520 j=1,nc
      l=l+1
      call loc (i,j,ir,nc,nc,1)
      f(l)=a(l)
  520 a(l)=r(ir)
      l=0
      do 530 i=1,nc
      do 530 j=1,nc
      l=l+1
      call loc (i,j,ir,nc,nc,1)
  530 r(l)=phi(ir)
      write (unit=ish1) (r(i),i=1,nc2)
      write (unit=ish1) (a(i),i=1,nc2)
      write (unit=ish1) (u(i),i=1,nc)
      rewind (unit=ish1)
      read (unit=ish1) (g(i),i=1,nc2)
c     gamma transpose is in g.
      read (unit=ish1) (o(i),i=1,nc2)
c     i - gamma is in o.
      read (unit=ish1) (d(i),i=1,nc)
      read (unit=ish1) (s(i),i=1,nc2)
c     lambda is in s.
      read (unit=ish1) (phi(i),i=1,nc2)
c     phi is in phi.
      read (unit=ish1) (r(i),i=1,nc2)
c     phi inverse is in r.
c     estimated mu from initial obs is in d.
c     lambda inverse is in a.
      rewind (unit=ish1)
c     transformation matrix is stored columnwise in tr.
      ntrow=nc*(nc+1)/2
      l=ntrow*nc2
      do 540 i=1,l
  540 tr(i)=0
      lll=1
      do 580 j=1,nc
      do 580 i=1,nc
      l=(lll-1)*nc*(nc+1)/2
      do 560 jj=1,nc
      do 560 ii=1,jj
      ll=l+jj*(jj-1)/2+ii
      if ((jj.eq.j.and.ii.eq.i).or.(jj.eq.i.and.ii.eq.j)) go to 550
      go to 560
  550 tr(ll)=1.
      go to 570
  560 continue
  570 lll=lll+1
  580 continue
      inw=ind
c     if (inw.eq.0.or.inw.eq.4) inw=2
      if (inw.eq.4) inw=2
      call likftn (nc,nc2,fn,nburst,n,g,o,s,phi,r,a,b,c,d,e,u,sl,sinn,
     1    crval,ntrow,tr,crifn,inw,storxy,nfixes,iolp)
      nbrach=1
      if (noitco.eq.0.and.ind.ge.0) then
         call dphead
         write (unit=iolp,fmt=1320)
         call linect (6)
      endif
  590 continue
      call krprod (nc,nc,g,nc,nc,g,n1,n2,e)
      nn=n1*n2
      do 600 i=1,nn
  600 e(i)=-e(i)
      do 610 i=1,nc2
      l=(i-1)*nc2+i
  610 e(l)=e(l)+1.
      call krprod (nc,nc,r,nc,nc,r,n1,n2,c)
      call dgmprd (e,c,d,nc2,nc2,nc2)
      call dgmtra (e,c,nc2,nc2)
      call dgmprd (d,c,b,nc2,nc2,nc2)
      nn=nc2*nc2
      do 620 i=1,nn
  620 b(i)=.5*b(i)
      if (nbrach.ne.10) go to 640
      do 630 i=1,nn
  630 e(i)=d(i)
  640 continue
      call dgmprd (tr,b,c,ntrow,nc2,nc2)
      call dgmtra (tr,b,ntrow,nc2)
      call dgmprd (c,b,o,ntrow,nc2,ntrow)
      call krprod (nc,nc,a,nc,nc,a,n1,n2,b)
      nn=n1*n2
      do 650 i=1,nn
  650 b(i)=.5*b(i)
      call dgmprd (tr,b,c,ntrow,nc2,nc2)
      call dgmtra (tr,b,ntrow,nc2)
      call dgmprd (c,b,sl,ntrow,nc2,ntrow)
      if (ind.ge.0) then
         write (unit=iolp,fmt=1330)
         call linect (3)
      endif
      do 660 i=1,ntrow
      l=(i-1)*ntrow
      if (ind.ge.0) then
         write (unit=iolp,fmt=1190) (sl(l+j),j=1,ntrow)
         call linect (2+nc/6)
      endif
  660 continue
      if (ind.ge.0) then
         write (unit=iolp,fmt=1340)
         call linect (3)
      endif
      do 670 i=1,ntrow
      l=(i-1)*ntrow
      if (ind.ge.0) then
         write (unit=iolp,fmt=1190) (o(l+j),j=1,ntrow)
         call linect (2+nc/6)
      endif
  670 continue
      if (noitco.eq.1) go to 780
      l=0
      do 680 i=1,nc
      k=(i-1)*nc
      do 680 j=1,i
      l=l+1
      b(l)=sini(k+j)
  680 c(l)=f(k+j)
      call dgmprd (sl,b,d,ntrow,ntrow,1)
      call dgmprd (o,c,b,ntrow,ntrow,1)
      do 690 i=1,ntrow
  690 b(i)=float(nburst)*d(i)+fn*b(i)
      l=0
      do 700 i=1,ntrow
      k=(i-1)*ntrow
      do 700 j=1,i
      l=l+1
  700 r(l)=float(nburst)*sl(k+j)+fn*o(k+j)
      call dgels (b,r,ntrow,1,eps,ier,s)
c     weighted estimate of lambda is in b.
      nn=nc*(nc+1)/2
      do 710 i=1,nn
  710 r(i)=b(i)
      call dsinv (r,nc,eps,ier)
c     inverse of weighted estimate of lambda is in r.
      l=0
      do 720 i=1,nc
      do 720 j=1,nc
      l=l+1
      call loc (i,j,ir,nc,nc,1)
      s(l)=b(ir)
  720 a(l)=r(ir)
c     weighted lambda in s, inverse in a.
      if (ind.ge.0) then
         write (unit=iolp,fmt=1350) nbrach
         call linect (3)
      endif
      do 730 i=1,nc
      l=(i-1)*nc
      if (ind.ge.0) then
         write (unit=iolp,fmt=1190) (s(l+j),j=1,nc)
         call linect (2+nc/6)
      endif
  730 continue
      call dgmprd (s,g,b,nc,nc,nc)
      call dgmtra (g,c,nc,nc)
      call dgmprd (c,b,d,nc,nc,nc)
      do 740 i=1,nc2
  740 phi(i)=s(i)-d(i)
      if (ind.ge.0) then
         write (unit=iolp,fmt=1360)
         call linect (3)
      endif
      do 750 i=1,nc
      l=(i-1)*nc
      if (ind.ge.0) then
         write (unit=iolp,fmt=1190) (phi(l+j),j=1,nc)
         call linect (2+nc/6)
      endif
  750 continue
      l=0
      do 760 i=1,nc
      k=(i-1)*nc
      do 760 j=1,i
      l=l+1
  760 b(l)=phi(k+j)
      call dsinv (b,nc,eps,ier)
      l=0
      do 770 i=1,nc
      do 770 j=1,nc
      l=l+1
      call loc (i,j,ir,nc,nc,1)
  770 r(l)=b(ir)
c     phi inverse is in r.
  780 read (unit=ish1) (g(i),i=1,nc2)
      read (unit=ish1) (o(i),i=1,nc2)
      read (unit=ish1) (d(i),i=1,nc)
      call dgmprd (r,o,b,nc,nc,nc)
      call dgmtra (b,sl,nc,nc)
      call dgmtra (o,c,nc,nc)
      call dgmprd (c,b,o,nc,nc,nc)
      call dgmprd (sl,z,b,nc,nc,1)
      call dgmprd (s,b,c,nc,nc,1)
      do 790 i=1,nc
  790 sinn(i)=float(nburst)*d(i)+fn*c(i)
      call dgmprd (s,o,d,nc,nc,nc)
      do 800 i=1,nc2
  800 d(i)=fn*d(i)
      do 810 i=1,nc
      l=(i-1)*nc+i
  810 d(l)=d(l)+float(nburst)
      call dgelg (sinn,d,nc,1,eps,ier)
      do 820 i=1,nc
  820 u(i)=sinn(i)
      if (ind.ge.0) then
         write (unit=iolp,fmt=1370) (u(i),i=1,nc)
         call linect (5+nc/7)
      endif
      inz=0
      if (ind.ge.0) inz=2
      call likftn (nc,nc2,fn,nburst,n,g,o,s,phi,r,a,b,c,d,e,u,sl,sinn,
     1 crval,ntrow,tr,crifn1,inz,storxy,nfixes,iolp)
      if (crifn1.ge.crifn.and.noitco.eq.1) ncode=1
      if (crifn1.lt.crifn.and.ind.ne.4) go to 830
      crifn=crifn1
      if (nbrach.eq.10) go to 830
      write (unit=ish1) (s(i),i=1,nc2)
c     new lambda is in s. old lambda from succ. obs. is in f.
      write (unit=ish1) (phi(i),i=1,nc2)
c     new phi is in phi.
      write (unit=ish1) (r(i),i=1,nc2)
c     new phi inverse in r.
      write (unit=ish1) (a(i),i=1,nc2)
c     new lambda inverse is in a.
      write (unit=ish1) (u(i),i=1,nc)
      rewind (unit=ish1)
      if (noitco.ne.1) nbrach=nbrach+1
      if (noitco.ne.1) go to 590
  830 rewind (unit=ish1)
      read (unit=ish1) (g(i),i=1,nc2)
c     gamma transpose is in g.
      read (unit=ish1) (o(i),i=1,nc2)
c     i - gamma is in o.
      read (unit=ish1) (d(i),i=1,nc)
c     estimated mu from initial obs is in d.
      read (unit=ish1) (s(i),i=1,nc2)
c     lambda is in s.
      read (unit=ish1) (phi(i),i=1,nc2)
c     phi is in phi.
      read (unit=ish1) (r(i),i=1,nc2)
c     phi inverse is in r.
      read (unit=ish1) (a(i),i=1,nc2)
c     lambda inverse is in a.
      read (unit=ish1) (z(i),i=1,nc)
      if (nbrach.ne.1) go to 840
      if (ind.ge.0) then
         write (unit=iolp,fmt=1380)
         call linect (3)
      endif
      if (ncode.eq.1.and.ind.ge.0) then
         write (unit=iolp,fmt=1390)
         call linect (1)
      endif
  840 if (ind.ge.0) then
         write (unit=iolp,fmt=1400) (z(i),i=1,nc)
         call linect (6+nc/6)
      endif
      davex=z(1)
      davey=z(2)
      if (ncoor.ne.2) go to 870
      do 850 ll=1,na
      i1=2*(ll-1)*nc+(ll-1)*2+1
      i2=i1+1
      i4=((ll-1)*2+1)*nc+2*ll
      area=dsqrt(phi(i1)*phi(i4)-phi(i2)**2)*18.82272672
      if (ind.ge.0) then
         write (unit=iolp,fmt=1260) ll,area
         call linect (2)
      endif
  850 continue
      do 860 ll=1,na
      i1=2*(ll-1)*nc+(ll-1)*2+1
      i2=i1+1
      i4=((ll-1)*2+1)*nc+2*ll
      area=dsqrt(s(i1)*s(i4)-s(i2)**2)*18.82272672
      if (ind.ge.0) then
         write (unit=iolp,fmt=1290) ll,area
         call linect (2)
      endif
      dvarx=s(i1)
      dvary=s(i4)
      dcov=s(i2)
c   convert area from square meters to hectares
      darea=area*0.0001
  860 continue
  870 continue
      call dgmprd (r,o,b,nc,nc,nc)
      call dgmtra (o,c,nc,nc)
      call dgmprd (c,b,o,nc,nc,nc)
      rewind (unit=ish1)
      if (ind.ge.0) then
         call dphead
         write (unit=iolp,fmt=1410)
         call linect (7)
         write (unit=iolp,fmt=1420)
         call linect (3)
      endif
      do 880 i=1,nc
      l=(i-1)*nc
      if (ind.ge.0) then
         write (unit=iolp,fmt=1190) (a(l+j),j=1,nc)
         call linect (2+nc/6)
      endif
  880 continue
      if (ind.ge.0) then
         write (unit=iolp,fmt=1430)
         call linect (3)
      endif
      do 890 i=1,nc
      l=(i-1)*nc
      if (ind.ge.0) then
         write (unit=iolp,fmt=1190) (o(l+j),j=1,nc)
         call linect (2+nc/6)
      endif
  890 continue
      l=0
      do 900 i=1,nc
      k=(i-1)*nc
      do 900 j=1,i
      l=l+1
  900 b(l)=fn*o(k+j)+float(nburst)*a(k+j)
      if (ind.ge.0) then
         write (unit=iolp,fmt=1440)
         call linect (3)
      endif
      do 910 i=1,nc
      l=i*(i-1)/2
      if (ind.ge.0) then
         write (unit=iolp,fmt=1190) (b(l+j),j=1,i)
         call linect (2+nc/6)
      endif
  910 continue
      call dsinv (b,nc,eps,ier)
      l=0
      do 920 i=1,nc
      do 920 j=1,nc
      l=l+1
      call loc (i,j,ir,nc,nc,1)
  920 o(l)=b(ir)
      if (ind.ge.0) then
         write (unit=iolp,fmt=1450)
         call linect (3)
      endif
      do 930 i=1,nc
      l=(i-1)*nc
      if (ind.ge.0) then
         write (unit=iolp,fmt=1190) (o(l+j),j=1,nc)
         call linect (2+nc/6)
      endif
  930 continue
      call dgmprd (s,g,b,nc,nc,nc)
      call dgmprd (b,r,c,nc,nc,nc)
      call dgmtra (c,sl,nc,nc)
      call krprod (nc,nc,sl,nc,nc,c,n1,n2,d)
      call dgmprd (b,sl,c,nc,nc,nc)
      call krprod (nc,nc,c,nc,nc,r,n1,n2,sl)
      call krprod (nc,nc,s,nc,nc,r,n1,n2,c)
      l=0
      do 940 i=1,nc2
      k=(i-1)*nc2
      do 940 ii=1,nc
      do 940 jj=1,nc
      l=l+1
      ll=ii+(jj-1)*nc+k
  940 b(l)=d(ll)
      nn=n1*n2
      do 950 i=1,nn
  950 d(i)=c(i)+sl(i)+b(i)
c     f(2,2) is in d
      if (ind.ge.0) then
         write (unit=iolp,fmt=1460)
         call linect (3)
      endif
      do 960 i=1,nc2
      l=(i-1)*nc2
      if (ind.ge.0) then
         write (unit=iolp,fmt=1190) (d(l+j),j=1,nc2)
         call linect (2+nc/6)
      endif
  960 continue
      ll=0
      do 970 i=1,nc2
      l=(i-1)*nc2
      do 970 j=1,i
      ll=ll+1
  970 bb(ll)=fn*d(l+j)
      call dgmprd (s,g,b,nc,nc,nc)
      call dgmtra (b,c,nc,nc)
      do 980 i=1,nc2
  980 sl(i)=0
      do 990 i=1,nc
      l=(i-1)*nc+i
  990 sl(l)=1.
c     regenerate e.
      call krprod (nc,nc,g,nc,nc,g,n1,n2,b)
      nn=n1*n2
      do 1000 i=1,nn
 1000 b(i)=-b(i)
      do 1010 i=1,nc2
      l=(i-1)*nc2+i
 1010 b(l)=b(l)+1.
      call krprod (nc,nc,r,nc,nc,r,n1,n2,sinn)
      call dgmprd (b,sinn,e,nc2,nc2,nc2)
      call krprod (nc,nc,r,nc,nc,r,n1,n2,sinn)
      call dgmprd (b,sinn,sini,nc2,nc2,nc2)
      call dgmtra (b,sinn,nc2,nc2)
      call dgmprd (sini,sinn,b,nc2,nc2,nc2)
      nn=nc2*nc2
      do 1020 i=1,nn
 1020 b(i)=.5*b(i)
      call dgmprd (tr,b,sinn,ntrow,nc2,nc2)
      call dgmtra (tr,b,ntrow,nc2)
      call dgmprd (sinn,b,sini,ntrow,nc2,ntrow)
      call krprod (nc,nc,a,nc,nc,a,n1,n2,b)
      nn=n1*n2
      do 1030 i=1,nn
 1030 b(i)=.5*b(i)
      call dgmprd (tr,b,sinn,ntrow,nc2,nc2)
      call dgmtra (tr,b,ntrow,nc2)
      call dgmprd (sinn,b,d,ntrow,nc2,ntrow)
      nn=ntrow**2
      do 1040 i=1,nn
 1040 sini(i)=float(nburst)*d(i)+fn*sini(i)
      call krprod (nc,nc,c,nc,nc,sl,n1,n2,b)
      call dgmprd (e,b,c,nc2,nc2,nc2)
      l=0
      do 1050 i=1,nc2
      k=(i-1)*nc2
      do 1050 ii=1,nc
         do 1050 jj=1,nc
            l=l+1
            ll=ii+(jj-1)*nc+k
 1050       b(l)=c(ll)
      nn=nc2*nc2
      do 1060 i=1,nn
 1060    d(i)=(-b(i)-c(i))/2.
      call dgmprd (tr,d,c,ntrow,nc2,nc2)
      if (ind.ge.0) then
         write (unit=iolp,fmt=1470)
         call linect (3)
      endif
      do 1070 i=1,nc2
         l=(i-1)*ntrow
         if (ind.ge.0) then
            write (unit=iolp,fmt=1190) (c(l+j),j=1,ntrow)
            call linect (2+nc/6)
         endif
 1070    continue
c     will print out as f(2)3) even though it resides as f(3,2)
      nn=nc2*(nc2+1)/2
      do 1080 i=1,nn
 1080    b(i)=bb(i)
      call dgmtra (c,d,ntrow,nc2)

      do 1110 i=1,ntrow
         l=(i-1)*nc2
         do 1090 j=1,nc2
            nn=nn+1
            if (nn.gt.maxcr4) then
               write (unit=iolp,fmt='(1x,a)')
     1            ' Exceed dimensions on array b in HOMEST.FOR'
               stop 'Error in code, call Gary White'
            endif
 1090       b(nn)=d(l+j)*fn
         l=(i-1)*ntrow
         do 1100 j=1,i
            nn=nn+1
            if (nn.gt.maxcr4) then
               write (unit=iolp,fmt='(1x,a)')
     1            ' Exceed dimensions on array b in HOMEST.FOR'
               stop 'Error in code, call Gary White'
            endif
 1100       b(nn)=sini(l+j)
 1110 continue
      n1=ntrow+nc2
      k=0
      if (ind.ge.0) then
         write (unit=iolp,fmt=1480)
         call linect (3)
      endif
      do 1120 i=1,n1
         k=i*(i-1)/2
         if (ind.ge.0) then
            write (unit=iolp,fmt=1190) (b(k+j),j=1,i)
            call linect (2+nc/6)
         endif
 1120    continue
      call dsinv (b,n1,eps,ier)
      if (ind.ge.0) then
         write (unit=iolp,fmt=1490)
         call linect (3)
      endif
      do 1130 i=1,n1
         k=i*(i-1)/2
         if (ind.ge.0) then
            write (unit=iolp,fmt=1190) (b(k+j),j=1,i)
            call linect (2+nc/6)
         endif
 1130    continue
      if (nbrach.eq.1) go to 1140
         read (unit=ish1) (g(i),i=1,nc2)
         read (unit=ish1) (o(i),i=1,nc2)
         read (unit=ish1) (d(i),i=1,nc)
         read (unit=ish1) (s(i),i=1,nc2)
         read (unit=ish1) (phi(i),i=1,nc2)
         read (unit=ish1) (r(i),i=1,nc2)
         read (unit=ish1) (a(i),i=1,nc2)
         call likftn (nc,nc2,fn,nburst,n,g,o,s,phi,r,a,b,c,d,e,z,
     1      sl,sinn,crval,ntrow,tr,crifn,inw,storxy,nfixes,iolp)
 1140 close (unit=ish1)
      return
c
 1160 format (1x,2(i2.2,1h/),i2.2,1x,2i2.2,6g14.8/(1x,6g14.8))
 1170 format (/11x,'number of data bursts =',i5//11x,'number of successi
     1ve observations =',i5//11x,'sampling interval =',f10.2//)
 1180 format (1h ,'estimated gamma matrix')
 1190 format (1h0,6d16.8/(1h ,6d16.8))
 1200 format (/4x,'eigenvalues of gamma',15x,'absolute value',10x,'eigen
     1values of b')
 1210 format (1h0,g12.5,' + ',g12.5,' i',10x,f10.5,7x,g12.5,' + ',g12.5,
     1 ' i')
 1220 format (1h0,'the process does not appear to be stationary. executi
     1on is terminated.')
 1230 format (//1h0,'estimated vector nu'/1h0,7d16.8/(1h ,7d16.8))
 1240 format (//1h0,'estimated vector mu from successive observations'/1
     1 h0,7d16.8/(1h ,7d16.8))
 1250 format (//1h0,'estimate of phi matrix'/)
 1260 format (/1h0,'area of 0.95 conditional confidence ellipse for anim
     1al',1x,i3,' = ',g12.5)
 1270 format (//1h0,'estimated matrix lambda from successive observation
     1s'/)
 1280 format (/1h0,'the above matrix is not positive definite. it is rep
     1laced by'/1h ,'the following estimate: lambda(s) = phi + g*lambda(
     2i)*g.'/1h ,'note: no iterated solution will be obtained.')
 1290 format (/1h0,'area of 0.95 confidence ellipse for animal',1x,i3,'
     1= ',g12.5)
 1300 format (//1h0,'estimated vector mu from initial observations'/1h0,
     1 7d16.8/(1h ,7d16.8))
 1310 format (//1h0,'estimated matrix lambda from initial observations'/
     1 )
 1320 format (1h0,'******************************************'/1h0,'iter
     1ated solution'/1h0,'******************************************')
 1330 format (/1h0,'information on lambda/initial observation')
 1340 format (/1h0,'information on lambda/successive observation')
 1350 format (/1h0,'iteration ',i4,'  estimated matrix lambda')
 1360 format (/1h0,'estimated phi matrix')
 1370 format (/1h0,'estimated vector mu'/1h0,7d16.8/(1h ,7d16.8))
 1380 format (1h0,'no iterative improvement is possible.  suggest using'
     1 /1h ,'the previously listed estimates from successive observation
     2s only')
 1390 format (1h ,'and the following:')
 1400 format (1h0,'******************************************'/1h0,'comb
     1ined estimate of mu'/1h0,6d16.8/(1h ,6d16.8))
 1410 format (/1h0,'******************************************'/1h0,'est
     1imated information and covariances'/1h0,'*************************
     2*****************')
 1420 format (/1h0,'information on mu/initial observation')
 1430 format (/1h0,'information on mu/successive observation')
 1440 format (/1h0,'total information on mu')
 1450 format (/1h0,'covariance for mu')
 1460 format (/1h0,'information on gamma/successive observation')
 1470 format (/1h0,'joint information on gamma and lambda/successive obs
     1ervation')
 1480 format (/1h0,'total information on gamma and lambda')
 1490 format (/1h0,'covariance of gamma and lambda')
      end