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

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