C PROJEC    SOURCE    PV        20/03/30    21:22:57     10567          
      SUBROUTINE PROJEC(ICPR,XPROJ,IOEIL,CGRAV,axez)
      IMPLICIT INTEGER(I-N)
      REAL*8 XO,XG,XP,XN,SN,XV,SV,UI,UJ
      DIMENSION XO(3),XP(3),XN(3),XG(3),XV(3),UI(3),UJ(3),CGRAV(*),
     >  axez(*)

-INC PPARAM
-INC CCOPTIO
-INC SMCOORD
      SEGMENT ICPR(nbpts)
      SEGMENT XPROJ(3,ITE)
      SEGACT MCOORD

      IF (IDIM .NE. 2) GOTO 5500
C       REPERE LOCAL SUR LE PLAN
        XN(1)=0.D0
        XN(2)=0.D0
        XN(3)=1.D0
        
        UJ(1)=0.D0
        UJ(2)=0.D0
        UJ(3)=1.D0

      DO 5501 I=1,ICPR(/1)
        IF (ICPR(I).EQ.0) GOTO 5501
        XPROJ(1,ICPR(I))=XCOOR(I*3-2)
        XPROJ(2,ICPR(I))=XCOOR(I*3-1)
 5501 CONTINUE
      GOTO 5502

C Seulement si IDIM .NE. 2
 5500 CONTINUE
      IREF=(IOEIL-1)*4
      XO(1)=XCOOR(IREF+1)
      XO(2)=XCOOR(IREF+2)
      XO(3)=XCOOR(IREF+3)
C POINT MOYEN
      DO 1 I=1,3
      XG(I)=0.D0
    1 CONTINUE
      DO 2 I=1,ICPR(/1)
      IF (ICPR(I).EQ.0) GOTO 2
      DO 3 J=1,3
      XG(J)=XG(J)+XCOOR(I*4-4+J)
   3  CONTINUE
    2 CONTINUE
      NBPOIN=XPROJ(/2)
      DO 4 J=1,3
      XG(J)=XG(J)/NBPOIN
      XN(J)=XO(J)-XG(J)
    4 CONTINUE
C  NORMALE AU PLAN
      SN=SQRT(XN(1)**2+XN(2)**2+XN(3)**2)
      IF (SN.EQ.0.D0) CALL ERREUR(21)
      IF (IERR.NE.0) RETURN
      DO 5 J=1,3
       XN(J)=XN(J)/SN
  5   CONTINUE
C  REPERE LOCAL SUR LE PLAN
      UJ(1)=0.D0
      UJ(2)=0.D0
      UJ(3)=1.D0
  21  CONTINUE
      SV=UJ(1)*XN(1)+UJ(2)*XN(2)+UJ(3)*XN(3)
      DO 20 J=1,3
      UJ(J)=UJ(J)-SV*XN(J)
  20  CONTINUE
      SV=UJ(1)**2+UJ(2)**2+UJ(3)**2
      IF (ABS(SV).LT.1.D-2) THEN
      UJ(1)=0.D0
      UJ(2)=1.D0
      UJ(3)=1.D0
      GOTO 21
      ENDIF
      SV=SQRT(SV)
      UJ(1)=UJ(1)/SV
      UJ(2)=UJ(2)/SV
      UJ(3)=UJ(3)/SV
      UI(1)=UJ(2)*XN(3)-UJ(3)*XN(2)
      UI(2)=UJ(3)*XN(1)-UJ(1)*XN(3)
      UI(3)=UJ(1)*XN(2)-UJ(2)*XN(1)
C PROJECTION CONIQUE SUR LE PLAN
      DO 12 I=1,ICPR(/1)
      IF (ICPR(I).EQ.0) GOTO 12
      DO 13 J=1,3
      XP(J)=XCOOR(I*4-4+J)
      XV(J)=XP(J)-XO(J)
  13  CONTINUE
*     XPROJ(3,ICPR(I))=SQRT(XV(1)**2+XV(2)**2+XV(3)**2)
      SV=XV(1)*XN(1)+XV(2)*XN(2)+XV(3)*XN(3)
      SN=(XP(1)-XG(1))*XN(1)+(XP(2)-XG(2))*XN(2)+(XP(3)-XG(3))*XN(3)
      XPROJ(3,ICPR(I))=-SN
      DO 14 J=1,3
      XP(J)=XP(J)-(SN/SV)*XV(J)-XG(J)
  14  CONTINUE
      XPROJ(1,ICPR(I))=XP(1)*UI(1)+XP(2)*UI(2)+XP(3)*UI(3)
      XPROJ(2,ICPR(I))=XP(1)*UJ(1)+XP(2)*UJ(2)+XP(3)*UJ(3)
  12  CONTINUE
*  rendre le centre de gravite pour eventuelle rotation
      cgrav(1)=xg(1)
      cgrav(2)=xg(2)
      cgrav(3)=xg(3)

 5502 CONTINUE
*  axez pour tourner correctement avec opengl
      axez(1)=0
      axez(2)=uj(3)
      axez(3)=sqrt(ABS(1-uj(3)**2))
      if (xn(3).lt.0.D0) axez(3)=-axez(3)
*     write (6,*) ' axez ',axez(1),axez(2),axez(3)
      RETURN
      END