text:
	PROGRAM GDSTAT
C************************************************************************
C* GDSTAT								*
C*									*
C* This program calculates an average, standard deviation and number	*
C* of points for grids of different times.				*
C*									*
C* Log:									*
C* I. Graffman/RDS	 4/87						*
C* M. desJardins/GSFC	11/88	GEMPAK 4.1				*
C* M. desJardins/GSFC	 1/89	Added call to IP_EXIT			*
C* M. desJardins/GSFC	11/89	Changed GR_FILE to DG_OFIL		*
C* M. desJardins/GSFC	 8/90	GEMPAK 5.0				*
C* M. desJardins/GSFC	10/90	Fix initialization of internal arrays	*
C* S. Jacobs/EAI         2/94   Added COLADD flag to DG_OFIL            *
C* L. Williams/EAI	 7/94	Removed call to GDYUPD			*
C* P. Bruehl/COMET	 8/95	Check for missing data in averaging (DK)*
C* K. Tyle/GSC		 8/96	Added FL_MFIL to search for file type	*
C* S. Maxwell/GSC        7/97   Increased input character length        *
C* T. Lee/GSC		 1/99	Increased LLMXTM to LLMXGT		*
C* S. Jacobs/NCEP	 3/01	Replaced DG_OFIL with DG_MFIL 		*
C* T. Lee/GSC		 7/01	Processed multiple files		*
C* T. Lee/SAIC		10/01	Reopened output file after DG_CLAL call	*
C* S. Jacobs/NCEP	 1/02	Fixed DG_MFIL to not add a col to grid	*
C* K. Brill/HPC		11/02	Added MAX and MIN grids			*
C* K. Brill/HPC		 4/03	CALL DG_INTL				*
C* K. Brill/HPC		 5/03	Check for change in grid size (kx,ky)	*
C* R. Tian/SAIC         11/03   Added nuflg to DG_INTL call             *
C* R. Tian/SAIC          2/04   Removed nuflg from DG_INTL call         *
C* A. Hardy/NCEP	 4/04	Removed check on ntime > MXLOOP		*
C* R. Tian/SAIC          3/05   Changes for time/file mngmnt            *
C* T. Piper/SAIC        01/08   Added GD_INIT; removed from IN_BDTA     *
C* B. Hebbard/SAIC      03/08   Remove LLMXGD size restriction on grids *
C************************************************************************
	INCLUDE 	'GEMPRM.PRM'
C*
	CHARACTER	gdfile*(LLMXLN), gfunc*(LLMXLN),
     +			gdatim*(LLMXLN), glevel*(LLMXLN),
     +			gvcord*(LLMXLN), gdoutf*(LLMXLN),
     +			grdnam*(LLMXLN)
        CHARACTER       dattim (2)*20
C*
	LOGICAL		respnd, done, proces, gottm
C*
	INCLUDE		'ERMISS.FNC'
C------------------------------------------------------------------------
C*  Initialize TAE.
C
	CALL IP_INIT  ( respnd, iperr )
	IF  ( iperr .eq. 0 )  THEN
	    CALL IP_IDNT  ( 'GDSTAT', ier )
C
C*  Initialize GEMPLT.
C
	    mode = 1
	    CALL GG_INIT  ( mode, ier )
	    IF  ( ier .eq. 0 )  THEN
C
C*  Initialize grid library common area grdcmn.cmn.
C
		CALL GD_INIT  ( ier )
C
C*  Initialize the DG library.
C
		CALL DG_INTL ( ier )
		done = .false.
	    ELSE
		done = .true.
	    END IF
	ELSE
	    done = .true.
	END IF
C
C*	Main loop for user input and computations.
C
	DO WHILE  ( .not. done )
C
C*	  Set flag to indicate processing will be done.
C
	  proces = .true.
	  CALL GDYINP  ( gdfile, gdoutf, gfunc, gdatim, glevel,
     +			 gvcord, grdnam, iperr )
C
C*	  Exit if there is an error.
C
	  IF  ( iperr .ne. 0 )  THEN
	    done   = .true.
	  ELSE
C
C*          Process the GDFILE and GDOUTF input.
C
            CALL DG_NFIL ( gdfile, gdoutf, ier )
            IF ( ier .ne. 0 ) THEN
                CALL ER_WMSG ( 'DG', ier, ' ', irr )
                proces = .false.
            END IF
C
C*          Process the GDATTIM input.
C
            CALL DG_NDTM ( gdatim, ier )
            IF ( ier .ne. 0 ) THEN
                CALL ER_WMSG ( 'DG', ier, gdatim, irr )
                proces = .false.
            END IF
C
C           Get the grid dimensions:
C
C           (a)  Before we can call DG_KXKY to get grid dimensions,
C                we make a call to DG_NTIM to advance to the next
C                (in this case first) time in the grid file...
C
            IF (proces) THEN
              CALL DG_NTIM ( .false., .false., dattim, gottm, ier )
              IF ( ier .ne. 0 ) THEN
                CALL ER_WMSG ( 'DG', ier, ' ', irr )
                proces = .false.
              END IF
            END IF
C
C           (b)  Now get the grid dimensions.  Note that this program
C                (GDSTAT) will error if all grids (across all times)
C                do not have the same dimensions, so it does not matter
C                that we query just the first grid in the file...
C
            IF (proces) THEN
              CALL DG_KXKY ( kx, ky, ier )
              IF ( ier .ne. 0 ) THEN
                CALL ER_WMSG ( 'DG', ier, ' ', irr )
                proces = .false.
              END IF
            END IF
C
C           (c) However, we need to reset the grid time back to the
C               first grid in the file, so that it will be selected
C               when the real work begins later...
C
            IF (proces) THEN
              CALL DG_RSTM ( ier )
              IF ( ier .ne. 0 ) THEN
                CALL ER_WMSG ( 'DG', ier, ' ', irr )
                proces = .false.
              END IF
            END IF
C
C*          Call the wrapper routine to allocate storage and do the work
C
            IF (proces) THEN
              CALL GDSTAA ( kx, ky, glevel, gvcord, gfunc, grdnam, ier )
            END IF
C
C
C*	    Prompt for next action to be done.
C
	    CALL IP_DYNM  ( done, ier )
	  END IF
	END DO
C
C*	Print general error messages if necessary.
C
        IF ( iperr .ne. 0 ) CALL ER_WMSG ( 'GDSTAT', iperr, ' ', ier )
C
C*      Exit from the TAE.
C
	CALL IP_EXIT  ( iret )
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
program P2622
    implicit none
    integer :: length_sum
    integer :: a, b, c
    integer :: count = 0
    read *, length_sum
    do a = 1, length_sum / 3
        do b = a, length_sum
            c = length_sum - a - b
            if (b <= c .and. (a + b) > c) then
                count = count + 1
            end if
        end do
    end do
    write (*, '(I0)') count
end program P2622

max_stars_repo_name: daily-boj/kiwiyou

max_issues_repo_name: daily-boj/kiwiyou

max_forks_repo_name: daily-boj/kiwiyou


text:
!-------------------------------------------------------------------------------
! Copyright (c) 2021, Whitman T. Dailey
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
! 1. Redistributions of source code must retain the above copyright notice, this
!    list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright notice,
!    this list of conditions and the following disclaimer in the documentation
!    and/or other materials provided with the distribution.
! 3. Neither the name of the copyright holder nor the names of its
!    contributors may be used to endorse or promote products derived from
!    this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
! FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!-------------------------------------------------------------------------------
Module Quadratures

    Use Kinds, Only: dp
    Implicit None
    Private
    Public :: Romberg_Quad
    Public :: Composite_Quad
    Public :: Composite_Trapezoid
    Public :: Composite_Simpson
    Public :: Adaptive_Quad
    Public :: Adaptive_Trapezoid
    Public :: Adaptive_Simpson
    Public :: MidPoint_Quad
    Public :: Trapezoid_Quad
    Public :: Simpson_Quad
    Public :: GaussLegendre5
    Public :: GaussLegendre16
    Public :: GaussLegendre96
    Public :: GaussLegendreN
    Public :: Progressive_GaussLegendre

    Real(dp), Parameter :: one_third = 1._dp / 3._dp
    Real(dp), Parameter :: one_sixth = 1._dp / 6._dp
    Integer, Parameter :: bisect_limit = 53

    Interface Romberg_Quad
        Module Procedure Romb_Quad
        Module Procedure Romb_Quad_ranges
    End Interface

Contains

Function Romb_Quad_ranges(f,ab,aTol,rTol,p) Result(q)
    !Integrates f(x) on the intervals between entries in ab(:) by extrapolation on successive composite trapezoid
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: q  !the result of the integration
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds,Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: ab(:)  !list of limits of integration
    Real(dp), Intent(In) :: aTol,rTol  !absolute and relative tolerances for convergence
    Real(dp), Intent(Out), Optional :: p  !estimated precision achived in the final estimate
    Real(dp) :: h  !the total width of the integration window from ab(1) to ab(n)
    Real(dp) :: hi  !the width of the i-th integration window
    Integer :: n  !length of the list of integration limits, resulting in n-1 integrations to be performed
    Integer :: i  !loop counter
    Real(dp) :: p_i  !estimated precision of the integral from the i-th interval

    q = 0._dp
    n = Size(ab)
    h = ab(n) - ab(1)
    Do i = 1,n-1
        hi = ab(i+1) - ab(i)
        If (Present(p)) Then
            q = q + Romb_Quad(f,ab(i),ab(i+1),(hi/h)*aTol,rTol,p=p_i)
            p = Min(p,p_i)
        Else
            q = q + Romb_Quad(f,ab(i),ab(i+1),(hi/h)*aTol,rTol)
        End If
    End Do
End Function Romb_Quad_ranges

Function Romb_Quad(f,a,b,aTol,rTol,allow_GL_fallback,p,n_ord,n_ext) Result(q)
    !Integrates f(x) on (a,b) by extrapolation on successive composite trapezoid
    Use Kinds, Only: dp
    Use Utilities, Only: Converged
    Use Utilities, Only: Prec
    Implicit None
    Real(dp) :: q    !the result of the integration
    Interface
        Function f(x)    !the function to be integrated
            Use Kinds,Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !limits of integration
    Real(dp), Intent(In) :: aTol,rTol  !absolute and relative tolerances for convergence
    Logical, Intent(In), Optional :: allow_GL_fallback  !flag signaling whether fallback to progressive G-L should be used in the
                                                        !event of a failure to converge during extrapolation stages
    Real(dp), Intent(Out), Optional :: p  !estimated precision achived in the final estimate
    Integer, Intent(Out), Optional :: n_ord  !total number of function evaluations (ordinates)
    Integer, Intent(Out), Optional :: n_ext  !total number of extrapolation stages (rows in the table)
    Integer, Parameter :: Tmax = 21  !maximum number of extrapolations in the table
    Real(dp) :: T(0:Tmax)  !Extrapolation table previous row
    Real(dp) :: Tk0,Tk  !Extrapolation table current row values
    Integer :: i,j,k  !counters: i for table row, j for quadrature ordinates, k for table column
    Integer :: n      !number of intervals
    Real(dp) :: h0,h  !spacing between quadrature ordinates
    Real(dp) :: fk    !multiplier for extrapolation steps
    Real(dp) :: s     !sum of function values at quadrature ordinates
    Real(dp) :: Prec0 !precision of the previous extrapolation (used to monitor convergence)
    Integer :: ifin   !final number of extrapolations performed
    Logical :: GL_fallback
    !    If the preprocessor directive ROMB_TABLES is set at compile-time,
    ! then the tables of computed values are dumped to file in the working
    ! directory as the routine runs.
#   if ROMB_TABLES
        Integer :: unit
#   endif

    If (a .EQ. b) Then !zero-width interval for integration
        q = 0._dp
        If (Present(p)) p = Prec(0._dp,q)
        If (Present(n_ord)) n_ord = 0
        If (Present(n_ext)) n_ext = 0
        Return
    End If
    !Initial trapezoid estimate
    n = 1
    s = 0.5_dp * (f(a) + f(b))
    h0 = b - a
    T(0) = h0 * s
    Prec0 = -1._dp
#   if ROMB_TABLES
        Open(NEWUNIT=unit,FILE='Romb_Tables.tst',ACTION='WRITE',STATUS='UNKNOWN',POSITION='APPEND')
        Do k = 0,Tmax
            Write(unit,'(I24)',ADVANCE='NO') k
        End Do
        Write(unit,'(ES24.15)') T(0)
        Close(unit)
