C BIOSAV    SOURCE    CHAT      05/11/29    21:15:02     5153           
      SUBROUTINE BIOSAV
*----------------------------------------------------------
*
*           CALCUL DES CHAMPS DE BIOT ET SAVART
*
*----------------------------------------------------------
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
-INC CCOPTIO
-INC CCREEL
-INC SMELEME
-INC SMCOORD
-INC TMTRAV
-INC SMCHPOI
      INTEGER ENT 
      DIMENSION C(3),PN1(3),RESP(3)
      DIMENSION VX(3),VZ(3),VY(3),VI(3),PP3(3)
      DIMENSION ZF1(1),ZALF(1),ZDENS(1),ZBET(1)
      DIMENSION P1(3),P2(3),P3(3),ROTA(3,3)
      SEGMENT ICPR(XCOOR(/1)/(IDIM+1))
      SEGMENT SANGLE 
      REAL*8 TETM(NT4)
      REAL*8 TETI(NT4)
      REAL*8 DTE(NT4)
      ENDSEGMENT
*
      CHARACTER*4 MCLE(4),CHSECT(1),CHOIX(2)
      LOGICAL INDU,POTE
*
      DATA EPSIL3/1.D-12/,EPSIL2/1.D-8/
*
      DATA MCLE/'CERC','ARC ','BARR','FIL '/
      DATA CHSECT/'TRAP'/
      DATA CHOIX/'INDU','POTE'/
*
      CALL LIROBJ('CHPOINT  ',MCHPO1,0,IRETOU)
      IF (IERR.NE.0) RETURN
      IF (MCHPO1.GT.0) THEN
        CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU)
        IF (IERR.NE.0) RETURN

* OPTION CALCUL PAR INTEGRALE ELLIPTIQUE (PHR 2002)

        CALL BIOKAM (MCHPO1,IPT1,ENT)
        IF (IERR.NE.0) RETURN
        MCHPOI = ENT
        GOTO 9990
      ENDIF
*   CHOIX DE L'OPTION DE CALCUL
*   (INDUCTION PAR DEFAUT)
*
      INDU=.TRUE.
      POTE=.FALSE.
      CALL LIRMOT(CHOIX,2,IMC,0)
      IF (IERR.NE.0) RETURN
      IF(IMC.EQ.2) THEN
         POTE=.TRUE.
         INDU=.FALSE.
      ENDIF
*
*   LECTURE D'UN MAILLAGE
*
      CALL LIROBJ('MAILLAGE',MELEME,1,IRETOU)
      IF (IERR.NE.0) RETURN
*
*   LECTURE D'UN MOT CLE
*
      CALL LIRMOT(MCLE,4,IMC,1)
      IF (IERR.NE.0) RETURN
*
*         IMC = 1  ON A LU CERCLE
*         IMC = 2  ON A LU ARC
*         IMC = 3  ON A LU BARRE
*         IMC = 4  ON A LU FIL
*
*
*   LECTURE DE 2 ou 3 POINTS
*
      CALL LIROBJ('POINT ',IP1,1,IRETOU)
      CALL LIROBJ('POINT ',IP2,1,IRETOU)
      IP3=0
      IF(IMC.LT.4) THEN
         CALL LIROBJ('POINT ',IP3,1,IRETOU)
      ELSE
C        CAS DU FIL RECTILIGNE
         GO TO 10
      ENDIF
      IF (IERR.NE.0) RETURN
*
*  LECTURE DE 4 OU 5 FLOTTANTS SELON LES CAS
*
      CALL LIRREE(F1,1,IRETOU)
      CALL LIRREE(F2,1,IRETOU)
      F3=0.D0
      IF(IMC.LT.3) CALL LIRREE(F3,1,IRETOU)
