C IPPID     SOURCE    BP208322  21/01/28    21:15:14     10868          
        SUBROUTINE IPPID(INUA,NBC,MTR1,X1,Y1,N,INDJ,NN,MM,NNMM,XP,XELIM)
C-----------------------------------------------------------------------
C NOM         : IPPID
C DESCRIPTION : Interpolation dans un NUAGE par Ponderation Inverse a la
C               Distance
C LANGAGE     : ESOPE
C AUTEUR      : Francois DI PAOLA
C-----------------------------------------------------------------------
C APPELE PAR  : IPLNU2
C APPELE      :
C-----------------------------------------------------------------------
C ENTREES     : INUA (Objet de type NUAGE)
C               NBC nombe de points (cad de NNMM-uplets) du nuage
C               NNMM nombe de composantes du nuage
C               NN nombre de composantes connues du point cible
C               MM nombre de composantes a calculer du point cible
C               MTR1.IADD contient la correspondance entre les numeros
C               des composantes du point cible et celles du nuage :
C               - La composante connue I (entre 1 et NN) est a la
C                 position IADD(I) dans le nuage
C               - La composante a calculer J (entre 1 et MM) est a la
C                 position IADD(J) dans le nuage
C               X1 tableau contenant les valeurs des NN composantes
C               connues au point cible
C               XP est l'exposant choisit pour la ponderation
C SORTIES     : Y1 tableau contenant les valeurs des MM composantes
C               interpolees au point cible
C-----------------------------------------------------------------------
C VERSION    : v1, 27/05/2016, version initiale
C HISTORIQUE : v1, 27/05/2016, creation
C HISTORIQUE :
C HISTORIQUE :
C-----------------------------------------------------------------------
C Priere de PRENDRE LE TEMPS de completer les commentaires
C en cas de modification de ce sous-programme afin de faciliter
C la maintenance !
C-----------------------------------------------------------------------
C REMARQUES : - L'interpolation est exacte, c'est-a-dire que l'on
C               retrouve les valeurs du nuage si l'on interpole en point
C               du nuage
C             - Pour plus d'information consultez :
C               * Shepard, Donald (1968). « A two-dimensional interpolation
C                 function for irregularly-spaced data ». Proceedings of
C                 the 1968 ACM National Conference: 517–524.
C               https://fr.wikipedia.org/wiki/Pond%C3%A9ration_inverse_%C3%A0_la_distance
C               https://en.wikipedia.org/wiki/Inverse_distance_weighting
C-----------------------------------------------------------------------
C
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT    REAL*8(A-H,O-Z)
-INC CCREEL

-INC PPARAM
-INC CCOPTIO
-INC SMNUAGE
      SEGMENT MTR1
        INTEGER IADD(NNMM)
      ENDSEGMENT
      REAL*8  X1,Y1
      DIMENSION X1(N,*)
      DIMENSION Y1(N,*)
      REAL*8  DIS(NBC),W(NBC)
C
C     Critere pour le calcul des poids avec des distance trop grandes
      XGRANDR=XGRAND**(1.D0/ABS(XP))
C
C     Calcul du nombre de points dans le nuage
      MNUAG1=INUA
C
C     Si nuage vide, operation impossible
      IF (NBC.EQ.0) THEN
        CALL ERREUR(680)
        RETURN
      ENDIF
C
C     Calcul preliminaire des distances euclidiennes des points du nuage
C     au point XI (sur les NN composantes connues)
      NZER=0
C     Boucle sur les points du nuage
      DO I=1,NBC
        DI=XZERO
C       Boucle sur les NN composantes connues du champ et calcul de la
C       distance euclidienne du point I
        DO J=1,NN
          X1J=X1(INDJ,J)
          JJ=IADD(J)
          NUAVF1=MNUAG1.NUAPOI(JJ)
          XNJ=NUAVF1.NUAFLO(I)
          DXJ=X1J-XNJ
          DI=DI+(DXJ*DXJ)
        ENDDO
        DI=SQRT(DI)
C       Si le point I du nuage est a une distance nulle, on arrete le
C       calcul des distances car il pese pour la totalite de
C       l'interpolation, ce qui permet d'avoir une interpolation exacte 
c       -> facile : on fait le travail ici !
        IF (ABS(DI).LT.XELIM) THEN
c         boucle sur les MM composantes a traiter
          DO J=1,MM 
            II=IADD(J+NN)
            NUAVF1=MNUAG1.NUAPOI(II)
            Y1(INDJ,J)=NUAVF1.NUAFLO(I)
          ENDDO
          RETURN
        ENDIF
c       si le point n'est pas a une distance nulle on stocke la distance
        DIS(I)=DI
      ENDDO
      
C     Calcul des ponderations associees a chaque point du nuage
C     (inverses a la distance)
C     Si un point du nuage est a une distance nulle, 
      WSOM=XZERO
      DO I=1,NBC
C       Si la distance est trop elevee, erreur division par zero
        IF (DIS(I).GT.XGRANDR) THEN
          CALL ERREUR(835)
          RETURN
        ENDIF
C       On peut diviser plus sereinement
        W(I)=(1.D0/DIS(I))**(XP)
        WSOM=WSOM+W(I)
      ENDDO
C
C     Calcul des valeurs des MM composantes inconnues au point X1
      DO J=1,MM
        FJ=XZERO
        II=IADD(J+NN)
        NUAVF1=MNUAG1.NUAPOI(II)
C       Boucle sur les points du nuage
        DO I=1,NBC
          WI=W(I)
          FI=NUAVF1.NUAFLO(I)
          FJ=FJ+(WI*FI)
        ENDDO
        FJ=FJ/WSOM
        Y1(INDJ,J)=FJ
      ENDDO
      RETURN
      END