#   endif
    Do i = 1,Tmax !up to Tmax rows in the table
        !Trapezoid estimate for the 0-th column of the i-th row of table
        n = n * 2
        h = h0 / Real(n,dp)
        Do j = 1,n-1,2  !Odd values of j are NEW points at which to evaluate f
            s = s + f(a + Real(j,dp)*h)
        End Do
        Tk0 = h * s
        !Fill i-th row, columns k = 1:i, with extrapolated estimates
        fk = 1._dp
        Do k = 1,i  !up to i columns this row
            fk = fk * 4._dp
            Tk = (fk * Tk0 - T(k-1)) / (fk - 1._dp)
            If (k .LT. i) Then
                T(k-1) = Tk0  !store Tk0 for next i
                Tk0 = Tk  !store Tk for next k
            End If !otherwise, skip storage steps if working final column
        End Do
        !Check for convergence
        If (Converged(T(i-1),Tk,rTol=Real(2**k-1)*rTol,aTol=aTol)) Then
            q = Tk
#           if ROMB_TABLES
                T(i-1) = Tk0
                T(i) = Tk
                Open(NEWUNIT=unit,FILE='Romb_Tables.tst',ACTION='WRITE',STATUS='UNKNOWN',POSITION='APPEND')
                Do k = 0,i
                    Write(unit,'(ES24.15)',ADVANCE='NO') T(k)
                End Do
                Write(unit,'(A)') '*'
                Write(unit,*)
                Close(unit)
#           endif
            If (Present(p)) p = Prec(T(i-1),Tk)
            If (Present(n_ord)) n_ord = n
            If (Present(n_ext)) n_ext = i
            Return  !Normal exit
        Else  !Check for failures and prep for the next time though the loop
            !check for exit conditions other than convergence (failures)
            If (i .EQ. Tmax) Then  !extrapolation limit has been reached without convergence
                ifin = Tmax
                Exit
            Else If (i .GT. Tmax/2) Then
                If (Prec(T(i-1),Tk) .LT. Prec0) Then  !precision was LOST on this extrapolation, convergence will not occur
                    ifin = i
                    Exit
                End If
            End If
            !store Tk0 and Tk for next i, and update precision monitor
            Prec0 = Prec(T(i-1),Tk)
            T(i-1) = Tk0
            T(i) = Tk
#           if ROMB_TABLES
                Open(NEWUNIT=unit,FILE='Romb_Tables.tst',ACTION='WRITE',STATUS='UNKNOWN',POSITION='APPEND')
                Do k = 0,i
                    Write(unit,'(ES24.15)',ADVANCE='NO') T(k)
                End Do
                Write(unit,*)
                Close(unit)
#           endif
        End If
    End Do
    !If we get this far, we did not converge
    If (Present(allow_GL_fallback)) Then
        GL_fallback = allow_GL_fallback
    Else
        GL_fallback = .TRUE.
    End If
    If (GL_fallback) Then
        If (Present(p)) Then
            q = Progressive_GaussLegendre(f,a,b,aTol=aTol,rTol=rTol,n_start=5,n_stride=3,n_done=n,p=p)
        Else
            q = Progressive_GaussLegendre(f,a,b,aTol=aTol,rTol=rTol,n_start=5,n_stride=3,n_done=n)
        End If
        If (Present(n_ord)) n_ord = n
        If (Present(n_ext)) n_ext = -1
    Else
        Write(*,*)
        Write(*,'(A,I0,A)')          'WARNING:  Quadratures: Romberg_Quad:  Failed to converge in ',ifin,' extrapolations.'
        If (ifin .LT. Tmax) Then
            Write(*,'(A10,A)')              '','Extrapolation was terminated for loss of precision.'
        End If
        Write(*,'(A10,A,ES23.15E2,A,F0.5,A)') '','Final estimated value: ',Tk,' (~',Prec(T(ifin-1),Tk),' good digits)'
        Write(*,'(A10,A,2(ES10.3E2,A))')      '','Final Extrapolation Error: ',Abs(Tk-T(ifin-1)),   ' (abs), ', &
                                       &                                       Abs(Tk-T(ifin-1))/Tk,' (rel)'
        Write(*,'(A10,A,2(ES10.3E2,A))')      '','Convergence Criteria:      ',atol,                ' (abs), ', &
                                       &                                       rtol,                ' (rel)'
        q = Tk
#       if ROMB_TABLES
            T(i-1) = Tk0
            T(i) = Tk
            Open(NEWUNIT=unit,FILE='Romb_Tables.tst',ACTION='WRITE',STATUS='UNKNOWN',POSITION='APPEND')
            Do k = 0,i
                Write(unit,'(ES24.15)',ADVANCE='NO') T(k)
            End Do
            Write(unit,'(A)') '*'
            Write(unit,*)
            Close(unit)
#       endif
        If (Present(p)) p = Prec(T(ifin-1),Tk)
        If (Present(n_ord)) n_ord = n
        If (Present(n_ext)) n_ext = ifin
    End If
    ! ERROR STOP
    !Currently, no STOP is implemented for a failure to converge.  The result of the quadrature is set to the final value from the
    !Romberg table, and additional required values are set, and execution returns to the calling routine after the above warnings
    !are sent to the screen.
End Function Romb_Quad

Function Richardson_Extrapolation(fh,f2h,p) Result(r)
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: r
    Real(dp), Intent(In) :: fh,f2h
    Integer, Intent(In) :: p

    r = ( fh*Real(2**p,dp) - f2h ) / Real(2**p-1,dp)
End Function Richardson_Extrapolation

Function Adaptive_Quad(f,a0,b0,Quad,aTol,rTol) Result(q)
    !Integrates f(x) on (a,b) by adaptively applying composite quadrature
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: q  !the result of the integration
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds,Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
        Function Quad(f,a,b,p)  !the quadrature rule to use
            Use Kinds,Only: dp
            Implicit None
            Real(dp) :: Quad  !the result of the quadrature
            Interface
                Function f(x)    !the function to be integrated
                    Use Kinds,Only: dp
                    Implicit None
                    Real(dp) :: f
                    Real(dp), Intent(In) :: x
                End Function f
            End Interface
            Real(dp),Intent(In) :: a,b  !limits of integration
            Integer, Intent(Out), Optional :: p
        End Function Quad
    End Interface
    Real(dp), Intent(In) :: a0,b0  !limits of integration
    Real(dp), Intent(In) :: aTol,rTol  !absolute and relative tolerance to be used as criteria for convergence
    Integer :: n,i
    Real(dp), Allocatable :: s(:,:)
    integer, Allocatable :: sb(:)
    Real(dp) :: a,b,ab
    Real(dp) :: q0,q1,q2
    Real(dp) :: h
    Real(dp) :: rTol_per_h,aTol_per_h
    Real(dp) :: r,r0,dr
    Integer :: p
    Integer :: bisects

    n = 16
    i = 0
    Allocate(s(1:3,1:n))
    s = 0._dp
    Allocate(sb(1:n))
    sb = 0
    q0 = Quad(f,a0,b0,p=p)
    Call Put_s(a0,b0,q0,0)
    rTol_per_h = rTol * Real(2**p-1,dp)
    aTol_per_h = aTol / (b0 - a0)
    q = 0._dp
    Do
        If (Get_s(a,b,q0,bisects)) Then  !got an interval to do
            ab = 0.5_dp * (a + b)
            q1 = Quad(f,a,ab)
            q2 = Quad(f,ab,b)
            If (p .GT. 1) Then
                r = Richardson_Extrapolation(q1+q2,q0,p)
                r0 = q1 + q2
            Else
                r = q1 + q2
                r0 = q0
            End If
            h = b - a
            dr = Abs( r - r0 )
            If (      dr .LT. h * rTol_per_h    & !relative convergence
               & .OR. dr .LT. h * aTol_per_h    & !absolute convergence
               & .OR. bisects .GT. bisect_limit & !bisection limit exceeded
               &      ) Then
                q = q + r  !add to result of the quadrature
            Else  !add the subintervals to the list to be refined again
                Call Put_s(a,ab,q1,bisects+1)
                Call Put_s(ab,b,q2,bisects+1)
            End If
        Else !no intervals left to quadrature
            Exit
        End If
    End Do
Contains
    Subroutine Put_s(c1,c2,c3,bs)
        Use Kinds, Only: dp
        Implicit None
        Real(dp), Intent(In) :: c1,c2,c3
        Integer, Intent(In) :: bs
        Real(dp), Allocatable :: swaps(:,:)
        Integer, Allocatable :: swapsb(:)

        i = i + 1
        If (i .GT. n) Then !need more space
            Allocate(swaps(1:3,1:n))
            swaps = s
            Deallocate(s)
            Allocate(s(1:3,1:2*n))
            s = 0._dp
            s(:,1:n) = swaps
            Deallocate(swaps)
            Allocate(swapsb(1:n))
            swapsb = sb
            Deallocate(sb)
            Allocate(sb(1:2*n))
            sb = 0
            sb(1:n) = swapsb
            Deallocate(swapsb)
            n = n * 2
        End If
        s(:,i) = (/ c1 , c2 , c3 /)
        sb(i) = bs
    End Subroutine Put_s
    Function Get_s(c1,c2,c3,bs) Result(bingo)
        Use Kinds, Only: dp
        Implicit None
        Logical :: bingo
        Real(dp), Intent(Out) :: c1,c2,c3
        Integer, Intent(Out) :: bs

        If (i .LT. 1) Then !there was no s to get
            bingo = .FALSE.
        Else
            bingo = .TRUE.
            c1 = s(1,i)
            c2 = s(2,i)
            c3 = s(3,i)
            bs = sb(i)
            i = i - 1
        End If
    End Function Get_s
End Function Adaptive_Quad

Function Adaptive_Trapezoid(f,a0,b0,aTol,rTol) Result(q)
    !Integrates f(x) on (a,b) by adaptively applying composite trapezoid
    Use Kinds,Only: dp
    Implicit None
    Real(dp) :: q  !the result of the integration
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds,Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a0,b0  !limits of integration
    Real(dp), Intent(In) :: aTol,rTol  !absolute and relative tolerance to be used as criteria for convergence
    Integer :: n,i
    Real(dp), Allocatable :: s(:,:)
    Integer, Allocatable :: sb(:)
    Real(dp) :: a,b,ab
    Real(dp) :: fa,fb,fab
    Real(dp) :: q0,q1,q2,r,dr
    Real(dp) :: h
    Real(dp) :: rTol_per_h,aTol_per_h
    Integer :: bisects

    n = 16
    i = 0
    Allocate(s(1:5,1:n))
    s = 0._dp
    Allocate(sb(1:n))
    sb = 0
    i = 0
    !Initial Trapezoid estimate
    fa = f(a0)
    fb = f(b0)
    q0 = (b0 - a0) * (fa + fb) * 0.5_dp
    Call Put_s(a0,b0,fa,fb,q0,0)
    rTol_per_h = rTol * 3._dp
    aTol_per_h = aTol / (b0 - a0)
    q = 0._dp
    Do
        If (Get_s(a,b,fa,fb,q0,bisects)) Then  !got an interval to do
            ab = 0.5_dp * (a + b)
            fab = f(ab)
            h = b - a
            q1 = 0.25_dp * h * (fa + fab)
            q2 = 0.25_dp * h * (fab + fb)
            r = Richardson_Extrapolation(q1+q2,q0,2)
            dr = Abs( q1 + q2 - r )
            If (      dr .LT. h * rTol_per_h    & !relative convergence
               & .OR. dr .LT. h * aTol_per_h    & !absolute convergence
               & .OR. bisects .GT. bisect_limit & !bisection limit exceeded
               &      ) Then
                q = q + r
            Else  !add the subintervals to the list to be quadratured
                Call Put_s(a,ab,fa,fab,q1,bisects+1)
                Call Put_s(ab,b,fab,fb,q2,bisects+1)
            End If
        Else !no intervals left to quadrature
            Exit
        End If
    End Do
Contains
      Subroutine Put_s(c1,c2,c3,c4,c5,bs)
        Use Kinds, Only: dp
        Implicit None
        Real(dp), Intent(In) :: c1,c2,c3,c4,c5
        Integer, Intent(In) :: bs
        Real(dp), Allocatable :: swaps(:,:)
        Integer, Allocatable :: swapsb(:)

        i = i + 1
        If (i .GT. n) Then !need more space
            Allocate(swaps(1:5,1:n))
            swaps = s
            Deallocate(s)
            Allocate(s(1:5,1:2*n))
            s = 0._dp
            s(:,1:n) = swaps
            Deallocate(swaps)
            Allocate(swapsb(1:n))
            swapsb = sb
            Deallocate(sb)
            Allocate(sb(1:2*n))
            sb = 0
            sb(1:n) = swapsb
            Deallocate(swapsb)
            n = n * 2
        End If
        s(:,i) = (/ c1 , c2 , c3 , c4, c5 /)
        sb(i) = bs
    End Subroutine Put_s
    Function Get_s(c1,c2,c3,c4,c5,bs) Result(bingo)
        Use Kinds, Only: dp
        Implicit None
        Logical :: bingo
        Real(dp), Intent(Out) :: c1,c2,c3,c4,c5
        Integer, Intent(Out) :: bs

        If (i .LT. 1) Then !there was no s to get
            bingo = .FALSE.
        Else
            bingo = .TRUE.
            c1 = s(1,i)
            c2 = s(2,i)
            c3 = s(3,i)
            c4 = s(4,i)
            c5 = s(5,i)
            bs = sb(i)
            i = i - 1
        End If
    End Function Get_s
