***   This program estimates the parameters of the trouser trawl   ***
***   model proposed by Millar and Walsh (1990). The selection     ***
***   curve is assumed to be the 2 parameter logistic function.    ***
***   Two models can be fitted: the first assumes that the split   ***
***   of fish in the legs of the trawl is 50:50.  The second       ***
***   estimates the split.                                         ***
***   The user may specify that either or both models be fitted.   ***
***                                                                ***
***   The user is prompted for the name of the file containing the ***
***   trouser trawl length frequency data (a maximum of 100 length ***
***   classes).  Each row contains 3 observations -                ***
***   the length class, the number of fish in the wide (test) mesh ***
***   leg, and the number of fish in the fine mesh leg.            ***
***   Computed statistics include the parameters of the fitted     ***
***   logistic selectivity curve, estimates of L.25, L.50 and L.75 ***
***   and their (asymptotic) variances.  These statistics are      ***
***   written to the file params.out.  The fitted selectivity      *** 
***   curve (retention probabilities) and the deviance residuals   ***    
***   are written to rprobs.out.                                   ***
***                                                                ***
***   The maximization algorithm is Newton-Raphson. The user must  ***
***   specify the maximum number of iterations and the convergence ***
***   tolerance (relative absolute difference between successive   ***
***   parameter estimates), epsilon.  Initial parameter values     ***
***   are derived from user supplied initial guesses of L.25 and   ***
***   L.50                                                         ***
***                                                                *** 
***   This program's name is ttrawl.f. It was written by N. Cad-   ***
***   igan, Dept. of Fisheries and Oceans, Government of Canada,   ***
***   St. John's, Newfoundland, Canada, in the fall of 1990. The   ***
***   program is written in standard fortran and was developed on  ***
***   a SUN workstation (UNIX) using a SUN compiler. The program   ***
***   has been successfully compiled and executed on VAX-VMS main- ***
***   frame and on a PC-XT using a Microsoft Fortran compiler.     ***

	real*8 n1(100),n2(100),np(100),l(100),l25,l50,xeps
	integer nlclass,miter,izero

	real*8 param(3),vparam(3,3),lp(3),vlp(3)
	real*8 aeps,rll,fll,dll(100),r,c
	integer niter,model
        character*50 infile
	common /data/ n1,n2,np,xeps,l25,l50,miter,nlclass,l

        call spaces(15)
        write(6,*)'        Input name of data file'
        call spaces(15)
        read(5,*) infile
	open(1,file = infile,status='old')
	open(2,file = 'params.out')
	open(3,file = 'rprobs.out')

***     l(i) is the vector of length classes.
***     n1(i) and n2(i) contain the length frequency data; np(i) is
***     their sum (i = 1(1)nlclass).
***     nlclass is the number of length classes.
***     izero records the number of groups with np(i) = 0.


	izero = 0
	nlclass = 1
99	read(1,*,end=100)l(nlclass),n1(nlclass),n2(nlclass)
	np(nlclass) = 0.0d0
	np(nlclass) = n1(nlclass) + n2(nlclass)
	if (np(nlclass).eq.0.0) izero = izero + 1
	nlclass = nlclass + 1
        goto 99
100     continue
	nlclass = nlclass - 1
	
	call spaces(15)
500	write(6,*)'         Input initial guesses of L.25 and L.50'
	call spaces(15)
	read(5,*)l25,l50
	if (l25 .gt. l50) then
	  write(6,*)
	  write(6,*)'         Initial guesses are incorrect'
	  write(6,*)'         L.25 must be less than L.50'
	  write(6,*)'         Re-enter guesses'
	  write(6,*)
	  goto 500
	end if

***     starting values of a and b are computed from l.25 and l.50   

	call spaces(15)
501	write(6,*)'      Input convergence tolerance and maximum'
        write(6,*)'      number of iterations'
	call spaces(15)
	read(5,*)xeps,miter
	if (xeps .le. 0) then
	  write(6,*)
	  write(6,*)'         Epsilon must be > 0'
	  write(6,*)'         Re-enter guesses'
	  write(6,*)
	  goto 501
	end if 
	if (miter .lt. 0) then
	  write(6,*)'         Maximum number of iterations must be >= 0'
	  write(6,*)'         Re-enter guesses'
	  write(6,*)
	  goto 501
	end if
	call spaces(15)
