!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
!     Software for computation of 23 quality factors
!     for a given tetrahedron.
!
!     The 2 following subroutines can be used:
!
!          cgeompara computes all the geometric parameters of the
!                    tetrahedron (distances, angles, surfaces, volume..)
!          cgeomcrit computes the 17 geometric criterions, all
!                    ranged beetween (0.-1.) and regarded as different
!                    quality factors for the measurement of vectorial
!                    parameters such as current density.
!
!     Program qualifaclist, given in a separed file, provide an example
!     of use, and compute and print the value of the Q.F. for a given
!     configuration of the tetrahedron.
!
!     Criterions 1 to 13 as described in the following papers:
!
!     Patrick ROBERT and Alain ROUX
!     INFLUENCE OF THE SHAPE OF THE TETRAHEDRON ON THE ACCURACY
!     OF THE ESTIMATE OF THE CURRENT DENSITY
!     CRPE/CNET-CNRS, 92131 Issy les Moulineaux, France
!     Proceedings of the ESA Conference on 'Spatio-Temporel Analysis
!     for Resolving plasma Turbulence' (START), 'Method for Analysing
!     Plasma Turbulence', ESA SP, Aussois, France, Jan 31-Feb 5, 1993.
!
!     Criterions 14 to 17  has been computed by J. Schoenmakers
!     ESOC/OAD, Robert Bosch Str. 5, D-64293 Darmstadt
!                       _________________________
!
!     This library has been updated with subroutine allowing computation
!     of reciprocal vector of the barycentric coordinates, delivered
!     by G. Chanteur, CETP, september 1994. These vectors are used to
!     compute easily curl(B) and div(B), and more recently (2004) the
!     radius of curvatur and the normal vector to the osculating plane.
!     *** routines not tested!! problems ****
!
!
!                                        Patrick ROBERT
!                                        CETP/CNRS-UVSQ
!                                        10-12 avenue de l'Europe
!                                        F-78140 Velizy
!
!                                        May 1994
!                         addendum september 1994 (csomdebdl)
!                   Revision for E et P, May 2001
!                            addendum March  2004 (ccurvbary)
!                            addendum March  2020 (csomdebds)
!                       _________________________
!
! For computation of Curl and div, prefer the following routines
! (test succesfull):
!         call cnormales
!         call csurfaces
!         call csomdebdl(satpos,satmag)
!         call ccurlinco(1,gjx1,gjy1,gjz1)
!         call ccurlinco(2,gjx2,gjy2,gjz2)
!         call ccurlinco(3,gjx3,gjy3,gjz3)
!         call ccurlinco(4,gjx4,gjy4,gjz4)

!         gjx=0.25*(gjx1 +gjx2 +gjx3 +gjx4)
!         gjy=0.25*(gjy1 +gjy2 +gjy3 +gjy4)
!         gjz=0.25*(gjz1 +gjz2 +gjz3 +gjz4)

!         call csomdebds(satmag)
!         call cdivlinco(1,div1)
!         call cdivlinco(2,div2)
!         call cdivlinco(3,div3)
!         call cdivlinco(4,div4)
!         div=(div1+div2+div3+div4)/4.
!
!   P. Robert, March 2020
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine cocurldiv(satpos,satmag,vjx,vjy,vjz,curl,div,ajb,elong,plana)
      
      real  satmag(3,4), satpos(3,4)
      
      real*8 dirax(3,3)
      real*8 a,b,c,elong8,plana8
      
!     ---------------------------------------------------------------+--
!     Computation of curl(B) and div(B)
!     P. Robert, revised March 2020
!
!     input:  satpos en km, satmag en nT
!     output: J en nA/m2 et div en nT/m, ajb en degres
!     ---------------------------------------------------------------+--

!     ------------------------------------------------------------------   
      common /vdensi/  gjx134 , gjx142 , gjx123 , gjx432 , &
                       gjy134 , gjy142 , gjy123 , gjy432 , &
                       gjz134 , gjz142 , gjz123 , gjz432 , &
                       gjm134 , gjm142 , gjm123 , gjm432
!     ---------------------------------------------------------------+--  !
!     conversion de div ou curl de nT/Re en nT/m
!     conv_div=1.5693E-7 ! 1./6372.2e3
  
!     conversion de curl en nT/m en J nA/m2  
!     conv_curl=0.124857 ! 1.56e-7 nT/m  /4pi e-7 Tm/A
      
      
! *   passage dans le centre de masse
!
      call  cbaryclus(satpos,pgx,pgy,pgz)
      call  transclus(satpos,pgx,pgy,pgz)
!
! *   calcul des courants etc.
!
      call  cdistance(satpos)
      call  csurfaces
      call  cnormales
      call  cvolumeto(satpos)

      call  csomdebdl(satpos,satmag)
      call  ccurlinco(1,vjx134,vjy134,vjz134)
      call  ccurlinco(2,vjx142,vjy142,vjz142)
      call  ccurlinco(3,vjx123,vjy123,vjz123)
      call  ccurlinco(4,vjx432,vjy432,vjz432)

      vjx=0.25*(vjx134 +vjx142  +vjx123 +vjx432)*1.e3  ! pour passer de A/km2 a nA/m2
      vjy=0.25*(vjy134 +vjy142  +vjy123 +vjy432)*1.e3
      vjz=0.25*(vjz134 +vjz142  +vjz123 +vjz432)*1.e3

      vjmod=sqrt(vjx**2 +vjy**2 + vjz**2)

!     rot(B)=mu0 J
      curl=4.*3.14159*1.e-7*vjmod
!
! *   calcul du champ moyen
!
      bmx=  (satmag(1,1) +satmag(1,2) +satmag(1,3) +satmag(1,4))/4.
      bmy=  (satmag(2,1) +satmag(2,2) +satmag(2,3) +satmag(2,4))/4.
      bmz=  (satmag(3,1) +satmag(3,2) +satmag(3,3) +satmag(3,4))/4.

      bmod  =sqrt(bmx**2 +bmy**2 + bmz**2)

! *   calcul de la normale N= J X Bmoy
!
      rnx=  vjy*bmz - vjz*bmy
      rny=  vjz*bmx - vjx*bmz
      rnz=  vjx*bmy - vjy*bmx
!
! *   calcul de l'angle entre J et Bmoy
!
      if(bmod < 1.e-30  .or. vjmod < 1.e-30) then
        ajb=0.
        ratio=0.
                        else
        call cangratbis(vjx,vjy,vjz,bmx,bmy,bmz,ajb,ratio)
        ajb= ajb*180./3.1415927
      endif

! calcul de la divergence

      call  csomdebds(satmag)
      call  cdivlinco(1,div1)
      call  cdivlinco(2,div2)
      call  cdivlinco(3,div3)
      call  cdivlinco(4,div4)   

      div=(div1+div2+div3+div4)/4.e3 ! pour passer de nt/km a nT/m

! *   calcul elongation et planarity

      call  cellicrit(satpos,a,b,c,dirax,elong8,plana8)
      elong=sngl(elong8)
      plana=sngl(plana8)
      
! retour aux coordonnees initiales

      call  transclus(satpos,-pgx,-pgy,-pgz)
      
      return
      end
      
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine cangratbis(ux,uy,uz,vx,vy,vz,angle,ratio) 
!                                                                       
! ----------------------------------------------------------------------
! *   Class  : basic compute modules of Rocotlib Software               
! *   Object : compute_angle_and_ratio beetween U and V vectors         
! *   Author : P. Robert, CRPE, 1992                                    
!                                                                       
! *   Input  : ux,uy,uz                                                 
!              vx,vy,vz                                                 
!                                                                       
! *   Output : sp=U.V                                                   
!              angle=angle beetween U and V (radians)                   
!              ratio= mod(U)/mod(V)                                     
! ----------------------------------------------------------------------
!                                                                       
!                                                                       
      double precision u1,u2,u3,v1,v2,v3,dp,ru,rv,cot 
!                                                                       
      u1=dble(ux) 
      u2=dble(uy) 
      u3=dble(uz) 
      v1=dble(vx) 
      v2=dble(vy) 
      v3=dble(vz) 
!                                                                       
      dp= u1*v1 + u2*v2 + u3*v3 
      ru= dsqrt(u1*u1 + u2*u2 + u3*u3) 
      rv= dsqrt(v1*v1 + v2*v2 + v3*v3) 
      cot=dp/(ru*rv) 
      cot=cot -sign(1.d-13,cot) 
      ratio=real(ru/rv) 
                                                                        
      if(cot.ge.1.d0)  then 
                 print*, ' *** Rocotlib/cp_angle_and_ratio: cos > 1 !!!' 
                 print*, '                       angle set to 0.' 
                 angle=0. 
                 return 
                 endif 
                                                                        
      if(cot.lt.-1.d0) then 
                 print*, ' *** Rocotlib/cp_angle_and_ratio: cos < 1 !!!' 
                 print*, '                       angle set to 0.' 
                 angle=0. 
                 return 
                 endif 
!                                                                       
      angle=real(dacos(cot)) 
!                                                                       
      return 
      END                                           
!   
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine cgeompara(s)
!
!----------------------------------------------------------------------!
! *   Object: compute geometric parameters of a given tetrahedron
! *   Author: P. Robert, CRPE, 1994
!----------------------------------------------------------------------!

      real s(3,4)
!
      call cdistance(s)
      call csurfaces
      call cnormales
      call canglefac
      call cvolumeto(s)
      call cspherins(s)