End Function Adaptive_Trapezoid

Function Adaptive_Simpson(f,a0,b0,aTol,rTol) Result(q)
    !Integrates f(x) on (a,b) by adaptively applying composite Simpson
    Use Kinds,Only: dp
    Implicit None
    Real(dp) :: q  !the result of the integration
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds,Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a0,b0  !limits of integration
    Real(dp), Intent(In) :: aTol,rTol  !absolute and relative tolerance to be used as criteria for convergence
    Integer :: n,i
    Real(dp), Allocatable :: s(:,:)
    Integer, Allocatable :: sb(:)
    Real(dp) :: x1,x2,x3,x4,x5
    Real(dp) :: y1,y2,y3,y4,y5
    Real(dp) :: q0,q1,q2,r,dr
    Real(dp) :: h
    Real(dp) :: rTol_per_h,aTol_per_h
    Integer :: bisects

    n = 16
    i = 0
    Allocate(s(1:7,1:n))
    s = 0._dp
    Allocate(sb(1:n))
    sb = 0
    i = 0
    !Initial Simpson estimate
    x3 = 0.5_dp * (a0 + b0)
    y1 = f(a0)
    y3 = f(x3)
    y5 = f(b0)
    q0 = (b0 - a0) * (y1 + 4._dp*y3 + y5) * one_sixth
    Call Put_s(a0,x3,b0,y1,y3,y5,q0,0)
    rTol_per_h = rTol * 15._dp
    aTol_per_h = aTol / (b0 - a0)
    q = 0._dp
    Do
        If (Get_s(x1,x3,x5,y1,y3,y5,q0,bisects)) Then  !got an interval to do
            x2 = 0.5_dp * (x1 + x3)
            x4 = 0.5_dp * (x3 + x5)
            y2 = f(x2)
            y4 = f(x4)
            q1 = (x3-x1) * (y1 + 4._dp*y2 + y3) * one_sixth
            q2 = (x5-x3) * (y3 + 4._dp*y4 + y5) * one_sixth
            r = Richardson_Extrapolation(q1+q2,q0,4)
            h = x5 - x1
            dr = Abs( q1 + q2 - r )
            If (      dr .LT. h * rTol_per_h    & !relative convergence
               & .OR. dr .LT. h * aTol_per_h    & !absolute convergence
               & .OR. bisects .GT. bisect_limit & !bisection limit exceeded
               &      ) Then
                q = q + r
            Else  !add the subintervals to the list to be quadratured
                Call Put_s(x1,x2,x3,y1,y2,y3,q1,bisects+1)
                Call Put_s(x3,x4,x5,y3,y4,y5,q2,bisects+1)
            End If
        Else !no intervals left to quadrature
            Exit
        End If
    End Do
Contains
      Subroutine Put_s(c1,c2,c3,c4,c5,c6,c7,bs)
        Use Kinds, Only: dp
        Implicit None
        Real(dp), Intent(In) :: c1,c2,c3,c4,c5,c6,c7
        Integer, Intent(In) :: bs
        Real(dp), Allocatable :: swaps(:,:)
        Integer, Allocatable :: swapsb(:)

        i = i + 1
        If (i .GT. n) Then !need more space
            Allocate(swaps(1:7,1:n))
            swaps = s
            Deallocate(s)
            Allocate(s(1:7,1:2*n))
            s = 0._dp
            s(:,1:n) = swaps
            Deallocate(swaps)
            Allocate(swapsb(1:n))
            swapsb = sb
            Deallocate(sb)
            Allocate(sb(1:2*n))
            sb = 0
            sb(1:n) = swapsb
            Deallocate(swapsb)
            n = n * 2
        End If
        s(:,i) = (/ c1 , c2 , c3 , c4, c5 , c6 , c7 /)
        sb(i) = bs
    End Subroutine Put_s
    Function Get_s(c1,c2,c3,c4,c5,c6,c7,bs) Result(bingo)
        Use Kinds, Only: dp
        Implicit None
        Logical :: bingo
        Real(dp), Intent(Out) :: c1,c2,c3,c4,c5,c6,c7
        Integer, Intent(Out) :: bs

        If (i .LT. 1) Then !there was no s to get
            bingo = .FALSE.
        Else
            bingo = .TRUE.
            c1 = s(1,i)
            c2 = s(2,i)
            c3 = s(3,i)
            c4 = s(4,i)
            c5 = s(5,i)
            c6 = s(6,i)
            c7 = s(7,i)
            bs = sb(i)
            i = i - 1
        End If
    End Function Get_s
End Function Adaptive_Simpson

Function Composite_Quad(f,a,b,Quad,aTol,rTol) Result(q)
    !Integrates f(x) on (a,b) by composite quadrature
    Use Kinds, Only: dp
    Use Utilities, Only: Converged
    Implicit None
    Real(dp) :: q    !the result of the integration
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
        Function Quad(f,a,b,p)  !the quadrature rule to use
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: Quad  !the result of the quadrature
            Interface
                Function f(x)    !the function to be integrated
                    Use Kinds ,Only: dp
                    Implicit None
                    Real(dp) :: f
                    Real(dp), Intent(In) :: x
                End Function f
            End Interface
            Real(dp), Intent(In) :: a,b  !interval on which to apply the rule
            Integer, Intent(Out), Optional :: p
        End Function Quad
    End Interface
    Real(dp), Intent(In) :: a,b  !limits of integration
    Real(dp), Intent(In) :: aTol,rTol  !absolute and relative tolerance to be used as criteria for convergence
    Integer :: i,j,n
    Real(dp) :: h,h0
    Real(dp) :: q0,r,dr
    Real(dp) :: eps
    Integer :: p

    n = 1
    q0 = Quad(f,a,b,p=p)
    h0 = b - a
    eps = Real(2**p-1,dp) * h0 * rTol
    Do i = 1,bisect_limit+1  !max intervals is bounded by splitting interval bisect_limit+1 times
        q = 0._dp
        n = n * 2
        h = h0 / Real(n,dp)
        Do j = 1,n
            q = q + Quad(f , a+h*Real((j-1),dp) , a+h*Real(j,dp))
        End Do
        If (p .GT. 1) Then !apply a stage of Richardson extrapolation for a low-cost refinement of the estimate
            r = Richardson_Extrapolation(q,q0,p)
            dr = Abs( r - q )
            If (      dr .LT. eps  & !relative convergence
               & .OR. dr .LT. aTol & !absolute convergence
               &      ) Then
                q = r
                Return
            End If
        Else !Richardson extrapolation does not provide a benefit for methods with p<2 (i.e. Midpoint, etc.)
            If (Converged(q0,q,rTol=rTol,aTol=aTol)) Return
        End If
        q0 = q
    End Do
    !if we get this far, we failed to converge
    Write(*,'(A,I0,A)') 'ERROR:  Quadratures: Composite_Quad:  Failed to converge on ',n,' intervals.'
    ERROR STOP
End Function Composite_Quad

Function MidPoint_Quad(f,a,b,p) Result(q)
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: q  !the result of the quadrature
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !interval on which to apply the rule
    Integer, Intent(Out), Optional :: p

    q = (b - a) * f( 0.5_dp * (a + b) )
    If (Present(p)) p = 2
End Function MidPoint_Quad

Function Trapezoid_Quad(f,a,b,p) Result(q)
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: q  !the result of the quadrature
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !interval on which to apply the rule
    Integer, Intent(Out), Optional :: p

    q = 0.5_dp * (b - a) * ( f(a) + f(b) )
    If (Present(p)) p = 2
End Function Trapezoid_Quad

Function Simpson_Quad(f,a,b,p) Result(q)
    Use Kinds, Only: dp
    Implicit None
    Real(dp):: q  !the result of the quadrature
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !interval on which to apply the rule
    Integer, Intent(Out), Optional :: p

    q = one_sixth * (b - a) * ( f(a) + 4._dp * f( 0.5_dp * (a + b) ) + f(b) )
    If (Present(p)) p = 4
End Function Simpson_Quad

Function Composite_Trapezoid(f,a,b,aTol,rTol) Result(q)
    !Integrates f(x) on (a,b) by composite trapezoid
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: q  !the result of the integration
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !limits of integration
    Real(dp), Intent(In) :: aTol,rTol  !absolute and relative tolerance to be used as criteria for convergence
    Integer :: i,j,n
    Real(dp) :: h,h0
    Real(dp) :: q0,s
    Real(dp) :: r,dr
    Real(dp) :: eps

    n = 1
    s = 0.5_dp * ( f(a) + f(b) )
    h0 = b - a
    eps = 3._dp * h0 * rTol
    q0 = h0 * s
    Do i = 1,bisect_limit+1  !max intervals is bounded by splitting interval bisect_limit+1 times
        n = n * 2
        h = h0 / Real(n,dp)
        Do j = 1,n-1,2  !only odd values of j, these are the NEW points at which to evaluate f
            s = s + f( a + h*Real(j,dp) )
        End Do
        q = h * s
        r = Richardson_Extrapolation(q,q0,2)
        dr = Abs( r - q )
        If (      dr .LT. eps  & !relative convergence
           & .OR. dr .LT. aTol & !absolute convergence
           &      ) Then
             q = r
            Return
        End If
        q0 = q
    End Do
    !if we get this far, we failed to converge
    Write(*,'(A,I0,A)') 'ERROR:  Quadratures: Progressive_Trapezoid:  Failed to converge on ',n,' intervals.'
    ERROR STOP
End Function Composite_Trapezoid

Function Composite_Simpson(f,a,b,aTol,rTol) Result(q)
    !Integrates f(x) on (a,b) by composite Simpson
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: q  !the result of the integration
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !limits of integration
    Real(dp), Intent(In) :: aTol,rTol  !absolute and relative tolerance to be used as criteria for convergence
    Integer :: i,j,n
    Real(dp) :: h,h0
    Real(dp) :: q0,s1,s2,s3
    Real(dp) :: r,dr,eps

    n = 2
    s1 = f(a) + f(b)
    s2 = 0._dp
    s3 = f( 0.5_dp * (a + b) )
    h0 = b - a
    eps = 15._dp * h0 * rTol
    q0 = one_sixth * h0 * (s1 + 4._dp*s3)
    Do i = 1,bisect_limit+1  !max intervals is bounded by splitting interval bisect_limit+1 times
        n = n * 2
        h = h0 / Real(n,dp)
        s2 = s2 + s3
        s3 = 0._dp
        Do j = 1,n-1,2  !only odd values of j, these are the NEW points at which to evaluate f
            s3 = s3 + f( a + h*Real(j,dp) )
        End Do
        q = one_third * h * (s1 + 2._dp*s2 + 4._dp*s3)
        r = Richardson_Extrapolation(q,q0,4)
        dr = Abs( r - q )
        If (      dr .LT. eps  & !relative convergence
           & .OR. dr .LT. aTol & !absolute convergence
           &      ) Then
             q = r
            Return
        End If
        q0 = q
    End Do
    !if we get this far, we failed to converge
    Write(*,'(A)') 'ERROR:  Quadratures: Progressive_Simpson:  Failed to converge on ',n,' intervals.'
    ERROR STOP
End Function Composite_Simpson

Function GaussLegendre(f,a,b,n,wi,xi) Result(q)
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: q  !result of the quadrature
    Interface
        Function f(x)  !the function to integrate
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !interval on which to apply the rule
    Integer, Intent(In) :: n
    Real(dp), Intent(In) :: wi(1:n),xi(1:n)
    Real(dp) :: c1,c2  !changes interval from (a,b) to (-1,1)
    Integer :: i
    Real(dp) :: fi(1:n)

    c1 = 0.5_dp * (b - a)
    c2 = 0.5_dp * (b + a)
    Do i = 1,n
        fi(i) = f(c1 * xi(i) + c2)
    End Do
    q = c1 * Dot_Product(wi,fi)
End Function GaussLegendre

Function GaussLegendre5(f,a,b) Result(q)
    Use Kinds, Only: dp
    Use GaussLegendre_w_a, Only: n => n005
    Use GaussLegendre_w_a, Only: wi => w005
    Use GaussLegendre_w_a, Only: xi => a005
    Implicit None
    Real(dp) :: q  !result of the quadrature
    Interface
        Function f(x)  !the function to integrate
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !interval on which to apply the rule

    q = GaussLegendre(f,a,b,n,wi,xi)
End Function GaussLegendre5

