C INTERP    SOURCE    GOUNAND   12/09/11    21:15:03     7496           
      SUBROUTINE INTERP
C
C=======================================================================
C
C     Opérateur IPOL
C
C----------------------------
C Phrases d'appel (GIBIANE) :
C----------------------------
C
C         OBJET2 = 'IPOL'  OBJET1  | LISTREE1  LISTREE2  | ;
C                                  | EVOL1               |
C
C         OBJET3 = 'IPOL'  TAB1    FLOT1 ;
C
C
C         OBJET2 = 'IPOL'  OBJET1  NUA1  | ('GAUSS') | ;
C                                        |  'RATIO'  |
C
C------------------------
C Opérandes et résultat :
C------------------------
C
C     Syntaxe 1
C
C       LISTREE1 : liste des abscisses T    de la fonction
C
C       LISTREE2 : liste des ordonnées F(T) de la fonction
C
C       EVOL1    : Evolution représentant la fonction
C
C       OBJET1   : objet contenant les valeurs à interpoler sur F
C                  (type FLOTTANT, CHPOINT, LISTREEL, MCHAML)
C
C       OBJET2   : valeurs de la fonction interpolée aux points donnés
C                  (même type que OBJET1)
C
C     Syntaxe 2
C
C       TAB1     : table de sous-type 'RESULTAT'
C
C       FLOT1    : valeur à interpoler sur F
C
C       OBJET3   : valeurs de la fonction interpolée au point donné
C                  (type CHPOINT ou MCHAML)
C
C     Syntaxe 3
C
C       NUA1     : objet nuage représentant la fonction à interpoler
C
C=======================================================================
C
C    Remarques
C
C       Les listes LISTREE1 et LISTREE2 doivent se correspondre
C
C       L'évolution EVOL1 doit être élémentaire
C
C       ATTENTION : la liste des abscisses donnéee est supposée monotone
C
C=======================================================================
C
C    Auteur   : 
C    Création : 
C    Modifications : 
C       PM le 03/10/2006 : lecture d'une évolution en entrée
C             19/03/2007 : vérification que les listes ont une longueur > 1
C       SG    22/10/2009 : option SPLINE (biblio : numerical recipes)
C=======================================================================

      IMPLICIT INTEGER(I-N)
      IMPLICIT    REAL*8(A-H,O-Z)
      CHARACTER*8 TYPE1,TYPE2,TYPEI,TYPO1,TYPO2,TYPOBJ
      CHARACTER*8 CHARIN,CHARRE
      LOGICAL     LOGIN,LOGRE
      CHARACTER*4 CLE(1)
      DATA CLE/'SPLI'/
      CHARACTER*4 CDER(2)
      DATA CDER/'DGAU','DDRO'/
      LOGICAL LCDER(2)
      REAL*8 XCDER(2)
      INTEGER OOOVAL

-INC SMLREEL
      POINTEUR MLREE4.MLREEL,MLDERS.MLREEL,MLDER.MLREEL
-INC CCREEL
-INC SMCHPOI
-INC SMELEME
-INC CCOPTIO
-INC SMCHAML
-INC SMTABLE
-INC SMEVOLL
      SEGMENT TABR
         REAL*8   TEMA(LOR)
      ENDSEGMENT

*     syntaxe 3 ?
      CALL LIROBJ('NUAGE',INUA,0,IRETOU)
      IF(IRETOU.EQ.1) THEN
        CALL IPLNU1(INUA)
        RETURN
      ENDIF
*     syntaxe 2 ?
      CALL LIROBJ('TABLE',IPOT,0,IRETOU)
      IF(IRETOU.EQ.1) THEN
        CALL  LIRREE(TEMPS,1,IRETOU)
        IF(IERR.NE.0) RETURN
        GOTO 50
      ENDIF

*     syntaxe 1 (INDIC=1 à 4)
      CALL LIRREE(TEMPS,0,IRETOU)
      IF(IRETOU.EQ.1) THEN
        INDIC=1
      ELSE
        CALL LIROBJ('CHPOINT ',MTEMP,0,IRETOU)
        IF(IRETOU.EQ.1) THEN
          INDIC=2
        ELSE
          CALL LIROBJ('MCHAML',IPO1,0,IRETOU)
          IF(IRETOU.EQ.1) THEN
            INDIC=4
          ELSE
            CALL LIROBJ('LISTREEL',MLIST,0,IRETOU)
            IF(IRETOU.EQ.0) THEN
