anchor.diff:        
   c
         subroutine logit(theta,n,ltheta)
  -c Truncates theta to between 0 and 1.
  +c Maps (0,1) -> R.
   cf2py intent(hide) n
   cf2py intent(out) ltheta
         DOUBLE PRECISION theta(n), ltheta(n)
  
positive.code:       end do
        RETURN
        END
  c
        subroutine invlogit(ltheta,n,theta)
  c Truncates theta to between 0 and 1.
  cf2py intent(hide) n
  cf2py intent(out) theta
        DOUBLE PRECISION theta(n), ltheta(n)
  
negative.code:       end do
        RETURN
        END      
  
  c
  
        subroutine stukel_logit(theta,n,ltheta,a1,a2,na1,na2)
  !
  ! Reference: Therese A. Stukel, 'Generalized Logistic Models',
  ! JASA vol 83 no 402, pp.426-431 (June 1988)
  !
  cf2py intent(hide) n, na1, na2
  cf2py intent(out) ltheta
  cf2py intent(copy) theta
        DOUBLE PRECISION theta(n), ltheta(n)
        DOUBLE PRECISION a1(na1), a2(na2), a1t, a2t
        LOGICAL a1_isscalar, a2_isscalar      
        DOUBLE PRECISION infinity
        PARAMETER (infinity = 1.7976931348623157d308)      
        INTEGER n, i, na1, na2
  
        a1t = a1(1)
        a2t = a2(1)
  
        CALL logit(theta,n,ltheta)
  
        a1_isscalar = (na1.LT.n)
        a2_isscalar = (na2.LT.n)      
  
        do i=1,n
            
            if (ltheta(i).GT.0.0D0) then
                if (.NOT.a1_isscalar) then
                    a1t = a1(i)
                end if
                if (a1t.GT.0.0D0) then
                    ltheta(i)=dlog(ltheta(i)*a1t+1.0D0)/a1t
                else if (a1t.LT.0.0D0) then
                    ltheta(i) = (1.0D0-dexp(-ltheta(i)*a1t))/a1t
                end if
                
            else if (ltheta(i).LT.0.0D0) then
                if (.NOT.a2_isscalar) then
                    a2t = a2(i)
                end if
                if (a2t.GT.0.0D0) then
                    ltheta(i)=-dlog(-ltheta(i)*a2t+1.0D0)/a2t
                else if (a2t.LT.0.0D0) then
                    ltheta(i)=-(1.0D0-dexp(ltheta(i)*a2t))/a2t
                end if
                
            else
                ltheta(i) = 0.0D0
            end if
  
        end do
        
        RETURN
        END
  
  c
  
        subroutine stukel_invlogit(ltheta,n,theta,a1,a2,na1,na2)
  !
  ! Reference: Therese A. Stukel, 'Generalized Logistic Models',
  ! JASA vol 83 no 402, pp.426-431 (June 1988)
  !
  cf2py intent(hide) n, na1, na2
  cf2py intent(out) theta
  cf2py intent(copy) ltheta
        DOUBLE PRECISION theta(n), ltheta(n)
        DOUBLE PRECISION a1(na1), a2(na2), a1t, a2t
        LOGICAL a1_isscalar, a2_isscalar
        DOUBLE PRECISION infinity
        PARAMETER (infinity = 1.7976931348623157d308)      
        INTEGER n, i, na1, na2
  
        a1t = a1(1)
        a2t = a2(1)
        
        a1_isscalar = (na1.LT.n)
        a2_isscalar = (na2.LT.n)      
        
        do i=1,n
            if (ltheta(i).GT.0.0D0) then
                if (.NOT.a1_isscalar) then
                    a1t = a1(i)
                end if
                if (a1t.GT.0.0D0) then
                    ltheta(i) = (dexp(a1t*ltheta(i))-1.0D0)/a1t
                else if (a1t.LT.0.0D0) then
                    ltheta(i) = -dlog(1.0D0-a1t*ltheta(i))/a1t
                end if
                
            else if (ltheta(i).LT.0.0D0) then
                if (.NOT.a2_isscalar) then
                    a2t = a2(i)
                end if
                if (a2t.GT.0.0D0) then
                    ltheta(i) = -(dexp(-a2t*ltheta(i))-1.0D0)/a2t
                else if (a2t.LT.0.0D0) then
                    ltheta(i) = dlog(1.0D0+a2t*ltheta(i))/a2t
                end if
                
            else
                ltheta(i) = 0.5D0
            end if
  
        end do
  
        
        CALL invlogit(ltheta,n,theta)
        
        RETURN
        END      
  