Sources Esope | Cast3M

brui
C BRUI      SOURCE    CB215821  25/04/23    21:15:02     12247          
      SUBROUTINE BRUI
C-----------------------------------------------------------------------
C Génération d'un bruit blanc obéissant à une loi statistique décrite
C via les arguments transmis. Ce bruit est utilisée pour créer :
C  1) Un LISTREEL lorsque on donne le nombre de valeurs à générer ;
C  2) Un objet EVOLUTION si un LISTREEL de temps est fourni ;
C  3) Un CHAMPOINT si le maillage GEO1 est précisé ;
C  4) Un LISTENTI lorsqu'on tire des variables entiere selon un
C     processus de Poisson.
C-----------------------------------------------------------------------
C
C---------------------------
C Phrase d'appel (GIBIANE) :
C---------------------------
C
C                                         | ENTI2         |
C RES1 = 'BRUI' 'BLAN' MOT1 FLOT1 (FLOT2) | LREEL1 (COUL) | (ENTI3) ;
C                                         | GEO1          |
C ou
C
C RES1 = 'BRUI' 'BLAN' 'POIS' ENTI1 ENTI2 (ENTI3)
C
C------------------------
C Opérandes et résultat :
C------------------------
C
C   1e Syntaxe :
C   ------------
C
C BLAN   : Mot indiquant qu'il s'agit d'un bruit blanc
C MOT1   : Mot indiquant la loi statistique suivi par le bruit :
C              MOT1 = 'GAUS' : Distribution gaussienne,
C              MOT1 = 'UNIF' : Distribution uniforme,
C              MOT1 = 'EXPO' : Distribution exponentielle.
C FLOT1  : Moyenne statistique du bruit à créer. Ne sert à rien
C          si MOT1='EXPO'.
C FLOT2  : Ecart type du bruit à créer. Ne sert à rien si MOT1='POIS'.
C ENTI2  : Nombre de valeurs du LISTREEL à générer.
C LREEL1 : LISTREEL contenant la liste des temps pour l'EVOLUTION.
C COUL   : Mot clef indiquant la couleur associée à l'EVOLUTION RES1.
C GEO1   : Maillage contenant le support géométrique du CHAMPOINT.
C ENTI3  : Entier positif ou nul. Modifie l'initialisation du bruit.
C
C RES1   : LISTREEL, EVOLUTION, CHAMPOINT selon la syntaxe utilisée.
C
C   2e Syntaxe :
C   ------------
C
C BLAN   : Mot indiquant qu'il s'agit d'un bruit blanc
C 'POIS' : Mot-cle que les valeurs suivent une distribution de Poisson.
C ENTI1  : Valeur moyenne de la distribution.
C ENTI2  : Nombre de valeurs du LISTENTI à générer.
C ENTI3  : Entier positif ou nul. Modifie l'initialisation du bruit.
C
C RES1   : Resultat, LISTENTI de valeurs aleatoires.
C
C----------------------
C Variables en COMMON :
C----------------------
C
C NBCOUL : Nombre de couleurs admises par CASTEM (in CCGEOME)
C NCOUL  : Tableau de CHAR*4, dim NBCOUL, Noms des couleurs (in CCGEOME)
C IDCOUL : Valeur de la couleurs par défaut (in CCOPTIO)
C IERR   : Numéro de l'erreur détectée (in CCOPTIO)
C IFOUR  : Indique le type de calcul (in CCOPTIO)
C NIFOUR : Numéro de l'harmonique de fourier si IFOUR=1 (cf CCOPTIO)
C TITREE : CHAR*72, titre des tracés (cf CCOPTIO)
C
C-----------------------------------------------------------------------
C
C Langage : ESOPE + FORTRAN77
C
C Modifs  : F.DABBENE 06/95 (Extension LISTREEL et CHAMPOINT)
C Modifs  : S.PASCAL  06/06 (Extension distribution de Poisson)
C
C-----------------------------------------------------------------------
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C
-INC CCGEOME
 
-INC PPARAM
-INC CCOPTIO
-INC SMEVOLL
-INC SMLENTI
-INC SMLREEL
-INC SMCHPOI
-INC SMCOORD
-INC SMELEME
C
      CHARACTER*4 MOTYP(4),MOTB(1)
C
      DATA MOTB  /'BLAN'/
      DATA MOTYP /'GAUS','UNIF','EXPO','POIS'/
C
C- Lecture du type de bruit
C
      CALL LIRMOT(MOTB,1,IVAL,1)
      IF (IVAL.EQ.0) RETURN