502	write(6,*)'                Choose which model(s) to fit'
	write(6,*)
	write(6,*)
	write(6,*)'                Enter 1,2 or 3'
	write(6,*)
	write(6,*)'                1 : 50:50 split'
	write(6,*)
	write(6,*)'                2 : Estimated split'
	write(6,*)
	write(6,*)'                3 : Both 1 and 2'
	call spaces(5)
	read(5,*)model
	if ((model .lt. 1).or.(model .gt. 3)) then
	  write(6,*)
	  write(6,*)'         Invalid choice'
	  write(6,*)'         enter a 1, 2 or 3'
	  write(6,*)
	  goto 502
	end if

	if ((model.eq.1).or.(model.eq.3)) then
          goto 101
	else
	  goto 102
	end if
	  
***     This portion fits the 2 parameter model.

101	call twop(param,vparam,niter,aeps,rll,fll,dll)
	  
***     this portion computes the 25%, 50% and 75% retention lengths
***     lp(i) (i = 1(1)3) and their variances vlp(i) (i = 1(1)3).

	r = 0.25d0
	do 10 i = 1,3
	  c = dlog(r/(1.0d0-r))
	  vlp(i) = (param(2)**(-4))*(param(2)*param(2)*vparam(1,1) - 
     a    (2.0d0*param(2)*(param(1)-c)*vparam(1,2))+
     a    ((param(1)-c)**2)*vparam(2,2))
	  lp(i) = (c-param(1))/param(2)
  	  r = r + 0.25d0
10	continue

	write(2,*)
	write(2,*)'         Statistics for the two parameter model'
	write(2,*)
	if (aeps.gt.xeps) then
	  write(2,*)'         Maximum number of iterations (',miter,')'
     a              ,'reached'
	else
	  write(2,*)'         Convergence criteria met in ',niter-1,
     a              ' iterations'
	end if 
	write(2,*)
	write(2,*)'         Parameter Estimates and Covariances'
	write(2,1002)'a         = ', param(1)
	write(2,1002)'b         = ', param(2)
	write(2,1002)'var(a)    = ',vparam(1,1)
	write(2,1002)'var(b)    = ',vparam(2,2)
	write(2,1002)'cov(a,b)  = ',vparam(1,2)
	write(2,1002)'L.25      = ',lp(1)
	write(2,1002)'L.50      = ',lp(2)
	write(2,1002)'L.75      = ',lp(3)
	write(2,1002)'var(L.25) = ', vlp(1)
	write(2,1002)'var(L.50) = ', vlp(2)
	write(2,1002)'var(L.75) = ', vlp(3)
	write(2,1002)'Loglikelihood of full model    = ',fll
	write(2,1002)'Loglikelihood of reduced model = ',rll
	write(2,*)
	write(2,1003)'Model deviance (',nlclass-2-izero,'df) =',
     a  2.0d0*(fll-rll)
	write(3,*)
	write(3,*)'             Retention probabilities for 2',
     a            '-parameter model'
	write(3,*)
	write(3,*)'             Length                            Deviance'
	write(3,*)'             Group           r(l)              Residual'
	write(3,*)

***     this portion computes the retention probabilities for each
***     length group.

	do 11 i = 1,nlclass
	  r = dexp(param(1) + param(2)*l(i))/(1.0d0 + dexp(param(1) 
     a    + param(2)*l(i)))
	  write(3,1000)l(i),r,dll(i)
11	continue

	if (model.eq.3) then
	  goto 102
	else
	  goto 103
	end if

***     This portion fits the three parameter model

102	call threep(param,vparam,niter,aeps,rll,fll,dll)
	  
***     this portion computes the 25%, 50% and 75% retention lengths
***     lp(i) (i = 1(1)3) and their variances vlp(i) (i = 1(1)3).
		  
	r = 0.25d0
	do 12 i = 1,3
	  c = dlog(r/(1.0d0 - r))
	  vlp(i) = (param(2)**(-4))*(param(2)*param(2)*vparam(1,1) - 
     a    (2.0d0*param(2)*(param(1)-c)*vparam(1,2))+
     a    ((param(1)-c)**2)*vparam(2,2))
	  lp(i) = (c-param(1))/param(2)
  	  r = r + 0.25d0