*
*  LECTURE DU TYPE DE SECTION
*
      ISECT=0
      CALL LIRMOT(CHSECT,1,ISECT,0)
      IF(ISECT.EQ.1) THEN
         CALL LIRREE(PENT1,1,IRETOU)
         CALL LIRREE(PENT2,1,IRETOU)
         IF(IERR.NE.0) RETURN
      ELSE
         ISECT=3
      ENDIF
 10   CONTINUE
      CALL LIRREE(DENS,1,IRETOU)
      CALL LIRREE(XIU,1,IRETOU)
C   ON VIENT DE LIRE  MU0 DANS L UNITE UTILISATEUR
      RUMKSA= 4.D0 * XPI * 1.D-7
      RAP=(XIU / RUMKSA)
      IU= 1
      IF (IERR.NE.0) RETURN
      ZF1(1)=F1
      ZDENS(1)=DENS
*
*  EXEMPLE DE RECUPERATION DES COORDONNEES DU POINT P1
*
      IREF=(IDIM+1)*IP1
      P1(1)=XCOOR(IREF-IDIM)
      P1(2)=XCOOR(IREF-IDIM+1)
      P1(3)=0.D0
      IF (IDIM.EQ.3) P1(3)=XCOOR(IREF-IDIM+2)
      IREF=(IDIM+1)*IP2
      P2(1)=XCOOR(IREF-IDIM)
      P2(2)=XCOOR(IREF-IDIM+1)
      P2(3)=0.D0
      IF (IDIM.EQ.3) P2(3)=XCOOR(IREF-IDIM+2)
      IF(IP3.NE.0) THEN
         IREF=(IDIM+1)*IP3
         P3(1)=XCOOR(IREF-IDIM)
         P3(2)=XCOOR(IREF-IDIM+1)
         P3(3)=0.D0
         IF (IDIM.EQ.3) P3(3)=XCOOR(IREF-IDIM+2)
      ENDIF
*
*                -----------------------
*                |   LES CALCULS ....  |
*                -----------------------
*
*    ON COMMENCE PAR EXPLORER LE MAILLAGE
*
      SEGACT MELEME
      IF (ITYPEL.NE.1.OR.LISOUS(/1).NE.0) CALL CHANGE(MELEME,1)
*
*    ON VERIFIE QUE CHAQUE POINT N'APPARAIT QU'UNE SEULE FOIS
*
      NNNOE=NUM(/2)
      SEGINI ICPR
      ICON=0
      DO 12 I=1,NNNOE
         IKI=NUM(1,I)
         IF(ICPR(IKI).NE.0) GO TO 12
         ICON=ICON+1
         ICPR(IKI)=ICON
 12   CONTINUE
      SEGSUP ICPR
*
      IF(ICON.NE.NNNOE) THEN
         CALL ERREUR(75)
         SEGDES MELEME
         RETURN
      ENDIF
*
*    ON INITIALISE LES NOMS DES COMPOSANTES ET LES HARMONIQUES
*           POUR L'INSTANT  LIMITE AU 3D ET 2D  PLAN
*
*
      NNIN=IDIM
      SEGINI MTRAV
      NHRM=NIFOUR
      IF(NIFOUR.EQ.1) NHRM=IFOUR
      INCO(1)='BX  '
      IF(POTE) INCO(1)='AX  '
      INCO(2)='BY  '
      IF(POTE) INCO(2)='AY  '
      NHAR(1)=NHRM
      NHAR(2)=NHRM
      IF(IDIM.EQ.3) THEN
         INCO(3)='BZ  '
         IF(POTE) INCO(3)='AZ  '
         NHAR(3)=NHRM
      ENDIF
C
      GOTO (100,100,200,300),IMC
 100  CONTINUE
C          BOBINE CIRCULAIRE

