*...:....1....:....2....:....3....:....4....:....5....:....6....:....7....:....8
*    
*     density.f
*
*     Copyright (C) 1991.
*
*     Marie Davidian and A. Ronald Gallant 
*     P.O. Box 5513, Raleigh NC 27650-5513, USA   
*
*     Permission to use, copy, modify, and distribute this
*     software and its documentation for any purpose and
*     without fee is hereby granted, provided that the
*     above copyright notice appear in all copies and that
*     both that copyright notice and this permission notice
*     appear in supporting documentation.
*
*     This software is provided "as is" without any
*     expressed or implied warranty.
*
*...:....1....:....2....:....3....:....4....:....5....:....6....:....7....:....8
*
*     This file can be used as a templete in many pharmacokinetic 
*     applications; edit the lines marked with an X in column 73.

      subroutine density(isw)
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      save

      include "global.f"

      real*8 time,conc
      real*8 data
      real*8 Cl,V
      real*8 P,PwCl,PwV,PwS,Plag,PwClag,PwVlag,PwSlag
      real*8 E,O,OwCl,OwV,S,f,fwO,fwS
      real*8 Clwta1,Clwta2,Clwta3                                       X
      real*8 Vwtau1,Vwtau2,Vwtau3                                       X
      real*8 Swtau1,Swtau2,Swtau3                                       X
      real*8 Clwz1,Clwz2                                                X
      real*8 Vwz1,Vwz2                                                  X
      
      integer*4 i,isw
      integer*4 id0,id0lag,iptr
      integer*4 count,point

      real*8 vpow
      integer*4 mrows,mcols

      parameter (vpow=2.d0)                                             X
      parameter (mrows=2,mcols=2000)                                    X

      dimension count(mcols),point(mcols)
      dimension data(mrows,mcols)


*     If isw=1, read data; compute nobs, pointers, and counts.

      if (isw.eq.1) then

        open (unit=14,file='monoexp.dat',status='old',form='formatted') X

        nobs=0
        iptr=1
        count(iptr)=0
        point(iptr)=1
        id0lag=1

        do i=1,mcols
        
          read(14,14001,end=10) id0,time,conc                           X

          data(1,i)=time                                                X
          data(2,i)=conc                                                X

          nobs=nobs+1                                                   X

          if ((id0lag.eq.id0).or.(i.eq.1)) then
            count(iptr)=count(iptr)+1
          else
            iptr=iptr+1
            count(iptr)=1
            point(iptr)=i
            id0lag=id0
          end if
        end do
        
10      continue

        if (iptr.ne.ncases) then
          write(3,3901)
          stop
        end if

        close(14)
        return
      end if
      
      
*     If isw=2, compute the density.

      if (isw.eq.2) then

        Cl = dexp(tau(1)+z(1))                                          X
        V  = dexp(tau(2)+z(2))                                          X

        Clwz1 = Cl                                                      X
        Clwz2 = 0.d0                                                    X
        
        Vwz1  = 0.d0                                                    X
        Vwz2  = V                                                       X
 
        Clwta1 = Cl                                                     X
        Clwta2 = 0.d0                                                   X
        Clwta3 = 0.d0                                                   X

        Vwtau1 = 0.d0                                                   X
        Vwtau2 = V                                                      X
        Vwtau3 = 0.d0                                                   X

        Plag   = 1.d0
        PwClag = 0.d0
        PwVlag = 0.d0
        PwSlag = 0.d0

        do i=point(id),point(id)-1+count(id)
          time=data(1,i)                                                X
          conc=data(2,i)                                                X
          E = dexp(-(Cl/V)*time)
          O = (1.d0/V)*E
          OwV = (-1.d0/V**2)*E + O*(Cl/V**2)*time
          OwCl = O*(-1.d0/V)*time
          S = tau(3)                                                    X
          Swtau1 = 0.d0                                                 X
          Swtau2 = 0.d0                                                 X
          Swtau3 = 1.d0                                                 X
          call normalh(conc,O,S,vpow,f,fwO,fwS)
          P    = Plag*f
          PwS  = Plag*fwS + PwSlag*f
          PwCl = Plag*fwO*OwCl + PwClag*f
          PwV  = Plag*fwO*OwV + PwVlag*f
          Plag   = P
          PwSlag = PwS
          PwClag = PwCl
          PwVlag = PwV
        end do

        d = P

        dwtau(1) = PwCl*Clwta1 + PwV*Vwtau1 + PwS*Swtau1                X
        dwtau(2) = PwCl*Clwta2 + PwV*Vwtau2 + PwS*Swtau2                X
        dwtau(3) = PwCl*Clwta3 + PwV*Vwtau3 + PwS*Swtau3                X

        dwz(1) = PwCl*Clwz1 + PwV*Vwz1                                  X
        dwz(2) = PwCl*Clwz2 + PwV*Vwz2                                  X
        
        return
      end if


*     If isw=3, write random coefficients to unit 9.

      if (isw.eq.3) then
        write(9,9001) id,Cl,V
        return
      end if
        
        
*     If isw=4, print additional information if any.

      if (isw.eq.4) then
        return
      end if
      
      return
      
3901  format(1x,'error, density, error reading data')
9001  format(i5,2e30.18)
14001 format(i5,f8.0,d30.18)

      end