12	continue

	write(2,*)
	write(2,*)
        write(2,*)'______________________________________________'
	write(2,*)
	write(2,*)'         Statistics for the three parameter model'
	write(2,*)
	if (aeps.gt.xeps) then
	  write(2,*)'         Maximum number of iterations (',miter,')'
     a             ,' reached' 	  
	else
	  write(2,*)'         Convergence criteria met in ',niter-1,
     a              ' iterations'
	end if 
	write(2,*)
	write(2,*)'         Parameter Estimates and Covariances'
	write(2,1002)'a         = ', param(1)
	write(2,1002)'b         = ', param(2)
	write(2,1002)'p         = ', param(3)
	write(2,1002)'var(a)    = ',vparam(1,1)
	write(2,1002)'var(b)    = ',vparam(2,2)
	write(2,1002)'var(p)    = ',vparam(3,3)
	write(2,1002)'cov(a,b)  = ',vparam(1,2)
	write(2,1002)'cov(a,p)  = ',vparam(1,3)
	write(2,1002)'cov(b,p)  = ',vparam(2,3)
	write(2,1002)'L.25      = ',lp(1)
	write(2,1002)'L.50      = ',lp(2)
	write(2,1002)'L.75      = ',lp(3)
	write(2,1002)'var(L.25) = ', vlp(1)
	write(2,1002)'var(L.50) = ', vlp(2)
	write(2,1002)'var(L.75) = ', vlp(3)
	write(2,1002)'Loglikelihood of full model    = ',fll
	write(2,1002)'Loglikelihood of reduced model = ',rll
	write(2,*)
	write(2,1003)'Model deviance (',nlclass-3-izero,'df) =',
     a  2.0d0*(fll-rll)
	write(3,*)
	write(3,*)'             Retention probabilties for 3',
     a            '-parameter model'
	write(3,*)
	write(3,*)'           Length                              Deviance'
	write(3,*)'           Group             r(l)              Residual'
	write(3,*)

***     this portion computes the retention probabilities for each
***     length group.

	do 13 i = 1,nlclass
	  r = dexp(param(1) + param(2)*l(i))/(1.0d0 + 
     a       dexp(param(1) + param(2)*l(i)))
	  write(3,1000)l(i),r,dll(i)
13	continue

        write(6,*)
        write(6,*)' Program finished:'
        write(6,*)'  - estimated statistics are in file params.out'
        write(6,*)'  - estimated selectivity curve(s) and deviance ',
     a            'residuals'
        write(6,*)'    are in file rprobs.out'

1000    format(11x,f7.3,5x,e14.6,5x,e14.6)
1001    format(/,a37)
1002    format(/,a40,3x,e14.6)
1003    format(a28,1x,i3,1x,a6,4x,e14.6)
103	stop
	end

***____________________________________________________________________

	subroutine twop(p,v,iter,eps,rl,fl,dll)

***     This routine computes the estimates of a and b (p - fixed at
***      0.5) for the two parameter model. Parameter estimates are 
***     returned invector p, their variances in v. The loglikelihood
***     of the fitted model is returned in rl; that of the full model 
***     in fl. The number of iterations to convergence and the last  
***     computed epsilon are returned in iter and eps respectively.
***     Note: dl a 2x1 vector of 1st order partial derivatives of the
***	      likelihood function.(regardless of dimension).
***           xj a 2x2 matrix of 2nd order partial derivatives of the
***           likelihood function.(regardless of dimension).
***           t's are temporary variables.
***           nd is the dimension of the problem i.e. 2.

	real*8 n1(100),n2(100),np(100),l(100),l25,l50,xeps
	integer nlclass,miter

	real*8 p(3),dl(3),xj(3,3),v(3,3)
	real*8 eps,rl,fl,t1,t2,t3,t4,t5,det,dll(100),trl,tfl
	integer iter,nd
	common /data/ n1,n2,np,xeps,l25,l50,miter,nlclass,l

	nd = 2
	eps = 1.0d0
	iter = 0
	p(2) = -1.0986123d0/(l25-l50)
	p(1) = -p(2)*l50

