Télécharger la plaquette Cast3M


fr en

Télécharger inter9.eso

Retour à la liste

Numérotation des lignes : 

inter9

  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