Sources Esope | Cast3M

tfrinv
C TFRINV    SOURCE    CB215821  19/07/30    21:18:14     10273          
      SUBROUTINE TFRINV
C
C=======================================================================
C  =  CALCUL DE LA TRANSFORMEE DE FOURIER INVERSE D'UN SIGNAL
C  =
C  =  SYNTAXE : FFTR = TFRI EVOL    ;
C  =
C  =  EVOL    : OBJET DE TYPE EVOLUTIO CONTENANT L'OBJET EVOLUTION
C  =                                                                  =
C  =  CREATION    : 26/04/88                                          =
C  =  F. ROULLIER
C=======================================================================
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT  REAL*8(A-H,O-Z)
      CHARACTER *72 TI
      CHARACTER *4  ITYP1
c     RLTOCX : fonction de conversion real*8 -> complex*16  
      COMPLEX*16 RLTOCX
      LOGICAL INV
 
-INC PPARAM
-INC CCOPTIO
-INC SMEVOLL
-INC SMLREEL
-INC CCREEL
C
c       SEGMENT COMPL
c           COMPLEX*16 XX(NPT),YY(NPT),W(NEXP)
c           REAL*8     AR(NPOIN),AI(NPOIN)
c           REAL*8     DM(NPTD),DP(NPTD)
c       ENDSEGMENT
 
c     segment pour appel a fftpack5.1      
      segment wfft51
      real*8 work(lenwrk)
      real*8 wsave(lensav)
      real*8 XX51(NCOU)
      endsegment
C
C==== LECTURES =========================================================
c
C     LECTURE DE L'OBJET EVOLUTION COMPLEXE A TRANSFORMER
C
      CALL LIROBJ('EVOLUTIO',IEV1,1,IRETOU)
      IF(IERR.NE.0) RETURN
      CALL ACTOBJ('EVOLUTIO',IEV1,1)
C
      MEVOLL=IEV1
      IF (ITYEVO.NE.'COMPLEXE') THEN
        MOTERR(1:8)='EVOLUTION'
        CALL ERREUR (302)
        RETURN
      ENDIF
      KEVOL1=IEVOLL(1)
      KEVOL2=IEVOLL(2)
      ITYP1=KEVOL1.NUMEVY
      ICOU1=KEVOL1.NUMEVX
      MLREE3=KEVOL1.IPROGX
      MLREE1=KEVOL1.IPROGY
      MLREE2=KEVOL2.IPROGY
C
c     on recupere le pas en frequence, la dimension et les ordonnees
      L1=MLREE3.PROG(/1)
      DF=MLREE3.PROG(2)-MLREE3.PROG(1)
C
c     on verifie le nombre de points
      NPT=(L1-1)*2
      NP=NPT
      NN=0
      DO 1 I=1,1000
      NP=NP/2
      NN=NN+1
      IF (NP.LE.1) GO TO 2
    1 CONTINUE
    2 NP=2**NN
      IF (NP.NE.NPT) THEN
         MOTERR(1:8)='EVOLUTIO'
         CALL ERREUR (1083)
         RETURN
      ENDIF
c     on verifie -que la 1ere frequence soit 0,
      IF(MLREE3.PROG(1).NE.0.D0) THEN
c La premiere frequence doit etre egale a 0.
         REAERR(1)=0.d0
         CALL ERREUR (1084)
         RETURN
      ENDIF
C                -que la partie imag de la 1ere soit 0 ...
      IF(ITYP1.EQ.'PREE') THEN
        DR=MLREE1.PROG(1)
        DI=MLREE2.PROG(1)
        DM=SQRT(DR**2+DI**2)
      ELSE
        DM = MLREE1.PROG(1)
        DP = MLREE2.PROG(1) * XPI/180.     
        DI = DM*SIN(DP)
      ENDIF       
      IF(ABS(DI).GT.(XZPREC*DM)) THEN
c La 1 eme valeur imaginaire doit etre egale a 0.
         INTERR(1)=1
         REAERR(1)=0.d0
         CALL ERREUR (1086)        
         RETURN
      ENDIF
C
c       NEXP=NPT/2
c       NPOIN=L1
c       IF (ITYP1.NE.'PREE') THEN
c         NPTD=L1
c       ELSE
c         NPTD=0
c       ENDIF
c       SEGINI COMPL
c       write(*,*) 'TFR en entree de dim:',L1
c       write(*,*) 'nombre de point en sortie',NPT,NEXP
c       write(*,*) 'dim des tableaux de travail',DM(/1),AR(/1),W(/1)
 
C     creation du tableau de travail pour FFTPACK5.1
      NCOU=NPT
      lenwrk = 2 * NPT
      lensav = NPT + int(log ( real (NPT) ) / log (2.0) ) + 4
      segini,wfft51      
