text:
C$Procedure ZZGFOCU ( GF, occultation utilities )

      SUBROUTINE ZZGFOCU ( OCCTYP, FRONT,  FSHAPE, FFRAME,
     .                     BACK,   BSHAPE, BFRAME, OBSRVR,
     .                     ABCORR, TIME,   OCSTAT          )

C$ Abstract
C
C     SPICE Private routine intended solely for the support of SPICE
C     routines.  Users should not call this routine directly due
C     to the volatile nature of this routine.
C
C     This routine contains the entry points that produce the
C     computations needed for solving for occultation states
C     in the geometry finding routines.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     GF
C     NAIF_IDS
C     PCK
C     SPK
C     TIME
C
C$ Keywords
C
C     SEARCH
C     GEOMETRY
C     OCCULTATION
C
C$ Declarations

      IMPLICIT NONE

      INCLUDE 'dsk.inc'
      INCLUDE 'gf.inc'
      INCLUDE 'zzabcorr.inc'
      INCLUDE 'zzdsk.inc'
      INCLUDE 'zzocced.inc'

      CHARACTER*(*)         OCCTYP
      CHARACTER*(*)         FRONT
      CHARACTER*(*)         FSHAPE
      CHARACTER*(*)         FFRAME
      CHARACTER*(*)         BACK
      CHARACTER*(*)         BSHAPE
      CHARACTER*(*)         BFRAME
      CHARACTER*(*)         OBSRVR
      CHARACTER*(*)         ABCORR
      DOUBLE PRECISION      TIME
      LOGICAL               OCSTAT

C$ Brief_I/O
C
C     VARIABLE  I/O  Entry points
C     --------  ---  --------------------------------------------------
C     OCCTYP     I   ZZGFOCIN
C     FRONT      I   ZZGFOCIN
C     FSHAPE     I   ZZGFOCIN
C     FFRAME     I   ZZGFOCIN
C     BACK       I   ZZGFOCIN
C     BSHAPE     I   ZZGFOCIN
C     BFRAME     I   ZZGFOCIN
C     OBSRVR     I   ZZGFOCIN
C     ABCORR     I   ZZGFOCIN
C     TIME       I   ZZGFOCST
C     OCSTAT     O   ZZGFOCST
C
C$ Detailed_Input
C
C     See entry points.
C
C$ Detailed_Output
C
C     See entry points.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     See entry points.
C
C$ Files
C
C
C     Appropriate SPICE kernels must be loaded by the calling program
C     before this routine is called.
C
C     The following data are required:
C
C        - SPK data: the calling application must load ephemeris data
C          for the target, source and observer that cover the time
C          period specified by the window CNFINE. If aberration
C          corrections are used, the states of the target bodies and of
C          the observer relative to the solar system barycenter must be
C          calculable from the available ephemeris data. Typically
C          ephemeris data
C          are made available by loading one or more SPK files via
C          FURNSH.
C
C        - PCK data: bodies modeled as triaxial ellipsoids must have
C          semi-axis lengths provided by variables in the kernel pool.
C          Typically these data are made available by loading a text
C          PCK file via FURNSH.
C
C        - FK data: if either of the reference frames designated by
C          BFRAME or FFRAME are not built in to the SPICE system,
C          one or more FKs specifying these frames must be loaded.
C
C     The following data may be required:
C
C        - DSK data: if either FSHAPE or BSHAPE indicates that DSK
C          data are to be used, DSK files containing topographic data
C          for the target body must be loaded. If a surface list is
C          specified, data for at least one of the listed surfaces must
C          be loaded.
C
C        - Surface name-ID associations: if surface names are specified
C          in FSHAPE or BSHAPE, the association of these names with
C          their corresponding surface ID codes must be established by
C          assignments of the kernel variables
C
C             NAIF_SURFACE_NAME
C             NAIF_SURFACE_CODE
C             NAIF_SURFACE_BODY
C
C          Normally these associations are made by loading a text
C          kernel containing the necessary assignments. An example
C          of such a set of assignments is
C
C             NAIF_SURFACE_NAME += 'Mars MEGDR 128 PIXEL/DEG'
C             NAIF_SURFACE_CODE += 1
C             NAIF_SURFACE_BODY += 499
C
C        - CK data: either of the body-fixed frames to which FFRAME or
C          BFRAME refer might be a CK frame. If so, at least one CK
C          file will be needed to permit transformation of vectors
C          between that frame and the J2000 frame.
C
C        - SCLK data: if a CK file is needed, an associated SCLK
C          kernel is required to enable conversion between encoded SCLK
C          (used to time-tag CK data) and barycentric dynamical time
C          (TDB).
C
C     Kernel data are normally loaded once per program run, NOT every
C     time this routine is called.
C
C$ Particulars
C
C     This routine is designed to determine whether a specified
C     type of occultation or transit is in progress at a specified
C     epoch. Three methods of modeling the shapes of the target
C     bodies are supported:
C
C        1)  Model both target bodies as triaxial ellipsoids. For this
C            case, the user may choose between occultations that are
C            partial, full or annular. See the entry header for
C            ZZGFOCIN for an explanation of these terms.
C
C        2)  Treat one target body as a point object and the other
C            target body as a triaxial ellipsoid. The only supported
C            occultation type is "ANY" for this case.
C
C        3)  Treat one target body as a point object and the other
C            target body as an extended shape modeled by DSK data.
C            The only supported occultation type is "ANY" for this
C            case.
C
C     This routine contains two entry points that support searches
C     for occultations performed using ZZGFSOLV:
C
C        ZZGFOCIN   Saves the user-supplied inputs defining the
C                   occultation computation to be performed.
C                   Initializes the occultation search.
C
C        ZZGFOCST   Returns the occultation state for a specified
C                   time.
C
C$ Examples
C
C     See GFOCCE.
C
C$ Restrictions
C
C     This is a SPICELIB private routine; it should not be called by
C     user applications.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     N.J. Bachman   (JPL)
C     L.S. Elson     (JPL)
C     W.L. Taber     (JPL)
C     I.M. Underwood (JPL)
C     E.D. Wright    (JPL)
C
C$ Version
C
C-    SPICELIB Version 2.0.0, 21-FEB-2017 (NJB)
C
C        Modified FAILED tests.
C
C        31-DEC-2016 (NJB)
C
C        Corrected long error message for blank frame string.
C        Previously the order of ERRCH and SIGERR calls was
C        inverted.
C
C        Updated long error message for incorrect body-fixed
C        frame center, so the ID of the target body is included.
C
C        24-FEB-2016 (NJB)
C
C        DSK capability was integrated.
C
C        04-MAR-2015 (NJB)
C
C        Initial updates to support surfaces represented by DSKs.
C
C-    SPICELIB Version 1.1.0, 18-MAY-2014 (NJB)
C
C        Bug fix: in entry point ZZGFOCIN, corrected long error message
C        for the case in which a body-fixed reference frame is not
C        centered on the correct body.
C
C-    SPICELIB Version 1.0.0, 05-MAR-2009 (NJB) (LSE) (WLT) (IMU) (EDW)
C
C-&

C
C     SPICELIB functions
C
      DOUBLE PRECISION      DASINE
      DOUBLE PRECISION      HALFPI
      DOUBLE PRECISION      VNORM
      DOUBLE PRECISION      VSEP

      INTEGER               ISRCHC
      INTEGER               ZZOCCED

      LOGICAL               EQSTR
      LOGICAL               FAILED
      LOGICAL               RETURN

C
C     Local parameters
C
C

C     ALPHA is a bound for the fraction of the speed of light
C     at which target body may move, relative to the solar
C     system barycenter.
C
      DOUBLE PRECISION      ALPHA
      PARAMETER           ( ALPHA = 1.D-2 )

C
C     ATOL is a tolerance value for computing arc sine.
C
      DOUBLE PRECISION      ATOL
      PARAMETER           ( ATOL  = 1.D-12 )

      INTEGER               BDNMLN
      PARAMETER           ( BDNMLN = 36 )

      INTEGER               FRNMLN
      PARAMETER           ( FRNMLN = 32 )

      INTEGER               POSLEN
      PARAMETER           ( POSLEN = 10 )

C
C     Local variables
C
      CHARACTER*(FRNMLN)    FIXFRM
      CHARACTER*(CVTLEN)    PNTDEF
      CHARACTER*(POSLEN)    POSNAM
      CHARACTER*(SHPLEN)    SHAPE
      CHARACTER*(SHPLEN)    SHPSTR
      CHARACTER*(SUBLEN)    SUBTYP
      CHARACTER*(FRNMLN)    SVBFRM
      CHARACTER*(MTHLEN)    SVBMTH
      CHARACTER*(BDNMLN)    SVBNAM
      CHARACTER*(SHPLEN)    SVBSHP
      CHARACTER*(CORLEN)    SVCORR
      CHARACTER*(FRNMLN)    SVFFRM
      CHARACTER*(MTHLEN)    SVFMTH
      CHARACTER*(BDNMLN)    SVFNAM
      CHARACTER*(SHPLEN)    SVFSHP
      CHARACTER*(BDNMLN)    SVONAM
      CHARACTER*(OCLLN)     SVTYPE
      CHARACTER*(OCLLN)     SVTYPS ( NOCTYP )
      CHARACTER*(TMTLEN)    TRMTYP

      DOUBLE PRECISION      BCKFRT ( 3 )
      DOUBLE PRECISION      BCKOBS ( 3 )
      DOUBLE PRECISION      BCKPOS ( 3 )
      DOUBLE PRECISION      BDIST
      DOUBLE PRECISION      BSMAXS ( 3, 3 )
      DOUBLE PRECISION      ETBCOR
      DOUBLE PRECISION      ETFCOR
      DOUBLE PRECISION      FDIST
      DOUBLE PRECISION      FRTBCK ( 3 )
      DOUBLE PRECISION      FRTOBS ( 3 )
      DOUBLE PRECISION      FRTPOS ( 3 )
      DOUBLE PRECISION      FSMAXS ( 3, 3 )
      DOUBLE PRECISION      LTBACK
      DOUBLE PRECISION      LTFRNT
      DOUBLE PRECISION      MAXANG
      DOUBLE PRECISION      MINANG
      DOUBLE PRECISION      MTEMP  ( 3, 3 )
      DOUBLE PRECISION      SVORIG ( 3 )
      DOUBLE PRECISION      RADII  ( 3 )
      DOUBLE PRECISION      SPOINT ( 3 )
      DOUBLE PRECISION      SRFVEC ( 3 )
      DOUBLE PRECISION      SRAD
      DOUBLE PRECISION      SVBRAD ( 3 )
      DOUBLE PRECISION      SVFRAD ( 3 )
      DOUBLE PRECISION      SVMNBR
      DOUBLE PRECISION      SVMNFR
      DOUBLE PRECISION      SVMXBR
      DOUBLE PRECISION      SVMXFR
      DOUBLE PRECISION      T2SEP
      DOUBLE PRECISION      TDIST
      DOUBLE PRECISION      TRGEPC
      DOUBLE PRECISION      TRGSEP

      INTEGER               CENTER
      INTEGER               CLSSID
      INTEGER               FFRMID
      INTEGER               FRCLSS
      INTEGER               I
      INTEGER               IDBACK
      INTEGER               IDFRNT
      INTEGER               IDOBS
      INTEGER               LOC
      INTEGER               N
      INTEGER               NSURF
      INTEGER               OCCNUM
      INTEGER               OCCODE
      INTEGER               SRFLST ( MAXSRF )
      INTEGER               SVBACK
      INTEGER               SVFRNT
      INTEGER               SVOBS
      INTEGER               TRGID

      LOGICAL               ATTBLK ( NABCOR )
      LOGICAL               FOUND
      LOGICAL               PNTOCC
      LOGICAL               PRI

      INTEGER               NCALLS
      SAVE                  NCALLS

