Sources Esope | Cast3M

dspr
C DSPR      SOURCE    BP208322  16/11/18    21:16:32     9177           
      SUBROUTINE DSPR
      IMPLICIT INTEGER(I-N)
      IMPLICIT  REAL*8(A-H,O-V)
      IMPLICIT COMPLEX*16 (W-Z)
      CHARACTER *72 TI
      CHARACTER*12 MOTX,MOTY
      CHARACTER*4 MCLE(2)
      LOGICAL INV
C
C=======================================================================
C  =  CALCUL DE LA DENSITE SPECTRALE DE PUISSANCE D'UN SIGNAL         =
C  =                                                                  =
C  =  SYNTAXE :                                                       =
C  =                                                                  =
C  =    SPEC = DSPR EXP2 EVO1 FMIN V1 FMAX V2 (COUL) ;                =
C  =                                                                  =
C  =                                                                  =
C  =    EXP2    : EXPOSANT DANS NPOINT=2^EXP2, NPOINT ETANT           =
C  =              LE NOMBRE DE REELS DANS LISTREEL                    =
C  =    EVO1    : OBJET DE TYPE EVOLUTIO CONTENANT LE SIGNAL A TRAITER=
C  =              ( UNE COURBE SEULEMENT )                            =
C  =    FMIN    : MOT-CLE                                             =
C  =    V1      : FREQUENCE MINIE A VISUALISER                        =
C  =    FMAX    : MOT-CLE                                             =
C  =    V2      : FREQUENCE MAXIE A VISUALISER                        =
C  =    COUL    : COULEUR ATTRIBUEE A L'OBJET CREE (FACULTATIF)       =
C  =                                                                  =
C  =                                                                  =
C  =  CREATION    : 01/04/87                                          =
C  =  PROGRAMMEUR : BEAUFILS                                          =
C  =  MODIFICATIO : PEG WARNING LE IRET2.... 22/2/90                  =
C  =  MODIFICATIO : PEG FREQUENCES CORRECTES  1/6/90                  =
C  =  MODIFICATIO : PEG PERIODOGRAM CORRECT   1/6/90                  =
C=======================================================================
C
-INC CCGEOME
 
-INC PPARAM
-INC CCOPTIO
-INC SMEVOLL
-INC SMLREEL
C
      SEGMENT MTRAV1
          IMPLIED XX(NCOU),YY(NCOU)
      ENDSEGMENT
C
      SEGMENT MTRAV2
          IMPLIED W(NEXP)
      ENDSEGMENT
C
      DATA MOTX/'FREQ      HZ'/
      DATA MOTY/'DSP UN**2/HZ'/
      DATA MCLE(1),MCLE(2)/'FMIN','FMAX'/
C
C     LECTURE DE EXP2
C
      CALL LIRENT(N2,1,IRET1)
      IF(IRET1.EQ.0)GOTO 666
C
C     LECTURE DE L'OBJET EVOLUTIO CONTENANT LE SIGNAL
C
      CALL LIROBJ('EVOLUTIO',IPSIG,1,IRET2)
C PEG IF(IRET1.EQ.0)GOTO 666
      IF(IRET2.EQ.0)GOTO 666
C
C     LECTURE DE LA FREQUENCE MINIE
C
      CALL LIRMOT(MCLE(1),1,IRET,0)
      IF(IRET.EQ.1) THEN
          CALL LIRREE(FRMI,0,IRET1)
          IF(IRET1.EQ.0) THEN
              MOTERR(1:4)=MCLE(1)
              CALL ERREUR(166)
              GOTO 666
          ENDIF
          IF(FRMI.LT.(-1.D-20)) FRMI=-FRMI
          IF(ABS(FRMI).LT.1.D-20) FRMI=-1.D0
      ELSE
          FRMI=-1.D0
      ENDIF
C
C     LECTURE DE LA FREQUENCE MAXIE
C
      CALL LIRMOT(MCLE(2),1,IRET,0)
      IF(IRET.EQ.1) THEN
          CALL LIRREE(FRMA,0,IRET1)
          IF(IRET1.EQ.0) THEN
              MOTERR(1:4)=MCLE(2)
              CALL ERREUR(166)
              GOTO 666
          ENDIF
          IF(FRMA.LT.(-1.D-20)) FRMA=-FRMA
          IF(ABS(FRMA).LT.1.D-20) FRMA=-1.D0
      ELSE
          FRMA=-1.D0
      ENDIF
C
C     LECTURE DE LA COULEUR
C
      CALL LIRMOT(NCOUL,NBCOUL,ICOUL,0)
      IF(ICOUL.EQ.0) ICOUL=IDCOUL+1
      ICOUL=ICOUL-1
C
      IF(IERR.NE.0) GOTO 666
C
      MEVOL1=IPSIG
      SEGACT MEVOL1
      KEVOL1=MEVOL1.IEVOLL(1)
      SEGACT KEVOL1
      MLREE1=KEVOL1.IPROGX
      MLREE2=KEVOL1.IPROGY
      SEGACT MLREE1,MLREE2
C
      L1=MLREE1.PROG(/1)
      DT=MLREE1.PROG(2)-MLREE1.PROG(1)
      SEGDES MLREE1
