C INTER9 SOURCE PV 11/03/07 21:17:00 6885 C C======================================================================= C SUBROUTINE INTER9(T0,TP0,ZA0,VOIS2,COEF2,ilent1) C C======================================================================= C C Calcul de transformations de phases C appelee par TRPHA2 C C interpolation lineaire dans un espace 3D C C T0,TP0,ZA0 coordonnees du points a interpoler C VOIS2 points supports de l'interpolation C COEF2 coef de l'interpolation C (coordonnees barycentriques) C Ilent1 suite de nuages ens des pts experimentaux C C La routine gere les cas ou C * le systeme est degenere (ie) les points de VOIS2 generent C un plan ou une droite C * le point a interpoler est en dehors du quadrilatere forme C par les points de VOIS2 : on extrapole C C routines appelees C PROJ projection C DET3 calcul de determinant C EQHOM3 resolution d'une equation sans second membre C EQPAR3 resolution d'une equation avec second membre C INTE3P calcul de coordonnees barycentriques dans un plan C C Michael Martinez 07/98 C C======================================================================= C C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C -INC SMLENTI -INC SMLREEL C REAL*8 VOIS2(4,3),A(3,3),COEF2(4),X(3),Y(3) REAL*8 M0(3),M1(3),M2(3),M3(3),M4(3),P0(3) LOGICAL DIM1,DIM2,DIM3,DETNUL C----------------------------------------------------------------------- C------ RECUPERATION DES 4 PLUS PROCHES VOISINS STOCKES DANS VOIS2 ----- C----------------------------------------------------------------------- C mlent1=ilent1 segact mlent1 C IHIST=nint(VOIS2(1,1)) mlenti = mlent1.lect(ihist) segact mlenti ITEMP=nint(VOIS2(1,2)) mlreel = lect(ITEMP) segact mlreel T1=prog(1) TP1=prog(2) ZA1=prog(3) segdes mlreel,mlenti C IHIST=nint(VOIS2(2,1)) mlenti = mlent1.lect(ihist) segact mlenti ITEMP=nint(VOIS2(2,2)) mlreel = lect(ITEMP) segact mlreel T2=prog(1) TP2=prog(2) ZA2=prog(3) segdes mlreel,mlenti C IHIST=nint(VOIS2(3,1)) mlenti = mlent1.lect(ihist) segact mlenti ITEMP=nint(VOIS2(3,2)) mlreel = lect(ITEMP) segact mlreel T3=prog(1) TP3=prog(2) ZA3=prog(3) segdes mlreel,mlenti C IHIST=nint(VOIS2(4,1)) mlenti = mlent1.lect(ihist) segact mlenti ITEMP=nint(VOIS2(4,2)) mlreel = lect(ITEMP) segact mlreel T4=prog(1) TP4=prog(2) ZA4=prog(3) segdes mlreel,mlenti segdes mlent1 C M0(1)=T0 M0(2)=TP0 M0(3)=ZA0 M1(1)=T1 M1(2)=TP1 M1(3)=ZA1 M2(1)=T2 M2(2)=TP2 M2(3)=ZA2 M3(1)=T3 M3(2)=TP3 M3(3)=ZA3 M4(1)=T4 M4(2)=TP4 M4(3)=ZA4 C DIM1=.FALSE. DIM2=.FALSE. DIM3=.FALSE. C C----------------------------------------------------------------------- C---------- CAS OU M1,2,3,4 GENERENT L'ESPACE -------------------------- C---------- LES SOLUTIONS LANDA1,2,3,4 SONT UNIQUES -------------------- C---------- A UN FACTEUR MULT. PRES ------------------------------------ C----------------------------------------------------------------------- C A(1,1)=T0-T1 A(1,2)=T0-T2 A(1,3)=T0-T3 A(2,1)=TP0-TP1 A(2,2)=TP0-TP2 A(2,3)=TP0-TP3 A(3,1)=ZA0-ZA1 A(3,2)=ZA0-ZA2 A(3,3)=ZA0-ZA3 C Y(1)=T4-T0 Y(2)=TP4-TP0 Y(3)=ZA4-ZA0 C C--------------- RESOLUTION DE AX = Y ---------------------------------- C CALL EQPAR3(X,A,Y,DETNUL) C C--------------- DET NON NUL -> M0 N'EST PAS COPLANAIRE AVEC M1,2,3 ---- C IF (DETNUL.EQV..FALSE.) THEN SLANDA=X(1)+X(2)+X(3)+1.D0 C C--------------- SLANDA NON NUL -> M1,2,3,4 GENERENT L'ESPACE ---------- C--------------- DANS CE CAS LE CALCUL S'ARRETE ICI -------------------- C IF (ABS(SLANDA).GT.1.0D-12) THEN DIM3=.TRUE. COEF2(1)=X(1)/SLANDA COEF2(2)=X(2)/SLANDA COEF2(3)=X(3)/SLANDA COEF2(4)=1.D0/SLANDA C C--------------- SLANDA NUL -> M1,2,3,4 GENERENT UN PLAN OU UNE DROITE - C--------------- M0 EST PROJETE SUR CE PLAN (DROITE -> VOIR PLUS LOIN) - C ELSE DIM3=.FALSE. CALL PROJ(P0,M0,M1,M2,M3,DIM2) IF (DIM2.EQV..FALSE.) THEN CALL PROJ(P0,M0,M1,M2,M4,DIM2) ENDIF IF (DIM2.EQV..FALSE.) THEN DIM1=.TRUE. ELSE DIM1=.FALSE. T0=P0(1) TP0=P0(2) ZA0=P0(3) ENDIF ENDIF ENDIF C C--------------- CAS OU M0 EST COPLANAIRE AVEC M1,2,3 ------------------ C IF (DIM3.EQV..FALSE..AND.DIM2.EQV..FALSE..AND.DIM1.EQV..FALSE.) & THEN C A(1,1)=T1-T2 A(1,2)=T1-T3 A(1,3)=T1-T4 A(2,1)=TP1-TP2 A(2,2)=TP1-TP3 A(2,3)=TP1-TP4 A(3,1)=ZA1-ZA2 A(3,2)=ZA1-ZA3 A(3,3)=ZA1-ZA4 C CALL DET3(DETA,A) C C--------------- CAS OU M0,1,2,3,4 SONT COPLANAIRES --------------------- C IF (ABS(DETA).LE.1.D-12) THEN DIM2=.TRUE. C C--------------- M4 NON COPLANAIRE AVEC M0,1,2,3 ------------------------ C ELSE CALL INTE3P(T0,TP0,ZA0,T1,TP1,ZA1,T2,TP2,ZA2,T3,TP3,ZA3, & COEF2,DIM1,DIM2,DIM3) COEF2(4)=0.D0 ENDIF ENDIF C C C IF (DIM3.EQV..FALSE..AND.DIM2.EQV..FALSE..AND.DIM1.EQV..FALSE.) C & THEN C DIM1=.TRUE. C ENDIF C C----------------------------------------------------------------------- C---------- CAS OU M0,1,2,3,4 GENERENT UN PLAN ------------------------- C---------- LES SOLUTIONS COEF2(1,2,3,4) DE NORME MINI SONT CHOISIES --- C----------------------------------------------------------------------- C C IF (DIM3.EQV..FALSE..AND.DIM1.EQV..FALSE.) THEN C C--------------- CALCUL DE LA 1e SOLUTION ------------------------------ C A(1,1)=T0-T1 A(1,2)=T0-T2 A(1,3)=T0-T3 A(2,1)=TP0-TP1 A(2,2)=TP0-TP2 A(2,3)=TP0-TP3 A(3,1)=ZA0-ZA1 A(3,2)=ZA0-ZA2 A(3,3)=ZA0-ZA3 C C--------------- RESOLUTION DE AX=0 ------------------------------------ C CALL EQHOM3 (X,A) C COEF11=X(1) COEF12=X(2) COEF13=X(3) COEF14=0.D0 SLAND1=COEF11+COEF12+COEF13+COEF14 C C--------------- SI M1,2,3 SONT ALIGNES -------------------------------- C IF (ABS(SLAND1).LT.1.D-12) THEN A(1,1)=T0-T1 A(1,2)=T0-T4 A(1,3)=T0-T3 A(2,1)=TP0-TP1 A(2,2)=TP0-TP4 A(2,3)=TP0-TP3 A(3,1)=ZA0-ZA1 A(3,2)=ZA0-ZA4 A(3,3)=ZA0-ZA3 C CALL EQHOM3 (X,A) C COEF11=X(1) COEF12=0.D0 COEF13=X(3) COEF14=X(2) SLAND1=COEF11+COEF12+COEF13+COEF14 ENDIF C C C C--------------- CALCUL DE LA 2e SOLUTION ------------------------------ C C A(1,1)=T0-T1 A(1,2)=T0-T2 A(1,3)=T0-T4 A(2,1)=TP0-TP1 A(2,2)=TP0-TP2 A(2,3)=TP0-TP4 A(3,1)=ZA0-ZA1 A(3,2)=ZA0-ZA2 A(3,3)=ZA0-ZA4 C CALL EQHOM3 (X,A) C COEF21=X(1) COEF22=X(2) COEF23=0.D0 COEF24=X(3) SLAND2=COEF21+COEF22+COEF23+COEF24 C C--------------- SI M1,2,4 SONT ALIGNES -------------------------------- C IF (ABS(SLAND2).LT.1.D-12) THEN A(1,1)=T0-T3 A(1,2)=T0-T2 A(1,3)=T0-T4 A(2,1)=TP0-TP3 A(2,2)=TP0-TP2 A(2,3)=TP0-TP4 A(3,1)=ZA0-ZA3 A(3,2)=ZA0-ZA2 A(3,3)=ZA0-ZA4 C CALL EQHOM3 (X,A) C COEF21=0.D0 COEF22=X(2) COEF23=X(1) COEF24=X(3) SLAND2=COEF21+COEF22+COEF23+COEF24 ENDIF C C C C------ CALCUL DU POINT DE NORME MINI -------------------------------- C IF (ABS(SLAND1).GE.1.D-12.AND.ABS(SLAND2).GE.1.D-12) THEN DIM2 = .TRUE. SLAND1=COEF11+COEF12+COEF13+COEF14 SLAND2=COEF21+COEF22+COEF23+COEF24 C COEF11=COEF11/SLAND1 COEF12=COEF12/SLAND1 COEF13=COEF13/SLAND1 COEF14=COEF14/SLAND1 COEF21=COEF21/SLAND2 COEF22=COEF22/SLAND2 COEF23=COEF23/SLAND2 COEF24=COEF24/SLAND2 C SCARR1=COEF11**2+COEF12**2+COEF13**2+COEF14**2 SCARR2=COEF21**2+COEF22**2+COEF23**2+COEF24**2 SDBPRO=COEF11*COEF21+COEF12*COEF22+COEF13*COEF23 . +COEF14*COEF24 C IF ((SCARR1-2*SDBPRO+SCARR2).GT.1.D-12) THEN ALPH1=(SCARR2-SDBPRO)/(SCARR1-2*SDBPRO+SCARR2) BETA1=(SCARR1-SDBPRO)/(SCARR1-2*SDBPRO+SCARR2) ELSE IF (SCARR2.GT.SDBPRO) THEN ALPH1=0.D0 BETA1=1.D0 ELSE ALPH1=1.D0 BETA1=0.D0 ENDIF ENDIF C COEFF1=ALPH1*COEF11+BETA1*COEF21 COEFF2=ALPH1*COEF12+BETA1*COEF22 COEFF3=ALPH1*COEF13+BETA1*COEF23 COEFF4=ALPH1*COEF14+BETA1*COEF24 C SLANDA=COEFF1+COEFF2+COEFF3+COEFF4 C COEF2(1)=COEFF1/SLANDA COEF2(2)=COEFF2/SLANDA COEF2(3)=COEFF3/SLANDA COEF2(4)=COEFF4/SLANDA C ELSE DIM2 = .FALSE. ENDIF C ENDIF C C C----------------------------------------------------------------------- C---------- CAS OU M0,1,2,3,4 GENERENT UNE DROITE ---------------------- C----------------------------------------------------------------------- C IF (DIM2.EQV..FALSE..AND.DIM3.EQV..FALSE.) THEN C C > LE PROJETE P0 DE M0 SUR LA DROITE M1,2,3 A POUR COORDONNEES C > BARYCENTRIQUES (M1,1.-SA2) (M2,SA2) OU SA2 = M1M0.M1M2 C SA2=(T0-T1)*(T2-T1)+(TP0-TP1)*(TP2-TP1)+(ZA0-ZA1)*(ZA2-ZA1) C IF (SA2.GT.0..AND.SA2.LE.1.) THEN COEF2(1)=1.-SA2 COEF2(2)=SA2 COEF2(3)=0. COEF2(4)=0. ELSE C C > LE PROJETE P0 DE M0 SUR LA DROITE M1,2,3 A POUR COORDONNEES C > BARYCENTRIQUES (M1,1.-SA3) (M3,SA3) OU SA3 = M1M0.M1M3 C SA3=(T0-T1)*(T3-T1)+(TP0-TP1)*(TP3-TP1)+(ZA0-ZA1)*(ZA3-ZA1) C IF (SA3.GT.0..AND.SA3.LE.1.) THEN COEF2(1)=1.-SA3 COEF2(2)=0. COEF2(3)=SA3 COEF2(4)=0. ELSE C C > LE PROJETE P0 DE M0 SUR LA DROITE M1,2,3 A POUR COORDONNEES C > BARYCENTRIQUES (M1,1.-SA4) (M4,SA4) OU SA4 = M1M0.M1M4 C SA4=(T0-T1)*(T4-T1)+(TP0-TP1)*(TP4-TP1) . +(ZA0-ZA1)*(ZA4-ZA1) C COEF2(1)=1.-SA4 COEF2(2)=0. COEF2(3)=SA4 COEF2(4)=0. ENDIF ENDIF ENDIF C C----------------------------------------------------------------------- C---------- SI L'UNE DES COORDONNEES BARYCENTRIQUES EST NEGATIVE ------- C---------- M0 SE TROUVE A L'EXTERIEUR DU QUADRILATERE M1,2,3,4 -------- C---------- ON LE PROJETE SUR L'UNE DES FACES DU QUADRILATERE ---------- C---------- L'INTERPOLATION SE FAIT ALORS SUR TROIS POINTS ------------- C----------------------------------------------------------------------- C C CMIN1=MIN(COEF2(1),COEF2(2),COEF2(3),COEF2(4)) IF (CMIN1.LT.-0.1D0) THEN CALL INTE3P(T0,TP0,ZA0,T1,TP1,ZA1,T2,TP2,ZA2,T3,TP3,ZA3,COEF2, & DIM1,DIM2,DIM3) COEF2(4)=0.D0 CMIN2=MIN(COEF2(1),COEF2(2),COEF2(3)) C IF (CMIN2.LT.-0.1D0) THEN CALL INTE3P(T0,TP0,ZA0,T4,TP4,ZA4,T2,TP2,ZA2,T3,TP3,ZA3, & COEF2,DIM1,DIM2,DIM3) COEF2(4)=COEF2(1) COEF2(1)=0.D0 CMIN3=MIN(COEF2(4),COEF2(2),COEF2(3)) C IF (CMIN3.LT.-0.1D0) THEN CALL INTE3P(T0,TP0,ZA0,T1,TP1,ZA1,T4,TP4,ZA4,T3,TP3,ZA3, & COEF2,DIM1,DIM2,DIM3) COEF2(4)=COEF2(2) COEF2(2)=0.D0 CMIN4=MIN(COEF2(1),COEF2(4),COEF2(3)) C IF (CMIN4.LT.-0.1D0) THEN CALL INTE3P(T0,TP0,ZA0,T1,TP1,ZA1,T2,TP2,ZA2,T4,TP4, & ZA4,COEF2,DIM1,DIM2,DIM3) COEF2(4)=COEF2(3) COEF2(3)=0.D0 CMIN5=MIN(COEF2(1),COEF2(2),COEF2(4)) C IF (CMIN5.LT.-0.1D0) THEN COEF2(1)=1.D0 COEF2(2)=0.D0 COEF2(3)=0.D0 COEF2(4)=0.D0 ENDIF ENDIF ENDIF ENDIF ENDIF C C RETURN END