cdeck homer5
      subroutine homer5 (nfixes,ndim,storxy)
      integer nfixes,ndim
      real storxy(3,ndim)
c
c     homer5 computes home range estimates by the minimum area method.
c     note - x and y arrays are rearranged by subroutine homer5 - do not
c     attempt to use again in calling program
      include 'maincomo.for'
      include 'homerc.for'
      real ext(4)
      double precision sumx, sumy, sumsqx, sumsqy, sumsqxy
      double precision davex, davey, dvarx, dvary, dcov
      real darea
      double precision avex, avey, varx, vary, cov, uavex, uavey, uvarx,
     1 uvary, ucov, timval, temp, e1, e2
      real areajt, uareajt, area, arear, acres, acresr, xdmax, ydmax, df
     1 , chi025, chi975, ajtl, ajtu, uajtl, uajtu
      integer iobs(1400),i1,i2,i3,nrep
c   calculate value of pi
      pi=2.*asin(1.)
c
c   perform test for serial correlation of locations
c
      if (qoptsr(12)) then
      i1=nfixes+1
      i2=i1+nfixes/3+1
      nrep=1000
      i3=i2+nrep/3+1
      if (i3.lt.ndim) then
         call homr5a(storxy,nfixes,storxy(1,i1),storxy(1,i2),nrep)
      else
         write(unit=iolp,fmt='(/9x,a)')
     1      'too little memory for randomization test'
         call linect(2)
      endif
      endif
c
c nuf  - number of unique fixes
c
      nuf=0
      if (nfixes.lt.3) go to 90
      if (qoptsr(6)) then
         davex=0.
         davey=0.
         dvarx=0.
         dvary=0.
         dcov=0.
         darea=0.
         if (duntim(idic).eq.0.) then
             timval=(storxy(3,2)-storxy(3,1))*1.000001
         else
            timval=duntim(idic)
         endif
         call homest (storxy,nfixes,iolp,timval,dunprt,davex,davey,dvarx
     1   ,dvary,dcov,darea)
      endif
      if (qoptsr(3).and.qoptsr(4)) return
      sumx=dble(storxy(1,1))
      sumy=dble(storxy(2,1))
      sumsqx=sumx*sumx
      sumsqy=sumy*sumy
      sumsqxy=sumx*sumy
c     ext - (1) = maxy, (2) = minx, (3) = miny, (4) = maxx
      ext(1)=storxy(2,1)
      ext(3)=storxy(2,1)
      ext(2)=storxy(1,1)
      ext(4)=storxy(1,1)
      if (qoptsr(13)) then
c
c        Write raw locations to unit 21
c
         open (unit=21,file='homer.loc', status='unknown')
      endif
      do 10 i=2,nfixes
         if (qoptsr(13)) then
            write (21,'(2(1x,g14.8))') storxy(1,i), storxy(2,i)
         endif
c     store minimum and maximum x and y
         if (storxy(2,i).gt.ext(1)) then
            ext(1)=storxy(2,i)
         elseif (storxy(2,i).lt.ext(3)) then
            ext(3)=storxy(2,i)
         endif
         if (storxy(1,i).gt.ext(4)) then
            ext(4)=storxy(1,i)
         elseif (storxy(1,i).lt.ext(2)) then
            ext(2)=storxy(1,i)
         endif
c         cumulate sums and sums of squares and cross-products
         if (.not.qoptsr(4)) then
         sumx=sumx+dble(storxy(1,i))
         sumy=sumy+dble(storxy(2,i))
         sumsqx=sumsqx+dble(storxy(1,i))**2
         sumsqy=sumsqy+dble(storxy(2,i))**2
         sumsqxy=sumsqxy+dble(storxy(1,i))*dble(storxy(2,i))
         endif
   10    continue
      if (qoptsr(13)) then
         close(unit=21)
      endif
