Télécharger histog.eso

Retour à la liste

histog
C HISTOG    SOURCE    CB215821  24/04/12    21:16:14     11897          
      SUBROUTINE HISTOG
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
C----------------------------------------------------------------------C
C                                                                      C
C SUBROUTINE DE L'OPERATEUR "HIST" :                                   C
C                                                                      C
C----------------------------------------------------------------------C
C
C SYNTAXE 1 :
C
C      LENT1 LENT2  = HIST  LVAL1 LCLAS1 (LVAL2 LCLAS2 ( ... )) ...
C                           ... ('CLAS' 'OCCU') ;
C OBJET :
C
C      Etant donne :
C      - les N evenements definis par les m-uplets {X1 X2 ... Xm}
C        fournis sous la forme de m LISTREELS : LVAL1, ... LVALm
C      - les classes correspondantes LCLAS1 (de dime N1+1), ...
C        LCLASm (de dime Nm+1) de type LISTREEL egalement,
C      L'operateur 'HIST' renvoie le LISTENTI correspondant :
C      - a la classe de chaque evenement (option 'CLAS') --> LENT1
C      - au nombre d'occurences des evenements dans chacune des
C        classes (option 'OCCU')  --> LENT2
C      Par defaut (aucune option), on renvoie les 2 LISTENTI.
C
C REMARQUE :
C
C      On numerote de maniere globale les classes de telle sorte que
C      la k^eme classe (numero global) renvoie aux classes k1, k2 ...
C      avec : k = k1 + N1*(k2-1) + N1*N2*(k3-1) + ...
C
C EXEMPLE :
C
C      Soit la suite de 4 evenements :
C        {0.2 4} {0.1 14} {0.5 10} {0.4 1}
C      definie par :
C        x1 = prog  0.2   0.1  1.1   0.5  0.4 ;
C        x2 = prog  4.0  14.0  9.0  10.0  1.0 ;
C      et les classes associees :
C        y1 = prog  0.  0.5  1. ;
C        y2 = prog  0.  5.  10.  15.  20.;
C
C      la numerotation globale des classes est :
C            0      0.5     1.0
C          0 +-------+-------+----->y1
C            |   1   |   2   |
C          5 +-------+-------+
C            |   3   |   4   |
C         10 +-------+-------+
C            |   5   |   6   |
C         15 +-------+-------+
C            |   7   |   8   |
C         20 +-------+-------+
C            |
C         y2 v
C
C      lclass loccu = HIST x1 y1 x2 y2 'CLAS' 'OCCU';
C      -->  lclass contient la suite de 5 entiers :
C             1  5  6  1
C      -->  loccu contient la suite de 8 entiers :
C             2  0  0  0  1  1  0  0
C
C----------------------------------------------------------------------C
C                                                                      C
C SYNTAXE 2 :                                                          C
C                                                                      C
C       EV2 = HIST (COUL) MOD1 CHAM1 (ABS) LRE1 |(MOT1 )| ;            C
C                                               |(LMOT1)|              C
C                                                                      C
C OBJET :                                                              C
C       L'OPERATEUR 'HIST' RENVOIE UNE EVOLUTION DE TYPE HISTOGRAMME   C
C       REPRESENTANT LA DENSITE DE DISTRIBUTION DES VALEURS D'UN       C
C       MCHAML SUR UN MAILLAGE.                                        C
C                                                                      C
C PRINCIPE :                                                           C
C                                                                      C
C     LES INTERVALLES D'ECHANTILLONNAGE DE LA DENSITE DE DISTRIBUTION  C
C     SONT FOURNIS PAR LE 'LISTREEL' LRE1. DANS LE CAS OU CES INTERVAL-C
C     LES NE SONT CONSTANTS, ON REALISE UN HACHAGE EN INTERVALLES REGU-C
C     -LIERS DE TAILLE INFERIEURE AU + PETIT INTERVALLE DE LRE1. POUR  C
C     CHAQUE INTERVALLE, ON MESURE LA LONGUEUR, LA SURFACE OU LE VOLU- C
C     -ME DU MAILLAGE ASSOCIE AUX POINTS SUPPORTS DU CHAMP DONT LA VA- C
C     -LEUR EST DANS L'INTERVALLE. CETTE MESURE EST DONNEE PAR LE PRO- C
C     -DUIT DU POIDS DE GAUSS DU POINT CONSIDERE AVEC LA VALEUR DU JA- C
C     -COBIEN EN CE POINT.                                             C
C                                                                      C
C REMARQUE :                                                           C
C                                                                      C
C     UN CHAMP EXPRIME AU NOEUDS EST INTERPOLE AUX POINTS SUPPORTS DES C
C     RIGIDITES. POUR LES AUTRES POINTS SUPPORTS, C'EST LE JACOBIEN    C
C     QUI EST INTERPOLE AUX POINT SUPPORTS DU CHAMP.                   C
C                                                                      C
C DESCRIPTION DES VARIABLES D'INTERET DU PROGRAMME :                   C
C                                                                      C
C     IPMOD1   : POINTEUR SUR MMODEL                                   C
C     IPCHE1   : POINTEUR SUR CHAM1                                    C
C     MLREE1   : LISTREEL DES INTERVALLES DE MESURE                    C
C     MLREE2   : ORDONNEES DE LA DENSITE DE DISTRIBUTION               C
C     LMOT1    : LISTMOTS CONTENANT LE OU LES NOMS DES COMPOSANTES DE  C
C                CHAM1 A TRAITER                                       C
C     N1       : NOMBRE D'INTERVALLES DE MLREE1                        C
C     X1,X2    : BORNES INF ET SUP DES INTERVALLES D'ECHANTILLONNAGE   C
C     XPASH    : PAS DU HACHAGE                                        C
C     NH       : NOMBRE D'INTERVALLES DU HACHAGE                       C
C     XINDIC   : SEGMENT CONTENANT LES RELATIONS ENTRE LES INTERVALLES C
C                DE MLREE1 ET DU HACHAGE
C     IPCHE2   : POINTEUR SUR MCHAML DU JACOBIEN ASSOCIE AU MAILLAGE.  C
C                OBTENU PAR APPEL A JACOPO                             C
C     TZC1     : TABLEAU D'ENTIER, DE DIMENSION LE NOMBRE DE SOUS-ZONE C
C                DE CHAM1 x NOMBRE DE COMPOSANTES A TRAITER. CONTIENT  C
C                LES POINTEURS SUR LES LISTREELS RELATIFS AUX DISTRI-  C
C                -BUTIONS POUR CHAQUE SOUS-ZONE ET CHAQUE COMPOSANTE   C
C     NCP2     : NOMBRE DE COMPOSANTES EFFECTIVEMENT TRAITEES          C
C     LMOT2    : LISTMOTS, NOMS DES COMPOSANTES EFFECTIVEMENT TRAITEES C
C                CERTAINS NOMS DE COMPOSANTES SPECIFIES DANS LISMOT1   C
C                PEUVENT NE PAS EXISTER DANS CHAM1                     C
C     MLENT1   : LISTE D'ENTIER DONNANT LES CORRESPONDANCES ENTRE LA   C
C                Ie COMPOSANTE DU CHAMP ET LA Je COMPOSANTE A TRAITER  C
C     MLNORM   : LISTREEL : MESURE DE LA SOUS-ZONE ASSOCIEE A LA Je    C
C                COMPOSANTE A TRAITER                                  C
C                                                                      C
C----------------------------------------------------------------------C
 
