;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function Calc_PAmp, LAmp, PFrac, dPAmp_dLAmp  = dPAmp_dLAmp, $
                                 dPAmp_dPFrac = dPAmp_dPFrac
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;       Calculate Positronium amplitude from the line amplitude
;       and the positronium fraction according to the equation:
;
;       Posit Frac = 2.0 / (2.25*(Line Amp/Posit Amp)+ 1.5)
;
;
;       Parameters:   LAmp  -  511 keV line amplitude
;                     PFrac -  Positronium FRACTION
;
;  Modifications:
;
;       3-Dec-96 (wrp) - New routine.
;
;
if( PFrac eq 2.0/1.5) then begin
  PAmp         = 9.0d9
  dPAmp_dLAmp  = 9.0d9
  dPAmp_dPFrac = 9.0d9
endif else begin
  PAmp         = 2.25d0*PFrac*LAmp / (2.0d0 - 1.5d0*PFrac)
  dPAmp_dLAmp  = 2.25d0*PFrac / (2.0d0 - 1.5d0*PFrac)
  dPAmp_dPFrac = 4.50d0*LAmp / (2.0d0 - 1.5d0*PFrac)^2.0
endelse


return, PAmp
end



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function posit_function, energy, epeak
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;       Positronium Continuum Function
;
;       Reference:  Ore and Powell, Phys. Rev. 75, p. 1696 (1949)
;
;       Parameters:   ENERGY -  Energy array for which positronium
;                               continuum function is to be evaluated.
;                     EPEAK  -  Energy of annihilation line (scalar).
;
;  Modifications:
;
;       29-Nov-93 (wrp) - New routine.  Note that the term 0.869598d0 is
;                         a normalization term so that the integral over
;                         the range 0.0 - 0.511 MeV is unity.
;
;
e0         = double( epeak)
e          = double( energy)
continuum  = dblarr( n_elements( energy))

indx = where( e lt e0, nindx)
if( nindx gt 0) then begin
  e = e( indx)
  continuum(indx) = (2.0/e0) * ( (e*(e0-e)/(2*e0-e)^2) + (2*e0-e)/e + $
             ((2*e0*(e0-e)/e^2) - (2*e0*(e0-e)^2/(2*e0-e)^3)) * alog((e0-e)/e0))
endif

return, float( continuum / 0.869598d0)
end


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function posit_continuum, low_edg, high_edg, epeak
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;       Positronium Continuum Model
;
;       Reference:  Ore and Powell, Phys. Rev. 75, p. 1696 (1949)
;
;       Parameters:   LOW_EDG  -  Array of channel lower energies for
;                                 integrating the positronium continuum
;                     HIGH_EDG -  Array of channel upper energies for
;                                 integrating the positronium continuum
;                     EPEAK    -  Energy of annihilation line (scalar).
;
;  Modifications:
;
;       29-Nov-93 (wrp) - New routine.  Note that the positronium continuum
;                         function is integrated over the channels using
;                         Simpson's Rule.
;
;
e0         = double( epeak)
wid        = double( high_edg) - double( low_edg)
continuum  = dblarr( n_elements( low_edg))

nloop = 11
dwid  = double( wid / float( nloop))
for iloop=1,nloop-1 do begin
  continuum = continuum + 4.0*dwid*posit_function(low_edg+(iloop-0.5)*dwid, e0)
  continuum = continuum + 2.0*dwid*posit_function(low_edg+(iloop)*dwid, e0)
endfor
continuum = continuum + 4.0*dwid*posit_function(low_edg+(nloop-0.5)*dwid, e0)

continuum = continuum + dwid*posit_function( low_edg, e0)
continuum = continuum + dwid*posit_function( high_edg, e0)

continuum = continuum / 6.0d0


indx1 = where( wid gt 0.0d0, nindx1)
indx2 = where( wid le 0.0d0, nindx2)
if( nindx1 gt 0) then continuum(indx1) = continuum(indx1) / wid(indx1)
if( nindx2 gt 0) then continuum(indx2) = 0.0

