PROGRAM quadpdf

c     using the quadrature model, compute PDF(t,tbar,sigup,sigdn), the
c     probability of the error distribution as a function of the
c     independent variable T for a measurement with mean value TBAR
c     having asymmetric errors SIGUP and SIGDN.  The computation is
c     performed as a function of X, where X=T-TBAR

c     USER INPUT
c     N     = number of "pixels" in the output function (parameter)
c     tbar  = measured mean value
c     sigup = downward error bar in units of T
c     sigup = upward error bar in units of T
c     xmax  = interval to each side of TBAR over which to compute the PDF

c     OUTPUT
c     x     = independent variable relative to the mean, X=T-TBAR
c     f1    = dimidated Gaussian distribution (for comparison only)
c     f2    = the desried PDF based upon the quadrature model

      implicit none
      logical             error
      integer       i,j,k,n
      parameter     (n=1800)
      double precision sigma,alpha,tpi,sigup,sigdn
      double precision agy,y,r
      double precision u,upos,uneg,arg,value2
      double precision dx,xmax,x,gx,gu,dxdu1,dxdu2,f1,f_pos,f_neg
      double precision t,tbar
      character*30       in1,in2,in3

c     constant 2*pi

      tpi   = 2.0d0*3.14159d0

      CALL getarg(1,in1)
      CALL getarg(2,in2)
      CALL getarg(3,in3)
      READ(in1,*) tbar
      READ(in2,*) sigdn
      READ(in3,*) sigup


      sigup = sigup - tbar
      sigdn = tbar - sigdn

      if (sigdn.eq.sigup) sigup = sigup + 0.001

c     ----------------------------------------------
c     USER INPUTS (can be modified to read in a file)
c     the mean, and the upward/downward "sigmas" relative to the mean

C      tbar  = 45.    ! measured value
C      sigup = 1.9    ! upward error bar
C      sigdn = 2.1    ! downward error bar

c     maximum +/- range of |T-TBAR| over which to compute the PDF

      xmax  = 180.0d0+tbar

c     ----------------------------------------------
c     BEGIN COMPUTATIONS

c     DX = "pixel size" at which to output the PDF; N is declared as a
c     parameter (above); factor of 2 accounts for XMAX range to both
c     sides about the mean
      
      dx    = 0.2

c     quadrature parameters for PDF calculation

      sigma = 0.50d0*(sigup+sigdn)       ! average of sigma
      alpha = 0.50d0*(sigup-sigdn)       ! half difference of sigma

      DO 11 i=1,n                        ! loop over x
       x   = -xmax + float(i-1)*dx       ! run of independent variable
       arg = (sigma**2) + 4.0d0*alpha*x  ! bias/skew term 

c     compute the dimidated Gaussian (for comparisons only)

       IF (x.le.0.0d0) then  ! below the mean
        gx    = exp(-0.50d0*x*x/(sigdn*sigdn))/sqrt(tpi)
        dxdu1 = sigdn
        f1    = gx/dxdu1
       ELSE                  ! above the mean
        gx    = exp(-0.50d0*x*x/(sigup*sigup))/sqrt(tpi)
        dxdu1 = sigup
        f1    = gx/dxdu1
       END IF

c     compute the Barstow quadraditc asymmetric PDF; arg<0 invalid

       f_pos = 0.0d0
       f_neg = 0.0d0

       IF (arg.gt.0.0) then 

c     psotive root
        upos  = (sqrt(arg)-sigma)/(2.0d0*alpha)
        u     = upos
        gu    = exp(-0.50d0*u*u/(sigma*sigma))/sqrt(tpi)
        dxdu2 = abs(sigma + 2.0d0*alpha*u)
        f_pos = gu/dxdu2

c     negative root
        uneg  = (-sqrt(arg)-sigma)/(2.0d0*alpha)
        u     = uneg
        gu    = exp(-0.50d0*u*u/(sigma*sigma))/sqrt(tpi)
        dxdu2 = abs(sigma + 2.0d0*alpha*u)
        f_neg = gu/dxdu2

       END IF

c     compute Kato, Omachi, Aso univartiate asymmetric Gaussian
c     http://www.iic.ecei.tohoku.ac.jp/~kato/papers/t.kato_spr2002a.pdf
c     http://www.springerlink.com/content/yjulk4559lg3k3g9/

       r = sigup/sigdn
       y = x**2/(2.0d0*sigup*sigup)
       if (x.le.0.0) y = y/(r*r) 
       agy = 2.0d0/sqrt(tpi)/sigup/(r+1.0d0) * exp(-y)

c     write to standard out

       IF (f_pos.le.1.0E-20) f_pos = 1.0d-10
       IF (f_neg.le.1.0E-20) f_neg = 1.0d-10
       IF (agy.le.1.0E-20) agy = 1.0d-10

       t = tbar + x
       WRITE(6,*) t,x,f_pos,agy

 11   CONTINUE  ! next x


      STOP
      END