C
C     Saved variables
C
      SAVE                  SVBACK
      SAVE                  SVBFRM
      SAVE                  SVBMTH
      SAVE                  SVBNAM
      SAVE                  SVBRAD
      SAVE                  SVBSHP
      SAVE                  SVCORR
      SAVE                  SVFFRM
      SAVE                  SVFMTH
      SAVE                  SVFNAM
      SAVE                  SVFRAD
      SAVE                  SVFRNT
      SAVE                  SVFSHP
      SAVE                  SVMNBR
      SAVE                  SVMNFR
      SAVE                  SVMXBR
      SAVE                  SVMXFR
      SAVE                  SVOBS
      SAVE                  SVONAM
      SAVE                  SVORIG
      SAVE                  SVTYPE
      SAVE                  SVTYPS

      DATA NCALLS / 0 /

C
C     Initial values
C
      DATA                  SVORIG / 3 * 0.D0 /

C
C     Below we initialize the list of occultation types.
C
      DATA                  SVTYPS / ANNULR, ANY, PARTL, FULL /

C
C     This routine should never be called directly.
C

      CALL CHKIN  ( 'ZZGFOCU' )

      CALL SIGERR ( 'SPICE(BOGUSENTRY)' )

      CALL CHKOUT ( 'ZZGFOCU' )
      RETURN

C$Procedure  ZZGFOCIN ( GF, occultation initialization )

      ENTRY ZZGFOCIN  ( OCCTYP, FRONT,  FSHAPE, FFRAME,
     .                  BACK,   BSHAPE, BFRAME, OBSRVR,
     .                  ABCORR                         )

C$ Abstract
C
C    Perform initialization functions for occultation state
C    determination.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     NAIF_IDS
C     FRAMES
C     PCK
C     SPK
C     TIME
C
C$ Keywords
C
C     SEARCH
C     GEOMETRY
C     OCCULTATION
C
C$ Declarations
C
C     CHARACTER*(*)         OCCTYP
C     CHARACTER*(*)         FRONT
C     CHARACTER*(*)         FSHAPE
C     CHARACTER*(*)         FFRAME
C     CHARACTER*(*)         BACK
C     CHARACTER*(*)         BSHAPE
C     CHARACTER*(*)         BFRAME
C     CHARACTER*(*)         OBSRVR
C     CHARACTER*(*)         ABCORR
C
C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     OCCTYP     I   Type of occultation.
C     FRONT      I   Name of body occulting the other.
C     FSHAPE     I   Type of shape model used for front body.
C     FFRAME     I   Body-fixed, body-centered frame for front body.
C     BACK       I   Name of body occulted by the other.
C     BSHAPE     I   Type of shape model used for back body.
C     BFRAME     I   Body-fixed, body-centered frame for back body.
C     OBSRVR     I   Name of the observing body.
C     ABCORR     I   Aberration correction flag.
C
C$ Detailed_Input
C
C
C     OCCTYP     indicates the type of occultation that is to be found.
C                The full set of possible values of OCCTYP may be used
C                when both target bodies are modeled as ellipsoids.
C                When either target is modeled as a point, OCCTYP must
C                be set to 'ANY' (see description below).
C
C                Supported values of OCCTYP and corresponding
C                definitions are:
C
C
C                   'FULL'      denotes the full occultation of the
C                               body designated by BACK by the body
C                               designated by FRONT, as seen from the
C                               location of the observer. In other
C                               words, the occulted body is completely
C                               invisible as seen from the observer's
C                               location.
C
C                   'ANNULAR'   denotes an annular occultation: the
C                               body designated by FRONT blocks part
C                               of, but not the limb of, the body
C                               designated by BACK, as seen from the
C                               location of the observer.
C
C                   'PARTIAL'   denotes a partial, non-annular
C                               occultation: the body designated by
C                               FRONT blocks part, but not all, of the
C                               limb of the body designated by BACK, as
C                               seen from the location of the observer.
C
C                   'ANY'       denotes any of the above three types of
C                               occultations: 'PARTIAL', 'ANNULAR', or
C                               'FULL'.
C
C                               'ANY' should be used to search for
C                               times when the body designated by FRONT
C                               blocks any part of the body designated
C                               by BACK.
C
C                               The option 'ANY' must be used if either
C                               the front or back target body is
C                               modeled as a point.
C
C                Case and leading or trailing blanks are not
C                significant in the string OCCTYP.
C
C
C     FRONT      is the name of the target body that occults---that is,
C                passes in front of---the other. Optionally, you may
C                supply the integer NAIF ID code for the body as a
C                string. For example both 'MOON' and '301' are
C                legitimate strings that designate the Moon.
C
C                Case and leading or trailing blanks are not
C                significant in the string FRONT.
C
C
C
C     FSHAPE     is a string indicating the geometric model used to
C                represent the shape of the front target body. The
C                supported options are:
C
C                   'ELLIPSOID'
C
C                       Use a triaxial ellipsoid model with radius
C                       values provided via the kernel pool. A kernel
C                       variable having a name of the form
C
C                          'BODYnnn_RADII'
C
C                       where nnn represents the NAIF integer code
C                       associated with the body, must be present in
C                       the kernel pool. This variable must be
C                       associated with three numeric values giving the
C                       lengths of the ellipsoid's X, Y, and Z
C                       semi-axes.
C
C                   'POINT'
C
C                       Treat the body as a single point. When a point
C                       target is specified, the occultation type must
C                       be set to 'ANY'.
C
C                   'DSK/UNPRIORITIZED[/SURFACES = <surface list>]'
C
C                       Use topographic data provided by DSK files to
C                       model the body's shape. These data must be
C                       provided by loaded DSK files.
C
C                       The surface list specification is optional. The
C                       syntax of the list is
C
C                          <surface 1> [, <surface 2>...]
C
C                       If present, it indicates that data only for the
C                       listed surfaces are to be used; however, data
C                       need not be available for all surfaces in the
C                       list. If absent, loaded DSK data for any surface
C                       associated with the target body are used.
C
C                       The surface list may contain surface names or
C                       surface ID codes. Names containing blanks must
C                       be delimited by double quotes, for example
C
C                          SURFACES = "Mars MEGDR 128 PIXEL/DEG"
C
C                       If multiple surfaces are specified, their names
C                       or IDs must be separated by commas.
C
C                       See the Particulars section below for details
C                       concerning use of DSK data.
C
C                The combinations of the shapes of the target bodies
C                FRONT and BACK must be one of:
C
C                   One ELLIPSOID, one POINT
C                   Two ELLIPSOIDs
C                   One DSK, one POINT
C
C                Case and leading or trailing blanks are not
C                significant in the string FSHAPE.
C
C
C     FFRAME     is the name of the body-fixed, body-centered reference
C                frame associated with the front target body. Examples
C                of such names are 'IAU_SATURN' (for Saturn) and
C                'ITRF93' (for the Earth).
C
C                If the front target body is modeled as a point, FFRAME
C                should be left blank.
C
C                Case and leading or trailing blanks bracketing a
C                non-blank frame name are not significant in the string
C                FFRAME.
C
C
C     BACK       is the name of the target body that is occulted
C                by---that is, passes in back of---the other.
C                Optionally, you may supply the integer NAIF ID code
C                for the body as a string. For example both 'MOON' and
C                '301' are legitimate strings that designate the Moon.
C
C                Case and leading or trailing blanks are not
C                significant in the string BACK.
C
C
C     BSHAPE     is the shape specification for the body designated
C                by BACK. See the description of FSHAPE above for
C                details.
C
C
C     BFRAME     is the name of the body-fixed, body-centered reference
C                frame associated with the ``back'' target body.
C                Examples of such names are 'IAU_SATURN' (for Saturn)
C                and 'ITRF93' (for the Earth).
C
C                If the back target body is modeled as a point, BFRAME
C                should be left blank.
C
C                Case and leading or trailing blanks bracketing a
C                non-blank frame name are not significant in the string
C                BFRAME.
C
C
C     OBSRVR     is the name of the body from which the occultation is
C                observed. Optionally, you may supply the integer NAIF
C                ID code for the body as a string.
C
C                Case and leading or trailing blanks are not
C                significant in the string OBSRVR.
C
C
C     ABCORR     indicates the aberration corrections to be applied to
C                the state of the target body to account for one-way
C                light time.  Stellar aberration corrections are
C                ignored if specified, since these corrections don't
C                improve the accuracy of the occultation determination.
C
C                See the header of the SPICE routine SPKEZR for a
C                detailed description of the aberration correction
C                options. For convenience, the options supported by
C                this routine are listed below:
C
C                   'NONE'     Apply no correction.
C
C                   'LT'       "Reception" case:  correct for
C                              one-way light time using a Newtonian
C                              formulation.
C
C                   'CN'       "Reception" case:  converged
C                              Newtonian light time correction.
C
C                   'XLT'      "Transmission" case:  correct for
C                              one-way light time using a Newtonian
C                              formulation.
C
C                   'XCN'      "Transmission" case:  converged
C                              Newtonian light time correction.
C
C                Case and blanks are not significant in the string
C                ABCORR.
C
C$ Detailed_Output
C
C     None.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1)  If name of either target or the observer cannot be translated
C         to a NAIF ID code, the error SPICE(IDCODENOTFOUND) is
C         signaled.
C
C     2)  If either of the target bodies FRONT or BACK coincides with
C         the observer body OBSRVR, or if the targets coincide,
C         the error SPICE(BODIESNOTDISTINCT) will be signaled.
C
C     3)  If either of the body model specifiers FSHAPE or BSHAPE
C         is not recognized, the error SPICE(INVALIDSHAPE) will be
C         signaled.
C
C     4)  If both of the body model specifiers FSHAPE and BSHAPE
C         specify point targets, the error SPICE(INVALIDSHAPECOMBO)
C         will be signaled.
C
C     5)  If one of the target shape arguments FSHAPE and BSHAPE
C         specifies a DSK model, and the other argument does not
C         specify a point target, the error SPICE(INVALIDSHAPECOMBO)
C         will be signaled.
C
C     6)  If an unrecognized value of OCCTYP is seen, the error
C         SPICE(INVALIDOCCTYPE) is signaled.
C
C     7)  If one target body is modeled as a point and OCCTYP is not
C         set to 'ANY', the error SPICE(BADTYPESHAPECOMBO) is signaled.
C
C     8)  If a target indicated to be an ellipsoid by its shape
C         specification argument does not have three associated
C         positive radii, the error SPICE(DEGENERATECASE) will be
C         signaled.
C
C     9)  If the number of radii associated with a target body is
C         not three, the error SPICE(BADRADIUSCOUNT) will be
C         signaled.
C
C     10) If a target body-fixed reference frame associated with a
C         non-point target is not recognized, the error
C         SPICE(INVALIDFRAME) will be signaled.
C
C     11) If a target body-fixed reference frame is not centered at
C         the corresponding target body, the error
C         SPICE(INVALIDFRAME) will be signaled.
C
C     12) If the aberration correction string is invalid, the error
C         will be diagnosed by a routine in the call tree of this
C         routine.
C
C     13) If either FSHAPE or BSHAPE specifies that the target surface
C         is represented by DSK data, and no DSK files are loaded for
C         the specified target, the error is signaled by a routine in
C         the call tree of this routine.
C
C     14) If either FSHAPE or BSHAPE specifies that the target surface
C         is represented by DSK data, but the shape specification is
C         invalid, the error is signaled by a routine in the call tree
C         of this routine.
C$ Files
C
C     See the header of the umbrella routine ZZGFOCU.
C
C$ Particulars
C
C     This entry point initializes the parameters needed by the
C     occultation state determination entry point ZZGFOCST.
C
C$ Examples
C
C     See implementation of GFOCCE.
C
C$ Restrictions
C
C     This is a SPICELIB private routine; it should not be called by
C     user applications.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     L.S. Elson     (JPL)
C     W.L. Taber     (JPL)
C     N.J. Bachman   (JPL)
C     E.D. Wright    (JPL)
C
C$ Version
C
C-    SPICELIB Version 2.0.0, 21-FEB-2017 (NJB)
C
C        Modified FAILED tests.
C
C        20-FEB-2016 (NJB)
C
C        DSK capability was integrated.
C
C        04-MAR-2015 (NJB)
C
C        Initial updates to support surfaces represented by DSKs.
C
C-    SPICELIB Version 1.1.0, 18-MAY-2014 (NJB)
C
C        Bug fix: corrected long error message for the case in which a
C        body-fixed reference frame is not centered on the correct
C        body.
C
C-    SPICELIB Version 1.0.0, 15-APR-2009 (LSE) (WLT) (NJB) (EDW)
C
C-&

