Télécharger ipol_pid.dgibi

Retour à la liste

Numérotation des lignes :

  1. OPTI 'DIME' 3 'ELEM' 'CUB8' 'ECHO' 0 ;
  2.  
  3.  
  4.  
  5. ************************************************************************
  6. * Test de l'operateur IPOL option 'PID' *
  7. * interpolation multi-dimensionnelle *
  8. * par Ponderation Inverse a la Distance *
  9. * *
  10. * Application a l'interpolation d'une fonction de 1 variable, puis *
  11. * une fonction de 2 variables d'espace (sinus cardinal) *
  12. * - test avec des CHPOINTs et des MCHAMLs *
  13. * - test sur des points aleatoires *
  14. * - test sur des points identiques au nuage *
  15. ************************************************************************
  16.  
  17.  
  18.  
  19. ** Indicateur pour le trace des champs interpoles
  20. itrac = FAUX ;
  21.  
  22. ** Tolerance pour les tests des valeurs interpolees
  23. tol1 = 1.E-12 ;
  24.  
  25.  
  26.  
  27.  
  28.  
  29. ************************************************************************
  30. * PETIT TEST EN 1D *
  31. ************************************************************************
  32.  
  33. ** Maillage d'une ligne
  34. p1 = 0. 0. 0. ;
  35. p2 = 10. 0. 0. ;
  36. l1 = DROI 3 p1 p2 ;
  37.  
  38. ** Fonction sinus cardinal
  39. x = COOR 1 l1 ;
  40. p0 = x POIN 'EGALE' 0. ;
  41. mp1 = CHAN 'POI1' l1 ;
  42. mp1b = DIFF mp1 p0 ;
  43. xb = REDU x mp1b ;
  44. f1 = (SIN (xb*180./ pi)) / xb ;
  45. f1b = MANU 'CHPO' p0 1 'SCAL' 1. 'NATURE' 'DIFFUS' ;
  46. f1 = f1 ET f1b ;
  47. ev0 = EVOL 'CHPO' f1 l1 ;
  48.  
  49. ** Nuage contenant l'information de depart
  50. ch1 = (NOMC 'X' x) ET (NOMC 'F1' f1) ;
  51. nu1 = NUAG ch1 ;
  52.  
  53. ** Points ou effectuer l'interpolation
  54. p3 = -5. 0. 0. ;
  55. p4 = 15. 0. 0. ;
  56. l2 = DROI 10000 p3 p4 ;
  57.  
  58. ** Interpolation
  59. x2 = NOMC (COOR 1 l2) 'X' ;
  60. np = 2. ;
  61. f2p = IPOL nu1 x2 'PID' np ;
  62. f2g = IPOL nu1 x2 'GAUSS' ;
  63. f2r = IPOL nu1 x2 'RATIO' ;
  64.  
  65. ** Trace des graphes
  66. SI itrac ;
  67. lx2 = EXTR (EVOL 'CHPO' x2 l2) 'ORDO' ;
  68. lp = EXTR (EVOL 'CHPO' f2p l2) 'ORDO' ;
  69. lg = EXTR (EVOL 'CHPO' f2g l2) 'ORDO' ;
  70. lr = EXTR (EVOL 'CHPO' f2r l2) 'ORDO' ;
  71. evp = EVOL 'ROUG' 'MANU' lx2 lp ;
  72. evg = EVOL 'VERT' 'MANU' lx2 lg ;
  73. evr = EVOL 'JAUN' 'MANU' lx2 lr ;
  74. tit1 = CHAI 'IPOL par PID (rouge), par GAUSS (vert) '
  75. 'et par RATIO (jaune)' ;
  76. tleg = TABL ;
  77. tleg . 1 = MOT 'MARQ CARR NOLI' ;
  78. tleg . TITRE = TABL ;
  79. tleg . TITRE . 1 = 'Sinus cardinal' ;
  80. tleg . TITRE . 2 = CHAI 'PID avec p = ' np ;
  81. tleg . TITRE . 3 = 'GAUSS' ;
  82. tleg . TITRE . 4 = 'RATIO' ;
  83. DESS (ev0 ET evp ET evg ET evr) 'TITR' tit1 'LEGE' tleg ;
  84. FINSI ;
  85.  
  86.  
  87.  
  88.  
  89.  
  90. ************************************************************************
  91. * PETIT TEST EN 2D *
  92. ************************************************************************
  93.  
  94. ** CREATION DU NUAGE
  95. ** -----------------
  96. ** Un petit maillage
  97. ne = 30 ;
  98. xmax = 10. ;
  99. mxmax = -1. * xmax ;
  100. l1 = D ne (mxmax mxmax 0.) (xmax mxmax 0.) ;
  101. s1 = l1 TRAN ne (0. (2. * xmax) 0.) ;
  102. mp1 = CHAN 'POI1' s1 ;
  103.  
  104. ** Coordonnees x et y puis fonction (ici on prend un sinus cardinal
  105. * un peu tordu)
  106. x y z = COOR s1 ;
  107. r = ((x*x) + (y*y))**0.5 ;
  108. p0 = r POIN 'EGALE' 0. ;
  109. mp1b = DIFF mp1 p0 ;
  110. rb = REDU r mp1b ;
  111. xb = REDU x mp1b ;
  112. f1 = (10.*(SIN (rb*180./ pi))/rb) + (0.001*(xb**4)) ;
  113. f1b = MANU 'CHPO' p0 1 'SCAL' 10. 'NATURE' 'DIFFUS' ;
  114. f1 = f1 ET f1b ;
  115. ** Nuage de valeurs : ensemble de n-uplets (x,y,f1)
  116. ch1 = (NOMC 'X' x) ET (NOMC 'Y' y) ET (NOMC 'F1' f1) ;
  117. nu1 = NUAG ch1 ;
  118.  
  119.  
  120.  
  121. ** CREATION DES POINTS A CALCULER
  122. ** ------------------------------
  123. ** On tire nbp points dans les environs du nuage
  124. nbp = 100 ;
  125. lx = BRUI 'BLAN' 'GAUS' 0. (0.6 * xmax) nbp ;
  126. ly = BRUI 'BLAN' 'GAUS' 0. (0.6 * xmax) nbp ;
  127. lz = PROG nbp*0. ;
  128. mp2 = POIN lx ly lz ;
  129.  
  130. ** On y ajoute quelques points "singuliers" : presents dans le nuage
  131. mp3 = (s1 POIN 'PROC' (0. 0. 0.)) ET
  132. (s1 POIN 'PROC' (0. xmax 0.)) ET
  133. (s1 POIN 'PROC' (mxmax mxmax 0.)) ET
  134. (s1 POIN 'PROC' ((-0.5 * xmax) (0.2 * xmax) 0.)) ET
  135. (s1 POIN 'PROC' ((0.3 * xmax) (0.3 * xmax) 0.)) ;
  136. mp2 = mp2 ET mp3 ;
  137.  
  138.  
  139.  
  140. ** INTERPOLATION
  141. ** -------------
  142. ** Interpolation (dans un CHPOINT et un MCHAML)
  143. p = 3 ;
  144. xmp2 ymp2 zmp2 = COOR mp2 ;
  145. chp2 = (NOMC 'X' xmp2) ET (NOMC 'Y' ymp2) ;
  146. f2p = IPOL nu1 chp2 'PID' p ;
  147. chm2 = CHAN 'CHAM' chp2 mp2 ;
  148. f2m = IPOL nu1 chm2 'PID' p ;
  149.  
  150.  
  151.  
  152. ** COMPARAISON VALEURS CALCULEES / THEORIE
  153. ** ---------------------------------------
  154. maxerr = 0. ;
  155. MESS 'Numero | Valeur | Valeurs | Erreur ' ;
  156. MESS 'du point | theorique | interpolees | relative' ;
  157. MESS '--------------------------------------------------------' ;
  158. REPE b0 (NBEL mp2) ;
  159. p0 = mp2 POIN &b0 ;
  160. * valeurs calculees (CHPOINT et MCHAML)
  161. f0ca1 = EXTR f2p 'F1' p0 ;
  162. f0ca2 = EXTR f2m 'F1' 1 &b0 1 ;
  163. * valeur theorique
  164. x0 y0 z0 = COOR p0 ;
  165. lw = PROG ;
  166. f0th = 0. ;
  167. REPE b1 (NBEL mp1) ;
  168. p1 = mp1 POIN &b1 ;
  169. x1 = EXTR ch1 'X' p1 ;
  170. y1 = EXTR ch1 'Y' p1 ;
  171. d1 = (((x1 - x0)**2) + ((y1 - y0)**2))**0.5 ;
  172. SI (d1 <EG 1.E-15) ;
  173. lw = PROG 1. ;
  174. f0th = (EXTR ch1 'F1' p1) ;
  175. QUIT b1 ;
  176. FINSI ;
  177. w1 = 1. / (d1**p) ;
  178. lw = lw ET w1 ;
  179. f0th = f0th + (w1 * (EXTR ch1 'F1' p1)) ;
  180. FIN b1 ;
  181. wtot = SOMM lw ;
  182. f0th = f0th / wtot ;
  183. * comparaison
  184. err1 = ABS ((f0ca1 - f0th) / f0th) ;
  185. err2 = ABS ((f0ca2 - f0th) / f0th) ;
  186. MESS &b0 ' |' f0th '|' f0ca1 '|' err1 ;
  187. MESS ' | |' f0ca2 '|' err2 ;
  188. MESS '--------------------------------------------------------' ;
  189. maxerr = MAXI (PROG maxerr err1 err2) ;
  190. FIN b0 ;
  191.  
  192.  
  193.  
  194. ** TEST FINAL
  195. ** ----------
  196. SI (maxerr <EG tol1) ;
  197. MESS '***** Succes du cas test ! *****' ;
  198. SINON ;
  199. MESS '***** Echec du cas test ! *****' ; MESS ; MESS ;
  200. ERRE 5 ;
  201. FINSI ;
  202.  
  203.  
  204.  
  205. ** VISUALISATION DES POINTS CALCULES SUR LE NUAGE INITIAL
  206. ** ------------------------------------------------------
  207. SI itrac ;
  208. u1 = (MANU 'CHPO' s1 2 'UX' 0. 'UY' 0. 'NATURE' 'DIFFUS') ET
  209. (NOMC 'UZ' f1 'NATURE' 'DIFFUS') ;
  210. s1b = s1 PLUS u1 ;
  211. u2p = (MANU 'CHPO' mp2 2 'UX' 0. 'UY' 0. 'NATURE' 'DIFFUS') ET
  212. (NOMC 'UZ' f2p 'NATURE' 'DIFFUS') ;
  213. mp2p = (mp2 PLUS u2p) COUL 'ROUG' ;
  214. f2g = IPOL nu1 chp2 'GAUSS' ;
  215. u2g = (MANU 'CHPO' mp2 2 'UX' 0. 'UY' 0. 'NATURE' 'DIFFUS') ET
  216. (NOMC 'UZ' f2g 'NATURE' 'DIFFUS') ;
  217. mp2g = (mp2 PLUS u2g) COUL 'VERT' ;
  218. f2r = IPOL nu1 chp2 'RATIO' ;
  219. u2r = (MANU 'CHPO' mp2 2 'UX' 0. 'UY' 0. 'NATURE' 'DIFFUS') ET
  220. (NOMC 'UZ' f2r 'NATURE' 'DIFFUS') ;
  221. mp2r = (mp2 PLUS u2r) COUL 'JAUN' ;
  222. tit1 = CHAI 'IPOL par PID (rouge), par GAUSS (vert) '
  223. 'et par RATIO (jaune)' ;
  224. TRAC (s1b ET mp2p et mp2g et mp2r) 'TITR' tit1 ;
  225. FINSI ;
  226.  
  227.  
  228.  
  229. FIN ;
  230.  
  231.  
  232.  

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