c     compute area of enclosing rectangle, maximum x and y distances,
c     center of activity and its variances.
      xdmax=ext(4)-ext(2)
      ydmax=ext(1)-ext(3)
      arear=xdmax*ydmax
      if (.not.qoptsr(4)) then
      avex=sumx/dble(nfixes)
      avey=sumy/dble(nfixes)
      varx=(sumsqx-sumx**2/dble(nfixes))/dble(nfixes-1)
      vary=(sumsqy-sumy**2/dble(nfixes))/dble(nfixes-1)
      cov=(sumsqxy-sumx*sumy/dble(nfixes))/dble(nfixes-1)
      areajt=dble(pi)*sqrt(varx*vary-cov*cov)*5.99d0
c   convert from square meters to hectares
      areajt=areajt*0.0001
c   generate confidence bounds
      df=2.*real(nfixes)-4.
      call mdchi (0.025,df,chi025,ier)
      ajtu=df*areajt/chi025
      call mdchi (0.975,df,chi975,ier)
      ajtl=df*areajt/chi975
c   cumulate statistics for unique fixes
      sumx=0.d0
      sumy=0.d0
      sumsqx=0.d0
      sumsqy=0.d0
      sumsqxy=0.d0
      nuf=nfixes
      do 50 i=1,nfixes
         if (i.gt.nuf) go to 60
         sumx=sumx+dble(storxy(1,i))
         sumy=sumy+dble(storxy(2,i))
         sumsqx=sumsqx+dble(storxy(1,i))**2
         sumsqy=sumsqy+dble(storxy(2,i))**2
         sumsqxy=sumsqxy+dble(storxy(1,i))*dble(storxy(2,i))
         if (i.eq.nuf) go to 60
         do 40 k=i+1,nfixes
   20       if (k.gt.nuf) go to 50
c           check for redundant fixes.
            if (storxy(1,i).ne.storxy(1,k).or.storxy(2,i).ne.storxy(2,k)
     1      ) go to 40
c      swap redundant fix and the last unique fix
            do 30 ii=1,2
               temp=storxy(ii,k)
               storxy(ii,k)=storxy(ii,nuf)
               storxy(ii,nuf)=temp
   30          continue
            nuf=nuf-1
            if (k.gt.nuf) go to 50
            go to 20
   40       continue
   50    continue
c     compute center of range and its variance - unique fixes only.
   60 if (nuf.le.3) then
         uavex=0.d0
         uavey=0.d0
         uvarx=0.d0
         uvary=0.d0
         ucov=0.d0
         uareajt=0.
         uajtl=0.
         uajtu=0.
         area=0.
         npt=0
            else
         uavex=sumx/dble(nuf)
         uavey=sumy/dble(nuf)
         uvarx=(sumsqx-sumx**2/dble(nuf))/dble(nuf-1)
         uvary=(sumsqy-sumy**2/dble(nuf))/dble(nuf-1)
         ucov=(sumsqxy-sumx*sumy/dble(nuf))/dble(nuf-1)
         uareajt=dble(pi)*sqrt(uvarx*uvary-ucov*ucov)*5.99d0
c      convert from square meters to hectares
         uareajt=uareajt*0.0001
c   generate confidence bounds
      df=2.*real(nuf)-4.
      call mdchi (0.025,df,chi025,ier)
      uajtu=df*uareajt/chi025
      call mdchi (0.975,df,chi975,ier)
      uajtl=df*uareajt/chi975
      endif
      endif
      call dphead
      if (.not.qoptsr(3)) then
c         compute area of minimum convex polygon
            if (nuf.gt.280) then
c              homr10 is a faster method of calculating the minimum
c              area polygon, but is limited by storage space.
c              the array iobs has only 1400 storage spaces.
c              iobs must be of length 5*nuf, hence only use homr10 if
c              nuf*5 is less than 1400.
               call homer6 (storxy,nuf,npt,area)
               do 70 i=1,npt
   70             iobs(i)=i
            else
              call homr10 (storxy,nuf,npt,area,iobs)
            endif
         acres=area*0.0002471054
         acresr=arear*0.0002471054
c        convert square meters to hectares
         area=area*0.0001
         arear=arear*0.0001
            write (iolp,100) area,acres
            call linect (5)
            if (area.gt.arear) then
               write (iolp,110)
               call linect (1)
            endif
         write (iolp,150) npt
         call linect (1)