!
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cgeomcrit(s,valcrit,titcrit,n)
!
!----------------------------------------------------------------------!
! *   Object: compute geometric criterions and give their title
! *   Author: P. Robert, CRPE, 1994
!----------------------------------------------------------------------!
!
!
      real s(3,4)
      real*8 valcrit(n),dirabc(3,3)
      character*(*) titcrit(n)
      real*8 a,b,c,gl,gf
!
!     ------------------------------------------------------------------
!
      common /distan/     d12 ,    d13 ,    d14 ,    d23 ,  d24,  d34, &
                         dx12 ,   dx13 ,   dx14 ,   dx21 , dx31, dx41, &
                         dy12 ,   dy13 ,   dy14 ,   dy21 , dy31, dy41, &
                         dz12 ,   dz13 ,   dz14 ,   dz21 , dz31, dz41
      common /surfac/    s134 ,   s142 ,   s123 ,   s432
      common /angles/      a1 ,     a2 ,     a3 ,     a4 ,  a5 ,  a6
      common /volume/    v134 ,   v142 ,   v123 ,   v432
      common /sphere/   rsphe ,  vsphe
!
!     ------------------------------------------------------------------
!
!
      if(n.lt.25) then
      stop 'CGEOMCRIT      *** ABORTED ! CRIT. ARRAY MUST BE DIM 25'
      endif
!
      r2=sqrt(2.)
      r3=sqrt(3.)
!
      dmin=min(d12,d13,d14,d23,d24,d34)
      dmax=max(d12,d13,d14,d23,d24,d34)
      dmoy=(d12 + d13 + d14 + d23 + d24 + d34)/6.
!
      amin=min(a1,a2,a3,a4,a5,a6)
      amax=max(a1,a2,a3,a4,a5,a6)
      amoy=(a1 + a2 + a3 + a4 + a5 + a6)/6.
!
      smin=min(s134,s142,s123,s432)
      smax=max(s134,s142,s123,s432)
      smoy=r3*(dmoy**2)
!
      vmoy=(dmoy**3)*r2/12.
!
      surf=s134 + s142 + s123 + s432
      volu=(v134 + v142 + v123 + v432)/4.
!
      call inipobary(s)
      call calbabary(qfbary)
!
      rr9=r2*sqrt(r3)/36.
      rr10=2.*r3/(9.*3.1415927)
!
      valcrit( 1)=dble(dmin/dmax)
      valcrit( 2)=dble(amin/amax)
      valcrit( 3)=dble(smin/smax)
      valcrit( 4)=dble(dmin/dmoy)
      valcrit( 5)=dble(amin/amoy)
      valcrit( 6)=dble(smin*4./smoy)
      valcrit( 7)=dble(surf/smoy)
      valcrit( 8)=dble(volu/vmoy)
      valcrit( 9)=dble((volu/((sqrt(surf))**3))/rr9)
      valcrit(10)=dble(((volu/vsphe)/rr10)**(1./3.))
      valcrit(11)=dble(volu/vmoy + surf/smoy + 1.)/3.d0
      valcrit(12)=dble(volu/vmoy + surf/smoy + dmin/dmax)/3.d0
      valcrit(13)=((valcrit(11)*3.d0)-1.d0)/2.d0
!
! *** computation of crt # 14,15,16,17 and 18
!
      call cellicrit(s,a,b,c,dirabc,gl,gf)
!
      valcrit(14)=b/a
      valcrit(15)=c/a
      valcrit(16)=c/b
      valcrit(17)=(a+b+c)/(3.d0*a)
      valcrit(18)=(valcrit(17)-1.d0/3.d0)*3.d0/2.d0
!
      valcrit(19)=dble(qfbary)
      valcrit(20)=1.d0 -sqrt((gl**2 + gf**2)/2.d0)
      valcrit(21)=1.d0 -((gl+gf)**2)/4.d0
      valcrit(22)=sqrt((1.d0-gl)*(1.d0-gf))
      valcrit(23)=1.d0 -sqrt(((gl**2 + gf**2)**1.5d0)/(gl+gf))
      valcrit(24)= gl
      valcrit(25)= gf
!
!
      titcrit( 1) ='Dmin/Dmax '
      titcrit( 2) ='Angle min/Angle max '
      titcrit( 3) ='Surf. min/Surf. max '
      titcrit( 4) ='Dmin/Dmoy '
      titcrit( 5) ='Angle min/Angle moy '
      titcrit( 6) ='4*Surf. min/Surf. moy '
      titcrit( 7) ='Surf. tot/Surf. moy '
      titcrit( 8) ='Volume/Vol.  moy '
      titcrit( 9) ='(Volume/Surf. tot. **3/2)*norm '
      titcrit(10)='(Volume/vol. sphere)*norm '
      titcrit(11)='(Vol./Vol.moy + Surf./Surf.moy + 1)/3 '
      titcrit(12)='(Vol./Vol.moy + Surf./Surf.moy + Dmin/Dmax)/3 '
      titcrit(13)='(Qgm-1)/2.'
      titcrit(14)='b/a'
      titcrit(15)='c/a'
      titcrit(16)='c/b'
      titcrit(17)='(a+b+c)/3a'
      titcrit(18)='0.5(a+b+c/a -1)'
      titcrit(19)='K2min/K2max base bary'
      titcrit(20)='1.-sqrt((L2+F2)/2.)'
      titcrit(21)='1.-0.25(L+F)**2'
      titcrit(22)='sqrt[(1-L)(1-F)]'
      titcrit(23)='1-sqrt[(L2+F2)**3/2]/(L+F)'
      titcrit(24)='Elongation'
      titcrit(25)='Planarity'
!
!
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cdistance(s)
!
!----------------------------------------------------------------------!
! *   Object:  compute all the distances of the tetrahedron
! *   Author: P. Robert, CRPE, 1994
!----------------------------------------------------------------------!
!
!
      real s(3,4)
!
!     ------------------------------------------------------------------
      common /distan/     d12 ,    d13 ,    d14 ,    d23 ,  d24,  d34, &
                         dx12 ,   dx13 ,   dx14 ,   dx21 , dx31, dx41, &
                         dy12 ,   dy13 ,   dy14 ,   dy21 , dy31, dy41, &
                         dz12 ,   dz13 ,   dz14 ,   dz21 , dz31, dz41
!     ------------------------------------------------------------------
!
      dx12= s(1,2)-s(1,1)
      dy12= s(2,2)-s(2,1)
      dz12= s(3,2)-s(3,1)
!
      dx13= s(1,3)-s(1,1)
      dy13= s(2,3)-s(2,1)
      dz13= s(3,3)-s(3,1)
!
      dx14= s(1,4)-s(1,1)
      dy14= s(2,4)-s(2,1)
      dz14= s(3,4)-s(3,1)
!
      dx23= s(1,3)-s(1,2)
      dy23= s(2,3)-s(2,2)
      dz23= s(3,3)-s(3,2)
!
      dx24= s(1,4)-s(1,2)
      dy24= s(2,4)-s(2,2)
      dz24= s(3,4)-s(3,2)
!
      dx34= s(1,4)-s(1,3)
      dy34= s(2,4)-s(2,3)
      dz34= s(3,4)-s(3,3)
!
!
      dx21= s(1,1)-s(1,2)
      dy21= s(2,1)-s(2,2)
      dz21= s(3,1)-s(3,2)
!
      dx31= s(1,1)-s(1,3)
      dy31= s(2,1)-s(2,3)
      dz31= s(3,1)-s(3,3)
!
      dx41= s(1,1)-s(1,4)
      dy41= s(2,1)-s(2,4)
      dz41= s(3,1)-s(3,4)
!
!
      d12= sqrt(dx12**2 + dy12**2 + dz12**2)
      d13= sqrt(dx13**2 + dy13**2 + dz13**2)
      d14= sqrt(dx14**2 + dy14**2 + dz14**2)
      d23= sqrt(dx23**2 + dy23**2 + dz23**2)
      d24= sqrt(dx24**2 + dy24**2 + dz24**2)
      d34= sqrt(dx34**2 + dy34**2 + dz34**2)
!
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine csurfaces
!
!----------------------------------------------------------------------!
! *   Object: compute all the surfaces of the tetrahedron
! *   Author: P. Robert, CRPE, 1994
!----------------------------------------------------------------------!
!
!
!     ------------------------------------------------------------------
      common /distan/     d12 ,    d13 ,    d14 ,    d23 ,  d24,  d34, &
                         dx12 ,   dx13 ,   dx14 ,   dx21 , dx31, dx41, &
                         dy12 ,   dy13 ,   dy14 ,   dy21 , dy31, dy41, &
                         dz12 ,   dz13 ,   dz14 ,   dz21 , dz31, dz41
      common /surfac/    s134 ,   s142 ,   s123 ,   s432
!     ------------------------------------------------------------------
!
!     calcul des surfaces
!     -------------------
!
      p134= (d13 + d34 + d14)/2.
      p142= (d14 + d24 + d12)/2.
      p123= (d12 + d23 + d13)/2.
      p432= (d34 + d23 + d24)/2.
