Télécharger ippid.eso

Retour à la liste

Numérotation des lignes :

ippid
  1. C IPPID SOURCE BP208322 21/01/28 21:15:14 10868
  2. SUBROUTINE IPPID(INUA,NBC,MTR1,X1,Y1,N,INDJ,NN,MM,NNMM,XP,XELIM)
  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.  
  55. -INC PPARAM
  56. -INC CCOPTIO
  57. -INC SMNUAGE
  58. SEGMENT MTR1
  59. INTEGER IADD(NNMM)
  60. ENDSEGMENT
  61. REAL*8 X1,Y1
  62. DIMENSION X1(N,*)
  63. DIMENSION Y1(N,*)
  64. REAL*8 DIS(NBC),W(NBC)
  65. C
  66. C Critere pour le calcul des poids avec des distance trop grandes
  67. XGRANDR=XGRAND**(1.D0/ABS(XP))
  68. C
  69. C Calcul du nombre de points dans le nuage
  70. MNUAG1=INUA
  71. C
  72. C Si nuage vide, operation impossible
  73. IF (NBC.EQ.0) THEN
  74. CALL ERREUR(680)
  75. RETURN
  76. ENDIF
  77. C
  78. C Calcul preliminaire des distances euclidiennes des points du nuage
  79. C au point XI (sur les NN composantes connues)
  80. NZER=0
  81. C Boucle sur les points du nuage
  82. DO I=1,NBC
  83. DI=XZERO
  84. C Boucle sur les NN composantes connues du champ et calcul de la
  85. C distance euclidienne du point I
  86. DO J=1,NN
  87. X1J=X1(INDJ,J)
  88. JJ=IADD(J)
  89. NUAVF1=MNUAG1.NUAPOI(JJ)
  90. XNJ=NUAVF1.NUAFLO(I)
  91. DXJ=X1J-XNJ
  92. DI=DI+(DXJ*DXJ)
  93. ENDDO
  94. DI=SQRT(DI)
  95. C Si le point I du nuage est a une distance nulle, on arrete le
  96. C calcul des distances car il pese pour la totalite de
  97. C l'interpolation, ce qui permet d'avoir une interpolation exacte
  98. c -> facile : on fait le travail ici !
  99. IF (ABS(DI).LT.XELIM) THEN
  100. c boucle sur les MM composantes a traiter
  101. DO J=1,MM
  102. II=IADD(J+NN)
  103. NUAVF1=MNUAG1.NUAPOI(II)
  104. Y1(INDJ,J)=NUAVF1.NUAFLO(I)
  105. ENDDO
  106. RETURN
  107. ENDIF
  108. c si le point n'est pas a une distance nulle on stocke la distance
  109. DIS(I)=DI
  110. ENDDO
  111.  
  112. C Calcul des ponderations associees a chaque point du nuage
  113. C (inverses a la distance)
  114. C Si un point du nuage est a une distance nulle,
  115. WSOM=XZERO
  116. DO I=1,NBC
  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. WSOM=WSOM+W(I)
  125. ENDDO
  126. C
  127. C Calcul des valeurs des MM composantes inconnues au point X1
  128. DO J=1,MM
  129. FJ=XZERO
  130. II=IADD(J+NN)
  131. NUAVF1=MNUAG1.NUAPOI(II)
  132. C Boucle sur les points du nuage
  133. DO I=1,NBC
  134. WI=W(I)
  135. FI=NUAVF1.NUAFLO(I)
  136. FJ=FJ+(WI*FI)
  137. ENDDO
  138. FJ=FJ/WSOM
  139. Y1(INDJ,J)=FJ
  140. ENDDO
  141. RETURN
  142. END
  143.  
  144.  
  145.  
  146.  

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