return, float( continuum)
end



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function positf, edg, wid, epeak, fwhm, lamp, pfrac
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;       Positronium Model using positronium FRACTION
;
;       Form:   F(E) = P(E) + G(E) where P(E) is the positronium continuum
;               component and G(E) is the integrated gaussian line component.
;
;       Reference:  Ore and Powell, Phys. Rev. 75, p. 1696 (1949)
;
;       Parameters:     edg     energy edge array
;                       wid     energy width array
;                       epeak   annihilation line energy
;                       fwhm    full-width, half max
;                       LAmp    total line flux
;                       PFrac   positronium FRACTION
;
;  Modifications:
;
;       13-Jan-93 (wrp) - Corrected error in normalization of positronium
;                         continuum and added use of double precision when
;                         evaluating continuum function.
;
;       29-Nov-93 (wrp) - Modified to call POSIT_CONTINUUM which does numerical
;                         integration of positronium continuum function over
;                         the specified channels.
;
;
e0 = double( epeak)
gausscmp = new_gauss0( edg, edg+wid, 1.d0, fwhm, e0) / wid
positcmp = posit_continuum( edg, edg+wid, e0)

PAmp = Calc_PAmp(LAmp,PFrac,dPAmp_dLAmp=dPAmp_dLAmp,dPAmp_dPFrac=dPAmp_dPFrac)


return, float( (LAmp*gausscmp) + (PAmp * positcmp))
end


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function dpositf_dlamp, edg, wid, epeak, fwhm, lamp, pFrac
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;       Positronium Model
;
;       Form:   F(E) = P(E) + G(E) where P(E) is the positronium continuum
;               component and G(E) is the integrated gaussian line component.
;
;       Reference:  Ore and Powell, Phys. Rev. 75, p. 1696 (1949)
;
;       Parameters:     edg     energy edge array
;                       epeak   annihilation line energy
;                       fwhm    full-width, half max
;                       LAmp    total line flux
;                       PFrac   positronium FRACTION
;
e0 = double( epeak)
gausscmp = new_gauss0( edg, edg+wid, 1.0d0, fwhm, e0) / wid

positcmp = posit_continuum( edg, edg+wid, e0)

PAmp = Calc_PAmp(LAmp,PFrac,dPAmp_dLAmp=dPAmp_dLAmp,dPAmp_dPFrac=dPAmp_dPFrac)


return, float( gausscmp + dPAmp_dLAmp*positcmp)
end


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function dpositf_dfwhm, edg, wid, epeak, fwhm, lamp, pfrac
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;       Positronium Model - Derivative wrt FWHM
;
;       Form:   F(E) = P(E) + G(E) where P(E) is the positronium continuum
;               component and G(E) is the integrated gaussian line component.
;
;       Reference:  Ore and Powell, Phys. Rev. 75, p. 1696 (1949)
;
;       Parameters:     edg     energy edge array
;                       epeak   annihilation line energy
;                       fwhm    full-width, half max
;                       LAmp    total line flux
;                       PAmp    total positronium flux
;
;       29-Nov-93 (wrp) - New routine.  Returns analytic derivative of the
;                         positronium model with respect to the FWHM parameter.
;
e0 = double( epeak)
dgauss_dfwhm = partial_ng0_fwhm( edg, edg+wid, 1.0d0, fwhm, e0) / wid

return, float( LAmp*dgauss_dfwhm)
end


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function dpositf_dpfrac, edg, wid, epeak, fwhm, lamp, pfrac
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;       Positronium Model
;
;       Form:   F(E) = P(E) + G(E) where P(E) is the positronium continuum
;               component and G(E) is the integrated gaussian line component.
;
;       Reference:  Ore and Powell, Phys. Rev. 75, p. 1696 (1949)
;
;       Parameters:     edg     energy edge array
;                       wid     energy width array
;                       epeak   annihilation line energy
;                       fwhm    full-width, half max
;                       LAmp    total line flux
;                       PAmp    total positronium flux
;
;  Modifications:
;
;       13-Jan-93 (wrp) - Corrected error in normalization of positronium
;                         continuum and added use of double precision when
;                         evaluating continuum function.
;
;       29-Nov-93 (wrp) - Modified to call POSIT_CONTINUUM which does numerical
;                         integration of positronium continuum function over
;                         the specified channels.
;
;
e0 = double( epeak)
positcmp = posit_continuum( edg, edg+wid, e0)