!
      s134c= p134*(p134-d13)*(p134-d34)*(p134-d14)
      s142c= p142*(p142-d14)*(p142-d24)*(p142-d12)
      s123c= p123*(p123-d12)*(p123-d23)*(p123-d13)
      s432c= p432*(p432-d34)*(p432-d23)*(p432-d24)
      
      if(s134c < 0.) then
                     s134=0.
                     else
                     s134=sqrt(s134c)
      endif
      
      
      if(s142c < 0.) then
                     s142=0.
                     else
                     s142=sqrt(s142c)
      endif      
      
      if(s123c < 0.) then
                     s123=0.
                     else
                     s123=sqrt(s123c)
      endif
      
      if(s432c < 0.) then
                     s432=0.
                     else
                     s432=sqrt(s432c)
      endif
      
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cnormales
!
!----------------------------------------------------------------------!
! *   Object: compute all the normales of the tetrahedron
! *   Author: P. Robert, CRPE, 1994
!----------------------------------------------------------------------!
!
!     ------------------------------------------------------------------
      common /distan/     d12 ,    d13 ,    d14 ,    d23 ,  d24,  d34, &
                         dx12 ,   dx13 ,   dx14 ,   dx21 , dx31, dx41, &
                         dy12 ,   dy13 ,   dy14 ,   dy21 , dy31, dy41, &
                         dz12 ,   dz13 ,   dz14 ,   dz21 , dz31, dz41
      common /normal/   rx134 ,  rx142 ,  rx123 ,  rx432 , &
                        ry134 ,  ry142 ,  ry123 ,  ry432 , &
                        rz134 ,  rz142 ,  rz123 ,  rz432 , &
                        rn134 ,  rn142 ,  rn123 ,  rn432
!     ------------------------------------------------------------------
!
!
!     calcul des normales aux faces
!     -----------------------------
!
      rx134= dy13*dz14 - dz13*dy14
      ry134= dz13*dx14 - dx13*dz14
      rz134= dx13*dy14 - dy13*dx14
!
      rx142= dy14*dz12 - dz14*dy12
      ry142= dz14*dx12 - dx14*dz12
      rz142= dx14*dy12 - dy14*dx12
!
      rx123= dy12*dz13 - dz12*dy13
      ry123= dz12*dx13 - dx12*dz13
      rz123= dx12*dy13 - dy12*dx13
!
      rx432= (dy14-dy12)*(dz13-dz12) - (dz14-dz12)*(dy13-dy12)
      ry432= (dz14-dz12)*(dx13-dx12) - (dx14-dx12)*(dz13-dz12)
      rz432= (dx14-dx12)*(dy13-dy12) - (dy14-dy12)*(dx13-dx12)
!
! *** normalisation des normales
!
      rn134=sqrt(rx134**2 + ry134**2 + rz134**2)
      rn142=sqrt(rx142**2 + ry142**2 + rz142**2)
      rn123=sqrt(rx123**2 + ry123**2 + rz123**2)
      rn432=sqrt(rx432**2 + ry432**2 + rz432**2)

      if(rn134 .lt. 1.e-30) then
         rx134=0.
         ry134=0.
         rz134=0.
         else
         rx134=rx134/rn134
         ry134=ry134/rn134
         rz134=rz134/rn134
      endif

      if(rn142 .lt. 1.e-30) then
         rx142=0.
         ry142=0.
         rz142=0.
         else
         rx142=rx142/rn142
         ry142=ry142/rn142
         rz142=rz142/rn142
      endif

      if(rn123 .lt. 1.e-30) then
         rx123=0.
         ry123=0.
         rz123=0.
         else
         rx123=rx123/rn123
         ry123=ry123/rn123
         rz123=rz123/rn123
      endif

      if(rn432 .lt. 1.e-30) then
         rx432=0.
         ry432=0.
         rz432=0.
         else
         rx432=rx432/rn432
         ry432=ry432/rn432
         rz432=rz432/rn432
      endif
!
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine canglefac
!
!----------------------------------------------------------------------!
! *   Object: compute all the angles of the tetrahedron
! *   Author: P. Robert, CRPE, 1994
!----------------------------------------------------------------------!
!
!
!     ------------------------------------------------------------------
      common /normal/   rx134 ,  rx142 ,  rx123 ,  rx432 , &
                        ry134 ,  ry142 ,  ry123 ,  ry432 , &
                        rz134 ,  rz142 ,  rz123 ,  rz432 , &
                        rn134 ,  rn142 ,  rn123 ,  rn432
      common /angles/      a1 ,     a2 ,     a3 ,     a4 , a5 , a6
!     ------------------------------------------------------------------
!
!     angles entre faces
!     ------------------
!
      a1= rx134*rx142 + ry134*ry142 + rz134*rz142
      a2= rx134*rx123 + ry134*ry123 + rz134*rz123
      a3= rx134*rx432 + ry134*ry432 + rz134*rz432
      a4= rx142*rx123 + ry142*ry123 + rz142*rz123
      a5= rx142*rx432 + ry142*ry432 + rz142*rz432
      a6= rx123*rx432 + ry123*ry432 + rz123*rz432

      if(abs(a1).gt.1.) a1= sign(1.,a1)
      if(abs(a2).gt.1.) a2= sign(1.,a2)
      if(abs(a3).gt.1.) a3= sign(1.,a3)
      if(abs(a4).gt.1.) a4= sign(1.,a4)
      if(abs(a5).gt.1.) a5= sign(1.,a5)
      if(abs(a6).gt.1.) a6= sign(1.,a6)

!      print 100, a1,a2,a3,a4,a5,a6
!  100 format(6F12.7)

      a1= acos(a1)
      a2= acos(a2)
      a3= acos(a3)
      a4= acos(a4)
      a5= acos(a5)
      a6= acos(a6)

      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cvolumeto(s)
!
!----------------------------------------------------------------------!
! *   Object: compute the volume of the tetrahedron
! *   Author: P. Robert, CRPE, 1994
!----------------------------------------------------------------------!
!
!
      real s(3,4)
!
!     ------------------------------------------------------------------
      common /normal/   rx134 ,  rx142 ,  rx123 ,  rx432 , &
                        ry134 ,  ry142 ,  ry123 ,  ry432 , &
                        rz134 ,  rz142 ,  rz123 ,  rz432 , &
                        rn134 ,  rn142 ,  rn123 ,  rn432
      common /surfac/    s134 ,   s142 ,   s123 ,   s432
      common /volume/    v134 ,   v142 ,   v123 ,   v432
!     ------------------------------------------------------------------
!
!     volume du tetraedre
!     -------------------

      tol=0.01
!
! *** face de base 134
!
      p134 =s(1,1)*rx134 + s(2,1)*ry134 + s(3,1)*rz134
      p134b=s(1,3)*rx134 + s(2,3)*ry134 + s(3,3)*rz134
      p134c=s(1,4)*rx134 + s(2,4)*ry134 + s(3,4)*rz134
!
      h134=abs(s(1,2)*rx134 + s(2,2)*ry134 + s(3,2)*rz134 -p134)
      v134=s134*h134/3.
!
      if(h134 .gt. 1.e-10) then
         erb=abs((p134b-p134)/h134)
         erc=abs((p134c-p134)/h134)
         if(erb.gt.tol) print 100, &
                         'lib_CLU_tools: unprecise p134-b :',erb
         if(erc.gt.tol) print 100, &
                         'lib_CLU_tools: unprecise p134-c :',erc
      endif
!
! *** face de base 142
!
      p142 =s(1,1)*rx142 + s(2,1)*ry142 + s(3,1)*rz142
      p142b=s(1,4)*rx142 + s(2,4)*ry142 + s(3,4)*rz142
      p142c=s(1,2)*rx142 + s(2,2)*ry142 + s(3,2)*rz142
!
      h142=abs(s(1,3)*rx142 + s(2,3)*ry142 + s(3,3)*rz142 -p142)
      v142=s142*h142/3.
!
      if(h142 .gt. 1.e-10) then
         erb=abs((p142b-p142)/h142)
         erc=abs((p142c-p142)/h142)
         if(erb.gt.tol) print 100, &
                         'lib_CLU_tools: unprecise p142-b :',erb
         if(erc.gt.tol) print 100, &
                         'lib_CLU_tools: unprecise p142-c :',erc
      endif
!
! *** face de base 123
!
      p123 =s(1,1)*rx123 + s(2,1)*ry123 + s(3,1)*rz123
      p123b=s(1,2)*rx123 + s(2,2)*ry123 + s(3,2)*rz123
      p123c=s(1,3)*rx123 + s(2,3)*ry123 + s(3,3)*rz123
!
      h123=abs(s(1,4)*rx123 + s(2,4)*ry123 + s(3,4)*rz123 -p123)
      v123=s123*h123/3.
!
      if(h123 .gt. 1.e-10) then
         erb=abs((p123b-p123)/h123)
         erc=abs((p123c-p123)/h123)
         if(erb.gt.tol) print 100, &
                         'lib_CLU_tools: unprecise p123-b :',erb
         if(erc.gt.tol) print 100, &
                         'lib_CLU_tools: unprecise p123-c :',erc
      endif
!
! *** face de base 432
!
      p432 =s(1,4)*rx432 + s(2,4)*ry432 + s(3,4)*rz432
      p432b=s(1,3)*rx432 + s(2,3)*ry432 + s(3,3)*rz432
      p432c=s(1,2)*rx432 + s(2,2)*ry432 + s(3,2)*rz432
!
      h432=abs(s(1,1)*rx432 + s(2,1)*ry432 + s(3,1)*rz432 -p432)
      v432=s432*h432/3.
!
      if(h432 .gt. 1.e-10) then
         erb=abs((p432b-p432)/h432)
         erc=abs((p432c-p432)/h432)
         if(erb.gt.tol) print 100, &
                         'lib_CLU_tools: unprecise p432-b :',erb
         if(erc.gt.tol) print 100, &
                         'lib_CLU_tools: unprecise p432-c :',erc
      endif
!
  100 format(a,f13.6)
!
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cspherins(s)
!
!----------------------------------------------------------------------!
! *   Object: compute radius&volume of sphere including tetrahedron
! *   Author: P. Robert, CRPE, 1994
!----------------------------------------------------------------------!
!
!
      real s(3,4)