C                  ETABLISSEMENT DU REPERE LOCAL
C   AXE DES Z   (P2 - P1 ) * (P3 - P1)
      DO  801 I=1,3
         VX(I)=P2(I)- P1(I)
         VI(I)=P3(I)- P1(I)
 801  CONTINUE
      CALL PROVEC(VX,VI,VZ)
      CALL PROVEC(VZ,VX,VY)
      CALL ROTREP (ROTA,VX,VY,VZ)
      IF(IERR.NE.0) THEN
         SEGDES MELEME
         SEGSUP MTRAV
         RETURN
      ENDIF
      NB=1
      IF(F1.EQ.0.D0)  RETURN
      ALF=F2/F1
      BET=F3/(2.D0*F1)
      ZALF(1)=ALF
      ZBET(1)=BET
      IC= 3
      DO  802 J=1,3
         PP3(J)= P3(J) - P1(J)
 802  CONTINUE
      RAYO2=SDOT(3,VX,1,VX,1)
      PROD=SDOT(3,VX,1,VI,1)
      RAYO2=RAYO2+EPSIL3
      TM=ACOS(PROD/RAYO2)
      GOTO  350
 200  CONTINUE
C   CAS   DES BARRES ON VA TOURNER / OX POUR AMENER OY //  P1P2
      GY=0.D0
      DO  803 I=1,3
         VX(I)=P2(I)- P1(I)
         VI(I)=P3(I)- P1(I)
         GY=GY +( VX(I)*VX(I))
         C(I)= (P2(I)+ P1(I))/2.D0
 803  CONTINUE
      GY= SQRT( GY)
      CALL PROVEC(VX,VI,VZ)
      CALL PROVEC(VZ,VX,VY)
      CALL ROTREP (ROTA,VX,VY,VZ)
      IF(IERR.NE.0) THEN
         SEGDES MELEME
         SEGSUP MTRAV
         RETURN
      ENDIF
      DO 804 I= 1,3
         VZ(I)=0.D0
         DO 8041 J= 1,3
            VZ(I)= VZ(I) + ROTA(I,J) * C(J)
 8041    CONTINUE
 804  CONTINUE
*      CALL ROTATI (C,1,1,1)
      GO TO 350
C
C     CAS DU FIL
 300  CONTINUE
C
 350  CONTINUE
      UPSI =1.D-9
*
*    ON BOUCLE SUR LES NOEUDS
*
      NT4=8000
      SEGINI SANGLE
      DO 30 IPT=1,NNNOE
         L=NUM(1,IPT)
         IREF= (IDIM+1)*L
*   L EST LE NUMERO DU NOEUD  DANS LE TABLEAU DES COORDONNEES
         DO  805 J=1,3
            PN1(J)=XCOOR(IREF-IDIM-1+J)
 805     CONTINUE
*
         GOTO (400,400,500,600),IMC
 400     CONTINUE
*   CAS DES CERCLES ET DES ARC
C   MESURE DU POINT DANS LE REPERE LOCAL
         DO  806 J=1,3
            PN1(J)=PN1(J) - P1(J)
 806     CONTINUE
         DO 807 I= 1,3
            VZ(I)=0.D0
            VY(I)=0.D0
            DO 8071 J= 1,3
               VY(I)= VY(I) + ROTA(I,J) * PN1(J)
               VZ(I)= VZ(I) + ROTA(I,J) * PP3(J)
 8071       CONTINUE
 807     CONTINUE
*      CALL ROTATI (PP3,1,1,1)
*      CALL ROTATI (PN1,1,1,1)
C                       ON  CHANGE DE RERERENTIEL
C     MESURE DANS LE REPERE LOCAL DES POINTS DE CALCUL DE HS

         ZP= VY(3)
         RP= SQRT ( VY(1)**2 + VY(2)**2 )
         IF(RP.EQ.0.D0) THEN
            TMIN=0.D0
            COST=1.D0
            SINT=0.D0
            RP=EPSIL2
         ELSE
            TMIN=ACOS(VY(1)/(RP+EPSIL3))
            COST=COS(TMIN)
            SINT=VY(2)/(RP+EPSIL3)
         ENDIF
         DO 808 J= 1,3
            C(J)=0.D0
 808     CONTINUE
C
         IF(BET.EQ.0.D0.AND.ALF.EQ.1.D0) THEN
C       CAS DU FIL CIRCULAIRE
            RAYON=F1
            COUR=DENS
