Esope source files

   1 : C LAPLIN    SOURCE    CHAT      05/01/13    01:12:59     5004           
   2 :       SUBROUTINE LAPLIN
   3 : C     
   4 :       IMPLICIT INTEGER(I-N)
   5 :       IMPLICIT REAL*8     (A-H,O-Y)
   6 :       IMPLICIT COMPLEX*16 (Z)
   7 : C     
   8 : C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
   9 : C     
  10 : C   OPERATEUR  LAPL
  11 : C     
  12 : C   CALCUL LA TRANSFORMEE DE LAPLACE INVERSE
  13 : C     
  14 : C     
  15 : C   SORTIES :
  16 : C     
  17 : C     
  18 : C   AUTEUR : SAINT-DIZIER ET DE LANGRE
  19 : C   DATE   : 05 DECEMBRE 1989
  20 : C     
  21 : C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  22 : C     
  23 :       LOGICAL INV
  24 : 
  25 : -INC CCREEL
  26 : -INC SMLREEL
  27 : -INC SMEVOLL
  28 : -INC CCOPTIO
  29 : C     
  30 :       SEGMENT MTRAV
  31 :       COMPLEX*16 ZW(NEXP)
  32 :       COMPLEX*16 ZXX(NPT)
  33 :       COMPLEX*16 ZYY(NPT)
  34 :       ENDSEGMENT
  35 : C     
  35 :       CALL LIROBJ('LISTREEL',MLREE1,1,IRETOU)
  36 :       IF (IERR.NE.0) RETURN
  36 :       CALL LIROBJ('LISTREEL',MLREE2,1,IRETOU)
  37 :       IF (IERR.NE.0) RETURN
  37 :       CALL LIROBJ('LISTREEL',MLREE3,1,IRETOU)
  38 :       IF (IERR.NE.0) RETURN
  38 :       CALL LIRREE(AA,1,IRETOU)
  39 :       IF (IERR.NE.0) RETURN
  39 :       CALL LIRENT(LL,1,IRETOU)
  40 :       IF (IERR.NE.0) RETURN
  41 : C     
  42 :       SEGACT MLREE1
  43 :       SEGACT MLREE2
  44 :       SEGACT MLREE3
  45 : C     
  46 :       NSUM = MLREE1.PROG(/1)
  47 :       NN = NSUM / LL
  48 :       JG = NN
  49 :       POINTEUR MLREE4.MLREEL
  50 :       POINTEUR MLREE5.MLREEL
  51 :       SEGINI MLREE4
  52 :       SEGINI MLREE5
  53 : C     
  54 :       NEXP = NN/2
  55 :       NPT  = NN
  56 :       SEGINI MTRAV
  57 : C     
  58 :       N  = 1
  59 :       SEGINI MEVOL1
  60 :       SEGINI KEVOL1
  61 : C     
  62 : C -- DETERMINATION DES CARACTERISTIQUES DU CALCUL
  63 : C     
  64 :       ZI = (0.D0,1.D0)
  65 :       T = 2.D0 * XPI / (MLREE1.PROG (2))
  66 :       DELTAT = T / NN
  67 : C     
  68 : C -- ALGORITHME DE DURBIN AVEC FFT
  69 : C     
  70 :       DO 20 K = 0 , NN - 1
  71 :          XA = 0D0
  72 :          XB = 0D0
  73 :          DO 21 KL= 0 , LL -1
  74 :             IKL =K + KL * NN
  75 :             XA = XA + MLREE2.PROG(IKL+1)
  76 :             XB = XB + MLREE3.PROG(IKL+1)
  77 :  21      CONTINUE
  78 :          ZXX(K + 1) = XA + ZI* XB
  79 :  20   CONTINUE
  80 : C     
  81 :       INV = .TRUE.
  81 :       CALL WEXP(INV,NPT,ZW)
  81 :       CALL FFTL(ZXX,ZYY,ZW,NPT)
  82 : C     
  83 : C     
  84 :       DO 30 J = 0 , NN-1
  85 :          TJ = J*DELTAT
  86 :          MLREE4.PROG(J+1)= TJ
  87 :          CJ = EXP(AA*TJ)
  88 :          FF =CXTORE(ZXX(J+1))
  89 :          MLREE5.PROG(J+1) = (2.D0*CJ/T)*(-.5D0*MLREE2.PROG(1)+FF)
  90 : C           MLREE5.PROG(J+1) = (2.D0*CJ/T)*(                   FF)
  91 : C           MLREE5.PROG(J+1) =                               FF
  92 :  30   CONTINUE
  93 : C     
  94 :       MEVOL1.ITYEVO = 'REEL'
  95 :       MEVOL1.IEVTEX = 'TRANSFORMEE DE LAPLACE'
  96 :       MEVOL1.IEVOLL(1) = KEVOL1
  97 : C     
  98 :       KEVOL1.IPROGX = MLREE4
  99 :       KEVOL1.IPROGY = MLREE5
 100 :       KEVOL1.NUMEVX = 0
 101 :       KEVOL1.TYPX   = 'LISTREEL'
 102 :       KEVOL1.TYPY   = 'LISTREEL'
 103 :       KEVOL1.NOMEVX = 'TEMPS'
 104 :       KEVOL1.NOMEVY = 'F(T)'
 105 :       KEVOL1.KEVTEX = 'LAPLACE INVERSE'
 106 : C     
 106 :       CALL ECROBJ('EVOLUTION',MEVOL1)
 107 : C     
 108 :       SEGDES MLREE1
 109 :       SEGDES MLREE2
 110 :       SEGDES MLREE3
 111 :       SEGDES MLREE4
 112 :       SEGDES MLREE5
 113 :       SEGDES MEVOL1
 114 :       SEGDES KEVOL1
 115 :       END
 116 : 
 117 : 
 118 :  
 119 :  
 120 :