-INC PPARAM
-INC CCOPTIO
-INC CCGEOME
-INC SMELEME
-INC SMCHAML
-INC SMMODEL
-INC SMLMOTS
-INC SMLENTI
-INC SMLREEL
-INC SMEVOLL
-INC SMINTE
C                                                                      C
C---------------------------- DECLARATIONS ----------------------------C
C                                                                      C
      CHARACTER*4 MOTCLE(1)
      CHARACTER*(LOCOMP) MCOMPJ,MCOMP
      MOTCLE(1)='ABS '
C                                                                      C
      SEGMENT XINDIC
        INTEGER IIND1(IND1,2)
        REAL*8  XIND1(IND1)
      ENDSEGMENT
      SEGMENT TZC
        INTEGER ITZC1(NZ,NC)
      ENDSEGMENT
C
      POINTEUR MLNORM.MLREEL
 
C=======================================================================
C     SYNTAXE 2 PAR DEFAUT
C=======================================================================
C                                                                      C
C----------------------- ACQUISITION DES ENTREES ----------------------C
C                                                                      C
C LECTURE OPTIONNELLE DE LA COULEUR :
      CALL LIRMOT(NCOUL,NBCOUL,ICOUL,0)
      IF (ICOUL.EQ.0) ICOUL=IDCOUL+1
      ICOUL=ICOUL-1
