!*****************************************************************************************
!>
!  For reading the JPL planetary and lunar ephemerides.
!  This is an extensively modified version of the original FORTRAN 77 code from JPL.
!
!### Ephemeris files
!  Note that this module uses the JPL binary ephemeris files, which
!  can be obtained using the instructions
!  [here](ftp://ssd.jpl.nasa.gov/pub/eph/planets/fortran/userguide.txt).
!  See also the comments in [[ephemeris_test]] for more details.
!
!### License
!
!  ***Original JPL License***
!
!  THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
!  CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
!  GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
!  ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
!  PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
!  TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
!  WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
!  PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
!  SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
!  SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
!
!  IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
!  BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
!  LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
!  INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
!  REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
!  REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
!
!  RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
!  THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
!  CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
!  ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
!
!  ***Modifications***
!
!  Modifications for the Fortran Astrodynamics Toolkit are covered
!  under the [following license](https://github.com/jacobwilliams/Fortran-Astrodynamics-Toolkit/blob/master/LICENSE).
!
!### History
!  * [Original code from JPL](ftp://ssd.jpl.nasa.gov/pub/eph/planets/fortran/) Version : March 25, 2013
!  * Extensive modifications by Jacob Williams for the Fortran Astrodynamics Toolkit.

    module jpl_ephemeris_module

    use, intrinsic :: iso_fortran_env, only: real64,error_unit
    use ephemeris_module

    implicit none

    private

    integer,parameter :: wp = real64  !! using double precision reals since that is how the ephemeris is stored.

    integer,parameter,public :: nmax   = 1000   !! Current maximum number of ephemeris constants
                                                !! used in the integration and listed
                                                !! in the `header.xxx` file for the ephemeris.
    integer,parameter        :: oldmax = 400    !! For earlier versions of the code, the maximum
                                                !! number of ephemeris constants used in the
                                                !! integration and listed in the `header.xxx`
                                                !! file for the ephemeris.

    integer,parameter :: nrecl = 4  !! `nrecl=1` if `recl` in the open statement is the record length in s.p. words
                                    !! `nrecl=4` if `recl` in the open statement is the record length in bytes

    type,extends(ephemeris_class),public :: jpl_ephemeris

        !! Main class for accessing a JPL ephemeris file.

        character(len=:),allocatable :: namfil  !! name of the binary ephemeris file

        integer :: ksize = 2036  !! ksize must be set by the user according to the ephemeris to be read:
                                 !!   for ***de200***, set `ksize=1652`,
                                 !!   for ***de405***, set `ksize=2036`,
                                 !!   for ***de406***, set `ksize=1456`,
                                 !!   for ***de414***, set `ksize=2036`,
                                 !!   for ***de418***, set `ksize=2036`,
                                 !!   for ***de421***, set `ksize=2036`,
                                 !!   for ***de422***, set `ksize=2036`,
                                 !!   for ***de423***, set `ksize=2036`,
                                 !!   for ***de424***, set `ksize=2036`,
                                 !!   for ***de430***, set `ksize=2036`.

        !ephhdr
        integer,dimension(3,13)  :: ipt   = 0       !! ipt(39)
        real(wp),dimension(nmax) :: cval  = 0.0_wp
        real(wp),dimension(3)    :: ss    = 0.0_wp
        real(wp)                 :: au    = 0.0_wp
        real(wp)                 :: emrat = 0.0_wp
        integer                  :: ncon  = 0
        integer                  :: numde = 0

        !stcomx
        real(wp),dimension(6) :: pvsun = 0.0_wp   !! dp 6-word array containing the barycentric position and
                                                  !! velocity of the sun.
        logical               :: km    = .true.   !! logical flag defining physical units of the output states.
                                                  !!   km = .true.  : km and km/sec
                                                  !!   km = .false. : au and au/day
                                                  !! for nutations and librations.  angle unit is always radians.
        logical               :: bary  = .false.  !! logical flag defining output center.
                                                  !! only the 9 planets are affected.
                                                  !!   bary = .true.  : center is solar-system barycenter
                                                  !!   bary = .false. : center is sun

        !chrhdr
        character(len=6),dimension(14,3) :: ttl  = ''
        character(len=6),dimension(nmax) :: cnam = ''

        logical :: initialized = .false.  !! is the ephemeris initialized?
        integer :: nrfile      = 0        !! file unit for the ephemeris file
        integer :: nrl         = -1       !! this was formerly in state
        integer :: ncoeffs     = 0

        ! formerly in state:
        real(wp),dimension(1500) :: buf = 0.0_wp

        ! formerly in interp:
        real(wp),dimension(18) :: pc = 0.0_wp
        real(wp),dimension(18) :: vc = 0.0_wp
        integer                :: np = 2
        integer                :: nv = 3
        real(wp)               :: twot = 0.0_wp

    contains

        procedure,public :: get_rv => get_rv_from_jpl_ephemeris

        procedure,public :: initialize => initialize_ephemeris
        procedure,public :: get_state
        procedure,public :: get_constants
        procedure,public :: close => close_ephemeris

        procedure :: interp
        procedure :: state

    end type jpl_ephemeris

    !>
    !  These correspond to the numbering convention for 'ntarg' and 'ncent' in [[get_state]]:
    character(len=*),dimension(15),parameter :: list_of_bodies = [ &
                                                        'mercury                ',&
                                                        'venus                  ',&
                                                        'earth                  ',&
                                                        'mars                   ',&
                                                        'jupiter                ',&
                                                        'saturn                 ',&
                                                        'uranus                 ',&
                                                        'neptune                ',&
                                                        'pluto                  ',&
                                                        'moon                   ',&
                                                        'sun                    ',&
                                                        'solar-system barycenter',&
                                                        'earth-moon barycenter  ',&
                                                        'nutations              ',&
                                                        'librations             ']

    !public routines:
    public :: ephemeris_test

    contains
!*****************************************************************************************

!*****************************************************************************************
!>
!  Interface for the [[ephemeris_module]].

    subroutine get_rv_from_jpl_ephemeris(me,et,targ,obs,rv,status_ok)

    use time_module,           only: et_to_jd
    use conversion_module,     only: day2sec
    use celestial_body_module, only: celestial_body
    use numbers_module,        only: zero

    implicit none

    class(jpl_ephemeris),intent(inout) :: me
    real(wp),intent(in)                :: et         !! ephemeris time [sec]
    type(celestial_body),intent(in)    :: targ       !! target body
    type(celestial_body),intent(in)    :: obs        !! observer body
    real(wp),dimension(6),intent(out)  :: rv         !! state of targ w.r.t. obs [km,km/s] in ICRF frame
    logical,intent(out)                :: status_ok  !! true if there were no problems

    real(wp) :: jd     !! julian date for input to [[get_state]].
    integer :: ntarg   !! id code for target body
    integer :: ncent   !! id code for observer body

    if (targ==obs) then
        !don't bother if target and observer are the same body
        rv = zero
        status_ok = .true.
    else

        !convert to expected inputs:
        jd    = et_to_jd(et)
        ntarg = spice_id_to_old_id(targ%id)
        ncent = spice_id_to_old_id(obs%id)

        if (ntarg>0 .and. ncent>0) then
            call me%get_state(jd,ntarg,ncent,rv,status_ok)
            if (status_ok) then
                if (.not. me%km) then
                    !we must return in units of km/s
                    !so, convert from AU, AU/day to km, km/s
                    rv      = rv * me%au
                    rv(4:6) = rv(4:6) / day2sec
                end if
            else
                write(error_unit,'(A)') 'Error in get_rv_from_jpl_ephemeris: '//&
                            'Error calling ephemeris.'
            end if
        else
            write(error_unit,'(A)') 'Error in get_rv_from_jpl_ephemeris: '//&
                        'No ephemeris for this body.'
            status_ok = .false.
        end if

    end if

    end subroutine get_rv_from_jpl_ephemeris
!*****************************************************************************************

!*****************************************************************************************
!> author: Jacob Williams
!  date: 3/20/2016
!
!  Convert the NAIF SPICE ID code to the old one used by the JPL ephemeris.
!  Returns `0` if the body was not found.

    pure function spice_id_to_old_id(spice_id) result(old_id)

    implicit none

    integer,intent(in) :: spice_id !! the ID code used by SPICE
    integer            :: old_id   !! the ID code used by this module (old JPL ephemeris code)

    integer :: i !! counter

    !>
    !  The index of this array is the old ID code. The value is the new code.
    !  See: [NAIF Integer ID codes](http://naif.jpl.nasa.gov/pub/naif/toolkit_docs/FORTRAN/req/naif_ids.html)
    integer,parameter,dimension(13) :: new_ids = &
        [   199,&  ! mercury
            299,&  ! venus
            399,&  ! earth
            499,&  ! mars
            599,&  ! jupiter
            699,&  ! saturn
            799,&  ! uranus
            899,&  ! neptune
            999,&  ! pluto
            301,&  ! moon
            10, &  ! sun
            0,  &  ! solar-system barycenter
            3   ]  ! earth-moon barycenter

    !just a simple search of the list:
    ! [small enough that bisection search probably not worth it]
    do i=1,size(new_ids)
        if (new_ids(i)==spice_id) then
            old_id = i
            return
        end if
    end do

    !not found:
    old_id = 0

    end function spice_id_to_old_id
!*****************************************************************************************

!*****************************************************************************************
!>
!  This subroutine reads the JPL planetary ephemeris
!  and gives the position and velocity of the point `ntarg`
!  with respect to `ncent`.
!
!### Notes
!
!  The numbering convention for `ntarg` and `ncent` is:
!
!    1 = mercury      8 = neptune
!    2 = venus        9 = pluto
!    3 = earth       10 = moon
!    4 = mars        11 = sun
!    5 = jupiter     12 = solar-system barycenter
!    6 = saturn      13 = earth-moon barycenter
!    7 = uranus      14 = nutations (longitude and obliq)
!                    15 = librations, if on eph file
!
!  (if nutations are wanted, set ntarg = 14.
!  for librations, set ntarg = 15. set ncent=0.)

    subroutine get_state(me,jd,ntarg,ncent,rrd,status_ok)

    implicit none

    class(jpl_ephemeris),intent(inout) :: me
    real(wp),intent(in)                :: jd         !! d.p. Julian ephemeris date at which interpolation is wanted.
    integer,intent(in)                 :: ntarg      !! integer number of 'target' point.
    integer,intent(in)                 :: ncent      !! integer number of 'center' point.
    real(wp),dimension(6),intent(out)  :: rrd        !! output 6-word d.p. array containing position and velocity
                                                     !! of point `ntarg` relative to `ncent`.
                                                     !! the units are AU and AU/day (or km and km/sec if `me%km=.true.`).
                                                     !! For librations the units are radians and radians
                                                     !! per day. In the case of nutations the first four words of
                                                     !! `rrd` will be set to nutations and rates, having units of
                                                     !! radians and radians/day.
    logical,intent(out)                :: status_ok  !! true if there were no problems

    real(wp),dimension(2)    :: et2
    real(wp),dimension(6,13) :: pv
    real(wp),dimension(6,11) :: pvst
    real(wp),dimension(4)    :: pnut
    integer,dimension(12)    :: list
    integer                  :: i,j,k
    logical                  :: bsave

    status_ok = .false.

    if (me%initialized) then

        ! initialize et2 for 'state' and set up component count
        et2(1) = jd
        et2(2) = 0.0_wp
        list   = 0
        rrd    = 0.0_wp

        if (ntarg /= ncent) then

            select case (ntarg)

            case (14)    !nutation

                if (me%ipt(2,12)>0) then            !me%ipt(35)
                    list(11) = 2
                    call me%state(et2,list,pvst,pnut,status_ok)
                    if (.not. status_ok) then
                        write(error_unit,'(A)') 'error interpolating nutations in get_state.'
                        return
                    end if
                    rrd(1:4) = pnut(1:4)
                    rrd(5) = 0.0_wp
                    rrd(6) = 0.0_wp
                    return
                else
                    rrd(1:4) = 0.0_wp
                    write(error_unit,'(A)') 'error in get_state: the ephemeris file does not contain nutations.'
                    return
                endif

            case (15)    !librations

                if (me%ipt(2,13)>0) then            !me%ipt(38)
                    list(12) = 2
                    call me%state(et2,list,pvst,pnut,status_ok)
                    if (.not. status_ok) then
                        write(error_unit,'(A)') 'error interpolating librations in get_state.'
                        return
                    end if
                    rrd = pvst(:,11)
                    return
                else
                    write(error_unit,'(A)') 'error in get_state: the ephemeris file does not contain librations.'
                    return
                endif

               case default

                ! force barycentric output by 'state'

                bsave = me%bary
                me%bary  = .true.

                ! set up proper entries in 'list' array for state call

                do i=1,2
                    k=ntarg
                    if (i == 2)  k = ncent
                    if (k <= 10) list(k)  = 2
                    if (k == 10) list(3)  = 2
                    if (k == 3)  list(10) = 2
                    if (k == 13) list(3)  = 2
                enddo

                ! make call to state

                call me%state(et2,list,pvst,pnut,status_ok)
                if (.not. status_ok) then
                    write(error_unit,'(A)') 'error interpolating state in get_state.'
                    return
                end if
                do i=1,10
                    do j = 1,6
                        pv(j,i) = pvst(j,i)
                    end do
                enddo

                if (ntarg == 11 .or. ncent == 11) pv(:,11) = me%pvsun
                if (ntarg == 12 .or. ncent == 12) pv(:,12) = 0.0_wp
                if (ntarg == 13 .or. ncent == 13) pv(:,13) = pvst(:,3)

                if (ntarg*ncent == 30 .and. ntarg+ncent == 13) then
                    pv(:,3) = 0.0_wp
                else
                    if (list(3) == 2)  pv(:,3)  = pvst(:,3) - pvst(:,10)/(1.0_wp+me%emrat)
                    if (list(10) == 2) pv(:,10) = pv(:,3) + pvst(:,10)
                end if

                rrd = pv(:,ntarg) - pv(:,ncent)

                me%bary = bsave

            end select

        end if

        status_ok = .true.

    else
        write(error_unit,'(A)') 'error in get_state: the ephemeris is not initialized.'
    end if

    end subroutine get_state
!*****************************************************************************************

!*****************************************************************************************
!>
!  this subroutine differentiates and interpolates a
!  set of chebyshev coefficients to give position and velocity.

    subroutine interp(me,buf,t,ncf,ncm,na,ifl,pv)

    implicit none

    class(jpl_ephemeris),intent(inout) :: me
    integer,intent(in)                 :: ncf  !! # of coefficients per component
    integer,intent(in)                 :: ncm  !! # of components per set of coefficients
    real(wp),dimension(ncf,ncm,*)      :: buf  !! 1st location of array of d.p. chebyshev coefficients of position
    real(wp),dimension(2),intent(in)   :: t
    integer,intent(in)                 :: na   !! # of sets of coefficients in full array
                                               !! (i.e., # of sub-intervals in full interval)
    integer,intent(in)                 :: ifl  !! integer flag
                                               !! = 1 for positions only
                                               !! = 2 for pos and vel
    real(wp),dimension(ncm,*)          :: pv   !! interpolated quantities requested.  dimension
                                               !! expected is pv(ncm,ifl), dp.

    real(wp) :: dna,dt1,temp,vfac,tc
    integer :: l,i,j

    ! entry point. get correct sub-interval number for this set
    ! of coefficients and then get normalized chebyshev time
    ! within that subinterval.

    dna  = dble(na)
    dt1  = int(t(1))
    temp = dna*t(1)
    l    = int(temp-dt1)+1

    ! tc is the normalized chebyshev time (-1 <= tc <= 1)

    tc=2.0_wp*(mod(temp,1.0_wp)+dt1)-1.0_wp

    ! check to see whether chebyshev time has changed,
    ! and compute new polynomial values if it has.
    ! (the element pc(2) is the value of t1(tc) and hence
    ! contains the value of tc on the previous call.)

    if (tc/=me%pc(2)) then
        me%np    = 2
        me%nv    = 3
        me%pc(2) = tc
        me%twot  = tc+tc
    endif

    ! be sure that at least 'ncf' polynomials have been evaluated
    ! and are stored in the array 'pc'.

    if (me%np<ncf) then
        do i=me%np+1,ncf
            me%pc(i)=me%twot*me%pc(i-1)-me%pc(i-2)
        end do
        me%np=ncf
    endif

    ! interpolate to get position for each component

    do i=1,ncm
        pv(i,1)=0.0_wp
        do j=ncf,1,-1
            pv(i,1)=pv(i,1)+me%pc(j)*buf(j,i,l)
        end do
    end do
    if (ifl<=1) return

    ! if velocity interpolation is wanted, be sure enough
    ! derivative polynomials have been generated and stored.

    vfac=(dna+dna)/t(2)
    me%vc(3)=me%twot+me%twot
    if (me%nv<ncf) then
        do i=me%nv+1,ncf
            me%vc(i)=me%twot*me%vc(i-1)+me%pc(i-1)+me%pc(i-1)-me%vc(i-2)
        end do
        me%nv=ncf
    endif

    ! interpolate to get velocity for each component

    do i=1,ncm
        pv(i,2)=0.0_wp
        do j=ncf,2,-1
            pv(i,2)=pv(i,2)+me%vc(j)*buf(j,i,l)
        end do
        pv(i,2)=pv(i,2)*vfac
    end do

    end subroutine interp
!*****************************************************************************************

!*****************************************************************************************
!>
!  this subroutine breaks a d.p. number into a d.p. integer
!  and a d.p. fractional part.

    subroutine split(tt,fr)

    implicit none

    real(wp),intent(in)               :: tt  !! d.p. input number
    real(wp),dimension(2),intent(out) :: fr  !! d.p. 2-word output array.
                                             !! fr(1) contains integer part
                                             !! fr(2) contains fractional part
                                             !! for negative input numbers, fr(1) contains the next
                                             !! more negative integer; fr(2) contains a positive fraction.

    ! get integer and fractional parts

    fr(1) = int(tt)
    fr(2) = tt-fr(1)

    if (tt>=0.0_wp .or. fr(2)==0.0_wp) return

    ! make adjustments for negative input number

    fr(1) = fr(1) - 1.0_wp
    fr(2) = fr(2) + 1.0_wp

    end subroutine split
!*****************************************************************************************

!*****************************************************************************************
!>
!  Initialize the ephemeris.
!  This routine may be called to load a different ephemeris file.
!  Otherwise, it is called on the first call to get_state, and loads
!  the file specified in the module header.
!
!# Note
!  * Based on code formerly in [[state]].

    subroutine initialize_ephemeris(me,filename,ksize,km,bary,status_ok)

    implicit none

    class(jpl_ephemeris),intent(inout) :: me
    character(len=*),intent(in)        :: filename   !! ephemeris file name
    integer,intent(in),optional        :: ksize      !! corresponding `ksize`
    logical,intent(in),optional        :: km         !! defining physical units of the output states.
                                                     !! `km = .true.`  : km and km/sec [default],
                                                     !! `km = .false.` : au and au/day.
    logical,intent(in),optional        :: bary       !! logical flag defining output center.
                                                     !! only the 9 planets are affected.
                                                     !! `bary = .true.`  : center is solar-system barycenter,
                                                     !! `bary = .false.` : center is sun [default].
    logical,intent(out) :: status_ok                 !! true if there were not problems.

    !local variables:
    integer :: irecsz,istat,i,j,k,l

    !just in case it was already open:
    call me%close()  ! clears everything in the class

    !ephemeris file name:
    me%namfil = trim(filename)

    !optional inputs:
    if (present(ksize)) me%ksize  = ksize
    if (present(km))    me%km     = km
    if (present(bary))  me%bary   = bary

    irecsz     = nrecl*me%ksize
    me%ncoeffs = me%ksize/2

    open(newunit = me%nrfile,     &
         file    = me%namfil,     &
         access  = 'DIRECT',      &
         form    = 'UNFORMATTED', &
         action  = 'READ',        &   !JW added
         recl    = irecsz,        &
         iostat  = istat,         &
         status  = 'OLD'          )

         !write(*,*) "istat=",istat

    status_ok = (istat==0)  !if there were no problems opening the file

    if (status_ok) then

        read(me%nrfile,rec=1,iostat=istat) &
                 me%ttl,(me%cnam(k),k=1,oldmax),me%ss,me%ncon,me%au,me%emrat,&
                 ((me%ipt(i,j),i=1,3),j=1,12),me%numde,(me%ipt(i,13),i=1,3), &
                 (me%cnam(l),l=oldmax+1,me%ncon)

        if (istat==0) then
            if (me%ncon <= oldmax) then
                read(me%nrfile,rec=2,iostat=istat) (me%cval(i),i=1,oldmax)
            else
                read(me%nrfile,rec=2,iostat=istat) (me%cval(i),i=1,me%ncon)
            endif
            if (istat==0) then
                me%nrl = 0
                me%initialized = .true.
            end if
        end if

        ! check if the reads went OK:
        status_ok = me%initialized
        if (status_ok) then
            !initialize some of the class variables:
            ! [note: this was formerly done in the interp routine]
            me%pc(1) = 1.0_wp
            me%vc(2) = 1.0_wp
        else
            write(error_unit,'(A)') 'Error reading ephemeris file: '//trim(me%namfil)
        end if

    else
        write(error_unit,'(A)') 'Error opening ephemeris file: '//trim(me%namfil)
    end if

    end subroutine initialize_ephemeris
!*****************************************************************************************

!*****************************************************************************************
!>
!  Close the ephemeris.

    subroutine close_ephemeris(me)

    implicit none

    class(jpl_ephemeris),intent(inout) :: me

    logical :: is_open
    integer :: istat

    if (me%initialized) then

        !close the file:
        inquire(unit=me%nrfile,opened=is_open,iostat=istat)
        if (is_open) close(unit=me%nrfile,iostat=istat)

        !initialize all class variables to defaults:
        call clear(me)

    end if

    contains

        subroutine clear(eph)
            !! clear all the variables in the [[jpl_ephemeris]] class.
            implicit none
            class(jpl_ephemeris),intent(out) :: eph
        end subroutine clear

    end subroutine close_ephemeris
!*****************************************************************************************

!*****************************************************************************************
!>
!  This subroutine reads and interpolates the JPL planetary ephemeris file.
!
!@note The ephemeris is assumed to have been initialized.

    subroutine state(me,et2,list,pv,pnut,status_ok)

    implicit none

    class(jpl_ephemeris),intent(inout)   :: me
    real(wp),dimension(2),intent(in)     :: et2     !! 2-word Julian ephemeris epoch at which interpolation
                                                    !! is wanted.  any combination of `et2(1)+et2(2)` which falls
                                                    !! within the time span on the file is a permissible epoch.
                                                    !!
                                                    !! ***a.*** for ease in programming, the user may put the
                                                    !!    entire epoch in `et2(1)` and set `et2(2)=0`.
                                                    !! ***b.*** for maximum interpolation accuracy, set `et2(1)` =
                                                    !!    the most recent midnight at or before interpolation
                                                    !!    epoch and set `et2(2)` = fractional part of a day
                                                    !!    elapsed between `et2(1)` and epoch.
                                                    !! ***c.*** as an alternative, it may prove convenient to set
                                                    !!    `et2(1)` = some fixed epoch, such as start of integration,
                                                    !!    and `et2(2)` = elapsed interval between then and epoch.

    integer,dimension(12),intent(in)     :: list    !! 12-word integer array specifying what interpolation
                                                    !! is wanted for each of the bodies on the file.
                                                    !!    ***list(i) = 0*** : no interpolation for body i,
                                                    !!    ***list(i) = 1*** : position only,
                                                    !!    ***list(i) = 2*** : position and velocity.
                                                    !!
                                                    !! The designation of the astronomical bodies by i is:
                                                    !!    ***i = 1*** : mercury,
                                                    !!    ***i = 2*** : venus,
                                                    !!    ***i = 3*** : earth-moon barycenter,
                                                    !!    ***i = 4*** : mars,
                                                    !!    ***i = 5*** : jupiter,
                                                    !!    ***i = 6*** : saturn,
                                                    !!    ***i = 7*** : uranus,
                                                    !!    ***i = 8*** : neptune,
                                                    !!    ***i = 9*** : pluto,
                                                    !!    ***i =10*** : geocentric moon,
                                                    !!    ***i =11*** : nutations in longitude and obliquity,
                                                    !!    ***i =12*** : lunar librations (if on file).

    real(wp),dimension(6,11),intent(out) :: pv      !! dp 6 x 11 array that will contain requested interpolated
                                                    !! quantities (other than nutation, stored in `pnut`).
                                                    !! the body specified by `list(i)` will have its
                                                    !! state in the array starting at `pv(1,i)`.
                                                    !! (on any given call, only those words in `pv` which are
                                                    !! affected by the first 10 `list` entries, and by `list(12)`
                                                    !! if librations are on the file, are set.
                                                    !! the rest of the `pv` array is untouched.)
                                                    !! the order of components starting in `pv(1,i)` is:
                                                    !! `x`,`y`,`z`,`dx`,`dy`,`dz`.
                                                    !!
                                                    !! All output vectors are referenced to the earth mean
                                                    !! equator and equinox of J2000 if the DE number is 200 or
                                                    !! greater; of B1950 if the DE number is less than 200.
                                                    !!
                                                    !! The moon state is always geocentric; the other nine states
                                                    !! are either heliocentric or solar-system barycentric,
                                                    !! depending on the setting of the `bary` variable in the class.
                                                    !!
                                                    !! Lunar librations, if on file, are put into `pv(k,11)` if
                                                    !! `list(12)` is `1` or `2`.

    real(wp),dimension(4),intent(out)    :: pnut    !! dp 4-word array that will contain nutations and rates,
                                                    !! depending on the setting of `list(11)`.  the order of
                                                    !! quantities in `pnut` is:
                                                    !!
                                                    !! * `d psi`  (nutation in longitude),
                                                    !! * `d epsilon` (nutation in obliquity),
                                                    !! * `d psi dot`,
                                                    !! * `d epsilon dot`.
    logical,intent(out) :: status_ok !! true if there were no problems

    real(wp),dimension(2)    :: t
    real(wp),dimension(4)    :: pjd
    real(wp) :: aufac,s,tmp1,tmp2
    integer :: istat,i,j,k,nr

    status_ok = .true.

    s = et2(1)-0.5_wp
    call split(s,pjd(1:2))
    call split(et2(2),pjd(3:4))
    pjd(1) = pjd(1)+pjd(3)+0.5_wp
    pjd(2) = pjd(2)+pjd(4)
    call split(pjd(2),pjd(3:4))
    pjd(1) = pjd(1)+pjd(3)

    ! error return for epoch out of range

    if (pjd(1)+pjd(4)<me%ss(1) .or. pjd(1)+pjd(4)>me%ss(2)) then
        write(error_unit,'(A,F12.2,A,2F12.2)') &
                'Error: requested jed,',&
                et2(1)+et2(2),&
                ' not within ephemeris limits,',&
                me%ss(1),me%ss(2)
        status_ok = .false.
        return
    end if

    ! calculate record # and relative time in interval

    nr = int((pjd(1)-me%ss(1))/me%ss(3))+3
    if (pjd(1)==me%ss(2)) nr = nr-1

    tmp1 = dble(nr-3)*me%ss(3) + me%ss(1)
    tmp2 = pjd(1) - tmp1
    t(1) = (tmp2 + pjd(4))/me%ss(3)

    ! read correct record if not in core

    if (nr/=me%nrl) then
        me%nrl = nr
        read(me%nrfile,rec=nr,iostat=istat) (me%buf(k),k=1,me%ncoeffs)
        if (istat/=0) then
            write(error_unit,'(2F12.2,A80)') et2,'Error return in state'
            status_ok = .false.
            return
        end if
    endif

    if (me%km) then
        t(2) = me%ss(3)*86400.0_wp
        aufac = 1.0_wp
    else
        t(2) = me%ss(3)
        aufac = 1.0_wp/me%au
    endif

    ! interpolate ssbary sun

    call me%interp(me%buf(me%ipt(1,11)),t,me%ipt(2,11),3,me%ipt(3,11),2,me%pvsun)

    me%pvsun = me%pvsun * aufac

    ! check and interpolate whichever bodies are requested

    do i=1,10
        if (list(i)==0) cycle
        call me%interp(me%buf(me%ipt(1,i)),t,me%ipt(2,i),3,me%ipt(3,i),list(i),pv(1,i))
        do j=1,6
            if (i<=9 .and. .not.me%bary) then
                pv(j,i) = pv(j,i)*aufac - me%pvsun(j)
            else
                pv(j,i) = pv(j,i)*aufac
            endif
        enddo
    end do

    ! do nutations if requested (and if on file)

    if (list(11)>0 .and. me%ipt(2,12)>0) &
        call me%interp(me%buf(me%ipt(1,12)),t,me%ipt(2,12),2,me%ipt(3,12),list(11),pnut)

    ! get librations if requested (and if on file)

    if (list(12)>0 .and. me%ipt(2,13)>0) &
        call me%interp(me%buf(me%ipt(1,13)),t,me%ipt(2,13),3,me%ipt(3,13),list(12),pv(1,11))

    end subroutine state
!*****************************************************************************************

!*****************************************************************************************
!>
!  Obtain the constants from the ephemeris file.

    subroutine get_constants(me,nam,val,sss,n)

    implicit none

    class(jpl_ephemeris),intent(inout)        :: me
    character(len=6),dimension(:),intent(out) :: nam  !! array of constant names
    real(wp),dimension(:),intent(out)         :: val  !! array of values of constants
    real(wp),dimension(3),intent(out)         :: sss  !! jd start, jd stop, step of ephemeris
    integer,intent(out)                       :: n    !! number of entries in `nam` and `val` arrays

    integer :: i

    if (me%initialized) then

        n   = me%ncon
        sss = me%ss

        do i=1,n
            nam(i) = me%cnam(i)
            val(i) = me%cval(i)
        enddo

    else
        write(error_unit,'(A)') 'error in get_constants: the ephemeris is not initialized.'
    end if

    end subroutine get_constants
!*****************************************************************************************

!*****************************************************************************************
!>
!  Ephemeris test routine.
!
!### Ephemeris files
!  This routine requires the DE405 and DE421 JPL binary ephemeris files
!  to be present in the `../eph` directory.
!  These can be built by using the instructions
!  [here](ftp://ssd.jpl.nasa.gov/pub/eph/planets/fortran/userguide.txt).
!
!  For example (on Linux):
!```bash
!  wget ftp://ssd.jpl.nasa.gov/pub/eph/planets/fortran/*
!  wget ftp://ssd.jpl.nasa.gov/pub/eph/planets/ascii/de405/*
!  #edit asc2eph.f file to set NRECL = 4:
!  sed -i '_original' '/^C.*PARAMETER ( NRECL = 4 )/s/^C//' asc2eph.f
!  gfortran asc2eph.f -o asc2eph
!  cat header.405 ascp*.405 | ./asc2eph
!  mkdir Fortran-Astrodynamics-Toolkit/eph
!  mv JPLEPH Fortran-Astrodynamics-Toolkit/eph/JPLEPH.405
!```

    subroutine ephemeris_test()

    implicit none

    character(len=6),dimension(nmax) :: nams
    real(wp) :: jd
    real(wp),dimension(6) :: rv,rv1,rv2,diffrv
    real(wp),dimension(3) :: ss
    real(wp),dimension(nmax) :: vals
    integer :: nvs,ntarg,nctr,i,j
    type(jpl_ephemeris) :: eph405, eph421
    logical :: status_ok_405,status_ok_421

    character(len=*),parameter :: ephemeris_file_405 = '../eph/JPLEPH.405' !! JPL DE405 ephemeris file
    character(len=*),parameter :: ephemeris_file_421 = '../eph/JPLEPH.421' !! JPL DE421 ephemeris file

    write(*,*) ''
    write(*,*) '---------------'
    write(*,*) ' ephemeris_test'
    write(*,*) '---------------'
    write(*,*) ''

    !initialize:
    call eph405%initialize(filename=ephemeris_file_405,status_ok=status_ok_405)
    call eph421%initialize(filename=ephemeris_file_421,status_ok=status_ok_421)

    if (status_ok_405) then

        !get some constants from the file:
        call eph405%get_constants(nams,vals,ss,nvs)

        write(*,'(A)') ''
        write(*,'(A)') 'Ephemeris initialized'
        write(*,'(A,1X,F15.3,1X,A,1X,F15.3)') 'JD range: ',ss(1), 'to ', ss(2)
        write(*,'(A)') ''
        do i=1,nvs
            write(*,'(A,1X,D25.16)') nams(i), vals(i)
        end do

        jd = 2451536.5d0       ! julian date

        if (jd < ss(1) .or. jd > ss(2)) then

            write(*,'(A)') ''
            write(*,*) 'error: jed out of bounds.'
            write(*,*) 'jed   = ', jd
            write(*,*) 'ss(1) = ', ss(1)
            write(*,*) 'ss(2) = ', ss(2)

        else

            !test DE405:
            do j=1,2

                if (j==1) then
                    ntarg = 3      !earth
                    nctr  = 11     !sun
                else
                    ntarg = 10     !moon
                    nctr  = 3      !earth
                end if

                write(*,*) ''
                write(*,*) 'DE405'
                write(*,*) 'state of "'//trim(list_of_bodies(ntarg))//&
                            '" wrt "'//trim(list_of_bodies(nctr))//'"'

                do i=1,10
                    call eph405%get_state( jd, ntarg, nctr, rv, status_ok_405)
                    write(*,'(F15.2,1X,*(E25.16,1X))') jd, norm2(rv(1:3)), rv
                    jd = jd + 10.0_wp
                end do

            end do
        end if

    else
        write(*,*) 'Error opening DE405 ephemeris file'
    end if

    if (status_ok_405 .and. status_ok_421) then

        !compare DE405 with DE421
        do j=1,2

            if (j==1) then
                ntarg = 3      !earth
                nctr  = 11     !sun
            else
                ntarg = 10     !moon
                nctr  = 3      !earth
            end if

            write(*,*) ''
            write(*,*) 'DE421 - DE405 Difference'
            write(*,*) 'state of "'//trim(list_of_bodies(ntarg))//&
                        '" wrt "'//trim(list_of_bodies(nctr))//'"'

            do i=1,10
                call eph405%get_state( jd, ntarg, nctr, rv1, status_ok_405)
                call eph421%get_state( jd, ntarg, nctr, rv2, status_ok_421)
                diffrv = rv2 - rv1
                write(*,'(F15.2,1X,*(E25.16,1X))') jd, norm2(diffrv(1:3)), norm2(diffrv(4:6))
                jd = jd + 10.0_wp
            end do

        end do

    else
        write(*,*) 'Error opening DE421 ephemeris file'
    end if

    !cleanup:
    call eph405%close()
    call eph421%close()

    end subroutine ephemeris_test
!*****************************************************************************************

!*****************************************************************************************
    end module jpl_ephemeris_module
!*****************************************************************************************


text:
  MODULE Parabolic_constants

  use clubb_precision, only: &
        dp ! double precision

  IMPLICIT NONE
  public ! Default
  INTEGER, PARAMETER  :: r8 = dp
! INCLUDE 'constants'
    ! Begin file constants
    REAL(r8), PARAMETER :: lnsqrttwopi=0.9189385332046727418_r8     ! log(sqrt(2*pi))
    REAL(r8), PARAMETER :: sqrttwopi=2.5066282746310005024_r8       ! sqrt(2*pi)
    REAL(r8), PARAMETER :: eulmasc=0.5772156649015328606_r8         ! Euler-Mascheroni constant
    REAL(r8), PARAMETER :: pihalf=1.5707963267948966192_r8          ! pi/2
    REAL(r8), PARAMETER :: pikwart=0.7853981633974483096_r8         ! pi/4
    REAL(r8), PARAMETER :: pidriekwart=2.3561944901923449288_r8     ! 3pi/4
    REAL(r8), PARAMETER :: sqrtpi=1.7724538509055160272_r8          ! sqrt(pi)
    REAL(r8), PARAMETER :: lnpi=1.1447298858494001741_r8            ! log(pi)
    REAL(r8), PARAMETER :: pi=3.1415926535897932385_r8              ! pi
    REAL(r8), PARAMETER :: onepi=.31830988618379067153_r8           ! 1/pi
    REAL(r8), PARAMETER :: sqrt2=1.4142135623730950488_r8           ! sqrt(2)
    REAL(r8), PARAMETER :: sqrt3=1.7320508075688772935_r8           ! sqrt(3)
    REAL(r8), PARAMETER :: sqrt2opi=0.7978845608028653559_r8        ! sqrt(2/pi)
    REAL(r8), PARAMETER :: twopi=6.2831853071795864769_r8           ! 2*pi
    REAL(r8), PARAMETER :: oneoversqrtpi=0.5641895835477562869_r8   ! 1/sqrt(pi)
    REAL(r8), PARAMETER :: twooversqrtpi=1.1283791670955125739_r8   ! 2/sqrt(pi)
    REAL(r8), PARAMETER :: onethird=0.33333333333333333_r8          ! 1/3
    REAL(r8), PARAMETER :: twothird=0.6666666666666666666667_r8     ! 2/3
    REAL(r8), PARAMETER :: onesix=.166666666666666666666666666667_r8! 1/6
    REAL(r8), PARAMETER :: piquart=0.78539816339744830962_r8        ! pi/4
    REAL(r8), PARAMETER :: twoexp14=1.18920711500272106671749997_r8 ! 2**(1/4)
    REAL(r8), PARAMETER :: dwarf=TINY(0.0_r8)*1000.0_r8             ! safe underflow limit
    REAL(r8), PARAMETER :: giant=HUGE(0.0_r8)/1000.0_r8             ! safe overflow limit
    REAL(r8), PARAMETER :: mac_tol=EPSILON(0.0_r8)                  ! machine-epsilon
    REAL(r8) :: epss                                                ! demanded accuracy
  ! End File constants
  END MODULE Parabolic_constants

max_stars_repo_licenses: Intel Unlicense NetCDF

max_issues_repo_licenses: Intel Unlicense NetCDF

max_forks_repo_licenses: Intel Unlicense NetCDF


!---------------------------------------------------------------------!
! OWNER: Ithaca Combustion Enterprise, LLC                            !
! COPYRIGHT: © 2012, Ithaca Combustion Enterprise, LLC                !
! LICENSE: BSD 3-Clause License (The complete text of the license can !
!  be found in the `LICENSE-ICE.txt' file included in the ISAT-CK7    !
!  source directory.)                                                 !
!---------------------------------------------------------------------!

module sell_bt

!  Module for searching a set of ellipsoids (SELLs) using a binary tree (BT).

!  This code is based on sell_ebt, where nodes are referred to as BEs.
!  Here, to avoid name conflicts, BE is changed to ND (for node).

!  A BT has nodes and leaves, with an ND associated with each.
!  The ND associated with a leaf is an ellipsoid (ELL) in the SELL.
!  Each node has two children (1 and 2).

!  S.B. Pope  6/3/2006

!  Data structures
!     bt_type    - for BTs
!     nd_type    - for bounding ellipsoids (ND)
!     nd_pointer - ponter to ND

!================================================

! Public soubroutines

! bt_initialize  ! allocate and initialize BT
! bt_destroy     ! destroy and deallocate BT
! bt_add         ! add leaf ELL to BT
! bt_remove      ! remove leaf ELL from BT
! bt_purge       ! remove all nodes (NDs) from BT
! bt_rebuild     ! rebuild BT
! bt_traverse    ! traverse BT based on query point x
! bt_ball_pts    ! find points in ball

! bt_param_set   ! set control parameters
! bt_param_get   ! get control parameters
! bt_prop        ! get BT properties
! bt_write       ! write BT to  logical unit LU (for checkpointing)
! bt_read        ! read BT from logical unit LU

! Private subroutines

! nd_create      ! allocate and initialize ND to null settings
! nd_destroy     ! deallocate ND
! nd_write       ! write ND to logical unit LU
! nd_read        ! create and read ND from logical unit LU
! bt_nd_pt_double! double pointer length of array

! bt_initialize( sell, bt, n_pt, check, idlist )
! bt_destroy( bt )
! bt_add( bt, id )
! bt_remove( bt, id )
! bt_purge( bt )
! bt_rebuild( bt, idlist, mode )
! bt_traverse( bt, x, id )
! bt_ball_pts( bt, c, r, npts_max, npts, ids, cp_tests, in_tests )
! bt_param_set( bt, check )
! bt_param_get( bt, check )
! bt_prop( bt, prop )
! bt_write( bt, lu, n_leaves )
! bt_read( sell, bt, idlist, lu, n_leaves )
! nd_create( bt, node_leaf, nd )
! nd_destroy( nd )
! nd_write( nd, lu )
! nd_read( bt, nd, lu )
! bt_nd_pt_double( bt )

!================================================

use sell_m
implicit none

   private nd_create, nd_destroy, nd_write, nd_read, bt_nd_pt_double

   type :: bt_type  !  Binary Tree  (BT) !--------
      type (sell_type), pointer    :: sell     ! SELL to which BT belongs
      type(nd_type),    pointer    :: root     ! root ND of BT
	  type(nd_pointer), pointer    :: nd_pt(:) ! pointer to leaves: nd_pt(id)
	  integer                      :: n_pt     ! dimension of nd_pt

      integer                      :: check    ! level of checking [0,1,2]
      type (id_list_type), pointer :: idlist   ! data structure used in id_list to generate IDs

      real(k_d)                    :: prop(10) ! properties of BT - see bt_prop
	  integer                      :: depth    ! depth of last traverse
   end type bt_type

   type :: nd_type   !  Bounding Ellipsoid (ND) ----------------
      type (bt_type),    pointer :: bt        ! BT to which ND belongs
      type (nd_type),    pointer :: parent    ! nd%parent  is the parent of ND
      type (nd_pointer), pointer :: child(:)  ! nd%child(k)%nd  is the k-th child of ND
	  integer                    :: id        ! nd%id is the ID of ND's leaf (if ND is a leaf of the bt)
	                                          ! nd%id < 0 indicates that ND is a node
	  integer                    :: k         ! ND is parent's k-th child (k = 1 or 2)

	  real(k_d), pointer         :: p(:)      ! point on cutting hyperplane
	  real(k_d), pointer         :: v(:)      ! cutting plane vector:  s(x) = v^T * (x-p)
	  integer                    :: k_first   ! first child to be processed
	  logical                    :: done(2)   ! flag for children having been treated
   end type nd_type  !------------------------------------------

   type :: nd_pointer
      type (nd_type), pointer :: nd
   end type nd_pointer

contains  !===================================================================

subroutine bt_initialize( sell, bt, n_pt, check, idlist )

!  Allocate and initialize (empty) BT.  Optionally set check.

   type (sell_type), pointer              :: sell    ! (in) SELL on which bt is to be formed
   type (bt_type),  pointer               :: bt     ! (out) data structure for BT

!  Optional
   integer, intent(in), optional          :: n_pt    ! initial dimension of pointer array
   type (id_list_type), pointer, optional :: idlist  ! (in) data structure used in id_list to generate IDs

   integer, intent(in), optional          :: check   ! level of checking  =0 minimal; =1 some; =2 maximal

   integer :: i

   allocate( bt )
   bt%sell => sell
   nullify(bt%root)

   bt%n_pt = n_pt_0
   if( present(n_pt) ) bt%n_pt = n_pt
   allocate( bt%nd_pt( bt%n_pt ) )

   do i = 1, bt%n_pt
      nullify( bt%nd_pt(i)%nd )
   end do

   if( present(idlist) ) then
      bt%idlist => idlist
   else
      nullify( bt%idlist)
   endif

   if( present(check) ) then
      call bt_param_set( bt, check )
   else
      bt%check = 0
   endif

   bt%prop  = 0.d0
   bt%depth = 0

   return
end subroutine bt_initialize  !------------------------------------------------------------

subroutine bt_destroy( bt )

!  Destroy and deallocate BT bt.
!  All NDs should first be destroyed by a call to bt_purge

   type (bt_type), pointer :: bt    ! (in) data structure for BT

   if( .not.associated(bt) ) then
      write(0,*) 'bt_destroy: bt not associated; stopping'
	  stop
   endif

   deallocate( bt%nd_pt )
   deallocate( bt )

   return
end subroutine bt_destroy  !------------------------------------------------------------

subroutine bt_add( bt, id )  !  add ELL (ID=id) to BT

! BT is traversed to identify ND_SIBLING.
! A new ND node is added (in place of ND_SIBLING) whose 1st and 2nd
! children are ND and ND_SIBLING, respectively.

   implicit none

   type (bt_type), pointer  :: bt         ! BT
   integer, intent(in)      :: id         ! ID of ELL

   type (nd_type), pointer  :: nd, nd_parent, nd_sibling
   integer                  :: id_status, ids_assigned, id_max, nx, ngeom, &
                               id_sibling, ksp, depth

   if( .not.associated(bt) ) then  !  check input
      write(0,*) 'bt_add: bt not associated; stopping'
	  stop
   endif

   if( bt%check > 1  .and.  associated( bt%idlist ) ) then
	  call id_query( bt%idlist, id, id_status, ids_assigned, id_max )
	  if( id_status /= 1 ) then
		 write(0,*) 'bt_add: id has not been assigned, id = ', id
	     stop
	  endif
   endif

   do while ( id > bt%n_pt )
      call bt_nd_pt_double( bt )
   end do

   nx    = bt%sell%nx
   ngeom = bt%sell%ngeom

!  initialize ND for new leaf
   call nd_create( bt, 2, nd )
   bt%prop(1) = bt%prop(1) + 1.d0

   nd%id  =  id
   bt%nd_pt(id)%nd => nd

! traverse bt to leaf = nd_sibling
   call bt_traverse( bt, bt%sell%ell_pt(id)%ell%geom(1:nx), id_sibling, depth )

! special treatment for first ND in BT
   if( id_sibling == 0 ) then
      bt%root => nd   !  ND is root of bt
	  return
   endif

   nd_sibling => bt%nd_pt(id_sibling)%nd

   bt%prop(2) = max( bt%prop(2), bt%depth+1.d0 )  !  update maximum depth

!  set  nd  to 1st child; nd_sibling to second child of nd_parent

!  allocate and initialize ND = nd_parent for new node
   call nd_create( bt, 1, nd_parent )

   if( associated( nd_sibling%parent ) ) then  !  nd_sibling was not root
      nd_parent%parent => nd_sibling%parent
	  ksp              =  nd_sibling%k
      nd_parent%k      =  ksp
	  nd_sibling%parent%child(ksp)%nd => nd_parent
   else                          !  nd_sibling was root
	  bt%root => nd_parent       !  nd_parent is now root
   endif

   nd_parent%child(1)%nd => nd
   nd_parent%child(2)%nd => nd_sibling

!  re-set nd_sibling

   nd_sibling%parent => nd_parent
   nd_sibling%k      =  2

!  set  nd

   nd%parent => nd_parent
   nd%k      =  1

!  set nd_parent's cutting plane (based on nd_sibling's ELL and nd's center)

   call ell_pt_hyper( nx, bt%sell%ell_pt(id_sibling)%ell%geom(1:nx), &
                          bt%sell%ell_pt(id_sibling)%ell%geom(nx+1:ngeom-2), &
                          bt%sell%ell_pt(id)%ell%geom(1:nx), nd_parent%p, nd_parent%v )

   nd_parent%v = - nd_parent%v ! s < 0 indicates parent%child(1) => nd

   return
end subroutine bt_add  !----------------------------------------------------------

subroutine bt_remove( bt, id )  !  remove leaf with ID=id from BT
   implicit none

   type (bt_type), pointer :: bt ! given BT
   integer, intent(in)     :: id ! ID of leaf to be removed

   integer                   :: k, ks, id_status, ids_assigned, id_max
   type (nd_type), pointer   :: nd, nd_sibling, nd_parent, nd_grandparent

   if( .not.associated(bt) ) then  !  check input
      write(0,*) 'bt_remove: bt not associated; stopping'
	  stop
   endif

   if( id < 1 ) then  !  check that id corresponds to an extant leaf
      write(0,*)'bt_remove: id is non-positive'
	  stop
   elseif( bt%check > 1 ) then
      call id_query( bt%idlist, id, id_status, ids_assigned, id_max )
	  if( id_status /= 1 ) then
	     write(0,*)'bt_remove: id is not assigned ', id, id_status, ids_assigned, id_max
	     stop
	  endif
   endif

   nd         => bt%nd_pt(id)%nd   !  leaf to be removed
   bt%prop(1) = bt%prop(1) - 1.d0

   if( .not.associated(nd%parent) ) then  !  ND is root of bt
      nullify( bt%root )                  !  BT is empty
	  call nd_destroy( nd )
	  nullify( bt%nd_pt( id )%nd )
	  return
   endif

   nd_parent  => nd%parent     ! identify ND's parent and sibling
   k          =  nd%k
   ks         =  3 - k
   nd_sibling => nd_parent%child(ks)%nd

   if( .not.associated(nd_parent%parent) ) then  !  nd_parent is root of BT
      nd%bt%root => nd_sibling                   !  nd's sibling is sole surviving ND
	  nd_sibling%k = 0
	  nullify( nd_sibling%parent )

   else  !  grandparent exists

      nd_grandparent    => nd_parent%parent   !  new parent of sibling
      nd_sibling%parent => nd_grandparent
	  k                 = nd_parent%k         !  index of parent's child (old and new)
      nd_sibling%k      = k
      nd_grandparent%child(k)%nd => nd_sibling
   endif

   call nd_destroy( nd )
   call nd_destroy( nd_parent )
   nullify( bt%nd_pt( id )%nd )

   return
end subroutine bt_remove   !-------------------------------------------------------

subroutine bt_purge( bt )   !  remove all NDs from BT

   implicit none
   type (bt_type), pointer  :: bt   ! given BT

   type (nd_type), pointer  :: nd, parent
   integer                  :: i

   if( .not.associated(bt) ) then  !  check input
      write(0,*) 'bt_purge: bt not associated; stopping'
	  stop
   endif

!  start at root (if it exists)

   if( .not.(associated(bt%root) ) ) return  !  already purged

   nd => bt%root  !  root
   nullify( parent )            !  ...has no parent

   do   !  on each iteration of loop,  ND  and  parent  must be defined (parent may be null)

      if( nd%id < 1  ) then  !  ND is a node

         if( associated(nd%child(1)%nd) ) then      !  go to child(1) if it exists
	        parent => nd
		    nd     => nd%child(1)%nd
		    nullify( parent%child(1)%nd )
		    cycle

         elseif( associated(nd%child(2)%nd) ) then  !  go to child(2) if it exists
	        parent => nd
		    nd     => nd%child(2)%nd
		    nullify( parent%child(2)%nd )
		    cycle
		 endif

	  endif

!  ND  is childless: destroy and go to parent (if it exists)

	  call nd_destroy( nd )                   !  destroy  nd
	  if( .not.associated( parent ) ) exit    !  all done
	  nd => parent                            !  go to parent

	  if( associated(nd%parent) ) then        !  identify parent (if any)
		 parent => nd%parent
	  else
		 nullify( parent )
	  endif

   end do

   nullify( bt%root )

   do i = 1, bt%n_pt
      nullify( bt%nd_pt(i)%nd )
   end do

   bt%prop = 0.d0  ! re-initialize bt%prop

   return

end subroutine bt_purge  !------------------------------------------------------------------

subroutine bt_rebuild( bt, idlist, mode )   !  rebuild BT

   implicit none
   type (bt_type), pointer      :: bt     !  BT to be rebuilt
   type (id_list_type), pointer :: idlist !  IDLIST
   integer, intent(in)          :: mode   !  method of rebuilding

! mode = 1 - purge BT, and then re-add ELLs in a random order
! (only mode currently supported)

   integer              :: i, ids_assigned, id_max
   integer, allocatable :: id_array(:)

   if( .not.associated(bt) ) then  !  check input
      write(0,*) 'bt_rebuild: bt not associated; stopping'
	  stop
   endif

   call id_stats( idlist, ids_assigned, id_max )
   if( ids_assigned < 3 ) return  !  unique configuration

   if( mode == 1 ) then  !  purge and re-add ELLs in random order

	  allocate( id_array(ids_assigned) )
	  call id_rand_array( idlist, id_array )
      call bt_purge( bt )

	  do i = 1, ids_assigned
	     call bt_add( bt, id_array(i) )
	  end do

	  deallocate( id_array )

   else
      write(0,*)'bt_rebuild: invalid mode = ', mode
	  stop
   endif

   return
end subroutine bt_rebuild  !------------------------------------------------------------------

subroutine bt_traverse( bt, x, id, depth )  !  traverse BT based on x

!  Traverse BT based on x to indentify leaf with ID=id.
!  If the root of the BT does not exist, id=0 is returned.
!  bt%depth is set to the depth of the traverse.

   implicit none

   type (bt_type), pointer   :: bt            ! given BT
   real(k_d), intent(in)     :: x(bt%sell%nx) ! query point
   integer, intent(out)      :: id            ! ID of leaf
   integer, intent(out)      :: depth         ! depth of traverse path

   type (nd_type), pointer   :: nd
   real(k_d)                 :: s

   if( .not.associated( bt ) ) then  ! check that bt exists
      write(0,*)'bt_traverse; bt not associated'
	  stop
   endif

   if( .not.associated( bt%root) ) then
      id       = 0 !  indicate that root does not exist
	  depth    = 0
	  bt%depth = 0
	  return
   endif

   nd         => bt%root            !  start traverse at root
   bt%prop(3) =  bt%prop(3) + 1.d0  !  increment number of traverses
   depth      =  1

   do  !  loop over NDs  ------------------------------------------

      if( nd%id > 0 ) then  !  leaf has been reached - all done
	     id         = nd%id
	     bt%prop(4) = bt%prop(4) + depth  !  increment sum of depths
		 bt%depth   = depth
	     return
      endif

	  s = sum( nd%v * (x-nd%p) )  !  compute s
      if( s < 0.d0 ) then   !  go to child 1 or 2 depending on s
	     nd => nd%child(1)%nd
	  else
	     nd => nd%child(2)%nd
	  endif

	  depth = depth + 1

   end do

end subroutine bt_traverse  !---------------------------------------------------------------

subroutine bt_ball_pts( bt, c, r, npts_max, npts, ids, dsq, cp_tests, in_tests )

!  Given a ball, centered at  c  of radius r, traverse the binary tree BT
!  to identify up to npts_max points covered by the ball.

   implicit none

   type (bt_type), pointer   :: bt            ! given BT
   real(k_d), intent(in)     :: c(bt%sell%nx) ! center of ball
   real(k_d), intent(in)     :: r             ! radius of ball
   integer, intent(in)       :: npts_max      ! max. no. of points to be returned
   integer, intent(out)      :: npts          ! no. of points found
   integer, intent(out)      :: ids(npts_max) ! IDs of points
   real(k_d), intent(out)    :: dsq(npts_max) ! square of distance of point from center
   integer, intent(out)      :: cp_tests      ! number of cutting plane tests
   integer, intent(out)      :: in_tests      ! number of in-ball tests

   type (nd_type), pointer   :: nd
   real(k_d)                 :: s, x(bt%sell%nx), rsq, xsq
   logical                   :: from_above
   integer                   :: nx, id

   if( .not.associated( bt ) ) then  ! check that bt exists
      write(0,*)'bt_traverse; bt not associated'
	  stop
   endif

   npts     = 0  !  null initialization
   cp_tests = 0
   in_tests = 0

   if( npts_max < 1 ) return

   if( .not.associated( bt%root) ) return

   rsq        = r*r
   nx         = bt%sell%nx
   nd         => bt%root     !  start traverse at root
   from_above = .true.

   do  !  loop over NDs  ------------------------------------------
      if( from_above ) then

	     if( nd%id > 0 ) then  !  ND is a leaf - test for in ball
	        id  = nd%id
		    x   = bt%sell%ell_pt(id)%ell%geom(1:nx) - c
		    xsq = sum( x*x )
            in_tests = in_tests + 1

            if( xsq <= rsq ) then
               npts = npts + 1  !  in ball
			   ids(npts) = id
			   dsq(npts) = xsq
			   if( npts == npts_max ) return
		    endif

            if( .not.associated( nd%parent ) ) return  ! BT consists of single leaf

			nd => nd%parent  !  go to parent
			from_above = .false.
			cycle

		 else  ! node - for the first time: determine which children need to be examined
		    nd%done = .false.

	        s = sum( nd%v * (c-nd%p) ) !  compute signed distance to center of ball
			cp_tests = cp_tests + 1

			if( s > r ) then
			   nd%done(1) = .true.  ! no need to examine child 1

			elseif( s < -r ) then
			   nd%done(2) = .true.  ! no need to examine child 2
			endif
		 endif

	  endif

! at node; node%done has been set: go to child or parent

      if( .not.nd%done(1) ) then
	     nd%done(1) = .true.
	     nd         => nd%child(1)%nd
		 from_above = .true.

	  elseif( .not.nd%done(2) ) then
	     nd%done(2) = .true.
	     nd         => nd%child(2)%nd
		 from_above = .true.

	  else
	     if( .not.associated( nd%parent ) ) return  ! at root: all done
	     nd => nd%parent  !  go to parent
		 from_above = .false.
	  endif

   end do

end subroutine bt_ball_pts  !---------------------------------------------------------------

subroutine bt_param_set( bt, check )

! set parameters

   type (bt_type), pointer         :: bt    ! data structure for BT.
   integer, intent(in), optional   :: check  ! level of checking [0,1,2]

   if( .not.associated(bt) ) then
      write(0,*) 'bt_param_set: bt not associated; stopping'
	  stop
   endif

   if( present(check) ) then
      if( check < 0  .or.  check >2 ) then
	     write(0,*) 'bt_param_set: invalid  check = ', check
	     stop
	  endif

      bt%check = check
   endif

   return
end subroutine bt_param_set  !--------------------------------------------------------------

subroutine bt_param_get( bt, check )

! get parameters

   type (bt_type), pointer         :: bt   ! data structure for bt.
   integer, intent(out), optional   :: check  ! level of checking

   if( .not.associated(bt) ) then
      write(0,*) 'bt_param_get: bt not associated; stopping'
	  stop
   endif

   if( present(check) ) check = bt%check

   return
end subroutine bt_param_get  !--------------------------------------------------------------

subroutine bt_prop( bt, prop )

! get properties of BT

   type (bt_type), pointer :: bt       ! data structure for BT.
   real(k_d), intent(out)  :: prop(10) ! properties

! 1 - number of leaves in tree
! 2 - maximum depth of tree
! 3 - number of traverses performed
! 4 - sum of depths of traverses
! 5 - average depth of traverses

   prop = 0.d0

   if( .not.associated(bt) ) then
      write(0,*) 'bt_prop: bt not associated; stopping'
	  stop
   endif

   prop(1:4) = bt%prop(1:4)
   prop(5)   = prop(4) / max( prop(3), 1.d0 )

   return
end subroutine bt_prop  !--------------------------------------------------------------

subroutine bt_write( bt, lu, n_leaves )  !  write bt to  logical unit LU

   type (bt_type), pointer   :: bt        !  BT to be written
   integer, intent(in)       :: lu        !  LU: must be open for writing
   integer, intent(out)      :: n_leaves  ! number of leaves written

   type (nd_type), pointer   :: nd
   logical                   :: opened, from_above
   integer                   :: i, id, id_parent

   if( .not.associated(bt) ) then
      write(0,*)'bt_write: bt not associated'
	  stop
   endif

   inquire( lu, opened=opened )
   if( .not.opened ) then
      write(0,*)'bt_write: lu not opened ', lu
	  stop
   endif

   write(lu, err=100 ) bt%n_pt, bt%check, bt%prop

   n_leaves = 0
   id       = 0
   do i = 1, bt%n_pt  !  write leaf NDs
      if( associated(bt%nd_pt(i)%nd) ) then
	     n_leaves = n_leaves + 1
		 id       = bt%nd_pt(i)%nd%id
         write(lu, err=110 ) id
		 call nd_write( bt%nd_pt(i)%nd, lu )
	  else
	     write(lu) 0
	  endif
   end do

   if( n_leaves < 2 ) return  !  all done if no nodes

!  first traverse of bt to assign unique negative id's to nodes, and to write node NDs
!  note that the root node has id = -1

   id         = 0
   from_above = .true.
   nd         => bt%root

   do
      if( from_above ) then  ! entering node or leaf for the first time

         if( nd%id < 0 ) then   !  first time at node
	        id = id - 1
	        nd%id = id          ! assign ID and write node ND
	        call nd_write( nd, lu )

			nd%done(2) = .false.
		    nd => nd%child(1)%nd  !  go to child 1, cycle
         else  !  leaf:  do nothing; return to parent
			from_above = .false.
			if( .not.associated(nd%parent) ) exit
			nd => nd%parent  !  go to parent, cycle
		 endif

	  else ! from_below
	     if( .not.nd%done(2) ) then
		    nd%done(2) = .true.
			from_above = .true.
		    nd => nd%child(2)%nd  !  go to child 2, cycle
         else  ! go to parent
			from_above = .false.
			if( .not.associated(nd%parent) ) exit
			nd => nd%parent  !  go to parent, cycle
		 endif
	  endif
   end do

   if( -id /= n_leaves-1 ) then
      write(0,*)'bt_write: node miscount, ', -id, n_leaves-1
	  stop
   endif

! second traverse of bt nodes to write ID of ND, parent and children
   from_above = .true.
   nd         => bt%root

   do
      if( from_above ) then  !first time at node
	     if( associated(nd%parent) ) then
		    id_parent = nd%parent%id
		 else
		    id_parent = 0
		 endif

		 write(lu, err=120 ) nd%id, id_parent, nd%child(1)%nd%id, nd%child(2)%nd%id

		 nd%done(2) = .false.
		 if( nd%child(2)%nd%id > 0 ) nd%done(2) = .true.

		 if( nd%child(1)%nd%id < 0 ) then
		    nd         => nd%child(1)%nd
		 elseif( nd%child(2)%nd%id < 0 ) then
		    nd%done(2) = .true.
		    nd         => nd%child(2)%nd
		 else
		    nd%done = .true.
		    from_above = .false.
			if( .not.associated(nd%parent) ) exit
			nd => nd%parent
		 endif

      else !  from below

         if( .not.nd%done(2) ) then
		    from_above = .true.
			nd%done(2) = .true.
		    nd => nd%child(2)%nd
		 else
			if( .not.associated(nd%parent) ) exit
			nd => nd%parent
		 endif

      endif
   end do

   return

100   write(0,*)'bt_write: error writing n_pt'
      stop

110   write(0,*)'bt_write: error writing bt%nd_pt(i)%nd%id'
      stop

120   write(0,*)'bt_write: error writing nd%id'
      stop

end subroutine bt_write    !-----------------------------------------------

subroutine bt_read( sell, bt, idlist, lu, n_leaves )  !  read BT from logical unit LU

!  BT should not be allocated prior to call

   type (sell_type), pointer    :: sell
   type (bt_type), pointer      :: bt
   type (id_list_type), pointer :: idlist
   integer, intent(in)          :: lu
   integer, intent(out)         :: n_leaves  ! number of leaves read

   type (nd_type), pointer      :: nd
   type (nd_pointer), pointer   :: node_pt(:)
   real(k_d) :: prop(10)
   integer   :: n_pt, check

   integer            :: i, id, n_nodes, id_node, id_parent, id_child1, id_child2
   logical            :: opened, from_above

   inquire( lu, opened=opened )
   if( .not.opened ) then
      write(0,*)'bt_write: lu not opened ', lu
	  stop
   endif

   read(lu, err=100) n_pt, check, prop

   call bt_initialize( sell, bt, n_pt=n_pt, check=check, idlist=idlist )
   bt%prop = prop

   n_leaves = 0
   do i = 1, bt%n_pt  !  read leaf NDs
      read(lu, err=110 ) id
      if( id > 0 ) then
	     n_leaves = n_leaves + 1
		 call nd_read( bt, nd, lu )
		 if( nd%id /= id ) then
		    write(0,*)'bt_read: id mismatch ', nd%id, id
			stop
		 endif
		 bt%nd_pt(i)%nd => nd
	  else
	     nullify( bt%nd_pt(i)%nd )
	  endif
   end do

   if( n_leaves < 2 ) then
      if( n_leaves == 1 ) bt%root => nd
	  return  !  all done if no nodes
   endif

   n_nodes = n_leaves - 1
   allocate( node_pt(n_nodes) )

   do id = 1, n_nodes  !  read node NDs
      call nd_read( bt, nd, lu )
	  if( nd%id /= -id ) then
	     write(0,*)'bt_read: node ID mismatch ', nd%id, -id
		 stop
	  endif
	  node_pt(id)%nd => nd
   end do

   bt%root => node_pt(1)%nd

!  traverse bt to assign child(1:2) and parent

   from_above = .true.
   nd         => bt%root

   do
      if( from_above ) then  ! entering node or leaf for the first time
	     read(lu, err=120 ) id_node, id_parent, id_child1, id_child2

		 nd => node_pt(-id_node)%nd

		 if( id_parent < 0 ) then  !  set parent
		    nd%parent => node_pt(-id_parent)%nd
		 else
		    nullify( nd%parent )
		 endif

         nd%done = .false.
		 if( id_child1 < 0 ) then  !  set child 1
		    nd%child(1)%nd => node_pt(-id_child1)%nd
		 else
            nd%child(1)%nd => bt%nd_pt(id_child1)%nd
			bt%nd_pt(id_child1)%nd%parent => nd
			nd%done(1)     =  .true.
		 endif

		 if( id_child2 < 0 ) then  !  set child 2
		    nd%child(2)%nd => node_pt(-id_child2)%nd
		 else
            nd%child(2)%nd => bt%nd_pt(id_child2)%nd
			bt%nd_pt(id_child2)%nd%parent => nd
			nd%done(2)     =  .true.
		 endif

		 if( .not.nd%done(1) ) then
		    nd%done(1) = .true.
		    nd         => nd%child(1)%nd
		 elseif( .not.nd%done(2) ) then
		    nd%done(2) = .true.
		    nd         => nd%child(2)%nd
		 else
		    from_above = .false.
			if( .not.associated(nd%parent) ) exit
			nd => nd%parent
		 endif

      else !  from below

         if( .not.nd%done(2) ) then
		 	nd%done(2) = .true.
		    from_above = .true.
		    nd => nd%child(2)%nd
		 else
			if( .not.associated(nd%parent) ) exit
			nd => nd%parent
		 endif

      endif
   end do

   deallocate( node_pt )

   return

100   write(0,*)'bt_read:  error reading n_pt'
      stop

110   write(0,*)'bt_read:  error reading id'
      stop

120   write(0,*)'bt_read:  error reading id_node'
      stop

end subroutine bt_read    !-----------------------------------------------

subroutine nd_create( bt, node_leaf, nd )  !  allocate and initialize ND to null settings

   type (bt_type), pointer  :: bt        !  (in)  BT to which ND belongs
   integer, intent(in)      :: node_leaf ! =1 for node; =2 for leaf
   type (nd_type),  pointer :: nd        !  (out) ND

   if( .not.associated( bt ) ) then  ! check that BT exists
      write(0,*)'nd_create; BT not associated'
	  stop
   endif

   if( node_leaf < 1  .or.  node_leaf > 2 ) then
      write(0,*)'nd_create; bad value of node_leaf = ', node_leaf
	  stop
   endif

   allocate( nd )
   nd%bt => bt
   nullify(  nd%parent )

   nd%id      = -1
   nd%k       = -1
   nd%k_first = -1
   nd%done    = .false.

   if( node_leaf == 1 ) then  ! node
      allocate( nd%child(2) )
      nullify(  nd%child(1)%nd )
      nullify(  nd%child(2)%nd )
      allocate( nd%p(bt%sell%nx) )
      allocate( nd%v(bt%sell%nx) )
	  nd%p       = -1.d0
      nd%v       = -1.d0
   else                        ! leaf
      nullify( nd%child )
      nullify( nd%p )
      nullify( nd%v )
   endif

   return
end subroutine nd_create  !------------------------------------------------

subroutine nd_destroy( nd )  !   deallocate and nullify ND

   type (nd_type), pointer :: nd

   if( .not.associated( nd ) ) then  ! check that ND exists
      write(0,*)'nd_destroy; nd not associated'
	  stop
   endif

   if( nd%id < 1 ) then  !  node
      deallocate( nd%child )
      deallocate( nd%p )
      deallocate( nd%v )
   endif

   deallocate( nd )

   return
end subroutine nd_destroy  !------------------------------------------------

subroutine nd_write( nd, lu )  !  write ND to logical unit LU

   type (nd_type),  pointer :: nd
   integer, intent(in)      :: lu

   integer :: nx

   nx    = nd%bt%sell%nx

! Note: nd%k_first, nd%done(1:2) are not written (a) because they are temporary, and not needed
!       for a perfect restart (b) to ensure repeatability

   write(lu, err=100 ) nd%id
   write(lu, err=110 ) nd%k

   if( nd%id < 0 ) write(lu, err=120) nd%p(1:nx), nd%v(1:nx)

   return

100  write(0,*)'nd_write: error writing id'
     stop

110  write(0,*)'nd_write: error writing k'
     stop

120  write(0,*)'nd_write: error writing p, v'
     stop

end subroutine nd_write  !------------------------------------------------

subroutine nd_read( bt, nd, lu )  !  create and read ND from logical unit LU

   type (bt_type), pointer  :: bt        !  (in)  BT to which ND belongs
   type (nd_type),  pointer :: nd        !  (out) ND
   integer, intent(in)      :: lu

   integer :: id, node_leaf, nx

   nx    = bt%sell%nx

   read(lu, err=100 ) id

   node_leaf = 1
   if( id > 0 ) node_leaf = 2

   call nd_create( bt, node_leaf, nd )

   nd%id = id
   read(lu, err = 110 ) nd%k

   if( nd%id < 0 ) read(lu, err = 120 ) nd%p(1:nx), nd%v(1:nx)

   return

100  write(0,*)'nd_read: error reading id'
     stop

110  write(0,*)'nd_read: error reading k'
     stop

120  write(0,*)'nd_read: error reading p, v'
     stop

end subroutine nd_read  !------------------------------------------------

subroutine bt_nd_pt_double( bt )

!  double size of array bt%nd_pt

      type (bt_type),   pointer :: bt

      type (nd_pointer), pointer :: pt_temp(:)
	  integer                    :: id

!  allocate temporary array pt_temp of twice size
   	  allocate( pt_temp(2*bt%n_pt) )

	  do id = 1, 2*bt%n_pt
	     nullify( pt_temp(id)%nd )
      end do

!  copy ell%pt to pt_temp
	  do id = 1, bt%n_pt
	     if( associated( bt%nd_pt(id)%nd ) ) pt_temp(id)%nd => bt%nd_pt(id)%nd
		 nullify( bt%nd_pt(id)%nd )
	  end do

!  double size of bt%nd_pt
	  deallocate( bt%nd_pt )

      bt%n_pt = 2 * bt%n_pt
	  allocate( bt%nd_pt(bt%n_pt) )

	  do id = 1, bt%n_pt

	     if( associated( pt_temp(id)%nd ) ) then
		    bt%nd_pt(id)%nd => pt_temp(id)%nd
	     else
		    nullify( bt%nd_pt(id)%nd )
		 endif

		 nullify(pt_temp(id)%nd)
	  end do

	  deallocate( pt_temp )

	  return

end subroutine bt_nd_pt_double   !-----------------------------------------------

end module sell_bt


text:
!	Get integer weights for weighted sum.
!	skillname is the array of factor names.
!	numcatobs is the array of numbers of categories in items.
!	slopedim is the array of dimensions for each factor and item.
!	weightnames are names of weights and weights are weights.
subroutine getweights(skillname,numcatobs,slopedim,weightnames,weights)
implicit none
interface
    subroutine readweights(weightname,weight)
        implicit none
        character(len=32),intent(inout)::weightname
        integer,intent(inout)::weight(:)
    end subroutine readweights
end interface
character(len=32),intent(in)::skillname(:)
integer,intent(in)::numcatobs(:),slopedim(:,:)
character(len=32),intent(out)::weightnames(:)
integer,intent(out)::weights(:,:)
!	buff is a buffer.
!	weightname is used for a nameless group.
character(len=4)::buff
character(len=32)::weightname
!	cat counts categories.
!	count weights.
!	counter finds position in weight.
!	io is error flag.
!	item counts items.
!
integer::cat,count,counter,io,item
!	weight to be read.
integer,allocatable::weight(:)

allocate(weight(size(weights,1)),stat=io)
if(io/=0)stop 'Allocation for category weights failed.'

weights=0
if(size(skillname)>=size(weightnames))then
	weightnames=skillname(1:size(weightnames))
	if(size(skillname)>size(weightnames))then
		if(size(weightnames)>1)then
			weightnames(size(weightnames))='Remainder'
		else
			weightnames(1)='Total'
		end if
	end if
	counter=0
	do item=1,size(numcatobs)
		do count=1,size(skillname)
			do cat=1,numcatobs(item)
				if(slopedim(count,item)>0)weights(counter+cat,min(count,size(weightnames)))=cat-1
			end do

		end do
		counter=counter+numcatobs(item)
	end do
else
	weightnames(1:size(skillname))=skillname
	counter=0
	do item=1,size(numcatobs)
		do count=1,size(skillname)
			do cat=1,numcatobs(item)
				if(slopedim(count,item)>0)weights(counter+cat,count)=cat-1
			end do

		end do
		counter=counter+numcatobs(item)
	end do
	do count=size(skillname)+1,size(weightnames)
		if(count==size(skillname)+1.and.count>2)then
			weightnames(count)='Total'
			counter=0
			do item=1,size(numcatobs)
				do cat=1,numcatobs(item)
					weights(counter+cat,count)=cat-1
				end do
				counter=counter+numcatobs(item)
			end do
		else
			write(buff,'(i4)') count
			weightnames(count)='Sum_'//trim(adjustl(buff))
		end if
	end do
end if

do count=1,size(weightnames)
	weightname=weightnames(count)
	weight=weights(:,count)
	call readweights(weightname,weight)
	weightnames(count)=weightname
	weights(:,count)=weight
end do

return
end subroutine getweights
subroutine readweights(weightname,weight)
    implicit none
    character(len=32),intent(inout)::weightname
    integer,intent(inout)::weight(:)
    integer::io
    namelist/readweight/weightname,weight
    read(*,nml=readweight,iostat=io)
    if(io/=0)stop "Weights not read successfully."

    return
end subroutine readweights

max_stars_repo_name: EducationalTestingService/MIRT

max_issues_repo_name: EducationalTestingService/MIRT

max_forks_repo_name: EducationalTestingService/MIRT


text:
	PROGRAM XdmfCreate
        !--------------------------------------------------------------------------------
        ! .. Purpose ..
        ! XdmfCreate is a postprocessing program which creates header files for
        ! Paraview and Visit
        ! It uses the INPUTFILE and assumes that the filenames in the program
        ! are
        ! consistent with those in the current file.
        !
        ! .. PARAMETERS .. INITIALIZED IN INPUTFILE
        ! time = start time of the simulation
        ! Nx = power of two, number of modes in the x direction
        ! Ny = power of two, number of modes in the y direction
        ! Nz = power of two, number of modes in the z direction
        ! Nt = the number of timesteps
        ! plotgap = the number of timesteps to take before plotting
        ! Lx = definition of the periodic domain of computation in x direction
        ! Ly = definition of the periodic domain of computation in y direction
        ! Lz = definition of the periodic domain of computation in z direction
        ! Dt = the time step
        ! Es = focusing or defocusing
        !
        ! REFERENCES
        !
        ! www.xdmf.org
        ! Paraview users mailing list
        ! VisIt users mailing list
        !
        ! ACCURACY
        !
        ! ERROR INDICATORS AND WARNINGS
        !
        ! FURTHER COMMENTS
        !------------------------------------------------------------------------------------
        ! EXTERNAL ROUTINES REQUIRED
        IMPLICIT NONE
        ! .. Scalar Arguments ..
        INTEGER(KIND=4) :: Nx, Ny, Nz, Nt, plotgap
        REAL(KIND=8)    :: Lx, Ly, Lz, DT, Es, time
        ! .. Local scalars ..
        INTEGER(KIND=4) :: stat,plotnum,ind,n,numplots
        REAL(KIND=8), PARAMETER :: pi=3.14159265358979323846264338327950288419716939937510d0
        ! .. Local Arrays ..
        CHARACTER*50    :: InputFileName, OutputFileName, OutputFileName2
        CHARACTER*10    :: number_file
        InputFileName='inputfile'
        OPEN(unit=11,FILE=trim(InputFileName),status="OLD")
        REWIND(11)
        READ(11,*) Nx, Ny, Nz, Nt, plotgap, Lx, Ly, Lz, Es, DT, time
        CLOSE(11)
        time=0.0
        plotnum=1
        numplots=1+Nt/plotgap
        OutputFileName='udata.xmf'
        OPEN(unit=11,FILE=trim(OutputFileName),status="UNKNOWN")
        REWIND(11)
        WRITE(11,'(a)') "<?xml version=""1.0"" ?>"
        WRITE(11,'(a)') "<!DOCTYPE Xdmf SYSTEM ""Xdmf.dtd"" []>"
        WRITE(11,'(a)') "<Xdmf xmlns:xi=""http://www.w3.org/2001/XInclude"" Version=""2.0"">"
        WRITE(11,'(a)') "<Domain>"
        WRITE(11,'(a)') " <Topology name=""topo"" TopologyType=""3DCoRectMesh"""
        WRITE(11,*) "  Dimensions=""",Nx,Ny,Nz,""">"
        WRITE(11,'(a)') " </Topology>"
        WRITE(11,'(a)') " <Geometry name=""geo"" Type=""ORIGIN_DXDYDZ"">"
        WRITE(11,'(a)') "  <!-- Origin -->"
        WRITE(11,'(a)') "  <DataItem Format=""XML"" Dimensions=""3"">"
        WRITE(11,*) -pi*Lx, -pi*Ly, -pi*Lz
        WRITE(11,'(a)') "  </DataItem>"
        WRITE(11,'(a)') "  <!-- DxDyDz -->"
        WRITE(11,'(a)') "  <DataItem Format=""XML"" Dimensions=""3"">"
        WRITE(11,*) REAL(2.0*pi*Lx/Nx,kind(0d0)), REAL(2.0*pi*Ly/Ny,kind(0d0)), REAL(2.0*pi*Lz/Nz,kind(0d0))
        WRITE(11,'(a)') "   </DataItem>"
        WRITE(11,'(a)') "  </Geometry>"
        WRITE(11,'(a)')
        WRITE(11,'(a)') "       <Grid Name=""TimeSeries"" GridType=""Collection"" CollectionType=""Temporal"">"
        WRITE(11,'(a)') "               <Time TimeType=""HyperSlab"">"
        WRITE(11,'(a)') "                       <DataItem Format=""XML"" NumberType=""Float"" Dimensions=""3"">"
        WRITE(11,*) "                   0.0", dt," 2"
        WRITE(11,'(a)') "                       </DataItem>"
        WRITE(11,'(a)') "               </Time>"

        DO n=1,numplots

                OutputFileName = 'data/u'
                ind = index(OutputFileName,' ') - 1
                WRITE(number_file,'(i0)') 10000000+plotnum
                OutputFileName = OutputFileName(1:ind)//number_file
                ind = index(OutputFileName,' ') - 1
                OutputFileName = OutputFileName(1:ind)//'.datbin'

                OutputFileName2 = "<Grid Name=""T"
                WRITE(number_file,'(i0)') n
                OutputFileName2 = OutputFileName2(1:13)//number_file
                ind =  index(number_file,' ') - 1
                OutputFileName2 = OutputFileName2(1:13+ind)//'" GridType="Uniform">'
                plotnum=plotnum+1

                WRITE(11,'(a)')
                WRITE(11,'(3X,a)') "   ",OutputFileName2
                WRITE(11,'(a)') "       <Topology Reference=""/Xdmf/Domain/Topology[1]""/>"
                WRITE(11,'(a)') "       <Geometry Reference=""/Xdmf/Domain/Geometry[1]""/>"
                WRITE(11,'(a)') "        <Attribute Name=""Absolute Value of Field"" Center=""Node"">"
                WRITE(11,'(a)') "               <DataItem Format=""Binary"" "
                WRITE(11,'(a)') "                DataType=""Float"" Precision=""8"" Endian=""Little"" "
                WRITE(11,*) "            Dimensions=""",Nx, Ny, Nz,""">"
                WRITE(11,'(16X,a)') "           ",OutputFileName
                WRITE(11,'(a)') "               </DataItem>"
                WRITE(11,'(a)') "       </Attribute>"
                WRITE(11,'(3X,a)') "</Grid>"
        END DO
        WRITE(11,'(a)')"        </Grid>"
        WRITE(11,'(a)') "</Domain>"
        WRITE(11,'(a)') "</Xdmf>"
        CLOSE(11)
        END PROGRAM XdmfCreate

max_stars_repo_name: bcloutier/PSNM

max_issues_repo_name: bcloutier/PSNM

max_forks_repo_name: bcloutier/PSNM


text:
*----------------------------------------------------------------------*
      subroutine import_hamint_gamess(hlist,str_info,orb_info)
*----------------------------------------------------------------------*
*     import one- and two-electron integrals from GAMESS
*     environment
*     we need:
*      MOINTS for fock matrix
*      MOINTS for two-electron integrals
*----------------------------------------------------------------------*
      implicit none

      include 'stdunit.h'
      include 'def_graph.h'
      include 'def_strinf.h'
      include 'def_orbinf.h'
      include 'def_filinf.h'
      include 'def_operator.h'
      include 'def_me_list.h'
      include 'par_gamess.h'

      type(me_list), intent(inout) ::
     &     hlist
      type(strinf), intent(in) ::
     &     str_info
      type(orbinf), intent(in) ::
     &     orb_info

      ! read reference energy and fock matrix from MOINTS
      ! and sort fock matrix into operator file
      call import_fock_dalton(hlist,str_info,orb_info,.true.,moints)

      ! get 2-electron integrals and sort them as well
      call import_h2_gamess(hlist,str_info,orb_info)

      return
      end

ext: f

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


program krackentest
use M_kracken95
implicit none
integer :: ierror

! call with an arbitrary verb name and a string that defines
! allowable switch names and default values

      call kracken('mycommand',&
     &' -i 10 -r 200.3 -par1 "#N#" -par2 -par3 10 20 30 -files ', &
     & ierror)

! that's it. Now you can just retrieve the values as strings
! using names of the form VERB_SWITCHNAME anywhere in your program
! the special name VERB_oo is for the string before a switch

      call samples()
contains

subroutine samples
use M_kracken95
implicit none
character(len=255) :: value
integer :: iilen, ier
real    :: anumber

! Getting everything before a switch
      call retrev('mycommand_oo',value,iilen,ier)
      write(*,*)'before any switch=',value(:iilen)

! Getting an integer from -i value
      call retrev('mycommand_i',value,iilen,ier)
      write(*,*)'value for -i =',value(:iilen)
      call string_to_real(value,anumber,ier)
      write(*,*)'divided by 2 is ',nint(anumber)/2

! Getting a real number from -r value
      call retrev('mycommand_r',value,iilen,ier)
      write(*,*)'value for -r =',value(:iilen)
      call string_to_real(value,anumber,ier)
      write(*,*)'divided by 2 is ',anumber/2

! Getting a string (with a default) from -par1
      call retrev('mycommand_par1',value,iilen,ier)
      write(*,*)'value for -par1 =',value(:iilen)

! Getting a string (with no default) from -par2
      call retrev('mycommand_par2',value,iilen,ier)
      write(*,*)'value for -par2 =',value(:iilen)

! Getting a string with a multi-word default and splitting it from -par3
      call listem('mycommand_par3','-par3',.true.)

! a list of files
      call listem('mycommand_files','-files',.false.)

end subroutine samples

subroutine listem(keyword,label,toreal)
use M_kracken95
implicit none
!     Just getting fancy and showing the use of DELIM()
!     SAMPLE that decomposes list of strings and optionally, numbers
!     and prints it.
!     delimit with space, comma, or colon
logical toreal
character(len=*)    :: keyword
character(len=*)    :: label
character(len=255)  :: value

!     for DELIM call if you want to break down a list
character(len=255)  ::  array(132)
integer ibegin(132)
integer iterm(132)
integer :: i10, igot, iilen, ier
real    :: anumber

!     get the value of the keyword and the length
      call retrev(keyword,value,iilen,ier)
      write(*,*)'value for ',label,'=',value(:iilen)
!     split the list into one word per array
      call delim(value,array,132,igot,ibegin,iterm,iilen,' ,:')
!     convert each word into a numeric value
      do i10=1,igot
         write(*,*)i10,') ',array(i10)(:len_trim(array(i10)))
         if(toreal)then
            call string_to_real(array(i10),anumber,ier)
            if(ier.eq.0)then
               write(*,*)' which is a number'
            else
               write(*,*) ' which is not a number'
            endif
         endif
      enddo
end subroutine listem
end program krackentest