C
      NPOINT=2**N2
      IF(NPOINT.GT.L1) THEN
          IF(IIMPI.EQ.1) WRITE(IOIMP,1000) L1,N2,NPOINT
 1000     FORMAT(1H ,'LE NOMBRE DE POINTS DANS LISTEMPS ',I6, ' EST ',
     &    'INFERIEURE @ 2**',I5,
     &    /' ON PRENDRA UNE LONGUEUR DE LISTEMPS DE ',I6,' MOTS     ',
     &    /' ET ON COMPLETERA PAR DES ZEROS')
      ELSE
           IF(NPOINT.LT.L1) THEN
               IF(IIMPI.EQ.1) WRITE(IOIMP,1010) N2
 1010          FORMAT(1H ,'ON TRONQUE LE SIGNAL A 2**',I5,' MOTS',/)
           ELSE
               IF(IIMPI.EQ.1) WRITE(IOIMP,1030)N2,NPOINT
 1030          FORMAT(1H ,'LA LONGUEUR DU LISTEMP EST EGALE A 2**',I5,
     &        ' = ',I6,/)
           ENDIF
      ENDIF
C
      NCOU=NPOINT
      SEGINI MTRAV1
C
      NEXP=NPOINT/2
      SEGINI MTRAV2
C
C     CHARGEMENT DES TABLEAUX DE TRAVAIL
C
      IND1=L1
      IF(NPOINT.LT.L1)IND1 = NPOINT
        DO 10 I=1,IND1
        XX(I)=MLREE2.PROG(I)
C       IF(IIMPI.EQ.1)WRITE(IOIMP,*) ' XX(',I,') = ',XX(I)
  10    CONTINUE
      IF(NPOINT.GT.L1) THEN
          L2=L1+1
            DO 11 I=L2,NPOINT
            XX(I)=0.D0
  11        CONTINUE
      ENDIF
C
C     CALCUL DE LA FFT
C
      IF(IIMPI.EQ.1) WRITE(IOIMP,*)' CALCUL DE LA FFT DU SIGNAL '
      INV=.FALSE.
      CALL WEXP(INV,NPOINT,W(1))
      CALL FFTL(XX(1),YY(1),W(1),NPOINT)
C
      IF(IIMPI.EQ.1) WRITE(IOIMP,*)' FFT CALCULEE '
      SEGDES MLREE2
C
C     CREATION ET CALCUL DES LISTES DE LA DSP
C
      IF(IIMPI.EQ.1) WRITE(IOIMP,*)' CALCUL DE LA DSP DU SIGNAL '
      DUREE=DT*DBLE(NPOINT)
      DFREQ=1.D0/DUREE
      KHALF=NEXP+1
CPP +
      KHALF1=KHALF
CPP +
      KMIN=INT(FRMI/DFREQ)+1
      KMAX=INT(FRMA/DFREQ)+1
      KDEBU=1
      IF(KMIN.GT.0) KDEBU=KMIN
      IF((KMAX.LT.KHALF).AND.(FRMA.GT.0.)) KHALF=KMAX
      JG=KHALF-KDEBU+1
      SEGINI MLREE1
      SEGINI MLREE2
CPP ! COEF=2.D0*DFREQ*DT*DT
      COEF=     DFREQ*DT*DT
        DO 20 I=KDEBU,KHALF
        FREQ=DBLE(I-1)*DFREQ
C       IF(IIMPI.EQ.1) WRITE(IOIMP,*)I,FREQ,' XX = ',XX(I)
        MLREE1.PROG(I-KDEBU+1)=FREQ
        YCOMP=CONJG(XX(I))
        PMODU2=XX(I)*YCOMP
cbp: on pourrait aussi ecrire   PMODU2=ABS(XX(I))**2
CPP +
          IF(I.NE.1.AND.I.NE.KHALF1)THEN
            YCOMP=CONJG(XX(NPOINT-I+2))
            PMODU2=PMODU2+XX(NPOINT-I+2)*YCOMP
cbp: on pourrait aussi ecrire   PMODU2=PMODU2+ABS(XX(NPOINT-I+2))**2
          ENDIF
CPP +
        MLREE2.PROG(I-KDEBU+1)=COEF*PMODU2
C       IF(IIMPI.EQ.1) WRITE(IOIMP,*)'DSP(',I,')=',MLREE2.PROG(I)
  20    CONTINUE
C
      SEGDES MLREE1,MLREE2
C
C     CREATION DE L'OBJET EVOLUTIO DSP
C
      N=1
      SEGINI MEVOLL
      IPVO=MEVOLL
      TI(1:72)=TITREE
      IEVTEX=TI
      ITYEVO='REEL'
      SEGINI KEVOLL
c       KEVTEX=TI
      KEVTEX=KEVOL1.KEVTEX
      IEVOLL(1)=KEVOLL
      TYPX='LISTREEL'
      TYPY='LISTREEL'
C
      IPROGX=MLREE1
      NOMEVX=MOTX(1:12)
C
      IPROGY=MLREE2
      NOMEVY=MOTY(1:12)
C
      NUMEVX=ICOUL
      NUMEVY='REEL'
C
      SEGSUP MTRAV1,MTRAV2
      SEGDES KEVOL1
      SEGDES MEVOL1
C
      SEGDES KEVOLL,MEVOLL
      CALL ECROBJ('EVOLUTIO',IPVO)
 666  CONTINUE
      RETURN
      END