C                                                                      C
C LECTURE DU MODELE (OBLIGATOIRE pour la syntaxe 2):
      MOTERR(1:8)=' MODELE '
      CALL MESLIR(-137)
      CALL LIROBJ('MMODEL  ',IPMOD1,0,IRET)
      IF(IERR.NE.0) RETURN
      IF(IRET.EQ.0) GOTO 100
c     ==> PAS DE MODELE, ON TENTE LA SYNTAXE 1
 
      CALL ACTOBJ('MMODEL  ',IPMOD1,1)
C                                                                      C
C LECTURE OBLIGATOIRE DU CHAMP :
      MOTERR(1:8)=' MCHAML '
      CALL MESLIR(-137)
      CALL LIROBJ('MCHAML',IPCHE1,1,IRET)
      CALL ACTOBJ('MCHAML',IPCHE1,1)
      IF (IERR.NE.0) RETURN
C     REDUCTION IPCHE1 SUR IPMOD1
      CALL REDUAF(IPCHE1,IPMOD1,IPCHE2,0,IRET,KERRE)
      IF (IRET.NE.1) THEN
        CALL ERREUR(KERRE)
        RETURN
      ENDIF
      IPCHE1=IPCHE2
C     VERIFICATION QUE LES SOUS-ZONES ONT TOUTES MEME SUPPORT
      CALL QUESUP(IPMOD1,IPCHE1,0,0,IRET,ISUP1)
      IF (ISUP1.EQ.9999) THEN
        CALL ERREUR(560)
        RETURN
      ENDIF
C                                                                      C
C LECTURE OPTIONNELLE DU MOT-CLE "ABS" :
      CALL LIRMOT(MOTCLE,1,IABS,0)
C                                                                      C
C LECTURE OBLIGATOIRE DU LISTREEL (MLREE1 = INTERVALLES DE MESURE) :
      MOTERR(1:8)='LISTREEL'
      CALL MESLIR(-137)
      CALL LIROBJ('LISTREEL',MLREE1,1,IRET)
      IF (IERR.NE.0) RETURN
      SEGACT,MLREE1
      N1=MLREE1.PROG(/1)
      IF (N1.LE.1) THEN
        CALL ERREUR(249)
        RETURN
      ENDIF
      X1=MLREE1.PROG(1)
      X2=MLREE1.PROG(2)
      IF (X1.GE.X2) THEN
        CALL ERREUR(249)
        RETURN
      ENDIF
C     CONSTRUCTION D'UN HACHAGE DANS LE CAS OU LONG. INTERV. PAS CONST.:
C     1. RECHERCHE DU + PETIT PAS DANS LISTREEL FOURNI (XPAS)
      JG=N1-1
      SEGINI,MLREE3
      XPAS=X2-X1
      MLREE3.PROG(1)=XPAS
      DO IH=3,N1
        X1=X2
        X2=MLREE1.PROG(IH)
        XPASI=X2-X1
        MLREE3.PROG(IH-1)=XPASI
        IF (XPAS.GT.XPASI) XPAS=XPASI
      ENDDO
      IF (XPAS.LE.0.) THEN
        CALL ERREUR(249)
        RETURN
      ENDIF