Function GaussLegendre16(f,a,b) Result(q)
    Use Kinds, Only: dp
    Use GaussLegendre_w_a, Only: n => n016
    Use GaussLegendre_w_a, Only: wi => w016
    Use GaussLegendre_w_a, Only: xi => a016
    Implicit None
    Real(dp) :: q  !result of the quadrature
    Interface
        Function f(x)  !the function to integrate
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !interval on which to apply the rule

    q = GaussLegendre(f,a,b,n,wi,xi)
End Function GaussLegendre16

Function GaussLegendre96(f,a,b) Result(q)
    Use Kinds, Only: dp
    Use GaussLegendre_w_a, Only: n => n096
    Use GaussLegendre_w_a, Only: wi => w096
    Use GaussLegendre_w_a, Only: xi => a096
    Implicit None
    Real(dp) :: q  !result of the quadrature
    Interface
        Function f(x)  !the function to integrate
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !interval on which to apply the rule

    q = GaussLegendre(f,a,b,n,wi,xi)
End Function GaussLegendre96

Function GaussLegendreN(n,f,a,b,p) Result(q)
    Use Kinds, Only: dp
    Use GaussLegendre_w_a
    Implicit None
    Real(dp) :: q  !result of the quadrature
    Integer, Intent(In) :: n  !number of quadrature points
    Interface
        Function f(x)  !the function to integrate
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !interval on which to apply the rule
    Integer, Intent(Out), Optional :: p
    Real(dp), Allocatable :: wi(:),ai(:)  !weights and abscissa

    Allocate(wi(1:n))
    Allocate(ai(1:n))
    Select Case (N)
        Case (1)
            wi = w001
            ai = a001
        Case (2)
            wi = w002
            ai = a002
        Case (3)
            wi = w003
            ai = a003
        Case (4)
            wi = w004
            ai = a004
        Case (5)
            wi = w005
            ai = a005
        Case (6)
            wi = w006
            ai = a006
        Case (7)
            wi = w007
            ai = a007
        Case (8)
            wi = w008
            ai = a008
        Case (9)
            wi = w009
            ai = a009
        Case (10)
            wi = w010
            ai = a010
        Case (11)
            wi = w011
            ai = a011
        Case (12)
            wi = w012
            ai = a012
        Case (13)
            wi = w013
            ai = a013
        Case (14)
            wi = w014
            ai = a014
        Case (15)
            wi = w015
            ai = a015
        Case (16)
            wi = w016
            ai = a016
        Case (17)
            wi = w017
            ai = a017
        Case (18)
            wi = w018
            ai = a018
        Case (19)
            wi = w019
            ai = a019
        Case (20)
            wi = w020
            ai = a020
        Case (21)
            wi = w021
            ai = a021
        Case (22)
            wi = w022
            ai = a022
        Case (23)
            wi = w023
            ai = a023
        Case (24)
            wi = w024
            ai = a024
        Case (25)
            wi = w025
            ai = a025
        Case (26)
            wi = w026
            ai = a026
        Case (27)
            wi = w027
            ai = a027
        Case (28)
            wi = w028
            ai = a028
        Case (29)
            wi = w029
            ai = a029
        Case (30)
            wi = w030
            ai = a030
        Case (31)
            wi = w031
            ai = a031
        Case (32)
            wi = w032
            ai = a032
        Case (33)
            wi = w033
            ai = a033
        Case (34)
            wi = w034
            ai = a034
        Case (35)
            wi = w035
            ai = a035
        Case (36)
            wi = w036
            ai = a036
        Case (37)
            wi = w037
            ai = a037
        Case (38)
            wi = w038
            ai = a038
        Case (39)
            wi = w039
            ai = a039
        Case (40)
            wi = w040
            ai = a040
        Case (41)
            wi = w041
            ai = a041
        Case (42)
            wi = w042
            ai = a042
        Case (43)
            wi = w043
            ai = a043
        Case (44)
            wi = w044
            ai = a044
        Case (45)
            wi = w045
            ai = a045
        Case (46)
            wi = w046
            ai = a046
        Case (47)
            wi = w047
            ai = a047
        Case (48)
            wi = w048
            ai = a048
        Case (49)
            wi = w049
            ai = a049
        Case (50)
            wi = w050
            ai = a050
        Case (51)
            wi = w051
            ai = a051
        Case (52)
            wi = w052
            ai = a052
        Case (53)
            wi = w053
            ai = a053
        Case (54)
            wi = w054
            ai = a054
        Case (55)
            wi = w055
            ai = a055
        Case (56)
            wi = w056
            ai = a056
        Case (57)
            wi = w057
            ai = a057
        Case (58)
            wi = w058
            ai = a058
        Case (59)
            wi = w059
            ai = a059
        Case (60)
            wi = w060
            ai = a060
        Case (61)
            wi = w061
            ai = a061
        Case (62)
            wi = w062
            ai = a062
        Case (63)
            wi = w063
            ai = a063
        Case (64)
            wi = w064
            ai = a064
        Case (65)
            wi = w065
            ai = a065
        Case (66)
            wi = w066
            ai = a066
        Case (67)
            wi = w067
            ai = a067
        Case (68)
            wi = w068
            ai = a068
        Case (69)
            wi = w069
            ai = a069
        Case (70)
            wi = w070
            ai = a070
        Case (71)
            wi = w071
            ai = a071
        Case (72)
            wi = w072
            ai = a072
        Case (73)
            wi = w073
            ai = a073
        Case (74)
            wi = w074
            ai = a074
        Case (75)
            wi = w075
            ai = a075
        Case (76)
            wi = w076
            ai = a076
        Case (77)
            wi = w077
            ai = a077
        Case (78)
            wi = w078
            ai = a078
        Case (79)
            wi = w079
            ai = a079
        Case (80)
            wi = w080
            ai = a080
        Case (81)
            wi = w081
            ai = a081
        Case (82)
            wi = w082
            ai = a082
        Case (83)
            wi = w083
            ai = a083
        Case (84)
            wi = w084
            ai = a084
        Case (85)
            wi = w085
            ai = a085
        Case (86)
            wi = w086
            ai = a086
        Case (87)
            wi = w087
            ai = a087
        Case (88)
            wi = w088
            ai = a088
        Case (89)
            wi = w089
            ai = a089
        Case (90)
            wi = w090
            ai = a090
        Case (91)
            wi = w091
            ai = a091
        Case (92)
            wi = w092
            ai = a092
        Case (93)
            wi = w093
            ai = a093
        Case (94)
            wi = w094
            ai = a094
        Case (95)
            wi = w095
            ai = a095
        Case (96)
            wi = w096
            ai = a096
        Case (97)
            wi = w097
            ai = a097
        Case (98)
            wi = w098
            ai = a098
        Case (99)
            wi = w099
            ai = a099
        Case (100)
            wi = w100
            ai = a100
        Case (101)
            wi = w101
            ai = a101
        Case (102)
            wi = w102
            ai = a102
        Case (103)
            wi = w103
            ai = a103
        Case (104)
            wi = w104
            ai = a104
        Case (105)
            wi = w105
            ai = a105
        Case (106)
            wi = w106
            ai = a106
        Case (107)
            wi = w107
            ai = a107
        Case (108)
            wi = w108
            ai = a108
        Case (109)
            wi = w109
            ai = a109
        Case (110)
            wi = w110
            ai = a110
        Case (111)
            wi = w111
            ai = a111
        Case (112)
            wi = w112
            ai = a112
        Case (113)
            wi = w113
            ai = a113
        Case (114)
            wi = w114
            ai = a114
        Case (115)
            wi = w115
            ai = a115
        Case (116)
            wi = w116
            ai = a116
        Case (117)
            wi = w117
            ai = a117
        Case (118)
            wi = w118
            ai = a118
        Case (119)
            wi = w119
            ai = a119
        Case (120)
            wi = w120
            ai = a120
        Case (121)
            wi = w121
            ai = a121
        Case (122)
            wi = w122
            ai = a122
        Case (123)
            wi = w123
            ai = a123
        Case (124)
            wi = w124
            ai = a124
        Case (125)
            wi = w125
            ai = a125
        Case (126)
            wi = w126
            ai = a126
        Case (127)
            wi = w127
            ai = a127
        Case (128)
            wi = w128
            ai = a128
        Case (129)
            wi = w129
            ai = a129
        Case (130)
            wi = w130
            ai = a130
        Case (131)
            wi = w131
            ai = a131
        Case (132)
            wi = w132
            ai = a132
        Case (133)
            wi = w133
            ai = a133
        Case (134)
            wi = w134
            ai = a134
        Case (135)
            wi = w135
            ai = a135
        Case (136)
            wi = w136
            ai = a136
        Case (137)
            wi = w137
            ai = a137
        Case (138)
            wi = w138
            ai = a138
        Case (139)
            wi = w139
            ai = a139
        Case (140)
            wi = w140
            ai = a140
        Case (141)
            wi = w141
            ai = a141
        Case (142)
            wi = w142
            ai = a142
        Case (143)
            wi = w143
            ai = a143
        Case (144)
            wi = w144
            ai = a144
        Case (145)
            wi = w145
            ai = a145
        Case (146)
            wi = w146
            ai = a146
        Case (147)
            wi = w147
            ai = a147
        Case (148)
            wi = w148
            ai = a148
        Case (149)
            wi = w149
            ai = a149
        Case (150)
            wi = w150
            ai = a150
        Case Default
            Write(*,'(A,I0)') 'ERROR:  Quadratures: GaussLegendreN:  Unsupported N, n = ',N
            ERROR STOP
    End Select
    q = GaussLegendre(f,a,b,n,wi,ai)
    If (Present(p)) p = 2*n - 1
End Function GaussLegendreN

Function Progressive_GaussLegendre(f,a,b,aTol,rTol,n_start,n_stride,n_done,p) Result(q)
    !Integrates f(x) on (a,b) by successive quadrature via Gauss-Legendre n-point rules
    Use Kinds, Only: dp
    Use Utilities, Only: Converged
    Use Utilities, Only: Prec
    Implicit None
    Real(dp) :: q  !result of the integration
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function
    End Interface
    Real(dp), Intent(In) :: a,b  !limits of integration
    Real(dp), Intent(In) :: aTol,rTol  !absolute and relative tolerance to be used as criteria for convergence
    Integer, Intent(In), Optional :: n_start,n_stride
    Integer, Intent(Out), Optional :: n_done
    Real(dp), Intent(Out), Optional :: p
    Real(dp) :: q0
    Integer :: n,dn
    Integer :: i

    If (Present(n_start)) Then
        n = n_start
    Else
        n = 1
    End If
    If (Present(n_stride)) Then
        dn = n_stride
    Else
        dn = 1
    End If
    q0 = GaussLegendreN(n,f,a,b)
    Do i = n+dn,150,dn !Gauss-Legendre weights and abscissa are available up to 150 points
        q = GaussLegendreN(i,f,a,b)
        If (Converged(q,q0,rTol=rTol,aTol=aTol)) Then
            If (Present(n_done)) n_done = i
            If (Present(p)) Then
                If (i .LT. 150) Then !progress one more stage of GL points to estimate precision
                    q0 = q
                    q = GaussLegendreN(i+1,f,a,b)
                    p = Prec(q0,q)
                Else
                    p = Prec(q0,q)
                End If
            End If
            RETURN !NORMAL EXIT
        End If
        q0 = q
    End Do
    !If we get this far, we failed to converge on 150-point Gauss
    q0 = GaussLegendreN(149,f,a,b) !regress one stage of GL points to estimate precision
    Write(*,*)
    Write(*,'(A,I0,A)') 'WARNING:  Quadratures: Progressive_GaussLegendre:  Failed to converge with ',150,' Gauss-points.'
    Write(*,'(A10,A,ES23.15E2,A,F0.5,A)') '','Final estimated value: ',q,' (~',Prec(q,q0),' good digits)'
    Write(*,'(A10,A,2(ES10.3E2,A))')      '','Final estimated error:     ',Abs(q-q0),  ' (abs), ', &
                                   &                                       Abs(q-q0)/q,' (rel)'
    Write(*,'(A10,A,2(ES10.3E2,A))')      '','Convergence Criteria:      ',atol,       ' (abs), ', &
                                   &                                       rtol,       ' (rel)'
    If (Present(n_done)) n_done = -1
    If (Present(p)) p = Prec(q0,q)
    ! ERROR STOP
    !Currently, no STOP is implemented for a failure to converge.  The result of the quadrature is set to the final value from the
    !quadrature, and additional required values are set, and execution returns to the calling routine after the above warnings
    !are sent to the screen.
End Function Progressive_GaussLegendre

End Module Quadratures

max_stars_repo_name: wtdailey/commonTools

max_issues_repo_name: wtdailey/commonTools

max_forks_repo_name: wtdailey/commonTools


text:
!
!  file BDENSLIB5.f90
!
! utilities for data structures for two-body densities
!

