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

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