!
!     ------------------------------------------------------------------
      common /distan/     d12 ,    d13 ,    d14 ,    d23 ,  d24,  d34, &
                         dx12 ,   dx13 ,   dx14 ,   dx21 , dx31, dx41, &
                         dy12 ,   dy13 ,   dy14 ,   dy21 , dy31, dy41, &
                         dz12 ,   dz13 ,   dz14 ,   dz21 , dz31, dz41
      common /sphere/   rsphe ,  vsphe
!     ------------------------------------------------------------------
!
!     sphere inscrite
!     ---------------
!
      gx=s(1,2)*dx21 + s(2,2)*dy21 + s(3,2)*dz21
      gy=s(1,3)*dx31 + s(2,3)*dy31 + s(3,3)*dz31
      gz=s(1,4)*dx41 + s(2,4)*dy41 + s(3,4)*dz41
!
      call setmatrix(dx21,dy21,dz21, dx31,dy31,dz31, dx41,dy41,dz41)
      call detmatrix(det)
      call invmatrix(gx,gy,gz,px,py,pz)
!
      rsphe=sqrt((s(1,1)-px)**2 + (s(2,1)-py)**2 + (s(3,1)-pz)**2)/2.
      vsphe=4.*3.1415927*(rsphe**3)/3.
!
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine csomdebdl(s,bmag)
!
!----------------------------------------------------------------------!
! *   Object: computation of J vector inside the 4 S/C tetrahedron
!             by the means of the ingrations over each of the 4 faces
! *   Author: P. Robert, CRPE, 1994
!----------------------------------------------------------------------!
!
!
      real s(3,4)
      real bmag(3,4)
!
!     ------------------------------------------------------------------
      common /distan/     d12 ,    d13 ,    d14 ,    d23 ,  d24,  d34, &
                         dx12 ,   dx13 ,   dx14 ,   dx21 , dx31, dx41, &
                         dy12 ,   dy13 ,   dy14 ,   dy21 , dy31, dy41, &
                         dz12 ,   dz13 ,   dz14 ,   dz21 , dz31, dz41
      common /surfac/    s134 ,   s142 ,   s123 ,   s432
      common /normal/   rx134 ,  rx142 ,  rx123 ,  rx432 , &
                        ry134 ,  ry142 ,  ry123 ,  ry432 , &
                        rz134 ,  rz142 ,  rz123 ,  rz432 , &
                        rn134 ,  rn142 ,  rn123 ,  rn432
      common /determ/    d134 ,   d142 ,   d123 ,   d432
      common /densit/   gj134 ,  gj142 ,  gj123 ,  gj432
      
      common /vdensi/  gjx134 , gjx142 , gjx123 , gjx432 , &
                       gjy134 , gjy142 , gjy123 , gjy432 , &
                       gjz134 , gjz142 , gjz123 , gjz432 , &
                       gjm134 , gjm142 , gjm123 , gjm432
!     ------------------------------------------------------------------
!     valeur de Mu0 en nT*km/A
      rmu=4.*acos(-1.)*(1.e-7)*(1.e6)
!
!     les unites d'entree sont en nT et en km
!     en sortie: J en A/km2
!     1 A/km2 = 1000 nA/m2
!     ------------------------------------------------------------------
!
! *** calcul des vecteurs pour la face 1-3-4
!
      dx13= s(1,3)-s(1,1)
      dy13= s(2,3)-s(2,1)
      dz13= s(3,3)-s(3,1)
!
      dx34= s(1,4)-s(1,3)
      dy34= s(2,4)-s(2,3)
      dz34= s(3,4)-s(3,3)
!
      dx41= s(1,1)-s(1,4)
      dy41= s(2,1)-s(2,4)
      dz41= s(3,1)-s(3,4)
!
! *** calcul des vecteurs pour la face 1-4-2
!
      dx14= s(1,4)-s(1,1)
      dy14= s(2,4)-s(2,1)
      dz14= s(3,4)-s(3,1)
!
      dx42= s(1,2)-s(1,4)
      dy42= s(2,2)-s(2,4)
      dz42= s(3,2)-s(3,4)
!
      dx21= s(1,1)-s(1,2)
      dy21= s(2,1)-s(2,2)
      dz21= s(3,1)-s(3,2)
!
! *** calcul des vecteurs pour la face 1-2-3
!
      dx12= s(1,2)-s(1,1)
      dy12= s(2,2)-s(2,1)
      dz12= s(3,2)-s(3,1)
!
      dx23= s(1,3)-s(1,2)
      dy23= s(2,3)-s(2,2)
      dz23= s(3,3)-s(3,2)
!
      dx31= s(1,1)-s(1,3)
      dy31= s(2,1)-s(2,3)
      dz31= s(3,1)-s(3,3)
!
! *** calcul des vecteurs pour la face 4-3-2
!
      dx43= s(1,3)-s(1,4)
      dy43= s(2,3)-s(2,4)
      dz43= s(3,3)-s(3,4)
!
      dx32= s(1,2)-s(1,3)
      dy32= s(2,2)-s(2,3)
      dz32= s(3,2)-s(3,3)
!
      dx24= s(1,4)-s(1,2)
      dy24= s(2,4)-s(2,2)
      dz24= s(3,4)-s(3,2)
!
!
!     calculs des champs moyens au milieu de chaque vecteur
!     -----------------------------------------------------
!
      bx1=bmag(1,1)
      bx2=bmag(1,2)
      bx3=bmag(1,3)
      bx4=bmag(1,4)
!
      by1=bmag(2,1)
      by2=bmag(2,2)
      by3=bmag(2,3)
      by4=bmag(2,4)
!
      bz1=bmag(3,1)
      bz2=bmag(3,2)
      bz3=bmag(3,3)
      bz4=bmag(3,4)
!
!
! *** calcul des champs au milieu des vecteurs pour la face 1-3-4
!
      bx13= (bx1 + bx3) /2.
      by13= (by1 + by3) /2.
      bz13= (bz1 + bz3) /2.
!
      bx34= (bx3 + bx4) /2.
      by34= (by3 + by4) /2.
      bz34= (bz3 + bz4) /2.
!
      bx41= (bx4 + bx1) /2.
      by41= (by4 + by1) /2.
      bz41= (bz4 + bz1) /2.
!
! *** calcul des champs au milieu des vecteurs pour la face 1-4-2
!
      bx14= (bx1 + bx4) /2.
      by14= (by1 + by4) /2.
      bz14= (bz1 + bz4) /2.
!
      bx42= (bx4 + bx2) /2.
      by42= (by4 + by2) /2.
      bz42= (bz4 + bz2) /2.
!
      bx21= (bx2 + bx1) /2.
      by21= (by2 + by1) /2.
      bz21= (bz2 + bz1) /2.
!
! *** calcul des champs au milieu des vecteurs pour la face 1-2-3
!
      bx12= (bx1 + bx2) /2.
      by12= (by1 + by2) /2.
      bz12= (bz1 + bz2) /2.
!
      bx23= (bx2 + bx3) /2.
      by23= (by2 + by3) /2.
      bz23= (bz2 + bz3) /2.
!
      bx31= (bx3 + bx1) /2.
      by31= (by3 + by1) /2.
      bz31= (bz3 + bz1) /2.
!
! *** calcul des champs au milieu des vecteurs pour la face 4-3-2
!
      bx43= (bx4 + bx3) /2.
      by43= (by4 + by3) /2.
      bz43= (bz4 + bz3) /2.
!
      bx32= (bx3 + bx2) /2.
      by32= (by3 + by2) /2.
      bz32= (bz3 + bz2) /2.
!
      bx24= (bx2 + bx4) /2.
      by24= (by2 + by4) /2.
      bz24= (bz2 + bz4) /2.
!
!
!     calcul des integrales de contour
!     --------------------------------
!
!
! *** calcul de l'integrale sur le contour correspondant a la face 1-3-4
!
      sb134= bx13*dx13 + by13*dy13 + bz13*dz13 &
           + bx34*dx34 + by34*dy34 + bz34*dz34 &
           + bx41*dx41 + by41*dy41 + bz41*dz41
!
! *** calcul de l'integrale sur le contour correspondant a la face 1-4-2
!
      sb142= bx14*dx14 + by14*dy14 + bz14*dz14 &
           + bx42*dx42 + by42*dy42 + bz42*dz42 &
           + bx21*dx21 + by21*dy21 + bz21*dz21
!
! *** calcul de l'integrale sur le contour correspondant a la face 1-2-3
!
      sb123= bx12*dx12 + by12*dy12 + bz12*dz12 &
           + bx23*dx23 + by23*dy23 + bz23*dz23 &
           + bx31*dx31 + by31*dy31 + bz31*dz31
!
! *** calcul de l'integrale sur le contour correspondant a la face 4-3-2
!
      sb432= bx43*dx43 + by43*dy43 + bz43*dz43 &
           + bx32*dx32 + by32*dy32 + bz32*dz32 &
           + bx24*dx24 + by24*dy24 + bz24*dz24
!
!
!     calcul des intensites sur chaque face
!     -------------------------------------
!
!     avec Mu0 en nT*km/A
      gi134=sb134/rmu
      gi142=sb142/rmu
      gi123=sb123/rmu
      gi432=sb432/rmu
