subroutine likftn (nc,nc2,fn,nburst,n,g,o,s,phi,r,a,b,c,d,e,z,sl,
     1 sini,crval,ntrow,tr,sum,ind,storxy,nfixes,iolp)
      double precision g(1), o(1), s(1), phi(1), r(1), a(1), b(1), c(1),
     1 d(1), e(1), z(1), sl(1), sini(1), sum, tot, tr(1)
      integer nfixes, iolp
      real storxy(3,nfixes)
      dimension n(1)
      character*4 yes, no, hit
      data yes /' yes'/, no /'  no'/
c     gamma transpose in g.
c     i - gamma in o.
c     lambda in s.
c     phi in phi.
c     phi inverse in r.
c     lambda inverse in a.
c     mu in z.
      sum=0
      miss=0
      eps=1.e-14
      do 10 i=1,nc2
      d(i)=0
      sini(i)=0
      b(i)=0
      c(i)=0
   10 e(i)=0
      if (ind.eq.1.or.ind.eq.3) then
            call dphead
            write (unit=iolp,fmt=210)
         call linect (2)
      endif
      iobs=0
      do 130 i=1,nburst
      iobs=iobs+1
      call readxy (storxy(1,iobs),o,nc)
      do 20 j=1,nc
      o(j)=o(j)-z(j)
   20 b(j)=b(j)+o(j)
      do 30 j=1,nc
      do 30 k=1,nc
      l=(k-1)*nc+j
   30 sum=sum+a(l)*o(j)*o(k)
      l=0
      do 40 j=1,nc
      do 40 jj=1,nc
      l=l+1
   40 sini(l)=sini(l)+o(j)*o(jj)
      if (ind.eq.1.or.ind.eq.3.and.n(i).ne.0) then
            write (unit=iolp,fmt=220) n(i)
            call linect (4)
      endif
      if (n(i).eq.0) go to 130
      ll=n(i)
      do 120 ij=1,ll
      do 50 k=1,nc
   50 o(nc+k)=0
      do 70 ii=1,nc
      k=(ii-1)*nc
      do 60 jj=1,nc
   60 o(nc+ii)=o(nc+ii)+g(k+jj)*o(jj)
      o(nc+ii)=o(nc+ii)+z(ii)
   70 sl(ii)=o(nc+ii)
c     predicted resides in rhs of o.
      iobs=iobs+1
      call readxy (storxy(1,iobs),o(2*nc+1),nc)
      do 80 ii=1,nc
      o(nc+ii)=o(2*nc+ii)-o(nc+ii)
   80 c(ii)=c(ii)+o(nc+ii)
      tot=0
      do 90 j=1,nc
      do 90 k=1,nc
      l=(k-1)*nc+j
   90 tot=tot+r(l)*o(nc+j)*o(nc+k)
c     write (unit=iolp,fmt=21)tot
      if (tot.gt.crval) miss=miss+1
      hit=no
      if (tot.le.crval) hit=yes
      sum=sum+tot
      l=0
      do 100 j=1,nc
      do 100 jj=1,nc
      l=l+1
      d(l)=d(l)+o(nc+j)*o(nc+jj)
  100 e(l)=e(l)+o(j)*o(nc+jj)
      if (ind.eq.1.or.ind.eq.3) then
         write (unit=iolp,fmt=230) i,ij,(o(2*nc+j),j=1,nc),(sl(j),j=1,nc
     1   ),(o(nc+j),j=1,nc)
         call linect (3)
      endif
      if (ind.eq.1.or.ind.eq.3) then
         write (unit=iolp,fmt=240) nc,tot,hit
         call linect (2)
      endif