C       ON CONSTRUIT LE REPERE LOCAL
            sigt=1.d0
            if( abs(sint).gt.1.e-6)SIGT=SIGN(1.D0,-SINT)
            TMIN=SIGT*(TMIN-XPI)+XPI
            IF(IMC.EQ.1) TM=2.D0*XPI
            TMAX=TMIN+TM
            CALL TEGRAL(SANGLE)
            IF(INDU) THEN
               CALL BFILCI(SANGLE,COUR,RAYON,TMIN,TMAX,RP,ZP,BX,BY,BZ)
            ENDIF
            IF(POTE) THEN
               CALL AFILCI(SANGLE,COUR,RAYON,TMIN,TMAX,RP,ZP,BX,BY,BZ)
            ENDIF
            RESP(1)=BX*COST-BY*SINT
            RESP(2)=BX*SINT+BY*COST
            RESP(3)=BZ
         ELSE IF(BET.EQ.0.D0.AND.ISECT.EQ.3) THEN
C       CAS DE LA PLAQUE CIRCULAIRE
            RINT=F1
            REXT=F2
            DR=REXT-RINT
            COUR=DENS*DR
C       ON CONSTRUIT LE REPERE LOCAL
            sigt=1.d0
            if( abs(sint).gt.1.e-6)  SIGT=SIGN(1.D0,-SINT)
            TMIN=SIGT*(TMIN-XPI)+XPI
            IF(IMC.EQ.1) TM=2.D0*XPI
            TMAX=TMIN+TM
            CALL TEGRAL(SANGLE)
            IF(INDU) THEN
               CALL BPLQCI(SANGLE,COUR,RINT,REXT,TMIN,TMAX,RP,ZP,
     &              BX,BY,BZ)
            ENDIF
            IF(POTE) THEN
               CALL APLQCI(SANGLE,COUR,RINT,REXT,TMIN,TMAX,RP,ZP,
     &              BX,BY,BZ)
            ENDIF
            RESP(1)=BX*COST-BY*SINT
            RESP(2)=BX*SINT+BY*COST
            RESP(3)=BZ
         ELSE IF(ALF.EQ.1D0) THEN
C       CAS DU CYLINDRE
            RAYON=F1
            H=F3
            COUR=DENS*H
C       ON CONSTRUIT LE REPERE LOCAL
            sigt=1.d0
            if( abs(sint).gt.1.e-6)SIGT=SIGN(1.D0,-SINT)
            TMIN=SIGT*(TMIN-XPI)+XPI
            IF(IMC.EQ.1) TM=2.D0*XPI
            TMAX=TMIN+TM
            CALL TEGRAL(SANGLE)
            IF(INDU) THEN
               CALL BCYLCI(SANGLE,COUR,RAYON,H,TMIN,TMAX,RP,ZP,
     &              BX,BY,BZ)
            ENDIF
            IF(POTE) THEN
               CALL ACYLCI(SANGLE,COUR,RAYON,H,TMIN,TMAX,RP,ZP,
     &              BX,BY,BZ)
            ENDIF
            RESP(1)=BX*COST-BY*SINT
            RESP(2)=BX*SINT+BY*COST
            RESP(3)=BZ
         ELSE
            IF (IMC.EQ.1.AND.ISECT.EQ.3) THEN
               IF(INDU) THEN
                  CALL CHAMBO(IU,IC,NB,C(3),ZF1,ZALF,ZBET,ZDENS,ZP,RP,BZ
     $                 ,BR,BP)
                  IF (RP.LE.UPSI) THEN
                     RESP(1)=0.D0
                     RESP(2)=0.D0
                  ELSE
                     RESP(1)=BR * VY(1) / RP
                     RESP(2)=BR * VY(2) / RP
                  END IF
                  RESP(3)=BZ
                  GO TO 700
               ENDIF
               IF(POTE) THEN
                  RINT=F1
                  REXT=F2
                  HCEN=F3
                  WIDTH=F2-F1
                  AIRE=F3*WIDTH
                  COUR=DENS*AIRE
                  DELTA=0.D0