C
C- Lecture du type de distribution
C
      CALL LIRMOT(MOTYP,4,IVAB,1)
      IF (IVAB.EQ.0) RETURN
C
C- Lecture de la moyenne et de l'écart type
C
      IF (IVAB.EQ.4) THEN
        CALL LIRENT(NMOYE,1,IRET1)
        IF (IERR.NE.0) RETURN
      ELSE
        CALL LIRREE(VMOYE,1,IRET1)
        IF (IERR.NE.0) RETURN
        CALL LIRREE(ECAR,1,IRET1)
        IF (IERR.NE.0) RETURN
        IF (IVAB.EQ.3) THEN
           REAERR(1) = REAL(VMOYE)
           CALL ERREUR(-304)
           VMOYE = 0.D0
        ENDIF
        IF (ECAR.LE.0.D0) THEN
          REAERR(1) = REAL(0.D0)
          REAERR(2) = REAL(ECAR)
          CALL ERREUR(191)
          RETURN
        ENDIF
      ENDIF
C
C- Lecture d'un LISTREEL, d'un MELEME ou d'un ENTIER
C
      CALL LIROBJ('LISTREEL',IPT1,0,IRET1)
      IF (IERR.NE.0) RETURN
      IK1 = 0
      IF (IRET1.EQ.0) THEN
         CALL LIROBJ('MAILLAGE',IPT1,0,IRET1)
         IF (IERR.NE.0) RETURN
         IK1 = 1
      ENDIF
      IF (IRET1.EQ.0) THEN
         CALL LIRENT(IPT1,1,IRET1)
         IF (IERR.NE.0) RETURN
         IK1 = 2
      ENDIF
C
C- Lecture facultative pour l'initialisation du générateur
C
      CALL LIRENT(NSTRT,0,IRET1)
      IF (IRET1.EQ.0) THEN
         NSTRT = 0
      ELSEIF (NSTRT.LT.0) THEN
         INTERR(1) = 0
         INTERR(2) = NSTRT
         CALL ERREUR(190)
         RETURN
      ENDIF
C
C- Lecture facultative de la couleur si RES1 est une évolution
C
      IF (IK1.EQ.0) THEN
         CALL LIRMOT(NCOUL,NBCOUL,ICOUL,0)
         IF (ICOUL.EQ.0) ICOUL=IDCOUL+1
         ICOUL=ICOUL-1
      ENDIF
C
C------------------------------------------
C Génération du LISTREEL de NPTBLO valeurs
C------------------------------------------
C
C- Initialisation de NPTBLO, nombre de valeurs à générer.
C- Les éléments du MELEME sont transformés en POI1 si nécessaire.
C- Le maillage de pointeur IPT1 est ACTIF en sortie de CHANGE.
C
      IF (IK1.EQ.0) THEN
         MLREEL = IPT1
         SEGACT MLREEL
         NPTBLO = PROG(/1)
         SEGDES MLREEL
      ELSEIF (IK1.EQ.1) THEN
         MELEME = IPT1
         SEGACT MELEME
         NBSOUS = LISOUS(/1)
         IF ((NBSOUS.NE.0).OR.(ITYPEL.NE.1)) THEN
            CALL CHANGE(IPT1,1)
            IF (IERR.NE.0) RETURN
            MELEME = IPT1
         ENDIF
         NPTBLO = NUM(/2)
         SEGDES MELEME
      ELSE
         NPTBLO = IPT1
      ENDIF
C
C- Création du LISTREEL ou du LISTENTI qui va contenir les valeurs
C- créées.
C
      JG = NPTBLO
      IF (JG .GE. 0) THEN
        IF (IVAB.EQ.4) THEN
          SEGINI MLENTI
        ELSE
          SEGINI MLREEL
        ENDIF
      ELSE
C       Cas ou la taille donnee est negative
        INTERR = NPTBLO
        CALL ERREUR(36)
        RETURN
      ENDIF
C
C- Initialisation du générateur TDRAND
C
      DO 10 I=1,NSTRT
         CALL TDRAND(XRAN)
  10  CONTINUE
