Télécharger ippid.eso

Retour à la liste

Numérotation des lignes :

  1. C IPPID SOURCE PV 16/09/13 13:11:50 9084
  2. SUBROUTINE IPPID(INUA,NBC,MTR1,X1,Y1,N,INDJ,NN,MM,NNMM,XP)
  3. C-----------------------------------------------------------------------
  4. C NOM : IPPID
  5. C DESCRIPTION : Interpolation dans un NUAGE par Ponderation Inverse a la
  6. C Distance
  7. C LANGAGE : ESOPE
  8. C AUTEUR : Francois DI PAOLA
  9. C-----------------------------------------------------------------------
  10. C APPELE PAR : IPLNU2
  11. C APPELE :
  12. C-----------------------------------------------------------------------
  13. C ENTREES : INUA (Objet de type NUAGE)
  14. C NBC nombe de points (cad de NNMM-uplets) du nuage
  15. C NNMM nombe de composantes du nuage
  16. C NN nombre de composantes connues du point cible
  17. C MM nombre de composantes a calculer du point cible
  18. C MTR1.IADD contient la correspondance entre les numeros
  19. C des composantes du point cible et celles du nuage :
  20. C - La composante connue I (entre 1 et NN) est a la
  21. C position IADD(I) dans le nuage
  22. C - La composante a calculer J (entre 1 et MM) est a la
  23. C position IADD(J) dans le nuage
  24. C X1 tableau contenant les valeurs des NN composantes
  25. C connues au point cible
  26. C XP est l'exposant choisit pour la ponderation
  27. C SORTIES : Y1 tableau contenant les valeurs des MM composantes
  28. C interpolees au point cible
  29. C-----------------------------------------------------------------------
  30. C VERSION : v1, 27/05/2016, version initiale
  31. C HISTORIQUE : v1, 27/05/2016, creation
  32. C HISTORIQUE :
  33. C HISTORIQUE :
  34. C-----------------------------------------------------------------------
  35. C Priere de PRENDRE LE TEMPS de completer les commentaires
  36. C en cas de modification de ce sous-programme afin de faciliter
  37. C la maintenance !
  38. C-----------------------------------------------------------------------
  39. C REMARQUES : - L'interpolation est exacte, c'est-a-dire que l'on
  40. C retrouve les valeurs du nuage si l'on interpole en point
  41. C du nuage
  42. C - Pour plus d'information consultez :
  43. C * Shepard, Donald (1968). « A two-dimensional interpolation
  44. C function for irregularly-spaced data ». Proceedings of
  45. C the 1968 ACM National Conference: 517–524.
  46. C https://fr.wikipedia.org/wiki/Pond%C3%A9ration_inverse_%C3%A0_la_distance
  47. C https://en.wikipedia.org/wiki/Inverse_distance_weighting
  48. C-----------------------------------------------------------------------
  49. C
  50. C
  51. IMPLICIT INTEGER(I-N)
  52. IMPLICIT REAL*8(A-H,O-Z)
  53. -INC CCREEL
  54. -INC CCOPTIO
  55. -INC SMNUAGE
  56. SEGMENT MTR1
  57. INTEGER IADD(NNMM)
  58. ENDSEGMENT
  59. REAL*8 X1,Y1
  60. DIMENSION X1(N,*)
  61. DIMENSION Y1(N,*)
  62. REAL*8 DIS(NBC),W(NBC)
  63. C
  64. C Critere pour le calcul des poids avec des distance trop grandes
  65. XGRANDR=XGRAND**(1.D0/ABS(XP))
  66. C
  67. C Calcul du nombre de points dans le nuage
  68. MNUAG1=INUA
  69. C
  70. C Si nuage vide, operation impossible
  71. IF (NBC.EQ.0) THEN
  72. CALL ERREUR(680)
  73. RETURN
  74. ENDIF
  75. C
  76. C Calcul preliminaire des distances euclidiennes des points du nuage
  77. C au point XI (sur les NN composantes connues)
  78. NZER=0
  79. C Boucle sur les points du nuage
  80. DO I=1,NBC
  81. DI=XZERO
  82. C Boucle sur les NN composantes connues du champ et calcul de la
  83. C distance euclidienne du point I
  84. DO J=1,NN
  85. X1J=X1(INDJ,J)
  86. JJ=IADD(J)
  87. NUAVF1=MNUAG1.NUAPOI(JJ)
  88. XNJ=NUAVF1.NUAFLO(I)
  89. DXJ=X1J-XNJ
  90. DI=DI+(DXJ*DXJ)
  91. ENDDO
  92. DI=SQRT(DI)
  93. C Si le point I du nuage est a une distance nulle, on arrete le
  94. C calcul des distances
  95. IF (ABS(DI).LT.XPETIT) THEN
  96. NZER=I
  97. GOTO 10
  98. ENDIF
  99. DIS(I)=DI
  100. ENDDO
  101. 10 CONTINUE
  102. C
  103. C Calcul des ponderations associees a chaque point du nuage
  104. C (inverses a la distance)
  105. C Si un point du nuage est a une distance nulle, il pese pour la
  106. C totalite de l'interpolation, ce qui permet d'avoir une
  107. C interpolation exacte !
  108. WSOM=XZERO
  109. DO I=1,NBC
  110. IF (NZER.NE.0) THEN
  111. IF (I.EQ.NZER) THEN
  112. W(I)=1.D0
  113. ELSE
  114. W(I)=XZERO
  115. ENDIF
  116. ELSE
  117. C Si la distance est trop elevee, erreur division par zero
  118. IF (DIS(I).GT.XGRANDR) THEN
  119. CALL ERREUR(835)
  120. RETURN
  121. ENDIF
  122. C On peut diviser plus sereinement
  123. W(I)=(1.D0/DIS(I))**(XP)
  124. ENDIF
  125. WSOM=WSOM+W(I)
  126. ENDDO
  127. C
  128. C Calcul des valeurs des MM composantes inconnues au point X1
  129. DO I=1,MM
  130. FI=XZERO
  131. II=IADD(I+NN)
  132. NUAVF1=MNUAG1.NUAPOI(II)
  133. C Boucle sur les points du nuage
  134. DO J=1,NBC
  135. WJ=W(J)
  136. FJ=NUAVF1.NUAFLO(J)
  137. FI=FI+(WJ*FJ)
  138. ENDDO
  139. FI=FI/WSOM
  140. Y1(INDJ,I)=FI
  141. ENDDO
  142. RETURN
  143. END
  144.  
  145.  
  146.  

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