Sources Esope | Cast3M

dertra
C DERTRA    SOURCE    PV090527  23/07/13    21:15:03     11708          
      SUBROUTINE DERTRA(NPTS,TRAC,X,Y,YPRIM,XINF,XSUP)
C
C--------------------------------------------------------------------
C Objet: calcul de la valeur d'une fonction F et de sa derivee
C        en un point X donne.Cette fonction est discretisee
C        en NPTS/2 points et tabulee dans un tableau de forme
C        suivante Trac(X1,F(X1),X2,F(X2),.........) avec X1&lt;X2&lt;...
C--------------------------------------------------------------------
C
C--------------------------------------------------------------------
C Entree: NPTS dimension du tableau TRAC
C         TRAC tableau de la fonction tabulee
C         X point ou sera calcule la fonction et sa derivee
C--------------------------------------------------------------------
C
C--------------------------------------------------------------------
C Sortie: Y valeur de la fonction en X
C         YPRIM valeur de la derivee en X
C         XINF,XSUP bornes entre lesquelles est compris X
C--------------------------------------------------------------------
      IMPLICIT INTEGER(I-N)
       IMPLICIT REAL*8 (A-H,O-Z)
       DIMENSION TRAC(*)
-INC CCREEL
C
       EPS1=1.D-14
       LIMIT=NPTS-1
C----------------------------------------------------------------
C la fonction est constante ou X est inferieur a la plus petite
C des abscisses ou X est superieur a la plus grande des abscisses
C----------------------------------------------------------------
       IF (NPTS.EQ.1) THEN
                                  Y=TRAC(1)
                                  YPRIM=0.D0
                                  XINF1=X
                                  XSUP1=X
                                  GOTO 250
       ENDIF
       IF (X.LT.TRAC(1)) THEN
                                  Y=TRAC(2)
                                  YPRIM=0.D0
                                  XINF1=X
                                  XSUP1=X
                                  GOTO 250
       ENDIF
       IF (X.GT.TRAC(LIMIT)) THEN
                                  Y=TRAC(NPTS)
                                  YPRIM=0.D0
                                  XINF1=X
                                  XSUP1=X
                                  GOTO 250
       ENDIF
C
C-------------------------------------------------------------------
C Recherche de l'intervalle [A,B] comprenant X
C-------------------------------------------------------------------
       EPS2=ABS(X)*EPS1
       XPLUS=X+EPS2
       XMOINS=X-EPS2
*  recherche par dichotomie
       NPT=NPTS/2
       IHAU=NPT
       IDEP=1
       IBAS=2*IDEP-IHAU
  15   CONTINUE
       IDEP=(IBAS+IHAU)/2
       IF (IDEP*2+1.GT.LIMIT) THEN
         I=LIMIT
         GOTO 10
       ENDIF
       IF (IDEP.LE.1) THEN
         I=1
         GOTO 10
       ENDIF
       IF (XMOINS.GT.TRAC(IDEP*2+1)) THEN
* on est trop petit
         IF (IBAS.EQ.IDEP) GOTO 17
         IBAS=IDEP
         GOTO 15
       ELSEIF(XPLUS.LT.TRAC(IDEP*2-1)) THEN
* on est trop grand
         IF (IHAU.EQ.IDEP) GOTO 17
         IHAU=IDEP
         GOTO 15
       ELSE
* on est bon
         I=IDEP*2-1
         GOTO 10
       ENDIF
17     write(6,*) ' erreur dans la recherche dichotomique '
       I=1
10     CONTINUE
       IF (XMOINS.GT.TRAC(I+2)) THEN
         I=I+2
         IF (I.EQ.LIMIT) THEN
          XINF1=TRAC(LIMIT-2)
          FXINF=TRAC(LIMIT-1)
          XSUP1=TRAC(LIMIT)
          FXSUP=TRAC(LIMIT+1)
          GOTO 100
         ENDIF
         GOTO 10
       ENDIF
       IF (XMOINS.GT.TRAC(I).AND.XPLUS.LT.TRAC(I+2)) THEN
           XINF1=TRAC(I)
           FXINF=TRAC(I+1)
           XSUP1=TRAC(I+2)
           FXSUP=TRAC(I+3)
           GOTO 100
       ENDIF
       IF (XMOINS.LE.TRAC(I)) THEN
           IF ( I.EQ.1 ) THEN
               XINF1=TRAC(1)
               FXINF=TRAC(2)
               XSUP1=TRAC(3)
               FXSUP=TRAC(4)
               GOTO 100
           ELSE
               XINF1=TRAC(I-2)
               FXINF=TRAC(I-1)
               XSUP1=TRAC(I+2)
               FXSUP=TRAC(I+3)
               Y=TRAC(I+1)
              XDEN=sign(max(xpetit/xzprec,abs(xsup1-xinf1)),xsup1-xinf1)
               YPRIM=(FXSUP-FXINF)/XDEN
               GOTO 250
           ENDIF
       ENDIF
       IF (XPLUS.GE.TRAC(I+2)) THEN
             IF ( I+2.EQ.LIMIT ) THEN
               XINF1=TRAC(LIMIT-2)
               FXINF=TRAC(LIMIT-1)
               XSUP1=TRAC(LIMIT)
               FXSUP=TRAC(LIMIT+1)
               GOTO 100
             ELSE
               XINF1=TRAC(I)
               FXINF=TRAC(I+1)
               XSUP1=TRAC(I+4)
               FXSUP=TRAC(I+5)
               Y=TRAC(I+3)
               YPRIM=(FXSUP-FXINF)/(XSUP1-XINF1)
               GOTO 250
           ENDIF
       ENDIF
       write (6,*) ' on ne devrait pas passer la (dans dertra) '
       write (6,*) ' X LIMIT I ',x,limit,i
       STOP 12
C
C------------------------------------------------------------------
C Calcul de la derivee de F et de la valeur de F en X
C------------------------------------------------------------------
100    CONTINUE                                                    
       XDEN=sign(max(xpetit/xzprec,abs(xsup1-xinf1)),xsup1-xinf1)
       YPRIM=(FXSUP-FXINF)/XDEN
       Y=FXINF+YPRIM*(X-XINF1)
250    CONTINUE
       XSUP=XSUP1
       XINF=XINF1
C
       RETURN
       END