C
C- Génération du bruit selon le type de loi repéré par IVAB
C- 1 - Distribution Gaussienne
C- 2 - Distribution Uniforme
C- 3 - Distribution Exponentielle
C- 4 - Distribution de Poisson
C
      IF (IVAB.EQ.1) THEN
         DO 20 I=1,NPTBLO
            AK = ECAR
            CALL TDRAND(XRAN)
            IF (XRAN.GT.0.5D0) THEN
               AK   = -ECAR
               XRAN = 1.D0 - XRAN
            ENDIF
            IF (XRAN.LT.1.D-6) XRAN=1.D-6
            T  = SQRT( LOG(1.D0 / (XRAN*XRAN)) )
            YY = VMOYE + AK * ( T - (2.30753D0 + 0.27061D0*T) /
     #                   (1.0D0 + T * (0.99229D0 + 0.04481D0*T)))
            PROG(I) = YY
  20     CONTINUE
      ELSEIF (IVAB.EQ.2) THEN
         DO 30 I=1,NPTBLO
            CALL TDRAND(XRAN)
            YY =VMOYE + (XRAN - 0.5D0) * 2.D0 * ECAR
            PROG(I)=YY
  30     CONTINUE
      ELSEIF (IVAB.EQ.3) THEN
         DO 40 I=1,NPTBLO
            CALL TDRAND(XRAN)
            IF (XRAN.LT.1.D-6) XRAN = 1.D-6
            YY = -LOG(XRAN) * ECAR
            PROG(I)= YY
  40     CONTINUE
      ELSE
C Pour generer des variables selon une distrib. de Poisson, on emploi 2
C methodes differentes selon que la moyenne de la distrib. est sup. ou
C non a la valeur 50 :
C  - Si sup. a 50 : approximation par une gaussienne ;
C  - Sinon : methode directe.
        XMOYE=FLOAT(NMOYE)
        IF (NMOYE.GE.50) THEN
          DO 50 I=1,NPTBLO
            AK = SQRT(XMOYE)
            CALL TDRAND(XRAN)
            IF (XRAN.GT.0.5D0) THEN
               AK   = -1.D0*AK
               XRAN = 1.D0 - XRAN
            ENDIF
            IF (XRAN.LT.1.D-6) XRAN=1.D-6
            T  = SQRT( LOG(1.D0 / (XRAN*XRAN)) )
            YY = XMOYE + AK * ( T - (2.30753D0 + 0.27061D0*T) /
     #                   (1.0D0 + T * (0.99229D0 + 0.04481D0*T)))
            LECT(I) = INT(YY)
  50      CONTINUE
        ELSE
          XMOYE = EXP(-1.D0*XMOYE)
          DO 60 I=1,NPTBLO
            CALL TDRAND(XRAN)
            A=XRAN
            DO 61 J=1,(10*NPTBLO)
              CALL TDRAND(XRAN)
              A=A*XRAN
              IF (A.LT.XMOYE) THEN
                LECT(I) = J
                GOTO 60
              ENDIF
  61        CONTINUE
  60      CONTINUE
        ENDIF
      ENDIF
C
      IF (IVAB.EQ.4) THEN
        SEGDES MLENTI
      ELSE
        SEGDES MLREEL
      ENDIF
C
C- Création des objets EVOLUTION, CHAMPOINT, LISTREEL ou LISTENTI
C
      IF (IK1.EQ.0) THEN
         N = 1
         SEGINI MEVOLL
         IEVTEX   = TITREE
         ITYEVO   = 'REEL'
         SEGINI KEVOLL
         KEVTEX    = TITREE
         IEVOLL(1) = KEVOLL
         NUMEVX    = ICOUL
         NUMEVY    = 'REEL'
         TYPX      = 'LISTREEL'
         TYPY      = 'LISTREEL'
         IPROGX = IPT1
         NOMEVX = 'TEMPS'
         IPROGY = MLREEL
         NOMEVY = 'SIGNAL'
         SEGDES KEVOLL,MEVOLL
         CALL ECROBJ('EVOLUTIO',MEVOLL)
      ELSEIF (IK1.EQ.1) THEN
         NAT    = 1
         NSOUPO = 1
         SEGINI MCHPOI
         MTYPOI    = ' '
         MOCHDE    = ' '
         JATTRI(1) = 2
         IFOPOI    = IFOUR
         NC        = 1
         SEGINI MSOUPO
         IPCHP(1)  = MSOUPO
         NOCOMP(1) = 'SCAL'
         IGEOC     = MELEME
         NOHARM(1) = NIFOUR
         N = NPTBLO
         SEGINI MPOVAL
         IPOVAL = MPOVAL
         SEGACT MLREEL
         DO 70 I=1,NPTBLO
            VPOCHA(I,1) = PROG(I)
  70     CONTINUE
         SEGSUP MLREEL
         CALL ACTOBJ('CHPOINT ',MCHPOI,1)
         CALL ECROBJ('CHPOINT ',MCHPOI)
      ELSE
        IF (IVAB.EQ.4) THEN
          CALL ECROBJ('LISTENTI',MLENTI)
        ELSE
          CALL ECROBJ('LISTREEL',MLREEL)
        ENDIF
      ENDIF
C
      END