module twobodydens_util
	use interaction
	use coupledmatrixelements

	implicit none

	integer :: dens2bfile = 47

	integer :: npnsection
    integer, allocatable :: pnpairstart(:), pnpairstop(:)
    integer(8), allocatable :: npnmesum(:)

contains

! ROUTINE FOR COUPLING TOGETHER pp/nn TWO-BODY DENSITIES
! here it is important to distinguish between ITBME, which labels uncoupled (m-scheme) matrix elements
! which are in turn stored in hmatpp/hmatnn
! and INDEX, which labels coupled (J-scheme) matrix elements
! stored in pp2bden and nn2bden
!
	subroutine couple_2bdensXX(it,JJi,JJf,zerome)

      use verbosity
      use spstate
      use haiku_info
      use system_parameters
      use interaction
      use coupledmatrixelements
      use butil_mod

      implicit none
      integer, intent(in)::  it      ! which species
	  integer, intent(in) :: JJi,JJf  ! 2 x initial, final J of states
      logical, intent(in) ::zerome  ! needed to reset
      integer(kind=8) :: dpair,cpair
	  integer Jabmin,Jabmax,Jcdmin,Jcdmax,Jab,Jcd
      integer Jtot,jmin,jmax
      integer T
      integer M
      integer par
      integer ik,il
      integer ii,ij

      integer i,j,k,l
      integer ith,jth,kth,lth
      integer isps,jsps,ksps,lsps

      integer nk,nl,ni,nj
      integer jk,jl,ji,jj
      integer mk,ml,mi,mj

      integer(kind=8) :: iref,istart
      integer(kind=8) :: itbme

      integer ia,ib,ic,id
	  integer :: ccouple,dcouple
	  integer :: abcouple,cdcouple

      integer phasekl,phaseij

      integer indx
      integer itmp
      real vtmp

      real cleb
	  real cleb4reduce
      real hfactor  ! Accounts for ordering for identical orbits
      real(8),pointer :: v(:)

	  type (vjs2b), pointer :: x2bden(:)
      integer pcref, pcstart, ppar
      integer,pointer :: xxmap(:)
      logical :: foundit
      integer :: aerr
	  integer(8) :: coupleblocksize,pairblocksize

	  logical :: absame,cdsame   ! if pairs have identical orbits, this restricts Js

      if(Np(it) < 2)return

      if(it ==1)then

        v => dmatpp
      else

        v => dmatnn

      endif

      if(it == 1)then
        x2bden => pp2bden
        xxmap => ppcouplemap

      else
        x2bden =>   nn2bden
        xxmap => nncouplemap

      endif

!----------------------------- ZERO OUT ARRAY-----------------

      T = 1        ! for pp, nn matrix elements

      do dpair = 1,npairXX(it)        ! loop over destruction m-scheme pair

        m = XX2(it)%pair(dpair)%m     ! find Jz, parity of pair
        par = XX2(it)%pair(dpair)%par

        iref = XX2(it)%meref(m,par)      ! iref, istart used to encode index for
        istart = XX2(it)%mestart(m,par)  ! the uncoupled matrix element
		pairblocksize= XX2(it)%block(m,par)

        if(iref > dpair)then
          print*,' problem with XX iref (d) ',it
          print*,dpair,iref
          print*,m,par
          print*,xx2(it)%mestart(:,1)
          stop
        endif

!---------------------- GET Q#s FOR PARTICLES IN DESTRUCTION PAIR --------
        k = XX2(it)%pair(dpair)%ia
        l = XX2(it)%pair(dpair)%ib
        if( k < 0)then
          ksps = -k
          kth = -it
        else
          ksps = k
          kth  = it
        endif
        ic = hspsqn(kth,ksps)%orb
        jk = hspsqn(kth,ksps)%j
        mk = hspsqn(kth,ksps)%m

        if( l < 0)then
          lsps = -l
          lth = -it
        else
          lsps = l
          lth  = it
        endif
        id = hspsqn(lth,lsps)%orb
        jl = hspsqn(lth,lsps)%j
        ml = hspsqn(lth,lsps)%m

		cdsame =.false.
		if(ic == id)cdsame = .true.

        phasekl = 1
!-------------- MUST HAVE ia > ib, else swap and pick up a phase ---

        if(ic < id)then
             itmp = ic
             ic = id
             id = itmp

             itmp = jk
             jk = jl
             jl = itmp

             itmp = mk
             mk   = ml
             ml   = itmp

             phasekl = -1  ! the rest of the phase comes from the Clebsch Gordan

        endif

        do cpair = 1,npairXX(it)       ! loop over creation m-scheme pair

          if(XX2(it)%pair(cpair)%m /=m)cycle       ! enforce quantum numbers
          if(XX2(it)%pair(cpair)%par /= par)cycle

          if(iref > cpair)then
            print*,' problem with XX iref (c) '
            print*,cpair,iref
            stop
          endif

!--------------- HOW INDEX FOR UNCOUPLED (m-scheme) tbme IS COMPUTED---------
          itbme = istart + pairblocksize*(dpair-iref-1)+cpair-iref

          if(itbme > nmatXX(it))then          ! error trap
             print*,' me label too large'
             print*,itbme,nmatXX(it)
             print*,cpair,dpair,iref,istart
             stop
          endif
          if(itbme <= 0)then
               print*,' problem itbme ',itbme
               stop
          endif

		  if(v(itbme)==0.0 .and. .not.zerome)cycle

!--------------- GET Q#s FOR CREATION PAIRS ---------------

          i = XX2(it)%pair(cpair)%ia
          j = XX2(it)%pair(cpair)%ib

!          print*,' operators ',i,j,k,l
          if( i < 0)then
            isps = -i
            ith = -it
          else
            isps = i
            ith  = it
          endif
          ia = hspsqn(ith,isps)%orb
          ji = hspsqn(ith,isps)%j
          mi = hspsqn(ith,isps)%m

          if( j < 0)then
            jsps = -j
            jth = -it
          else
            jsps = j
            jth  = it
          endif
          ib = hspsqn(jth,jsps)%orb
          jj = hspsqn(jth,jsps)%j
          mj = hspsqn(jth,jsps)%m

!-------------- account for double counting
          hfactor = 1.0
          !if(cpair == dpair)then
          !   hfactor = 0.5
          !else
          !   hfactor = 1.0
          !endif

		  if(ia==ib)hfactor = hfactor*2. ! this is because   a^+_i  a^+j is ordered
		  if(ic==id)hfactor = hfactor*2.

!------------- PUT INTO "STANDARD" ORDER; MAY GET PHASE -------

          phaseij = 1

		  absame = .false.
		  if(ia==ib)absame =.true.

          if(ia < ib) then

             itmp = ia
             ia = ib
             ib = itmp

             itmp = ji
             ji = jj
             jj = itmp
             itmp = mi
             mi   = mj
             mj   = itmp

             phaseij = -1  ! the rest of the phase comes from the Clebsch Gordan

          endif
          dcouple = ic*(ic-1)/2 + id  ! used to find index of coupled TBME

          ccouple = ia*(ia-1)/2 + ib

!-------------- ALSO MUST HAVE STANDARD ORDERING OF PAIRS -----------

          ccouple = xxmap(ccouple)
          dcouple = xxmap(dcouple)
          if(ccouple==-1 .or. dcouple == -1)then
             print*,' hmm strange ... (A)',it
             print*,ccouple,dcouple
             print*,ia,ib,ic,id
             print*,xxmap
             stop
          endif
          ppar = XXcouples(it)%pairc(ccouple)%par
          if(ppar /= XXcouples(it)%pairc(dcouple)%par)then
              print*,' sigh parity problem '
              print*,ia,ib,ic,id,ppar,par
              print*,ccouple,dcouple
              stop
          endif

          pcref = XXcouples(it)%meref(ppar)
          pcstart = XXcouples(it)%mestart(ppar)

          if(ccouple <= XXcouples(it)%meref(2))then
			coupleblocksize = XXcouples(it)%meref(2)
		  else
			coupleblocksize = ncouplesXX(it)-XXcouples(it)%meref(2)

		  end if
!.......... GET INDEX OF COUPLED (J-scheme) MATRIX ELEMENT.....

           indx = (dcouple - pcref -1 )*coupleblocksize +  		 &
		          ccouple - pcref + pcstart
!				  write(6,*)' INDEX AGAIN ',ccouple,dcouple,indx
          if(indx <= 0)then
             print*,' problem indx ',indx
             stop
          endif

!-------------FIND JMIN,JMAX ALLOWED ----------------
          Jabmax = (ji+jj)/2
		  Jabmin =abs(ji-jj)/2
		  Jcdmax = (jk+jl)/2
		  Jcdmin = abs(jk-jl)/2

!.... CHECK THESE AGAINST PRE-DETERMINED LIMITS

          if(Jabmax /= x2bden(indx)%Jabmax)then
			  print*,' mismatch in Jabmax ',Jabmax, x2bden(indx)%Jabmax,' indx ',indx
			  print*,ia,ib,ic,id
			  print*,dcouple,pcref,coupleblocksize,ccouple,pcstart

			  stop
		  end if

          if(Jabmin /= x2bden(indx)%Jabmin)then
			  print*,' mismatch in Jabmin ',Jabmin, x2bden(indx)%Jabmin
			  print*,i,j,ji,jj,indx
			  print*,coupleblocksize
			  stop
		  end if
          if(Jcdmax /= x2bden(indx)%Jcdmax)then
			  print*,' mismatch in Jcdmax ',Jcdmax, x2bden(indx)%Jcdmax
			  stop
		  end if
          if(Jcdmin /= x2bden(indx)%Jcdmin)then
			  print*,' mismatch in Jcdmin ',Jcdmin, x2bden(indx)%Jcdmin
			  stop
		  end if

!-------------NOW EXTRACT FROM TBMEs ---------------
!-------      NOTE FACTOR ZETA = SQRT(1+DELTA(A,B)) INCLUDED
          vtmp = v(itbme)*phasekl*phaseij *hfactor

		  if(vtmp==0.0 .and. .not.zerome)cycle

do Jab = Jabmin,Jabmax
	if(absame .and. (mod(Jab,2)/=0))cycle
	do Jcd = Jcdmin,Jcdmax
		if(cdsame .and. (mod(Jcd,2)/=0))cycle

        Jmin = x2bden(indx)%Jmin
		Jmin = max(Jmin,abs(Jab-Jcd))
        Jmax = x2bden(indx)%Jmax
		Jmax = min(Jmax,Jab+Jcd)
		if(Jmax < Jmin)cycle
		if(zerome)then
			x2bden(indx)%v(:,:,:) = 0.0
			cycle
		end if
          do jtot = jmin,jmax

			  cleb4reduce = cleb(jtot*2,0,JJi,Jz,JJf,Jz)
			  if(abs(cleb4reduce) < 0.00001)then
				  x2bden(indx)%v(Jab,Jcd,Jtot)= -999.
				  cycle
			  end if
			  x2bden(indx)%v(Jab,Jcd,Jtot)=x2bden(indx)%v(Jab,Jcd,Jtot) + &
			  vtmp * &
                    cleb(jk,mk,jl,ml,2*jcd,2*m)*cleb(ji,mi,jj,mj,2*jab,2*m)*cleb(2*jab,2*m,2*jcd,-2*m,2*Jtot,0)  &
					* (-1)**(Jcd -m)   & ! PHASE FROM time-reversal of second pair
					* sqrt(JJf+1.0)*(-1)**(Jtot+ (JJf-JJi)/2)/ cleb4reduce &    ! reduce
					/sqrt(2.0*Jtot+1.0)    ! needed for proper definition of density matrix

          enddo

    end do  ! Jcd
end do   ! Jab

        enddo  ! loop over cpair

      enddo  ! loop over dpair

      return

	end subroutine couple_2bdensXX

!
!==================================================================================
subroutine setup_PNarrays

	use interaction
	implicit none
	integer :: m, par
	integer :: nme
	integer :: isection,dimsection

	integer :: xpair,ipair,xstart

	integer :: aerr

      nme = 0

      m = -9999
      par = 0
!---------------- THE FOLLOWING SETS UP ARRAYS TO MAKE POSSIBLE OpenMP PARALLELIZATION--
!-------- COUNT HOW MANY SECTIONS OF PAIRS by m, parity
      isection = 0
      do xpair = 1,npairpn       ! loop over destruction pairs
!---------- FIND OUT IF THERE IS A CHANGE IN M OR PARITY
        if(par /= PN2%pair(xpair)%par .or. m /= PN2%pair(xpair)%m)then
          isection = isection +1
          m = PN2%pair(xpair)%m
          par = PN2%pair(xpair)%par
        end if
      end do
!---------- NOW ALLOCATE---------
      npnsection = isection
      if(.not. allocated( pnpairstart) )then
          allocate(pnpairstart(npnsection),pnpairstop(npnsection), stat=aerr )
          if(aerr /= 0) call memerror("count_uncouplepn 1")
          allocate(npnmesum(npnsection), stat=aerr )
          if(aerr /= 0) call memerror("count_uncouplepn 2")
      end if

      nme = 0

      m = -9999
      par = 0