!
!
!     calcul des densites sur chaque face
!     -----------------------------------
!
      if(s134 < 1. .or. s142 < 1. .or. s123 < 1. .or. s432 < 1.) then
          print*, '*** csomdebdl: one face area < 1 km2, J set to zero'
          
          gjx134=0. 
          gjx142=0. 
          gjx123=0. 
          gjx432=0. 
 
          gjy134=0. 
          gjy142=0. 
          gjy123=0. 
          gjy432=0. 
 
          gjz134=0. 
          gjz142=0. 
          gjz123=0. 
          gjz432=0. 
 
          gjm134=0. 
          gjm142=0. 
          gjm123=0. 
          gjm432=0.
          
          return
      endif
          
      gj134=gi134/s134
      gj142=gi142/s142
      gj123=gi123/s123
      gj432=gi432/s432
!
!     ---------------------------------------------------
!     calcul des determinants et inversion des 4 matrices
!     pour passer du systeme non orthogal au systeme orthogonal
!     des mesures de B et Pos.
!     ---------------------------------------------------
!
      call setmatrix(rx134,ry134,rz134, &
                     rx142,ry142,rz142, &
                     rx123,ry123,rz123)
      call detmatrix(d432)
      call invmatrix(gj134,gj142,gj123,gjx432,gjy432,gjz432)
!
      call setmatrix(rx142,ry142,rz142, &
                     rx123,ry123,rz123, &
                     rx432,ry432,rz432)
      call detmatrix(d134)
      call invmatrix(gj142,gj123,gj432,gjx134,gjy134,gjz134)
!
      call setmatrix(rx123,ry123,rz123, &
                     rx432,ry432,rz432, &
                     rx134,ry134,rz134)
      call detmatrix(d142)
      call invmatrix(gj123,gj432,gj134,gjx142,gjy142,gjz142)
!
      call setmatrix(rx432,ry432,rz432, &
                     rx134,ry134,rz134, &
                     rx142,ry142,rz142)
      call detmatrix(d123)
      call invmatrix(gj432,gj134,gj142,gjx123,gjy123,gjz123)
      
!     precautions tetrahedre plat pour J

      valmin= -1.e15
      valmax=  1.e15

      if( gjx134 > valmax) gjx134= valmax
      if( gjx142 > valmax) gjx142= valmax
      if( gjx123 > valmax) gjx123= valmax
      if( gjx432 > valmax) gjx432= valmax
      
      if( gjy134 > valmax) gjy134= valmax
      if( gjy142 > valmax) gjy142= valmax
      if( gjy123 > valmax) gjy123= valmax
      if( gjy432 > valmax) gjy432= valmax
      
      if( gjz134 > valmax) gjz134= valmax
      if( gjz142 > valmax) gjz142= valmax
      if( gjz123 > valmax) gjz123= valmax
      if( gjz432 > valmax) gjz432= valmax
      
      if( gjx134 < valmin) gjx134= valmin
      if( gjx142 < valmin) gjx142= valmin
      if( gjx123 < valmin) gjx123= valmin
      if( gjx432 < valmin) gjx432= valmin
      
      if( gjy134 < valmin) gjy134= valmin
      if( gjy142 < valmin) gjy142= valmin
      if( gjy123 < valmin) gjy123= valmin
      if( gjy432 < valmin) gjy432= valmin
      
      if( gjz134 < valmin) gjz134= valmin
      if( gjz142 < valmin) gjz142= valmin
      if( gjz123 < valmin) gjz123= valmin
      if( gjz432 < valmin) gjz432= valmin

!
!     module des 4 J calcules
!     -----------------------
!
      gjm432=sqrt(gjx432**2 + gjy432**2 + gjz432**2)
      gjm134=sqrt(gjx134**2 + gjy134**2 + gjz134**2)
      gjm142=sqrt(gjx142**2 + gjy142**2 + gjz142**2)
      gjm123=sqrt(gjx123**2 + gjy123**2 + gjz123**2)

      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine ccurlinco(iface,denjx,denjy,denjz)
!
!----------------------------------------------------------------------!
! *   Object: compute J from a given face of the tetrahedron
! *   Author: P. Robert, CRPE, 1994
!----------------------------------------------------------------------!
!
!
!     ------------------------------------------------------------------
      common /vdensi/  gjx134 , gjx142 , gjx123 , gjx432 , &
                       gjy134 , gjy142 , gjy123 , gjy432 , &
                       gjz134 , gjz142 , gjz123 , gjz432 , &
                       gjm134 , gjm142 , gjm123 , gjm432
!     ------------------------------------------------------------------
!
!
      if(iface.eq.1) then
                     denjx=gjx134
                     denjy=gjy134
                     denjz=gjz134
                     endif
!
      if(iface.eq.2) then
                     denjx=gjx142
                     denjy=gjy142
                     denjz=gjz142
                     endif
!
      if(iface.eq.3) then
                     denjx=gjx123
                     denjy=gjy123
                     denjz=gjz123
                     endif
!
      if(iface.eq.4) then
                     denjx=gjx432
                     denjy=gjy432
                     denjz=gjz432
                     endif
!
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine csomdebds(bmag)
!
!----------------------------------------------------------------------!
! *   Object: computation of dib(B)
!             by the means of the ingrations Som B.DS
! *   Author: P. Robert, ScientiDev, 2020
!----------------------------------------------------------------------!
!
      real bmag(3,4)
!
!     ------------------------------------------------------------------
      common /surfac/    s134 ,   s142 ,   s123 ,   s432
      common /normal/   rx134 ,  rx142 ,  rx123 ,  rx432 , &
                        ry134 ,  ry142 ,  ry123 ,  ry432 , &
                        rz134 ,  rz142 ,  rz123 ,  rz432 , &
                        rn134 ,  rn142 ,  rn123 ,  rn432
      common /volume/    v134 ,   v142 ,   v123 ,   v432
      common /determ/    d134 ,   d142 ,   d123 ,   d432
      common /divdeb/  div134 , div142 , div123 , div432
!     ------------------------------------------------------------------
!

!     champ moyen pour la face 134

      bx134= (bmag(1,1) +bmag(1,3) +bmag(1,4))/3.
      by134= (bmag(2,1) +bmag(2,3) +bmag(2,4))/3.
      bz134= (bmag(3,1) +bmag(3,3) +bmag(3,4))/3.

!     champ moyen pour la face 142

      bx142= (bmag(1,1) +bmag(1,4) +bmag(1,2))/3.
      by142= (bmag(2,1) +bmag(2,4) +bmag(2,2))/3.
      bz142= (bmag(3,1) +bmag(3,4) +bmag(3,2))/3.


!     champ moyen pour la face 123

      bx123= (bmag(1,1) +bmag(1,2) +bmag(1,3))/3.
      by123= (bmag(2,1) +bmag(2,2) +bmag(2,3))/3.
      bz123= (bmag(3,1) +bmag(3,2) +bmag(3,3))/3.


!     champ moyen pour la face 432

      bx432= (bmag(1,4) +bmag(1,3) +bmag(1,2))/3.
      by432= (bmag(2,4) +bmag(2,3) +bmag(2,2))/3.
      bz432= (bmag(3,4) +bmag(3,3) +bmag(3,2))/3.

!     calcul de som B.ds

      sb134= (bx134*rx134 +by134*ry134 +bz134*rz134)*s134
      sb142= (bx142*rx142 +by142*ry142 +bz142*rz142)*s142
      sb123= (bx123*rx123 +by123*ry123 +bz123*rz123)*s123
      sb432= (bx432*rx432 +by432*ry432 +bz432*rz432)*s432

!     calcul des divergences sur chaque face

      vol=0.25*(v134 +v142 +v123 +v432)

      div134=sb134/vol
      div142=sb142/vol
      div123=sb123/vol
      div432=sb432/vol

      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cdivlinco(iface,divB)
!
!----------------------------------------------------------------------!
! *   Object: compute div(B) from a given face of the tetrahedron
! *   Author: P. Robert, Scientidev 2020
!----------------------------------------------------------------------!
!
!
!     ------------------------------------------------------------------
      common /divdeb/  div134 , div142 , div123 , div432
!     ------------------------------------------------------------------
!
!
      if(iface.eq.1) divB=div134
      if(iface.eq.2) divB=div142
      if(iface.eq.3) divB=div123
      if(iface.eq.4) divB=div432
!
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine setmatrix(tax,tay,taz, tbx,tby,tbz, tcx,tcy,tcz)
!
!----------------------------------------------------------------------!
! *   Object: set matrices for further linear resolution system
! *   Author: P. Robert, CRPE, 1994
!----------------------------------------------------------------------!
!
!
!     ------------------------------------------------------------------
      common/imat1/ ax,ay,az, bx,by,bz, cx,cy,cz, dd
!     ------------------------------------------------------------------
!
!     le systeme a inverser est de la forme:
!
!     gx= ax*px + ay*py + az*cz
!     gy= bx*px + by*py + bz*pz
!     gz= cx*px + cy*py + cz*pz
!
!     connaissant gx,gy,gz l'inversion permet de calculer px,py,pz
!     les termes ai,bi,ci etant connus et constants
!
      ax=tax
      ay=tay
      az=taz
!
      bx=tbx
      by=tby
      bz=tbz
!
      cx=tcx
      cy=tcy
      cz=tcz
!
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine detmatrix(det)
!
!----------------------------------------------------------------------!
! *   Object: compute determinant of the matrix defined by setmat
! *   Author: P. Robert, CRPE, 1994
!----------------------------------------------------------------------!
!
!
!     ------------------------------------------------------------------
      common/imat1/ ax,ay,az, bx,by,bz, cx,cy,cz, dd
!     ------------------------------------------------------------------
!
      det= ax*by*cz + bx*cy*az + cx*ay*bz &
         - cx*by*az - ax*cy*bz - bx*ay*cz
!
      if (abs(det) .gt. 1.e-30) then
                                dd=det
                                else
                                dd= 1.e-30
                                endif
