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