!-------- COUNT HOW MANY SECTIONS OF PAIRS by m, parity
      dimsection = 0
	  isection = 0
      do xpair = 1,npairpn       ! loop over destruction pairs
!---------- FIND OUT IF THERE IS A CHANGE IN M OR PARITY
        if(par /= PN2%pair(xpair)%par .or. m /= PN2%pair(xpair)%m)then
          isection = isection +1

!------------- ADD TO TOTAL THOSE FROM LAST TIME
          nme = nme + (dimsection)*(dimsection)
          npnmesum(isection) = nme

!---------IF SO, THEN COUNT UP HOW MANY HAVE THE SAME M, PARITY
          m = PN2%pair(xpair)%m
          par = PN2%pair(xpair)%par
          xstart = xpair
          do ipair = xstart,npairpn  ! loop upwards until M, PARITY change
                                      ! pairs have previously been sorted
             if(par == PN2%pair(ipair)%par .and. m == PN2%pair(ipair)%m)then
!------------ NSECTION IS HOW MANY PAIRS HAVE SAME M, PARITY
                dimsection = ipair -xstart+1
             else
                exit
             endif
          enddo
          pnpairstart(isection) = xpair
          pnpairstop(isection)  = xpair-1+dimsection
          if(isection > 1)then
             if(pnpairstart(isection)/=pnpairstop(isection-1)+1)then
                print*,' Mismatch in pnpairstart/stop '
                print*,isection, pnpairstop(isection-1),pnpairstart(isection)
                stop
             end if
          end if
        endif

      end do

	return

end subroutine setup_PNarrays

!==================================================================================
! ROUTINE FOR COUPLING TOGETHER pn TWO-BODY DENSITIES

	subroutine couple_2bdensPN(JJi,JJf,zerome)
        use verbosity
        use spstate
        use haiku_info
        use system_parameters
        use interaction
        use coupledmatrixelements
        use butil_mod
		use sporbit

		implicit none
        integer, intent(in) :: JJi,JJf  ! 2 x initial, final J of states
        logical, intent(in) ::zerome  ! needed to reset

!...................... INTERNAL..........................
        real cleb

        integer :: dpair,cpair
	    integer :: m, par
	    integer :: ia,ib,ic,id
	    integer :: i,j,k,l
	  	integer :: isps,ith,jsps,jth,ksps,kth,lsps,lth
		integer :: ji,jj,jk,jl,mi,mj,mk,ml
		integer :: phaseij,phasekl
		integer :: pair1,pair2,pair11,pair22
		integer :: itbme,indx
		integer :: pcref,pcstart
		integer :: cindx,dindx
		integer :: iref, istart

		logical :: first
		integer :: ppar
		integer :: nsection, nme

		integer :: Jab,Jcd,Jtot
		integer :: Jabmin,Jabmax,Jcdmin,Jcdmax,Jmin,Jmax

		real vtmp

		integer :: pairblocksize

		integer :: isection
		real :: cleb4reduce

        do isection = 1,npnsection
	    dindx = 1
        do dpair = pnpairstart(isection),pnpairstop(isection) ! loop over destruction pairs
            m = PN2%pair(dpair)%m
            par = PN2%pair(dpair)%par

!-------------- EXTRACT QUANTUM NUMBERS OF DESTRUCTION PAIR--------
           k = PN2%pair(dpair)%ia
           l = PN2%pair(dpair)%ib
           if( k < 0)then
              ksps = -k
              kth = -1
           else
              ksps = k
              kth  = 1
           endif
           ic = hspsqn(kth,ksps)%orb
           jk = hspsqn(kth,ksps)%j
           mk = hspsqn(kth,ksps)%m
!           gk = hspsqn(kth,ksps)%group

           if( l < 0)then
             lsps = -l
             lth = -2
           else
             lsps = l
             lth  = 2
           endif
           id = hspsqn(lth,lsps)%orb
           jl = hspsqn(lth,lsps)%j
           ml = hspsqn(lth,lsps)%m
!           gl = hspsqn(lth,lsps)%group
           phasekl = 1

!-------------- FIND INDEX "PAIR1" used to find coupled tbme

           pair1 = numorb(2)*(ic-1)+id
!		   print*,ic,id,pair1

           pair1 = PNcouplemap(pair1)
           if(pair1 == -1)then
              print*,' uh problem pn map boss 1 '
              stop
           endif

		   Jcdmax = (jk+jl)/2
		   Jcdmin = abs(jk-jl)/2
		   Jcdmin = max(Jcdmin,abs(m))

	  	   if(par ==1)then
	           pairblocksize = PNcouples%meref(2)
	  	   else
	           pairblocksize = ncouplesPN-PNcouples%meref(2)

	  	   end if

!----------- SET UP MAPPING ARRAYS----------------------
!           dpnpair(l,k) = dindx + nme
!------------- TIME REVERSE -----------
!           if(useTR .and. m== 0)then
!             phasepnpair(l,k) = 1
!           endif
!          if(useTR .and. m/= 0)then
!            phasepnpair(l,k) = 1
!            ktr = hspsqn(kth,ksps)%tr
!            if(kth >0 .and. mk /= 0)ktr = -ktr
!            ltr = hspsqn(lth,lsps)%tr
!            if(lth >0 .and. ml /= 0)ltr = -ltr
!            phasepnpair(ltr,ktr) = int((-1)**( (jl+jk)/2),1)
!            dpnpair(ltr,ktr) = dindx + nme
  !          print*,l,k,dpnpair(l,k),ltr,ktr,dpnpair(ltr,ktr)
!          endif
!----------------- LOOP OVER CREATION PAIRS--------------
           cindx = 0
           do cpair = pnpairstart(isection),pnpairstop(isection)
			   if(m/= PN2%pair(cpair)%m)cycle
			   if(par/= PN2%pair(cpair)%par)cycle

              cindx = cindx + 1

              itbme = dindx + (cindx-1)*(pnpairstop(isection)+1-pnpairstart(isection)) + npnmesum(isection)
!              itbmetr = cindx + (dindx-1)*nsection +nme ! "time-reversed"

              if(itbme > nmatpn)then ! .or. itbmetr > nmatpn)then
                 print*,' me label too large'
                print*,itbme, nmatpn
                print*,cindx,dindx,nme,m,par
                print*,cpair,dpair,iref,istart
                stop
              endif
              if(itbme <= 0)then
                 print*,' problem itbme ',itbme
                 stop
              endif
!-------------- EXTRACT QUANTUM NUMBERS OF CREATION PAIR--------

              i = PN2%pair(cpair)%ia
              j = PN2%pair(cpair)%ib

              if( i < 0)then
                isps = -i
                ith = -1
              else
                isps = i
                ith  = 1
              endif
              ia = hspsqn(ith,isps)%orb
              ji = hspsqn(ith,isps)%j
              mi = hspsqn(ith,isps)%m
!              gi = hspsqn(ith,isps)%group
              if( j < 0)then
                jsps = -j
                jth = -2
              else
               jsps = j
               jth  = 2
              endif
              ib = hspsqn(jth,jsps)%orb
              jj = hspsqn(jth,jsps)%j
              mj = hspsqn(jth,jsps)%m
!              gj = hspsqn(jth,jsps)%group

!------------------ SET MAPPING ARRAY
!              cpnpair(j,i) = (cindx-1)*nsection
!            if(useTR .and. m/= 0)then   ! NOTE FOR NOW, TR NOT USED
!              itr = hspsqn(ith,isps)%tr
!              if(ith >0 .and. mi/=0)itr = -itr
!              jtr = hspsqn(jth,jsps)%tr
!              if(jth >0 .and. mj /= 0)jtr = -jtr
!              cpnpair(jtr,itr) = (cindx-1)*nsection
!            endif

!---------------- BECAUSE I HAVE TIME REVERSAL, I CAN SKIP A LOT
!            if(cpair > dpair)cycle

!------------ IF TWO-BODY DENSITIES, DO NOT UNCOUPLE MATRIX ELEMENTS

!            if(menu_char=='2b')cycle

!------------- PUT INTO "STANDARD" ORDER; MAY GET PHASE -------

             phaseij = 1

             pair2 = numorb(2)*(ia-1)+ib
             pair2 = PNcouplemap(pair2)
             if(pair2==-1)then
                print*,' oh boy should not have gotten that '
                stop
             endif

!            ppar = PNcouples%pairc(pair1)%par
             pcref = PNcouples%meref(par)
             pcstart = PNcouples%mestart(par)
!  MUST FIX
             indx = pairblocksize*(pair1-pcref-1)+pair2-pcref +pcstart

!             indx = istart + (dpair-iref-1)*pairblocksize+cpair-iref

             if(indx <= 0 .or. indx > nv2bmedim(0))then
                print*,' problem indx ',indx,pair1,pair2
                print*,ia,ib,ic,id
                stop
             endif

			 if(indx > nv2bmedim(0))then
				 print*,'error ',pairblocksize,dpair,cpair,indx,nv2bmedim(0)

				 print*,PNcouples%meref(2) , ncouplesPN - PNcouples%meref(2)

!				 print*,npnsection,pnpairstart,pnpairstop

			 end if

!-------------FIND JMIN,JMAX ALLOWED ----------------

!             jmax = bmin( jk+jl,ji+jj) /2
!             jmin = bmax( abs(jk-jl),abs(ji-jj))/2
!             jmin = bmax(jmin,abs(m))

!-------------NOW EXTRACT FROM TBMEs ---------------
!-------      NOTE FACTOR ZETA = SQRT(1+DELTA(A,B)) INCLUDED
!if(.not.zerome)print*,itbme,hmatpn(itbme),i,j,k,l
              vtmp = dmatpn(itbme) !*phasekl*phaseij !*hfactor

              if(vtmp==0.0 .and. .not.zerome)cycle

  		     Jabmax = (ji+jj)/2
  		     Jabmin = abs(ji-jj)/2

			 Jabmin = max(Jabmin,abs(m))

!			 print*,'test ', indx,ic,id,jcdmin,jcdmax

			 if(jcdmax > pn2bden(indx)%Jcdmax)then

				 print*,' whoops ',indx
				 print*,pair1,pair2,pairblocksize
				 print*,pair1,ic,id,Jcdmax,pair2,pairblocksize
				 stop
			 end if
!			 if(.not.zerome)print*,' test ',vtmp,i,j,k,l
			 do Jab = Jabmin,Jabmax
			 	do Jcd = Jcdmin,Jcdmax

			         Jmin = pn2bden(indx)%Jmin
			 		Jmin = max(Jmin,abs(Jab-Jcd))
			         Jmax = pn2bden(indx)%Jmax
			 		Jmax = min(Jmax,Jab+Jcd)
			 		if(Jmax < Jmin)cycle
			 		if(zerome)then
			 			pn2bden(indx)%v(:,:,:) = 0.0
			 			cycle
			 		end if
			           do jtot = jmin,jmax

			 			  cleb4reduce = cleb(jtot*2,0,JJi,Jz,JJf,Jz)
			 			  if(abs(cleb4reduce) < 0.00001)then
			 				  pn2bden(indx)%v(Jab,Jcd,Jtot)= -999.
			 				  cycle
			 			  end if
			 			  pn2bden(indx)%v(Jab,Jcd,Jtot)=pn2bden(indx)%v(Jab,Jcd,Jtot) + &
			 			  vtmp * &
			                     cleb(jk,mk,jl,ml,2*jcd,2*m)*cleb(ji,mi,jj,mj,2*jab,2*m)*cleb(2*jab,2*m,2*jcd,-2*m,2*Jtot,0)  &
			 					* (-1)**(Jcd -m)   & ! PHASE FROM time-reversal of second pair
			 					* sqrt(JJf+1.0)*(-1)**(Jtot+ (JJf-JJi)/2)/ cleb4reduce &    ! reduce
			 					/sqrt(2.0*Jtot+1.0)    ! needed for proper definition of density matrix

			           enddo

			     end do  ! Jcd
			 end do   ! Jab

          enddo  ! loop over cpair
            dindx = dindx + 1
        enddo  ! loop over dpair
	end do ! isection

		return
	end subroutine couple_2bdensPN
!==================================================================================

!........ SET UP TWO-BODY DENSITY MATRIX......
!         because these files can be large, create a separate file

    subroutine setup_2bdens_output

		use nodeinfo
		use io
		use program_info
		use sporbit
		use densities,only:pndensities
		use system_parameters
		implicit none
		integer :: ilast
	    character(10) :: date,time,zone
	    integer :: datetimeval(8)
		integer :: iorb

		if(iproc/=0)return

