Télécharger ipol_pid.dgibi

Retour à la liste

Numérotation des lignes :

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

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