c
c      print peripheral points of minimum area convex polygon
c
         write (iolp,200)
         write (iolp,210) (iobs(i),storxy(1,iobs(i)),storxy(2,iobs(i)),i
     1   =1,npt)
         call linect (npt)
         if (qoptsr(8)) then
c
c               Print perimeter of convex polygon to output
c                       unit 21, homer.con
c
            open (unit=21, file='homer.con', status='unknown')
            write (21,'(2(2x,g14.8))')
     1         (storxy(1,iobs(j)),storxy(2,iobs(j)),j=1,npt),
     2          storxy(1,iobs(1)),storxy(2,iobs(1))
            close(unit=21)
         endif
         write (iolp,120) arear,acresr
         call linect (4)
         write (iolp,130) ext(2),ext(3)
         call linect (2)
         write (iolp,140) xdmax,ydmax
         call linect (2)
      endif
      if (.not.qoptsr(4)) then
         write (unit=iolp,fmt=160)
         call linect (4)
         write (iolp,190) avex,uavex,avey,uavey,varx,uvarx,vary,uvary
     1   ,cov,ucov,areajt,uareajt,ajtl,uajtl,ajtu,uajtu,areajt*2.471054
     2   ,uareajt*2.471054,ajtl*2.471054,uajtl*2.471054,ajtu*2.471054
     3   ,uajtu*2.471054,nfixes,nuf
         call linect (13)
c
c   generate lotus file for 95% ellipses if options set
c
         if (qoptsr(9)) then
            temp=sqrt((varx+vary)**2-4.d0*(varx*vary-cov**2))
            e1=(varx+vary+temp)/2.d0
            e2=(varx+vary-temp)/2.d0
            theta=atan2(cov,e1-vary)*180./pi
            if (theta.lt.0.) theta=theta+360.
            rl=sqrt(5.99*e1)
            rs=sqrt(5.99*e2)
c                         Calculate 200 perimeter points for ellipse,
c                         send to output file
            rotat = theta * pi / 180.
            sinrot = sin( rotat )
            cosrot = cos( rotat )
            open (unit=21,file='homer.elp', status='unknown')
            do 1000 theta2 = 0., 6.2831853, 0.0314159
              xelip = rl * cos( theta2 )
              yelip = rs * sin( theta2 )
              xpt = xelip*cosrot - yelip*sinrot + avex
              ypt = xelip*sinrot + yelip*cosrot + avey
              write (21,'(1x,g14.8,3x,g14.8)') xpt,ypt
 1000         continue
            close(unit=21)
         endif
         if (nuf.le.3) go to 90
         if (qoptsr(10)) then
            temp=sqrt((uvarx+uvary)**2-4.d0*(uvarx*uvary-ucov**2))
            e1=(uvarx+uvary+temp)/2.d0
            e2=(uvarx+uvary-temp)/2.d0
            theta=atan2(ucov,e1-uvary)*180./pi
            if (theta.lt.0.) theta=theta+360.
            rl=sqrt(5.99*e1)
            rs=sqrt(5.99*e2)
c               Calculate 200 perimeter points for ellipse,
c                   send to output file
            rotat = theta * pi / 180.
            sinrot = sin( rotat )
            cosrot = cos( rotat )
            open (unit=21,file='homer.elu', status='unknown')
            do 1001 theta2 = 0., 6.2831853, 0.0314159
              xelip = rl * cos( theta2 )
              yelip = rs * sin( theta2 )
              xpt = xelip*cosrot - yelip*sinrot + avex
              ypt = xelip*sinrot + yelip*cosrot + avey
              write (21,'(1x,g14.8,3x,g14.8)') xpt,ypt
 1001         continue
            close(unit=21)
         endif
      endif
      if (qoptsr(6) .and. davex.ne.0. .and. davey.ne.0.) then
         write (unit=iolp,fmt=170)
         call linect (3)
         write (unit=iolp,fmt=180) davex,davey,dvarx,dvary,dcov,darea
     1   ,darea*2.471054
         call linect (6)
         if (davex.ne.0..and.davey.ne.0.) then
            if (qoptsr(11)) then
               temp=sqrt((dvarx+dvary)**2-4.d0*(dvarx*dvary-dcov**2))
               e1=(dvarx+dvary+temp)/2.d0
               e2=(dvarx+dvary-temp)/2.d0
               theta=atan2(dcov,e1-dvary)*180./pi
               if (theta.lt.0.) theta=theta+360.
               rl=sqrt(5.99*e1)
               rs=sqrt(5.99*e2)