*             Pas d opérande correcte trouvée
              CALL QUETYP(MOTERR(1:8),0,IRETOU)
              IF(IRETOU.NE.0) THEN
*               On ne veut pas d'objet de type %m1:8
                CALL ERREUR(39)
              ELSE
*               Cet opérateur a encore besoin d'un opérande
                CALL ERREUR(533)
              ENDIF
              RETURN
            ENDIF
            INDIC=3
          ENDIF
        ENDIF
      ENDIF

*     Lecture de la fonction à interpoler
      CALL LIROBJ('EVOLUTIO',IEVL,0,IRETOU)
      IF(IRETOU.EQ.1) THEN
        MEVOLL=IEVL
        SEGACT, MEVOLL
*       Vérification que l'évolution est élémentaire
        IF (IEVOLL(/1).NE.1) THEN
*         Opération interdite sur un objet complexe
          CALL ERREUR(25)
          RETURN 
        ENDIF
        KEVOLL=IEVOLL(1)
        SEGACT, KEVOLL
*       de type scalaire
        IF (ITYEVO.NE.'REEL') THEN
*         Opération interdite sur un objet complexe
          CALL ERREUR(25)
          RETURN 
        ENDIF
*       peuplée de flottants
        IF ((TYPX.NE.'LISTREEL').OR.(TYPY.NE.'LISTREEL')) THEN
*         Certaines courbes ne sont pas du bon type
          CALL ERREUR(871)
          RETURN 
        ENDIF
        KTE=IPROGX
        KFT=IPROGY
        SEGDES,KEVOLL,MEVOLL
      ELSE
        CALL LIROBJ('LISTREEL',KTE,1,IRETOU)
        IF(IERR.NE.0) RETURN
        CALL LIROBJ('LISTREEL',KFT,1,IRETOU)
        IF(IERR.NE.0) RETURN
      ENDIF

*     Longueur des listes décrivant la fonction
      MLREE1=KTE
      MLREE2=KFT
      SEGACT MLREE1,MLREE2
      IF(MLREE1.PROG(/1).NE.MLREE2.PROG(/1)) THEN
*       Les suites n ont pas la même longueur
        CALL ERREUR(212)
        RETURN
      ENDIF
      LON=MLREE1.PROG(/1)
C 897 2
C La dimension des LISTREEL doit etre plus grande que 1
      IF (LON.LT.1) THEN
         CALL ERREUR(897)
         RETURN 
      ENDIF 
      TDEB=MLREE1.PROG(1)
      TFIN=MLREE1.PROG(LON)
*     Les x sont-ils croissants ou décroissants ?
      IF (TFIN.GE.TDEB) THEN
         ISENS=0
      ELSE 
         ISENS=1
         JG=LON
         SEGINI,MLREE3
         SEGINI,MLREE4
         DO ILON=1,LON
            MLREE3.PROG(ILON)=MLREE1.PROG(LON-ILON+1)
            MLREE4.PROG(ILON)=MLREE2.PROG(LON-ILON+1)
         ENDDO
         SEGDES,MLREE1
         SEGDES,MLREE2
         MLREE1=MLREE3
         MLREE2=MLREE4
         TDEB=MLREE1.PROG(1)
         TFIN=MLREE1.PROG(LON)
      ENDIF
* On va vérifier que la liste est ordonnée
      IVERI=1
C
C Vérification que la liste est ordonnée
C      
      IF (IVERI.EQ.1) THEN
         TPRE=TDEB
         DO ILON=2,LON
            TCOU=MLREE1.PROG(ILON)
            IF (TCOU.LT.TPRE) GOTO 6661
            TPRE=TCOU
         ENDDO
      ENDIF
      GOTO 1
 6661 CONTINUE
C     "La suite de réels doit être croissante"
C  "monotone" !
      CALL ERREUR(249)
      RETURN
 1    CONTINUE
*
* Y a-t-il une option SPLINE ?
*
      CALL LIRMOT(CLE,1,ISPLIN,0)
      IF (ISPLIN.EQ.1) THEN
         LCDER(1)=.FALSE.
         LCDER(2)=.FALSE.