10	if ((eps.gt.xeps).and.(iter.le.miter)) then

 	  do 11 i = 1,3
	    dl(i)=0.0d0
	    xj(1,i)=0.0d0
	    xj(2,i)=0.0d0	
	    xj(3,i)=0.0d0
11	  continue

	  do 12 i = 1,nlclass
	    t1 = 1.0/(2.0d0 + dexp(-p(1)-p(2)*l(i)))
	    t2 = 1.0d0 - t1
	    t3 = 1.0d0 - 2.0d0*t1
	    t4 = t1*t3/(t2**2)
	    t5 = 2.0d0*(t1**2) - 4.0d0*t1 + 1.0d0

	    dl(1) = dl(1) + t3*(n1(i) - n2(i)*t1/t2)

	    dl(2) = dl(2) + l(i)*t3*(n1(i) - n2(i)*t1/t2)

	    xj(1,1) = xj(1,1) - t4*(n1(i) + np(i)*t5)

	    xj(2,2) = xj(2,2) - l(i)*l(i)*t4*(n1(i) + np(i)*t5)

	    xj(1,2) = xj(1,2) - l(i)*t4*(n1(i) + np(i)*t5)

12	  continue
	  xj(2,1) = xj(1,2)

	  call mxint(xj,nd,3,det)
	  if(det.eq.0.0d0) then
             write(6,*)'solution diverges, try different'
	     write(6,*)'starting values of L.25 and L.50'
	     stop
	  end if
	  p(1) = p(1) - xj(1,1)*dl(1) - xj(1,2)*dl(2)
	  p(2) = p(2) - xj(2,1)*dl(1) - xj(2,2)*dl(2)

	  eps = dabs(max((xj(1,1)*dl(1) + xj(1,2)*dl(2))/p(1),
     a    (xj(2,1)*dl(1) + xj(2,2)*dl(2))/p(2)))

	  iter = iter + 1
	goto 10
	endif
	
	rl = 0.0d0
	fl = 0.0d0
	do 14 i = 1,3
	  do 13 j = 1,3
	    v(i,j) = 0.0d0
13	  continue
14	continue

	do 15 i = 1,nlclass
	    dll(i) = 0
	    t1 = 1.0d0/(2.0d0 + dexp(-p(1)-p(2)*l(i)))
	    rl = rl + n1(i)*dlog(t1) + n2(i)*dlog(1.0d0-t1)
	    trl = n1(i)*dlog(t1) + n2(i)*dlog(1.0d0-t1)
	    tfl = 0
	    t2 = 0
	    if (np(i).gt.0.0) t2 = n1(i)/np(i)
	    if ((n1(i).gt.0.0).and.(n2(i).gt.0.0)) then
	      t2 = n1(i)/np(i) 
	      fl = fl + n1(i)*dlog(t2) + n2(i)*dlog(1.0d0-t2)
	      tfl = n1(i)*dlog(t2) + n2(i)*dlog(1.0d0-t2)
	    end if
	    dll(i) = sign(1,t2-t1)*dsqrt(2.0d0*(tfl-trl))

	  t3 = np(i)*dexp(p(1)+p(2)*l(i))/((1.0d0+dexp(p(1)
     a    +p(2)*l(i)))*((1.0d0 + 2.0d0*dexp(p(1)+p(2)*l(i)))**2))
	  v(1,1) = v(1,1) + t3
	  v(1,2) = v(1,2) + l(i)*t3
	  v(2,2) = v(2,2) + l(i)*l(i)*t3
15	continue
	v(2,1) = v(1,2)

	call mxint(v,nd,3,det)

	return
	end

*** __________________________________________________________________

	subroutine threep(p,v,iter,eps,rl,fl,dll)

***     This routine computes the estimates of a, b and p (the 
***     proportion of fish entering the wide mesh).
***     Parameter estimates are returned in vector p, their 
***     variances in v. The loglikelihood of the estimated model
***     is returned in rl; that of the full model in fl. The number of
***     iterations to convergence and the last computed epsilon are
***     returned in iter and eps respectively.

