Télécharger la plaquette Cast3M


fr en

Télécharger inters.eso

Retour à la liste

Numérotation des lignes : 

inters

  1. C INTERS    SOURCE    PASCAL    21/07/15    21:15:01     11076          
  2.       SUBROUTINE INTERS
  3. C-----------------------------------------------------------------------
  4. C
  5. C     INTERSECTIONS DE SURFACES (PLAN,SPHERE,CYLINDRE,CONE,TORE)
  6. C
  7. C-----------------------------------------------------------------------
  8. C
  9. C     IOB1   EVENTUELLEMENT  OBJET QUE L'ON PROLONGE
  10. C     IOB2   EVENTUELLEMENT  OBJET QUE L'ON AVANCE
  11. C     IPOIN1 POINT INITIAL
  12. C     IPOIN2 POINT FINAL
  13. C     ICLE   TYPE DE LA PREMIERE SURFACE
  14. C     JCLE   TYPE DE LA DEUXIEME SURFACE
  15. C     XPAUX  TABLEAU DE POINTS AUXILIAIRES SERVANT A RECTIFIER L'INTERSE
  16. C             LA QUATRIEME COMPOSANTE EST L'ABCISSE CURVILIGNE
  17. C
  18. C-----------------------------------------------------------------------
  19. C
  20. C     17/09/93 : intersection de deux MELEMEs (sens ensembliste)
  21. C
  22. C  PM 30/08/05 : Erreur si intersection vide, sauf si mot-clef 'NOVERIF'
  23. C
  24. C  BP 28/08/2012 : Ajout de l'intersection geometrique de 2 maillages
  25. C
  26. C-----------------------------------------------------------------------
  27. c
  28.       IMPLICIT INTEGER(I-N)
  29.       IMPLICIT REAL*8 (A-H,O-Z)
  30.  
  31. -INC PPARAM
  32. -INC CCOPTIO
  33. -INC CCGEOME
  34.  
  35. -INC SMCOORD
  36. -INC SMELEME
  37.  
  38.       SEGMENT TABPAR(NPOIN)
  39.  
  40.       LOGICAL VERIF,ltelq
  41.       CHARACTER*4 MCLE(8),MCLE2(1),MCLE3(1)
  42. C     INITIALISATION A L'USAGE DE L'OPTIMISEUR
  43.       PARAMETER ( NBT = 22 )
  44.       DIMENSION XPAUX(4,NBT)
  45.  
  46.       DATA RAYI/1./
  47.       DATA MCLE/'PLAN','SPHE','CYLI','CONI','TORI','DINI','DFIN','LONG'/
  48.       DATA MCLE2/'NOVE'/
  49.       DATA MCLE3/'GEOM'/
  50. C
  51. C-----------------------------------------------------------------------
  52. C     INTERSECTION DE DEUX LISTREELS
  53. C     Dans ce cas, le resultat est directment ecrit dans la pile.
  54.  
  55.       CALL LIROBJ('LISTREEL',IPROG1,0,IRETOU)
  56.       IF (IERR.NE.0) RETURN
  57.  
  58.       IF (IRETOU.EQ.1) THEN
  59.         CALL LIROBJ('LISTREEL',IPROG2,1,IRETOU)
  60.         IF (IERR.NE.0) RETURN
  61.  
  62.         CALL INTLRE(IPROG1,IPROG2)
  63.         RETURN
  64.  
  65.       ENDIF
  66. C
  67. C-----------------------------------------------------------------------
  68. C     INTERSECTION DE DEUX RIGIDITES
  69. C     Dans ce cas, le resultat est directment ecrit dans la pile.
  70.  
  71.       CALL LIROBJ('RIGIDITE',IPRIG1,0,IRETOU)
  72.       IF (IERR.NE.0) RETURN
  73.  
  74.       IF (IRETOU.EQ.1) THEN
  75.         CALL LIROBJ('RIGIDITE',IPRIG2,1,IRETOU)
  76.         IF (IERR.NE.0) RETURN
  77.  
  78.         CALL INTRIG(IPRIG1,IPRIG2)
  79.         RETURN
  80.  
  81.       ENDIF
  82. C
  83. C-----------------------------------------------------------------------
  84. C     INTERSECTION DE DEUX MODELES
  85. C     Dans ce cas, le resultat est directment ecrit dans la pile.
  86.  
  87.       CALL LIROBJ('MMODEL',IPMOD1,0,IRETOU)
  88.       IF (IERR.NE.0) RETURN
  89.  
  90.       IF (IRETOU.EQ.1) THEN
  91.         CALL LIROBJ('MMODEL',IPMOD2,1,IRETOU)
  92.         IF (IERR.NE.0) RETURN
  93.  
  94.         CALL INTMOD(IPMOD1,IPMOD2)
  95.         RETURN
  96.  
  97.       ENDIF
  98. C
  99. C-----------------------------------------------------------------------
  100. C     INTERSECTION DE DEUX MAILLAGES
  101. C
  102. C     sens ensembliste ou geometrique?
  103.       CALL LIRMOT(MCLE3,1,IGEO3,0)
  104. c       write(ioimp,*) 'option GEOM ? ',IGEO3
  105. C
  106. c     lecture des maillages
  107.       CALL LIROBJ('MAILLAGE',IOB1,0,IRETOU)
  108.       IF (IRETOU.EQ.1) THEN
  109.         CALL LIROBJ('MAILLAGE',IOB2,0,IRETOV)
  110. C       On laisse la syntaxe NOVERIF pour l'intersection Géométrique
  111.         CALL LIRMOT(MCLE2,1,ICLE,0)
  112.         VERIF=(ICLE.EQ.0)
  113.  
  114. C------ Intersection ensembliste ------------------
  115.         if(IGEO3.eq.0) then
  116.           IF (IRETOV.EQ.1) THEN
  117. C           Intersection de deux maillages
  118.             CALL INTERB(IOB1,IOB2,IRETIB,IPT3)
  119.             IF (IRETIB.EQ.1) THEN
  120. C             On fournit un maillage vide
  121.               NBNN=0
  122.               NBELEM=0
  123.               NBSOUS=0
  124.               NBREF=0
  125.               SEGINI,IPT3
  126.               IPT3.ITYPEL=0
  127.               SEGDES,IPT3
  128.             ENDIF
  129.             CALL ECROBJ('MAILLAGE',IPT3)
  130.             RETURN
  131.           ELSE
  132.             CALL REFUS
  133.           ENDIF
  134.  
  135. C------ Intersection geometrique ------------------
  136.         else
  137.           IOB3=0
  138.           IOB4=0
  139.           CALL INTGEO(IOB1,IOB2,IOB3,IOB4,ICLE)
  140. C         Cas ou n'a pas pu calculer l'intersection
  141.           if(IOB3.eq.0) then
  142. C           On fournit 2 maillage vide + les 2 initiaux
  143.             NBNN=0
  144.             NBELEM=0
  145.             NBSOUS=0
  146.             NBREF=0
  147.             SEGINI,IPT3
  148.             IPT3.ITYPEL=0
  149.             IOB3=IPT3
  150.             SEGDES,IPT3
  151.             IF(IRETOV.EQ.1)THEN
  152.               SEGINI,IPT3,IPT4
  153.               IPT4.ITYPEL=0
  154.               IOB4=IPT4
  155.               SEGDES,IPT4
  156.             ENDIF
  157.           endif
  158.           IF(IRETOV.EQ.1) CALL ECROBJ('MAILLAGE',IOB4)
  159.                           CALL ECROBJ('MAILLAGE',IOB3)
  160.           IF(IRETOV.EQ.1) CALL ECROBJ('MAILLAGE',IOB2)
  161.                           CALL ECROBJ('MAILLAGE',IOB1)
  162.           RETURN
  163.         endif
  164.  
  165.       ENDIF
  166.  
  167.  
  168. C-----------------------------------------------------------------------
  169. C     CREATION D'UNE LIGNE RESULTAT D'UNE INTERSECTION DE SURFACES
  170. C     DEFINIES ANALYTIQUEMENT (PLAN, SPHERE, ...)
  171. C
  172.       IF (IDIM.NE.3) CALL ERREUR(14)
  173.       IF (IERR.NE.0) RETURN
  174.  
  175. C  Y A T'IL UN FACTEUR DE DECOUPAGE
  176.       INBR=0
  177.       CALL LIRENT(INBR,0,IRETOU)
  178.       IMPOI=0
  179.       IMPOF=0
  180.       IMPLON=0
  181.       SEGACT MCOORD*mod
  182. C  LECTURE DU PREMIER POINT . SI C'EST UN ELEMENT ON EN EXTRAIT LE
  183. C  DERNIER POINT .
  184.       IOB1=0
  185.       CALL LIROBJ('MAILLAGE',IOB1,0,IRETOU)
  186.       IF (IRETOU.EQ.0) GOTO 1
  187.       IPT6=IOB1
  188.       SEGACT IPT6
  189.       CALL EXTRPO(IOB1,1,IPOIN1)
  190.       IF(IERR .NE. 0)RETURN
  191.       GOTO 2
  192.    1  CALL LIROBJ('POINT   ',IPOIN1,1,IRETOU)
  193.    2  IF (IERR.NE.0) RETURN
  194. C  PREMIERE SURFACE
  195.       CALL LIRMOT(MCLE,8,ICLE,1)
  196.       IF (IERR.NE.0) RETURN
  197. C  Y A T-IL DES DENSITES IMPOSEES
  198.   82  IF (IMPOI.EQ.1.OR.ICLE.NE.6) GOTO 80
  199.       IMPOI=1
  200.       CALL LIRREE(TPOIN1,1,IRETOU)
  201.       CALL LIRMOT(MCLE,8,ICLE,1)
  202.       IF (IERR.NE.0) RETURN
  203.   80  IF (IMPOF.EQ.1.OR.ICLE.NE.7) GOTO 81
  204.       IMPOF=1
  205.       CALL LIRREE(TPOIN2,1,IRETOU)
  206.       CALL LIRMOT(MCLE,8,ICLE,1)
  207.       IF (IERR.NE.0) RETURN
  208.       IF (IMPOI.EQ.0) GOTO 82
  209.   81  IF (IMPLON.EQ.1.OR.ICLE.NE.8) GOTO 800
  210.       IMPLON=1
  211.       CALL LIRREE(DLONG1,1,IRETOU)
  212.       CALL LIRMOT(MCLE,8,ICLE,1)
  213.       IF (IERR.NE.0) RETURN
  214.       goto 82
  215.  800  continue
  216. C  ACQUISITION DES DONNEES DE LA PREMIERE SURFACE
  217.       GOTO (20,21,22,23,24),ICLE
  218.   20  CALL LIROBJ('POINT   ',IPL,1,IRETOU)
  219.       IF (IERR.NE.0) RETURN
  220.       IREF=4*(IPL-1)
  221.       XIPL=XCOOR(IREF+1)
  222.       YIPL=XCOOR(IREF+2)
  223.       ZIPL=XCOOR(IREF+3)
  224.       GOTO 25
  225.   21  CALL LIROBJ('POINT   ',IPL,1,IRETOU)
  226.       IF (IERR.NE.0) RETURN
  227.       IREF=4*(IPL-1)
  228.       XIPL=XCOOR(IREF+1)
  229.       YIPL=XCOOR(IREF+2)
  230.       ZIPL=XCOOR(IREF+3)
  231.       GOTO 25
  232.   22  CALL LIROBJ('POINT   ',IPL,1,IRETOU)
  233.       IF (IERR.NE.0) RETURN
  234.       IREF=4*(IPL-1)
  235.       XIPL=XCOOR(IREF+1)
  236.       YIPL=XCOOR(IREF+2)
  237.       ZIPL=XCOOR(IREF+3)
  238.       CALL LIROBJ('POINT   ',IPL2,1,IRETOU)
  239.       IF (IERR.NE.0) RETURN
  240.       IREF=4*(IPL2-1)
  241.       XIV=XCOOR(IREF+1)-XIPL
  242.       YIV=XCOOR(IREF+2)-YIPL
  243.       ZIV=XCOOR(IREF+3)-ZIPL
  244.       S=SQRT(XIV*XIV+YIV*YIV+ZIV*ZIV)
  245.       IF (S.EQ.0.) THEN
  246.         CALL ERREUR(21)
  247.         RETURN
  248.       ENDIF
  249.       XIV=XIV/S
  250.       YIV=YIV/S
  251.       ZIV=ZIV/S
  252.       GOTO 25
  253.   23  CALL LIROBJ('POINT   ',IPL,1,IRETOU)
  254.       IF (IERR.NE.0) RETURN
  255.       IREF=4*(IPL-1)
  256.       XIPL=XCOOR(IREF+1)
  257.       YIPL=XCOOR(IREF+2)
  258.       ZIPL=XCOOR(IREF+3)
  259.       CALL LIROBJ('POINT   ',IPL2,1,IRETOU)
  260.       IF (IERR.NE.0) RETURN
  261.       IREF=4*(IPL2-1)
  262.       XIV=XCOOR(IREF+1)-XIPL
  263.       YIV=XCOOR(IREF+2)-YIPL
  264.       ZIV=XCOOR(IREF+3)-ZIPL
  265.       S=SQRT(XIV**2+YIV**2+ZIV**2)
  266.       IF (S.EQ.0.) THEN
  267.         CALL ERREUR(21)
  268.         RETURN
  269.       ENDIF
  270.       XIV=XIV/S
  271.       YIV=YIV/S
  272.       ZIV=ZIV/S
  273.       GOTO 25
  274.   24  CALL LIROBJ('POINT   ',IPL,1,IRETOU)
  275.       IF (IERR.NE.0) RETURN
  276.       IREF=4*(IPL-1)
  277.       XIPL=XCOOR(IREF+1)
  278.       YIPL=XCOOR(IREF+2)
  279.       ZIPL=XCOOR(IREF+3)
  280.       CALL LIROBJ('POINT   ',IAXE,1,IRETOU)
  281.       IF (IERR.NE.0) RETURN
  282.       IREF=4*(IAXE-1)
  283.       XIAXE=XCOOR(IREF+1)-XIPL
  284.       YIAXE=XCOOR(IREF+2)-YIPL
  285.       ZIAXE=XCOOR(IREF+3)-ZIPL
  286.       SIAXE=SQRT(XIAXE**2+YIAXE**2+ZIAXE**2)
  287.       IF (SIAXE.EQ.0.) THEN
  288.         CALL ERREUR(21)
  289.         RETURN
  290.       ENDIF
  291.       XIAXE=XIAXE/SIAXE
  292.       YIAXE=YIAXE/SIAXE
  293.       ZIAXE=ZIAXE/SIAXE
  294.       CALL LIROBJ('POINT   ',IP1,1,IRETOU)
  295.       IF (IERR.NE.0) RETURN
  296.       IREF=4*(IP1-1)
  297.       XIP1=XCOOR(IREF+1)
  298.       YIP1=XCOOR(IREF+2)
  299.       ZIP1=XCOOR(IREF+3)
  300.       GRAYI=SQRT((XIPL-XIP1)**2+(YIPL-YIP1)**2+(ZIPL-ZIP1)**2)
  301.       IF (GRAYI.EQ.0.) THEN
  302.         CALL ERREUR(21)
  303.         RETURN
  304.       ENDIF
  305.   25  CONTINUE
  306. C  DEUXIEME SURFACE
  307.       CALL LIRMOT(MCLE,8,JCLE,1)
  308.       IF (IERR.NE.0) RETURN
  309. C  Y A T-IL DES DENSITES IMPOSEES
  310.   85  IF (IMPOI.EQ.1.OR.JCLE.NE.6) GOTO 83
  311.       IMPOI=1
  312.       CALL LIRREE(TPOIN1,1,IRETOU)
  313.       CALL LIRMOT(MCLE,8,JCLE,1)
  314.       IF (IERR.NE.0) RETURN
  315.   83  IF (IMPOF.EQ.1.OR.JCLE.NE.7) GOTO 84
  316.       IMPOF=1
  317.       CALL LIRREE(TPOIN2,1,IRETOU)
  318.       CALL LIRMOT(MCLE,8,JCLE,1)
  319.       IF (IERR.NE.0) RETURN
  320.       IF (IMPOI.EQ.0) GOTO 85
  321.   84  IF (IMPLON.EQ.1.OR.ICLE.NE.8) GOTO 801
  322.       IMPLON=1
  323.       CALL LIRREE(DLONG1,1,IRETOU)
  324.       CALL LIRMOT(MCLE,8,ICLE,1)
  325.       IF (IERR.NE.0) RETURN
  326.       goto 85
  327.  801  continue
  328. C  AQUISITION DES DONNEES DE LA DEUXIEME SURFACE
  329.       GOTO (40,41,42,43,44),JCLE
  330.   40  CALL LIROBJ('POINT   ',JPL,1,IRETOU)
  331.       IF (IERR.NE.0) RETURN
  332.       IREF=4*(JPL-1)
  333.       XJPL=XCOOR(IREF+1)
  334.       YJPL=XCOOR(IREF+2)
  335.       ZJPL=XCOOR(IREF+3)
  336.       GOTO 45
  337.   41  CALL LIROBJ('POINT   ',JPL,1,IRETOU)
  338.       IF (IERR.NE.0) RETURN
  339.       IREF=4*(JPL-1)
  340.       XJPL=XCOOR(IREF+1)
  341.       YJPL=XCOOR(IREF+2)
  342.       ZJPL=XCOOR(IREF+3)
  343.       GOTO 45
  344.   42  CALL LIROBJ('POINT   ',JPL,1,IRETOU)
  345.       IF (IERR.NE.0) RETURN
  346.       IREF=4*(JPL-1)
  347.       XJPL=XCOOR(IREF+1)
  348.       YJPL=XCOOR(IREF+2)
  349.       ZJPL=XCOOR(IREF+3)
  350.       CALL LIROBJ('POINT   ',JPL2,1,IRETOU)
  351.       IF (IERR.NE.0) RETURN
  352.       IREF=4*(JPL2-1)
  353.       XJV=XCOOR(IREF+1)-XJPL
  354.       YJV=XCOOR(IREF+2)-YJPL
  355.       ZJV=XCOOR(IREF+3)-ZJPL
  356.       S=SQRT(XJV**2+YJV**2+ZJV**2)
  357.       IF(S.EQ.0.) THEN
  358.         CALL ERREUR(21)
  359.         RETURN
  360.       ENDIF
  361.       XJV=XJV/S
  362.       YJV=YJV/S
  363.       ZJV=ZJV/S
  364.       GOTO 45
  365.   43  CALL LIROBJ('POINT   ',JPL,1,IRETOU)
  366.       IF (IERR.NE.0) RETURN
  367.       IREF=4*(JPL-1)
  368.       XJPL=XCOOR(IREF+1)
  369.       YJPL=XCOOR(IREF+2)
  370.       ZJPL=XCOOR(IREF+3)
  371.       CALL LIROBJ('POINT   ',JPL2,1,IRETOU)
  372.       IF (IERR.NE.0) RETURN
  373.       IREF=4*(JPL2-1)
  374.       XJV=XCOOR(IREF+1)-XJPL
  375.       YJV=XCOOR(IREF+2)-YJPL
  376.       ZJV=XCOOR(IREF+3)-ZJPL
  377.       S=SQRT(XJV**2+YJV**2+ZJV**2)
  378.       IF (S.EQ.0.) THEN
  379.         CALL ERREUR(21)
  380.         RETURN
  381.       ENDIF
  382.       XJV=XJV/S
  383.       YJV=YJV/S
  384.       ZJV=ZJV/S
  385.       GOTO 45
  386.   44  CALL LIROBJ('POINT   ',JPL,1,IRETOU)
  387.       IF (IERR.NE.0) RETURN
  388.       IREF=4*(JPL-1)
  389.       XJPL=XCOOR(IREF+1)
  390.       YJPL=XCOOR(IREF+2)
  391.       ZJPL=XCOOR(IREF+3)
  392.       CALL LIROBJ('POINT   ',JAXE,1,IRETOU)
  393.       IF (IERR.NE.0) RETURN
  394.       IREF=4*(JAXE-1)
  395.       XJAXE=XCOOR(IREF+1)-XJPL
  396.       YJAXE=XCOOR(IREF+2)-YJPL
  397.       ZJAXE=XCOOR(IREF+3)-ZJPL
  398.       SJAXE=SQRT(XJAXE**2+YJAXE**2+ZJAXE**2)
  399.       IF (SJAXE.EQ.0.) THEN
  400.         CALL ERREUR(21)
  401.         RETURN
  402.       ENDIF
  403.       XJAXE=XJAXE/SJAXE
  404.       YJAXE=YJAXE/SJAXE
  405.       ZJAXE=ZJAXE/SJAXE
  406.       CALL LIROBJ('POINT   ',JP1,1,IRETOU)
  407.       IF (IERR.NE.0) RETURN
  408.       IREF=4*(JP1-1)
  409.       XJP1=XCOOR(IREF+1)
  410.       YJP1=XCOOR(IREF+2)
  411.       ZJP1=XCOOR(IREF+3)
  412.       GRAYJ=SQRT((XJPL-XJP1)**2+(YJPL-YJP1)**2+(ZJPL-ZJP1)**2)
  413.       IF (GRAYJ.EQ.0.) THEN
  414.         CALL ERREUR(21)
  415.         RETURN
  416.       ENDIF
  417.   45  CONTINUE
  418. C  Y A T-IL DES DENSITES IMPOSEES
  419.       IRETOU=0
  420.   88  IF (IMPOI*IMPOF*implon.EQ.0) CALL LIRMOT(MCLE(6),3,MDPOS,0)
  421.       IF (MDPOS.EQ.0) GOTO 89
  422.       IF (IMPOI.EQ.1.OR.MDPOS.NE.1) GOTO 86
  423.       IMPOI=1
  424.       CALL LIRREE(TPOIN1,1,IRETOU)
  425.       goto 88
  426.   86  IF (IMPOF.EQ.1.OR.MDPOS.NE.2) GOTO 87
  427.       IMPOF=1
  428.       CALL LIRREE(TPOIN2,1,IRETOU)
  429.       GOTO 88
  430.   87  if (implon.EQ.1.OR.MDPOS.NE.3) GOTO 802
  431.       implon=1
  432.       CALL LIRREE(dlong1,1,IRETOU)
  433.       goto 88
  434.  802  CALL REFUS
  435.   89  CONTINUE
  436.       if (ierr.ne.0) return
  437. C  LECTURE DU DEUXIEME POINT . SI C'EST UN ELEMENT ON EN EXTRAIT LE
  438. C  PREMIER POINT .
  439.       IOB2=0
  440.       iretou=0
  441.       if (implon.eq.0) CALL LIROBJ('MAILLAGE',IOB2,0,IRETOU)
  442.       IF (IRETOU.EQ.0) GOTO 50
  443.       IPT6=IOB2
  444.       SEGACT IPT6
  445.       CALL EXTRPO(IOB2,2,IPOIN2)
  446.       IF(IERR .NE. 0)RETURN
  447.       GOTO 51
  448.   50  CALL LIROBJ('POINT   ',IPOIN2,1,IRETOU)
  449.   51  IF (IERR.NE.0) RETURN
  450.       NBP=NBT-2
  451.       IREF1=(IPOIN1-1)*4
  452.       IREF2=(IPOIN2-1)*4
  453.       XPOIN1=XCOOR(IREF1+1)
  454.       YPOIN1=XCOOR(IREF1+2)
  455.       ZPOIN1=XCOOR(IREF1+3)
  456.       IF (IMPOI.EQ.0) TPOIN1=XCOOR(IREF1+4)
  457.       XPOIN2=XCOOR(IREF2+1)
  458.       YPOIN2=XCOOR(IREF2+2)
  459.       ZPOIN2=XCOOR(IREF2+3)
  460.       IF (IMPOF.EQ.0) TPOIN2=XCOOR(IREF2+4)
  461.       XVEC=XPOIN2-XPOIN1
  462.       YVEC=YPOIN2-YPOIN1
  463.       ZVEC=ZPOIN2-ZPOIN1
  464.       SVEC=SQRT(XVEC**2+YVEC**2+ZVEC**2)
  465.       IF (SVEC.EQ.0.) THEN
  466.         CALL ERREUR(21)
  467.         RETURN
  468.       ENDIF
  469.       EPS2=1E-10*SVEC**2
  470. C  IL FAUT VERIFIER LA COMPATIBILITE DES DONNEES
  471.       GOTO (60,61,62,63,64),ICLE
  472.   60  CONTINUE
  473.       XPV1=XPOIN1-XIPL
  474.       YPV1=YPOIN1-YIPL
  475.       ZPV1=ZPOIN1-ZIPL
  476.       XPV2=XPOIN2-XIPL
  477.       YPV2=YPOIN2-YIPL
  478.       ZPV2=ZPOIN2-ZIPL
  479.       XIPN=YPV1*ZPV2-ZPV1*YPV2
  480.       YIPN=ZPV1*XPV2-XPV1*ZPV2
  481.       ZIPN=XPV1*YPV2-YPV1*XPV2
  482.       SIPN=SQRT(XIPN**2+YIPN**2+ZIPN**2)
  483.       IF (SIPN.EQ.0.) CALL ERREUR(21)
  484.       IF (IERR.NE.0) RETURN
  485.       XIPN=XIPN/SIPN
  486.       YIPN=YIPN/SIPN
  487.       ZIPN=ZIPN/SIPN
  488.       GOTO 65
  489.   61  CONTINUE
  490.       RAYI=SQRT((XPOIN1-XIPL)**2+(YPOIN1-YIPL)**2+(ZPOIN1-ZIPL)**2)
  491.       RAYB=SQRT((XPOIN2-XIPL)**2+(YPOIN2-YIPL)**2+(ZPOIN2-ZIPL)**2)
  492.       if (implon.eq.0) then
  493.        IF (ABS((RAYI-RAYB)/RAYI).GE.1E-2) THEN
  494.          CALL ERREUR(21)
  495.          RETURN
  496.        ENDIF
  497.        RAYI=SQRT(RAYI*RAYB)
  498.       endif
  499.       GOTO 65
  500.   62  CONTINUE
  501.       U=XPOIN1-XIPL
  502.       V=YPOIN1-YIPL
  503.       W=ZPOIN1-ZIPL
  504.       SCA=U*XIV+V*YIV+W*ZIV
  505.       U=U-SCA*XIV
  506.       V=V-SCA*YIV
  507.       W=W-SCA*ZIV
  508.       RAYI=SQRT(U**2+V**2+W**2)
  509.       U=XPOIN2-XIPL
  510.       V=YPOIN2-YIPL
  511.       W=ZPOIN2-ZIPL
  512.       SCA=U*XIV+V*YIV+W*ZIV
  513.       U=U-SCA*XIV
  514.       V=V-SCA*YIV
  515.       W=W-SCA*ZIV
  516.       RAYB=SQRT(U**2+V**2+W**2)
  517.       if (implon.eq.0) then
  518.        IF (ABS(RAYI-RAYB)/RAYI.GT.1E-2) CALL ERREUR(21)
  519.        RAYI=SQRT(RAYI*RAYB)
  520.       endif
  521.       GOTO 65
  522.   63  CONTINUE
  523.       U=XPOIN1-XIPL
  524.       V=YPOIN1-YIPL
  525.       W=ZPOIN1-ZIPL
  526.       SCA=ABS(U*XIV+V*YIV+W*ZIV)
  527.       U1=V*ZIV-W*YIV
  528.       V1=W*XIV-U*ZIV
  529.       W1=U*YIV-V*XIV
  530.       VEC=SQRT(U1**2+V1**2+W1**2)
  531.       IF (SCA.EQ.0.) CALL ERREUR(21)
  532.       IF (IERR.NE.0) RETURN
  533.       TANGI=VEC/SCA
  534.       U=XPOIN2-XIPL
  535.       V=YPOIN2-YIPL
  536.       W=ZPOIN2-ZIPL
  537.       SCA=ABS(U*XIV+V*YIV+W*ZIV)
  538.       U1=V*ZIV-W*YIV
  539.       V1=W*XIV-U*ZIV
  540.       W1=U*YIV-V*XIV
  541.       VEC=SQRT(U1**2+V1**2+W1**2)
  542.       IF (SCA.EQ.0.) CALL ERREUR(21)
  543.       IF (IERR.NE.0) RETURN
  544.       TANGB=VEC/SCA
  545.       if (implon.eq.0) then
  546.        IF (ABS(TANGI-TANGB)/TANGI.GE.1E-2) CALL ERREUR(21)
  547.        TANGI=SQRT(TANGI*TANGB)
  548.       endif
  549.       GOTO 65
  550.   64  CONTINUE
  551.       U=XPOIN1-XIPL
  552.       V=YPOIN1-YIPL
  553.       W=ZPOIN1-ZIPL
  554.       SCA=U*XIAXE+V*YIAXE+W*ZIAXE
  555.       U=U-SCA*XIAXE
  556.       V=V-SCA*YIAXE
  557.       W=W-SCA*ZIAXE
  558.       S=SQRT(U**2+V**2+W**2)
  559.       IF (S.EQ.0.) CALL ERREUR(21)
  560.       IF (IERR.NE.0) RETURN
  561.       XC=XIPL+U*GRAYI/S
  562.       YC=YIPL+V*GRAYI/S
  563.       ZC=ZIPL+W*GRAYI/S
  564.       PRAYI=SQRT((XPOIN1-XC)**2+(YPOIN1-YC)**2+(ZPOIN1-ZC)**2)
  565.       U=XPOIN2-XIPL
  566.       V=YPOIN2-YIPL
  567.       W=ZPOIN2-ZIPL
  568.       SCA=U*XIAXE+V*YIAXE+W*ZIAXE
  569.       U=U-SCA*XIAXE
  570.       V=V-SCA*YIAXE
  571.       W=W-SCA*ZIAXE
  572.       S=SQRT(U**2+V**2+W**2)
  573.       IF (S.EQ.0.) CALL ERREUR(21)
  574.       IF (IERR.NE.0) RETURN
  575.       XC=XIPL+U*GRAYI/S
  576.       YC=YIPL+V*GRAYI/S
  577.       ZC=ZIPL+W*GRAYI/S
  578.       PRAYB=SQRT((XPOIN2-XC)**2+(YPOIN2-YC)**2+(ZPOIN2-ZC)**2)
  579.       if (implon.eq.0) then
  580.        IF (ABS(PRAYI-PRAYB)/PRAYI.GE.1E-2) CALL ERREUR(21)
  581.        PRAYI=SQRT(PRAYI*PRAYB)
  582.       endif
  583.   65  IF (IERR.NE.0) RETURN
  584.       GOTO (70,71,72,73,74),JCLE
  585.   70  CONTINUE
  586.       XPV1=XPOIN1-XJPL
  587.       YPV1=YPOIN1-YJPL
  588.       ZPV1=ZPOIN1-ZJPL
  589.       XPV2=XPOIN2-XJPL
  590.       YPV2=YPOIN2-YJPL
  591.       ZPV2=ZPOIN2-ZJPL
  592.       XJPN=YPV1*ZPV2-ZPV1*YPV2
  593.       YJPN=ZPV1*XPV2-XPV1*ZPV2
  594.       ZJPN=XPV1*YPV2-YPV1*XPV2
  595.       SJPN=SQRT(XJPN**2+YJPN**2+ZJPN**2)
  596.       IF (SJPN.EQ.0.) CALL ERREUR(21)
  597.       IF (IERR.NE.0) RETURN
  598.       XJPN=XJPN/SJPN
  599.       YJPN=YJPN/SJPN
  600.       ZJPN=ZJPN/SJPN
  601.       GOTO 75
  602.   71  CONTINUE
  603.       RAYJ=SQRT((XPOIN1-XJPL)**2+(YPOIN1-YJPL)**2+(ZPOIN1-ZJPL)**2)
  604.       RAYB=SQRT((XPOIN2-XJPL)**2+(YPOIN2-YJPL)**2+(ZPOIN2-ZJPL)**2)
  605.       if (implon.eq.0) then
  606.        IF (ABS((RAYJ-RAYB)/RAYJ).GE.1E-2) CALL ERREUR(21)
  607.        RAYJ=SQRT(RAYJ*RAYB)
  608.       endif
  609.       GOTO 75
  610.   72  CONTINUE
  611.       U=XPOIN1-XJPL
  612.       V=YPOIN1-YJPL
  613.       W=ZPOIN1-ZJPL
  614.       SCA=U*XJV+V*YJV+W*ZJV
  615.       U=U-SCA*XJV
  616.       V=V-SCA*YJV
  617.       W=W-SCA*ZJV
  618.       RAYJ=SQRT(U**2+V**2+W**2)
  619.       U=XPOIN2-XJPL
  620.       V=YPOIN2-YJPL
  621.       W=ZPOIN2-ZJPL
  622.       SCA=U*XJV+V*YJV+W*ZJV
  623.       U=U-SCA*XJV
  624.       V=V-SCA*YJV
  625.       W=W-SCA*ZJV
  626.       RAYB=SQRT(U**2+V**2+W**2)
  627.       if (implon.eq.0) then
  628.        IF (ABS(RAYJ-RAYB)/RAYJ.GT.1E-2) CALL ERREUR(21)
  629.        RAYJ=SQRT(RAYJ*RAYB)
  630.       endif
  631.       GOTO 75
  632.   73  CONTINUE
  633.       U=XPOIN1-XJPL
  634.       V=YPOIN1-YJPL
  635.       W=ZPOIN1-ZJPL
  636.       SCA=ABS(U*XJV+V*YJV+W*ZJV)
  637.       U1=V*ZJV-W*YJV
  638.       V1=W*XJV-U*ZJV
  639.       W1=U*YJV-V*XJV
  640.       VEC=SQRT(U1**2+V1**2+W1**2)
  641.       IF (SCA.EQ.0.) CALL ERREUR(21)
  642.       IF (IERR.NE.0) RETURN
  643.       TANGJ=VEC/SCA
  644.       U=XPOIN2-XJPL
  645.       V=YPOIN2-YJPL
  646.       W=ZPOIN2-ZJPL
  647.       SCA=ABS(U*XJV+V*YJV+W*ZJV)
  648.       U1=V*ZJV-W*YJV
  649.       V1=W*XJV-U*ZJV
  650.       W1=U*YJV-V*XJV
  651.       VEC=SQRT(U1**2+V1**2+W1**2)
  652.       IF (SCA.EQ.0.) CALL ERREUR(21)
  653.       IF (IERR.NE.0) RETURN
  654.       TANGB=VEC/SCA
  655.       if (implon.eq.0) then
  656.        IF (ABS(TANGJ-TANGB)/TANGJ.GE.1E-2) CALL ERREUR(21)
  657.        TANGJ=SQRT(TANGJ*TANGB)
  658.       endif
  659.       GOTO 75
  660.   74  CONTINUE
  661.       U=XPOIN1-XJPL
  662.       V=YPOIN1-YJPL
  663.       W=ZPOIN1-ZJPL
  664.       SCA=U*XJAXE+V*YJAXE+W*ZJAXE
  665.       U=U-SCA*XJAXE
  666.       V=V-SCA*YJAXE
  667.       W=W-SCA*ZJAXE
  668.       S=SQRT(U**2+V**2+W**2)
  669.       IF (S.EQ.0.) CALL ERREUR(21)
  670.       IF (IERR.NE.0) RETURN
  671.       XC=XJPL+U*GRAYJ/S
  672.       YC=YJPL+V*GRAYJ/S
  673.       ZC=ZJPL+W*GRAYJ/S
  674.       PRAYJ=SQRT((XPOIN1-XC)**2+(YPOIN1-YC)**2+(ZPOIN1-ZC)**2)
  675.       U=XPOIN2-XJPL
  676.       V=YPOIN2-YJPL
  677.       W=ZPOIN2-ZJPL
  678.       SCA=U*XJAXE+V*YJAXE+W*ZJAXE
  679.       U=U-SCA*XJAXE
  680.       V=V-SCA*YJAXE
  681.       W=W-SCA*ZJAXE
  682.       S=SQRT(U**2+V**2+W**2)
  683.       IF (S.EQ.0.) CALL ERREUR(21)
  684.       IF (IERR.NE.0) RETURN
  685.       XC=XJPL+U*GRAYJ/S
  686.       YC=YJPL+V*GRAYJ/S
  687.       ZC=ZJPL+W*GRAYJ/S
  688.       PRAYB=SQRT((XPOIN2-XC)**2+(YPOIN2-YC)**2+(ZPOIN2-ZC)**2)
  689.       if (implon.eq.0) then
  690.        IF (ABS(PRAYJ-PRAYB)/PRAYJ.GE.1E-2) CALL ERREUR(21)
  691.        PRAYJ=SQRT(PRAYJ*PRAYB)
  692.       endif
  693.   75  IF (IERR.NE.0) RETURN
  694. C  ON PLACE NBP POINTS SUR L'INTERSECTION POUR CALCULER DES ABCISSES
  695. C  CURVILIGNES
  696.       XVECR=XVEC/(NBP+1)
  697.       YVECR=YVEC/(NBP+1)
  698.       ZVECR=ZVEC/(NBP+1)
  699.       XVEC=XVEC/SVEC
  700.       YVEC=YVEC/SVEC
  701.       ZVEC=ZVEC/SVEC
  702.       XPAUX(1,1)=XPOIN1
  703.       XPAUX(2,1)=YPOIN1
  704.       XPAUX(3,1)=ZPOIN1
  705.       XPAUX(4,1)=0.
  706.       IBP=0
  707.  102  IBP=IBP+1
  708.       IF (IBP.GT.NBP) GOTO 100
  709.       XP=XPAUX(1,IBP)+XVECR
  710.       YP=XPAUX(2,IBP)+YVECR
  711.       ZP=XPAUX(3,IBP)+ZVECR
  712.       XPA=XP
  713.       YPA=YP
  714.       ZPA=ZP
  715.       IRET = 101
  716.       GOTO 400
  717.  101  CONTINUE
  718.       IF (IIMPI.EQ.1) WRITE (IOIMP,8880) IBP,XPA,YPA,ZPA,XP1,YP1,ZP1
  719.  8880 FORMAT(I4,' POINT A PROJETER ',3G13.5,' POINT PROJETE ',3G13.5)
  720. C  ON RANGE LE POINT D'INTERSECTION
  721.       XPAUX(1,IBP+1)=XP1
  722.       XPAUX(2,IBP+1)=YP1
  723.       XPAUX(3,IBP+1)=ZP1
  724.       XPAUX(4,IBP+1)=XPAUX(4,IBP)+SQRT((XP1-XPAUX(1,IBP))**2+
  725.      # (YP1-XPAUX(2,IBP))**2+(ZP1-XPAUX(3,IBP))**2)
  726.       GOTO 102
  727.  100  CONTINUE
  728.       XPAUX(1,NBT)=XPOIN2
  729.       XPAUX(2,NBT)=YPOIN2
  730.       XPAUX(3,NBT)=ZPOIN2
  731.       XPAUX(4,NBT)=XPAUX(4,NBT-1)+SQRT((XP1-XPOIN2)**2+(YP1-YPOIN2)**2+
  732.      # (ZP1-ZPOIN2)**2)
  733.       DLONG=XPAUX(4,NBT)
  734.       if (implon.eq.1) then
  735.        if ( iimpi.eq.1) write(ioimp,8888) dlong1 ,dlong
  736.  8888 format ( ' longueur demandee et longueur donnee',2g13.5)
  737.        if (dlong1.gt.dlong) call erreur(21)
  738.        if (ierr.ne.0) return
  739.        dlong=dlong1
  740.       endif
  741.       DEN1=TPOIN1/DLONG
  742.       DEN2=TPOIN2/DLONG
  743.       CALL DECOUP(INBR,DEN1,DEN2,APROG,NBELEM,DENI,DECA,DLONG)
  744.       NX=NBELEM-1
  745.       DEN1=DEN1*DLONG
  746.       IF (IIMPI.EQ.1) WRITE (IOIMP,9900) DLONG,NBELEM,APROG
  747.  9900 FORMAT (' LONGUEUR ',G12.5,' ELEMENT ',I5,' RAISON ',G12.5)
  748. C  CREATION DU SEGMENT
  749.       NBNN=NBNNE(KDEGRE(ILCOUR))
  750.       IF (ILCOUR.EQ.0) CALL ERREUR(16)
  751.       IF (KDEGRE(ILCOUR).EQ.0) CALL ERREUR(16)
  752.       IF (NBNN.EQ.0) CALL ERREUR(16)
  753.       IF (IERR.NE.0) RETURN
  754.       NBSOUS=0
  755.       NBREF=0
  756.       SEGINI MELEME
  757.       ITYPEL=KDEGRE(ILCOUR)
  758.       NPOIN=(ITYPEL-1)*NBELEM-1
  759.       if (implon.eq.1) npoin=npoin+1
  760.       SEGINI TABPAR
  761. C  DEFINITION DE LA TOPOLOGIE
  762.       IPOO=nbpts
  763.       IADR=IPOO
  764.       DIN=DEN1
  765.       DL=0.
  766.       NUM(1,1)=IPOIN1
  767.       IF (NX.EQ.0) GOTO 510
  768.       DO 512 I=1,NX
  769.       DIN=DIN*APROG
  770.       IF (ITYPEL.EQ.2) GOTO 513
  771.       IPOO=IPOO+1
  772.       NUM(2,I)=IPOO
  773.       TABPAR(IPOO-IADR)=DIN/2.+DL
  774.  513  IPOO=IPOO+1
  775.       NUM(ITYPEL,I)=IPOO
  776.       NUM(1,I+1)=IPOO
  777.       DL=DL+DIN
  778.       TABPAR(IPOO-IADR)=DL
  779.  512  CONTINUE
  780.  510  CONTINUE
  781.       IF (ITYPEL.NE.3) GOTO 514
  782.       DIN=DIN*APROG
  783.       IPOO=IPOO+1
  784.       TABPAR(IPOO-IADR)=DL+DIN/2
  785.       NUM(2,NBELEM)=IPOO
  786.  514  NUM(ITYPEL,NBELEM)=IPOIN2
  787.       if (implon.eq.1) then
  788.        NUM(ITYPEL,NBELEM)=ipoo+1
  789.        TABPAR(IPOO+1-IADR)=DLong
  790.       endif
  791. C   CREATION DES POINTS  LA LONGUEUR A RESPECTER EST DANS TABPAR
  792. C   LA LONGUEUR DES POINTS CONNUS EST DANS XPAUX
  793.       IP=1
  794.       T=0.
  795.       DLR=0.
  796.       NBPTA=nbpts
  797.       NBPTS=NBPTA+NPOIN
  798.       SEGADJ MCOORD
  799.       I=0
  800.  523  I=I+1
  801.       IF (I.GT.NPOIN) GOTO 520
  802.       DLV=TABPAR(I)
  803.  522  IF (DLV.LE.DLR) GOTO 521
  804.       IP=IP+1
  805.       DLR=XPAUX(4,IP)
  806.       GOTO 522
  807.  521  CONTINUE
  808.       DLA=XPAUX(4,IP-1)
  809.       COF=(DLV-DLA)/(DLR-DLA)
  810.       XP=XPAUX(1,IP-1)+COF*(XPAUX(1,IP)-XPAUX(1,IP-1))
  811.       YP=XPAUX(2,IP-1)+COF*(XPAUX(2,IP)-XPAUX(2,IP-1))
  812.       ZP=XPAUX(3,IP-1)+COF*(XPAUX(3,IP)-XPAUX(3,IP-1))
  813.       IF (IIMPI.NE.0) WRITE (IOIMP,8878) I,IP,XP,YP,ZP
  814.  8878 FORMAT(' AVANT ',I4,' INTERPOLE ',I4,3G13.5)
  815.       IRET = 525
  816.       GOTO 400
  817.  525  CONTINUE
  818.       IF (IIMPI.NE.0) WRITE (IOIMP,8876) I,IP,XP1,YP1,ZP1
  819.  8876 FORMAT(' APRES ',I4,' INTERPOLE ',I4,3G13.5)
  820.       XCOOR(NBPTA*4+1)=XP1
  821.       XCOOR(NBPTA*4+2)=YP1
  822.       XCOOR(NBPTA*4+3)=ZP1
  823.       XCOOR(NBPTA*4+4)=TABPAR(I)-T
  824.       NBPTA=NBPTA+1
  825.       T=TABPAR(I)
  826.       IADR=IADR+1
  827.       GOTO 523
  828.  520  CONTINUE
  829.       SEGSUP TABPAR
  830.       IPT2=MELEME
  831.       DO 527 I=1,NUM(/2)
  832.         ICOLOR(I)=IDCOUL
  833.  527  CONTINUE
  834.       IF (IOB1.NE.0) THEN
  835.         ltelq=.false.
  836.         CALL FUSE(IOB1,MELEME,IPT2,ltelq)
  837.         IPT6=IOB1
  838.         SEGDES IPT6
  839.       ENDIF
  840.       MELEME=IPT2
  841.       IF (IOB2.NE.0) THEN
  842.         ltelq=.false.
  843.         CALL FUSE(MELEME,IOB2,IPT2,ltelq)
  844.         IPT6=IOB2
  845.         SEGDES IPT6
  846.       ENDIF
  847.       MELEME=IPT2
  848.       SEGDES MELEME
  849.       CALL ECROBJ('MAILLAGE',MELEME)
  850.       RETURN
  851. C  SOUS-PROGRAMME DE CALCUL ITERATIF D'INTERSECTION
  852.  400  CONTINUE
  853.       ICYCL=0
  854.  402  CONTINUE
  855. C  IL FAUT TROUVER L'INTERSECTION DU PLAN PERPENDICULAIRE A XVEC YVEC
  856. C  ZVEC EN XPOIN YPOIN ZPOIN ET DES DEUX SURFACES
  857. C  ON UTILISE UNE METHODE ITERATIVE
  858. C  ON VA CHERCHER LE POINT DE PROJECTION DE XP SUR CHAQUE SURFACE DANS
  859. C  LE PLAN ET UN VECTEUR NORMAL
  860. C  XP YP ZP  POINT A PROJETER
  861. C  XP1 YP1 ZP1 INTERSECTION EN RETOUR
  862.       GOTO (200,210,220,230,240),ICLE
  863.  250  CONTINUE
  864.       GOTO (300,310,320,330,340),JCLE
  865.  350  CONTINUE
  866. C   LES PTS ET LES VECTS SONT XPI... XPJ... XNI... XNJ...
  867. C  ON CALCULE L'INTERSECTION DES DEUX VECTEURS (AFFINES)
  868.       XDIF=XPJ-XPI
  869.       YDIF=YPJ-YPI
  870.       ZDIF=ZPJ-ZPI
  871. C  COMPOSANTES SUR CHAQUE VECTEUR
  872.       DET=(YNI*ZNJ-ZNI*YNJ)*XVEC
  873.      #   +(ZNI*XNJ-XNI*ZNJ)*YVEC
  874.      #   +(XNI*YNJ-YNI*XNJ)*ZVEC
  875.       IF (DET**2.LT.EPS2) CALL ERREUR(40)
  876.       A=((YDIF*ZNJ-ZDIF*YNJ)*XVEC
  877.      # +(ZDIF*XNJ-XDIF*ZNJ)*YVEC
  878.      # +(XDIF*YNJ-YDIF*XNJ)*ZVEC)/DET
  879.       XP1=XPI+A*XNI
  880.       YP1=YPI+A*YNI
  881.       ZP1=ZPI+A*ZNI
  882. C  TEST CONVERGENCE
  883.       ECART=(XP-XP1)**2+(YP-YP1)**2+(ZP-ZP1)**2
  884.       IF (ECART.LE.EPS2) THEN
  885.         IF (IRET.EQ.101) GOTO 101
  886.         IF (IRET.EQ.525) GOTO 525
  887.       ENDIF
  888.       XP=XP1
  889.       YP=YP1
  890.       ZP=ZP1
  891.       ICYCL=ICYCL+1
  892.       IF (IIMPI.EQ.1) WRITE (IOIMP,4478) ICYCL,ECART,EPS2,XPI,YPI,ZPI,
  893.      # XNI,YNI,ZNI,XPJ,YPJ,ZPJ,XNJ,YNJ,ZNJ,XP,YP,ZP
  894.  4478 FORMAT(' ICYCL ',I4,' ECART ',G12.5,' EPS2 ',G12.5,/,' PI ',
  895.      # 3G13.5,' NI ',3G13.5,/,' PJ ',3G13.5,' NJ ',3G13.5,' P ',3G13.5)
  896.       IF (ICYCL.LE.25) GOTO 402
  897.       CALL ERREUR(40)
  898.       RETURN
  899. C  PROJECTION SUR UN PLAN XIPL XIPN EN RESTANT DANS XVEC
  900.  200  CONTINUE
  901. C  VECT PERPENDICULAIRE A XVEC ET XIPN
  902.       U=YVEC*ZIPN-ZVEC*YIPN
  903.       V=ZVEC*XIPN-XVEC*ZIPN
  904.       W=XVEC*YIPN-YVEC*XIPN
  905.       S=SQRT(U**2+V**2+W**2)
  906.       IF (S.EQ.0.) CALL ERREUR(40)
  907.       U1=U/S
  908.       V1=V/S
  909.       W1=W/S
  910.       U=V1*ZVEC-W1*YVEC
  911.       V=W1*XVEC-U1*ZVEC
  912.       W=U1*YVEC-V1*XVEC
  913. C  VECT SUIVANT LEQUEL ON DEPLACE LE POINT
  914.       PAR=-((XP-XIPL)*XIPN+(YP-YIPL)*YIPN+(ZP-ZIPL)*ZIPN)/
  915.      # (U*XIPN+V*YIPN+W*ZIPN)
  916.       XPI=XP+PAR*U
  917.       YPI=YP+PAR*V
  918.       ZPI=ZP+PAR*W
  919.       XNI=U1
  920.       YNI=V1
  921.       ZNI=W1
  922.       GOTO 250
  923.  300  CONTINUE
  924.       U=YVEC*ZJPN-ZVEC*YJPN
  925.       V=ZVEC*XJPN-XVEC*ZJPN
  926.       W=XVEC*YJPN-YVEC*XJPN
  927.       S=SQRT(U**2+V**2+W**2)
  928.       IF (S.EQ.0.) CALL ERREUR(40)
  929.       U1=U/S
  930.       V1=V/S
  931.       W1=W/S
  932.       U=V1*ZVEC-W1*YVEC
  933.       V=W1*XVEC-U1*ZVEC
  934.       W=U1*YVEC-V1*XVEC
  935.       PAR=-((XP-XJPL)*XJPN+(YP-YJPL)*YJPN+(ZP-ZJPL)*ZJPN)/
  936.      # (U*XJPN+V*YJPN+W*ZJPN)
  937.       XPJ=XP+PAR*U
  938.       YPJ=YP+PAR*V
  939.       ZPJ=ZP+PAR*W
  940.       XNJ=U1
  941.       YNJ=V1
  942.       ZNJ=W1
  943.       GOTO 350
  944. C  PROJECTION SUR SPHERE CENTRE XIPL RAYON RAYI EN RESTANT SUR LE PLAN X
  945.  210  CONTINUE
  946.       U1=XIPL-XP
  947.       V1=YIPL-YP
  948.       W1=ZIPL-ZP
  949.       SCA=U1*XVEC+V1*YVEC+W1*ZVEC
  950.       XC=XIPL-SCA*XVEC
  951.       YC=YIPL-SCA*YVEC
  952.       ZC=ZIPL-SCA*ZVEC
  953.       IF (SCA.GT.RAYI) CALL ERREUR(40)
  954.       IF (IERR.NE.0) RETURN
  955.       RAYC=SQRT(RAYI**2-SCA**2)
  956.       U=XC-XP
  957.       V=YC-YP
  958.       W=ZC-ZP
  959.       S=SQRT(U**2+V**2+W**2)
  960.       IF (S.EQ.0.) CALL ERREUR(40)
  961.       IF (IERR.NE.0) RETURN
  962.       U=U/S
  963.       V=V/S
  964.       W=W/S
  965.       XPI=XC-RAYC*U
  966.       YPI=YC-RAYC*V
  967.       ZPI=ZC-RAYC*W
  968.       XNI=YVEC*W-ZVEC*V
  969.       YNI=ZVEC*U-XVEC*W
  970.       ZNI=XVEC*V-YVEC*U
  971.       GOTO 250
  972.  310  CONTINUE
  973.       U1=XJPL-XP
  974.       V1=YJPL-YP
  975.       W1=ZJPL-ZP
  976.       SCA=U1*XVEC+V1*YVEC+W1*ZVEC
  977.       XC=XJPL-SCA*XVEC
  978.       YC=YJPL-SCA*YVEC
  979.       ZC=ZJPL-SCA*ZVEC
  980.       IF (SCA.GT.RAYJ) CALL ERREUR(40)
  981.       IF (IERR.NE.0) RETURN
  982.       RAYC=SQRT(RAYJ**2-SCA**2)
  983.       U=XC-XP
  984.       V=YC-YP
  985.       W=ZC-ZP
  986.       S=SQRT(U**2+V**2+W**2)
  987.       IF (S.EQ.0.) CALL ERREUR(40)
  988.       IF (IERR.NE.0) RETURN
  989.       U=U/S
  990.       V=V/S
  991.       W=W/S
  992.       XPJ=XC-RAYC*U
  993.       YPJ=YC-RAYC*V
  994.       ZPJ=ZC-RAYC*W
  995.       XNJ=YVEC*W-ZVEC*V
  996.       YNJ=ZVEC*U-XVEC*W
  997.       ZNJ=XVEC*V-YVEC*U
  998.       GOTO 350
  999. C  PROJECTION SUR CYLINDRE AXE XIPL DIR XIV EN  RESTANT DANS XVEC
  1000.  220  CONTINUE
  1001.       U=XIPL-XP
  1002.       V=YIPL-YP
  1003.       W=ZIPL-ZP
  1004.       SCA=U*XIV+V*YIV+W*ZIV
  1005.       U=U-SCA*XIV
  1006.       V=V-SCA*YIV
  1007.       W=W-SCA*ZIV
  1008.       S=SQRT(U**2+V**2+W**2)
  1009.       IF (S.EQ.0.) CALL ERREUR(40)
  1010.       IF (IERR.NE.0) RETURN
  1011.       U1=U-U*RAYI/S
  1012.       V1=V-V*RAYI/S
  1013.       W1=W-W*RAYI/S
  1014.       DNUM=U1*XVEC+V1*YVEC+W1*ZVEC
  1015.       DNOM=XIV*XVEC+YIV*YVEC+ZIV*ZVEC
  1016.       IF (DNOM.EQ.0.) GOTO 221
  1017.       IF (IERR.NE.0) RETURN
  1018.       DPAR=-DNUM/DNOM
  1019.       XPI=XP+U1+DPAR*XIV
  1020.       YPI=YP+V1+DPAR*YIV
  1021.       ZPI=ZP+W1+DPAR*ZIV
  1022.       XNI=(V*ZVEC-W*YVEC)/S
  1023.       YNI=(W*XVEC-U*ZVEC)/S
  1024.       ZNI=(U*YVEC-V*XVEC)/S
  1025.       GOTO 250
  1026. C  L'AXE DU CYLINDRE EST PARALLELE AU PLAN XVEC
  1027.  221  SCA=U*XVEC+V*YVEC+W*ZVEC
  1028.       IF (SCA.GT.RAYI) CALL ERREUR(40)
  1029.       IF (IERR.NE.0) RETURN
  1030.       RAYC=SQRT(RAYI**2-SCA**2)
  1031.       U=U-SCA*XVEC
  1032.       V=V-SCA*YVEC
  1033.       W=W-SCA*ZVEC
  1034.       S=SQRT(U**2+V**2+W**2)
  1035.       IF (S.EQ.0.) CALL ERREUR(40)
  1036.       IF (IERR.NE.0) RETURN
  1037.       XPI=XP+U-U*RAYC/S
  1038.       YPI=YP+V-V*RAYC/S
  1039.       ZPI=ZP+W-W*RAYC/S
  1040.       XNI=XIV
  1041.       YNI=YIV
  1042.       ZNI=ZIV
  1043.       GOTO 250
  1044.  320  CONTINUE
  1045.       U=XJPL-XP
  1046.       V=YJPL-YP
  1047.       W=ZJPL-ZP
  1048.       SCA=U*XJV+V*YJV+W*ZJV
  1049.       U=U-SCA*XJV
  1050.       V=V-SCA*YJV
  1051.       W=W-SCA*ZJV
  1052.       S=SQRT(U**2+V**2+W**2)
  1053.       IF (S.EQ.0.) CALL ERREUR(40)
  1054.       IF (IERR.NE.0) RETURN
  1055.       U1=U-U*RAYJ/S
  1056.       V1=V-V*RAYJ/S
  1057.       W1=W-W*RAYJ/S
  1058.       DNUM=U1*XVEC+V1*YVEC+W1*ZVEC
  1059.       DNOM=XJV*XVEC+YJV*YVEC+ZJV*ZVEC
  1060.       IF (DNOM.EQ.0.) GOTO 321
  1061.       IF (IERR.NE.0) RETURN
  1062.       DPAR=-DNUM/DNOM
  1063.       XPJ=XP+U1+DPAR*XJV
  1064.       YPJ=YP+V1+DPAR*YJV
  1065.       ZPJ=ZP+W1+DPAR*ZJV
  1066.       XNJ=(V*ZVEC-W*YVEC)/S
  1067.       YNJ=(W*XVEC-U*ZVEC)/S
  1068.       ZNJ=(U*YVEC-V*XVEC)/S
  1069.       GOTO 350
  1070. C  L'AXE DU CYLINDRE EST PARALLELE AU PLAN XVEC
  1071.  321  SCA=U*XVEC+V*YVEC+W*ZVEC
  1072.       IF (SCA.GT.RAYJ) CALL ERREUR(40)
  1073.       IF (IERR.NE.0) RETURN
  1074.       RAYC=SQRT(RAYJ**2-SCA**2)
  1075.       U=U-SCA*XVEC
  1076.       V=V-SCA*YVEC
  1077.       W=W-SCA*ZVEC
  1078.       S=SQRT(U**2+V**2+W**2)
  1079.       IF (S.EQ.0.) CALL ERREUR(40)
  1080.       IF (IERR.NE.0) RETURN
  1081.       XPJ=XP+U-U*RAYC/S
  1082.       YPJ=YP+V-V*RAYC/S
  1083.       ZPJ=ZP+W-W*RAYC/S
  1084.       XNJ=XJV
  1085.       YNJ=YJV
  1086.       ZNJ=ZJV
  1087.       GOTO 350
  1088. C  PROJECTION SUR CONE SOMMET XIPL DIR XIV TANG ANGLE SOMM TANGI EN
  1089. C  RESTANT DANS XVEC
  1090.  230  CONTINUE
  1091.       U=XIPL-XP
  1092.       V=YIPL-YP
  1093.       W=ZIPL-ZP
  1094.       SCA=U*XIV+V*YIV+W*ZIV
  1095.       U=U-SCA*XIV
  1096.       V=V-SCA*YIV
  1097.       W=W-SCA*ZIV
  1098.       XC=XP+U
  1099.       YC=YP+U
  1100.       ZC=ZP+U
  1101.       RL2=(SCA*TANGI)**2
  1102.       SCA=U*XVEC+V*YVEC+W*ZVEC
  1103.       DL2=SCA**2
  1104.       IF (DL2.GT.RL2) CALL ERREUR(40)
  1105.       IF(IERR.NE.0) RETURN
  1106.       U=U-SCA*XVEC
  1107.       V=V-SCA*YVEC
  1108.       W=W-SCA*ZVEC
  1109.       S=SQRT(U**2+V**2+W**2)
  1110.       IF (S.EQ.0.) CALL ERREUR(40)
  1111.       IF (IERR.NE.0) RETURN
  1112.       U1=U/S
  1113.       V1=V/S
  1114.       W1=W/S
  1115.       VEC=SQRT(RL2-DL2)
  1116.       U=U-VEC*U1
  1117.       V=V-VEC*V1
  1118.       W=W-VEC*W1
  1119.       XPI=XP+U
  1120.       YPI=YP+V
  1121.       ZPI=ZP+W
  1122.       U=XIPL-XPI
  1123.       V=YIPL-YPI
  1124.       W=ZIPL-ZPI
  1125.       U1=V*ZIV-W*YIV
  1126.       V1=W*XIV-U*ZIV
  1127.       W1=U*YIV-V*XIV
  1128.       U2=V*W1-W*V1
  1129.       V2=W*U1-U*W1
  1130.       W2=U*V1-V*U1
  1131.       XNI=YVEC*W2-ZVEC*V2
  1132.       YNI=ZVEC*U2-XVEC*W2
  1133.       ZNI=XVEC*V2-YVEC*U2
  1134.       S=SQRT(XNI**2+YNI**2+ZNI**2)
  1135.       IF (S.EQ.0.) CALL ERREUR(40)
  1136.       IF (IERR.NE.0) RETURN
  1137.       XNI=XNI/S
  1138.       YNI=YNI/S
  1139.       ZNI=ZNI/S
  1140.       GOTO 250
  1141.  330  CONTINUE
  1142.       U=XJPL-XP
  1143.       V=YJPL-YP
  1144.       W=ZJPL-ZP
  1145.       SCA=U*XJV+V*YJV+W*ZJV
  1146.       U=U-SCA*XJV
  1147.       V=V-SCA*YJV
  1148.       W=W-SCA*ZJV
  1149.       XC=XP+U
  1150.       YC=YP+U
  1151.       ZC=ZP+U
  1152.       RL2=(SCA*TANGJ)**2
  1153.       SCA=U*XVEC+V*YVEC+W*ZVEC
  1154.       DL2=SCA**2
  1155.       IF (DL2.GT.RL2) CALL ERREUR(40)
  1156.       IF(IERR.NE.0) RETURN
  1157.       U=U-SCA*XVEC
  1158.       V=V-SCA*YVEC
  1159.       W=W-SCA*ZVEC
  1160.       S=SQRT(U**2+V**2+W**2)
  1161.       IF (S.EQ.0.) CALL ERREUR(40)
  1162.       IF (IERR.NE.0) RETURN
  1163.       U1=U/S
  1164.       V1=V/S
  1165.       W1=W/S
  1166.       VEC=SQRT(RL2-DL2)
  1167.       U=U-VEC*U1
  1168.       V=V-VEC*V1
  1169.       W=W-VEC*W1
  1170.       XPJ=XP+U
  1171.       YPJ=YP+V
  1172.       ZPJ=ZP+W
  1173.       U=XJPL-XPJ
  1174.       V=YJPL-YPJ
  1175.       W=ZJPL-ZPJ
  1176.       U1=V*ZJV-W*YJV
  1177.       V1=W*XJV-U*ZJV
  1178.       W1=U*YJV-V*XJV
  1179.       U2=V*W1-W*V1
  1180.       V2=W*U1-U*W1
  1181.       W2=U*V1-V*U1
  1182.       XNJ=YVEC*W2-ZVEC*V2
  1183.       YNJ=ZVEC*U2-XVEC*W2
  1184.       ZNJ=XVEC*V2-YVEC*U2
  1185.       S=SQRT(XNJ**2+YNJ**2+ZNJ**2)
  1186.       IF (S.EQ.0.) CALL ERREUR(40)
  1187.       IF (IERR.NE.0) RETURN
  1188.       XNJ=XNJ/S
  1189.       YNJ=YNJ/S
  1190.       ZNJ=ZNJ/S
  1191.       GOTO 350
  1192. C  PROJECTION SUR UN TORE CENTRE XIPL AXE XIAXE CENTRE SECONDAIRE XIP1
  1193. C  GRAND RAYON GRAYI PETIT RAYON PRAYI EN RESTANT DANS XVEC
  1194.  240  CONTINUE
  1195.       U=XP-XIPL
  1196.       V=YP-YIPL
  1197.       W=ZP-ZIPL
  1198.       SCA=U*XIAXE+V*YIAXE+W*ZIAXE
  1199.       U=U-SCA*XIAXE
  1200.       V=V-SCA*YIAXE
  1201.       W=W-SCA*ZIAXE
  1202.       S=SQRT(U**2+V**2+W**2)
  1203.       IF (S.EQ.0.) CALL ERREUR(40)
  1204.       IF (IERR.NE.0) RETURN
  1205. C  CENTRE DU CERCLE INTERS
  1206.       XC=XIPL+U*GRAYI/S
  1207.       YC=YIPL+V*GRAYI/S
  1208.       ZC=ZIPL+W*GRAYI/S
  1209.       U=XP-XC
  1210.       V=YP-YC
  1211.       W=ZP-ZC
  1212.       S=SQRT(U**2+V**2+W**2)
  1213.       IF (S.EQ.0) CALL ERREUR(40)
  1214.       IF (IERR.NE.0) RETURN
  1215. C  PROJECTION SUR LE CERCLE
  1216.       X1=XC+U*PRAYI/S
  1217.       Y1=YC+V*PRAYI/S
  1218.       Z1=ZC+W*PRAYI/S
  1219. C  IL FAUT TOURNER SUR LE TORE JUSQU'AU PLAN XP XVEC
  1220. C  CENTRE DE ROTATION
  1221.       U=X1-XIPL
  1222.       V=Y1-YIPL
  1223.       W=Z1-ZIPL
  1224.       SCA=U*XIAXE+V*YIAXE+W*ZIAXE
  1225.       U=SCA*XIAXE
  1226.       V=SCA*YIAXE
  1227.       W=SCA*ZIAXE
  1228.       XC1=XIPL+U
  1229.       YC1=YIPL+V
  1230.       ZC1=ZIPL+W
  1231. C  RAYON DE CE CERCLE
  1232.       RAY=SQRT((X1-XC1)**2+(Y1-YC1)**2+(Z1-ZC1)**2)
  1233. C  DIRECTION DE L'INTERSECTION DU PLAN DU CERCLE ET DE XP XVEC
  1234.       U=YIAXE*ZVEC-ZIAXE*YVEC
  1235.       V=ZIAXE*XVEC-XIAXE*ZVEC
  1236.       W=XIAXE*YVEC-YIAXE*XVEC
  1237.       S=SQRT(U**2+V**2+W**2)
  1238.       IF (S.EQ.0.) CALL ERREUR(40)
  1239.       IF (IERR.NE.0) RETURN
  1240.       U=U/S
  1241.       V=V/S
  1242.       W=W/S
  1243.       U1=V*ZIAXE-W*YIAXE
  1244.       V1=W*XIAXE-U*ZIAXE
  1245.       W1=U*YIAXE-V*XIAXE
  1246.       U2=XC1-XP
  1247.       V2=YC1-YP
  1248.       W2=ZC1-ZP
  1249.       DNUM=U2*XVEC+V2*YVEC+W2*ZVEC
  1250.       DNOM=U1*XVEC+V1*YVEC+W1*ZVEC
  1251.       IF (DNOM.EQ.0.) CALL ERREUR(40)
  1252.       IF (IERR.NE.0.) RETURN
  1253.       PAR=-DNUM/DNOM
  1254.       XC2=XC1+PAR*U1
  1255.       YC2=YC1+PAR*V1
  1256.       ZC2=ZC1+PAR*W1
  1257.       S=SQRT((XC2-XC1)**2+(YC2-YC1)**2+(ZC2-ZC1)**2)
  1258.       IF (S.GT.RAY) CALL ERREUR(40)
  1259.       IF (IERR.NE.0) RETURN
  1260.       D=SQRT(RAY**2-S**2)
  1261.       IF ((X1-XC2)*U+(Y1-YC2)*V+(Z1-ZC2)*W.LE.0.) D=-D
  1262.       XPI=XC2+D*U
  1263.       YPI=YC2+D*V
  1264.       ZPI=ZC2+D*W
  1265. C  CENTRE DU PETIT CERCLE
  1266.       U=XPI-XIPL
  1267.       V=YPI-YIPL
  1268.       W=ZPI-ZIPL
  1269.       SCA=U*XIAXE+V*YIAXE+W*ZIAXE
  1270.       U=U-SCA*XIAXE
  1271.       V=V-SCA*YIAXE
  1272.       W=W-SCA*ZIAXE
  1273.       S=SQRT(U**2+V**2+W**2)
  1274.       XC1=XIPL+U*GRAYI/S
  1275.       YC1=YIPL+V*GRAYI/S
  1276.       ZC1=ZIPL+W*GRAYI/S
  1277.       U=XPI-XC1
  1278.       V=YPI-YC1
  1279.       W=ZPI-ZC1
  1280.       XNI=V*ZVEC-W*YVEC
  1281.       YNI=W*XVEC-U*ZVEC
  1282.       ZNI=U*YVEC-V*XVEC
  1283.       S=SQRT(XNI**2+YNI**2+ZNI**2)
  1284.       IF (S.EQ.0.) GOTO 250
  1285.       XNI=XNI/S
  1286.       YNI=YNI/S
  1287.       ZNI=ZNI/S
  1288.       GOTO 250
  1289.  340  CONTINUE
  1290.       U=XP-XJPL
  1291.       V=YP-YJPL
  1292.       W=ZP-ZJPL
  1293.       SCA=U*XJAXE+V*YJAXE+W*ZJAXE
  1294.       U=U-SCA*XJAXE
  1295.       V=V-SCA*YJAXE
  1296.       W=W-SCA*ZJAXE
  1297.       S=SQRT(U**2+V**2+W**2)
  1298.       IF (S.EQ.0.) CALL ERREUR(40)
  1299.       IF (IERR.NE.0) RETURN
  1300. C  CENTRE DU CERCLE INTERS
  1301.       XC=XJPL+U*GRAYJ/S
  1302.       YC=YJPL+V*GRAYJ/S
  1303.       ZC=ZJPL+W*GRAYJ/S
  1304.       U=XP-XC
  1305.       V=YP-YC
  1306.       W=ZP-ZC
  1307.       S=SQRT(U**2+V**2+W**2)
  1308.       IF (S.EQ.0) CALL ERREUR(40)
  1309.       IF (IERR.NE.0) RETURN
  1310. C  PROJECTION SUR LE CERCLE
  1311.       X1=XC+U*PRAYJ/S
  1312.       Y1=YC+V*PRAYJ/S
  1313.       Z1=ZC+W*PRAYJ/S
  1314. C  IL FAUT TOURNER SUR LE TORE JUSQU'AU PLAN XP XVEC
  1315. C  CENTRE DE ROTATION
  1316.       U=X1-XJPL
  1317.       V=Y1-YJPL
  1318.       W=Z1-ZJPL
  1319.       SCA=U*XJAXE+V*YJAXE+W*ZJAXE
  1320.       U=SCA*XJAXE
  1321.       V=SCA*YJAXE
  1322.       W=SCA*ZJAXE
  1323.       XC1=XJPL+U
  1324.       YC1=YJPL+V
  1325.       ZC1=ZJPL+W
  1326. C  RAYON DE CE CERCLE
  1327.       RAY=SQRT((X1-XC1)**2+(Y1-YC1)**2+(Z1-ZC1)**2)
  1328. C  DIRECTION DE L'INTERSECTION DU PLAN DU CERCLE ET DE XP XVEC
  1329.       U=YJAXE*ZVEC-ZJAXE*YVEC
  1330.       V=ZJAXE*XVEC-XJAXE*ZVEC
  1331.       W=XJAXE*YVEC-YJAXE*XVEC
  1332.       S=SQRT(U**2+V**2+W**2)
  1333.       IF (S.EQ.0.) CALL ERREUR(40)
  1334.       IF (IERR.NE.0) RETURN
  1335.       U=U/S
  1336.       V=V/S
  1337.       W=W/S
  1338.       U1=V*ZJAXE-W*YJAXE
  1339.       V1=W*XJAXE-U*ZJAXE
  1340.       W1=U*YJAXE-V*XJAXE
  1341.       U2=XC1-XP
  1342.       V2=YC1-YP
  1343.       W2=ZC1-ZP
  1344.       DNUM=U2*XVEC+V2*YVEC+W2*ZVEC
  1345.       DNOM=U1*XVEC+V1*YVEC+W1*ZVEC
  1346.       IF (DNOM.EQ.0.) CALL ERREUR(40)
  1347.       IF (IERR.NE.0.) RETURN
  1348.       PAR=-DNUM/DNOM
  1349.       XC2=XC1+PAR*U1
  1350.       YC2=YC1+PAR*V1
  1351.       ZC2=ZC1+PAR*W1
  1352.       S=SQRT((XC2-XC1)**2+(YC2-YC1)**2+(ZC2-ZC1)**2)
  1353.       IF (S.GT.RAY) CALL ERREUR(40)
  1354.       IF (IERR.NE.0) RETURN
  1355.       D=SQRT(RAY**2-S**2)
  1356.       IF ((X1-XC2)*U+(Y1-YC2)*V+(Z1-ZC2)*W.LE.0.) D=-D
  1357.       XPJ=XC2+D*U
  1358.       YPJ=YC2+D*V
  1359.       ZPJ=ZC2+D*W
  1360. C  CENTRE DU PETIT CERCLE
  1361.       U=XPJ-XJPL
  1362.       V=YPJ-YJPL
  1363.       W=ZPJ-ZJPL
  1364.       SCA=U*XJAXE+V*YJAXE+W*ZJAXE
  1365.       U=U-SCA*XJAXE
  1366.       V=V-SCA*YJAXE
  1367.       W=W-SCA*ZJAXE
  1368.       S=SQRT(U**2+V**2+W**2)
  1369.       XC1=XJPL+U*GRAYJ/S
  1370.       YC1=YJPL+V*GRAYJ/S
  1371.       ZC1=ZJPL+W*GRAYJ/S
  1372.       U=XPJ-XC1
  1373.       V=YPJ-YC1
  1374.       W=ZPJ-ZC1
  1375.       XNJ=V*ZVEC-W*YVEC
  1376.       YNJ=W*XVEC-U*ZVEC
  1377.       ZNJ=U*YVEC-V*XVEC
  1378.       S=SQRT(XNJ**2+YNJ**2+ZNJ**2)
  1379.       IF (S.EQ.0.) GOTO 250
  1380.       XNJ=XNJ/S
  1381.       YNJ=YNJ/S
  1382.       ZNJ=ZNJ/S
  1383.       GOTO 350
  1384.       END
  1385.  
  1386.  
  1387.  
  1388.  
  1389.  
  1390.  
  1391.  
  1392.  
  1393.  

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