*  Lecture des mots clés et valeurs associées
*  On lit les valeurs des dérivées premières à gauche et à droite
*  Si elles ne sont pas données, c'est la condition à la limite
*  naturelle qui s'applique
 77      CONTINUE
         CALL LIRMOT(CDER,2,ICDER,0)
         IF (ICDER.GT.0) THEN
            LCDER(ICDER)=.TRUE.
            CALL LIRREE(XCDER(ICDER),1,IRETOU)
            IF(IERR.NE.0) RETURN
            GOTO 77
         ENDIF
         JG=LON
         SEGINI MLDERS
         SEGINI MLDER
         IF (LCDER(1)) THEN
* Cas où on prescrit la dérivée première à gauche
            MLDERS.PROG(1)=-0.5D0
            DX=MLREE1.PROG(2)-MLREE1.PROG(1)
            DY=MLREE2.PROG(2)-MLREE2.PROG(1)
            MLDER.PROG(1)=(3.D0/DX)*((DY/DX)-XCDER(1))
         ELSE
* Condition de bord naturelle (dérivée seconde nulle)
            MLDERS.PROG(1)=XZERO
            MLDER.PROG(1)=XZERO
         ENDIF
         DO ILON=2,LON-1
            XIM=MLREE1.PROG(ILON-1)
            XI=MLREE1.PROG(ILON)
            XIP=MLREE1.PROG(ILON+1)
            YIM=MLREE2.PROG(ILON-1)
            YI=MLREE2.PROG(ILON)
            YIP=MLREE2.PROG(ILON+1)
            DXIM=XI-XIM
            DXI2=XIP-XIM
            DXIP=XIP-XI
            XRAP=DXIM/DXI2
            XP=XRAP*MLDERS.PROG(ILON-1)+2.D0
            MLDERS.PROG(ILON)=(XRAP-1.D0)/XP
            DYIP=YIP-YI
            DYIM=YI-YIM
            MLDER.PROG(ILON)=(6.D0*(DYIP/DXIP-DYIM/DXIM)/DXI2-XRAP
     $           *MLDER.PROG(ILON-1))/XP
         ENDDO
         IF (LCDER(2)) THEN
            XQN=0.5D0
            DX=MLREE1.PROG(LON)-MLREE1.PROG(LON-1)
            DY=MLREE2.PROG(LON)-MLREE2.PROG(LON-1)
            XUN=(3.D0/DX)*(XCDER(2)-(DY/DX))
         ELSE
* Condition de bord naturelle (dérivée seconde nulle)
            XQN=0.D0
            XUN=0.D0
         ENDIF
         MLDERS.PROG(LON)=(XUN-XQN*MLDER.PROG(LON-1))/
     $        (XQN*MLDERS.PROG(LON-1)+1.D0)
         DO ILON=LON-1,1,-1
            MLDERS.PROG(ILON)=MLDERS.PROG(ILON)*MLDERS.PROG(ILON+1)
     $           +MLDER.PROG(ILON)
         ENDDO
         SEGSUP MLDER
      ELSE
         MLDERS=0
      ENDIF
* On ne va pas renvoyer d'erreur si t est en-dehors des bornes 
* On renvoie la valeur de la fonction pour t le plus proche
      IVERI=0
*
*     Répartition suivant le type de l'objet fourni
*
      GOTO (10,20,30,40) INDIC

******************   T0  FLOTTANT  *******************************

10    CONTINUE
      CALL INTER5(TEMPS,MLREE1,MLREE2,FT0,ISPLIN,MLDERS,IVERI,IRET)
      IF (IRET.EQ.1) CALL ECRREE(FT0)
      GOTO 999

*********************   T0  CHPOINT  *****************************

