Télécharger inter9.eso

Retour à la liste

Numérotation des lignes :

  1. C INTER9 SOURCE PV 11/03/07 21:17:00 6885
  2. C
  3. C=======================================================================
  4. C
  5. SUBROUTINE INTER9(T0,TP0,ZA0,VOIS2,COEF2,ilent1)
  6. C
  7. C=======================================================================
  8. C
  9. C Calcul de transformations de phases
  10. C appelee par TRPHA2
  11. C
  12. C interpolation lineaire dans un espace 3D
  13. C
  14. C T0,TP0,ZA0 coordonnees du points a interpoler
  15. C VOIS2 points supports de l'interpolation
  16. C COEF2 coef de l'interpolation
  17. C (coordonnees barycentriques)
  18. C Ilent1 suite de nuages ens des pts experimentaux
  19. C
  20. C La routine gere les cas ou
  21. C * le systeme est degenere (ie) les points de VOIS2 generent
  22. C un plan ou une droite
  23. C * le point a interpoler est en dehors du quadrilatere forme
  24. C par les points de VOIS2 : on extrapole
  25. C
  26. C routines appelees
  27. C PROJ projection
  28. C DET3 calcul de determinant
  29. C EQHOM3 resolution d'une equation sans second membre
  30. C EQPAR3 resolution d'une equation avec second membre
  31. C INTE3P calcul de coordonnees barycentriques dans un plan
  32. C
  33. C Michael Martinez 07/98
  34. C
  35. C=======================================================================
  36. C
  37. C
  38. IMPLICIT INTEGER(I-N)
  39. IMPLICIT REAL*8 (A-H,O-Z)
  40. C
  41. -INC SMLENTI
  42. -INC SMLREEL
  43. C
  44. REAL*8 VOIS2(4,3),A(3,3),COEF2(4),X(3),Y(3)
  45. REAL*8 M0(3),M1(3),M2(3),M3(3),M4(3),P0(3)
  46. LOGICAL DIM1,DIM2,DIM3,DETNUL
  47. C-----------------------------------------------------------------------
  48. C------ RECUPERATION DES 4 PLUS PROCHES VOISINS STOCKES DANS VOIS2 -----
  49. C-----------------------------------------------------------------------
  50. C
  51. mlent1=ilent1
  52. segact mlent1
  53. C
  54. IHIST=nint(VOIS2(1,1))
  55. mlenti = mlent1.lect(ihist)
  56. segact mlenti
  57. ITEMP=nint(VOIS2(1,2))
  58. mlreel = lect(ITEMP)
  59. segact mlreel
  60. T1=prog(1)
  61. TP1=prog(2)
  62. ZA1=prog(3)
  63. segdes mlreel,mlenti
  64. C
  65. IHIST=nint(VOIS2(2,1))
  66. mlenti = mlent1.lect(ihist)
  67. segact mlenti
  68. ITEMP=nint(VOIS2(2,2))
  69. mlreel = lect(ITEMP)
  70. segact mlreel
  71. T2=prog(1)
  72. TP2=prog(2)
  73. ZA2=prog(3)
  74. segdes mlreel,mlenti
  75. C
  76. IHIST=nint(VOIS2(3,1))
  77. mlenti = mlent1.lect(ihist)
  78. segact mlenti
  79. ITEMP=nint(VOIS2(3,2))
  80. mlreel = lect(ITEMP)
  81. segact mlreel
  82. T3=prog(1)
  83. TP3=prog(2)
  84. ZA3=prog(3)
  85. segdes mlreel,mlenti
  86. C
  87. IHIST=nint(VOIS2(4,1))
  88. mlenti = mlent1.lect(ihist)
  89. segact mlenti
  90. ITEMP=nint(VOIS2(4,2))
  91. mlreel = lect(ITEMP)
  92. segact mlreel
  93. T4=prog(1)
  94. TP4=prog(2)
  95. ZA4=prog(3)
  96. segdes mlreel,mlenti
  97.  
  98. segdes mlent1
  99. C
  100. M0(1)=T0
  101. M0(2)=TP0
  102. M0(3)=ZA0
  103. M1(1)=T1
  104. M1(2)=TP1
  105. M1(3)=ZA1
  106. M2(1)=T2
  107. M2(2)=TP2
  108. M2(3)=ZA2
  109. M3(1)=T3
  110. M3(2)=TP3
  111. M3(3)=ZA3
  112. M4(1)=T4
  113. M4(2)=TP4
  114. M4(3)=ZA4
  115. C
  116. DIM1=.FALSE.
  117. DIM2=.FALSE.
  118. DIM3=.FALSE.
  119. C
  120. C-----------------------------------------------------------------------
  121. C---------- CAS OU M1,2,3,4 GENERENT L'ESPACE --------------------------
  122. C---------- LES SOLUTIONS LANDA1,2,3,4 SONT UNIQUES --------------------
  123. C---------- A UN FACTEUR MULT. PRES ------------------------------------
  124. C-----------------------------------------------------------------------
  125. C
  126. A(1,1)=T0-T1
  127. A(1,2)=T0-T2
  128. A(1,3)=T0-T3
  129. A(2,1)=TP0-TP1
  130. A(2,2)=TP0-TP2
  131. A(2,3)=TP0-TP3
  132. A(3,1)=ZA0-ZA1
  133. A(3,2)=ZA0-ZA2
  134. A(3,3)=ZA0-ZA3
  135. C
  136. Y(1)=T4-T0
  137. Y(2)=TP4-TP0
  138. Y(3)=ZA4-ZA0
  139. C
  140. C--------------- RESOLUTION DE AX = Y ----------------------------------
  141. C
  142. CALL EQPAR3(X,A,Y,DETNUL)
  143. C
  144. C--------------- DET NON NUL -> M0 N'EST PAS COPLANAIRE AVEC M1,2,3 ----
  145. C
  146. IF (DETNUL.EQV..FALSE.) THEN
  147. SLANDA=X(1)+X(2)+X(3)+1.D0
  148. C
  149. C--------------- SLANDA NON NUL -> M1,2,3,4 GENERENT L'ESPACE ----------
  150. C--------------- DANS CE CAS LE CALCUL S'ARRETE ICI --------------------
  151. C
  152. IF (ABS(SLANDA).GT.1.0D-12) THEN
  153. DIM3=.TRUE.
  154. COEF2(1)=X(1)/SLANDA
  155. COEF2(2)=X(2)/SLANDA
  156. COEF2(3)=X(3)/SLANDA
  157. COEF2(4)=1.D0/SLANDA
  158. C
  159. C--------------- SLANDA NUL -> M1,2,3,4 GENERENT UN PLAN OU UNE DROITE -
  160. C--------------- M0 EST PROJETE SUR CE PLAN (DROITE -> VOIR PLUS LOIN) -
  161. C
  162. ELSE
  163. DIM3=.FALSE.
  164. CALL PROJ(P0,M0,M1,M2,M3,DIM2)
  165. IF (DIM2.EQV..FALSE.) THEN
  166. CALL PROJ(P0,M0,M1,M2,M4,DIM2)
  167. ENDIF
  168. IF (DIM2.EQV..FALSE.) THEN
  169. DIM1=.TRUE.
  170. ELSE
  171. DIM1=.FALSE.
  172. T0=P0(1)
  173. TP0=P0(2)
  174. ZA0=P0(3)
  175. ENDIF
  176. ENDIF
  177. ENDIF
  178. C
  179. C--------------- CAS OU M0 EST COPLANAIRE AVEC M1,2,3 ------------------
  180. C
  181. IF (DIM3.EQV..FALSE..AND.DIM2.EQV..FALSE..AND.DIM1.EQV..FALSE.)
  182. & THEN
  183. C
  184. A(1,1)=T1-T2
  185. A(1,2)=T1-T3
  186. A(1,3)=T1-T4
  187. A(2,1)=TP1-TP2
  188. A(2,2)=TP1-TP3
  189. A(2,3)=TP1-TP4
  190. A(3,1)=ZA1-ZA2
  191. A(3,2)=ZA1-ZA3
  192. A(3,3)=ZA1-ZA4
  193. C
  194. CALL DET3(DETA,A)
  195. C
  196. C--------------- CAS OU M0,1,2,3,4 SONT COPLANAIRES ---------------------
  197. C
  198. IF (ABS(DETA).LE.1.D-12) THEN
  199. DIM2=.TRUE.
  200. C
  201. C--------------- M4 NON COPLANAIRE AVEC M0,1,2,3 ------------------------
  202. C
  203. ELSE
  204. CALL INTE3P(T0,TP0,ZA0,T1,TP1,ZA1,T2,TP2,ZA2,T3,TP3,ZA3,
  205. & COEF2,DIM1,DIM2,DIM3)
  206. COEF2(4)=0.D0
  207. ENDIF
  208. ENDIF
  209. C
  210. C
  211. C IF (DIM3.EQV..FALSE..AND.DIM2.EQV..FALSE..AND.DIM1.EQV..FALSE.)
  212. C & THEN
  213. C DIM1=.TRUE.
  214. C ENDIF
  215. C
  216. C-----------------------------------------------------------------------
  217. C---------- CAS OU M0,1,2,3,4 GENERENT UN PLAN -------------------------
  218. C---------- LES SOLUTIONS COEF2(1,2,3,4) DE NORME MINI SONT CHOISIES ---
  219. C-----------------------------------------------------------------------
  220. C
  221. C
  222. IF (DIM3.EQV..FALSE..AND.DIM1.EQV..FALSE.) THEN
  223. C
  224. C--------------- CALCUL DE LA 1e SOLUTION ------------------------------
  225. C
  226. A(1,1)=T0-T1
  227. A(1,2)=T0-T2
  228. A(1,3)=T0-T3
  229. A(2,1)=TP0-TP1
  230. A(2,2)=TP0-TP2
  231. A(2,3)=TP0-TP3
  232. A(3,1)=ZA0-ZA1
  233. A(3,2)=ZA0-ZA2
  234. A(3,3)=ZA0-ZA3
  235. C
  236. C--------------- RESOLUTION DE AX=0 ------------------------------------
  237. C
  238. CALL EQHOM3 (X,A)
  239. C
  240. COEF11=X(1)
  241. COEF12=X(2)
  242. COEF13=X(3)
  243. COEF14=0.D0
  244. SLAND1=COEF11+COEF12+COEF13+COEF14
  245. C
  246. C--------------- SI M1,2,3 SONT ALIGNES --------------------------------
  247. C
  248. IF (ABS(SLAND1).LT.1.D-12) THEN
  249. A(1,1)=T0-T1
  250. A(1,2)=T0-T4
  251. A(1,3)=T0-T3
  252. A(2,1)=TP0-TP1
  253. A(2,2)=TP0-TP4
  254. A(2,3)=TP0-TP3
  255. A(3,1)=ZA0-ZA1
  256. A(3,2)=ZA0-ZA4
  257. A(3,3)=ZA0-ZA3
  258. C
  259. CALL EQHOM3 (X,A)
  260. C
  261. COEF11=X(1)
  262. COEF12=0.D0
  263. COEF13=X(3)
  264. COEF14=X(2)
  265. SLAND1=COEF11+COEF12+COEF13+COEF14
  266. ENDIF
  267. C
  268. C
  269. C
  270. C--------------- CALCUL DE LA 2e SOLUTION ------------------------------
  271. C
  272. C
  273. A(1,1)=T0-T1
  274. A(1,2)=T0-T2
  275. A(1,3)=T0-T4
  276. A(2,1)=TP0-TP1
  277. A(2,2)=TP0-TP2
  278. A(2,3)=TP0-TP4
  279. A(3,1)=ZA0-ZA1
  280. A(3,2)=ZA0-ZA2
  281. A(3,3)=ZA0-ZA4
  282. C
  283. CALL EQHOM3 (X,A)
  284. C
  285. COEF21=X(1)
  286. COEF22=X(2)
  287. COEF23=0.D0
  288. COEF24=X(3)
  289. SLAND2=COEF21+COEF22+COEF23+COEF24
  290. C
  291. C--------------- SI M1,2,4 SONT ALIGNES --------------------------------
  292. C
  293. IF (ABS(SLAND2).LT.1.D-12) THEN
  294. A(1,1)=T0-T3
  295. A(1,2)=T0-T2
  296. A(1,3)=T0-T4
  297. A(2,1)=TP0-TP3
  298. A(2,2)=TP0-TP2
  299. A(2,3)=TP0-TP4
  300. A(3,1)=ZA0-ZA3
  301. A(3,2)=ZA0-ZA2
  302. A(3,3)=ZA0-ZA4
  303. C
  304. CALL EQHOM3 (X,A)
  305. C
  306. COEF21=0.D0
  307. COEF22=X(2)
  308. COEF23=X(1)
  309. COEF24=X(3)
  310. SLAND2=COEF21+COEF22+COEF23+COEF24
  311. ENDIF
  312. C
  313. C
  314. C
  315. C------ CALCUL DU POINT DE NORME MINI --------------------------------
  316. C
  317. IF (ABS(SLAND1).GE.1.D-12.AND.ABS(SLAND2).GE.1.D-12) THEN
  318. DIM2 = .TRUE.
  319. SLAND1=COEF11+COEF12+COEF13+COEF14
  320. SLAND2=COEF21+COEF22+COEF23+COEF24
  321. C
  322. COEF11=COEF11/SLAND1
  323. COEF12=COEF12/SLAND1
  324. COEF13=COEF13/SLAND1
  325. COEF14=COEF14/SLAND1
  326. COEF21=COEF21/SLAND2
  327. COEF22=COEF22/SLAND2
  328. COEF23=COEF23/SLAND2
  329. COEF24=COEF24/SLAND2
  330. C
  331. SCARR1=COEF11**2+COEF12**2+COEF13**2+COEF14**2
  332. SCARR2=COEF21**2+COEF22**2+COEF23**2+COEF24**2
  333. SDBPRO=COEF11*COEF21+COEF12*COEF22+COEF13*COEF23
  334. . +COEF14*COEF24
  335. C
  336. IF ((SCARR1-2*SDBPRO+SCARR2).GT.1.D-12) THEN
  337. ALPH1=(SCARR2-SDBPRO)/(SCARR1-2*SDBPRO+SCARR2)
  338. BETA1=(SCARR1-SDBPRO)/(SCARR1-2*SDBPRO+SCARR2)
  339. ELSE
  340. IF (SCARR2.GT.SDBPRO) THEN
  341. ALPH1=0.D0
  342. BETA1=1.D0
  343. ELSE
  344. ALPH1=1.D0
  345. BETA1=0.D0
  346. ENDIF
  347. ENDIF
  348. C
  349. COEFF1=ALPH1*COEF11+BETA1*COEF21
  350. COEFF2=ALPH1*COEF12+BETA1*COEF22
  351. COEFF3=ALPH1*COEF13+BETA1*COEF23
  352. COEFF4=ALPH1*COEF14+BETA1*COEF24
  353. C
  354. SLANDA=COEFF1+COEFF2+COEFF3+COEFF4
  355. C
  356. COEF2(1)=COEFF1/SLANDA
  357. COEF2(2)=COEFF2/SLANDA
  358. COEF2(3)=COEFF3/SLANDA
  359. COEF2(4)=COEFF4/SLANDA
  360. C
  361. ELSE
  362. DIM2 = .FALSE.
  363. ENDIF
  364. C
  365. ENDIF
  366. C
  367. C
  368. C-----------------------------------------------------------------------
  369. C---------- CAS OU M0,1,2,3,4 GENERENT UNE DROITE ----------------------
  370. C-----------------------------------------------------------------------
  371. C
  372. IF (DIM2.EQV..FALSE..AND.DIM3.EQV..FALSE.) THEN
  373. C
  374. C > LE PROJETE P0 DE M0 SUR LA DROITE M1,2,3 A POUR COORDONNEES
  375. C > BARYCENTRIQUES (M1,1.-SA2) (M2,SA2) OU SA2 = M1M0.M1M2
  376. C
  377. SA2=(T0-T1)*(T2-T1)+(TP0-TP1)*(TP2-TP1)+(ZA0-ZA1)*(ZA2-ZA1)
  378. C
  379. IF (SA2.GT.0..AND.SA2.LE.1.) THEN
  380. COEF2(1)=1.-SA2
  381. COEF2(2)=SA2
  382. COEF2(3)=0.
  383. COEF2(4)=0.
  384. ELSE
  385. C
  386. C > LE PROJETE P0 DE M0 SUR LA DROITE M1,2,3 A POUR COORDONNEES
  387. C > BARYCENTRIQUES (M1,1.-SA3) (M3,SA3) OU SA3 = M1M0.M1M3
  388. C
  389. SA3=(T0-T1)*(T3-T1)+(TP0-TP1)*(TP3-TP1)+(ZA0-ZA1)*(ZA3-ZA1)
  390. C
  391. IF (SA3.GT.0..AND.SA3.LE.1.) THEN
  392. COEF2(1)=1.-SA3
  393. COEF2(2)=0.
  394. COEF2(3)=SA3
  395. COEF2(4)=0.
  396. ELSE
  397. C
  398. C > LE PROJETE P0 DE M0 SUR LA DROITE M1,2,3 A POUR COORDONNEES
  399. C > BARYCENTRIQUES (M1,1.-SA4) (M4,SA4) OU SA4 = M1M0.M1M4
  400. C
  401. SA4=(T0-T1)*(T4-T1)+(TP0-TP1)*(TP4-TP1)
  402. . +(ZA0-ZA1)*(ZA4-ZA1)
  403. C
  404. COEF2(1)=1.-SA4
  405. COEF2(2)=0.
  406. COEF2(3)=SA4
  407. COEF2(4)=0.
  408. ENDIF
  409. ENDIF
  410. ENDIF
  411. C
  412. C-----------------------------------------------------------------------
  413. C---------- SI L'UNE DES COORDONNEES BARYCENTRIQUES EST NEGATIVE -------
  414. C---------- M0 SE TROUVE A L'EXTERIEUR DU QUADRILATERE M1,2,3,4 --------
  415. C---------- ON LE PROJETE SUR L'UNE DES FACES DU QUADRILATERE ----------
  416. C---------- L'INTERPOLATION SE FAIT ALORS SUR TROIS POINTS -------------
  417. C-----------------------------------------------------------------------
  418. C
  419. C
  420. CMIN1=MIN(COEF2(1),COEF2(2),COEF2(3),COEF2(4))
  421. IF (CMIN1.LT.-0.1D0) THEN
  422. CALL INTE3P(T0,TP0,ZA0,T1,TP1,ZA1,T2,TP2,ZA2,T3,TP3,ZA3,COEF2,
  423. & DIM1,DIM2,DIM3)
  424. COEF2(4)=0.D0
  425. CMIN2=MIN(COEF2(1),COEF2(2),COEF2(3))
  426. C
  427. IF (CMIN2.LT.-0.1D0) THEN
  428. CALL INTE3P(T0,TP0,ZA0,T4,TP4,ZA4,T2,TP2,ZA2,T3,TP3,ZA3,
  429. & COEF2,DIM1,DIM2,DIM3)
  430. COEF2(4)=COEF2(1)
  431. COEF2(1)=0.D0
  432. CMIN3=MIN(COEF2(4),COEF2(2),COEF2(3))
  433. C
  434. IF (CMIN3.LT.-0.1D0) THEN
  435. CALL INTE3P(T0,TP0,ZA0,T1,TP1,ZA1,T4,TP4,ZA4,T3,TP3,ZA3,
  436. & COEF2,DIM1,DIM2,DIM3)
  437. COEF2(4)=COEF2(2)
  438. COEF2(2)=0.D0
  439. CMIN4=MIN(COEF2(1),COEF2(4),COEF2(3))
  440. C
  441. IF (CMIN4.LT.-0.1D0) THEN
  442. CALL INTE3P(T0,TP0,ZA0,T1,TP1,ZA1,T2,TP2,ZA2,T4,TP4,
  443. & ZA4,COEF2,DIM1,DIM2,DIM3)
  444. COEF2(4)=COEF2(3)
  445. COEF2(3)=0.D0
  446. CMIN5=MIN(COEF2(1),COEF2(2),COEF2(4))
  447. C
  448. IF (CMIN5.LT.-0.1D0) THEN
  449. COEF2(1)=1.D0
  450. COEF2(2)=0.D0
  451. COEF2(3)=0.D0
  452. COEF2(4)=0.D0
  453. ENDIF
  454. ENDIF
  455. ENDIF
  456. ENDIF
  457. ENDIF
  458. C
  459. C
  460. RETURN
  461. END
  462.  
  463.  
  464.  
  465.  
  466.  
  467.  

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