c               Calculate 200 perimeter points for ellipse,
c                   send to output file
               rotat = theta * pi / 180.
               sinrot = sin( rotat )
               cosrot = cos( rotat )
               open (unit=21,file='homer.dun', status='unknown')
               do 1002 theta2 = 0., 6.2831853, 0.0314159
                  xelip = rl * cos( theta2 )
                  yelip = rs * sin( theta2 )
                  xpt = xelip*cosrot - yelip*sinrot + avex
                  ypt = xelip*sinrot + yelip*cosrot + avey
                  write (21,'(1x,g14.8,3x,g14.8)') xpt,ypt
 1002             continue
              close(unit=21)
            endif
         endif
      endif
      return
   90 write (iolp,220) nfixes,nuf
      call linect (3)
      return
c
  100 format (///9x,'home range by minimum area method',g14.8,
     1 'hectares'/41x,'(',g14.8,'acres)')
  110 format (9x,5(1h*),'disregard preceeding home range estimate becaus
     1e it exceeds the area of the rectangle',5(1h*))
  120 format (///9x,'area of home range rectangle',g14.8,'hectares'/
     1 37x,'(',g14.8,'acres)')
  130 format (9x,'coordinates of lower left-hand corner of rectangle   x
     1 = ',g14.8/62x,'y = ',g14.8)
  140 format (9x,'maximum x distance  ',g14.8,' meters'/
     1        9x,'maximum y distance  ',g14.8,' meters')
  150 format (9x,'number of sides on minimum area polygon',3x,i3)
  160 format (///9x,'jennrich-turner home range estimates')
  170 format (//9x,'dunn home range estimates')
  180 format (34x,'all fixes'/19x,'mean x',8x,g14.8/19x,'mean y',8x,g14.
     1 8/19x,'variance x',4x,g14.8/19x,'variance y',4x,g14.8/19x,'covari
     2ance',4x,g14.8/19x,'95% ellipse',3x,g14.8,'  (hectares)'/
     3 19x,'95% ellipse',3x,g14.8,'  (acres)')
  190 format (34x,'all fixes    unique fixes'/19x,'mean x',8x,g14.8,2x
     1 ,g14.8/19x,'mean y',8x,g14.8,2x,g14.8/19x,'variance x',4x,g14.8,2
     2 x,g14.8/19x,'variance y',4x,g14.8,2x,g14.8/19x,'covariance',4x
     3 ,g14.8,2x,g14.8/19x,'95% ellipse',3x,g14.8,2x,g14.8,'  (hectares)
     4 '/19x,'lower 95% ci',2x,g14.8,2x,g14.8,'  (hectares)'
     1   /19x,'upper 95% ci',2x,g14.8,2x,g14.8,'  (hectares)'
     2   /19x,'95% ellipse',3x,g14.
     6 8,2x,g14.8,'  (acres)'/19x,'lower 95% ci',2x,g14.8,2x,g14.8,'  (a
     7cres)'/19x,'upper 95% ci',2x,g14.8,2x,g14.8,'  (acres)'/19x,'no. f
     8ixes',5x,i5,11x,i5)
  200 format (//9x,'coordinates of convex polygon peripheral points'/18x
     1 ,'case number',9x,'x',14x,'y')
  210 format (22x,i3,8x,g14.8,1x,g14.8)
  220 format (/9x,'no minimum area or ellipse method calculations',4x,'n
     1o. of fixes = ',i5,'  no. of unique fixes = ',i5/)
      end