20    CONTINUE
      MCHPO1=MTEMP
      SEGINI,MCHPOI=MCHPO1
      MFT0=MCHPOI
      NS=IPCHP(/1)
      DO 21 IA=1,NS
        MSOUP1=IPCHP(IA)
        SEGINI,MSOUPO=MSOUP1
        NC=NOHARM(/1)
        IPCHP(IA)=MSOUPO
        IPT1=IGEOC
        SEGINI,IPT2=IPT1
        SEGDES IPT2
        IGEOC=IPT2
        MPOVA1=IPOVAL
        SEGINI,MPOVAL=MPOVA1
        IPOVAL=MPOVAL
        N=VPOCHA(/1)
        DO IB=1,N
          DO IC=1,NC
            TEMPS=VPOCHA(IB,IC)
            CALL INTER5(TEMPS,MLREE1,MLREE2,FT0,ISPLIN,MLDERS,
     $           IVERI,IRET)
            IF (IRET.EQ.0) THEN
              SEGSUP MCHPOI,MSOUPO,MPOVAL
              GOTO 999
            ENDIF
            VPOCHA(IB,IC)=FT0
          ENDDO
        ENDDO
        SEGDES MPOVAL
        SEGDES MSOUPO
21    CONTINUE
      SEGDES MCHPOI
      CALL ECROBJ('CHPOINT ',MFT0)
      GOTO 999

*******************   T0 EST UN LISTREEL  ***********************

30    CONTINUE
      MLREE3=MLIST
      SEGACT MLREE3
      LONG=MLREE3.PROG(/1)
      JG=LONG
      SEGINI MLREEL
      MSOL=MLREEL
      DO 31 ILOOP=1,LONG
        TEMPS=MLREE3.PROG(ILOOP)
        CALL INTER5(TEMPS,MLREE1,MLREE2,FTO,ISPLIN,MLDERS,
     $       IVERI,IRET)
        IF(IRET.EQ.0) GOTO 999
        PROG(ILOOP)=FTO
31    CONTINUE
      SEGDES MLREEL,MLREE3
      CALL ECROBJ('LISTREEL',MSOL)
      GOTO 999

***********************   T0  MCHAML   ***************************

40    CONTINUE
      IRET=0
      MCHEL1=IPO1
      SEGINI,MCHELM=MCHEL1
      IRET=MCHELM
      NSOUS=IMACHE(/1)
      DO 72 IA=1,NSOUS
        MCHAM1=ICHAML(IA)
        SEGINI,MCHAML=MCHAM1
        ICHAML(IA)=MCHAML
        DO 75 ICOMP=1,IELVAL(/1)
          MELVA1 = IELVAL(ICOMP)
          SEGINI,MELVAL=MELVA1
          IELVAL(ICOMP) = MELVAL
          SEGACT MELVA1
          IF (TYPCHE(ICOMP).NE.'REAL*8') GOTO 166
          N1PTEL=VELCHE(/1)
          N1EL  =VELCHE(/2)
          N2PTEL=0
          N2EL  =0
          DO IB=1,N1PTEL
            DO ID=1,N1EL
              TEMPS=MELVA1.VELCHE(IB,ID)
              CALL INTER5(TEMPS,MLREE1,MLREE2,FT0,ISPLIN,
     $             MLDERS,IVERI,IREE)
              IF (IREE.EQ.0) THEN
                SEGSUP MCHELM,MCHAML,MELVAL
                SEGDES MELVA1
                GOTO 999
              ENDIF
              VELCHE(IB,ID)=FT0
            ENDDO
          ENDDO
166       SEGDES MELVAL,MELVA1
75      CONTINUE
        SEGDES MCHAML
72    CONTINUE
      SEGDES MCHELM
      CALL ECROBJ('MCHAML',IRET)
      GOTO 999

************************  OBJET1  TABLE  ******************************

50    CONTINUE
      MTABLE = IPOT
      SEGACT MTABLE
      LO = MLOTAB

*--   Vérification du format de la table
      IF (LO.LE.2) THEN
*       La table n'a pas le format désiré
        CALL ERREUR(647)
        SEGDES MTABLE
        RETURN
      ENDIF

      LOR = LO
      SEGINI TABR