C
C==== CHARGEMENT DES TABLEAUX DE TRAVAIL ===============================
C
c   - cas Reel,Imag :
      IF (ITYP1.EQ.'PREE') THEN
 
c       on a forcement le terme constant (f=0Hz)
        XX51(1) = MLREE1.PROG(1)
c       la suite par paire
        ITY=1
        DO 10 I=2,L1-1
          ITY=ITY+1
          XX51(2*I-2)=2. * (-1)**(2*I)   * MLREE1.PROG(ITY)
          XX51(2*I-1)=2. * (-1)**(2*I+1) * MLREE2.PROG(ITY)
  10    CONTINUE
c       dernier terme : on suppose que la TF n'a pas ete tronquee
        XX51(NPT) =  MLREE1.PROG(L1)
 
c   - cas Module, Phase :
      ELSE 
 
        RAP=XPI/180.
c       on a forcement le terme constant (f=0Hz) 
c       XX51(1) =  MLREE1.PROG(1)
        DM = MLREE1.PROG(1)
        DP = MLREE2.PROG(1) * RAP
        DR = DM*COS(DP)
        DI = DM*SIN(DP)
c         write(*,*) '1: DM,DP,DR,DI=',DM,DP,DR,DI
        XX51(1) = DR
c       la suite par paire
        ITY=1
        DO 11 I=2,L1-1
          ITY=ITY+1
c         MODULE , PHASE -> REEL, IMAG
          DM = MLREE1.PROG(ITY)
          DP = MLREE2.PROG(ITY) * RAP
          DR = DM*COS(DP)
          DI = DM*SIN(DP)
c           write(*,*) ITY,I,': DM,DP,DR,DI=',DM,DP,DR,DI
          XX51(2*I-2)=2.0 * (-1)**(2*I)   * DR
          XX51(2*I-1)=2.0 * (-1)**(2*I+1) * DI
  11    CONTINUE
c       dernier terme : on suppose que la TF n'a pas ete tronquee
c       XX51(NPT) = 0.25d0 * MLREE1.PROG(L1)
        DM = MLREE1.PROG(L1)
        DP = MLREE2.PROG(L1) * RAP
        DR = DM*COS(DP)
        DI = DM*SIN(DP)
        XX51(NPT) =  DR
 
c c       si la TF est fournie sous la forme (module,phase), conversion
c         DO 8 I=1,L1
c           DM(I)=MLREE1.PROG(I)
c           DP(I)=MLREE2.PROG(I)
c     8   CONTINUE
c         CALL CONVCP(AR,AI,DM,DP,L1,-1)
c       ELSE
c c       sinon on recupere directement la TF (Re,Im)
c         DO 9 I=1,L1
c           AR(I)=MLREE1.PROG(I)
c           AI(I)=MLREE2.PROG(I)
c     9   CONTINUE
 
      ENDIF
 
 
C==== CALCUL DE LA FFT inverse =========================================
 
c +++ via FFTPACK5.1 +++
 
c     Initialize the WSAVE array.
      call rfft1i (NPT,wsave,lensav,ier)
      IF (IERR.ne.0) RETURN
c
c     Compute inverse FFT of coefficients.  Should get back the
c     original data.
      inc = 1
      lenr = NPT
      call rfft1b(NPT,inc,XX51,lenr,wsave,lensav,work,lenwrk,ier)
      IF (IERR.ne.0) RETURN
 
  19  IF(IIMPI.EQ.1) WRITE(IOIMP,*)' FFT INVERSE CALCULEE '
 
C
C==== CREATION ET CALCUL DES LISTES DE LA FFT ==========================
C
      JG=NPT
      SEGINI MLREE1
      SEGINI MLREE2
      FMAX=DF*REAL(NPT)
      DTEMP=1.D0/FMAX
C
      DO 20 I=1,NPT
        TEMP=REAL(I-1)*DTEMP
        MLREE1.PROG(I)=TEMP
        MLREE2.PROG(I)=DF*XX51(I)
  20  CONTINUE
      SEGSUP wfft51
C
C                CREATION DE L'OBJET EVOLUTION REEL FFT INVERSE
C
      N=1
      SEGINI MEVOLL
      IPVO=MEVOLL
      IEVTEX=TITREE
      ITYEVO='REEL    '
C
      SEGINI KEVOLL
c       KEVTEX=TITREE
      KEVTEX='TF^{-1}'
      TYPX='LISTREEL'
      TYPY='LISTREEL'
      IEVOLL(1)=KEVOLL
      IPROGX=MLREE1
      NOMEVX='TEMPS    SEC'
      IPROGY=MLREE2
      NOMEVY='TF INVERSE  '
      NUMEVX=ICOU1
      NUMEVY='REEL'
C
      CALL ACTOBJ('EVOLUTIO',IPVO,1)
      CALL ECROBJ('EVOLUTIO',IPVO)
  666 CONTINUE
 
      END