C$ Revisions
C
C     None.
C
C-&

C
C     Standard SPICE error handling.
C
      IF ( RETURN () ) THEN
         RETURN
      END IF

      CALL CHKIN  ( 'ZZGFOCIN' )

C
C     Find NAIF IDs for FRONT, BACK, and OBSRVR.
C
      CALL BODS2C ( FRONT, IDFRNT, FOUND )

      IF ( .NOT. FOUND ) THEN

         CALL SETMSG ( 'The front target object, '
     .   //            '''#'', is not a recognized name for an '
     .   //            'ephemeris object. The cause of this '
     .   //            'problem may be that you need an updated '
     .   //            'version of the SPICE toolkit. '           )
         CALL ERRCH  ( '#', FRONT                                 )
         CALL SIGERR ( 'SPICE(IDCODENOTFOUND)'                    )
         CALL CHKOUT ( 'ZZGFOCIN'                                 )
         RETURN

      END IF

      CALL BODS2C ( BACK, IDBACK, FOUND )

      IF ( .NOT. FOUND ) THEN

         CALL SETMSG ( 'The back target object, '
     .   //            '''#'', is not a recognized name for an '
     .   //            'ephemeris object. The cause of this '
     .   //            'problem may be that you need an updated '
     .   //            'version of the SPICE toolkit. '           )
         CALL ERRCH  ( '#', BACK                                  )
         CALL SIGERR ( 'SPICE(IDCODENOTFOUND)'                    )
         CALL CHKOUT ( 'ZZGFOCIN'                                 )
         RETURN

      END IF

      CALL BODS2C ( OBSRVR, IDOBS, FOUND )

      IF ( .NOT. FOUND ) THEN

         CALL SETMSG ( 'The observer, '
     .   //            '''#'', is not a recognized name for an '
     .   //            'ephemeris object. The cause of this '
     .   //            'problem may be that you need an updated '
     .   //            'version of the SPICE toolkit. '           )
         CALL ERRCH  ( '#', OBSRVR                                )
         CALL SIGERR ( 'SPICE(IDCODENOTFOUND)'                    )
         CALL CHKOUT ( 'ZZGFOCIN'                                 )
         RETURN

      END IF

C
C     Make sure the observer and both targets are distinct.
C
      IF (      ( IDFRNT .EQ. IDBACK )
     .     .OR. ( IDFRNT .EQ. IDOBS  )
     .     .OR. ( IDBACK .EQ. IDOBS  ) ) THEN

         CALL SETMSG ( 'The observer and both targets must be '
     .   //            'distinct objects, but are not: '
     .   //            'OBSRVR = #; FRONT = #; BACK = #.'      )
         CALL ERRCH  ( '#', OBSRVR                             )
         CALL ERRCH  ( '#', FRONT                              )
         CALL ERRCH  ( '#', BACK                               )
         CALL SIGERR ( 'SPICE(BODIESNOTDISTINCT)'              )
         CALL CHKOUT ( 'ZZGFOCIN'                              )
         RETURN

      END IF

C
C     Save the objects' names. We'll need these if
C     we need to call SINCPT.
C
      SVFNAM = FRONT
      SVBNAM = BACK
      SVONAM = OBSRVR

C
C     Store the ID codes, shape specifications, and body-fixed,
C     body-centered frame names of the objects involved in this event.
C     The shape arguments must be parsed in case they contain
C     DSK specifications.
C
      SVFRNT = IDFRNT
      SVFFRM = FFRAME

      SVBACK = IDBACK
      SVBFRM = BFRAME

      SVOBS  = IDOBS

C
C     Save the input shape strings. These will be examined later,
C     but we need them in their original form for computations
C     involving DSK data. In the variable names below, "MTH"
C     stands for "method"---the name used in SPICE geometry
C     APIs for this type of input string.
C
      SVFMTH = FSHAPE
      SVBMTH = BSHAPE

C
C     Parse the front body shape string.
C
      IF (  EQSTR( FSHAPE, 'POINT' )  ) THEN

         SVFSHP = PTSHAP

      ELSE IF (  EQSTR( FSHAPE, 'ELLIPSOID' )  ) THEN

         SVFSHP = EDSHAP

      ELSE

         CALL ZZPRSMET ( IDFRNT, SVFMTH, MAXSRF, SHPSTR, SUBTYP,
     .                   PRI,    NSURF,  SRFLST, PNTDEF, TRMTYP )

         IF ( FAILED() ) THEN
            CALL CHKOUT ( 'ZZGFOCIN' )
            RETURN
         END IF

         IF (  EQSTR( SHPSTR, 'DSK' )  ) THEN

            SVFSHP = DSSHAP

         ELSE

            CALL SETMSG ( 'Front target shape from FSHAPE '
     .      //            'string was <#>. Valid shapes are '
     .      //            'ELLIPSOID, POINT, and DSK.'       )
            CALL ERRCH  ( '#', FSHAPE                        )
            CALL SIGERR ( 'SPICE(INVALIDSHAPE)'              )
            CALL CHKOUT ( 'ZZGFOCIN'                         )
            RETURN

         END IF

      END IF

C
C     Parse the back body shape string.
C
      IF (  EQSTR( BSHAPE, 'POINT' )  ) THEN

         SVBSHP = PTSHAP

      ELSE IF (  EQSTR( BSHAPE, 'ELLIPSOID' )  ) THEN

         SVBSHP = EDSHAP

      ELSE

         CALL ZZPRSMET ( IDFRNT, SVBMTH, MAXSRF, SHPSTR, SUBTYP,
     .                   PRI,    NSURF,  SRFLST, PNTDEF, TRMTYP )

         IF ( FAILED() ) THEN
            CALL CHKOUT ( 'ZZGFOCIN' )
            RETURN
         END IF

         IF (  EQSTR( SHPSTR, 'DSK' )  ) THEN

            SVBSHP = DSSHAP

         ELSE

            CALL SETMSG ( 'Back target shape from BSHAPE '
     .      //            'string was <#>. Valid shapes are '
     .      //            'ELLIPSOID, POINT, and DSK.'       )
            CALL ERRCH  ( '#', BSHAPE                        )
            CALL SIGERR ( 'SPICE(INVALIDSHAPE)'              )
            CALL CHKOUT ( 'ZZGFOCIN'                         )
            RETURN

         END IF

      END IF

C
C     Check for invalid shape combinations.
C
      IF (       ( SVFSHP .EQ. PTSHAP )
     .     .AND. ( SVBSHP .EQ. PTSHAP ) ) THEN

         CALL SETMSG ( 'Both front and back objects have '
     .   //            'POINT shape specifications; only '
     .   //            'one point shape is allowed. The '
     .   //            'other shape must be ELLIPSOID or '
     .   //            'DSK.'                             )
         CALL SIGERR ( 'SPICE(INVALIDSHAPECOMBO)'         )
         CALL CHKOUT ( 'ZZGFOCIN'                         )
         RETURN

      ELSE IF (      (        ( SVFSHP .EQ. DSSHAP )
     .                  .AND. ( SVBSHP .NE. PTSHAP ) )
     .          .OR. (        ( SVBSHP .EQ. DSSHAP )
     .                  .AND. ( SVFSHP .NE. PTSHAP ) )  ) THEN

         CALL SETMSG ( 'Front target shape from FSHAPE '
     .   //            'string was <#>; back target shape '
     .   //            'from BSHAPE was <#>. When one shape '
     .   //            'is DSK, the other must be POINT.'   )
         CALL ERRCH  ( '#', SVFSHP                          )
         CALL ERRCH  ( '#', SVBSHP                          )
         CALL SIGERR ( 'SPICE(INVALIDSHAPECOMBO)'           )
         CALL CHKOUT ( 'ZZGFOCIN'                           )
         RETURN

       END IF

C
C     Save a single upper-case character representing the occultation
C     type string.
C
      CALL LJUST ( OCCTYP, SVTYPE )
      CALL UCASE ( SVTYPE, SVTYPE )

C
C     Check the occultation type.
C
      OCCNUM = ISRCHC ( SVTYPE, NOCTYP, SVTYPS )

      IF ( OCCNUM .EQ. 0 ) THEN

         CALL SETMSG ( 'The occultation type # is not '
     .   //            'recognized.  Supported types are: '
     .   //            '#, #, #,  #.'                      )
         CALL ERRCH  ('#', OCCTYP                          )

         DO I = 1, NOCTYP
            CALL ERRCH  ('#', SVTYPS(I)                    )
         END DO

         CALL SIGERR ('SPICE(INVALIDOCCTYPE)'              )
         CALL CHKOUT ('ZZGFOCIN'                           )
         RETURN

      END IF

C
C     If we have a point target, the occultation type must
C     be 'ANY'.
C
      IF (  ( SVFSHP .EQ. PTSHAP ) .OR. ( SVBSHP .EQ. PTSHAP )  ) THEN

         IF ( SVTYPE .NE. ANY ) THEN

            CALL SETMSG ( 'Occultation type # is not allowed when '
     .      //            'either target body is modeled as a '
     .      //            'point. Set OCCTYP to ANY for use with '
     .      //            'point targets.'                         )
            CALL ERRCH  ( '#',  OCCTYP                             )
            CALL SIGERR ( 'SPICE(BADTYPESHAPECOMBO)'               )
            CALL CHKOUT ( 'ZZGFOCIN'                               )
            RETURN

         END IF

      END IF

C
C     Check the aberration correction. If SPKEZR can't handle it,
C     neither can we.
C
      CALL ZZVALCOR ( ABCORR, ATTBLK )

      IF ( FAILED() ) THEN
         CALL CHKOUT ( 'ZZGFOCIN' )
         RETURN
      END IF

C
C     Create a local aberration correction string without
C     a stellar aberration correction specifier.
C
      IF ( ATTBLK(GEOIDX) ) THEN

         SVCORR = 'NONE'

      ELSE
C
C        The correction string specified either Newtonian or converged
C        light time correction.
C
         IF ( ATTBLK(XMTIDX) ) THEN
            SVCORR = 'X'
         ELSE
            SVCORR = ' '
         END IF

         IF ( ATTBLK(CNVIDX) ) THEN
            CALL SUFFIX ( 'CN', 0, SVCORR )
         ELSE
            CALL SUFFIX ( 'LT', 0, SVCORR )
         END IF

      END IF

C
C     Check the front and back targets' shapes, frames
C     and radii.
C
      DO I = 1, 2

         IF ( I .EQ. 1 ) THEN

            POSNAM = 'front'
            FIXFRM = FFRAME
            TRGID  = IDFRNT
            SHAPE  = SVFSHP

         ELSE

            POSNAM = 'back'
            FIXFRM = BFRAME
            TRGID  = IDBACK
            SHAPE  = SVBSHP

         END IF

         IF ( SHAPE .EQ. EDSHAP ) THEN
C
C           Fetch and check the radii.
C
            CALL BODVCD ( TRGID, 'RADII', 3, N, RADII )

            IF ( FAILED() ) THEN
               CALL CHKOUT ( 'ZZGFOCIN' )
               RETURN
            END IF

C
C           Check the count of the radii.
C
            IF ( N .NE. 3 )  THEN

               CALL SETMSG ( 'Target # should have 3 radii but '
     .         //            'actually has #. This may be due '
     .         //            'to an error in a PCK file used '
     .         //            'to provide the radii.'             )
               CALL ERRCH  ( '#', POSNAM                         )
               CALL ERRINT ( '#', N                              )
               CALL SIGERR ( 'SPICE(BADRADIUSCOUNT)'             )
               CALL CHKOUT ( 'ZZGFOCIN'                          )
               RETURN

            END IF

C
C           Check to make sure the current target has 3 positive
C           semi-axis lengths.
C
            IF (     ( RADII(1) .LE. 0.D0 )
     .          .OR. ( RADII(2) .LE. 0.D0 )
     .          .OR. ( RADII(3) .LE. 0.D0 ) ) THEN

               CALL SETMSG ( 'One or more semi-axis lengths of '
     .         //            'the # target body are non-'
     .         //            'positive: 1 = #, 2 = #, 3 = #. '  )
               CALL ERRCH  ( '#', POSNAM                        )
               CALL ERRDP  ( '#', RADII (1)                     )
               CALL ERRDP  ( '#', RADII (2)                     )
               CALL ERRDP  ( '#', RADII (3)                     )
               CALL SIGERR ( 'SPICE(DEGENERATECASE)'             )
               CALL CHKOUT ( 'ZZGFOCIN'                          )
               RETURN

            END IF

C
C           Checks of radii have been completed.
C
            IF ( I .EQ. 1 ) THEN

               CALL MOVED ( RADII,  3,  SVFRAD )
C
C              Select smallest and largest semi-axis lengths of body
C              for later tests.
C
               CALL MINAD ( SVFRAD, 3,  SVMNFR, LOC )
               CALL MAXAD ( SVFRAD, 3,  SVMXFR, LOC )

            ELSE

               CALL MOVED ( RADII,  3,  SVBRAD )

               CALL MINAD ( SVBRAD, 3,  SVMNBR, LOC )
               CALL MAXAD ( SVBRAD, 3,  SVMXBR, LOC )

            END IF

            IF ( FAILED() ) THEN
               CALL CHKOUT ( 'ZZGFOCIN' )
               RETURN
            END IF

         END IF
C
C        We've performed radii checks for an ellipsoidal target.
C        Minimum and maximum bounding radii are set, if the target
C        shape is modeled as an ellipsoid.
C
C
C        Check body-fixed frame for extended targets.
C
         IF (  ( SHAPE .EQ. EDSHAP ) .OR. ( SHAPE .EQ. DSSHAP )  ) THEN
C
C           The target is ellipsoidal or is modeled using DSK data;
C           there must be a target body-fixed frame associated with
C           this body.
C
            IF ( FIXFRM .EQ. ' ' ) THEN

               CALL SETMSG ( 'The # target shape is represented '
     .         //            'by an ellipsoid or by DSK data, '
     .        //             'but the associated '
     .         //            'body-fixed frame name is blank.' )
               CALL ERRCH  ( '#',  POSNAM                      )
               CALL SIGERR ( 'SPICE(INVALIDFRAME)'             )
               CALL CHKOUT ( 'ZZGFOCIN'                        )
               RETURN

            ELSE
C
C              Look up the target's body-fixed frame ID code.
C
               CALL NAMFRM ( FIXFRM, FFRMID )

               IF ( FAILED() ) THEN
                  CALL CHKOUT ( 'ZZGFOCIN' )
                  RETURN
               END IF

               IF ( FFRMID .EQ. 0 ) THEN

                  CALL SETMSG ( 'The # target''s body-fixed frame '
     .            //            'name # is not recognized.'        )
                  CALL ERRCH  ( '#',  POSNAM                       )
                  CALL ERRCH  ( '#',  FIXFRM                       )
                  CALL SIGERR ( 'SPICE(INVALIDFRAME)'              )
                  CALL CHKOUT ( 'ZZGFOCIN'                         )
                  RETURN

               END IF

C
C              Obtain the center of the frame and verify it's the
C              Ith target.
C
               CALL FRINFO ( FFRMID, CENTER, FRCLSS, CLSSID, FOUND )

               IF ( FAILED() ) THEN
                  CALL CHKOUT ( 'ZZGFOCIN' )
                  RETURN
               END IF

               IF ( .NOT. FOUND ) THEN
C
C                 Since we mapped the frame name to an ID code, we
C                 expect to find the frame info. So control should
C                 never reach this point.
C
                  CALL SETMSG ( 'Frame ID found for # body-fixed '
     .            //            'frame # but FRINFO couldn''t '
     .            //            'find frame info.'                 )
                  CALL ERRCH  ( '#',  POSNAM                       )
                  CALL ERRCH  ( '#',  FIXFRM                       )
                  CALL SIGERR ( 'SPICE(BUG)'                       )
                  CALL CHKOUT ( 'ZZGFOCIN'                         )
                  RETURN

               END IF

               IF ( CENTER .NE. TRGID ) THEN
C
C                 The body-fixed frame for the current target
C                 isn't actually centered on the body.
C
                  CALL SETMSG ( 'Supposed body-fixed frame # for '
     .            //            '# target # is actually centered '
     .            //            'on body #.'                       )
                  CALL ERRCH  ( '#',  FIXFRM                       )
                  CALL ERRCH  ( '#',  POSNAM                       )
                  CALL ERRINT ( '#',  TRGID                        )
                  CALL ERRINT ( '#',  CENTER                       )
                  CALL SIGERR ( 'SPICE(INVALIDFRAME)'              )
                  CALL CHKOUT ( 'ZZGFOCIN'                         )
                  RETURN

               END IF

            END IF

         END IF
C
C        We've performed frame checks for an extended target.
C
C
C        Obtain radii of inner and outer bounding spheres for
C        DSK targets.
C
         IF ( SHAPE .EQ. DSSHAP ) THEN
C
C           Note that TRGID and FFRMID refer to the current
C           target (out of two); "FFRMID" means "fixed frame ID."
C
            CALL ZZSUDSKI ( TRGID, NSURF, SRFLST, FFRMID )

            IF ( FAILED() ) THEN
               CALL CHKOUT ( 'ZZGFOCIN' )
               RETURN
            END IF

            IF ( I .EQ. 1 ) THEN

               CALL ZZMINRAD ( SVMNFR )
               CALL ZZMAXRAD ( SVMXFR )
            ELSE
               CALL ZZMINRAD ( SVMNBR )
               CALL ZZMAXRAD ( SVMXBR )
            END IF

         END IF

C
C        Initialize bounding radii and body-fixed frame
C        names for point targets.
C
         IF ( SHAPE .EQ. PTSHAP ) THEN
C
C           Zero out radius values for this target; set the
C           frame to blank.
C
            IF ( I .EQ. 1 ) THEN

               CALL CLEARD ( 3, SVFRAD )

               SVMNFR = 0.D0
               SVMXFR = 0.D0
               SVFFRM = ' '

            ELSE

               CALL CLEARD ( 3, SVBRAD )

               SVMNBR = 0.D0
               SVMXBR = 0.D0
               SVBFRM = ' '

            END IF

         END IF
C
C        We've performed shape, and if applicable, frame and radii
C        checks for the Ith target. Bounding radii have been obtained
C        for extended targets.
C
      END DO
C
C     We've performed shape, and if applicable, frame and radii
C     checks for both targets.
C
      CALL CHKOUT ( 'ZZGFOCIN' )
      RETURN

C$Procedure ZZGFOCST ( GF, "in occultation?"  )

      ENTRY ZZGFOCST ( TIME, OCSTAT )

C$ Abstract
C
C     See if the object is currently occulted.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     SPK
C     TIME
C
C$ Keywords
C
C     SEARCH
C     GEOMETRY
C     OCCULTATION
C
C$ Declarations
C
C     DOUBLE PRECISION      TIME
C     LOGICAL               OCSTAT
C
C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     TIME       I   TDB epoch (in seconds past J2000)
C     OCSTAT     O   .TRUE. if the object is occulted, .FALSE.
C                    otherwise.
C
C$ Detailed_Input
C
C     TIME       is the epoch of interest in TDB seconds past the
C                J2000 epoch.
C
C$ Detailed_Output
C
C     OCSTAT     is a logical flag indicating the state of
C                occultation. If the configuration initialized by
C                ZZGFOCIN is in occultation at the epoch TIME, OCSTAT is
C                returned with the value .TRUE. Otherwise it is
C                returned with the value .FALSE.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1) If any SPK lookup fails, the error will be diagnosed by
C        routines in the call tree of this routine.
C
C     2) If any frame transformation lookup fails, the error will be
C        diagnosed by routines in the call tree of this routine.
C
C     3) If any occultation computation is done for ellipsoidal
C        targets, and if either semi-axis matrix is invalid, the error
C        will be diagnosed by routines in the call tree of this
C        routine.
C
C     4) If any two of the bodies defining the occultation geometry
C        intersect, either error SPICE(NOTDISJOINT) will be
C        signaled by this routine, or the error will be diagnosed by
C        routines in the call tree of this routine.
C
C     5)  If the body model specifiers FSHAPE and BSHAPE don't specify
C         either two ellipsoidal targets or one ellipsoidal target and
C         one point target, the error SPICE(INVALIDSHAPECOMBO)
C         will be signaled.
C
C$ Files
C
C     See the Files header section of the umbrella routine ZZGFOCU.
C
C$ Particulars
C
C     This routine determines the occultation state of the
C     configuration specified by the last call to ZZGFOCIN and the
C     input time value.
C
C$ Examples
C
C     See the umbrella routine ZZGFOCU.
C
C$ Restrictions
C
C     This is a SPICELIB private routine; it should not be called by
C     user applications.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     N.J. Bachman   (JPL)
C     L.S. Elson     (JPL)
C     W.L. Taber     (JPL)
C     E.D. Wright    (JPL)
C
C$ Version
C
C-    SPICELIB Version 2.0.0, 21-FEB-2017 (NJB)
C
C        Modified FAILED tests.
C
C        20-FEB-2016 (NJB)
C
C        DSK capability was integrated.
C
C        04-MAR-2015 (NJB)
C
C        Initial updates to support surfaces represented by DSKs.
C
C-    SPICELIB Version 1.0.0, 30-DEC-2008 (NJB) (LSE) (WLT) (EDW)
C
C-&

C
C     Standard SPICE error handling.
C
      IF ( RETURN () ) THEN
         RETURN
      END IF

      CALL CHKIN  ( 'ZZGFOCST' )

C
C     Initialize the state output.
C
      OCSTAT = .FALSE.

C
C     Get the apparent positions of FRONT and BACK as seen from the
C     observer.
C
      CALL SPKEZP ( SVFRNT,  TIME,   'J2000',  SVCORR,
     .              SVOBS,   FRTPOS, LTFRNT           )

      CALL SPKEZP ( SVBACK,  TIME,   'J2000',  SVCORR,
     .              SVOBS,   BCKPOS, LTBACK           )

      IF ( FAILED() ) THEN
         CALL CHKOUT ( 'ZZGFOCST' )
         RETURN
      END IF

C
C     Handle the cases of one and two extended targets
C     separately.
C
      IF (       ( SVBSHP .EQ. EDSHAP )
     .     .AND. ( SVFSHP .EQ. EDSHAP )  ) THEN
C
C        The caller has selected a test for a partial, annular or full
C        occultation using ellipsoidal shape models.
C
C        Look up the axes of each target body in the J2000 frame at the
C        light time corrected epoch for that body.
C
         CALL ZZCOREPC ( SVCORR,  TIME,   LTBACK, ETBCOR )

         CALL PXFORM   ( SVBFRM, 'J2000', ETBCOR, MTEMP  )

         IF ( FAILED() ) THEN
            CALL CHKOUT ( 'ZZGFOCST' )
            RETURN
         END IF
C
C        Scale the columns of MTEMP by the axis lengths of the back
C        target.
C
         DO I = 1, 3
            CALL VSCL ( SVBRAD(I), MTEMP(1,I), BSMAXS(1,I) )
         END DO

         CALL ZZCOREPC ( SVCORR,  TIME,   LTFRNT, ETFCOR )
         CALL PXFORM   ( SVFFRM, 'J2000', ETFCOR, MTEMP  )

         IF ( FAILED() ) THEN
            CALL CHKOUT ( 'ZZGFOCST' )
            RETURN
         END IF
C
C        Scale the columns of MTEMP by the axis lengths of the second
C        target.
C
         DO I = 1, 3
            CALL VSCL ( SVFRAD(I), MTEMP(1,I), FSMAXS(1,I) )
         END DO

C
C        Classify the occultation state of BACK by FRONT as seen from
C        the observer.
C
         OCCODE = ZZOCCED ( SVORIG, BCKPOS, BSMAXS, FRTPOS, FSMAXS )

         IF ( FAILED() ) THEN
            CALL CHKOUT ( 'ZZGFOCST' )
            RETURN
         END IF

         IF ( OCCODE .EQ. NOOCC ) THEN
C
C           Neither body occults the other.
C
            OCSTAT = .FALSE.

         ELSE IF (       ( SVTYPE .EQ. ANY )
     .             .AND. ( OCCODE .LT. 0   ) ) THEN
C
C           The "of" body (target 1) is at least partially occulted by
C           the BY object.
C
            OCSTAT = .TRUE.

         ELSE IF (       ( SVTYPE .EQ.  FULL   )
     .             .AND. ( OCCODE .EQ.  TOTAL1 ) ) THEN
C
C           The BACK body is in total occultation.
C

            OCSTAT = .TRUE.

         ELSE IF (       ( SVTYPE .EQ. ANNULR )
     .             .AND. ( OCCODE .EQ. ANNLR1 ) ) THEN
C
C           The  BACK body is in annular occultation.
C
            OCSTAT = .TRUE.

         ELSE IF (       ( SVTYPE .EQ. PARTL  )
     .             .AND. ( OCCODE .EQ. PARTL1 ) ) THEN
C
C           The BACK body is partially occulted.
C
            OCSTAT = .TRUE.

         ELSE
C
C           The occultation state doesn't match the requested state.
C
            OCSTAT = .FALSE.

         END IF

         CALL CHKOUT ( 'ZZGFOCST' )
         RETURN

      ELSE IF (     (       ( SVFSHP .EQ. EDSHAP )
     .                .AND. ( SVBSHP .EQ. PTSHAP ) )
     .         .OR. (       ( SVFSHP .EQ. DSSHAP )
     .                .AND. ( SVBSHP .EQ. PTSHAP ) )
     .         .OR. (       ( SVFSHP .EQ. PTSHAP )
     .                .AND. ( SVBSHP .EQ. EDSHAP ) )
     .         .OR. (       ( SVFSHP .EQ. PTSHAP )
     .                .AND. ( SVBSHP .EQ. DSSHAP ) )  ) THEN
C
C        One of the targets is modeled as a point; the other is
C        modeled as an ellipsoid or a DSK shape.
C
C        If the front target is an ellipsoid or a DSK shape and the
C        back target is a point, we'll classify the geometry as a
C        "point occultation." Otherwise we have a "point transit" case.
C        We'll set the logical flag PNTOCC to .TRUE. to indicate a
C        point occultation.
C
         PNTOCC = SVBSHP .EQ. PTSHAP

C
C        We're going to start out by doing some error checking.
C        We're looking for intersections of the participating
C        objects: these should never occur.
C
C        Let BDIST, FDIST be the distances from the observer
C        to the back and front targets, respectively.
C
         BDIST  =  VNORM ( BCKPOS )
         FDIST  =  VNORM ( FRTPOS )

C
C        Find the vector from BACK to FRONT.  We'll use this later,
C        but we want it now in order to make sure that BACK doesn't
C        intersect FRONT.
C
         CALL VSUB ( FRTPOS, BCKPOS, BCKFRT )

         IF ( PNTOCC ) THEN
C
C           The front target is an extended shape.
C
            IF ( FDIST .LE. SVMNFR ) THEN
C
C              The observer is INSIDE the front target. We
C              treat this as an error.
C
               CALL SETMSG ( 'Observer is inside front target body.' )
               CALL SIGERR ( 'SPICE(NOTDISJOINT)'                    )
               CALL CHKOUT ( 'ZZGFOCST'                              )
               RETURN

            ELSE IF ( BDIST .EQ. 0.D0 ) THEN

               CALL SETMSG ( 'Back target coincides with observer.'  )
               CALL SIGERR ( 'SPICE(NOTDISJOINT)'                    )
               CALL CHKOUT ( 'ZZGFOCST'                              )
               RETURN

            ELSE IF ( VNORM(BCKFRT) .LE. SVMNFR ) THEN

               CALL SETMSG ( 'BACK target is inside FRONT target.' )
               CALL SIGERR ( 'SPICE(NOTDISJOINT)'                  )
               CALL CHKOUT ( 'ZZGFOCST'                            )
               RETURN

            END IF

         ELSE
C
C           The back target is an extended shape.
C
            IF ( BDIST .LE. SVMNBR ) THEN
C
C              The observer is INSIDE the back target. We
C              treat this as an error.
C
               CALL SETMSG ( 'Observer is inside back target body.' )
               CALL SIGERR ( 'SPICE(NOTDISJOINT)'                   )
               CALL CHKOUT ( 'ZZGFOCST'                             )
               RETURN

            ELSE IF ( FDIST .EQ. 0.D0 ) THEN

               CALL SETMSG ( 'Front target coincides with observer.'  )
               CALL SIGERR ( 'SPICE(NOTDISJOINT)'                     )
               CALL CHKOUT ( 'ZZGFOCST'                               )
               RETURN

            ELSE IF ( VNORM(BCKFRT) .LE. SVMNBR ) THEN

               CALL SETMSG ( 'FRONT target is inside BACK target.' )
               CALL SIGERR ( 'SPICE(NOTDISJOINT)'                  )
               CALL CHKOUT ( 'ZZGFOCST'                            )
               RETURN

            END IF

         END IF

C
C        Find angular separation of the target centers as
C        seen by the observer.
C
         TRGSEP = VSEP ( BCKPOS, FRTPOS )

C
C        Find angular radius of the outer bounding sphere of the
C        extended target, as seen by the observer.
C
C        In computing this angular radius, scale up the bounding
C        sphere to compensate for the light time error we've made
C        by computing light time to the target's center. The
C        correct value to use is light time to the limb point having
C        minimum angular separation from the point target.
C
C        Presuming the extended target can move no faster than
C        alpha*c (where c represents the speed of light in a vacuum),
C        and considering the fact that the light time error cannot
C        exceed r/c, where r is the radius of the outer bounding sphere
C        of the ellipsoid, we find that the magnitude of the position
C        error of the extended target cannot exceed alpha*r. Then the
C        correctly positioned target---that is, located at
C        the position corresponding to the correct light time
C        correction---must be contained in the outer bounding
C        sphere we've found, if we scale the sphere up by 1+alpha.
C
C        Perform the test only if the observer is outside the
C        outer bounding sphere of the extended target.
C
         IF ( PNTOCC ) THEN

            SRAD  = (1+ALPHA) * SVMXFR
            TDIST = FDIST
         ELSE
            SRAD  = (1+ALPHA) * SVMXBR
            TDIST = BDIST
         END IF

         IF ( SRAD .LT. TDIST ) THEN

            MAXANG = DASINE ( SRAD / TDIST,  ATOL )

            IF ( FAILED() ) THEN
               CALL CHKOUT ( 'ZZGFOCST' )
               RETURN
            END IF

            IF ( TRGSEP .GT. MAXANG ) THEN
C
C              No occultation is possible.
C
               OCSTAT = .FALSE.

               CALL CHKOUT ( 'ZZGFOCST' )
               RETURN

            END IF

         END IF

C
C        We'll need the negatives of the observer-target vectors in
C        several places later, so compute them now.
C
         CALL VMINUS ( FRTPOS, FRTOBS )
         CALL VMINUS ( BCKPOS, BCKOBS )

C
C        Now check for an occulted state assuming a spherical extended
C        body with radius equal to the minimum semi-axis. Again,
C        adjust the sphere for our light time error.
C
         IF ( PNTOCC ) THEN
            MINANG = DASINE ( (1-ALPHA)*SVMNFR / FDIST, ATOL )
         ELSE
            MINANG = DASINE ( (1-ALPHA)*SVMNBR / BDIST, ATOL )
         END IF

         IF ( FAILED() ) THEN
            CALL CHKOUT ( 'ZZGFOCST' )
            RETURN
         END IF

         IF ( TRGSEP .LT. MINANG ) THEN
C
C           The targets must overlap as seen from the observer.
C
            IF ( PNTOCC ) THEN
C
C              Examine the angle between the vector from FRONT to the
C              observer and the vector from FRONT to BACK.  If that
C              angle is greater than or equal to the complement of the
C              angular radius of FRONT, then FRONT occults BACK. First
C              find the position of FRONT and BACK relative to each
C              other.
C
               CALL VMINUS  ( BCKFRT, FRTBCK )

               T2SEP = VSEP ( FRTOBS, FRTBCK )

               IF (  T2SEP  .GT.  ( HALFPI() - MINANG )  ) THEN
C
C                 There must be an occultation.
C
                  OCSTAT = .TRUE.

               ELSE
C
C                 There can't be an occultation: the "back" object
C                 is actually in transit across the "front" object.
C
                  OCSTAT = .FALSE.

               END IF

            ELSE
C
C              We're looking for a point transit condition.
C
               T2SEP = VSEP ( BCKOBS, BCKFRT )

               IF (  T2SEP  .LT.  ( HALFPI() - MINANG )  ) THEN
C
C                 There must be a transit.
C
                  OCSTAT = .TRUE.

               ELSE
C
C                 There can't be a transit: the "back" object
C                 actually occults the "front" object.
C
                  OCSTAT = .FALSE.

               END IF

            END IF

C
C           OCSTAT has been set.
C
            CALL CHKOUT ( 'ZZGFOCST' )
            RETURN

         END IF

C
C        If we've reached this point, we have a situation where we
C        can't classify the geometry using bounding spheres. Instead,
C        we'll see whether the observer-point target vector intersects
C        the extended body.
C
         IF ( PNTOCC ) THEN
C
C           The front body is the extended one.
C
            NCALLS = NCALLS + 1

            CALL SINCPT ( SVFMTH, SVFNAM,  TIME,    SVFFRM,
     .                    SVCORR, SVONAM,  'J2000', BCKPOS,
     .                    SPOINT, TRGEPC,  SRFVEC,  FOUND )

            IF ( FAILED()) THEN
               CALL CHKOUT ( 'ZZGFOCST' )
               RETURN
            END IF

            IF ( FOUND ) THEN
C
C              There's an intercept. If the distance from the observer
C              to the intercept is less than the distance from the
C              observer to the back target, then the back target is
C              occulted; otherwise there's a point transit, which is
C              not considered an occultation in this case.
C
               OCSTAT  =  VNORM(SRFVEC) .LT. BDIST

            ELSE
C
C              There's no overlap and hence no occultation.
C
               OCSTAT = .FALSE.

            END IF

         ELSE
C
C           The back body is the extended one.
C
            CALL SINCPT ( SVBMTH, SVBNAM,  TIME,    SVBFRM,
     .                    SVCORR, SVONAM,  'J2000', FRTPOS,
     .                    SPOINT, TRGEPC,  SRFVEC,  FOUND )

            IF ( FAILED()) THEN
               CALL CHKOUT ( 'ZZGFOCST' )
               RETURN
            END IF

            IF ( FOUND ) THEN
C
C              There's an intercept. If the distance from the observer
C              to the intercept is greater than the distance from the
C              observer to the front target, then the front target is
C              in transit across the back target; otherwise there's a
C              point occultation, which is not considered a transit in
C              this case.
C
               OCSTAT  =  VNORM(SRFVEC) .GT. FDIST

            ELSE
C
C              There's no overlap and hence no occultation.
C
               OCSTAT = .FALSE.

            END IF

         END IF

      ELSE
C
C        Bad combination of shapes. We expect this situation to have
C        been caught at initialization time, but make this check for
C        safety.
C
         CALL SETMSG ( 'The combination of shapes of front and '
     .   //            'back targets is not supported: front '
     .   //            'shape = #; back shape = #.'             )
         CALL ERRCH  ( '#', SVFSHP                              )
         CALL ERRCH  ( '#', SVBSHP                              )
         CALL SIGERR ( 'SPICE(INVALIDSHAPECOMBO)'               )
         CALL CHKOUT ( 'ZZGFOCST'                               )
         RETURN

      END IF

      CALL CHKOUT ( 'ZZGFOCST' )
      RETURN
      END


text:
Address(North El Macero Drive) is a residential street in El Macero.

 Intersecting Streets

  Mace Boulevard and across the intersection El Macero Drive
  Country Club Drive
  Middle Golf Drive
  Middle Golf Drive again
  Sandy Circle
  Country Club Drive again
  Caddy Court and across the intersection East El Macero Drive

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
program parser
use, intrinsic :: ISO_C_BINDING
implicit none

integer (kind = 4) :: physical_names_count;
integer (kind = 4) :: nodes_count;
integer (kind = 4) :: elements_count;

integer (kind = 4) :: i

type node_info
    integer (kind = 4) :: number;
    real    (kind = 4) :: x;
    real    (kind = 4) :: y;
    real    (kind = 4) :: z;
end type node_info
type(node_info), allocatable, dimension(:) :: nodes_info

type element_info
    integer (kind = 4) :: number;
    integer (kind = 4) :: type;
    integer (kind = 4) :: tag_count;
    integer (kind = 4), dimension(10) :: tags
    integer (kind = 4), dimension(8) :: nodes
end type element_info
type(element_info), allocatable, dimension(:) :: elements_info

type physical_name_info
    integer (kind = 4) :: dimension
    integer (kind = 4) :: physical_number
    character (len = 1), dimension(32) :: name
end type physical_name_info
type(physical_name_info), allocatable, dimension(:) :: physical_names_info

type ine_info

    integer (kind = 4) :: indicator;
    integer (kind = 4), dimension(2) :: nodes;
    integer (kind = 4) :: element;
    integer (kind = 4) :: edge_number;
end type ine_info
type(ine_info), allocatable, dimension(:) :: ines_info

print *, "Starting from fortran..."
print *, "Calling C code..."

call parse_input_file(physical_names_count, nodes_count, elements_count)
call label_the_edges()

allocate(elements_info(elements_count))
allocate(nodes_info(nodes_count))
allocate(ines_info(elements_count * 3))

do i = 1, elements_count
    call retrieve_element_info(i, elements_info(i)%number, elements_info(i)%type, &
     elements_info(i)%tag_count, elements_info(i)%tags, elements_info(i)%nodes)
    print *, i, "element ", elements_info(i)%number , "[", elements_info(i)%nodes(1), &
     ", ", elements_info(i)%nodes(2), ", ", elements_info(i)%nodes(3), "]"
end do

do i = 1, nodes_count
    call retrieve_node_info(i, nodes_info(i)%number, nodes_info(i)%x, nodes_info(i)%y, nodes_info(i)%z)
    print *, i, "node ", nodes_info(i)%number , "[", nodes_info(i)%x, ", ", nodes_info(i)%y, ", ", nodes_info(i)%z, "]"
end do

do i = 1, (elements_count * 3)
    call retrieve_ine_array(i, ines_info(i)%indicator, ines_info(i)%nodes, &
    ines_info(i)%element, ines_info(i)%edge_number)
    print *, i, "indicator ", ines_info(i)%indicator , "[", ines_info(i)%nodes(1), &
    ", ", ines_info(i)%nodes(2), ", ", ines_info(i)%element, ", ", ines_info(i)%edge_number,"]"
end do

end program parser

max_stars_repo_name: expee/MSH-File-Parser

max_issues_repo_name: expee/MSH-File-Parser

max_forks_repo_name: expee/MSH-File-Parser


text:
      subroutine ksinit (x,xlb,xub,scale,work,mdv,mcon,mobj,mside,
     1                   mscale,jprnt,jtmax,jgrad,jsdrst,rdfun,adfun,
     2                   fdelt,fdmin,rhomn,rhomx,rhodl,munit,ireq)
      implicit double precision (a-h,o-z)
      common /kscomm/ rdf   ,adf   ,fdl   ,fdm   ,rho   ,drho  ,rhomax,
     1                fun0  ,slope ,delx  ,alpha ,alpmax,a1    ,a2    ,
     2                a3    ,a4    ,f1    ,f2    ,f3    ,f4    ,alim  ,
     3                atest ,ftest ,ifscl ,ifoff ,isx   ,isx0  ,isxlb ,
     4                isxub ,iscl  ,ig0   ,idf   ,islp  ,iobj0 ,iy    ,
     5                ip    ,ih    ,ihess ,iside ,isact ,idobj ,idg   ,
     6                itmp1 ,itmp2 ,inext ,jnext ,jsel  ,itcnt ,icntr ,
     7                icnta ,isdflg,isdrst,ifncl ,nunit ,ndv   ,ncon  ,
     8                nobj  ,nside ,nscale,iprnt ,itmax ,igrad ,limit
      dimension x(*),xlb(*),xub(*),scale(*),work(*)
c
c          routine to initialize optimization parameters
c          for routine ksopt
c
c          input parameters are
c
c               x     -- initial design variable values
c
c               xlb   -- lower bounds on design variables
c
c               xub   -- upper bounds on design variables
c
c               scale -- scale factors for design variables
c
c               work  -- scratch array
c
c               mdv   -- number of design variables
c
c               mcon  -- number of constraints
c
c               mobj  -- number of objective functions
c
c               mside -- flag  =0 if no side constraints
c
c               mscale-- flag selecting design variable scaling
c                        =0 -- no scaling
c                        <0 -- user supplied scaling in vector scale
c                        >0 -- automatic scaling by ksopt every mscale
c                              iterations
c
c               jprnt -- print level flag (3 digit)
c                        each print level includes all lower print levels
c
c                        hundreds digit - one-dimensional search print
c                                         0 = no print
c                                         1 = alpha and k-s function
c                                         2 = proposed d.v. vector
c
c                        tens     digit - gradient print
c                                         0 = no print
c                                         1 = df and dg
c
c                        ones     digit - d.v. and constraint print
c                                         0 = no print
c                                         1 = initial and final iterations
c                                         2 = all iterations
c                                         3 = slope and search direction
c                                         4 = hessian matrix
c
c               jtmax -- maximum number of iterations
c                        (default is 20)
c
c               jgrad -- flag selecting user supplied information
c                        =0 -- user supplies function and constraints
c                              with gradients computed by finite
c                              differences from within ksopt
c                        =1 -- user supplies all function, constraint,
c                              and gradient information
c                        (default is 0)
c
c               jsdrst-- number of iterations before restarting search
c                        direction with steepest descent
c                        (default is ndv+1)
c
c               rdfun -- relative function change for termination
c                        (default is 0.01)
c
c               adfun -- absolute function change for termination
c                        (default is 0.0)
c
c               fdelt -- step size for computing finite differences
c                        derivatives
c                        (default is 0.01)
c
c               fdmin -- minimum difference for computing
c                        finite differences
c                        (default is 0.0001)
c
c               rhomn -- minimum multiplier for ks function
c                        (default is 5.0)
c
c               rhomx -- maximum multiplier for ks function
c                        (default is 100.0)
c
c               rhodl -- increment for rho
c                        (default is computed internally)
c
c               munit -- fortran unit number for output file
c                        (default is 0)
c
c          output parameters are
c
c               work  -- array containing information required
c                        by routine ksopt.  this array must not
c                        be altered during execution of ksopt.
c
c               ireq  -- required length of work array
c
c          author   - Gregory A. Wrenn
c          location - Lockheed Engineering and Sciences Co.
c                     144 Research Drive
c                     Hampton, Va. 23666
c
c          last modification - 19 July 1996
c
c          these subroutines are compatible with fortran 77
c
      jsel   = 1
      ndv    = mdv
      ncon   = mcon
      nobj   = mobj
      nside  = 0
      nscale = mscale
      iprnt  = 1
      itmax  = 20
      igrad  = 0
      isdrst = jsdrst
      rdf    = 0.01
      adf    = 0.0
      fdl    = 0.01
      fdm    = 0.0001
      rho    = 5.0
      rhomax = 100.0
      drho   = 0.0
      nunit  = 0
      itcnt  = 1
      icntr  = 0
      icnta  = 0
      isdflg = 0
      ifncl  = 0
c
      if (mdv   .le. 0)   ndv   = 0
      if (mcon  .le. 0)   ncon  = 0
      if (mobj  .le. 0)   nobj  = 0
      if (mside .ne. 0)   nside = 1
      if (jprnt .ge. 0)   iprnt = jprnt
      if (jtmax .gt. 0)   itmax = jtmax
      if (jgrad .ne. 0)   igrad = 1
      if (jsdrst.le. 0)   isdrst= ndv + 1
      if (rdfun .gt. 0.0) rdf   = rdfun
      if (adfun .gt. 0.0) adf   = adfun
      if (fdelt .gt. 0.0) fdl   = fdelt
      if (fdmin .gt. 0.0) fdm   = fdmin
      if (rhomn .gt. 0.0) rho   = rhomn
      if (rhomx .ge. rho) rhomax= rhomx
      if (rhodl .gt. 0.0) drho  = rhodl
      if (munit .gt. 0)   nunit = munit
c
      ntmp  = 2 * ndv
      if (ntmp .lt. (nobj + ncon)) ntmp = nobj + ncon
      ifscl = 64
      ifoff = ifscl + nobj
      isx   = ifoff + nobj
      isx0  = isx   + ndv
      isxlb = isx0  + ndv
      isxub = isxlb + ndv
      iscl  = isxub + ndv
      ig0   = iscl  + 2 * ndv
      idf   = ig0   + ncon
      islp  = idf   + ndv
      iobj0 = islp  + ndv
      iy    = iobj0 + nobj
      ip    = iy    + ndv
      ih    = ip    + ndv
      ihess = ih    + ndv * (ndv + 1) / 2
      iside = ihess + ndv * (ndv + 1) / 2
      isact = iside + ndv
      idobj = isact + ndv
      idg   = idobj + nobj * ndv
      itmp1 = idg   + ncon * ndv
      itmp2 = itmp1 + ntmp
      ireq  = itmp2 + ntmp - 1
c
      if (ndv .eq. 0) go to 40
      do 30 i = 1,ndv
        j = i - 1
        xx  = x(i)
        if (nside .eq. 0) go to 20
        xbl = xlb(i)
        xbu = xub(i)
        if (xbl .le. xbu) go to 10
        xb = (xbl + xbu) / 2.0
        xbl = xb
        xbu = xb
   10 continue
        if (xx .lt. xbl) xx = xbl
        if (xx .gt. xbu) xx = xbu
        x(i) = xx
   20 continue
        ss = 1.0
        if (nscale .gt. 0) ss = abs(xx)
        if (nscale .lt. 0) ss = abs(scale(i))
        if (ss .lt. 1.0e-04) ss = 1.0
        work(iscl  + j) = ss
        work(isx   + j) = xx / ss
        work(isx0  + j) = xx / ss
        if (nside .eq. 0) go to 30
        work(isxlb + j) = xbl / ss
        work(isxub + j) = xbu / ss
   30 continue
c
      if (drho .gt. 0.0) go to 40
      drho = (rhomax - rho) / 5.0
      if (drho .lt. 10.0) drho = 10.0
      if (drho .gt. 40.0) drho = 40.0
c
   40 continue
c
      if (iprnt .eq. 0 .and. ndv .gt. 0 .and. nobj .gt. 0) go to 60
      write (nunit,70)
      write (nunit,80) ireq
      if (ndv .eq. 0) write (nunit,90)
      if (nobj .eq. 0) write (nunit,100)
      write (nunit,110) ndv,ncon,nobj,nside,nscale,iprnt,itmax,igrad
      write (nunit,120) rdf,adf,fdl,fdm,rho,rhomax,drho
      if (nside .ne. 0) go to 50
      write (nunit,130)
      write (nunit,140) (i,x(i),work(iscl+i-1),i=1,ndv)
      go to 60
   50 continue
      write (nunit,150)
      write (nunit,160) (i,x(i),work(iscl+i-1),
     1                   work(isxlb+i-1)*work(iscl+i-1),
     2                   work(isxub+i-1)*work(iscl+i-1),i=1,ndv)
   60 continue
      if (ndv .eq. 0 .or. nobj .eq. 0) jsel = 0
      call kscomp (work(1))
      return
c
   70 format (1h1,
     1        20x,40h========================================/
     1        21x,40h=                                      =/
     2        21x,40h=    KSOPT Multiobjective Optimizer    =/
     3        21x,40h=                                      =/
     4        21x,40h=    Version 2.8         19 July 96    =/
     5        21x,40h=                                      =/
     6        21x,40h=    Written by    Gregory A. Wrenn    =/
     7        21x,40h=                                      =/
     8        21x,40h========================================///)
   80 format (10x,43hThe work array must be dimensioned at least,
     1        i10,12h words long.//)
   90 format (/53h   Number of design variables is zero -- optimization,
     1        16h is not possible//)
  100 format (/44h   Number of objective functions is zero -- ,
     1        28hoptimization is not possible//)
  110 format (/10h      ndv=,i5,10h     ncon=,i5,10h     nobj=,i5,
     1         10h    nside=,i5,//10h   nscale=,i5,10h   iprint=,i5,
     2         10h    itmax=,i5,10h    igrad=,i5)
  120 format (/10h    rdfun=,e14.7,2x,10h    adfun=,e14.7/
     1        /10h    fdelt=,e14.7,2x,10h    fdmin=,e14.7/
     2        /10h   rhomin=,e14.7,2x,10h   rhomax=,e14.7,2x,
     3         10h     drho=,e14.7)
  130 format (/42(1h-)/
     1        5x,6h d.v. ,6x,7hinitial,9x,6h scale/
     2        5x,6hnumber,6x,7h value ,9x,6hfactor/
     3        1x,42(1h-))
  140 format (4x,i5,4x,e12.5,4x,e12.5)
  150 format (/73(1h-)/
     1        5x,6h d.v. ,6x,7hinitial,9x,6h scale,10x,5hlower,
     2        11x,5hupper/
     3        5x,6hnumber,6x,7h value ,9x,6hfactor,10x,5hbound,
     4        11x,5hbound/
     5        1x,73(1h-))
  160 format (4x,i5,4x,e12.5,4x,e12.5,4x,e12.5,4x,e12.5)
      end

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


      SUBROUTINE PAIR_FORCE(RSQ,FORCEIJ)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Calculates the force between a pair of atoms
!! Currently supports 6 different forcelaws
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      USE header_file
      IMPLICIT NONE

      REAL(kind = double) :: FORCEIJ,RSQ
      REAL(kind = double) :: RRR,RR6,RR12
      REAL(kind = double) :: R,RDIF,RDIF2,RDIF3
      REAL(kind = double) :: RATIO
      REAL(kind = double) :: FACTOR

      SELECT CASE (IPOT_TYPE) ! choose option depending on potential type

        CASE(1)   ! Cut LJ potential

          RRR = 1.0D0/RSQ
          RR6 = RRR*RRR*RRR
          RR12 = RR6*RR6
          FORCEIJ = 48.0d0*EPS*(RR12 - 0.5d0*RR6)*RRR

        CASE(2)   ! Splined-LJ potential

          IF(RSQ.LE.RCUT_INNER_SQ) THEN
            RRR = 1.0D0/RSQ
            RR6 = RRR*RRR*RRR
            RR12 = RR6*RR6
            FORCEIJ = 48.0d0*EPS*(RR12 - 0.5d0*RR6)*RRR
          ELSE
            R = SQRT(RSQ)
            RDIF = R - RCUT
            RDIF2 = RDIF*RDIF
            RDIF3 = RDIF2*RDIF
            FORCEIJ = -(2.0d0*SPLINE_R2*RDIF + 3.0d0*SPLINE_R3*RDIF2)/R
          END IF

        CASE(3)   ! Smooth repulsive potential

          RATIO = RSQ/CUTSQ
          FORCEIJ = EPSIG8*(1.0D0 - RATIO)**3

        CASE(4)   ! m-n family of potentials

          FACTOR = 2.0D0 - RSQ
          FORCEIJ = FORCE_MN_FACT*(FACTOR**N_POWER_SUB1 - FACTOR**M_POWER_SUB1)

        CASE(5)   ! Lucy potential

          R = SQRT(RSQ)
          RATIO = R/RCUT
          FORCEIJ = FORCE_FACT_LUCY*(1.0D0 - RATIO)**2

        CASE(6)   ! Soft disk potential

          FORCEIJ = 2.0d0*SOFT_N/(RSQ**(SOFT_N+1))

        CASE(7)   ! WCA fluid
          RRR = 1.0D0/RSQ
          RR6 = RRR*RRR*RRR
          RR12 = RR6*RR6
          FORCEIJ = 48.0d0*(RR12 - 0.5d0*RR6)*RRR

      END SELECT

      END SUBROUTINE


text:
      FUNCTION fmin (ax, bx, f, tol)
      USE stel_kinds
      IMPLICIT NONE
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      REAL(rprec) ax, bx, f, tol
C-----------------------------------------------
C   L o c a l   P a r a m e t e r s
C-----------------------------------------------
      REAL(rprec),PARAMETER :: zero = 0,one = 1, two = 2,
     1   three = 3, five = 5, p5 = 0.5_dp
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      REAL(rprec):: a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,
     1   u,v,w,fu,fv,fw,fx,x,fmin
C-----------------------------------------------
C   E x t e r n a l   F u n c t i o n s
C-----------------------------------------------
      EXTERNAL f
C-----------------------------------------------
c
c  c is the squared inverse of the golden ratio
c
      c = p5*(three - SQRT(five))
c
c  eps is approximately the square root of the relative machine
c  precision.
c
      eps = one
      eps = eps/two
      tol1 = one + eps
 1003 CONTINUE
      IF (tol1 .le. one) GOTO 1002
      eps = eps/two
      tol1 = one + eps
      IF (tol1 .le. one) GOTO 1002
      eps = eps/two
      tol1 = one + eps
      IF (tol1 .le. one) GOTO 1002
      eps = eps/two
      tol1 = one + eps
      IF (tol1 .le. one) GOTO 1002
      eps = eps/two
      tol1 = one + eps
      GOTO 1003
 1002 CONTINUE
      eps = SQRT(eps)
c
c  initialization
c
      a = ax
      b = bx
      v = a + c*(b - a)
      w = v
      x = v
      e = zero
      fx = f(x)
      fv = fx
      fw = fx
c
c  main loop starts here
c
   20 CONTINUE
      xm = p5*(a + b)
      tol1 = eps*ABS(x) + tol/three
      tol2 = two*tol1
c
c  check STOPping criterion
c
      IF (ABS(x - xm) .le. tol2 - p5*(b - a)) GOTO 90
c
c is golden-section necessary
c
      IF (ABS(e) .le. tol1) GOTO 40
c
c  fit parabola
c
      r = (x - w)*(fx - fv)
      q = (x - v)*(fx - fw)
      p = (x - v)*q - (x - w)*r
      q = two*(q - r)
      IF (q .gt. zero) p = -p
      q = ABS(q)
      r = e
      e = d
c
c  is parabola acceptable
c
      IF (ABS(p) .ge. ABS(p5*q*r)) GOTO 40
      IF (p .le. q*(a - x)) GOTO 40
      IF (p .ge. q*(b - x)) GOTO 40
c
c  a parabolic interpolation step
c
      IF( q.ne.zero )d = p/q
      u = x + d
c
c  f must not be evaluated too CLOSE to ax or bx
c
      IF (u - a < tol2) d = SIGN(tol1,xm - x)
      IF (b - u < tol2) d = SIGN(tol1,xm - x)
      GOTO 50
c
c  a golden-section step
c
   40 CONTINUE
      IF (x .ge. xm) THEN
         e = a - x
      ELSE
         e = b - x
      END IF
      d = c*e
c
c  f must not be evaluated too CLOSE to x
c
   50 CONTINUE
      IF (ABS(d) .ge. tol1) THEN
         u = x + d
      ELSE
         u = x + SIGN(tol1,d)
      END IF
      fu = f(u)
c
c  update  a, b, v, w, and x
c
      IF (fu .gt. fx) GOTO 60
      IF (u .ge. x) THEN
         a = x
      ELSE
         b = x
      END IF
      v = w
      fv = fw
      w = x
      fw = fx
      x = u
      fx = fu
      GOTO 20
   60 CONTINUE
      IF (u .lt. x) THEN
         a = u
      ELSE
         b = u
      END IF
      IF (fu .le. fw) GOTO 70
      IF (w .eq. x) GOTO 70
      IF (fu .le. fv) GOTO 80
      IF (v .eq. x) GOTO 80
      IF (v .eq. w) GOTO 80
      GOTO 20
   70 CONTINUE
      v = w
      fv = fw
      w = u
      fw = fu
      GOTO 20
   80 CONTINUE
      v = u
      fv = fu
      GOTO 20
c
c  END of main loop
c
   90 CONTINUE
      fmin = x

      END FUNCTION fmin


! { dg-do compile }
!
! PR fortran/40873
!
! Failed to compile (segfault) with -fwhole-file.
! Cf. PR 40873 comment 24; test case taken from
! PR fortran/31867 comment 6.
!

pure integer function lensum (words, sep)
  character (len=*), intent(in)        :: words(:), sep
  lensum = (size (words)-1) * len (sep) + sum (len_trim (words))
end function

module util_mod
  implicit none
  interface
    pure integer function lensum (words, sep)
      character (len=*), intent(in)        :: words(:), sep
    end function
  end interface
  contains
  function join (words, sep) result(str)
! trim and concatenate a vector of character variables,
! inserting sep between them
    character (len=*), intent(in)        :: words(:), sep
    character (len=lensum (words, sep))  :: str
    integer                              :: i, nw
    nw  = size (words)
    str = ""
    if (nw < 1) then
      return
    else
      str = words(1)
    end if
    do i=2,nw
      str = trim (str) // sep // words(i)
    end do
  end function join
end module util_mod
!
program xjoin
  use util_mod, only: join
  implicit none
  character (len=5) :: words(2) = (/"two  ","three"/)
  write (*,"(1x,'words = ',a)") "'"//join (words, "&")//"'"
end program xjoin