*--   Vérification du sous-type de la table
*     IMOT est son indice dans la table
      IOK = 0
      DO 55 I=1,LO
        TYPE1 = MTABTI(I)
        IF(TYPE1.EQ.'MOT     ') THEN
          CHARIN = 'SOUSTYPE'
          TYPOBJ = ' '
          CALL ACCTAB(MTABLE,TYPE1,IVALIN,XVALIN,CHARIN,LOGIN,
     $               IOBIN,TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
          segact mtable
          IF(CHARIN.EQ.'SOUSTYPE') THEN
            IF(CHARRE.EQ.'RESULTAT') THEN
              IOK  = 1
              IMOT = I
            ENDIF
          ENDIF
        ENDIF
55    CONTINUE

      IF(IOK.EQ.0) THEN
*       Le sous-type de la table est incorrect
        CALL ERREUR(648)
        SEGDES MTABLE
        SEGSUP TABR
        RETURN
      ENDIF

*--   Vérification qu'on a bien des flottants en indice
*     en dehors du sous-type
      J = 0
      DO 56 I=1,LO
        IF (IMOT.EQ.I) GOTO 56
        J=J+1
        TYPEI = MTABTI(I)

        IF (TYPEI.NE.'FLOTTANT') THEN
*         La table n'a pas le format desire
          CALL ERREUR(647)
          SEGDES MTABLE
          SEGSUP TABR
          RETURN
        ENDIF

        TEMA(J) = RMTABI(I)
56    CONTINUE

*--   Vérification de l'ordonnancement des indices
      DO 57 I=1,LOR-2
        TEM1 = TEMA(I)
        TEM2 = TEMA(I+1)

        IF(TEM1.GT.TEM2) THEN
*         la liste des indices n'est pas ordonnee
          CALL ERREUR(211)
          SEGDES MTABLE
          SEGSUP TABR
          RETURN
        ENDIF

57    CONTINUE

      IF(IMOT.EQ.1) THEN
        TEM1   = RMTABI(2)
        TYPO1  = MTABTV(2)
        IVALO1 = MTABIV(2)
        NDEB   = 3
      ENDIF
      IF(IMOT.EQ.2) THEN
        TEM1   = RMTABI(1)
        TYPO1  = MTABTV(1)
        IVALO1 = MTABIV(1)
        NDEB   = 3
      ENDIF
      IF(IMOT.GE.3) THEN
        TEM1   = RMTABI(1)
        TYPO1  = MTABTV(1)
        IVALO1 = MTABIV(1)
        NDEB   = 2
      ENDIF

      DO 58 I=NDEB,LOR
        IF(IMOT.EQ.I) GOTO 58
        TEM2   = RMTABI(I)
        TYPO2  = MTABTV(I)
        IVALO2 = MTABIV(I)
        IF((TEM1.LE.TEMPS).AND.(TEMPS.LE.TEM2)) THEN
          DTEM = (TEMPS-TEM1)/(TEM2-TEM1)
          IF(TYPO1.EQ.'CHPOINT ') THEN
            IF(TYPO2.EQ.'CHPOINT ') THEN
              CALL ECROBJ('CHPOINT',IVALO2)
              CALL ECROBJ('CHPOINT',IVALO1)
            ELSE
              CALL ERREUR(647)
              SEGDES MTABLE
              SEGSUP TABR
              RETURN
            ENDIF
          ENDIF
          IF(TYPO1.EQ.'MCHAML  ') THEN
            IF(TYPO2.EQ.'MCHAML  ') THEN
              CALL ECROBJ('MCHAML',IVALO2)
              CALL ECROBJ('MCHAML',IVALO1)
            ELSE
              CALL ERREUR(647)
              SEGDES MTABLE
              SEGSUP TABR
              RETURN
            ENDIF
          ENDIF
          CALL ECRREE(DTEM)
          CALL ECRREE(1.D0-DTEM)
          CALL COLI
          SEGDES MTABLE
          SEGSUP TABR
          GOTO 500
        ENDIF
        TEM1   = TEM2
        TYPO1  = TYPO2
        IVALO1 = IVALO2
58    CONTINUE

*     Le temps est en dehors des limites de la table
      CALL ERREUR(210)
      SEGDES MTABLE
      SEGSUP TABR
      GOTO 500

999   CONTINUE
      IF (ISENS.EQ.0) THEN
         SEGDES MLREE1
         SEGDES MLREE2
      ELSE
         SEGSUP MLREE1
         SEGSUP MLREE2
      ENDIF
      IF (ISPLIN.EQ.1) THEN
         SEGSUP MLDERS
      ENDIF
*     Sortie
500   CONTINUE
      RETURN
      END