Télécharger la plaquette Cast3M


fr en

Télécharger tfor.eso

Retour à la liste

Numérotation des lignes : 

tfor

  1. C TFOR      SOURCE    GOUNAND   25/05/05    21:15:08     12259          
  2.  
  3.       SUBROUTINE TFOR(IOPT)
  4.  
  5.       IMPLICIT INTEGER(I-N)
  6.       IMPLICIT REAL*8(A-H,O-V)
  7.       IMPLICIT COMPLEX*16 (W-Z)
  8. C
  9. C=======================================================================
  10. C  CALCULS IMPLIQUANT UNE TRANSFORMEE DE FOURIER RAPIDE D'UN SIGNAL
  11. C=======================================================================
  12. C
  13. C  +-----------+------+----------------------------------------------+
  14. C  | OPERATEUR | IOPT |                   CALCUL                     |
  15. C  +-----------+------+----------------------------------------------+
  16. C  |   TFR     |  +1  | TRANSFORMEE DE FOURIER RAPIDE D'UN SIGNAL    |
  17. C  +-----------+------+----------------------------------------------+
  18. C  |   DSPR    |  +2  | DENSITE SPECTRALE DE PUISSANCE D'UN SIGNAL   |
  19. C  +-----------+------+----------------------------------------------+
  20. C
  21. C  SYNTAXE : voir notices des operateurs
  22. C  -------
  23. C
  24. C  NOMENCLATURE DE LA SOURCE :
  25. C  -------------------------
  26. C    EXP2    : EXPOSANT DANS NPOINT=2^EXP2, NPOINT ETANT
  27. C              LE NOMBRE DE REELS DANS LISTREEL
  28. C    EVO1    : OBJET DE TYPE EVOLUTIO CONTENANT LE SIGNAL A TRAITER
  29. C    SORT    : TYPE DE SORTIE ;
  30. C               = 'REIM'   PART REEL & PART IMAG / FREQ
  31. C               = 'MOPH'   MODULE    & PHASE     / FREQ
  32. C    FMIN    : MOT-CLE
  33. C    V1      : FREQUENCE MINIE A VISUALISER
  34. C    FMAX    : MOT-CLE
  35. C    V2      : FREQUENCE MAXIE A VISUALISER
  36. C    COU1    : COULEUR ATTRIBUEE A LA PREMIERE COURBE (FACULTATIF)
  37. C    COU2    : COULEUR ATTRIBUEE A LA DEUXIEME COURBE (FACULTATIF)
  38. C
  39. C  CREATION/MODIF :
  40. C  --------------
  41. C  - 16/04/87, BEAUFILS
  42. C  - 20/05/2015, Benoit Prabel : extension aux LISTCHPO,
  43. C  - 2018-10-01, Benoit Prabel : branchement de FFTPACK5.1
  44. c  - 2021-08-04, BP : aiguillage de DSPR vers tfor.eso pour mutualiser
  45. C
  46. c  TODO :
  47. c  ----
  48. c  - brancher aussi les lischpo a FFTPACK
  49. c    en inversant les boucles + //iser sur les inconnues
  50. C
  51. C=======================================================================
  52. C
  53. -INC CCGEOME
  54. -INC PPARAM
  55. -INC CCOPTIO
  56. -INC CCREEL
  57. -INC SMEVOLL
  58. -INC SMLREEL
  59. -INC SMLCHPO
  60. -INC SMCHPOI
  61. -INC SMCOORD
  62. C
  63. c     segments pour appel a WEXP + FFTL
  64.       SEGMENT MTRAV2
  65.           IMPLIED W(NEXP)
  66.       ENDSEGMENT
  67.  
  68.       SEGMENT MTRAV3
  69.           IMPLIED XXX(NCOU,NDDL),YYY(NCOU)
  70.       ENDSEGMENT
  71.  
  72.       SEGMENT,LIPOV(NIPOV)
  73.       SEGMENT,LIPOV2(NIPOV)
  74.       SEGMENT,LIPOV3(NIPOV)
  75.  
  76. c     segment pour appel a fftpack5.1
  77.       segment wfft51
  78.       real*8 work(lenwrk)
  79.       real*8 wsave(lensav)
  80.       real*8 XX51(NCOU)
  81.       endsegment
  82. C
  83.       CHARACTER *72 TI
  84.       CHARACTER *24 MOTY
  85.       CHARACTER *4  NSORT(2),MCLE(2)
  86.       LOGICAL INV
  87.  
  88.       DATA NSORT(1),NSORT(2)/'REIM','MOPH'/
  89.       DATA MCLE(1),MCLE(2)/'FMIN','FMAX'/
  90. C
  91. C==== INITIALISATIONS ==================================================
  92. C
  93.  
  94.       IPSIG=0
  95.       IPTAB=0
  96. c     FFT toujours directe ici
  97.       INV=.FALSE.
  98.  
  99. C
  100. C==== LECTURES =========================================================
  101. C
  102. C     LECTURE DE EXP2
  103. C
  104.       CALL LIRENT(N2,1,IRET1)
  105.       IF(IRET1.EQ.0)GOTO 666
  106.       NPOINT=2**N2
  107.       IF(IIMPI.EQ.1) write(ioimp,*) 'NPOINT=2**',N2,'=',NPOINT
  108.       IF(N2 .LT. 2)THEN
  109.         CALL ERREUR(21)
  110.         GOTO 666
  111.       ENDIF
  112. C
  113. C     LECTURE DE L'OBJET EVOLUTIO CONTENANT LE SIGNAL
  114. C     ou de L'OBJET LISTCHPO contenant la suite de chpoint + liste des temps
  115. C
  116.       CALL LIROBJ('EVOLUTIO',IPSIG,0,IRET2)
  117.       IF(IRET2.EQ.0) THEN
  118.         CALL LIROBJ ('LISTCHPO',IPLCH,1,IRET2)
  119.         IF(IRET2.EQ.0)GOTO 666
  120.         CALL LIROBJ ('LISTREEL',IPREE1,1,IRET2)
  121.         IF(IRET2.EQ.0)GOTO 666
  122.       ENDIF
  123.       IF(IIMPI.EQ.1) write(ioimp,*) 'signal=',IPSIG,IPLCH,IPREE1
  124. C
  125. C     TYPE DE SORTIE : LECTURE si TFR de 'REel/IMaginaire' ou 'MOdule/PHase'
  126. C                      necessairement MODULE**2 pour DSPR
  127.       IF(IOPT.EQ.1) THEN
  128.         CALL LIRMOT(NSORT,2,ISORT,0)
  129.         IF(ISORT.EQ.0) ISORT=1
  130.       ELSEIF(IOPT.EQ.2) THEN
  131.         ISORT=3
  132.       ELSE
  133.         WRITE(IOIMP,*) 'tfor: valeur de IOPT incorrect'
  134.         CALL ERREUR(5)
  135.       ENDIF
  136.       IF(IIMPI.EQ.1) write(ioimp,*) 'sortie :',ISORT
  137. C
  138. C     LECTURE DE LA FREQUENCE MINI FMIN -> FRMI
  139. C
  140.       CALL LIRMOT(MCLE(1),1,IRET,0)
  141.       IF(IRET.EQ.1) THEN
  142.           CALL LIRREE(FRMI,0,IRET1)
  143.           IF(IRET1.EQ.0) THEN
  144.               MOTERR=MCLE(1)(1:4)
  145.               CALL ERREUR(166)
  146.               GOTO 666
  147.           ENDIF
  148.           IF(FRMI.LT.(-1.D-20)) FRMI=-FRMI
  149.           IF(ABS(FRMI).LT.1.D-20) FRMA=-1.D0
  150.       ELSE
  151.           FRMI=-1.D0
  152.       ENDIF
  153. C
  154. C     LECTURE DE LA FREQUENCE MAXI FMAX -> FRMA
  155. C
  156.       CALL LIRMOT(MCLE(2),1,IRET,0)
  157.       IF(IRET.EQ.1) THEN
  158.           CALL LIRREE(FRMA,0,IRET1)
  159.           IF(IRET1.EQ.0) THEN
  160.               MOTERR=MCLE(2)(1:4)
  161.               CALL ERREUR(166)
  162.               GOTO 666
  163.           ENDIF
  164.           IF(FRMA.LT.(-1.D-20)) FRMA=-FRMA
  165.           IF(ABS(FRMA).LT.1.D-20) FRMA=-1.D0
  166.       ELSE
  167.           FRMA=-1.D0
  168.       ENDIF
  169.       IF(IIMPI.EQ.1) write(ioimp,*) FRMI,' < FREQ <',FRMA
  170. C
  171. C     LECTURE DE LA COULEUR
  172. C
  173.       CALL LIRMOT(NCOUL,NBCOUL,ICOU1,0)
  174.       IF (ICOU1.EQ.0) ICOU1=IDCOUL+1
  175.       ICOU1=ICOU1-1
  176.       CALL LIRMOT(NCOUL,NBCOUL,ICOU2,0)
  177.       IF (ICOU2.EQ.0) ICOU2=IDCOUL+1
  178.       ICOU2=ICOU2-1
  179.       IF(IIMPI.EQ.1) write(ioimp,*) 'couleurs=',ICOU1,ICOU2
  180. C
  181.       IF(IERR.NE.0) GOTO 666
  182. C
  183. C     VERIFICATION DU NOMBRE DE PAS DE TEMPS ...
  184. C
  185.       IF(IPSIG.NE.0) THEN
  186. C     ...POUR L'EVOLUTION EN ENTREE
  187.         MEVOL1=IPSIG
  188.         CALL ACTOBJ('EVOLUTIO',IPSIG,1)
  189.         IF(MEVOL1.IEVOLL(/1).NE.1) THEN
  190.           CALL ERREUR(687)
  191.         ENDIF
  192.         KEVOL1=MEVOL1.IEVOLL(1)
  193.         MLREE1=KEVOL1.IPROGX
  194.         MLREE2=KEVOL1.IPROGY
  195.         L1    =MLREE1.PROG(/1)
  196.         IF(L1 .GE. 2)THEN
  197.           DT=MLREE1.PROG(2)-MLREE1.PROG(1)
  198.         ELSE
  199.           CALL ERREUR(199)
  200.           RETURN
  201.         ENDIF
  202.         IF(DT.LE.XPETIT) THEN
  203.           CALL ERREUR(206)
  204.           GOTO 666
  205.         ENDIF
  206.       ELSE
  207. C     ...POUR LA LISTE DE CHPOINTs EN ENTREE
  208.         MLCHPO=IPLCH
  209.         SEGACT,MLCHPO
  210.         N1=ICHPOI(/1)
  211.         MLREE1=IPREE1
  212.         SEGACT MLREE1
  213.         L1=MLREE1.PROG(/1)
  214.         IF(L1 .GE. 2)THEN
  215.           DT=MLREE1.PROG(2)-MLREE1.PROG(1)
  216.         ELSE
  217.           CALL ERREUR(199)
  218.           RETURN
  219.         ENDIF
  220.         IF(IIMPI.EQ.1) write(ioimp,*) 'N1,L1,DT=',N1,L1,DT
  221.         IF(DT.LE.XPETIT) THEN
  222.           CALL ERREUR(206)
  223.           GO TO 666
  224.         ENDIF
  225.         IF(N1.LT.L1) THEN
  226.           CALL ERREUR(217)
  227.           goto 666
  228.         ENDIF
  229.       ENDIF
  230.  
  231. c     qq ecritures pour le debuggage
  232.       IF(IIMPI.EQ.1) THEN
  233.  
  234.         PRINT *,'NPOINT=',NPOINT,'L1=',L1
  235.         IF(NPOINT.GT.L1) THEN
  236.           WRITE(IOIMP,1000) L1,N2,NPOINT
  237.  1000     FORMAT(1H ,'LE NOMBRE DE POINTS DANS LISTEMPS ',I6, ' EST ',
  238.      &    'INFERIEURE @ 2**',I5,
  239.      &    /' ON PRENDRA UNE LONGUEUR DE LISTEMPS DE ',I6,' MOTS     ',
  240.      &    /' ET ON COMPLETERA PAR DES ZEROS')
  241.         ELSE
  242.            IF(NPOINT.LT.L1) THEN
  243.                WRITE(IOIMP,1010) N2
  244.  1010          FORMAT(1H ,'ON TRONQUE LE SIGNAL A 2**',I5,' MOTS',/)
  245.            ELSE
  246.                WRITE(IOIMP,1030)N2,NPOINT
  247.  1030          FORMAT(1H ,'LA LONGUEUR DU LISTEMP EST EGALE A 2**',I5,
  248.      &        ' = ',I6,/)
  249.            ENDIF
  250.         ENDIF
  251.       ENDIF
  252.  
  253. C     CALCUL DE QUELQUES VARIABLES UTILES
  254. cbp : inutile?      NDIBLK=NPOINT
  255. C     IND1 = indice max pour le remplissage temporel
  256.       IND1=L1
  257.       IF(NPOINT.LT.L1) IND1 = NPOINT
  258.       DUREE=DT*REAL(NPOINT)
  259.       DUR05=0.5D0*DUREE
  260.       DFREQ=1.D0/DUREE
  261.       KHALF=INT(NPOINT/2)+1
  262.       KMIN=INT(FRMI/DFREQ)+1
  263.       KMAX=INT(FRMA/DFREQ)+1
  264.       KDEBU=1
  265.       IF(KMIN.GT.0) KDEBU=KMIN
  266.       IF((KMAX.LT.KHALF).AND.(FRMA.GT.0.D0)) THEN
  267.         KHALF=KMAX
  268.         KFIN =KHALF
  269.       ELSE
  270.         KFIN =KHALF-1
  271.       ENDIF
  272.       NNN=KHALF-KDEBU+1
  273.  
  274. c     si pas d'evolution (alors listchpo) -> on go en 700
  275.       IF(IPSIG.eq.0) GOTO 700
  276.  
  277.  
  278. C=======================================================================
  279. C     CAS EVOLUTION
  280. C=======================================================================
  281.       IF(IIMPI.EQ.1) write(ioimp,*) '=== CAS EVOLUTION ==='
  282. C
  283. C==== CHARGEMENT DES TABLEAUX DE TRAVAIL ===============================
  284. C
  285. C     creation du tableau de travail pour FFTPACK5.1
  286.       NCOU=NPOINT
  287.       lenwrk = 2 * NPOINT
  288.       lensav = NPOINT + int(log ( real (NPOINT) ) / log (2.0) ) + 4
  289.       segini,wfft51
  290.       DO 10 I=1,IND1
  291.         XX51(I)=MLREE2.PROG(I)
  292. c         write(*,FMT="(I6,A,E24.16)") I,' ',XX51(I)
  293.   10  CONTINUE
  294. C     On complete avec des 0 (utile meme apres SEGINI)
  295.       IF(NPOINT.GT.L1) THEN
  296.           L2=L1+1
  297.           DO 11 I=L2,NPOINT
  298.              XX51(I)=0.d0
  299.   11      CONTINUE
  300.       ENDIF
  301.       do ii=1,lenwrk
  302.       work(ii)=0.d0
  303.       enddo
  304.       do iii=1,lensav
  305.       wsave(iii)=0.d0
  306.       enddo
  307. C
  308. C==== CALCUL DE LA FFT =================================================
  309. C
  310.       IF(IIMPI.EQ.1) WRITE(IOIMP,*)' CALCUL DE LA FFT DU SIGNAL '
  311.  
  312. c +++ via FFTPACK5.1 +++
  313.  
  314. c     Initialize the WSAVE array.
  315.       call rfft1i (NPOINT,wsave,lensav,ier)
  316.       IF (IERR.ne.0) RETURN
  317. c       do iou=1,lensav
  318. c         write(*,FMT="(I6,A,E24.16)") iou,' ',wsave(iou)
  319. c       enddo
  320. c
  321. c     Compute the FFT coefficients.
  322.       inc = 1
  323.       lenr = NPOINT
  324. c       do iou=1,NPOINT
  325. c         write(*,*) 'XX51_avant (',iou,')=',XX51(iou)
  326. c       enddo
  327.       call rfft1f (NPOINT,inc,XX51,lenr,wsave,lensav,work,lenwrk,ier)
  328.       IF (IERR.ne.0) RETURN
  329.       IF(IIMPI.EQ.1) WRITE(IOIMP,*)' FFT CALCULEE '
  330. C
  331. C==== CREATION ET REMPLISSAGE DES LISTES DE LA FFT =====================
  332. C
  333.       IF(IIMPI.EQ.1) WRITE(IOIMP,*)' CHARGEMENT DE LA FFT '
  334.       JG=NNN
  335.       SEGINI MLREE1,MLREE2,MLREE3
  336. C
  337. C==== PARTIE REELLE & IMAGINAIRE / FREQUENCE (par defaut) ==============
  338. C
  339.           ITY=0
  340. c         souhaite-t-on le terme constant (f=0Hz) ?
  341.           IF (KDEBU.EQ.1) THEN
  342. c             WRITE(*,*) 'on remplit: 1 (ordre 0)'
  343.             MLREE1.PROG(1)=0.d0
  344.             MLREE2.PROG(1)=DUREE*XX51(1)
  345.             MLREE3.PROG(1)=0.d0
  346.             KDEBU=KDEBU+1
  347.             ITY=ITY+1
  348.           ENDIF
  349. c           DO 20 I=KDEBU,(KHALF-2)
  350.           DO 20 I=KDEBU,KFIN
  351.             FREQ=REAL(I-1)*DFREQ
  352.             ITY=ITY+1
  353. c             WRITE(*,*) 'LOOP20: on remplit: ',ITY,'(ordre',I,')f=',FREQ
  354.             MLREE1.PROG(ITY)=FREQ
  355.             MLREE2.PROG(ITY)=(-1)**(2*I)  *DUR05*XX51(2*I-2)
  356.             MLREE3.PROG(ITY)=(-1)**(2*I+1)*DUR05*XX51(2*I-1)
  357.   20      CONTINUE
  358. c         derniere valeur seulement si :
  359. c         + on a un nombre pair en entree (tjrs vrai car NPOINT=2**N2)
  360. c         + on cherche la dernier valeur
  361.           IF (KFIN.EQ.(KHALF-1)) THEN
  362.             FREQ=FREQ+DFREQ
  363.             ITY=ITY+1
  364. c           WRITE(*,*) 'on remplit: ',ITY,'(ordre',I,')f=',FREQ
  365.             MLREE1.PROG(ITY)=FREQ
  366.             MLREE2.PROG(ITY)=DUREE*XX51(NPOINT)
  367.             MLREE3.PROG(ITY)=0.D0
  368.           ENDIF
  369. c           MOTY(1:24)='PART REEL   PART IMAG'
  370. C
  371. C==== SORTIE EN MODULE & PHASE / FREQUENCE : on retraite ===============
  372. C
  373.       IF (ISORT.EQ.2) THEN
  374.           DO 21 ITY=1,NNN
  375.             PREEL=MLREE2.PROG(ITY)
  376.             PIMAG=MLREE3.PROG(ITY)
  377. c           module
  378.             PMODU=SQRT(PREEL**2+PIMAG**2)
  379.             MLREE2.PROG(ITY)=PMODU
  380. c           phase
  381.             CALL FACOMP(PREEL,PIMAG,PMODU)
  382.             MLREE3.PROG(ITY)=PMODU
  383.   21      CONTINUE
  384. c           MOTY(1:8)='MODULE'
  385. c           MOTY(9:12)=' '
  386. c           MOTY(13:20)='PHASE'
  387. c           MOTY(21:24)=' '
  388.       ENDIF
  389. C
  390. C==== SORTIE DSPR : on retraite un peu différemment ===============
  391. C
  392.       IF (ISORT.EQ.3) THEN
  393.           DFREQ2=2.D0*DFREQ
  394.           DO 22 ITY=1,NNN
  395.             PREEL=MLREE2.PROG(ITY)
  396.             PIMAG=MLREE3.PROG(ITY)
  397.             PMODU=PREEL**2+PIMAG**2
  398.             MLREE2.PROG(ITY)=DFREQ2*PMODU
  399.   22      CONTINUE
  400.       ENDIF
  401.  
  402.       SEGSUP wfft51
  403.  
  404. C
  405. C==== CREATION DE L'OBJET EVOLUTIO FFT =================================
  406. C
  407. c     longueur du titre de l'evol en entree
  408.       CALL LENCHA(KEVOL1.KEVTEX,LEN1)
  409.       LEN1=min(LEN1,65)
  410.  
  411. c     evol de sortie
  412.       IF (ISORT.EQ.3) THEN
  413.         N=1
  414.       ELSE
  415.         N=2
  416.       ENDIF
  417.       SEGINI MEVOLL
  418.       IPVO=MEVOLL
  419.       TI(1:72)=TITREE
  420.       IEVTEX=TI
  421.       IF (ISORT.EQ.3) THEN
  422.         ITYEVO='REEL'
  423.       ELSE
  424.         ITYEVO='COMPLEXE'
  425.       ENDIF
  426. C
  427. C  -- PREMIERE SOUS-COURBE --
  428. c
  429.       SEGINI KEVOLL
  430.       IEVOLL(1)=KEVOLL
  431. C     abscisses
  432.       TYPX='LISTREEL'
  433.       IPROGX=MLREE1
  434.       NOMEVX='f (Hz)'
  435. C     ordonnees
  436.       TYPY='LISTREEL'
  437.       IPROGY=MLREE2
  438. c       NOMEVY=MOTY(1:12)
  439. C     couleur
  440.       NUMEVX=ICOU1
  441. c     titre des abscisses, type et titre de la courbe (legende)
  442.       IF(ISORT.EQ.1) THEN
  443.           NOMEVY='Re(X)'
  444.           NUMEVY='PREE'
  445.           KEVTEX='Re TF('//KEVOL1.KEVTEX(1:LEN1)//')'
  446.       ELSEIF(ISORT.EQ.2) THEN
  447.           NOMEVY='|X|'
  448.           NUMEVY='MODU'
  449.           KEVTEX='Amp TF('//KEVOL1.KEVTEX(1:LEN1)//')'
  450.       ELSEIF(ISORT.EQ.3) THEN
  451.           NOMEVY='S_{x}'
  452. c         bp : on omet l'unite = unite_de_x**2/Hz
  453.           NUMEVY='REEL'
  454.           KEVTEX='DSP('//KEVOL1.KEVTEX(1:LEN1)//')'
  455.       ENDIF
  456. C
  457. C  -- DEUXIEME SOUS-COURBE --
  458. c
  459.       IF (N.GE.2) THEN
  460.         SEGINI KEVOLL
  461.         IEVOLL(2)=KEVOLL
  462. C
  463.         TYPX='LISTREEL'
  464.         IPROGX=MLREE1
  465.         NOMEVX='f (Hz)'
  466. C
  467.         TYPY='LISTREEL'
  468.         IPROGY=MLREE3
  469. C
  470.         NUMEVX=ICOU2
  471.         IF(ISORT.EQ.1) THEN
  472.             NOMEVY='Im(X)'
  473.             NUMEVY='PIMA'
  474.             KEVTEX='Im TF('//KEVOL1.KEVTEX(1:LEN1)//')'
  475.         ELSEIF(ISORT.EQ.2) THEN
  476.             NOMEVY='\j(X)'
  477.             NUMEVY='PHAS'
  478.             KEVTEX='\j TF('//KEVOL1.KEVTEX(1:LEN1)//')'
  479.         ENDIF
  480.       ELSE
  481.         SEGSUP,MLREE3
  482.       ENDIF
  483. C
  484. C  -- ECRITURE --
  485.       CALL ACTOBJ('EVOLUTIO',IPVO,1)
  486.       CALL ECROBJ('EVOLUTIO',IPVO)
  487.  666  CONTINUE
  488.       RETURN
  489.  
  490.  
  491.  700  CONTINUE
  492.  
  493. C=======================================================================
  494. C     CAS LISTCHPO
  495. C=======================================================================
  496.       IF(IIMPI.EQ.1) write(ioimp,*) '=== CAS LISTCHPO ==='
  497. C
  498. C==== CALCUL DE L'EXPONENTIEL ==========================================
  499. C
  500.       IF(IIMPI.EQ.1) WRITE(IOIMP,*)' CALCUL DE L EXPONENTIEL'
  501.       NEXP=NPOINT/2
  502.       SEGINI MTRAV2
  503.       CALL WEXP(INV,NPOINT,W(1))
  504.       IF(IIMPI.EQ.1) write(ioimp,*) 'W=',(W(iou),iou=1,NEXP)
  505.  
  506. C     creation des listchpo de sortie
  507.       N1=NNN
  508.       SEGINI,MLCHP2,MLCHP3
  509.  
  510. C---- boucle sur les ddls des chpoints----------------------------------
  511. c     on suppose tous les chpoints de la meme forme que le 1er !
  512.  
  513. C     on stocke les MPOVAL dans le tableau LIPOV
  514.       NIPOV =IND1
  515.       segini,LIPOV
  516.       NIPOV =NNN
  517.       segini,LIPOV2,LIPOV3
  518.  
  519. C==== BOUCLE SUR LES MSOUPO DES CHPOINTS ===============================
  520.       IFOCHS = -99
  521.       ISOUPO=0
  522.  710  ISOUPO=ISOUPO+1
  523.       IF(IIMPI.EQ.1) write(ioimp,*) '-------- zone',ISOUPO,' --------'
  524.  
  525. c---- on ouvre les MSOUPO de tous les chpoints ----
  526. c       DO I1=1,L1
  527.       DO I1=1,IND1
  528.         MCHPOI=ICHPOI(I1)
  529.         CALL ACTOBJ('CHPOINT',MCHPOI,1)
  530.         IF(ISOUPO.EQ.1) THEN
  531. c           IF(IIMPI.EQ.1) write(ioimp,*) 'Time',I1,': on ouvre',MCHPOI
  532.           NSOUPO=IPCHP(/1)
  533.           IF(I1.EQ.1) NSOUP1=NSOUPO
  534.           IF(NSOUP1.NE.NSOUPO) THEN
  535.             CALL ERREUR(214)
  536.             return
  537.           ENDIF
  538.           IF (I1.EQ.1) IFOCHS = MCHPOI.IFOPOI
  539.           IF (MCHPOI.IFOPOI.NE.IFOCHS) THEN
  540.             moterr(1:8)='CHPOINT'
  541.             interr(1)=IFOCHS
  542.             interr(2)=MCHPOI.IFOPOI
  543.             interr(3)=IFOUR
  544. c-dbg           write(ioimp,*) '1132 TFOR',IFOCHS,MCHPOI.IFOPOI
  545.             call erreur(1132)
  546.             IFORIG = IFOUR
  547.           ENDIF
  548.         ENDIF
  549.         MSOUPO   = IPCHP(ISOUPO)
  550.         MPOVAL   = IPOVAL
  551.         LIPOV(I1)= MPOVAL
  552.         NC       = VPOCHA(/2)
  553.         N        = VPOCHA(/1)
  554.         IF(I1.EQ.1) THEN
  555.           NC1=NC
  556.           N1=N
  557.         ENDIF
  558.         IF(NC1.NE.NC.OR.N1.NE.N) THEN
  559.           CALL ERREUR(214)
  560.           return
  561.         ENDIF
  562. c       rem : on pourrait aussi verifier les noms de composantes, etc...
  563.       ENDDO
  564.  
  565. C---- creation du chpoint de sortie ----
  566.       DO I1=1,NNN
  567.         IF(ISOUPO.EQ.1) THEN
  568.           NAT=2
  569.           SEGINI,MCHPO2,MCHPO3
  570.           MLCHP2.ICHPOI(I1)=MCHPO2
  571.           MCHPO2.MTYPOI='REEL   '
  572.           MCHPO2.MOCHDE='CHAMP PAR POINTS CREE PAR TFR '
  573.           MCHPO2.IFOPOI=IFOCHS
  574.           MLCHP3.ICHPOI(I1)=MCHPO3
  575.           MCHPO3.MTYPOI='IMAG   '
  576.           MCHPO3.MOCHDE='CHAMP PAR POINTS CREE PAR TFR '
  577.           MCHPO3.IFOPOI=IFOCHS
  578.         ELSE
  579.           MCHPO2=MLCHP2.ICHPOI(I1)
  580.           MCHPO3=MLCHP3.ICHPOI(I1)
  581.         ENDIF
  582. C       creation du MSOUP2 et du MPOVA2 de sortie
  583.         SEGINI,MSOUP2=MSOUPO
  584.         MCHPO2.IPCHP(ISOUPO)=MSOUP2
  585.         SEGINI,MPOVA2=MPOVAL
  586.         MSOUP2.IPOVAL=MPOVA2
  587.         LIPOV2(I1)=MPOVA2
  588.         SEGINI,MSOUP3=MSOUPO
  589.         MCHPO3.IPCHP(ISOUPO)=MSOUP3
  590.         SEGINI,MPOVA3=MPOVAL
  591.         MSOUP3.IPOVAL=MPOVA3
  592.         LIPOV3(I1)=MPOVA3
  593.       ENDDO
  594.  
  595.  
  596. C==== CHARGEMENT DU TABLEAU DE TRAVAIL =================================
  597.  
  598. C     creation du tableau de travail pour ce MSOUPO
  599.       NCOU=NPOINT
  600.       NDDL=N*NC
  601.       SEGINI MTRAV3
  602.  
  603. C     boucle sur les pas de temps = boucle sur les chpoints
  604.       DO 720 I=1,IND1
  605. C     => pour ce MSOUPO, boucle sur les pas de temps = boucle sur les MPOVAL
  606.         MPOVAL=LIPOV(I)
  607. C       boucle sur les ddls = boucle sur sur les composantes + noeuds
  608.         IDDL=0
  609. C       boucle sur les composantes
  610.         DO 730 IC=1,NC
  611. C       boucle sur les noeuds
  612.         DO 731 IN=1,N
  613.             IDDL=IDDL+1
  614.             XXX(I,IDDL)= VPOCHA(IN,IC)
  615.  731    CONTINUE
  616.  730    CONTINUE
  617. c       IF(IIMPI.EQ.1)
  618. c      & write(ioimp,*) 'XXX(',I,')=',(XXX(I,iou),iou=1,IDDL)
  619.  720  CONTINUE
  620. C     On ne complete pas avec des 0 (inutile apres SEGINI MTRAV3)
  621.  
  622. C==== CALCUL DE LA FFT =================================================
  623.  
  624. c     boucle sur les ddls
  625.       IDDL=0
  626.       DO 740 IC=1,NC
  627.       DO 741 IN=1,N
  628.         IDDL=IDDL+1
  629.         IF(IIMPI.EQ.1) WRITE(IOIMP,*)' CALCUL DE LA FFT DU SIGNAL ',IDDL
  630.         CALL FFTL(XXX(1,IDDL),YYY(1),W(1),NPOINT)
  631.  741  CONTINUE
  632.  740  CONTINUE
  633.       IF(IIMPI.EQ.1) WRITE(IOIMP,*)' FFT CALCULEE '
  634. C
  635. C==== REMPLISSAGE DES LISTES DE CHPOINTS DE LA FFT =====================
  636. C
  637. C     creation de la liste des frequences
  638.       IF(ISOUPO.eq.1) THEN
  639.         JG=NNN
  640.         SEGINI MLREE1
  641.       ENDIF
  642.  
  643.       IF(IIMPI.EQ.1) WRITE(IOIMP,*)' SORTIE DE ',NNN,' POINTS '
  644.  
  645.       IF(ISORT.EQ.1) THEN
  646. C
  647. C     SORTIE EN PARTIE REELLE & IMAGINAIRE / FREQUENCE
  648. C
  649. c       Boucle sur les frequences = sur les chpoints
  650.         DO 750 I=KDEBU,KHALF
  651.           FREQ=REAL(I-1)*DFREQ
  652. c           IF(IIMPI.EQ.1) write(ioimp,*) I,' FREQ=',FREQ
  653. c           IF(IIMPI.EQ.1)
  654. c      &    write(ioimp,*)'XXX(',I,')=',(XXX(I,iou),iou=1,IDDL)
  655.           ITY=I-KDEBU+1
  656.           IF(ISOUPO.eq.1) MLREE1.PROG(ITY)=FREQ
  657.           MPOVA2=LIPOV2(ITY)
  658.           MPOVA3=LIPOV3(ITY)
  659. C         boucle sur les ddls = boucle sur sur les composantes + noeuds
  660.           IDDL=0
  661. C         boucle sur les composantes
  662.           DO 760 IC=1,NC
  663. C         boucle sur les noeuds
  664.           DO 761 IN=1,N
  665.               IDDL=IDDL+1
  666.               MPOVA2.VPOCHA(IN,IC)=DT*CXTORE(XXX(I,IDDL))
  667.               MPOVA3.VPOCHA(IN,IC)=DT*CXTOIM(XXX(I,IDDL))
  668.  761      CONTINUE
  669.  760      CONTINUE
  670.  
  671.  750    CONTINUE
  672.  
  673.       ELSE
  674. C
  675. C     SORTIE EN MODULE & PHASE / FREQUENCE
  676. C
  677.       ENDIF
  678.  
  679. C     SUPPRESSION
  680.       SEGSUP,MTRAV3
  681.       IF(ISOUPO.LT.NSOUPO) GOTO 710
  682. C     ici, on boucle sur les SOUPO...
  683.  
  684. C
  685. C==== DESACTIVATION ET ECRITURE DES LISTCHPOINT ET LISTREEL ============
  686. C
  687.       IF(IIMPI.EQ.1) WRITE(IOIMP,*)'NETTOYAGE'
  688.       SEGSUP MTRAV2
  689.       SEGSUP LIPOV,LIPOV2,LIPOV3
  690.       SEGDES MLCHPO,MLCHP2,MLCHP3,MLREE1
  691.       CALL ECROBJ('LISTREEL',MLREE1)
  692.       CALL ECROBJ('LISTCHPO',MLCHP3)
  693.       CALL ECROBJ('LISTCHPO',MLCHP2)
  694.  
  695.  777  CONTINUE
  696.       RETURN
  697.       END
  698.  
  699.  

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