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 141 : CALL EQPAR3(X,A,Y,DETNUL) 142 : C 143 : C--------------- DET NON NUL -> M0 N'EST PAS COPLANAIRE AVEC M1,2,3 ---- 144 : C 145 : IF (DETNUL.EQV..FALSE.) THEN 146 : SLANDA=X(1)+X(2)+X(3)+1.D0 147 : C 148 : C--------------- SLANDA NON NUL -> M1,2,3,4 GENERENT L'ESPACE ---------- 149 : C--------------- DANS CE CAS LE CALCUL S'ARRETE ICI -------------------- 150 : C 151 : IF (ABS(SLANDA).GT.1.0D-12) THEN 152 : DIM3=.TRUE. 153 : COEF2(1)=X(1)/SLANDA 154 : COEF2(2)=X(2)/SLANDA 155 : COEF2(3)=X(3)/SLANDA 156 : COEF2(4)=1.D0/SLANDA 157 : C 158 : C--------------- SLANDA NUL -> M1,2,3,4 GENERENT UN PLAN OU UNE DROITE - 159 : C--------------- M0 EST PROJETE SUR CE PLAN (DROITE -> VOIR PLUS LOIN) - 160 : C 161 : ELSE 162 : DIM3=.FALSE. 162 : CALL PROJ(P0,M0,M1,M2,M3,DIM2) 163 : IF (DIM2.EQV..FALSE.) THEN 163 : CALL PROJ(P0,M0,M1,M2,M4,DIM2) 164 : ENDIF 165 : IF (DIM2.EQV..FALSE.) THEN 166 : DIM1=.TRUE. 167 : ELSE 168 : DIM1=.FALSE. 169 : T0=P0(1) 170 : TP0=P0(2) 171 : ZA0=P0(3) 172 : ENDIF 173 : ENDIF 174 : ENDIF 175 : C 176 : C--------------- CAS OU M0 EST COPLANAIRE AVEC M1,2,3 ------------------ 177 : C 178 : IF (DIM3.EQV..FALSE..AND.DIM2.EQV..FALSE..AND.DIM1.EQV..FALSE.) 179 : & THEN 180 : C 181 : A(1,1)=T1-T2 182 : A(1,2)=T1-T3 183 : A(1,3)=T1-T4 184 : A(2,1)=TP1-TP2 185 : A(2,2)=TP1-TP3 186 : A(2,3)=TP1-TP4 187 : A(3,1)=ZA1-ZA2 188 : A(3,2)=ZA1-ZA3 189 : A(3,3)=ZA1-ZA4 190 : C 190 : CALL DET3(DETA,A) 191 : C 192 : C--------------- CAS OU M0,1,2,3,4 SONT COPLANAIRES --------------------- 193 : C 194 : IF (ABS(DETA).LE.1.D-12) THEN 195 : DIM2=.TRUE. 196 : C 197 : C--------------- M4 NON COPLANAIRE AVEC M0,1,2,3 ------------------------ 198 : C 199 : ELSE 199 : CALL INTE3P(T0,TP0,ZA0,T1,TP1,ZA1,T2,TP2,ZA2,T3,TP3,ZA3, 200 : & COEF2,DIM1,DIM2,DIM3) 201 : COEF2(4)=0.D0 202 : ENDIF 203 : ENDIF 204 : C 205 : C 206 : C IF (DIM3.EQV..FALSE..AND.DIM2.EQV..FALSE..AND.DIM1.EQV..FALSE.) 207 : C & THEN 208 : C DIM1=.TRUE. 209 : C ENDIF 210 : C 211 : C----------------------------------------------------------------------- 212 : C---------- CAS OU M0,1,2,3,4 GENERENT UN PLAN ------------------------- 213 : C---------- LES SOLUTIONS COEF2(1,2,3,4) DE NORME MINI SONT CHOISIES --- 214 : C----------------------------------------------------------------------- 215 : C 216 : C 217 : IF (DIM3.EQV..FALSE..AND.DIM1.EQV..FALSE.) THEN 218 : C 219 : C--------------- CALCUL DE LA 1e SOLUTION ------------------------------ 220 : C 221 : A(1,1)=T0-T1 222 : A(1,2)=T0-T2 223 : A(1,3)=T0-T3 224 : A(2,1)=TP0-TP1 225 : A(2,2)=TP0-TP2 226 : A(2,3)=TP0-TP3 227 : A(3,1)=ZA0-ZA1 228 : A(3,2)=ZA0-ZA2 229 : A(3,3)=ZA0-ZA3 230 : C 231 : C--------------- RESOLUTION DE AX=0 ------------------------------------ 232 : C 232 : CALL EQHOM3 (X,A) 233 : C 234 : COEF11=X(1) 235 : COEF12=X(2) 236 : COEF13=X(3) 237 : COEF14=0.D0 238 : SLAND1=COEF11+COEF12+COEF13+COEF14 239 : C 240 : C--------------- SI M1,2,3 SONT ALIGNES -------------------------------- 241 : C 242 : IF (ABS(SLAND1).LT.1.D-12) THEN 243 : A(1,1)=T0-T1 244 : A(1,2)=T0-T4 245 : A(1,3)=T0-T3 246 : A(2,1)=TP0-TP1 247 : A(2,2)=TP0-TP4 248 : A(2,3)=TP0-TP3 249 : A(3,1)=ZA0-ZA1 250 : A(3,2)=ZA0-ZA4 251 : A(3,3)=ZA0-ZA3 252 : C 252 : CALL EQHOM3 (X,A) 253 : C 254 : COEF11=X(1) 255 : COEF12=0.D0 256 : COEF13=X(3) 257 : COEF14=X(2) 258 : SLAND1=COEF11+COEF12+COEF13+COEF14 259 : ENDIF 260 : C 261 : C 262 : C 263 : C--------------- CALCUL DE LA 2e SOLUTION ------------------------------ 264 : C 265 : C 266 : A(1,1)=T0-T1 267 : A(1,2)=T0-T2 268 : A(1,3)=T0-T4 269 : A(2,1)=TP0-TP1 270 : A(2,2)=TP0-TP2 271 : A(2,3)=TP0-TP4 272 : A(3,1)=ZA0-ZA1 273 : A(3,2)=ZA0-ZA2 274 : A(3,3)=ZA0-ZA4 275 : C 275 : CALL EQHOM3 (X,A) 276 : C 277 : COEF21=X(1) 278 : COEF22=X(2) 279 : COEF23=0.D0 280 : COEF24=X(3) 281 : SLAND2=COEF21+COEF22+COEF23+COEF24 282 : C 283 : C--------------- SI M1,2,4 SONT ALIGNES -------------------------------- 284 : C 285 : IF (ABS(SLAND2).LT.1.D-12) THEN 286 : A(1,1)=T0-T3 287 : A(1,2)=T0-T2 288 : A(1,3)=T0-T4 289 : A(2,1)=TP0-TP3 290 : A(2,2)=TP0-TP2 291 : A(2,3)=TP0-TP4 292 : A(3,1)=ZA0-ZA3 293 : A(3,2)=ZA0-ZA2 294 : A(3,3)=ZA0-ZA4 295 : C 295 : CALL EQHOM3 (X,A) 296 : C 297 : COEF21=0.D0 298 : COEF22=X(2) 299 : COEF23=X(1) 300 : COEF24=X(3) 301 : SLAND2=COEF21+COEF22+COEF23+COEF24 302 : ENDIF 303 : C 304 : C 305 : C 306 : C------ CALCUL DU POINT DE NORME MINI -------------------------------- 307 : C 308 : IF (ABS(SLAND1).GE.1.D-12.AND.ABS(SLAND2).GE.1.D-12) THEN 309 : DIM2 = .TRUE. 310 : SLAND1=COEF11+COEF12+COEF13+COEF14 311 : SLAND2=COEF21+COEF22+COEF23+COEF24 312 : C 313 : COEF11=COEF11/SLAND1 314 : COEF12=COEF12/SLAND1 315 : COEF13=COEF13/SLAND1 316 : COEF14=COEF14/SLAND1 317 : COEF21=COEF21/SLAND2 318 : COEF22=COEF22/SLAND2 319 : COEF23=COEF23/SLAND2 320 : COEF24=COEF24/SLAND2 321 : C 322 : SCARR1=COEF11**2+COEF12**2+COEF13**2+COEF14**2 323 : SCARR2=COEF21**2+COEF22**2+COEF23**2+COEF24**2 324 : SDBPRO=COEF11*COEF21+COEF12*COEF22+COEF13*COEF23 325 : . +COEF14*COEF24 326 : C 327 : IF ((SCARR1-2*SDBPRO+SCARR2).GT.1.D-12) THEN 328 : ALPH1=(SCARR2-SDBPRO)/(SCARR1-2*SDBPRO+SCARR2) 329 : BETA1=(SCARR1-SDBPRO)/(SCARR1-2*SDBPRO+SCARR2) 330 : ELSE 331 : IF (SCARR2.GT.SDBPRO) THEN 332 : ALPH1=0.D0 333 : BETA1=1.D0 334 : ELSE 335 : ALPH1=1.D0 336 : BETA1=0.D0 337 : ENDIF 338 : ENDIF 339 : C 340 : COEFF1=ALPH1*COEF11+BETA1*COEF21 341 : COEFF2=ALPH1*COEF12+BETA1*COEF22 342 : COEFF3=ALPH1*COEF13+BETA1*COEF23 343 : COEFF4=ALPH1*COEF14+BETA1*COEF24 344 : C 345 : SLANDA=COEFF1+COEFF2+COEFF3+COEFF4 346 : C 347 : COEF2(1)=COEFF1/SLANDA 348 : COEF2(2)=COEFF2/SLANDA 349 : COEF2(3)=COEFF3/SLANDA 350 : COEF2(4)=COEFF4/SLANDA 351 : C 352 : ELSE 353 : DIM2 = .FALSE. 354 : ENDIF 355 : C 356 : ENDIF 357 : C 358 : C 359 : C----------------------------------------------------------------------- 360 : C---------- CAS OU M0,1,2,3,4 GENERENT UNE DROITE ---------------------- 361 : C----------------------------------------------------------------------- 362 : C 363 : IF (DIM2.EQV..FALSE..AND.DIM3.EQV..FALSE.) THEN 364 : C 365 : C > LE PROJETE P0 DE M0 SUR LA DROITE M1,2,3 A POUR COORDONNEES 366 : C > BARYCENTRIQUES (M1,1.-SA2) (M2,SA2) OU SA2 = M1M0.M1M2 367 : C 368 : SA2=(T0-T1)*(T2-T1)+(TP0-TP1)*(TP2-TP1)+(ZA0-ZA1)*(ZA2-ZA1) 369 : C 370 : IF (SA2.GT.0..AND.SA2.LE.1.) THEN 371 : COEF2(1)=1.-SA2 372 : COEF2(2)=SA2 373 : COEF2(3)=0. 374 : COEF2(4)=0. 375 : ELSE 376 : C 377 : C > LE PROJETE P0 DE M0 SUR LA DROITE M1,2,3 A POUR COORDONNEES 378 : C > BARYCENTRIQUES (M1,1.-SA3) (M3,SA3) OU SA3 = M1M0.M1M3 379 : C 380 : SA3=(T0-T1)*(T3-T1)+(TP0-TP1)*(TP3-TP1)+(ZA0-ZA1)*(ZA3-ZA1) 381 : C 382 : IF (SA3.GT.0..AND.SA3.LE.1.) THEN 383 : COEF2(1)=1.-SA3 384 : COEF2(2)=0. 385 : COEF2(3)=SA3 386 : COEF2(4)=0. 387 : ELSE 388 : C 389 : C > LE PROJETE P0 DE M0 SUR LA DROITE M1,2,3 A POUR COORDONNEES 390 : C > BARYCENTRIQUES (M1,1.-SA4) (M4,SA4) OU SA4 = M1M0.M1M4 391 : C 392 : SA4=(T0-T1)*(T4-T1)+(TP0-TP1)*(TP4-TP1) 393 : . +(ZA0-ZA1)*(ZA4-ZA1) 394 : C 395 : COEF2(1)=1.-SA4 396 : COEF2(2)=0. 397 : COEF2(3)=SA4 398 : COEF2(4)=0. 399 : ENDIF 400 : ENDIF 401 : ENDIF 402 : C 403 : C----------------------------------------------------------------------- 404 : C---------- SI L'UNE DES COORDONNEES BARYCENTRIQUES EST NEGATIVE ------- 405 : C---------- M0 SE TROUVE A L'EXTERIEUR DU QUADRILATERE M1,2,3,4 -------- 406 : C---------- ON LE PROJETE SUR L'UNE DES FACES DU QUADRILATERE ---------- 407 : C---------- L'INTERPOLATION SE FAIT ALORS SUR TROIS POINTS ------------- 408 : C----------------------------------------------------------------------- 409 : C 410 : C 411 : CMIN1=MIN(COEF2(1),COEF2(2),COEF2(3),COEF2(4)) 412 : IF (CMIN1.LT.-0.1D0) THEN 412 : CALL INTE3P(T0,TP0,ZA0,T1,TP1,ZA1,T2,TP2,ZA2,T3,TP3,ZA3,COEF2, 413 : & DIM1,DIM2,DIM3) 414 : COEF2(4)=0.D0 415 : CMIN2=MIN(COEF2(1),COEF2(2),COEF2(3)) 416 : C 417 : IF (CMIN2.LT.-0.1D0) THEN 417 : CALL INTE3P(T0,TP0,ZA0,T4,TP4,ZA4,T2,TP2,ZA2,T3,TP3,ZA3, 418 : & COEF2,DIM1,DIM2,DIM3) 419 : COEF2(4)=COEF2(1) 420 : COEF2(1)=0.D0 421 : CMIN3=MIN(COEF2(4),COEF2(2),COEF2(3)) 422 : C 423 : IF (CMIN3.LT.-0.1D0) THEN 423 : CALL INTE3P(T0,TP0,ZA0,T1,TP1,ZA1,T4,TP4,ZA4,T3,TP3,ZA3, 424 : & COEF2,DIM1,DIM2,DIM3) 425 : COEF2(4)=COEF2(2) 426 : COEF2(2)=0.D0 427 : CMIN4=MIN(COEF2(1),COEF2(4),COEF2(3)) 428 : C 429 : IF (CMIN4.LT.-0.1D0) THEN 429 : CALL INTE3P(T0,TP0,ZA0,T1,TP1,ZA1,T2,TP2,ZA2,T4,TP4, 430 : & ZA4,COEF2,DIM1,DIM2,DIM3) 431 : COEF2(4)=COEF2(3) 432 : COEF2(3)=0.D0 433 : CMIN5=MIN(COEF2(1),COEF2(2),COEF2(4)) 434 : C 435 : IF (CMIN5.LT.-0.1D0) THEN 436 : COEF2(1)=1.D0 437 : COEF2(2)=0.D0 438 : COEF2(3)=0.D0 439 : COEF2(4)=0.D0 440 : ENDIF 441 : ENDIF 442 : ENDIF 443 : ENDIF 444 : ENDIF 445 : C 446 : C 447 : RETURN 448 : END 449 : 450 : 451 : 452 : 453 :
© Cast3M 2003 - All rights reserved.
Disclaimer