C     2. CONSTRUCTION DU HACHAGE (INTERVALLES REGULIERS DE TAILLE
C        INFERIEURE AU + PETIT INTERVALLE DE MLREE1)
C        XINDIC PERMET D'ASSOCIER LES INTERVALLES DE MLREE1 AU HACHAGE :
C     IIND1(I,1)=0 : Ie INTERV. HACHAGE A CHEVAL SUR 2 INTERV. MLREE1
C     IIND1(I,1)=1 : Ie INTERV. HACHAGE STRICT. DANS 1 INTERV. MLREE1
C     IIND1(I,2)=K : Ie INTERV. HACHAGE ASSOCIE A Ke INTERV. MLREE1
C     XIND1(I)= BORNE DES 2 INTERV. MLREE1 DANS Ie INTERV. HACHAGE
C               (CAS OU IIND1(I,1)=0)
      X1=MLREE1.PROG(1)
      NH=INT((X2-X1)/XPAS)+1
      IND1=NH
      SEGINI,XINDIC
      IREF=1
      IIND1(1,1)=1
      IIND1(1,2)=IREF
      XPASH=(X2-X1)/FLOAT(NH)
      XREH=X1+XPASH
      DO IIND=2,NH
        XREF=MLREE1.PROG(IREF+1)
        XREH=XREH+XPASH
        IIND1(IIND,2)=IREF
        XIND1(IIND)=XREF
        IF(XREH.GT.XREF) THEN
          IIND1(IIND,1)=0
          IREF=IREF+1
          IF (IREF.GT.N1) THEN
c           write(6,*) 'IREF > N1 !!'
            RETURN
          ENDIF
        ELSE
          IIND1(IIND,1)=1
        ENDIF
      ENDDO
C                                                                      C
C LECTURE OPTIONNELLE DU(DES) NOM(S) DE LA COMPOSANTE A TRAITER :
      NCP1=0
C     DONNEE D'UN MOT ?
      CALL LIRCHA(MCOMP,0,IMOT)
      IF (IMOT.NE.0) THEN
        NCP1=1
        JGM=NCP1
        JGN=LOCOMP
        SEGINI,MLMOT1
        MLMOT1.MOTS(1)=MCOMP(1:IMOT)
      ELSE
C     SINON, DONNEE D'UN LISTMOTS ?
        CALL LIROBJ('LISTMOTS',MLMOT1,0,IRET)
        IF (IRET.NE.0) THEN
          SEGACT,MLMOT1
          NCP1=MLMOT1.MOTS(/2)
        ELSE
C     SINON, TOUTES LES COMPOSANTES DU CHAMP EN ENTREE :
C     (ATTENTION : IL FAUT PARCOURIR TOUTES LES SOUS-ZONES...)
          IMOT=1
          MCHEL1=IPCHE1
          NZ1=MCHEL1.ICHAML(/1)
          MCHAM1=MCHEL1.ICHAML(1)
C         INITIALISATION MLMOTS SUR 1ERE SOUS-ZONE :
          JGN=LOCOMP
          JGM=MCHAM1.NOMCHE(/2)
          SEGINI,MLMOT1
          NCP1=JGM
          DO I=1,NCP1
            MLMOT1.MOTS(I)=MCHAM1.NOMCHE(I)
          ENDDO
C         VERIF. NOMS COMPOSANTES SUR LES AUTRES SOUS-ZONES
          DO I=2,NZ1
            MCHAM1=MCHEL1.ICHAML(I)
            DO J=1,MCHAM1.NOMCHE(/2)
              MCOMPJ=MCHAM1.NOMCHE(J)
              CALL PLACE(MLMOT1.MOTS,NCP1,IPLA,MCOMPJ)
              IF (IPLA.EQ.0) THEN
                JGM=NCP1+1
                SEGADJ,MLMOT1
                NCP1=NCP1+1
                MLMOT1.MOTS(NCP1)=MCOMPJ
              ENDIF
            ENDDO
          ENDDO
        ENDIF
      ENDIF
C                                                                      C
C----------------------------- TRAITEMENT -----------------------------C
C                                                                      C
C JACOBIEN ASSOCIE AU MODELE :
C     CALCUL DU JACOBIEN :
      CALL JACOPO(IPMOD1,IPCHE2,IRET)
      IF (IRET.EQ.0) RETURN
      MMODEL=IPMOD1
      IMODEL=KMODEL(1)
