C ONDE SOURCE OF166741 25/02/20 21:17:13 12165 C SUBROUTINE ONDE IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-V) IMPLICIT COMPLEX*16 (W-Z) CHARACTER *72 TI CHARACTER *4 NSORT(4),MCLE(5) LOGICAL INV C C======================================================================= C = CALCUL DE LA TRANSFORMEE EN ONDELETTE CONTINUE D'UN SIGNAL = C = = C = SYNTAXE : = C = = C = TRANSFORMEE DE FOURIER DIRECTE = C = = C = EVO = 'ONDE' EVO1 EXP2 SORT 'FMIN' V1 'FMAX' V2 ( 'NFRQ' V3 ) = C = ( 'PULS' V4 ) ( 'EPSI' V5 ) = C = OU = C = CHP1 MAIL1 = 'ONDE' EVO1 EXP2 SORT 'FMIN' V1 'FMAX' V2 = C = ( 'NFRQ' V3 ) ( 'PULS' V4 ) = 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 = SORT : TYPE DE SORTIE ; = C = = 'CRMO' EXTRACTION CRETE SUR CRITERE MODULE = C = ( PAR DEFAUT ) = C = = 'CRPH' EXTRACTION CRETE SUR CRITERE PHASE = C = = 'REIM' PART REEL & PART IMAG = C = = 'MOPH' MODULE & PHASE = C = V1 : FREQUENCE MINI A VISUALISER = C = V2 : FREQUENCE MAXI A VISUALISER = C = V3 : NOMBRE DE PAS EN FREQUENCE = C = V4 : PULSATION ONDELETTE MERE = C = V5 : CRITERE DE NULLITE = C = COU1 : COULEUR ATTRIBUEE A LA PREMIERE COURBE (FACULTATIF) = C = COU2 : COULEUR ATTRIBUEE A LA DEUXIEME COURBE (FACULTATIF) = C = = C = = C======================================================================= C -INC CCREEL -INC CCGEOME -INC PPARAM -INC CCOPTIO -INC SMEVOLL -INC SMLREEL -INC SMCHPOI -INC SMELEME -INC SMCOORD C SEGMENT MTRAV1 IMPLIED XX(NCOU),YY(NCOU),WOND(NCOU),WONC(NCOU) IMPLIED YY1(NCOU),YY2(NCOU),VM(NCOU),VP(NCOU) IMPLIED VDP(NCOU),VPL(NCOU) ENDSEGMENT C SEGMENT MTRAV2 IMPLIED W(NEXP) ENDSEGMENT SEGMENT MTRAV3 IMPLIED VML((NA+1)*NCOU),VML2((NA+1)*NCOU) IMPLIED ICP1((NA+1)*NCOU) ENDSEGMENT C DATA NSORT(1),NSORT(2),NSORT(3)/'CRMO','CRPH','REIM'/ DATA NSORT(4)/'MOPH'/ DATA MCLE(1),MCLE(2),MCLE(3),MCLE(4)/'FMIN','FMAX','NFRQ','EPSI'/ DATA MCLE(5)/'PULS'/ 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) IF(IRET2.EQ.0)GOTO 666 C C LECTURE DU TYPE DE SORTIE C CALL LIRMOT(NSORT,4,ISORT,0) IF(ISORT.EQ.0) ISORT=1 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)(1:4) CALL ERREUR(166) GOTO 666 ENDIF IF(FRMI.LT.(-1.D-20)) FRMI=-FRMI IF(ABS(FRMI).LT.1.D-20) FRMA=-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)(1:4) 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 DU NOMBRE DE DIVISION EN FREQUENCE C CALL LIRMOT(MCLE(3),1,IRET,0) IF(IRET.EQ.1) THEN CALL LIRENT(NA,0,IRET1) IF(IRET1.EQ.0) THEN MOTERR(1:4)=MCLE(3)(1:4) CALL ERREUR(166) GOTO 666 ENDIF ELSE NA=100 ENDIF C C LECTURE DU CRITERE C CALL LIRMOT(MCLE(4),1,IRET,0) IF(IRET.EQ.1) THEN CALL LIRREE(EPSI,0,IRET1) IF(IRET1.EQ.0) THEN MOTERR(1:4)=MCLE(4)(1:4) CALL ERREUR(166) GOTO 666 ENDIF ELSE EPSI=1.D-4 ENDIF C C LECTURE DE LA FREQUENCE ONDELETTE C CALL LIRMOT(MCLE(5),1,IRET,0) IF(IRET.EQ.1) THEN CALL LIRREE(V0,0,IRET1) IF(IRET1.EQ.0) THEN MOTERR(1:4)=MCLE(5)(1:4) CALL ERREUR(166) GOTO 666 ENDIF ELSE V0=5.D0 ENDIF C C LECTURE DE LA COULEUR C CALL LIRMOT(NCOUL,NBCOUL,ICOU1,0) IF(ICOU1.EQ.0) ICOU1=IDCOUL+1 ICOU1=ICOU1-1 CALL LIRMOT(NCOUL,NBCOUL,ICOU2,0) IF(ICOU2.EQ.0) ICOU2=IDCOUL+1 ICOU2=ICOU2-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 C SEGACT MLREE1,MLREE2 L1=MLREE1.PROG(/1) DT=MLREE1.PROG(2)-MLREE1.PROG(1) TDEB=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 TFIN=TDEB+(NPOINT*DT) C NCOU=NPOINT SEGINI MTRAV1 C NEXP=NPOINT/2 SEGINI MTRAV2 SEGINI MTRAV3 C C CHARGEMENT DES TABLEAUX DE TRAVAIL C IND1=L1 IF(NPOINT.LT.L1)IND1 = NPOINT DO 1100 I=1,IND1 XX(I)=MLREE2.PROG(I) C IF(IIMPI.EQ.1)WRITE(IOIMP,*) ' XX(',I,') = ',XX(I) 1100 CONTINUE IF(NPOINT.GT.L1) THEN L2=L1+1 DO 1110 I=L2,NPOINT XX(I)=0.D0 1110 CONTINUE ENDIF DUREE=DT*DBLE(NPOINT) DFREQ=1.D0/DUREE KHALF=INT(NPOINT/2)+1 KDEBU=1 NNN=KHALF-KDEBU+1 JG=NNN C C CALCUL DE LA FFT C IF(IIMPI.EQ.1) WRITE(IOIMP,*)' CALCUL DE LA FFT DU SIGNAL ' C C================= CALCUL TFR FONCTION F(T) ============= C C INV=.FALSE. CALL WEXP(INV,NPOINT,W(1)) CALL FFTL(XX(1),YY(1),W(1),NPOINT) SEGDES MLREE2 C C======================================================== C AMIN=V0/2.D0/XPI/FRMA AMAX=V0/2.D0/XPI/FRMI NPO=NPOINT*(NA+1) GOTO(1999,2999,4199,4199),ISORT RETURN C C CRITERE SUR LE MODULE C 1999 SEGINI MLREE1 SEGINI MLREE2 SEGINI MLREE3 JG=0 SEGADJ MLREE1 K=1 C================ BOUCLE SUR LES FREQUENCES ================= C DO 2000 I=0,NA A=AMIN+((AMAX-AMIN)*I/NA) V=V0/A C PRODUIT ONDELETTE TFR SIGNAL DO 1200 J=1,KHALF VJ=DBLE(J-1)*DFREQ*2.D0*XPI*A V2=(VJ-V0)**2.D0 PM=EXP(-0.5D0*V2)*((2.D0*XPI)**0.5D0) WOND(J)=PM*XX(J) 1200 CONTINUE J=NPOINT+1 KK=1 L11=KHALF-1 DO 1205 II=2,L11 KK=KK+1 J=J-1 VJ=DBLE(J-1-NPOINT)*DFREQ*2.D0*XPI*A V2=(VJ-V0)**2.D0 PM=EXP(-0.5D0*V2)*((2.D0*XPI)**0.5D0) WOND(J)=CONJG(XX(KK))*PM 1205 CONTINUE C C================= CALCUL TFR INVERSE ============= C C INV=.TRUE. CALL WEXP(INV,NPOINT,W(1)) CALL FFTL(WOND(1),YY1(1),W(1),NPOINT) C C======================================================== DO 1210 J=1,NPOINT R1=WOND(J) R2=CXTOIM(WOND(J))*(-1.D0) RT=SQRT(R1*R1+R2*R2) VML((I*NPOINT)+J)=RT/(DBLE(NPOINT)) 1210 CONTINUE 2000 CONTINUE C RECHERCHE CRETE JN=INT(((2.D0*AMAX)-TDEB)/DT-1.D0) JX=INT((TFIN-(2.D0*AMAX)-TDEB)/DT-1.D0) DO 2500 J=JN,JX VP=-1.D0 DO 2100 I=0,NA A=AMIN+((AMAX-AMIN)*I/NA) V=V0/A TCOU=TDEB+(DT*DBLE(J-1)) VMIJ=VML((I*NPOINT)+J)/((XPI/2.D0)**0.5D0) IF(VMIJ.GT.VP) THEN VP=VMIJ VF=V/2.D0/XPI ENDIF 2100 CONTINUE IF(VP.GT.0) THEN JG=K SEGADJ MLREE1 SEGADJ MLREE2 SEGADJ MLREE3 MLREE3.PROG(K)=VP MLREE2.PROG(K)=VF MLREE1.PROG(K)=TCOU K=K+1 ENDIF 2500 CONTINUE GOTO 4001 C C CRITERE SUR LA PHASE C 2999 SEGINI MLREE1 SEGINI MLREE2 SEGINI MLREE3 JG=0 SEGADJ MLREE1 K=1 C C================ BOUCLE SUR LES FREQUENCES ================= C DO 4000 I=0,NA A=AMIN+((AMAX-AMIN)*I/NA) V=V0/A C PRODUIT ONDELETTE TFR SIGNAL DO 3200 J=1,KHALF VJ=DBLE(J-1)*DFREQ*2.D0*XPI*A V2=(VJ-V0)**2.D0 PM=EXP(-0.5D0*V2)*((2.D0*XPI)**0.5D0) WOND(J)=PM*XX(J) 3200 CONTINUE J=NPOINT+1 KK=1 L11=KHALF-1 DO 3205 II=2,L11 KK=KK+1 J=J-1 VJ=DBLE(J-1-NPOINT)*DFREQ*2.D0*XPI*A V2=(VJ-V0)**2.D0 PM=EXP(-0.5D0*V2)*((2.D0*XPI)**0.5D0) WOND(J)=CONJG(XX(KK))*PM 3205 CONTINUE C C================= CALCUL TFR INVERSE ============= C C INV=.TRUE. CALL WEXP(INV,NPOINT,W(1)) CALL FFTL(WOND(1),YY1(1),W(1),NPOINT) C C======================================================== DO 3210 J=1,NPOINT R1=WOND(J) R2=CXTOIM(WOND(J))*(-1.D0) RT=SQRT(R1*R1+R2*R2) VM(J)=RT/(DBLE(NPOINT)) CALL FACOMP(R1,R2,RT) VP(J)=ATAN2(R2,R1) VPL(J)=VP(J)+0.D0 3210 CONTINUE C================= CALCUL DERIVEE PHASE ============= VDP(1)=(VPL(2)-VPL(1))/DT DO 3220 J=2,(NPOINT-1) VDP(J)=(VPL(J+1)-VPL(J-1))/2.D0/DT 3220 CONTINUE VDP(NPOINT)=(VPL(NPOINT)-VPL(NPOINT-1))/DT C CORRECTION DO 3230 J=1,NPOINT DPI=XPI*0.5D0/DT SPI=DPI*(-1.D0) IF(VDP(J).GT.DPI) VDP(J)=VDP(J)-(XPI/DT) IF(VDP(J).LT.SPI) VDP(J)=VDP(J)+(XPI/DT) VDP(J)=VDP(J)+(V0/A) C RECHERCHE CRETE AVDP=ABS(VDP(J))-(EPSI*V0/A) IF(AVDP.LT.0) THEN FCOU=0.5D0/A TCOU=TDEB+(DT*DBLE(J-1)) TFCO=TFIN-(1.D0/FCOU) TDCO=TDEB+(1.D0/FCOU) IF(TCOU.GT.TDCO) THEN IF(TCOU.LT.TFCO) THEN VMJ=VM(J)/((XPI/2.D0)**0.5D0) JG=K SEGADJ MLREE1 SEGADJ MLREE2 SEGADJ MLREE3 MLREE3.PROG(K)=VMJ MLREE2.PROG(K)=V/2.D0/XPI MLREE1.PROG(K)=TCOU K=K+1 ENDIF ENDIF ENDIF 3230 CONTINUE 4000 CONTINUE C----------- SORTIES ------------ 4001 SEGDES MLREE1,MLREE2,MLREE3 C C CREATION DE L'OBJET EVOLUTIO CRETE C N=2 SEGINI MEVOLL IPVO=MEVOLL TI(1:72)=TITREE IEVTEX=TI C C PREMIERE COURBE C SEGINI KEVOLL TYPX='LISTREEL' TYPY='LISTREEL' c KEVTEX=TI KEVTEX='Freq (Hz)' IEVOLL(1)=KEVOLL C IPROGX=MLREE1 NOMEVX='TEMPS S' C IPROGY=MLREE2 NOMEVY='FREQ HZ' C NUMEVX=ICOU1 NUMEVY=' ' SEGDES KEVOLL C C DEUXIEME COURBE C SEGINI KEVOLL TYPX='LISTREEL' TYPY='LISTREEL' c KEVTEX=TI KEVTEX='Amp' IEVOLL(2)=KEVOLL C IPROGX=MLREE1 NOMEVX='TEMPS S' C IPROGY=MLREE3 NOMEVY='MODULE ' C NUMEVX=ICOU2 NUMEVY=' ' C SEGDES KEVOLL C SEGSUP MTRAV1,MTRAV2,MTRAV3 SEGDES KEVOL1,MEVOL1 C SEGDES MEVOLL CALL ECROBJ('EVOLUTIO',IPVO) RETURN C=============== CHAMPOINT =========================== 4199 IF (IDIM.EQ.0) RETURN NBSOUS=0 NBREF=0 NBELEM=NA*(NPOINT-1) NBNN=4 SEGINI MELEME ITYPEL=8 segact mcoord*mod IDEB=NBPTS C C================ BOUCLE SUR LES FREQUENCES ================= C DO 5000 I=0,NA A=AMIN+((AMAX-AMIN)*I/NA) V=V0/A C PRODUIT ONDELETTE TFR SIGNAL DO 4200 J=1,KHALF VJ=DBLE(J-1)*DFREQ*2.D0*XPI*A V2=(VJ-V0)**2.D0 PM=EXP(-0.5D0*V2)*((2.D0*XPI)**0.5D0) WOND(J)=PM*XX(J) 4200 CONTINUE J=NPOINT+1 KK=1 L11=KHALF-1 DO 4205 II=2,L11 KK=KK+1 J=J-1 VJ=DBLE(J-1-NPOINT)*DFREQ*2.D0*XPI*A V2=(VJ-V0)**2.D0 PM=EXP(-0.5D0*V2)*((2.D0*XPI)**0.5D0) WOND(J)=CONJG(XX(KK))*PM 4205 CONTINUE C C================= CALCUL TFR INVERSE ============= C C INV=.TRUE. CALL WEXP(INV,NPOINT,W(1)) CALL FFTL(WOND(1),YY1(1),W(1),NPOINT) C C======================================================== DO 4210 J=1,NPOINT II=(I*NPOINT)+J IF(ISORT.EQ.3) THEN VML(II)=WOND(J)/(CMPLX(NPOINT))/((XPI/CMPLX(2.D0))** & CMPLX(0.5D0)) VML2(II)=CXTOIM(WOND(J))/(DBLE(NPOINT))/((XPI/2.D0)** & 0.5D0) ENDIF IF(ISORT.EQ.4) THEN R1=WOND(J) R2=CXTOIM(WOND(J)) VML(II)=SQRT(R1*R1+R2*R2)/(DBLE(NPOINT))/((XPI/2.D0)** & 0.5D0) VML2(II)=ATAN2(R2,R1)*360.D0/XPI ENDIF FCOU=0.5D0/A TCOU=TDEB+(DT*DBLE(J-1)) TFCO=TFIN-(1.D0/FCOU) TDCO=TDEB+(1.D0/FCOU) IF ((TCOU.LT.TDCO).OR.(TCOU.GT.TFCO)) THEN VML(II)=0.D0 VML2(II)=0.D0 ENDIF NBPTS=NBPTS+1 SEGADJ MCOORD XCOOR((NBPTS-1)*(IDIM+1)+1)=DBLE(J-1)/(NPOINT-1) XCOOR((NBPTS-1)*(IDIM+1)+2)=1.D0-(DBLE(I)/NA) IF(IDIM.NE.2) THEN XCOOR((NBPTS-1)*(IDIM+1)+3)=0.D0 ENDIF XCOOR(NBPTS*(IDIM+1))=0.D0 4210 CONTINUE 5000 CONTINUE DO 5150 I=0,(NA-1) DO 5100 J=1,(NPOINT-1) INE=I*(NPOINT-1)+J ICOLOR(INE)=IDCOUL NUM(1,INE)=IDEB+J+(I*NPOINT) NUM(2,INE)=IDEB+(J+1)+(I*NPOINT) NUM(3,INE)=IDEB+(J+1)+((I+1)*NPOINT) NUM(4,INE)=IDEB+J+((I+1)*NPOINT) 5100 CONTINUE 5150 CONTINUE C--- TABLEAU CONVERSION QUA4->POI1 ---- DO 5151 I=1,NPO ICP1(I)=0 5151 CONTINUE ICON=0 DO 5155 I=1,4 DO 51551 J=1,NBELEM IKI=NUM(I,J)-IDEB IF (ICP1(IKI).NE.0) GOTO 51551 ICON=ICON+1 ICP1(IKI)=ICON 51551 CONTINUE 5155 CONTINUE C GOTO 4310 SEGDES MELEME CALL ECROBJ('MAILLAGE',MELEME) SEGACT MELEME CALL CHANGE (MELEME,1) NBPOIN=NUM(/2) NSOUPO=1 NC=2 N=NBPOIN NAT=1 SEGINI MCHPOI MOCHDE=' chpoint de coordonnees ' MTYPOI=' ' JATTRI(1) = 1 IPPOI=MCHPOI SEGINI MSOUPO SEGINI MPOVAL IPCHP(1)=MSOUPO IFOPOI = IFOUR SEGDES MCHPOI IF(ISORT.EQ.3) THEN NOCOMP(1)='PREE' NOCOMP(2)='PIMA' ENDIF IF(ISORT.EQ.4) THEN NOCOMP(1)='MODU' NOCOMP(2)='PHAS' ENDIF NOHARM(1)=NIFOUR NOHARM(2)=NIFOUR IGEOC=MELEME IPOVAL=MPOVAL SEGDES MSOUPO DO 5160 I= 1 ,NPO II=ICP1(I) VPOCHA(II,1)=VML(I) VPOCHA(II,2)=VML2(I) 5160 CONTINUE SEGDES MPOVAL SEGDES MELEME CALL ECROBJ('CHPOINT ',IPPOI) SEGSUP MTRAV1,MTRAV2,MTRAV3 666 CONTINUE RETURN END