;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 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 ;; Paul Bilodeau, 27-nov-2002 ;; Input is sigma - convert to FWHM a2 = p(2) * 2. * Sqrt( 2. * Alog(2) ) 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