PAmp = Calc_PAmp(LAmp,PFrac,dPAmp_dLAmp=dPAmp_dLAmp,dPAmp_dPFrac=dPAmp_dPFrac)


return, float( dPAmp_dPFrac*positcmp)
end


PRO POSIT_WFRAC, X, NX, Y, NY, P, NPDP, W, NW

;
; ***********************************************************************
;+
;
;  TITLE:       POSITRONIUM model w/positronium fraction as parameter
;
;  AUTHOR:      K. McNaron-Brown
;               Space Science Division
;               Naval Research Laboratory
;               Washington DC
;
;  DATE:        Sept. 2, 1992
;
;  PROJECT:     GRO/OSSE
;
;  PURPOSE:     Produces a photon spectrum model based on Positronium function.
;
;  CATEGORY:    OSSE analysis utility
;
;  CALLING SEQUENCE:
;
;       POSIT_WFRAC, X, NX, F, NY, P, NPDP, W, NW
;
;  INPUTS:
;
;       X       - 2-D Fltarr containing vector of lower count edges in
;                 first dimension, and vector of count edge x(1,*)ths in
;                 second.
;       NX      - Integer of N_ELEMENTS(X).
;                 (**Do not need to evaluate.  This is a SUPERFIT variable
;                       that is not used in this procedure.**)
;       NY      - Number of count model spectrums.
;                 (**Do not need.  This is a SUPERFIT variable
;                       that is not used in this procedure.**)
;       P       - Fltarr containing model parameter values.
;       NPDP    - Two element vector with the following definitions:
;                       NPDP(0) = Number of parameters in fit.
;                       NPDP(1) = 0 to evaluate the model itself,
;                                 i to compute the derivative of P(i).
;       W       - Fltarr containing count model uncertainties.
;                 (**Do not need to evaluate.  This is a SUPERFIT variable
;                       that is not used in this procedure.**)
;       NW      - Indicator of weighting technique.
;                 (**Do not need to evaluate.  This is a SUPERFIT variable
;                       that is not used in this procedure.**)
;
;  OUTPUTS:
;
;       Y       - Fltarr containing count model.
;
;  MODULES CALLED BY THIS MODULE:       none
;
;  MODULES THAT CALL THIS MODULE:       CALL_PROCEDURE, User programs
;
;  COMMENTS:
;
;       4-Dec-96 (wrp) - Module created from original POSITRONIUM routine.
;                        Note that this routine has the Positronium FRACTION
;                        as a parameter rather than the positronium amplitude.
;
;  MODIFICATION HISTORY:
; 		14-June-2001, Paul BIlodeau - fixed calls to positf such that
; 		width, not upper edge limit, is passed.
;
;-
; ***********************************************************************
;

        nd = npdp(1)
        if nd lt 0 then nd = 0
        a0 = p(0)                                       ; Line Amp
        a1 = p(1)                                       ; Posit FRACTION
        a2 = p(2)                                       ; FWHM
        a3 = p(3)                                       ; Line Energy


		wid = x(1,*)-x(0,*)

        case nd of
        0:      y = positf( x(0,*), wid, a3, a2, a0, a1 )

        1:      y = dpositf_dlamp( x(0,*), wid, a3, a2, a0, a1 )

        2:      y = dpositf_dpfrac( x(0,*), wid, a3, a2, a0, a1 )

        3:      y = dpositf_dfwhm( x(0,*), wid, a3, a2, a0, a1 )

        4:      begin   ;; numeric derivative
                  print,'Positronium w/Fraction Model dYdP - error, numeric derivatives required'
                  stop
                end

        else:
        endcase

return
end