!.......... OPEN FILE...................

        ilast = index(outfile,' ')-1

        open(unit=dens2bfile,file=outfile(1:ilast)//".den2b",status = 'unknown')

		print*,' '
		print*,' Two-body densities will be written to: ',outfile(1:ilast),".den2b "
		if(diag_den2b)print*,' (Expectation values of two-body operators written to ',outfile(1:ilast),".res  )"
		print*, ' '

!...........WRITE HEADER WITH S.P. INFORMATION.....
        write(dens2bfile,*)'!#  Two-body densities from BIGSTICK version ',version,lastmodified
	    call date_and_time(date,time,zone,datetimeval)
	    write(dens2bfile,*)'!#  Run date: ',date(1:4),'-',date(5:6),'-',date(7:8)
			write(dens2bfile,*)'!# Densities written in explicit proton-neutron formalism '
			write(dens2bfile,*)'!# Single-particle orbits information follows '
			write(dens2bfile,*)'!# Number of single-particle orbits (same for both protons, neutrons ) '
			write(dens2bfile,*)numorb(1)
		    write(dens2bfile,*)'!#      '

			write(dens2bfile,*)'!# Proton, neutron    N      L   2xJ '
			do iorb = 1,numorb(1)
				write(dens2bfile,'(5x,2i5,2x,3i7)')iorb,iorb+numorb(1),orbqn(1,iorb)%nr,orbqn(1,iorb)%l,orbqn(1,iorb)%j
			end do

			if(diag_den2b .and. pndensities)then
			    write(resultfile,*)'!#    orbit labels '

				write(resultfile,*)'!# Proton, neutron    N      L   2xJ '
				do iorb = 1,numorb(1)
					write(resultfile,'(5x,2i5,2x,3i7)')iorb,iorb+numorb(1),orbqn(1,iorb)%nr,orbqn(1,iorb)%l,orbqn(1,iorb)%j
				end do

			end if
			if(diag_den2b .and. .not.pndensities)then
				write(resultfile,*)'!#    Orbit     N       L    2xJ '
				do iorb = 1,numorb(1)
					write(resultfile,'(5x,i5,2x,3i7)')iorb,orbqn(1,iorb)%nr,orbqn(1,iorb)%l,orbqn(1,iorb)%j
				end do

			end if

		write(dens2bfile,*)'!#   Z   N '
		write(dens2bfile,'(3x,2i4)')np(1),np(2)

		return

	end subroutine setup_2bdens_output

! ROUTINE FOR PRINTING OUT COUPLED TWO-BODY DENSITIES
! must have already called couple_2bdnesXX,PN

	subroutine print_out_2bdens(istate,fstate,Ei,Ef,xJi,xJf,xTi,xTf)

		use interaction
		use coupledmatrixelements
		use sporbit
		use densities,only:pndensities
		use nodeinfo
		use system_parameters
		use io
		implicit none
		integer,intent(in) :: istate,fstate  ! label of initial,final many-body state
		real(4),intent(in) :: Ei,Ef          ! energies of initial,final state
		real(4),intent(in) :: xJi,xJf,xTi,xTf  ! J,T (if not in pn-format) for initial,final state
		integer :: it
        integer(kind=4) :: abcouple,cdcouple  ! "couples" of orbits
		integer(kind=4) :: a,b,c,d
		integer(kind=4) :: Jab,Jcd
        integer Jtot,jmin,jmax
		integer(8) :: indx, parblocksize

		type (vjs2b), pointer :: x2bden(:)

		integer(4) :: ipar
		integer(4) :: couplestart,coupleend

		real(8) :: test,sumrule,zetafactor,trfactor,zetaab,zetacd
		real(8) :: sumrulepp,sumrulenn,sumrulepn  ! intermediate for testing
		integer :: i

		sumrule =0.d0
		sumrulepn = 0.d0
		sumrulepp = 0.d0
		if(iproc/=0)return
		write(dens2bfile,*)'!#  '
		write(dens2bfile,*)'!#  '

		if(diag_den2b)then
  		   write(dens2bfile,*)'!#   State      Energy         Jstate     '
  		   write(dens2bfile,'(i10,6x,f10.5,5x,f5.1)')istate,Ei,xJi
 		   write(dens2bfile,*)'!# a   b   c   d   Jpair   rho_(Jpair,Jpair:J=0)   '

  		   write(resultfile,*)'!#   State      Energy         Jstate     '
  		   write(resultfile,'(i10,6x,f10.5,5x,f5.1)')istate,Ei,xJi
 		   if(pndensities)then
			   write(resultfile,*)'!# a   b   c   d     J     X_J(ab,cd) '
		   else
			   write(resultfile,*)'!# a   b   c   d     J     T      X_JT(ab,cd) '
		   endif

	   else
		   write(dens2bfile,*)'!# Ini state      Energy         Ji     '
		   write(dens2bfile,'(i10,6x,f10.5,5x,f5.1)')istate,Ei,xJi
		   write(dens2bfile,*)'!# Fin state      Energy         Jf     '
		   write(dens2bfile,'(i10,6x,f10.5,5x,f5.1)')fstate,Ef,xJf
		   write(dens2bfile,*)'!# a   b    Jab    c   d    Jcd  Jmin Jmax  rho_(Jab,Jcd:J), J=Jmin,Jmax)'

		end if

		do it = 1,2
			if(np(it) < 2)cycle

			if(it==1)then
				x2bden=>pp2bden

			else
				x2bden=>nn2bden

			end if

		    do abcouple = 1,ncouplesXX(it)
				a = XXcouples(it)%pairc(abcouple)%ia
				b = XXcouples(it)%pairc(abcouple)%ib

				if(it==2)then
					a = a+numorb(1)
					b = b+numorb(1)

				end if

!--- GET PARITY

                if(abcouple <= XXcouples(it)%meref(2))then
					ipar = 1
					parblocksize = XXcouples(it)%meref(2)
					couplestart = 1
					coupleend   = XXcouples(it)%meref(2)
				else
					ipar = 2
					parblocksize = ncouplesXX(it)-XXcouples(it)%meref(2)
					couplestart = XXcouples(it)%meref(2)+1
					coupleend   = ncouplesXX(it)

				end if

				do cdcouple= couplestart,coupleend !1,ncouplesXX(it)

					c = XXcouples(it)%pairc(cdcouple)%ia
					d = XXcouples(it)%pairc(cdcouple)%ib

					zetafactor = 1.0
					if(a==b)then
						zetafactor = zetafactor/sqrt(2.)
					end if
					if(c==d)then
						zetafactor = zetafactor/sqrt(2.)
					end if
					if(it==2)then
						c = c+numorb(1)
						d = d+numorb(1)

					end if

					if(diag_den2b)then  ! prevent writing out time-reversed states
						if(a < c)cycle
						if(a==c .and. b < d)cycle
					end if

					if(a==c .and. b==d)then
						trfactor = 1
					else
						trfactor = 2   ! used to account for Time reversal  (ab,cd) = (cd,ab)
					end if

!--- COMPUTE INDEX ........

                    indx = (cdcouple - XXcouples(it)%meref(ipar) -1 )*parblocksize +  		 &
					       abcouple - XXcouples(it)%meref(ipar) + XXcouples(it)%mestart(ipar)

!----- LOOP OVER Js
                    do Jab = x2bden(indx)%Jabmin,x2bden(indx)%Jabmax

	                    do Jcd = x2bden(indx)%Jcdmin,x2bden(indx)%Jcdmax
							if(diag_den2b)then
								Jmin=0
								Jmax= 0
							else
							    Jmin = abs(Jab-Jcd)
							    Jmin = max(Jmin, int(abs(xJi-xJf)))
						 	    Jmax = (Jab + Jcd)
							    Jmax = min(Jmax,int(xJi+xJf))
							end if

							if(Jmin > Jmax)cycle
							if(Jmin < x2bden(indx)%Jmin .or.   Jmax > x2bden(indx)%Jmax )then  !ERROR TRAP
							   print*,' ALARUM!  Some problem with min, max Js '
							   print*,Jab,Jcd,Jmin,Jmax

							   stop
						    end if
							test =  0.d0
							do Jtot= Jmin,Jmax
								if(x2bden(indx)%v(Jab,Jcd,Jtot) < -900)cycle
								test = test + abs(x2bden(indx)%v(Jab,Jcd,Jtot))
							end do
							if(test < 1.0e-5)then
								cycle
							endif
							if(diag_den2b)then
								write(dens2bfile,'(x,4i4,i6,2x,f12.5)')a,b,c,d,Jab, &
								   (x2bden(indx)%v(Jab,Jcd,0))*zetafactor

   								if(pndensities)write(resultfile,'(x,4i4,i6,2x,f12.5)')a,b,c,d,Jab, &
   								   x2bden(indx)%v(Jab,Jcd,0)*sqrt(2*Jab+1.0)*zetafactor*trfactor/sqrt(2*xJi+1.0)
								   if(a==c .and. b==d)then
									   sumrule= sumrule + 2*sqrt(2*Jab+1.0)*zetafactor*x2bden(indx)%v(Jab,Jcd,0)/sqrt(2*xJi+1.0)
									   if(it==1)sumrulepp= sumrulepp + 2*sqrt(2*Jab+1.0)*zetafactor*x2bden(indx)%v(Jab,Jcd,0)/sqrt(2*xJi+1.0)

								   end if
							else

							write(dens2bfile,'(x,2i4,i6,2x,2i4,i6,2x,2i4,2x,20f12.5)')a,b,Jab,c,d,Jcd,Jmin,Jmax, &
							   (x2bden(indx)%v(Jab,Jcd,Jtot)*zetafactor,Jtot=Jmin,Jmax)
						    end if
						end do ! Jcd

					end do !  Jab

				end do

			end do    !ab couple
		end do   !it
!...................  NOW FOR PROTON-NEUTRON

        if(np(1)*np(2) > 0)then
		    do abcouple = 1,ncouplesPN
				a = PNcouples%pairc(abcouple)%ia
				b = PNcouples%pairc(abcouple)%ib

				zetafactor = 1.0

			    b = b+numorb(1)
!--- GET PARITY

                if(abcouple <= PNcouples%meref(2))then
					ipar = 1
					parblocksize = PNcouples%meref(2)
					couplestart = 1
					coupleend   = PNcouples%meref(2)
				else
					ipar = 2
					parblocksize = ncouplesPN-PNcouples%meref(2)
					couplestart = PNcouples%meref(2)+1
					coupleend   = ncouplesPN

				end if

				do cdcouple= couplestart,coupleend !1,ncouplesXX(it)

					c = PNcouples%pairc(cdcouple)%ia
					d = PNcouples%pairc(cdcouple)%ib

					d = d+numorb(1)

					if(diag_den2b)then  ! prevent writing out time-reversed states
						if(a < c)cycle
						if(a==c .and. b < d)cycle
					end if
					if(a==c .and. b==d)then
						trfactor = 1
					else
						trfactor = 2   ! used to account for Time reversal  (ab,cd) = (cd,ab)
					end if
!--- COMPUTE INDEX ........

                    indx = (abcouple - PNcouples%meref(ipar) -1 )*parblocksize +  		 &
					       cdcouple - PNcouples%meref(ipar) + PNcouples%mestart(ipar)

!----- LOOP OVER Js
                    do Jab = pn2bden(indx)%Jabmin,pn2bden(indx)%Jabmax

	                    do Jcd = pn2bden(indx)%Jcdmin,pn2bden(indx)%Jcdmax
							if(diag_den2b)then
							   Jmin=0
							   Jmax= 0
							else
							   Jmin = abs(Jab-Jcd)
							   Jmin = max(Jmin, int(abs(xJi-xJf)))
							   Jmax = (Jab + Jcd)
							   Jmax = min(Jmax,int(xJi+xJf))
						    end if
							if(Jmin > Jmax)cycle
							if(Jmin < pn2bden(indx)%Jmin .or.   Jmax > pn2bden(indx)%Jmax )then  !ERROR TRAP
							   print*,' ALARUM!  Some problem with min, max Js in pn'
							   print*,Jab,Jcd,Jmin,Jmax
							   stop
						    end if
							! check to see if nonzero densities
							test =  0.d0
							do Jtot= Jmin,Jmax
								if(pn2bden(indx)%v(Jab,Jcd,Jtot) < -900)cycle
								test = test + abs(pn2bden(indx)%v(Jab,Jcd,Jtot))
							end do
!							print*,test
							if(test < 1.0e-5)then
								cycle
							endif
							if(diag_den2b)then
								write(dens2bfile,'(x,4i4,i6,2x,f12.5)')a,b,c,d,Jab, &
								   (pn2bden(indx)%v(Jab,Jcd,0))
      							if(pndensities)write(resultfile,'(x,4i4,i6,2x,f12.5)')a,b,c,d,Jab, &
      								   pn2bden(indx)%v(Jab,Jcd,0)*sqrt(2*Jab+1.0)/sqrt(2*xJi+1.0)  *trfactor !*zetafactor
								   if(a==c .and. b==d)then
									   sumrule= sumrule + 2*sqrt(2.0*Jab+1.0)*pn2bden(indx)%v(Jab,Jcd,0)/sqrt(2*xJi+1.0)
									   sumrulepn= sumrulepn + 2*sqrt(2.0*Jab+1.0)*pn2bden(indx)%v(Jab,Jcd,0)/sqrt(2*xJi+1.0)

								   end if

							else
							  write(dens2bfile,'(x,2i4,i6,2x,2i4,i6,2x,2i4,2x,20f12.5)')a,b,Jab,c,d,Jcd,Jmin,Jmax, &
							   (pn2bden(indx)%v(Jab,Jcd,Jtot),Jtot=Jmin,Jmax)
						    endif
						end do ! Jcd

					end do !  Jab

				end do

			end do    !ab couple
		end if
		sumrule = sumrule
		if(diag_den2b)then
			print*,'!# Sum rule = ',sumrule
			write(dens2bfile,*)'!# Sum rule = ',sumrule
		end if
		write(dens2bfile,*)'!#  '
		write(dens2bfile,*)'!# Two-body density matrices defined in  '
		write(dens2bfile,*)'!# BIGSTICK Manual, section 5.2, eqn (5.19) '
		write(dens2bfile,*)'!# rho(Jab,Jcd:J) = '
		write(dens2bfile,*)'!# - < Jf || [ [a^+ x b^+]_Jab x [ ~c x ~d ]_Jcd ]_J || Ji> '
		write(dens2bfile,*)'!#  ---------------------------------------------------------'
		write(dens2bfile,*)'!#     sqrt[ (2J +1) ( 1+ delta_ab)(1+delta_cd ) ] '

		if(diag_den2b .and. .not.pndensities)call isospin_expectations(xJi)
		return
	end subroutine print_out_2bdens

!=============================================================
!
!  added 7.9.5 to allow us to write out diagonal (expectation value) matrix elements in isospin format
!

subroutine isospin_expectations(xJi)
	use interaction
	use coupledmatrixelements
	use sporbit
	use densities,only:pndensities
	use nodeinfo
	use system_parameters
	use io

	real,intent(in) :: xJi
	integer :: ntbmes,ncpair
	integer :: dpair,cpair
	integer :: dparity
	integer :: iref, istart
	integer :: a,b,c,d
	integer :: jmin,jmax
	integer :: ja,jb,jc,jd,jj
	integer :: parblocksize
	integer :: couplestart,coupleend
	integer :: pncouplestart,pncoupleend
	integer :: pnparblocksize
	integer :: ipar
	real :: trfactor,zetafactor
	logical :: same
	integer :: phase

	integer :: abcouple,cdcouple, pnabcouple,pncdcouple,pnabcoupleR,pncdcoupleR

	real(4) :: isoscale(0:30),isovec(0:30)  ! tmp storage for matrix element

	integer :: indx,pnindx

	integer :: aerr

	if(ncouplesXX(1) /= ncouplesXX(2))then
		print*,' Cannot do isospin format in this case '
		stop
	end if
    do abcouple = 1,ncouplesXX(1)
		a = XXcouples(1)%pairc(abcouple)%ia
		b = XXcouples(1)%pairc(abcouple)%ib

		ja = orbqn(1,a)%j
		jb = orbqn(1,b)%j

        if(abcouple <= XXcouples(1)%meref(2))then
			ipar = 1
			parblocksize = XXcouples(1)%meref(2)
			couplestart = 1
			coupleend   = XXcouples(1)%meref(2)
		else
			ipar = 2
			parblocksize = ncouplesXX(1)-XXcouples(1)%meref(2)
			couplestart = XXcouples(1)%meref(2)+1
			coupleend   = ncouplesXX(1)

		end if

		do cdcouple= couplestart,coupleend !1,ncouplesXX(it)
			isovec= 0.0
			isoscale = 0.0

			c = XXcouples(1)%pairc(cdcouple)%ia
			d = XXcouples(1)%pairc(cdcouple)%ib

			if(a==b .or. c==d)then
				same = .true.
			else
				same = .false.
			end if

			jc = orbqn(1,c)%j
			jd = orbqn(1,d)%j

			jmin = max(abs(ja-jb),abs(jc-jd))/2
			jmax = min(abs(ja+jb),abs(jc+jd))/2
			if(jmax < jmin )cycle

			zetafactor = 1.0
			if(a==b)then
				zetafactor = zetafactor/sqrt(2.)
			end if
			if(c==d)then
				zetafactor = zetafactor/sqrt(2.)
			end if

			if(a < c)cycle
			if(a==c .and. b < d)cycle

			if(a==c .and. b==d)then
				trfactor = 1
			else
				trfactor = 2   ! used to account for Time reversal  (ab,cd) = (cd,ab)
			end if
			isoscale = 0.0

            indx = (cdcouple - XXcouples(1)%meref(ipar) -1 )*parblocksize +  		 &
			       abcouple - XXcouples(1)%meref(ipar) + XXcouples(1)%mestart(ipar)

			do jj = jmin,jmax
				if(same .and. mod(jj,2)==1)cycle  ! no T=1 matrix elements

				isovec(jj)= (pp2bden(indx)%v(jj,jj,0) + nn2bden(indx)%v(jj,jj,0)) &
				          *sqrt(2*Jj+1.0)*zetafactor*trfactor/sqrt(2*xJi+1.0)

			end do
!..................... NOW NEED TO FETCH PROTON-NEUTRON INFORMATION.............
            if(ipar==1)then
	           pnparblocksize = PNcouples%meref(2)
	           pncouplestart = 1
	           pncoupleend   = PNcouples%meref(2)
            else
	           pnparblocksize = ncouplesPN-PNcouples%meref(2)
	           pncouplestart = PNcouples%meref(2)+1
	           pncoupleend   = ncouplesPN

            end if
			do pnabcouple = pncouplestart,pncoupleend
				if(a==PNcouples%pairc(pnabcouple)%ia .and. b==PNcouples%pairc(pnabcouple)%ib)goto 1
			end do ! pnabcouples
			print*,' (1) Could not find couple ',a,b
			stop

			1 continue
			do pncdcouple = pncouplestart,pncoupleend
				if(c==PNcouples%pairc(pncdcouple)%ia .and. d==PNcouples%pairc(pncdcouple)%ib)goto 2
			end do ! pncdcouples
			print*,' (2) Could not find couple ',c,d
			stop

			2 continue
            pnindx = (pnabcouple - PNcouples%meref(ipar) -1 )*pnparblocksize +  		 &
			       pncdcouple - PNcouples%meref(ipar) + PNcouples%mestart(ipar)

	   		do jj = jmin,jmax
	   			if(same .and. mod(jj,2)==1)cycle  ! no T=1 matrix elements
!				if(jj==0)print*,' TEST Y1 ', pn2bden(pnindx)%v(jj,jj,0) *sqrt(2*Jj+1.0)/sqrt(2*xJi+1.0)

				isovec(jj) = isovec(jj) + (1./2.)*pn2bden(pnindx)%v(Jj,Jj,0)*sqrt(2*Jj+1.0)/sqrt(2*xJi+1.0)  *trfactor *zetafactor

			end do
	   		do jj = jmin,jmax
	   			if(same .and. mod(jj,2)==0)cycle  ! no T=0 matrix elements
				isoscale(jj) = isoscale(jj) + (1./2.)*pn2bden(pnindx)%v(Jj,Jj,0)*sqrt(2*Jj+1.0)/sqrt(2*xJi+1.0)  *trfactor *zetafactor
			end do
!............... SWAP   c -- > d
            phase = (-1)**( (jc+jd)/2)
			do pncdcoupleR = pncouplestart,pncoupleend
				if(d==PNcouples%pairc(pncdcoupleR)%ia .and. c==PNcouples%pairc(pncdcoupleR)%ib)goto 3
			end do ! pncdcouples
			print*,' (3) Could not find couple ',d,c
			stop

			3 continue
            pnindx = (pnabcouple - PNcouples%meref(ipar) -1 )*pnparblocksize +  		 &
			       pncdcoupleR - PNcouples%meref(ipar) + PNcouples%mestart(ipar)

	   		do jj = jmin,jmax
	   			if(same .and. mod(jj,2)==1)cycle  ! no T=1 matrix elements
!				if(jj==0)print*,' TEST Y2 ', pn2bden(pnindx)%v(jj,jj,0) , sqrt(2*Jj+1.0)/sqrt(2*xJi+1.0),trfactor ,zetafactor, &
!				(-1)**JJ*phase*0.5
				isovec(jj) = isovec(jj) - (1./2.)*pn2bden(pnindx)%v(Jj,Jj,0)*sqrt(2*Jj+1.0)/sqrt(2*xJi+1.0)*trfactor *zetafactor*(-1)**JJ*phase
			end do
	   		do jj = jmin,jmax
	   			if(same .and. mod(jj,2)==0)cycle  ! no T=0 matrix elements
				isoscale(jj)=isoscale(jj)+(1./2.)*pn2bden(pnindx)%v(Jj,Jj,0)*sqrt(2*Jj+1.0)/sqrt(2*xJi+1.0)*trfactor *zetafactor*(-1)**JJ*phase
			end do
!............... SWAP   a -- > b, c <---> d
            phase = (-1)**( (ja+jb+jc+jd)/2)

			do pnabcoupleR = pncouplestart,pncoupleend
				if(b==PNcouples%pairc(pnabcoupleR)%ia .and. a==PNcouples%pairc(pnabcoupleR)%ib)goto 4
			end do ! pnabcouples
			print*,' (4) Could not find couple ',b,a
			stop

			4 continue
            pnindx = (pnabcoupleR - PNcouples%meref(ipar) -1 )*pnparblocksize +  		 &
			       pncdcoupleR - PNcouples%meref(ipar) + PNcouples%mestart(ipar)

	   		do jj = jmin,jmax
	   			if(same .and. mod(jj,2)==1)cycle  ! no T=1 matrix elements
!				if(jj==0)print*,' TEST Y3 ', pn2bden(pnindx)%v(jj,jj,0) *sqrt(2*Jj+1.0)/sqrt(2*xJi+1.0)*trfactor*zetafactor*(-1)**JJ*phase*0.5
				isovec(jj) = isovec(jj) + (1./2.)*pn2bden(pnindx)%v(Jj,Jj,0)*sqrt(2*Jj+1.0)/sqrt(2*xJi+1.0)  *trfactor *zetafactor*phase
			end do
	   		do jj = jmin,jmax
	   			if(same .and. mod(jj,2)==0)cycle  ! no T=0 matrix elements
				isoscale(jj)=isoscale(jj)+(1./2.)*pn2bden(pnindx)%v(Jj,Jj,0)*sqrt(2*Jj+1.0)/sqrt(2*xJi+1.0)  *trfactor *zetafactor*phase
			end do
!............... SWAP   a -- > b
			phase = (-1)**( (ja+jb)/2)
            pnindx = (pnabcoupleR - PNcouples%meref(ipar) -1 )*pnparblocksize +  		 &
			       pncdcouple - PNcouples%meref(ipar) + PNcouples%mestart(ipar)

	   		do jj = jmin,jmax
	   			if(same .and. mod(jj,2)==1)cycle  ! no T=1 matrix elements
!				if(jj==0)print*,' TEST Y4 ', pn2bden(pnindx)%v(jj,jj,0) *sqrt(2*Jj+1.0)/sqrt(2*xJi+1.0) *trfactor*zetafactor*(-1)**JJ*phase*0.5

				isovec(jj) = isovec(jj) - (1./2.)*pn2bden(pnindx)%v(Jj,Jj,0)*sqrt(2*Jj+1.0)/sqrt(2*xJi+1.0)*trfactor*zetafactor*(-1)**JJ*phase
			end do
	   		do jj = jmin,jmax
	   			if(same .and. mod(jj,2)==0)cycle  ! no T=0 matrix elements
				isoscale(jj)=isoscale(jj)+(1./2.)*pn2bden(pnindx)%v(Jj,Jj,0)*sqrt(2*Jj+1.0)/sqrt(2*xJi+1.0)*trfactor *zetafactor*(-1)**JJ*phase
			end do
!............... finally, write out
			do jj = jmin,jmax
				if(abs(isovec(jj)) >  1.e-5)write(resultfile,'(x,4i4,2i6,2x,f12.5)')a,b,c,d,jj,1,isovec(jj)
				if(abs(isoscale(jj)) > 1.e-5)write(resultfile,'(x,4i4,2i6,2x,f12.5)')a,b,c,d,jj,0,isoscale(jj)
			end do

		end do  ! cd couple
	end do  ! ab couple

end subroutine isospin_expectations

end module twobodydens_util

max_stars_repo_name: noctildon/pyBigstick

max_issues_repo_name: noctildon/pyBigstick

max_forks_repo_name: noctildon/pyBigstick