*...:....1....:....2....:....3....:....4....:....5....:....6....:....7....:....8
*    
*     normalh.f
*
*     Copyright (C) 1991,1992.
*
*     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
    
*     Normalh computes the normal density evaluated at x and its deriv-
*     atives with respect to O and S.  The density is returned in f and
*     the derivatives in fwO and fwS.  O and vpow must be positive.

*     The mean of the density is O.  The variance of the density is
*     (S**2)*(O**vpow); that is, the variance is proportional to the
*     mean raised to the vpow power.

*     Normalh.f was substantially revised 3/28/92 to robustify the code
*     against negative and small values of O.   

*     The parameter eps is machine dependent.

*     Set the parameter eps to about 10,000 times the smallest real*8
*     that the machine can represent by editing the statement
*     parameter (eps=1.d-xxx)
*     For a machine with IEEE floats, such as a Sun SparcStation2, set
*     eps=1.d-300.  For an IBM System/370 machine, such as an IBM 3081,
*     set eps=1.d-80.  (The code has not been tested on IBM equipment.)

*     The value 10,000 = 100**2 above represents an assumption that 
*     normalh will not be called with x exceeding 100 in absolute value
*     when O is small.  If this assumption is violated, the computation 
*     of RV*(x-OO)**2 in the statement
*     logf = - 0.5d0*(RV*(x-OO)**2 + logV + log2pi)
*     may fail.

*     If x and S are such that RV and RV*(x-OO)**2 can be computed, 
*     then the code should not fail.

      subroutine normalh(x,O,S,vpow,f,fwO,fwS)
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      save

      real*8 x,O,S,vpow,f,fwO,fwS
      real*8 OO,SS,tol
      real*8 RV,logV,logf
      real*8 D,RVD,RVDD,fRVD,fRVDD

      integer*4 OOwO

      integer*4 ieeer,ieee_flags
      character*32 out
      
      real*8 log2pi,eps,vpsave 

      parameter (eps=1.d-300)
      parameter (log2pi=0.1837877066409345D+01)

      data vpsave /-500.d0/

      ieeer = ieee_flags('clear','exception','common',out)

      SS=dabs(S)

      if (vpow.ne.vpsave) then
        vpsave=vpow
        if (vpow.gt.1.d0) then
          tol=dexp(dlog(eps)/vpow)
        else
          tol=eps
        end if
      end if

      if (O.gt.tol) then
        OO=O
        OOwO=1
      else
        OO=tol
        OOwO=0
      end if
        
      logV=2.d0*dlog(SS)+vpow*dlog(OO)

      RV=dexp(-logV)

*     These three lines order the computations to protect against large
*     values of O.

      D=(x-O)
      RVD=RV*D
      RVDD=RVD*D

      logf = - 0.5d0*(RVDD + logV + log2pi)

      f=dexp(logf)

*     These two lines use small f to reduce the magnitude of RVDD prior 
*     to the division by OO below which helps when OO is small.

      fRVD=f*RVD
      fRVDD=fRVD*D

      if (OOwO.eq.1) then
        fwO = fRVD + 0.5d0*vpow*(fRVDD - f)/OO
      else
        fwO = 0.d0
      end if

      if (S.ge.0.d0) then
        fwS = (fRVDD - f)/SS
      else
        fwS = - (fRVDD - f)/SS
      end if

      ieeer = ieee_flags('get','exception','invalid',out)
      if (out.eq.'invalid') write(6,6901) out,x,O,S,vpow,f,fwO,fwS
      ieeer = ieee_flags('get','exception','overflow',out)
      if (out.eq.'overflow') write(6,6901) out,x,O,S,vpow,f,fwO,fwS
      ieeer = ieee_flags('get','exception','division',out)
      if (out.eq.'division') write(6,6901) out,x,O,S,vpow,f,fwO,fwS
      ieeer = ieee_flags('clear','exception','common',out)
      return
6901  format(a32,' x,O,S,vpow,f,fwO,fwS'/7d12.3)
      end