!
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine invmatrix(gx,gy,gz,px,py,pz)
!
!----------------------------------------------------------------------!
! *   Object: compute inverse matrix determined by setmat and detmat
! *   Author: P. Robert, CRPE, 1994
!----------------------------------------------------------------------!
!
!
!     ------------------------------------------------------------------
      common/imat1/ ax,ay,az, bx,by,bz, cx,cy,cz, dd
!     ------------------------------------------------------------------
!
      px= gx*by*cz + gy*cy*az + gz*ay*bz &
        - gz*by*az - gx*cy*bz - gy*ay*cz
!
      py= ax*gy*cz + bx*gz*az + cx*gx*bz &
        - cx*gy*az - ax*gz*bz - bx*gx*cz
!
      pz= ax*by*gz + bx*cy*gx + cx*ay*gy &
        - cx*by*gx - ax*cy*gy - bx*ay*gz
!
      px=px/dd
      py=py/dd
      pz=pz/dd
!
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cellicrit(s,a,b,c,dirax,gl,gf)
!
!----------------------------------------------------------------------!
! *   Object: Compute criterions defined by J. Schoenmaker
! *   Author: P. Robert, CRPE, 1994
!----------------------------------------------------------------------!
!
      real s(3,4)
!
      real*8 pos(3,4)
      real*8 sax(3),dir(3,3)
      real*8 a,b,c,gl,gf
      real*8 dirax(3,3)
!
!
      do 10 i=1,3
      do 10 j=1,4
      pos(i,j)=dble(s(i,j))
   10 continue
!
      call quaellips(pos,sax,dir)
!
      a=sax(1)
      b=sax(2)
      c=sax(3)
!
      gl=1.d0 -b/a
      if(b.lt.1e-30) then
                     gf=1.
                     else
                     gf=1.d0 -c/b
      endif
!
      do 20 i=1,3
      do 20 j=1,3
      dirax(i,j)=dir(i,j)
   20 continue
!
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
     subroutine cellipara(s,scalax,vectax)
!
!----------------------------------------------------------------------!
! *   Object: computation of the inertial ellipsoide
! *   Author: P. Robert, CRPE, 1994
!----------------------------------------------------------------------!
!
!     input: position of 4 sat s(3,4)
!     output: module of a,b,c in scalax(3)
!             direction of a,b,c in vectax(3,3)
!                                with a=(x-z,1), etc.
!
      real     s(3,4) , scalax(3) , vectax(3,3)
      real*8 pos(3,4) ,    sax(3) ,  dirax(3,3)
!
!
      do 10 i=1,3
      do 10 j=1,4
      pos(i,j)=dble(s(i,j))
   10 continue
!
      call quaellips(pos,sax,dirax)
!
      do 20 i=1,3
      scalax(i)=real(sax(i))
      do 20 j=1,3
      vectax(i,j)=real(dirax(i,j))
   20 continue
!
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

     subroutine cbaryclus(s,gx,gy,gz)

!----------------------------------------------------------------------!
! *   Object: compute barycentre of 4 points of cluster
! *   Author: P. Robert, CRPE, 1994
!----------------------------------------------------------------------!

      real s(3,4)

      gx=(s(1,1)+s(1,2)+s(1,3)+s(1,4))/4.
      gy=(s(2,1)+s(2,2)+s(2,3)+s(2,4))/4.
      gz=(s(3,1)+s(3,2)+s(3,3)+s(3,4))/4.

      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine transclus(s,gx,gy,gz)

!----------------------------------------------------------------------!
! *   Object: tranlate position of the 4 points of cluster of gx,gy,gz
! *   Author: P. Robert, CRPE, 1994
!----------------------------------------------------------------------!

      real s(3,4)

      do 10 j=1,4
      s(1,j)=s(1,j)-gx
      s(2,j)=s(2,j)-gy
      s(3,j)=s(3,j)-gz
   10 continue

      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine quaellips(POS,SAX,DIR)
!
!----------------------------------------------------------------------!
! *   Object: compute the properties of the elipse for a tetrahedron
! *   Author: J. Schoenmaekers, ESOC, 1993
!----------------------------------------------------------------------!
!
!     ******************************************************************
!     Software provided by J. Schoenmakers
!     ESOC/OAD, Robert Bosch Str. 5, D-64293 Darmstadt
!     ******************************************************************
!
!C =====================================================================
!C COMPUTES THE PROPERTIES OF THE QUALITY ELLIPSE FOR A TETRAHEDRON
!C =====================================================================
!C PROJ=CLU, SUBJ=MAN, AUTH=J.SCHOENMAEKERS(ESOC), DATA=26/11/93
!C =====================================================================
!I POS(3,4) R*8 POS(1..3,J): POSITION VECTOR OF EDGE J (J=1..4)
!O SAX(3)   R*8 SAX(1): SEMI-MAJOR  AXIS
!O              SAX(2): SEMI-MIDDLE AXIS
!O              SAX(3): SEMI-MINOR  AXIS
!O DIR(3,3) R*8 DIR(1..3,1): DIRECTION COSINES OF SEMI-MAJOR  AXIS
!O              DIR(1..3,2): DIRECTION COSINES OF SEMI-MIDDLE AXIS
!O              DIR(1..3,3): DIRECTION COSINES OF SEMI-MINOR  AXIS
! ======================================================================
      IMPLICIT REAL*8 (A-H, O-Z)
      IMPLICIT INTEGER*4 (I-N)
! ======================================================================
      DIMENSION POS(3,4)
      DIMENSION SAX(3)
      DIMENSION DIR(3,3)
      DIMENSION POM(3)
      DIMENSION POC(3)
      DIMENSION SEP(3,3)
      DIMENSION SYM(6)
! ======================================================================
      POM(1) = 0.25D0 * (POS(1,1) + POS(1,2) + POS(1,3) + POS(1,4))
      POM(2) = 0.25D0 * (POS(2,1) + POS(2,2) + POS(2,3) + POS(2,4))
      POM(3) = 0.25D0 * (POS(3,1) + POS(3,2) + POS(3,3) + POS(3,4))
      DO 1000 J = 1, 3
      SEP(1,J) = 0.D0
      SEP(2,J) = 0.D0
      SEP(3,J) = 0.D0
 1000 CONTINUE
      DO 2000 I = 1, 4
      POC(1) = POS(1,I) - POM(1)
      POC(2) = POS(2,I) - POM(2)
      POC(3) = POS(3,I) - POM(3)
      DO 3000 J = 1, 3
      SEP(1,J) = SEP(1,J) + 2.D0 * POC(1) * POC(J)
      SEP(2,J) = SEP(2,J) + 2.D0 * POC(2) * POC(J)
      SEP(3,J) = SEP(3,J) + 2.D0 * POC(3) * POC(J)
 3000 CONTINUE
 2000 CONTINUE
! ----------------------------------------------------------------------
      SYM(1) = SEP(1,1)
      SYM(2) = SEP(1,2)
      SYM(3) = SEP(2,2)
      SYM(4) = SEP(1,3)
      SYM(5) = SEP(2,3)
      SYM(6) = SEP(3,3)
      CALL CEIGENVAL (SYM, DIR, 3, 0)
! ----------------------------------------------------------------------
      if(SYM(1) < 0.) SYM(1)=0. ! P. Robert
      if(SYM(3) < 0.) SYM(3)=0. ! P. Robert
      if(SYM(6) < 0.) SYM(6)=0. ! P. Robert
      
      SAX(1) = DSQRT (SYM(1))
      SAX(2) = DSQRT (SYM(3))
      SAX(3) = DSQRT (SYM(6))
! ----------------------------------------------------------------------
      DO 6000 J = 1, 3
      IF (DIR(1,J) .GE. 0.D0) GOTO 6000
      DIR(1,J) = - DIR(1,J)
      DIR(2,J) = - DIR(2,J)
      DIR(3,J) = - DIR(3,J)
 6000 CONTINUE
! ======================================================================
      RETURN
      END

! ======================================================================
      subroutine CEIGENVAL(A,R,N,MV)
!
!----------------------------------------------------------------------!
! *   Object: computes eigenvalues&eigenvectors of real symm. matrix
! *   Author: J. Schoenmaekers, ESOC, 1993
!----------------------------------------------------------------------!

! ======================================================================
!P COMPUTES EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC MATRIX
!I A      R8     NN REAL SYMMETRIC (N*N)-MATRIX IS STORED IN A VECTOR OF
!I                  DIMENSION NN=N*(N+1)/2. THE UPPER PART OF THE
!I                  SYMMETRIC MATRIX IS STORED COLUMNWISE AS  1  2  4  7
!I                  AND IS CONVERTED BY THE FORMULA              3  5  8
!I                  A(K*(K-1)/2+I)=AA(I,K) WITH I=1,N AND K=I,N.    6  9
!I N      I4        ORDER OF MATRICES A AND R                         10
!I MV     I4        INPUT CODE
!I                  MV=0  COMPUTE EIGENVALUES AND EIGENVECTORS
!I                  MV=1  COMPUTE EIGENVALUES ONLY (R NEED NOT TO BE
!I                        DIMENSIONED BUT MUST STILL APPEAR IN
!I                        CALLING SEQUENCE)
!O A      R8     NN EIGENVALUES ARE DEVELOPED IN DIAGONAL OF MATRIX A
!O                  IN DESCENDING ORDER
!O R      R8    N*N MATRIX OF EIGENVECTORS (STORED COLUMNWISE, IN SAME
!O                  SEQUENCE AS EIGENVALUES)
!M MATRIX A CANNOT BE IN THE SAME LOCATION AS MATRIX R
!M DIAGONALIZATION METHOD ORIGINATED BY JACOBI AND ADAPTED
!M BY VON NEUMANN FOR LARGE COMPUTERS AS FOUND IN 'MATHEMATICAL METHODS
!M FOR DIGITAL COMPUTERS', EDITED BY A. RALSTON AND H.S. WILF,
!M JOHN WILEY AND SONS, NEW YORK, 1962, CHAPTER 7 .
! ======================================================================
      DOUBLE PRECISION A,R,ANORM,ANRMX,THR,X,Y,SINX,SINX2,COSX,COSX2, &
      SINCS,RANGE