C         ON CONSTRUIT LE REPERE LOCAL
                  SIGT=1.D0
                  if(ABS(SINT).gt.1.e-8)SIGT=SIGN(1.D0,-SINT)
                  TMIN=SIGT*(TMIN-XPI)+XPI
                  TM=2.D0*XPI
                  TMAX=TMIN+TM
                  PENT1=0.D0
                  PENT2=0.D0
                  CALL TEGRAL(SANGLE)
                  CALL AARCTR(SANGLE,COUR,HCEN,HCEN,RINT,REXT,
     *                 PENT1,PENT2,TMIN,TMAX,
     *                 RP,ZP,BX,BY,BZ)
                  RESP(1)=BX*COST-BY*SINT
                  RESP(2)=BX*SINT+BY*COST
                  RESP(3)=BZ
                  GO TO 700
               ENDIF
            ENDIF
            IF(ISECT.EQ.3) THEN
C       CAS DES SECTIONS RECTANGULAIRES
C       ON CALCULE POUR UN ARC  ---> ANGLES
               A1=0.D0
               IF (ABS(VZ(1)).LT.UPSI) THEN
                  A2=XPI/2.D0
               ELSE
                  AA= VZ(2)/VZ(1)
                  A2=ATAN ( AA )
                  IF ( A2.LT.0.D0) A2 = XPI + A2
               END IF
               IF(INDU) THEN
                  CALL CHAMAR(IU,IC,C,A1,A2,F1,ALF,BET,DENS,
     *                 VY(1),VY(2),VY(3),BX,BY,BZ)
                  RESP(1)=BX
                  RESP(2)=BY
                  RESP(3)=BZ
               ENDIF
               IF(POTE) THEN
                  RINT=F1
                  REXT=F2
                  HCEN=F3
                  WIDTH=F2-F1
                  AIRE=F3*WIDTH
                  COUR=DENS*AIRE
                  PENT1=0.D0
                  PENT2=0.D0
C         ON CONSTRUIT LE REPERE LOCAL
                  sigt=1.d0
                  if( abs(sint).gt.1.e-6)SIGT=SIGN(1.d0,-SINT)
                  TMIN=SIGT*(TMIN-XPI)+XPI
                  TMAX=TMIN+TM
                  CALL TEGRAL(SANGLE)
                  CALL AARCTR(SANGLE,COUR,HCEN,HCEN,RINT,REXT,
     *                 PENT1,PENT2,TMIN,TMAX,
     *                 RP,ZP,BX,BY,BZ)
                  RESP(1)=BX*COST-BY*SINT
                  RESP(2)=BX*SINT+BY*COST
                  RESP(3)=BZ
               ENDIF
            ELSEIF(ISECT.EQ.1) THEN
C     CAS DES SECTIONS TRAPEZOIDALES
               RINT=F1
               REXT=F2
               HCEN=F3
               WIDTH=F2-F1
               AIRE=F3*WIDTH
               COUR=DENS*AIRE
               DELTA=0.5D0*WIDTH*(PENT2-PENT1)
               HINT=HCEN-DELTA
               HEXT=HCEN+DELTA
C       ON CONSTRUIT LE REPERE LOCAL
               SIGT=1.d0
               if(abs(sint).gt.1.e-06)SIGT=SIGN(1.D0,-SINT)      
               TMIN=SIGT*(TMIN-XPI)+XPI
               IF(IMC.EQ.1) TM=2.D0*XPI
               TMAX=TMIN+TM
               CALL TEGRAL(SANGLE)
               IF(HCEN.EQ.0.D0) THEN