c     b contains sum v(0)
c     c contains sum v(j)
c     sini contains v(0)*v(0)
c     d contains v(j)*v(j)
c     e contains v(j)*(z(j-1)-mu)
      do 110 j=1,nc
  110 o(j)=o(2*nc+j)-z(j)
  120 continue
  130 continue
      tot=miss/fn
      if (ind.eq.1.or.ind.eq.3) then
         write (unit=iolp,fmt=250) miss,fn,tot
         call linect (5)
      endif
      call dgmprd (a,b,sl,nc,nc,1)
      call dgmprd (r,c,b,nc,nc,1)
      do 140 j=1,nc
  140 sl(j)=sl(j)+b(j)
      call dgmprd (g,b,c,nc,nc,1)
      do 150 j=1,nc
  150 sl(j)=sl(j)-c(j)
      if (ind.eq.2.or.ind.eq.3) then
         write (unit=iolp,fmt=260) (sl(j),j=1,nc)
         call linect (4+nc/6)
      endif
      call dgmprd (s,g,c,nc,nc,nc)
      call dgmtra (c,b,nc,nc)
      call dgmprd (r,b,c,nc,nc,nc)
      call dgmprd (d,c,sl,nc,nc,nc)
      do 160 j=1,nc2
  160 sl(j)=fn*b(j)+e(j)-sl(j)
      if (ind.eq.2.or.ind.eq.3) then
         write (unit=iolp,fmt=270) (sl(j),j=1,nc2)
         call linect (4+nc/6)
      endif
      call dgmprd (g,r,c,nc,nc,nc)
      call dgmtra (c,b,nc,nc)
      call dgmprd (d,b,e,nc,nc,nc)
      call dgmprd (c,e,sl,nc,nc,nc)
      call dgmprd (g,b,e,nc,nc,nc)
      do 170 j=1,nc2
  170 sl(j)=fn*e(j)-sl(j)
      call dgmprd (r,d,c,nc,nc,nc)
      call dgmprd (c,r,b,nc,nc,nc)
      call dgmprd (a,sini,d,nc,nc,nc)
      call dgmprd (d,a,c,nc,nc,nc)
      do 180 j=1,nc2
  180 sl(j)=sl(j)-float(nburst)*a(j)+c(j)-fn*r(j)+b(j)
      call dgmprd (tr,sl,c,ntrow,nc2,1)
      if (ind.eq.2.or.ind.eq.3) then
         write (unit=iolp,fmt=280) (c(j),j=1,ntrow)
         call linect (4+nc/6)
      endif
      l=0
      do 190 i=1,nc
      k=(i-1)*nc
      do 190 j=1,i
      l=l+1
      b(l)=s(k+j)
  190 c(l)=phi(k+j)
      call dmfsd (b,nc,eps,ier)
      call dmfsd (c,nc,eps,ier)
      sum=-sum/2.
      l=0
      do 200 i=1,nc
      l=l+i
  200 sum=sum-dlog(b(l))-dlog(c(l))
      sum=sum-float(nc)*(fn+float(nburst))*0.918938534
      if (ind.eq.2.or.ind.eq.3) then
            write (unit=iolp,fmt=290) sum
            call linect (2)
      endif
      return
c
  210 format (1h0,'residual analysis'/)
  220 format (1h0,'number of observations =',i5/1h0,'burst','  obs.',5x,
     1 'observed/predicted/residual')
  230 format (1h0,2i5,10f8.2/1h ,10x,8f10.2)
  240 format (12x,'chi-square (',i2,' d.f.) =',f10.3/
     1 12x,'observation in predicted confidence region ?',2x,a4)
  250 format (/1h0,'observation outside of single-step confidence region
     1 ',i6,' times in ',f6.0,' trials'/1h0,'error rate is ',g12.5)
  260 format (/1h0,'likelihood equation for mu'/(1h ,1p6d16.6))
  270 format (/1h0,'likelihood equation for gamma'/(1h ,1p6d16.6))
  280 format (/1h0,'likelihood equation for lambda'/(1h ,1p6d16.8))
  290 format (1h0,'log-likelihood is',g15.8)
      end