Télécharger inter5.eso

Retour à la liste

inter5
C INTER5    SOURCE    PASCAL    21/02/24    21:15:10     10898          
      SUBROUTINE INTER5(TEMPS,KTE,KFT,FT0,ISPLIN,
     $     MLDERS,IHORS,IRET)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C   INTERPOLATION LINEAIRE ; recherche par binary search
C   Source : Donald Knuth, The Art of Computer Programming
C   Volume 3, Sorting and searching, Third edition, 1997
C   Addison-Wesley pp. 409-426
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  TEMPS  : valeur où interpoler la fonction
C  KTE    : abscisses de la fonction à interpoler
C  KFT    : ordonnées de la fonction à interpoler
C  LON    : longueur de ces deux précédentes listes
C  FT0    : valeur interpolée de la fonction à l'abscisse TEMPS
C  ISENS  : =0 KTE est croissante =1 KTE est décroissante
C  IHORS    =0 renvoie  une erreur si t est en-dehors des bornes
C         : =1 ne renvoie pas d'erreur si t est en-dehors des bornes
C              mais la valeur de la fonction pour t le plus proche
C         : =2 ne renvoie pas d'erreur si t est en-dehors des bornes
C              et on extrapole lineairement la valeur a partir des bords
C  IRET   : code retour =1 si succès, 0 sinon
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
-INC SMLREEL
      POINTEUR MLDERS.MLREEL
 
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
*
      FT0=0.D0
      MLREE1=KTE
      MLREE2=KFT
      LON=MLREE1.PROG(/1)
* Cas de bord (bourrin)
      TINF=MLREE1.PROG(1)
      TSUP=MLREE1.PROG(LON)
C      WRITE(IOIMP,*) 'IHORS=',IHORS
C      WRITE(IOIMP,*) 'ISENS=',ISENS
C      WRITE(IOIMP,*) 'TINF=',TINF
C      WRITE(IOIMP,*) 'TSUP=',TSUP
 
C---- TRAITEMENT CAS HORS INTERVALLE DE DEFINITION 
C
C Indicateur IHORS : 
C     IHORS = 0 : on doit etre dans l'intervalle de def., sinon erreur
C     IHORS = 1 : on renvoie la valeur aux bornes
C     IHORS = 2 : on extrapole lineairement en dehors des bornes
C
      IF (TEMPS.LE.TINF) THEN
         IF (IHORS.EQ.0) THEN
            GOTO 50
         ELSE IF (IHORS.EQ.1) THEN
            FT0 = MLREE2.PROG(1)
         ELSEIF (IHORS.EQ.2) THEN
            T1 = TINF
            T2 = MLREE1.PROG(2)
            F1 = MLREE2.PROG(1)
            F2 = MLREE2.PROG(2)
            DT = T2-T1
            TM = MAX(ABS(T1),ABS(T2))
            IF (ABS(DT).LE.(XPETIT+XZPREC*TM)) THEN
               CALL ERREUR(1124)
               RETURN
            ENDIF
            FT0 = (F2-F1)/DT*(TEMPS-T1)+F1
         ELSE
            CALL ERREUR(5)
            RETURN
         ENDIF
         GOTO 70
      ELSEIF (TEMPS.GE.TSUP) THEN
         IF (IHORS.EQ.0) THEN
            GOTO 50
         ELSE IF (IHORS.EQ.1) THEN
            FT0 = MLREE2.PROG(LON)
         ELSE IF (IHORS.EQ.2) THEN
            T1 = MLREE1.PROG(LON-1)
            T2 = TSUP
            F1 = MLREE2.PROG(LON-1)
            F2 = MLREE2.PROG(LON)
            DT = T2-T1
            TM = MAX(ABS(T1),ABS(T2))
            IF (ABS(DT).LE.(XPETIT+XZPREC*TM)) THEN
               CALL ERREUR(1124)
               RETURN
            ENDIF
            FT0 = (F2-F1)/DT*(TEMPS-T1)+F1
         ELSE
            CALL ERREUR(5)
            RETURN
         ENDIF
         GOTO 70
      ENDIF
 
C---- TRAITEMENT CAS DANS INTERVALLE DE DEFINITION 
C
* Recherche de l'index
      IDXINF=1
      IDXSUP=LON
 1    CONTINUE
C         WRITE(IOIMP,*) 'IDXINF=',IDXINF
C         WRITE(IOIMP,*) 'IDXSUP=',IDXSUP
      IF (IDXSUP.LT.IDXINF) GOTO 3
      IDXMIL=(IDXINF+IDXSUP)/2
      TEMMIL=MLREE1.PROG(IDXMIL)
C         WRITE(IOIMP,*) 'IDXMIL=',IDXMIL
C         WRITE(IOIMP,*) 'TEMMIL=',TEMMIL
      IF (TEMPS.LT.TEMMIL) THEN
         IDXSUP=IDXMIL-1
         GOTO 1
      ELSEIF (TEMPS.GT.TEMMIL) THEN
         IDXINF=IDXMIL+1
         GOTO 1
      ELSE
         FT0=MLREE2.PROG(IDXMIL)
         GOTO 70
      ENDIF
 
* Interpolation de la valeur :
 3    CONTINUE
      T1=MLREE1.PROG(IDXINF)
      T2=MLREE1.PROG(IDXSUP)
      DT=T2-T1
      XP2=(TEMPS-T1)/DT
      XP1=(T2-TEMPS)/DT
      FT1=MLREE2.PROG(IDXINF)
      FT2=MLREE2.PROG(IDXSUP)
      IF (ISPLIN.EQ.1) THEN
         DSFT1=MLDERS.PROG(IDXINF)
         DSFT2=MLDERS.PROG(IDXSUP)
         FT0=XP1*FT1+XP2*FT2+((XP1**3-XP1)*DSFT1+(XP2**3-XP2)*DSFT2)
     $        *(DT**2)/6.D0
      ELSE
         FT0=(XP1*FT1)+(XP2*FT2)
      ENDIF
 
* Sortie avec resultat
 70   CONTINUE
      IRET=1
      RETURN
 
* Sortie avec erreur si IHORS=0 et hors intervalle def. donnees
 50   CONTINUE
      CALL ERREUR(210)
      IRET=0
      RETURN
 
      END
 
 
 
 
 
 

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