C         CAS DU TRONC DE CONE
                  TETA=ATAN(PENT1)
                  COUR=DENS*WIDTH/COS(TETA)
                  IF(INDU) THEN
                     CALL BTRCON(SANGLE,COUR,RINT,REXT,PENT1,TMIN,TMAX,
     &                    RP,ZP,BX,BY,BZ)
                  ENDIF
                  IF(POTE) THEN
                     CALL ATRCON(SANGLE,COUR,RINT,REXT,PENT1,TMIN,TMAX,
     &                    RP,ZP,BX,BY,BZ)
                  ENDIF
               ELSE
                  IF(INDU) THEN
                     CALL ARCTRA(SANGLE,COUR,HINT,HEXT,RINT,REXT,
     *                    PENT1,PENT2,TMIN,TMAX,
     *                    RP,ZP,BX,BY,BZ)
                  ENDIF
                  IF(POTE) THEN
                     CALL AARCTR(SANGLE,COUR,HINT,HEXT,RINT,REXT,
     *                    PENT1,PENT2,TMIN,TMAX,
     *                    RP,ZP,BX,BY,BZ)
                  ENDIF
               ENDIF
               RESP(1)=BX*COST-BY*SINT
               RESP(2)=BX*SINT+BY*COST
               RESP(3)=BZ
            ENDIF
         END IF
         IF(IERR.NE.0) THEN
            SEGDES MELEME
            SEGSUP MTRAV
            RETURN
         ENDIF
C
         GOTO   700
*
 500     CONTINUE
C      CAS DES BARRES         ON  CHANGE DE RERERENTIEL
         DO 809 I= 1,3
            VY(I)=0.D0
            DO 8091 J= 1,3
               VY(I)= VY(I) + ROTA(I,J) * PN1(J)
 8091       CONTINUE
 809     CONTINUE
*      CALL ROTATI (PN1,1,1,1)
C     MESURE DANS LE REPERE LOCAL DES POINTS DE CALCUL DE HS
         IF(ISECT.EQ.3) THEN
C     CAS DES SECTIONS RETANGULAIRES
            IF(F2.EQ.0.D0) THEN
C         CAS DE LA PLAQUE RECTANGULAIRE
               DO 813 I=1,3
                  VY(I)=0.D0
                  DO 8131 J=1,3
                     VY(I)=VY(I)+ROTA(I,J)*(PN1(J)-P1(J))
 8131             CONTINUE
 813           CONTINUE
               DY=F1
               REXT=0.5D0*DY
               RINT=-REXT
               COUR=DENS*DY
               YMIN=-VY(1)
               YMAX=YMIN+GY
               XP=-VY(2)
               ZP=VY(3)
               IF(INDU) THEN
                  CALL BPLQDR(COUR,RINT,REXT,YMIN,YMAX,XP,ZP,BXP,BYP,BZP
     $                 )
               ENDIF
               IF(POTE) THEN
                  CALL APLQDR(COUR,RINT,REXT,YMIN,YMAX,XP,ZP,BXP,BYP,BZP
     $                 )
               ENDIF
               RESP(1)=BYP
               RESP(2)=-BXP
               RESP(3)= BZP
            ELSE
               IF(INDU) THEN
                  CALL CHABAR(IU,VZ,GY,F1,F2,DENS,VY(1),VY(2),VY(3),
     *                 BXP,BYP,BZP)
                  RESP(1)=BXP
                  RESP(2)=BYP
                  RESP(3)=BZP
               ENDIF
               IF(POTE) THEN
                  DO 814 I=1,3
                     VY(I)=0.D0
                     DO 8141 J=1,3
                        VY(I)=VY(I)+ROTA(I,J)*(PN1(J)-P1(J))
 8141                CONTINUE
 814              CONTINUE
                  DY=F1
                  DZ=F2
                  REXT=0.5D0*DY
                  RINT=-REXT
                  AIRE=DY*DZ
                  COUR=DENS*AIRE
                  PENT1=0.D0
                  PENT2=0.D0
                  YMIN=-VY(1)
                  YMAX=YMIN+GY
                  XP=-VY(2)
                  ZP=VY(3)
                  CALL ASEGTR(COUR,DZ,DZ,RINT,REXT,PENT1,PENT2,
     *                 YMIN,YMAX,XP,ZP,BXP,BYP,BZP)
                  RESP(1)=BYP
                  RESP(2)=-BXP
                  RESP(3)= BZP
               ENDIF
            ENDIF
         ELSE IF(ISECT.EQ.1) THEN