! ======================================================================
! *** modif P.Robert June 1995
!     les dimensions initiales A(1) et R(1) sont etendues pour eviter
!     les erreurs d'execution avec l'option Check subscrit
!     a la compilation
!
      DIMENSION A(6),R(9)
! ======================================================================
!     GENERATE IDENTITY MATRIX
!
      RANGE=1.0D-12
      IF(MV-1) 10,25,10
   10 IQ=-N
      DO 20 J=1,N
         IQ=IQ+N
         DO 20 I=1,N
            IJ=IQ+I
            R(IJ)=0.D0
            IF(I-J) 20,15,20
   15       R(IJ)=1.D0
   20 CONTINUE
!
!     COMPUTE INITIAL AND FINAL NORMS (ANORM AND ANRMX)
!
   25 ANORM=0.D0
      DO 35 I=1,N
         DO 35 J=I,N
            IF(I-J) 30,35,30
   30       IA=I+(J*J-J)/2
            ANORM=ANORM+A(IA)*A(IA)
   35 CONTINUE
      IF(ANORM) 165,165,40
   40 ANORM=1.414D0*DSQRT(ANORM)
      ANRMX=ANORM*RANGE/DBLE(FLOAT(N))
!
!     INITIALIZE INDICATORS AND COMPUTE THRESHOLD, THR
!
      IND=0
      THR=ANORM
   45 THR=THR/DBLE(FLOAT(N))
   50 L=1
   55 M=L+1
!
!     COMPUTE SIN AND COS
!
   60 MQ=(M*M-M)/2
      LQ=(L*L-L)/2
      LM=L+MQ
      IF(DABS(A(LM))-THR) 130,65,65
   65 IND=1
      LL=L+LQ
      MM=M+MQ
      X=0.5d0*(A(LL)-A(MM))
      Y=-A(LM)/DSQRT(A(LM)*A(LM)+X*X)
      IF(X) 70,75,75
   70 Y=-Y
   75 SINX=Y/DSQRT(2.D0*(1.D0+(DSQRT(1.D0-Y*Y))))
      SINX2=SINX*SINX
      COSX=DSQRT(1.D0-SINX2)
      COSX2=COSX*COSX
      SINCS=SINX*COSX
!
!     ROTATE L AND M COLUMNS
!
      ILQ=N*(L-1)
      IMQ=N*(M-1)
      DO 125 I=1,N
         IQ=(I*I-I)/2
         IF(I-L) 80,115,80
   80    IF(I-M) 85,115,90
   85    IM=I+MQ
         GOTO 95
   90    IM=M+IQ
   95    IF(I-L) 100,105,105
  100    IL=I+LQ
         GOTO 110
  105    IL=L+IQ
  110    X=A(IL)*COSX-A(IM)*SINX
         A(IM)=A(IL)*SINX+A(IM)*COSX
         A(IL)=X
  115    IF(MV-1) 120,125,120
  120    ILR=ILQ+I
         IMR=IMQ+I
         X=R(ILR)*COSX-R(IMR)*SINX
         R(IMR)=R(ILR)*SINX+R(IMR)*COSX
         R(ILR)=X
  125 CONTINUE
      X=2.d0*A(LM)*SINCS
      Y=A(LL)*COSX2+A(MM)*SINX2-X
      X=A(LL)*SINX2+A(MM)*COSX2+X
      A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2)
      A(LL)=Y
      A(MM)=X
!
!     TESTS FOR COMPLETION
!
!     TEST FOR M=LAST COLUMN
!
  130 IF(M-N) 135,140,135
  135 M=M+1
      GOTO 60
!
!     TEST FOR L=SECOND FROM LAST COLUMN
!
  140 IF(L-(N-1)) 145,150,145
  145 L=L+1
      GOTO 55
  150 IF(IND-1) 160,155,160
  155 IND=0
      GOTO 50
!
!     COMPARE THRESHOLD WITH FINAL NORM
!
  160 IF(THR-ANRMX) 165,165,45
!
!     SORT EIGENVALUES AND EIGENVECTORS
!
  165 IQ=-N
      DO 185 I=1,N
         IQ=IQ+N
         LL=(I*I+I)/2
         JQ=N*(I-2)
         DO 185 J=I,N
            JQ=JQ+N
            MM=(J*J+J)/2
            IF(A(LL)-A(MM)) 170,185,185
  170       X=A(LL)
            A(LL)=A(MM)
            A(MM)=X
            IF(MV-1) 175,185,175
  175       DO 180 K=1,N
               ILR=IQ+K
               IMR=JQ+K
               X=R(ILR)
               R(ILR)=R(IMR)
  180       R(IMR)=X
  185 CONTINUE
! ======================================================================
      RETURN
      END
!

! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
! G. Chanteur method - Not tested, problems with simulation-see top
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!

      subroutine inipobary(s)
!
!----------------------------------------------------------------------!
! *   Object: set the initial position of 4 S/C for barycent. coord.
! *   Author: G. Chanteur, P. Robert, CRPE, 1994
!----------------------------------------------------------------------!
!
!
      real s(3,4)
!
!     ------------------------------------------------------------------
      common/satel/ x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4
!     ------------------------------------------------------------------
!
      x1 = s(1,1)
      y1 = s(2,1)
      z1 = s(3,1)
!
      x2 = s(1,2)
      y2 = s(2,2)
      z2 = s(3,2)
!
      x3 = s(1,3)
      y3 = s(2,3)
      z3 = s(3,3)
!
      x4 = s(1,4)
      y4 = s(2,4)
      z4 = s(3,4)
!
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine inimabary(bmag)
!
!----------------------------------------------------------------------!
! *   Object: set the the magnetic field values for barycent. coord.
! *   Author: G. Chanteur, P. Robert, CRPE, 1994
!----------------------------------------------------------------------!
!
!
      real bmag(3,4)
!
!     ------------------------------------------------------------------
      common/field/bx1,by1,bz1,bx2,by2,bz2,bx3,by3,bz3,bx4,by4,bz4
!     ------------------------------------------------------------------
!
!
      bx1 = bmag(1,1)
      by1 = bmag(2,1)
      bz1 = bmag(3,1)
!
      bx2 = bmag(1,2)
      by2 = bmag(2,2)
      bz2 = bmag(3,2)
!
      bx3 = bmag(1,3)
      by3 = bmag(2,3)
      bz3 = bmag(3,3)
!
      bx4 = bmag(1,4)
      by4 = bmag(2,4)
      bz4 = bmag(3,4)
!
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine calbabary(qfbary)
!
!----------------------------------------------------------------------!
! *   Object:  compute barycentric coordinates
! *   Author: G. Chanteur, P. Robert, CRPE, 1994
!----------------------------------------------------------------------!
!
!     input : /satel/
!     output: /base/
!
      real*4      kx1,ky1,kz1,kx2,ky2,kz2,kx3,ky3,kz3,kx4,ky4,kz4
!
!     ------------------------------------------------------------------
      common/satel/ x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4
      common/base /kx1,ky1,kz1,kx2,ky2,kz2,kx3,ky3,kz3,kx4,ky4,kz4,Q
!     ------------------------------------------------------------------
!
!
      ax4 = (y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)
      ay4 = (z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)
      az4 = (x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)

       v4 = ax4*(x4-x1)+ay4*(y4-y1)+az4*(z4-z1)

      if(v4 .lt. 1.e-30) then
         kx4=0.
         ky4=0.
         kz4=0.
                         else
         kx4 = ax4/v4
         ky4 = ay4/v4
         kz4 = az4/v4
      endif
!
      ax3 = (y1-y4)*(z2-z4)-(z1-z4)*(y2-y4)
      ay3 = (z1-z4)*(x2-x4)-(x1-x4)*(z2-z4)
      az3 = (x1-x4)*(y2-y4)-(y1-y4)*(x2-x4)

       v3 = ax3*(x3-x4)+ay3*(y3-y4)+az3*(z3-z4)

      if(V3 .lt. 1.e-30) then
         kx3=0.
         ky3=0.
         kz3=0.
                         else
         kx3 = ax3/v3
         ky3 = ay3/v3
         kz3 = az3/v3
      endif
!
      ax2 = (y4-y3)*(z1-z3)-(z4-z3)*(y1-y3)
      ay2 = (z4-z3)*(x1-x3)-(x4-x3)*(z1-z3)
      az2 = (x4-x3)*(y1-y3)-(y4-y3)*(x1-x3)

       v2 = ax2*(x2-x3)+ay2*(y2-y3)+az2*(z2-z3)

      if(V2 .lt. 1.e-30) then
         kx2=0.
         ky2=0.
         kz2=0.
                         else
         kx2 = ax2/v2
         ky2 = ay2/v2
         kz2 = az2/v2
      endif
!
      ax1 = (y3-y2)*(z4-z2)-(z3-z2)*(y4-y2)
      ay1 = (z3-z2)*(x4-x2)-(x3-x2)*(z4-z2)
      az1 = (x3-x2)*(y4-y2)-(y3-y2)*(x4-x2)

       v1 = ax1*(x1-x2)+ay1*(y1-y2)+az1*(z1-z2)

      if(v1 .lt. 1.e-30) then
         kx1=0.
         ky1=0.
         kz1=0.
                         else
         kx1 = ax1/v1
         ky1 = ay1/v1
         kz1 = az1/v1
      endif
!
!   Facteur de qualite du tetraedre
!
      rk1= kx1**2 + ky1**2 + kz1**2
      rk2= kx2**2 + ky2**2 + kz2**2
      rk3= kx3**2 + ky3**2 + kz3**2
      rk4= kx4**2 + ky4**2 + kz4**2
!
      rk2min=min(rk1,rk2,rk3,rk4)
      rk2max=max(rk1,rk2,rk3,rk4)
!
      if(rk2max .lt. 1.e-30) then
         qfbary=0.
                             else
         qfbary=rk2min/rk2max
      endif
!
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine estibbary
!
!----------------------------------------------------------------------!
! *   Object: Estimate vectorial quantities of B by barycent. coord.
! *   Author: G.Chanteur, CETP, 1994
!----------------------------------------------------------------------!
!
!     ******************************************************************
!     Estimate vectorial quantities of B by the barycentric coordinates
!     0riginal sofware delivered by G. Chanteur CETP/IPSL - CNRS & UVSQ,
!     september 1994. Radius of curvature from gradB, G.C., 2003/02/26.
!
!     input : /satel/
!             /base /
!             /field/
!     output: /estim/
!     ******************************************************************
!
      real*4      kx1,ky1,kz1,kx2,ky2,kz2,kx3,ky3,kz3,kx4,ky4,kz4
!
!     ------------------------------------------------------------------
      common/base /kx1,ky1,kz1,kx2,ky2,kz2,kx3,ky3,kz3,kx4,ky4,kz4,Q
      common/field/bx1,by1,bz1,bx2,by2,bz2,bx3,by3,bz3,bx4,by4,bz4
      common/estim/gbxx,gbxy,gbxz,gbyx,gbyy,gbyz,gbzx,gbzy,gbzz &
                  ,divb,rotbx,rotby,rotbz,rnormx,rnormy,rnormz,radcurv
!     ------------------------------------------------------------------
!
!
! *** Estimation du tenseur gradB
!
      gbxx   = kx1*bx1 +kx2*bx2 +kx3*bx3 +kx4*bx4
      gbxy   = kx1*by1 +kx2*by2 +kx3*by3 +kx4*by4
      gbxz   = kx1*bz1 +kx2*bz2 +kx3*bz3 +kx4*bz4
!
      gbyx   = ky1*bx1 +ky2*bx2 +ky3*bx3 +ky4*bx4
      gbyy   = ky1*by1 +ky2*by2 +ky3*by3 +ky4*by4
      gbyz   = ky1*bz1 +ky2*bz2 +ky3*bz3 +ky4*bz4
!
      gbzx   = kz1*bx1 +kz2*bx2 +kz3*bx3 +kz4*bx4
      gbzy   = kz1*by1 +kz2*by2 +kz3*by3 +kz4*by4
      gbzz   = kz1*bz1 +kz2*bz2 +kz3*bz3 +kz4*bz4
!
! *** Estimation de divB
!
      divb = gbxx+gbyy+gbzz
!
! *** Estimation de rotB
!
      rotbx = gbyz-gbzy
      rotby = gbzx-gbxz
      rotbz = gbxy-gbyx

! *** Estimation of the unit vector normal to field lines through
!     CLUSTER  and of the radius of curvature from gradB
!     G. Chanteur CETP/IPSL - CNRS & UVSQ 2003/02/26 at IASB, Brussels
!
! *** Magnetic field at center of mass and unit tangent vector T

      bmean_x = 0.25*(bx1 +bx2 +bx3 +bx4)
      bmean_y = 0.25*(by1 +by2 +by3 +by4)
      bmean_z = 0.25*(bz1 +bz2 +bz3 +bz4)

      bmean_m = sqrt(bmean_x**2 +bmean_y**2 +bmean_z**2)
!
      if(bmean_m .lt. 1.e-30) then
         tx=0.
         ty=0.
         tz=0.
                          else
         tx = bmean_x/bmean_m
         ty = bmean_y/bmean_m
         tz = bmean_z/bmean_m
      endif

! *** Estimation of the unit normal vector N and curvature 1/R

      tk1 = tx*kx1 +ty*ky1 +tz*kz1
      tk2 = tx*kx2 +ty*ky2 +tz*kz2
      tk3 = tx*kx3 +ty*ky3 +tz*kz3
      tk4 = tx*kx4 +ty*ky4 +tz*kz4
!
      dbds   = (tx*bx1 +ty*by1 +tz*bz1)*tk1 &
              +(tx*bx2 +ty*by2 +tz*bz2)*tk2 &
              +(tx*bx3 +ty*by3 +tz*bz3)*tk3 &
              +(tx*bx4 +ty*by4 +tz*bz4)*tk4
!
      Ax = tk1*bx1 +tk2*bx2 +tk3*bx3 +tk4*bx4 - dbds*tx
      Ay = tk1*by1 +tk2*by2 +tk3*by3 +tk4*by4 - dbds*ty
      Az = tk1*bz1 +tk2*bz2 +tk3*bz3 +tk4*bz4 - dbds*tz
      Am = sqrt(Ax**2 +Ay**2 +Az**2)
!
      if(Am .lt. 1.e-30) then
        rnormx=0.
        rnormy=0.
        rnormz=0.
        radcurv=0.
                         else
        rnormx= Ax/Am
        rnormy= Ay/Am
        rnormz= Az/Am
        radcurv = bmean_m/Am
      endif

      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine ccurlbary(curx,cury,curz)
!
!----------------------------------------------------------------------!
! *   Object: compute curl(B) from barycentric method
! *   Author: G. Chanteur, P. Robert, CRPE, 1994
!----------------------------------------------------------------------!
!
!
!     ------------------------------------------------------------------
      common/estim/gbxx,gbxy,gbxz,gbyx,gbyy,gbyz,gbzx,gbzy,gbzz &
                  ,divb,rotbx,rotby,rotbz,rnormx,rnormy,rnormz,radcurv
!     ------------------------------------------------------------------
!
!
!     Estimation de rotB
!
      curx = rotbx
      cury = rotby
      curz = rotbz
!
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine ccurvbary(rnormax,rnormay,rnormaz,rcurv)
!
!----------------------------------------------------------------------!
! *   Object: compute radius of curvature of field lines
!             and normal to osculateur plane from barycentric method
! *   Author: G.Chanteur, CETP, 2003
!----------------------------------------------------------------------!
!
!
!     ------------------------------------------------------------------
      common/estim/gbxx,gbxy,gbxz,gbyx,gbyy,gbyz,gbzx,gbzy,gbzz &
                  ,divb,rotbx,rotby,rotbz,rnormx,rnormy,rnormz,radcurv
!     ------------------------------------------------------------------
!
!
!     Estimation de la normale
!
      rnormax = rnormx
      rnormay = rnormy
      rnormaz = rnormz
!
!     Estimation du rayon de courbure
!
      rcurv=radcurv
!
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cgrabary(gra1x,gra1y,gra1z, gra2x,gra2y,gra2z, &
                          gra3x,gra3y,gra3z)
!
!----------------------------------------------------------------------!
! *   Object: compute grad(B) matrix from barycentric metho
! *   Author: G.Chanteur, CETP, 1994
!----------------------------------------------------------------------!
!
!
!     ------------------------------------------------------------------
      common/estim/gbxx,gbxy,gbxz,gbyx,gbyy,gbyz,gbzx,gbzy,gbzz &
                  ,divb,rotbx,rotby,rotbz,rnormx,rnormy,rnormz,radcurv
!     ------------------------------------------------------------------
!
!
!     Estimation du tenseur gradB
!
      gra1x = gbxx
      gra1y = gbyx
      gra1z = gbzx

      gra2x = gbxy
      gra2y = gbyy
      gra2z = gbzy

      gra3x = gbxz
      gra3y = gbyz
      gra3z = gbzz
!
!
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine cdivbary(divdeb)
!
!----------------------------------------------------------------------!
! *   Object: compute div(B) from barycentric method
! *   Author: G. Chanteur, P. Robert, CRPE, 1994
!----------------------------------------------------------------------!
!
!
!     ------------------------------------------------------------------
      common/estim/gbxx,gbxy,gbxz,gbyx,gbyy,gbyz,gbzx,gbzy,gbzz &
                  ,divb,rotbx,rotby,rotbz,rnormx,rnormy,rnormz,radcurv
!     ------------------------------------------------------------------
!
!
!     Estimation de divB
!
      divdeb = divb
!
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine ctestbary(iok)
!
!----------------------------------------------------------------------!
! *   Object: compute test on gradient sign from barycentric method
! *   Author: G. Chanteur, P. Robert, CRPE, 1994
!----------------------------------------------------------------------!
!
!
!     ------------------------------------------------------------------
      common/estim/gbxx,gbxy,gbxz,gbyx,gbyy,gbyz,gbzx,gbzy,gbzz &
                  ,divb,rotbx,rotby,rotbz,rnormx,rnormy,rnormz,radcurv
!     ------------------------------------------------------------------
!
      iok=1
!
      if(gbyz*gbzy.gt.0.) iok=0
      if(gbzx*gbxz.gt.0.) iok=0
      if(gbxy*gbyx.gt.0.) iok=0
!
      return
      end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