C
C AJUSTEMENT DES POINTS SUPPORTS :
C ISUP1 = VALEUR INDICATEUR POINTS SUPPORTS IPCHE1
C ISUP2 = VALEUR INDICATEUR POINTS SUPPORTS IPCHE2 (JACOBIEN)
      IF (ISUP1.EQ.1) THEN
        CALL CHASUP(IPMOD1,IPCHE1,IPCHE3,IRET,5)
        IF (IRET.NE.0) THEN
          CALL ERREUR(IRET)
          RETURN
        ENDIF
        IPCHE1=IPCHE3
      ELSEIF (ISUP1.NE.5) THEN
        CALL CHASUP(IPMOD1,IPCHE2,IPCHE3,IRET,ISUP1)
        IF (IRET.NE.0) THEN
          CALL ERREUR(IRET)
          RETURN
        ENDIF
        IPCHE2=IPCHE3
      ENDIF
C
C ECHANTILLONNAGE :
      NCP2=0
      MCHEL1=IPCHE1
      MCHEL2=IPCHE2
      NZ1=MCHEL1.ICHAML(/1)
      NZ=NZ1
      NC=NCP1
      SEGINI,TZC
      JGN=LOCOMP
      JGM=NCP1
      SEGINI,MLMOT2
      JG=NCP1
      SEGINI,MLNORM
      SEGINI,MLENT1
      DO I=1,NZ1
        MCHAM1=MCHEL1.ICHAML(I)
        IPT1=MCHEL1.IMACHE(I)
        ICOMP=MCHAM1.NOMCHE(/2)
        MINTE2=MCHEL2.INFCHE(I,4)
        MCHAM2=MCHEL2.ICHAML(I)
        MELVA2=MCHAM2.IELVAL(1)
        DO J=1,ICOMP
          MCOMPJ=MCHAM1.NOMCHE(J)
          CALL PLACE(MLMOT1.MOTS,NCP1,IPLA,MCOMPJ)
          IF (IPLA.NE.0) THEN
            IF (MCHAM1.TYPCHE(J).NE.'REAL*8') THEN
               MOTERR(1:4)='HIST'
               MOTERR(5:8)=MCOMPJ(1:4)
               CALL ERREUR(335)
               RETURN
            ENDIF
            JREF=IPLA
            IPLA=0
            IF (NCP2.NE.0) CALL PLACE(MLMOT2.MOTS,NCP2,IPLA,MCOMPJ)
            IF (IPLA.EQ.0) THEN
              NCP2=NCP2+1
              MLMOT2.MOTS(NCP2)=MCOMPJ
              MLENT1.LECT(JREF)=NCP2
            ENDIF
            JG=N1
            SEGINI,MLREE2
            ITZC1(I,JREF)=MLREE2
            MELVA1=MCHAM1.IELVAL(J)
            NP1=MELVA1.VELCHE(/1)
            NE1=MELVA1.VELCHE(/2)
            NP2=MELVA2.VELCHE(/1)
            NE2=MELVA2.VELCHE(/2)
            DO K=1,IPT1.NUM(/2)
              DO L=1,MINTE2.POIGAU(/1)
                VIJKL2=MINTE2.POIGAU(L)*
     &                 (ABS(MELVA2.VELCHE(MIN(NP2,L),MIN(NE2,K))))
                XIJKL1=MELVA1.VELCHE(MIN(NP1,L),MIN(NE1,K))
                IF (IABS.EQ.1) XIJKL1=ABS(XIJKL1)
                NIJKL1=INT((XIJKL1-X1)/XPASH)+1
                IF (NIJKL1.GE.1.AND.NIJKL1.LE.NH.AND.XIJKL1.GE.X1) THEN
                  IIJKL1=IIND1(NIJKL1,2)
                  IF (IIND1(NIJKL1,1).EQ.1) THEN
                    XIJKL2=VIJKL2/MLREE3.PROG(IIJKL1)
                    MLREE2.PROG(IIJKL1)=MLREE2.PROG(IIJKL1)+XIJKL2
                  ELSE
                    IF (XIJKL1.LT.XIND1(NIJKL1)) THEN
                      XIJKL2=VIJKL2/MLREE3.PROG(IIJKL1)
                      MLREE2.PROG(IIJKL1)=MLREE2.PROG(IIJKL1)+XIJKL2
                    ELSE
                      XIJKL2=VIJKL2/MLREE3.PROG(IIJKL1+1)
                      MLREE2.PROG(IIJKL1+1)=MLREE2.PROG(IIJKL1+1)+XIJKL2
                    ENDIF
                  ENDIF
                ENDIF
                MLNORM.PROG(JREF)=MLNORM.PROG(JREF)+VIJKL2
              ENDDO
            ENDDO
          ENDIF
        ENDDO
      ENDDO
