C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-C
C======================================================================C
C                                                                      C
C     M. Petitjean, J. Comput. Chem. 1994,15[5],507-523.               C
C                                                                      C
C======================================================================C
C                                                                      C
C     DOUBLE PRECISION FUNCTION LOIUAB ( UA , UB )                     C
C     DOUBLE PRECISION                   UA , UB                       C
C     LOIUAB = UA + (UB-UA) * random_real_over_01                      C
C     RETURN                                                           C
C     END                                                              C
C                                                                      C
C======================================================================C
C                                                                      C
C     ****************************************************             C
C     MESURE LA SURFACE DE L'INTERSECTION ET DE LA REUNION             C
C       DE N HYPERSPHERES  PAR LA METHODE DE MONTE-CARLO               C
C     ****************************************************             C
C                                                                      C
C                                                                      C
C     ARGUMENTS D'ENTREE :                                             C
C     ------------------                                               C
C                                                                      C
C     D      : DIMENSION DE L'ESPACE                                   C
C     NOBST  : NOMBRE D'OBSERVATIONS A GENERER                         C
C     N      : NOMBRE DE D-SPHERES                                     C
C     X      : COORDONNEES DES CENTRES : X (N,D)                       C
C     R      : RAYONS DES D-SPHERES                                    C
C                                                                      C
C                                                                      C
C     ARGUMENTS DE SORTIE :                                            C
C     -------------------                                              C
C                                                                      C
C     RR     : PUISSANCES D-1 DES RAYONS DES D-SPHERES                 C
C     SURFT  : SOMME DES SURFACES DES SPHERES                          C
C     INOBSI : NOMBRE D'OBSERVATIONS INTERNES A L'INTERSECTION         C
C     ISUR   : SURFACE DE L'INTERSECTION                               C
C     IDELTS : ECART-TYPE DE LA SURFACE DE L'INTERSECTION              C
C     IRELPS : PRECISION ASSOCIEE A IDELTS ( 0% A 100% )               C
C     RNOBSI : NOMBRE D'OBSERVATIONS INTERNES A LA REUNION             C
C     RSUR   : SURFACE DE LA REUNION                                   C
C     RDELTS : ECART-TYPE DE LA SURFACE DE LA REUNION                  C
C     RRELPS : PRECISION ASSOCIEE A RDELTS ( 0% A 100% )               C
C     OBSD   : BUFFER DE D ELEMENTS REELS                              C
C                                                                      C
C                                                                      C
C     RETOURS MULTIPLES :                                              C
C     -----------------                                                C
C                                                                      C
C     * 1   : AUCUNE OBSERVATION DANS L'INTERSECTION                   C
C     * 2   : AUCUNE OBSERVATION DANS LA REUNION                       C
C                                                                      C
C======================================================================C
C
      SUBROUTINE MSIRNS ( D , NOBST , N , X , R , RR , SURFT ,
     ,                    INOBSI , ISUR , IDELTS , IRELPS ,
     ,                    RNOBSI , RSUR , RDELTS , RRELPS ,
     ,                    OBSD , * , * )
C
      IMPLICIT INTEGER  ( A - Z )
C
      LOGICAL             IOINT , ROINT , INTERN
C
      DOUBLE PRECISION    X (N,D) , R (N) , RR (N) , SURFT ,
     ,                    ISUR , IDELTS , IRELPS ,
     ,                    RSUR , RDELTS , RRELPS ,
     ,                    OBSD (D) ,
     ,                    LOIUAB , DFLOAT , DSQRT ,
     ,                    RRTOT , LR , Z , ZZ , ZN ,
     ,                    IZP , IZQ , ISIGMA ,
     ,                    RZP , RZQ , RSIGMA ,
     ,                    ZERO   , UN     , DEUX   , CENT   , INFINI ,
     ,                    PI , COEFF
C
      DATA                ZERO   , UN     , DEUX   , CENT   , INFINI
     ,                  / 0.D0   , 1.D0   , 2.D0   , 100.D0 , 1.D+30  /
C
      DATA                PI  / 3.141592653589793238462643383279D0 /
C
C
C
C ... COEFFICIENT DESTINE AU CALCUL DE LA SURFACE D'UNE HYPERSPHERE
C     -------------------------------------------------------------
C
      D2     = D / 2
C
      COEFF  = UN
C
      DO I = 1 , D2
          COEFF  = COEFF * PI
      END DO
C
      IF ( D .EQ. (2*D2) ) THEN
          Z      = UN
          DO I = D2 , 2 , -1
              Z      = Z * DFLOAT (I)
          END DO
          COEFF  = COEFF * DFLOAT (D) / Z
      ELSE
          Z      = UN
          DO I = D-2 , 3 , -2
              Z      = Z * DFLOAT (I)
          END DO
          COEFF  = DEUX * COEFF / Z
          DO I = 1 , D2
              COEFF  = COEFF + COEFF
          END DO
      ENDIF
C
C
C
C ... POIDS TOTAL ET SURFACE TOTALE
C     -----------------------------
C
      RRTOT  = ZERO
C
      DO I = 1 , N
          RR (I) = UN
          DO J = 2 , D
              RR (I) = RR (I) * R (I)
          END DO
          RRTOT  = RRTOT + RR (I)
      END DO
C
      SURFT  = COEFF * RRTOT
C
C
C
C ... COMPTAGE DES OBSERVATIONS INTERNES AUX SPHERES
C     ----------------------------------------------
C
      INOBSI  = 0
      RNOBSI  = 0
C
      DO IOBS = 1 , NOBST
C
C
C
C ... ... CHOIX DE LA SPHERE
C         ------------------
C
          Z      = LOIUAB ( ZERO , RRTOT )
          LR     = ZERO
C
          DO I = 1 , N
              LR     = LR + RR (I)
              IF ( Z .LE. LR ) THEN
                  S      = I
                  GOTO 100
              ENDIF
          END DO
          S      = N
C
C
C
C ... ... OBSERVATION ALEATOIRE SUR LA SURFACE DE LA SPHERE S
C         ---------------------------------------------------
C
100       ZZ     = ZERO
          DO J = 1 , D
              OBSD (J) = LOIUAB ( -UN , +UN )
              ZZ       = ZZ + OBSD (J) * OBSD (J)
          END DO
          IF ( ZZ .GT. UN  .OR.  ZZ .LE. ZERO ) GOTO 100
C
          Z      = R (S) / DSQRT (ZZ)
          DO J = 1 , D
              OBSD (J) = X (S,J) + OBSD (J) * Z
          END DO
C
C
C
C ... ... LOCALISATION PAR RAPPORT AUX AUTRES SPHERES
C         -------------------------------------------
C
          IOINT  = .TRUE.
          ROINT  = .FALSE.
C
          DO I = 1 , N
              IF ( I .NE. S ) THEN
                  ZZ     = ZERO
                  DO J = 1 , D
                      Z         = OBSD (J) - X (I,J)
                      ZZ        = ZZ + Z * Z
                  END DO
                  INTERN = ZZ .LE. R(I)*R(I)
                  ROINT  = ROINT  .OR.   INTERN
                  IOINT  = IOINT  .AND.  INTERN
              ENDIF
          END DO
C
          IF ( IOINT ) INOBSI = INOBSI + 1
          IF ( ROINT ) RNOBSI = RNOBSI + 1
C
      END DO
C
      RNOBSI = NOBST - RNOBSI
C
C
C ... CALCUL DES SURFACES ET DES PRECISIONS ASSOCIEES
C     -----------------------------------------------
C
      ZN     = DFLOAT (NOBST)
C
      IZP    = DFLOAT (INOBSI) / ZN
      IZQ    = DFLOAT (NOBST-INOBSI) / ZN
C
      ISUR   = SURFT * IZP
      ISIGMA = DSQRT ( IZP*IZQ / ZN )
      IDELTS = SURFT * ISIGMA
C
      IF ( IZP .NE. ZERO ) THEN
          IRELPS = CENT * ISIGMA  / IZP
      ELSE
          IRELPS = CENT
      ENDIF
C
C
      RZP    = DFLOAT (RNOBSI) / ZN
      RZQ    = DFLOAT (NOBST-RNOBSI) / ZN
C
      RSUR   = SURFT * RZP
      RSIGMA = DSQRT ( RZP*RZQ / ZN )
      RDELTS = SURFT * RSIGMA
C
      IF ( RZP .NE. ZERO ) THEN
          RRELPS = CENT * RSIGMA  / RZP
      ELSE
          RRELPS = CENT
      ENDIF
C
C
      IF ( RNOBSI .LE. 0 ) RETURN 2
      IF ( INOBSI .LE. 0 ) RETURN 1
C
      RETURN
C
      END