C     CAS DES SECTIONS TRAPEZOIDALES
            DO 812 I=1,3
               VY(I)=0.D0
               DO 8121 J=1,3
                  VY(I)=VY(I)+ROTA(I,J)*(PN1(J)-P1(J))
 8121          CONTINUE
 812        CONTINUE
            DY=F1
            DZ=F2
            REXT=0.5D0*DY
            RINT=-REXT
            AIRE=DY*DZ
            COUR=DENS*AIRE
            DELTA=0.5D0*DY*(PENT2-PENT1)
            HINT=DZ-DELTA
            HEXT=DZ+DELTA
            YMIN=-VY(1)
            YMAX=YMIN+GY
            XP=-VY(2)
            ZP=VY(3)
            IF(INDU) THEN
               CALL SEGTRA(COUR,HINT,HEXT,RINT,REXT,PENT1,PENT2,
     *              YMIN,YMAX,XP,ZP,BXP,BYP,BZP)
            ENDIF
            IF(POTE) THEN
               CALL ASEGTR(COUR,HINT,HEXT,RINT,REXT,PENT1,PENT2,
     *              YMIN,YMAX,XP,ZP,BXP,BYP,BZP)
            ENDIF
            RESP(1)=BYP
            RESP(2)=-BXP
            RESP(3)= BZP
         ENDIF
         GO TO 700
C
C     CAS DU FIL
 600     CONTINUE
         GY=0.D0
         DO  601 I=1,3
            VY(I)=P2(I)-P1(I)
            VI(I)=PN1(I)-P1(I)
            GY=GY+(VY(I)*VY(I))
 601     CONTINUE
         GY=SQRT(GY)
         CALL PROVEC(VI,VY,VZ)
         CALL PROVEC(VY,VZ,VX)
         CALL ROTREP (ROTA,VX,VY,VZ)
         IF(IERR.NE.0) THEN
            SEGDES MELEME
            SEGSUP MTRAV
            RETURN
         ENDIF
         XP=0.D0
         YMIN=0.D0
         DO 602 J= 1,3
            XP=XP+ROTA(1,J)*VI(J)
            YMIN=YMIN-ROTA(2,J)*VI(J)
 602     CONTINUE
         YMAX=YMIN+GY
         COUR=DENS
         IF(INDU) THEN
            CALL BFILDR(COUR,YMIN,YMAX,XP,BZP)
            RESP(1)=0.D0
            RESP(2)=0.D0
            RESP(3)=BZP
         ENDIF
         IF(POTE) THEN
            CALL AFILDR(COUR,YMIN,YMAX,XP,AYP)
            RESP(1)=0.D0
            RESP(2)=AYP
            RESP(3)=0.D0
         ENDIF
CC
 700     CONTINUE
*
C     RETOUR  DU CHAMP DANS LE REPERE  GENERAL
        
         DO 810 I= 1,3
            VY(I)=0.D0
            DO 8101 J= 1,3
               VY(I)= VY(I) + ROTA(J,I) * RESP(J) * RAP
 8101       CONTINUE
 810     CONTINUE
*      CALL RET_REF   (RESP)
*   REMPLISSAGE DU SEGMENT DE TRAVAIL MTRAV
*
         IGEO(IPT)=L
         DO 811 J= 1,3
            BB(J,IPT) = VY(J)
            IBIN(J,IPT) = 1
 811     CONTINUE
C          FIN BOUCLE SUR LES POINTS
*
 30   CONTINUE
      SEGSUP SANGLE
      SEGDES,MELEME
*
      CALL CRECHP(MTRAV,MCHPOI)
      SEGSUP MTRAV
*
*       ATTRIBUTION D'UNE NATURE DISCRETE AU CHAMP QUI SORT
*
      
 9990 SEGACT MCHPOI*MOD
      JATTRI(1) = 2
      SEGDES MCHPOI
*
*
      CALL ECROBJ('CHPOINT ',MCHPOI)
      RETURN
      END