Télécharger la plaquette Cast3M


fr en

Télécharger proobj.eso

Retour à la liste

Numérotation des lignes : 

proobj

  1. C PROOBJ    SOURCE    SP204843  24/08/05    21:15:03     11978          
  2. C   CE SOUS-PROGRAMME PREPARE LA PROJECTION D'UN OBJET SUR UNE
  3. C   SURFACE :
  4. C                            PLAN P1 P2 P3
  5. C                            SPHE C P
  6. C                            CYLI C1 C2   P
  7. C                            CONI C C1   P
  8. C                            TORI C A Cs P
  9. C    NOUVELLE POSSIBILITE EN 2D NOVEMBRE 1986
  10. C                            DROI P1 P2
  11. C                            CERCLE P1
  12. C
  13. C    JANVIER 1987  INTRODUCTION DE LA PROJECTION CONIQUE
  14. C
  15. C    SEPTEMBRE 1996 PROJECTION D'UN OBJET MAILLAGE SUR UN OBJET
  16. C    MAILLAGE.
  17.       SUBROUTINE PROOBJ
  18.       IMPLICIT INTEGER(I-N)
  19.       IMPLICIT REAL*8 (A-H,O-Z)
  20.  
  21. -INC PPARAM
  22. -INC CCREEL 
  23. -INC CCOPTIO
  24. -INC SMCOORD
  25. -INC SMELEME
  26.       COMMON /CTOURN/XPT1,YPT1,ZPT1,XV1,YV1,ZV1,XV2,YV2,ZV2,XVEC,YVEC,
  27.      # ZVEC,ANGLE,ICLE,XP1,YP1,ZP1
  28.       CHARACTER*4 MCLE(7),MTYPP(3)
  29.       DATA MTYPP/'DIRE','CONI','CYLI'/
  30.       DATA MCLE/'PLAN','SPHE','CYLI','CONI','TORI','DROI','CERC'/
  31.       SEGACT MCOORD
  32.       MELEME=0
  33.       CALL LIROBJ('MAILLAGE',MELEME,0,IRETOU)
  34.       IF (IRETOU.EQ.1) GOTO 1
  35. C  C'EST UN POINT QU'ON PROJETTE
  36.       CALL LIROBJ('POINT ',IPOIN,1,IRETOU)
  37.       IF (IERR.NE.0) RETURN
  38.    1  CONTINUE
  39. *  CONIQUE OU CYLINDRIQUE
  40.       CALL LIRMOT(MTYPP,3,IYYT,0)
  41.       CALL LIROBJ('POINT ',IP1,1,IRETOU)
  42.       IF (IERR.NE.0) RETURN
  43.       IF (IYYT.EQ.3) THEN
  44.         IP2 = 0
  45.         CALL LIROBJ('POINT ',IP2,1,IRETOU)
  46.         IF (IERR.NE.0) RETURN
  47.       ENDIF
  48.       CALL LIROBJ('MAILLAGE',IPPP,0,IRETO)
  49.       IF (IYYT.EQ.1) THEN
  50.       IF(IRETO.NE.0) THEN
  51.         if (meleme.ne.0) then
  52.           CALL PROOB1(MELEME,IPPP,IP1)
  53.         else
  54.           call erreur(821)
  55.           return
  56.         endif
  57.       RETURN
  58.       ENDIF
  59.       ENDIF
  60.       IF (IYYT.EQ.2) THEN
  61.       IF(IRETO.NE.0) THEN
  62.         CALL PROOB4(MELEME,IPPP,IP1)
  63.       RETURN
  64.       ENDIF
  65.       ENDIF
  66.       SEGACT MCOORD
  67.       IREF=(IDIM+1)*(IP1-1)
  68.       XVEC=XCOOR(IREF+1)
  69.       YVEC=XCOOR(IREF+2)
  70.       ZVEC=XCOOR(IREF+3)
  71.       IF (IDIM.EQ.2) ZVEC=0.
  72.       IF (IYYT.EQ.3) THEN
  73.         IF (IDIM.EQ.2) THEN
  74.           INTERR(1)=2
  75.           CALL ERREUR(709)
  76.           RETURN
  77.         ENDIF
  78. C       write(6,*) 'Option CYL2'
  79.         XP1 = XVEC
  80.         YP1 = YVEC
  81.         ZP1 = ZVEC
  82.         IREF=(IDIM+1)*(IP2-1)
  83.         XP2=XCOOR(IREF+1)
  84.         YP2=XCOOR(IREF+2)
  85.         ZP2=XCOOR(IREF+3)
  86.         IF (IDIM.EQ.2) ZP2=0.
  87.         XVEC=XP2-XP1
  88.         YVEC=YP2-YP1
  89.         ZVEC=ZP2-ZP1
  90.         SVEC=XVEC**2+YVEC**2+ZVEC**2
  91.         UMAX=max(ABS(XP1),ABS(XP2))
  92.         VMAX=max(ABS(YP1),ABS(YP2))
  93.         WMAX=max(ABS(ZP1),ABS(ZP2))
  94.         IF(VMAX.GT.UMAX) UMAX=VMAX
  95.         IF(WMAX.GT.UMAX) UMAX=WMAX
  96.         IF (SVEC.LT.(UMAX*XZPREC)) THEN
  97.           CALL ERREUR(237)
  98.           RETURN
  99.         ENDIF
  100.         SVEC=SQRT(SVEC)
  101.         XVEC=XVEC/SVEC
  102.         YVEC=YVEC/SVEC
  103.         ZVEC=ZVEC/SVEC
  104. C       write(6,*) 'XP1,YP1,ZP1',XP1,YP1,ZP1
  105. C       write(6,*) 'XP2,YP2,ZP2',XP2,YP2,ZP2
  106. C       write(6,*) 'XVEC,YVEC,ZVEC,SVEC',XVEC,YVEC,ZVEC,SVEC
  107.       ENDIF
  108.       IALIR=0
  109.       IF (IYYT.EQ.0) IALIR=1
  110.       CALL LIRMOT(MCLE,7,IVAL,IALIR)
  111.       IF (IVAL.EQ.0) THEN
  112.          IVAL=IYYT+2
  113.          IYYT=0
  114.       ENDIF
  115.       IF (IERR.NE.0) RETURN
  116.       IF (IDIM.EQ.2.AND.IVAL.LE.5) CALL ERREUR(34)
  117.       IF (IDIM.EQ.3.AND.IVAL.GE.6) CALL ERREUR(34)
  118.       IF (IERR.NE.0) RETURN
  119.       ICLE=IVAL+7
  120.       GOTO (10,20,30,40,50,60,70),IVAL
  121.   10  CONTINUE
  122. C  PLAN ON LIT TROIS POINTS
  123.       CALL LIROBJ('POINT ',IP1,1,IRETOU)
  124.       CALL LIROBJ('POINT ',IP2,1,IRETOU)
  125.       CALL LIROBJ('POINT ',IP3,1,IRETOU)
  126.       IF (IERR.NE.0) RETURN
  127.       IREF=(IDIM+1)*(IP1-1)
  128.       XPT1=XCOOR(IREF+1)
  129.       YPT1=XCOOR(IREF+2)
  130.       ZPT1=XCOOR(IREF+3)
  131.       IREF=(IDIM+1)*(IP2-1)
  132.       XV2=XCOOR(IREF+1)-XPT1
  133.       YV2=XCOOR(IREF+2)-YPT1
  134.       ZV2=XCOOR(IREF+3)-ZPT1
  135.       IREF=(IDIM+1)*(IP3-1)
  136.       XV3=XCOOR(IREF+1)-XPT1
  137.       YV3=XCOOR(IREF+2)-YPT1
  138.       ZV3=XCOOR(IREF+3)-ZPT1
  139. C  ON GARDE LE VECTEUR NORMAL NORMALISE
  140.       XV1=YV2*ZV3-ZV2*YV3
  141.       YV1=ZV2*XV3-XV2*ZV3
  142.       ZV1=XV2*YV3-YV2*XV3
  143.       SV1=SQRT(XV1**2+YV1**2+ZV1**2)
  144.       IF (SV1.EQ.0.) CALL ERREUR(21)
  145.       IF (IERR.NE.0) RETURN
  146.       XV1=XV1/SV1
  147.       YV1=YV1/SV1
  148.       ZV1=ZV1/SV1
  149.       GOTO 100
  150.   20  CONTINUE
  151. C   SPHERE ON LIT LE CENTRE ET UN POINT
  152.       CALL LIROBJ('POINT ',IPCEN,1,IRETOU)
  153.       CALL LIROBJ('POINT ',IP1,1,IRETOU)
  154.       IF (IERR.NE.0) RETURN
  155.       IREF=(IPCEN-1)*(IDIM+1)
  156.       XPT1=XCOOR(IREF+1)
  157.       YPT1=XCOOR(IREF+2)
  158.       ZPT1=XCOOR(IREF+3)
  159.       IREF=(IP1-1)*(IDIM+1)
  160.       XV=XCOOR(IREF+1)-XPT1
  161.       YV=XCOOR(IREF+2)-YPT1
  162.       ZV=XCOOR(IREF+3)-ZPT1
  163.       ANGLE=SQRT(XV**2+YV**2+ZV**2)
  164.       IF (ANGLE.EQ.0.) CALL ERREUR(21)
  165.       IF (IERR.NE.0) RETURN
  166.       GOTO 100
  167.   30  CONTINUE
  168. C  CYLINDRE ON LIT DEUX POINTS DE L'AXE ET UN POINT DU CYLINDRE
  169.       CALL LIROBJ('POINT ',IPOIN1,1,IRETOU)
  170.       CALL LIROBJ('POINT ',IPOIN2,1,IRETOU)
  171.       CALL LIROBJ('POINT ',IP1,1,IRETOU)
  172.       IF (IERR.NE.0) RETURN
  173.       IREF=(IDIM+1)*(IPOIN1-1)
  174.       XPT1=XCOOR(IREF+1)
  175.       YPT1=XCOOR(IREF+2)
  176.       ZPT1=XCOOR(IREF+3)
  177.       IREF=(IDIM+1)*(IPOIN2-1)
  178.       XV1=XCOOR(IREF+1)-XPT1
  179.       YV1=XCOOR(IREF+2)-YPT1
  180.       ZV1=XCOOR(IREF+3)-ZPT1
  181.       S=SQRT(XV1**2+YV1**2+ZV1**2)
  182.       IF (S.EQ.0.) CALL ERREUR(21)
  183.       IF (IERR.NE.0) RETURN
  184.       XV1=XV1/S
  185.       YV1=YV1/S
  186.       ZV1=ZV1/S
  187.       IREF=(IDIM+1)*(IP1-1)
  188.       XV2=XCOOR(IREF+1)-XPT1
  189.       YV2=XCOOR(IREF+2)-YPT1
  190.       ZV2=XCOOR(IREF+3)-ZPT1
  191.       XV3=YV1*ZV2-ZV1*YV2
  192.       YV3=ZV1*XV2-XV1*ZV2
  193.       ZV3=XV1*YV2-YV1*XV2
  194.       ANGLE=SQRT(XV3**2+YV3**2+ZV3**2)
  195. C  RAYON DU CYLINDRE
  196.       IF (ANGLE.EQ.0) CALL ERREUR(21)
  197.       IF (IERR.NE.0) RETURN
  198.       GOTO 100
  199.   40  CONTINUE
  200. C  CONE ON LIT LE SOMMET UN POINT DE L'AXE ET UN POINT DU CONE
  201.       CALL LIROBJ('POINT ',IPT1,1,IRETOU)
  202.       CALL LIROBJ('POINT ',IP1,1,IRETOU)
  203.       CALL LIROBJ('POINT ',IP2,1,IRETOU)
  204.       IF (IERR.NE.0) RETURN
  205.       IREF=(IDIM+1)*(IPT1-1)
  206.       XPT1=XCOOR(IREF+1)
  207.       YPT1=XCOOR(IREF+2)
  208.       ZPT1=XCOOR(IREF+3)
  209.       IREF=(IDIM+1)*(IP1-1)
  210.       XV1=XCOOR(IREF+1)-XPT1
  211.       YV1=XCOOR(IREF+2)-YPT1
  212.       ZV1=XCOOR(IREF+3)-ZPT1
  213.       SV1=SQRT(XV1**2+YV1**2+ZV1**2)
  214.       IF (SV1.EQ.0.) CALL ERREUR(21)
  215.       IF (IERR.NE.0) RETURN
  216.       XV1=XV1/SV1
  217.       YV1=YV1/SV1
  218.       ZV1=ZV1/SV1
  219.       IREF=(IDIM+1)*(IP2-1)
  220.       XV2=XCOOR(IREF+1)-XPT1
  221.       YV2=XCOOR(IREF+2)-YPT1
  222.       ZV2=XCOOR(IREF+3)-ZPT1
  223.       SV2=SQRT(XV2**2+YV2**2+ZV2**2)
  224.       IF (SV2.EQ.0.) CALL ERREUR(21)
  225.       IF (IERR.NE.0) RETURN
  226.       XV2=XV2/SV2
  227.       YV2=YV2/SV2
  228.       ZV2=ZV2/SV2
  229.       ANGLE=(XV1*XV2+YV1*YV2+ZV1*ZV2)**2
  230.       IF (IIMPI.NE.0) WRITE (IOIMP,1012) ANGLE
  231.  1012 FORMAT(' COS **2 DE ANGLE AU SOMMET ',G12.5)
  232.       GOTO 100
  233.   50  CONTINUE
  234. C  TORE ON LIT LE CENTRE UN POINT DE L'AXE UN CENTRE SECONDAIRE ET
  235. C  UN POINT
  236.       CALL LIROBJ('POINT ',IPT1,1,IRETOU)
  237.       CALL LIROBJ('POINT ',IPAX,1,IRETOU)
  238.       CALL LIROBJ('POINT ',IPCS,1,IRETOU)
  239.       CALL LIROBJ('POINT ',IP1 ,1,IRETOU)
  240.       IF (IERR.NE.0) RETURN
  241.       IREF=(IDIM+1)*(IPT1-1)
  242.       XPT1=XCOOR(IREF+1)
  243.       YPT1=XCOOR(IREF+2)
  244.       ZPT1=XCOOR(IREF+3)
  245.       IREF=(IDIM+1)*(IPAX-1)
  246.       XV1=XCOOR(IREF+1)-XPT1
  247.       YV1=XCOOR(IREF+2)-YPT1
  248.       ZV1=XCOOR(IREF+3)-ZPT1
  249.       SV1=XV1**2+YV1**2+ZV1**2
  250.       IF (SV1.EQ.0.) CALL ERREUR(21)
  251.       IF (IERR.NE.0) RETURN
  252.       SV1=SQRT(SV1)
  253.       XV1=XV1/SV1
  254.       YV1=YV1/SV1
  255.       ZV1=ZV1/SV1
  256.       IREF=(IDIM+1)*(IPCS-1)
  257.       XV2=XCOOR(IREF+1)-XPT1
  258.       YV2=XCOOR(IREF+2)-YPT1
  259.       ZV2=XCOOR(IREF+3)-ZPT1
  260.       SCA=XV2*XV1+YV2*YV1+ZV2*ZV1
  261.       XPT1=XPT1+SCA*XV1
  262.       YPT1=YPT1+SCA*YV1
  263.       ZPT1=ZPT1+SCA*ZV1
  264.       XV2=XV2-SCA*XV1
  265.       YV2=YV2-SCA*YV1
  266.       ZV2=ZV2-SCA*ZV1
  267.       GR2=XV2**2+YV2**2+ZV2**2
  268.       IF (GR2.EQ.0.) CALL ERREUR(21)
  269.       IF (IERR.NE.0) RETURN
  270.       IREF=(IDIM+1)*(IP1-1)
  271.       XV2=XCOOR(IREF+1)-XPT1
  272.       YV2=XCOOR(IREF+2)-YPT1
  273.       ZV2=XCOOR(IREF+3)-ZPT1
  274.       SCA=XV2*XV1+YV2*YV1+ZV2*ZV1
  275.       XC=XV2-SCA*XV1
  276.       YC=YV2-SCA*YV1
  277.       ZC=ZV2-SCA*ZV1
  278.       SC=XC**2+YC**2+ZC**2
  279.       IF (SC.EQ.0.) CALL ERREUR(21)
  280.       IF (IERR.NE.0) RETURN
  281.       RAP=SQRT(GR2/SC)
  282.       XC=XC*RAP
  283.       YC=YC*RAP
  284.       ZC=ZC*RAP
  285.       PR2=(XV2-XC)**2+(YV2-YC)**2+(ZV2-ZC)**2
  286.       IF (PR2.EQ.0.) CALL ERREUR(21)
  287.       IF (IERR.NE.0) RETURN
  288.       XV2=PR2
  289.       YV2=GR2
  290.       IF (IIMPI.NE.0) WRITE (IOIMP,1015) XV2,YV2,XPT1,YPT1,ZPT1,
  291.      # XV1,YV1,ZV1
  292.  1015 FORMAT(' CARRE PETIT RAYON ',G12.5,
  293.      # ' CARRE GRAND RAYON ',G12.5,/,' CENTRE ',3G13.5,' AXE ',
  294.      #  3G12.5)
  295.       GOTO 100
  296.   60  CONTINUE
  297. C  DROITE ON LIT DEUX POINTS
  298.       CALL LIROBJ('POINT ',IP1,1,IRETOU)
  299.       CALL LIROBJ('POINT ',IP2,1,IRETOU)
  300.       IF (IERR.NE.0) RETURN
  301.       IREF=(IDIM+1)*(IP1-1)
  302.       XPT1=XCOOR(IREF+1)
  303.       YPT1=XCOOR(IREF+2)
  304.       ZPT1=0.
  305.       IREF=(IDIM+1)*(IP2-1)
  306.       XV2=XCOOR(IREF+1)-XPT1
  307.       YV2=XCOOR(IREF+2)-YPT1
  308.       ZV2=0.
  309.       XV3=0.
  310.       YV3=0.
  311.       ZV3=1.
  312. C  ON GARDE LE VECTEUR NORMAL NORMALISE
  313.       XV1=YV2*ZV3-ZV2*YV3
  314.       YV1=ZV2*XV3-XV2*ZV3
  315.       ZV1=XV2*YV3-YV2*XV3
  316.       SV1=SQRT(XV1**2+YV1**2+ZV1**2)
  317.       IF (SV1.EQ.0.) CALL ERREUR(21)
  318.       IF (IERR.NE.0) RETURN
  319.       XV1=XV1/SV1
  320.       YV1=YV1/SV1
  321.       ZV1=ZV1/SV1
  322.       GOTO 100
  323.   70  CONTINUE
  324. C   CERCLE ON LIT LE CENTRE ET UN POINT
  325.       CALL LIROBJ('POINT ',IPCEN,1,IRETOU)
  326.       CALL LIROBJ('POINT ',IP1,1,IRETOU)
  327.       IF (IERR.NE.0) RETURN
  328.       IREF=(IPCEN-1)*(IDIM+1)
  329.       XPT1=XCOOR(IREF+1)
  330.       YPT1=XCOOR(IREF+2)
  331.       ZPT1=0.
  332.       IREF=(IP1-1)*(IDIM+1)
  333.       XV=XCOOR(IREF+1)-XPT1
  334.       YV=XCOOR(IREF+2)-YPT1
  335.       ZV=0.
  336.       ANGLE=SQRT(XV**2+YV**2+ZV**2)
  337.       IF (ANGLE.EQ.0.) CALL ERREUR(21)
  338.       IF (IERR.NE.0) RETURN
  339.       GOTO 100
  340.  100  CONTINUE
  341.       IF (MELEME.EQ.0) GOTO 101
  342.       IF (IYYT.EQ.2) THEN
  343.         IVAL=IVAL+10000
  344.         ICLE=ICLE+10000
  345.       ENDIF
  346.       IF (IYYT.EQ.3) THEN
  347.         IVAL=IVAL-10000
  348.         ICLE=ICLE-10000
  349.       ENDIF
  350.       CALL INTOPE(MELEME)
  351.       RETURN
  352.  101  CONTINUE
  353.       IREF=(IDIM+1)*(IPOIN-1)
  354.       XPOIN=XCOOR(IREF+1)
  355.       YPOIN=XCOOR(IREF+2)
  356.       ZPOIN=XCOOR(IREF+3)
  357.       TPOIN=XCOOR(IREF+4)
  358.       IF (IDIM.EQ.2) ZPOIN=0
  359.       IF (IDIM.EQ.2) TPOIN=XCOOR(IREF+3)
  360.       IF (IYYT.EQ.2) THEN
  361.         XVEC=XVEC-XPOIN
  362.         YVEC=YVEC-YPOIN
  363.         ZVEC=ZVEC-ZPOIN
  364.         IF (IDIM.EQ.2) ZVEC=0.
  365.       ENDIF
  366.       IF (IYYT.EQ.3) THEN
  367.         CALL ERREUR(251)
  368.         RETURN
  369.       ENDIF
  370.       SVEC=SQRT(XVEC**2+YVEC**2+ZVEC**2)
  371.       IF (SVEC.EQ.0.) CALL ERREUR(21)
  372.       IF (IERR.NE.0) RETURN
  373.       XVEC=XVEC/SVEC
  374.       YVEC=YVEC/SVEC
  375.       ZVEC=ZVEC/SVEC
  376.       GOTO (110,120,130,140,150,160,170),IVAL
  377.  110  CONTINUE
  378.       XV=XPOIN-XPT1
  379.       YV=YPOIN-YPT1
  380.       ZV=ZPOIN-ZPT1
  381.       DENUM=XV*XV1+YV*YV1+ZV*ZV1
  382.       DENOM=XVEC*XV1+YVEC*YV1+ZVEC*ZV1
  383.       IF (DENOM.EQ.0.) CALL ERREUR(21)
  384.       IF (IERR.NE.0) RETURN
  385.       RAP=-DENUM/DENOM
  386.       XPOIN=XPOIN+RAP*XVEC
  387.       YPOIN=YPOIN+RAP*YVEC
  388.       ZPOIN=ZPOIN+RAP*ZVEC
  389.       GOTO 200
  390.  120  CONTINUE
  391.       XV=XPT1-XPOIN
  392.       YV=YPT1-YPOIN
  393.       ZV=ZPT1-ZPOIN
  394.       SCA=XVEC*XV+YVEC*YV+ZVEC*ZV
  395.       XV=XVEC*SCA
  396.       YV=YVEC*SCA
  397.       ZV=ZVEC*SCA
  398.       S2=(XPOIN+XV-XPT1)**2+(YPOIN+YV-YPT1)**2+(ZPOIN+ZV-ZPT1)**2
  399.       IF (S2.GT.ANGLE**2) CALL ERREUR(21)
  400.       IF (IERR.NE.0) RETURN
  401.       C=SQRT(ANGLE**2-S2)
  402.       IF (SCA.LT.0.) C=-C
  403.       XPOIN=XPOIN+XV-C*XVEC
  404.       YPOIN=YPOIN+YV-C*YVEC
  405.       ZPOIN=ZPOIN+ZV-C*ZVEC
  406.       GOTO 200
  407.  130  CONTINUE
  408. C  V2 = VEC VECT V1
  409.       XV2=YVEC*ZV1-ZVEC*YV1
  410.       YV2=ZVEC*XV1-XVEC*ZV1
  411.       ZV2=XVEC*YV1-YVEC*XV1
  412.       S2V2=XV2**2+YV2**2+ZV2**2
  413.       IF (S2V2.EQ.0.) CALL ERREUR(21)
  414.       IF (IERR.NE.0) RETURN
  415.       C2=(XVEC*XV1+YVEC*YV1+ZVEC*ZV1)**2
  416.       IF (C2.EQ.1.) CALL ERREUR(21)
  417.       IF (IERR.NE.0) RETURN
  418. C  V3= POIN PT1
  419.       XV3=XPT1-XPOIN
  420.       YV3=YPT1-YPOIN
  421.       ZV3=ZPT1-ZPOIN
  422.       XV4=YV3*ZV1-ZV3*YV1
  423.       YV4=ZV3*XV1-XV3*ZV1
  424.       ZV4=XV3*YV1-YV3*XV1
  425.       DNUM=(XV4*XVEC+YV4*YVEC+ZV4*ZVEC)**2
  426.       DIS2=DNUM/S2V2
  427.       IF (IIMPI.NE.0) WRITE (IOIMP,1000) DIS2,ANGLE
  428.  1000 FORMAT (' DISTANCE AU CARRE DES AXES ',G12.5,
  429.      # 'RAYON DU CYLINDRE ',G12.5)
  430.       IF (DIS2.GT.ANGLE**2) CALL ERREUR(21)
  431.       IF (IERR.NE.0) RETURN
  432.       DMU=SQRT((ANGLE**2-DIS2)/(1.-C2))
  433.       DNUM=XV2*XV4+YV2*YV4+ZV2*ZV4
  434.       DLA=DNUM/S2V2
  435.       DMU=SIGN(DMU,DLA)
  436.       IF (IIMPI.NE.0) WRITE (IOIMP,1001) DLA,DMU
  437.  1001 FORMAT(' LAMBDA,MU ',2G13.5)
  438.       XPOIN=XPOIN+XVEC*(DLA-DMU)
  439.       YPOIN=YPOIN+YVEC*(DLA-DMU)
  440.       ZPOIN=ZPOIN+ZVEC*(DLA-DMU)
  441.       GOTO 200
  442.  140  CONTINUE
  443.       XV2=XPOIN-XPT1
  444.       YV2=YPOIN-YPT1
  445.       ZV2=ZPOIN-ZPT1
  446.       VECV1=XVEC*XV1+YVEC*YV1+ZVEC*ZV1
  447.       VEC2=XVEC**2+YVEC**2+ZVEC**2
  448.       V2V1=XV2*XV1+YV2*YV1+ZV2*ZV1
  449.       VECV2=XVEC*XV2+YVEC*YV2+ZVEC*ZV2
  450.       V22=XV2**2+YV2**2+ZV2**2
  451.       A=VECV1**2-ANGLE*VEC2
  452.       B=2*(V2V1*VECV1-ANGLE*VECV2)
  453.       C=V2V1**2-ANGLE*V22
  454.       DELTA=B**2-4*A*C
  455.       IF (DELTA.LT.0.) CALL ERREUR(21)
  456.       IF (IERR.NE.0) RETURN
  457.       DEL=SQRT(DELTA)
  458.       X1=(-B+DEL)/(2*A)
  459.       X2=(-B-DEL)/(2*A)
  460.       X=X2
  461.       IF (ABS(X1).LT.ABS(X2)) X=X1
  462.       XPOIN=XPOIN+X*XVEC
  463.       YPOIN=YPOIN+X*YVEC
  464.       ZPOIN=ZPOIN+X*ZVEC
  465.       GOTO 200
  466.  150  CONTINUE
  467.       PR2=XV2
  468.       GR2=YV2
  469.       XOP=XPOIN-XPT1
  470.       YOP=YPOIN-YPT1
  471.       ZOP=ZPOIN-ZPT1
  472.       OPV=XOP*XVEC+YOP*YVEC+ZOP*ZVEC
  473.       R2=XOP**2+YOP**2+ZOP**2-GR2-PR2
  474.       VA=XVEC*XV1+YVEC*YV1+ZVEC*ZV1
  475.       QR2VA2=4*GR2*VA**2
  476.       OPA=XOP*XV1+YOP*YV1+ZOP*ZV1
  477.       HR2PV=8*GR2*OPA*VA
  478.       R=4*GR2*OPA**2-4*PR2*GR2
  479.       RLMD=0
  480. C  RESOLUTION DE L'EQUATION DU 4EME DEGRE PAR ITERATION
  481.       DO 151 ITER=1,25
  482.       EXP1=RLMD*(RLMD+2*OPV)+R2
  483.       FDLM=EXP1**2+QR2VA2*RLMD**2+HR2PV*RLMD+R
  484.       FPDLM=4*EXP1*(RLMD+OPV)+QR2VA2*2*RLMD+HR2PV
  485.       IF (FPDLM.EQ.0.) CALL ERREUR(40)
  486.       IF (IERR.NE.0) RETURN
  487.       CORR=FDLM/FPDLM
  488.       IF (IIMPI.NE.0) WRITE(IOIMP,1016) ITER,RLMD,CORR
  489.  1016 FORMAT(' ITER ',I6,' LAMBDA ',G12.5,' CORR ',G12.5)
  490.       RLMD=RLMD-CORR
  491.       IF (RLMD.EQ.0.) GOTO 151
  492.       IF (ABS(CORR/RLMD).LT.1E-5) GOTO 152
  493.  151  CONTINUE
  494.       CALL ERREUR(40)
  495.       RETURN
  496.  152  CONTINUE
  497.       XPOIN=XPOIN+RLMD*XVEC
  498.       YPOIN=YPOIN+RLMD*YVEC
  499.       ZPOIN=ZPOIN+RLMD*ZVEC
  500.       GOTO 200
  501.  160  CONTINUE
  502.       XV=XPOIN-XPT1
  503.       YV=YPOIN-YPT1
  504.       ZV=ZPOIN-ZPT1
  505.       DENUM=XV*XV1+YV*YV1+ZV*ZV1
  506.       DENOM=XVEC*XV1+YVEC*YV1+ZVEC*ZV1
  507.       IF (DENOM.EQ.0.) CALL ERREUR(21)
  508.       IF (IERR.NE.0) RETURN
  509.       RAP=-DENUM/DENOM
  510.       XPOIN=XPOIN+RAP*XVEC
  511.       YPOIN=YPOIN+RAP*YVEC
  512.       ZPOIN=ZPOIN+RAP*ZVEC
  513.       GOTO 200
  514.  170  CONTINUE
  515.       XV=XPT1-XPOIN
  516.       YV=YPT1-YPOIN
  517.       ZV=ZPT1-ZPOIN
  518.       SCA=XVEC*XV+YVEC*YV+ZVEC*ZV
  519.       XV=XVEC*SCA
  520.       YV=YVEC*SCA
  521.       ZV=ZVEC*SCA
  522.       S2=(XPOIN+XV-XPT1)**2+(YPOIN+YV-YPT1)**2+(ZPOIN+ZV-ZPT1)**2
  523.       IF (S2.GT.ANGLE**2) CALL ERREUR(21)
  524.       IF (IERR.NE.0) RETURN
  525.       C=SQRT(ANGLE**2-S2)
  526.       IF (SCA.LT.0.) C=-C
  527.       XPOIN=XPOIN+XV-C*XVEC
  528.       YPOIN=YPOIN+YV-C*YVEC
  529.       ZPOIN=ZPOIN+ZV-C*ZVEC
  530.       GOTO 200
  531.  200  CONTINUE
  532.       SEGACT MCOORD*mod
  533.       NBPTS=nbpts+1         
  534.       SEGADJ MCOORD
  535.       XCOOR((NBPTS-1)*(IDIM+1)+1)=XPOIN
  536.       XCOOR((NBPTS-1)*(IDIM+1)+2)=YPOIN
  537.       IF (IDIM.EQ.3) XCOOR((NBPTS-1)*(IDIM+1)+3)=ZPOIN
  538.       XCOOR(NBPTS*(IDIM+1))=TPOIN
  539.       CALL ECROBJ('POINT ',NBPTS)
  540.       RETURN
  541.       END
  542.  
  543.  
  544.  
  545.  
  546.  
  547.  
  548.  

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