***     Note: dl a 3x1 vector of 1st order partial derivatives of the
***	      likelihood function.
***           xj a 3x3 matrix of 2nd order partial derivatives of the
***           likelihood function.
***           t's are temporary variables.
***           nd is the dimension of the problem i.e. 3.

	real*8 n1(100),n2(100),np(100),l(100),l25,l50,xeps
	integer nlclass,miter

	real*8 p(3),dl(3),xj(3,3),v(3,3)
	real*8 eps,rl,fl,t1,t2,t3,t4,det,dll(100),trl,tfl
	integer iter,nd
	common /data/ n1,n2,np,xeps,l25,l50,miter,nlclass,l

	nd = 3
	eps = 1.0d0
	iter = 0
	p(2) = -1.0986123d0/(l25-l50)
	p(1) = -p(2)*l50
	p(3) = 0.5

10	if ((eps.gt.xeps).and.(iter.le.miter)) then

 	  do 11 i = 1,3
	    dl(i)=0.0d0
	    xj(1,i)=0.0d0
	    xj(2,i)=0.0d0
	    xj(3,i)=0.0d0
11	  continue

	  do 12 i = 1,nlclass
	    t1 = p(3)/(1.0d0 + (1.0d0 - p(3))*dexp(-p(1)-p(2)*l(i)))
	    t2 = 1.0d0 - t1
	    t3 = p(3) - t1

	    t4 = t3/p(3)
	    dl(1) = dl(1) + t4*((n2(i)*t1/t2) - n1(i))
	    dl(2) = dl(2) + l(i)*t4*((n2(i)*t1/t2) - n1(i))

	    t4 = p(3)*(1.0d0 - p(3))
	    dl(3) = dl(3) + t4*((t1*np(i)) - n1(i))

	    t4 = t1*t3/(t2*p(3)*p(3))
	    xj(1,1) = xj(1,1) + t4*(n1(i) - np(i)*t1 + (n2(i)*t3/t2))
	    xj(2,2) = xj(2,2) + l(i)*l(i)*t4*(n1(i) - np(i)*t1 +
     a                (n2(i)*t3/t2))
	    xj(1,2) = xj(1,2) + l(i)*t4*(n1(i) - np(i)*t1 +
     a                (n2(i)*t3/t2))

	    t4 = p(3)*p(3)/((1.0d0-p(3))**2)
	    xj(3,3) = xj(3,3) + t4*((2.0d0*p(3)*(np(i)*t1 - n1(i)))
     a      + n1(i) - (np(i)*t1*t1))  
	    
	    t4 = p(3)*p(3)/(1.0d0-p(3))
	    xj(1,3) = xj(1,3) + t4*t1*t3*np(i)
	    xj(2,3) = xj(2,3) + l(i)*t4*t1*t3*np(i)
12	  continue
	  xj(2,1) = xj(1,2)
	  xj(3,1) = xj(1,3)
	  xj(3,2) = xj(2,3)

	  call mxint(xj,nd,nd,det)
	  if(det.eq.0.0d0) then
             write(6,*)'solution diverges, try different'
	     write(6,*)'starting values of L.25 and L.50'
	     stop
	  end if

	  p(1) = p(1) - xj(1,1)*dl(1) - xj(1,2)*dl(2) - xj(1,3)*dl(3)
	  p(2) = p(2) - xj(2,1)*dl(1) - xj(2,2)*dl(2) - xj(2,3)*dl(3)
	  p(3) = p(3) - xj(3,1)*dl(1) - xj(3,2)*dl(2) - xj(3,3)*dl(3)

	  eps = dabs(max((xj(1,1)*dl(1)+xj(1,2)*dl(2)+
     a   xj(1,3)*dl(3))/p(1),
     a    (xj(2,1)*dl(1)+xj(2,2)*dl(2)+xj(2,3)*dl(3))/p(2),
     a   (xj(3,1)*dl(1)+xj(3,2)*dl(2)+xj(3,3)*dl(3))/p(3)))

	  iter = iter + 1
	goto 10
	endif
	
	rl = 0.0d0
	fl = 0.0d0
	do 14 i = 1,3
	  do 13 j = 1,3
	    v(i,j) = 0.0d0