C
C VERIFICATION NOMBRE DE COMPOSANTES EFFECTIVEMENT TRAITEES
c     write(6,*) 'NCP2 =',NCP2
      IF (NCP2.EQ.0) THEN
        IF (IMOT.NE.0) THEN
          MOTERR(1:4)=MLMOT1.MOTS(1)
          MOTERR(5:8)='    '
          CALL ERREUR(77)
          RETURN
        ELSE
          CALL ERREUR(26)
          RETURN
        ENDIF
      ENDIF
      IF (NCP2.LT.NCP1) THEN
        JGN=LOCOMP
        JGM=NCP2
        SEGADJ,MLMOT2
      ENDIF
      SEGSUP,MLREE3
C
C------------------------------- SORTIE -------------------------------C
C                                                                      C
C CONSTRUCTION DE L'EVOLUTION :
C     SOMMATION DES ORDONNEES SUR CHAQUE SOUS-ZONE,
C     COMPOSANTE PAR COMPOSANTE :
      N=NCP2
      SEGINI,MEVOLL
      ITYEVO='REEL'
      IEVTEX=' '
      DO I=1,NCP1
        ICP=MLENT1.LECT(I)
        IF (ICP.NE.0) THEN
          MLREE2=ITZC1(1,I)
C         SOMME SUR LES SOUS-ZONES :
          DO J=2,ITZC1(/1)
            MLREE3=ITZC1(J,I)
            IF (MLREE2.EQ.0) THEN
              MLREE2=MLREE3
            ELSE
              IF (MLREE3.NE.0) THEN
                CALL ADLISR(MLREE2,MLREE3,IRET,1)
                IF (IRET.EQ.0) RETURN
                MLREE2=IRET
              ENDIF
            ENDIF
          ENDDO
C         DIVISION DES ORDONNEES PAR XNORM => DENSITE DE DISTRIB.
          XNORM=MLNORM.PROG(I)
          CALL MUFLIR(MLREE2,XNORM,IRET,-1)
          SEGINI,KEVOLL
          IPROGX=MLREE1
          IPROGY=IRET
          TYPX='LISTREEL'
          TYPY='LISTREEL'
C         COURBE DE TYPE "HIST" POUR TRAITEMENT DANS DESSIN :
          NUMEVX=ICOUL
          NUMEVY='HIST'
          NOMEVX=MLMOT2.MOTS(ICP)
          NOMEVY=' '
          KEVTEX=' '
          MEVOLL.IEVOLL(ICP)=KEVOLL
        ENDIF
      ENDDO
C                                                                      C
C SORTIE :
      CALL ACTOBJ('EVOLUTIO',MEVOLL,1)
      CALL ECROBJ('EVOLUTIO',MEVOLL)
C                                                                      C
C UN PEU DE MENAGE :
      IF (IMOT.NE.0) THEN
        SEGSUP,MLMOT1
      ENDIF
      SEGSUP,MLMOT2,TZC,XINDIC,MLNORM,MLENT1
 
      RETURN
 
 
C=======================================================================
C     SYNTAXE 1 : --> HISTO1
C=======================================================================
 100  CONTINUE
      CALL HISTO1
      END
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales