Télécharger ipmul.eso

Retour à la liste

Numérotation des lignes :

  1. C IPMUL SOURCE FD218221 15/10/05 21:15:06 8658
  2. C-----------------------------------------------------------------------
  3. C NOM : IPMUL
  4. C DESCRIPTION : Interpolation multi-lineaire dans un NUAGE represantant
  5. C une grille de valeurs
  6. C LANGAGE : ESOPE
  7. C AUTEUR : Francois DI PAOLA
  8. C-----------------------------------------------------------------------
  9. C APPELE PAR : IPMULI via IPGRIL
  10. C APPELE :
  11. C-----------------------------------------------------------------------
  12. C ENTREES
  13. C MVAL1 : Tableau des N coordonnees du point ou l'on souhaite faire
  14. C l'interpolation
  15. C MLENT1 : Tableau de (N+1) LISTREELs :
  16. C - N listes des coordonnees des noeuds de la grille
  17. C - 1 liste des valeurs de la fonction sur chaque noeud
  18. C Les N listes de coordonnees sont rangees dans le meme
  19. C ordre que MVAL1
  20. C SORTIES
  21. C XVAL : Valeur de la fonction interpolee au point MVAL1
  22. C-----------------------------------------------------------------------
  23. C VERSION : v1, 05/10/2015, version initiale
  24. C HISTORIQUE : v1, 05/10/2015, creation
  25. C HISTORIQUE :
  26. C HISTORIQUE :
  27. C-----------------------------------------------------------------------
  28. C Priere de PRENDRE LE TEMPS de completer les commentaires
  29. C en cas de modification de ce sous-programme afin de faciliter
  30. C la maintenance !
  31. C-----------------------------------------------------------------------
  32. C REMARQUES : - L'interpolation est exacte, c'est-a-dire que l'on
  33. C retrouve les valeurs de la grille si l'on interpole en
  34. C un noeud de la grille
  35. C - La grille peut contenir autant de dimensions que
  36. C souhaitees
  37. C-----------------------------------------------------------------------
  38. C
  39. SUBROUTINE IPMUL(VAL,NBC,MLENT1,XVAL)
  40. C
  41. IMPLICIT INTEGER(I-N)
  42. IMPLICIT REAL*8 (A-H,O-Z)
  43. -INC CCNOYAU
  44. -INC SMLENTI
  45. -INC SMLREEL
  46. C Le tableau VAL peremt de stocker les coordonnees du point
  47. C d'interpolation
  48. C VAL(i) --> valeur de la coordonnee i du point d'interpolation
  49. REAL*8 VAL(NBC)
  50. C
  51. C Deux tableaux pour reperer les 2^n noeuds encadrant le point ou l'on
  52. C veut interpoler
  53. C IND(i,1) --> indice dans la liste de la coordonne i du noeud situe
  54. C en dessous
  55. C IND(i,2) --> indice dans la liste de la coordonne i du noeud situe
  56. C au dessus
  57. C VAL(i,1) --> valeur de la coordonne i du noeud situe en dessous
  58. C VAL(i,2) --> valeur de la coordonne i du noeud situe au dessus
  59. INTEGER IND(NBC,2)
  60. REAL*8 FVA(NBC,2)
  61. C
  62. C Un dernier pour stocker la valeur, en base 2, des 2^n entiers
  63. C IBIT(i) --> valeur du ieme bit (0 ou 1)
  64. INTEGER IBIT(NBC)
  65. C
  66. C ----------------------------------------------------------------
  67. C PARTIE 1 : Initialisations
  68. C - determination de la boite des noeuds encadrant le
  69. C point d'interpolation
  70. C - stockage des coordonnes des noeuds encadrant et des
  71. C valeurs de la fonction
  72. C ----------------------------------------------------------------
  73. C Recherche, pour chaque coordonnee, des 2 noeuds de la grille
  74. C encadrant le point d'interpolation -> remplissage de IND et FVA
  75. WTOT=1D0
  76. DO I=1,NBC
  77. X=VAL(I)
  78. C Recuperation, du LISTREEL de la grille correspondant la
  79. C coordonnee I
  80. MLREE1=MLENT1.LECT(I)
  81. SEGACT,MLREE1
  82. N1=MLREE1.PROG(/1)
  83. C Si la liste n'est pas de dimension 2, erreur
  84. IF (N1.LT.2) THEN
  85. CALL ERREUR(897)
  86. RETURN
  87. ENDIF
  88. C Si le point est en dehors des noeuds, on le projette sur le
  89. C noeud le plus proche
  90. XMIN=MLREE1.PROG(1)
  91. IF (X.LE.XMIN) THEN
  92. I1=1
  93. I2=2
  94. X1=XMIN
  95. X2=MLREE1.PROG(I2)
  96. VAL(I)=XMIN
  97. GOTO 10
  98. ENDIF
  99. XMAX=MLREE1.PROG(N1)
  100. IF (X.GE.XMAX) THEN
  101. I1=N1-1
  102. I2=N1
  103. X1=MLREE1.PROG(I1)
  104. X2=XMAX
  105. VAL(I)=XMAX
  106. GOTO 10
  107. ENDIF
  108. C Parcours de la liste a la recherche des valeurs encadrantes
  109. DO J=1,(N1-1)
  110. X1=MLREE1.PROG(J)
  111. X2=MLREE1.PROG(J+1)
  112. IF ((X1.LE.X).AND.(X.LE.X2)) THEN
  113. I1=J
  114. I2=J+1
  115. GOTO 10
  116. ENDIF
  117. ENDDO
  118. 10 CONTINUE
  119. IND(I,1)=I1
  120. IND(I,2)=I2
  121. FVA(I,1)=X1
  122. FVA(I,2)=X2
  123. C Calcul du "volume" de la boite de noeuds encadrante
  124. WTOT=WTOT*(MLREE1.PROG(I2)-MLREE1.PROG(I1))
  125. ENDDO
  126. C
  127. C -------------------------------------------------------------
  128. C PARTIE 2 : Interpolation
  129. C - recuperation des valeurs de F aux 2^n noeuds encadrant
  130. C - calcul des ponderations associees a chaque noeud
  131. C -------------------------------------------------------------
  132. C Initialisation du resultat de l'interpolation
  133. XVAL=0D0
  134. N1=2**NBC
  135. DO I=0,N1-1
  136. C Conversion de I en base 2 (calcul des NBC bits)
  137. II=I
  138. DO J=1,NBC
  139. IQ=II/2
  140. IR=II-(IQ*2)
  141. IBIT(NBC+1-J)=IR
  142. II=IQ
  143. ENDDO
  144. C Recuperation de la valeur de la fonction au noeud encadrant I
  145. III=1
  146. DO J=1,NBC
  147. JB=IBIT(J)+1
  148. JL=IND(J,JB)
  149. JP=1
  150. IF (J.GT.1) THEN
  151. DO K=1,J-1
  152. MLREE1=MLENT1.LECT(K)
  153. JP=JP*(MLREE1.PROG(/1))
  154. ENDDO
  155. ENDIF
  156. III=III+((JL-1)*JP)
  157. ENDDO
  158. MLREE1=MLENT1.LECT(NBC+1)
  159. FI=MLREE1.PROG(III)
  160. C Calcul de la ponderation du noeud I, c'est le "volume" du
  161. C polytope diagonalement oppose au noeud I par rapport au point
  162. C d'interpolation
  163. WI=1D0
  164. DO J=1,NBC
  165. XJ=VAL(J)
  166. C sur la coord J, le noeud I est il a gauche ou a droite ?
  167. JB=IBIT(J)
  168. IF (JB.EQ.0) THEN
  169. XJD=FVA(J,2)
  170. WJ=XJD-XJ
  171. ELSE
  172. XJG=FVA(J,1)
  173. WJ=XJ-XJG
  174. ENDIF
  175. WI=WI*WJ
  176. ENDDO
  177. C La contribution du noeud I dans l'interpolation est FI*WI
  178. XVAL=XVAL+(FI*WI)
  179. ENDDO
  180. C Enfin, on divise par la taille de la boite encadrante
  181. XVAL=XVAL/WTOT
  182. C
  183. RETURN
  184. END
  185.  
  186.  
  187.  

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