13	  continue
14	continue

	do 15 i = 1,nlclass
	  dll(i) = 0
	  t1 = p(3)*dexp(p(1)+p(2)*l(i))/(1.0d0 - p(3) + dexp(p(1) 
     a         + p(2)*l(i)))
	  t2 = p(3)*p(3)*(1.0d0-t1)
	  t3 = p(3)*p(3)*(1.0d0-p(3))

	    rl = rl + n1(i)*dlog(t1) + n2(i)*dlog(1.0d0-t1)
	    trl = n1(i)*dlog(t1) + n2(i)*dlog(1.0d0-t1)
	    tfl = 0
	    t4 = 0
	    if (np(i).gt.0.0) t4 = n1(i)/np(i)
	    if ((n1(i).gt.0.0).and.(n2(i).gt.0.0)) then
	      t4 = n1(i)/np(i) 
	      fl = fl + n1(i)*dlog(t4) + n2(i)*dlog(1.0d0-t4)
	      tfl = n1(i)*dlog(t4) + n2(i)*dlog(1.0d0-t4)
	    end if
	    dll(i) = sign(1,t4-t1)*dsqrt(2.0d0*(tfl-trl))

	  v(1,1) = v(1,1) + np(i)*((p(3) - t1)**2)*t1/t2
	  v(2,2) = v(2,2) + l(i)*l(i)*np(i)*((p(3) - t1)**2)*t1/t2
	  v(2,1) = v(2,1) + l(i)*np(i)*((p(3)-t1)*(p(3)-t1))*t1/t2

	  v(3,3) = v(3,3) + np(i)*(1.0d0 - t1)*t1/(t3*(1.0d0 - p(3)))

	  v(3,1) = v(3,1) + np(i)*(p(3) - t1)*t1/t3
	  v(3,2) = v(3,2) + l(i)*np(i)*(p(3) - t1)*t1/t3
15	continue

	v(1,2) = v(2,1)
	v(1,3) = v(3,1)
	v(2,3) = v(3,2)

	call mxint(v,nd,nd,det)

	return
	end

*** ____________________________________________________________________

	subroutine  mxint(a,nrr,nra,det)
        dimension isw(96), a(*) 
        double precision a,det

***     this matrix inversion subroutine takes a  nra*nra
***     square matrix and inverts the upper  nrr*nrr
***     submatrix 
 
        na = nra+1
        nd = nrr*nra
        do 1 i=1,nrr
    1   isw(i) = i
        ndiag = 1
        nup = 1
        det = 1.
        do 9 i=1,nrr
           ka = i
           am = 0.
           la = ndiag 
           do 2 j=i,nrr
              if (dabs(a(la)).le.am) go to 2
              ka = j
              am = dabs(a(la))
    2      la = la+1
           if (am.gt.1.0e-20) go to 3
	   write(6,*)'pivot equal zero'
           det = 0.
           return
    3      if (ka.eq.i) go to 5 
           mup = i
           do 4 mlow=ka,nd,nra
              s = a(mup)
              a(mup) = a(mlow)
              a(mlow) = s
    4      mup = mup+nra
           is = isw(ka)
           isw(ka) = isw(i)
           isw(i) = is
    5      aa = a(ndiag)
           det = det*aa
           do 6 j=i,nd,nra
    6      a(j) = a(j)/aa
           a(ndiag) = 1./aa
           nel = nup
           do 8 j=1,nrr
              if (j.eq.i) go to 8
              aa = a(nel)
              a(nel) = 0.
              ka = i
              do 7 k=j,nd,nra
                  a(k) = a(k)-aa*a(ka)
    7         ka = ka+nra
    8      nel = nel+1
           nup = nup+nra
    9   ndiag = ndiag+na
        i = 1
   10   if (isw(i).ne.i) go to 11
        i = i+1
        if (i.gt.nrr) go to 13
   11   j = isw(i)
        isw(i) = isw(j)
        isw(j) = j
        mlow = (i-1)*nra+1
        mup = (j-1)*nra+1
        do 12 k=1,nrr 
           s = a(mlow)
           a(mlow) = a(mup)
           a(mup) = s 
           mlow = mlow+1
   12   mup = mup+1
        go to 10
 
   13   return
        end

***  ________________________________________________________________

	subroutine spaces(k)
	integer k
	  do 10 i = 1,k
	    print*,' ' 
10	  continue
	return
	end