      program rocot_check_mva
c
c     ***************************************************************0**
c *   Rocot_Check_MVA: program de test des modules varmin de le rocotlib
c
c     calcul la base (vecteurs propres) de la matrice de variance
c     d'un signal bx,by,bz
c
c     P. Robert, aout 2001
c     ***************************************************************0**
c
      parameter (nbp= 1000)
c
      dimension bx(nbp),by(nbp),bz(nbp)
      dimension var(3,3),valp(3),vecp(3,3)
c
      print 200
      print 200, '======================='
      print 200, 'rocot_check_mva program'
      print 200, '======================='
      print 200
c
c *** chargement d'un signal test
c
      call sigtest(bx,by,bz,nbp)
c
c *** ouverture du fichier resu et chargement du signal d'entree
c
      open(1,file='rocot_check_mva.resu')
c
      call writ_para_tube(1)
c
      write(1,'(a)') 'Input signal:'
      write(1,300) nbp
      write(1,100) (bx(i),by(i),bz(i),i=1,nbp)
c
c *** calcul de variance du signal:
c     calcul de la matrice de variance, des  valeurs propres et
c     des vecteurs propres ; verification de l'orthogonalite
c     et du sens de la base trouvee, et passage des donnees
c     dans cette nouvelle base si icrep NE 0.
c
      icrep=1
c
      call mat_cp_varmin(6,bx,by,bz,nbp,icrep,var,valp,vecp)
c
c *** ecriture du signal dans la base de la variance
c
      write(1,'(a)') '-------------------------------------'
      write(1,'(a)') 'MVA signal:'
      write(1,100) (bx(i),by(i),bz(i),i=1,nbp)
c
c *** fermeture du fichier resultat
c
      close(1)
c
      print 200, ' '
      print 200, 'Input signal and MVA coord. syst. signal are in '//
     &        'rocot_check_mva.resu file'
      print 200
      print 200, '======================='
      print 200, 'end of rocot_check_mva '
      print 200, '======================='
c
  100 format(3(1Pe14.5))
  200 format(a)
  300 format(i8)
c
      stop 'rocot_check_mva.exe: normal termination'
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine sigtest(bx,by,bz,nbp)
c
      real bx(nbp),by(nbp),bz(nbp)
c
      common /tube/ grtube,pari,gjtube,tettu,phitu,gjxtu,gjytu,gjztu
c
      call read_para_tube
c
      dx=grtube/60. 
c
      print 200
      print 200,'Computation of B vector along the trajectory'
      print 200,'--------------------------------------------'
      print 200
c
      print 200, 'nb. of points=',nbp
c
      call writ_para_tube(6)
c
      do i=1,nbp
c
         sx= float(i-1)*dx -10.*pari
         sy= -pari
         sz= 0.
c
         call cal_b_tube(sx,sy,sz,bx(i),by(i),bz(i),bmoy)
c
      enddo

  200 format(a,i9)
 
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
c     magsructlib: "magnetic structure libarary"
c                  bibliotheque de modules permettant le calcul
c                  du champ et de la densite de courant produit
c                  par une structure de courant telle un tube de courant
c                  a densite uniforme ou gaussienne.
c
c                  Biblioteque developpee pour les simulation CLUSTER
c                  Patrick ROBERT, CRPE, 1990-1995
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine read_para_tube
c
      common /tube/ grtube,pari,gjtube,tettu,phitu,gjxtu,gjytu,gjztu
      common /repaj/ modj
c
      character*1 modj
c
c     ***********************************************************
c     lecture des parametres du tube
c     ***********************************************************
c
      print 200, ' '
      print 200, 'Current tube characteristics:'
      print 200, '-----------------------------'
c
      print 200, ' '
      print 200, 'Radius of the current tube, in km'
      read *,   grtube
      print 300,grtube
c
      print 200, ' '
      print 200, 'impact parameter (min dist to the curr. tube), in km'
      read *,pari
      print 300,pari
c
      print 200, ' '
      print 200, 'Current density, in A/km2'
      read *,   gjtube
      print 300,gjtube
c
      print 200, ' '
      print 200, 'Theta angle of the current direction/ given frame (d)'
      print 200, 'teta ?'
      read *,   tettu
      print 300,tettu
c
      print 200, ' '
      print 200, 'Phi   angle of the current direction/ given frame (d)'
      read *,   phitu
      print 300,phitu
c
      call carsphe(gjtube,tettu,phitu,gjxtu,gjytu,gjztu)
c
c
      print 200, ' '
      print 200, 'Current density model ("u" or "g"/ uniform or gauss.)'
      read (5,100) modj
      print   100, modj
      print 200
c
  100 format(a1)
  200 format(a)
  300 format(f10.3)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine writ_para_tube(ifc)
c
      common /tube/ grtube,pari,gjtube,tettu,phitu,gjxtu,gjytu,gjztu
      common /repaj/ modj
c
      character*1 modj
c
c     ***********************************************************
c     ecriture des parametres du tube
c     ***********************************************************
c
      write(ifc,100)  grtube, pari, gjtube, tettu, phitu, modj
c
  100 format(' R=',f8.1,' km  IP=',f8.1,' km J=',f6.3,' A/km2  tet=', 
     &       f5.1,' d.  phi=',f5.1,' d. mod=',a1)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cal_b_tube(xs,ys,zs,bxs,bys,bzs,bmod)
c
      common /tube/ grtube,pari,gjtube,tettu,phitu,gjxtu,gjytu,gjztu
      common /repaj/ modj
c
      character*1 modj
c
c     ***********************************************************
c     calcul de B dans le repere structure
c     repere structure=repere ou le tube est oblique
c     common /tube/=axe et module de J dans le repere structure
c     ***********************************************************
c
      pi=acos(-1.)
c     valeur de Mu0 en nT*km/A
c     rmu=4.*pi*(1.e-7)*(1.e6)
c     c=Mu0*J/2
      c=0.2*pi*gjtube
c
c
c *** passage au repere tube incline par rapport au repere structure
c
      call xybvdh(xs,ys,zs,gjxtu,gjytu,gjztu,x,y,z)
c
      r2=x*x+y*y
      gr2=grtube*grtube
c
c
      if(modj.ne.'u' .and. modj.ne.'g')
     &         stop 'CAL_B_TUBE     *** ABORTED ! UNDEFINED MODEL OF J'

      if(modj.eq.'u') then
c
c *** densite uniforme
c
          if(r2.ge.gr2) then
                        bx=-c*gr2*y/r2
                        by= c*gr2*x/r2
                        bz=0.
c
                        else
                        bx=-c*y
                        by= c*x
                        bz= 0.
          endif
c
      else
c
c *** densite gaussienne
c
      rsr2=r2/gr2
c
      bx=-c*(1./rsr2)*(1.-exp(-rsr2))*y
      by= c*(1./rsr2)*(1.-exp(-rsr2))*x
      bz= 0.
c
c     um=1.06
c     rm=sqrt(um)*grtube
c     bmax=c*rm*(1.-exp(-um))/um
c
      endif
c
      bmod=sqrt(bx*bx+by*by+bz*bz)
c
c *** passage au repere structure
c
      call vdhxyb(bx,by,bz,gjxtu,gjytu,gjztu,bxs,bys,bzs)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine cal_j_tube(xs,ys,zs,xsj,ysj,zsj,rsj)
c
      common /tube/ grtube,pari,gjtube,tettu,phitu,gjxtu,gjytu,gjztu
      common /repaj/ modj
c
      character*1 modj
c
c     ***********************************************************
c     calcul de J dans le repere structure
c     repere structure=repere ou le tube est oblique
c     common /tube/=axe et module de J dans le repere structure
c     ***********************************************************
c
c
c *** passage au repere tube
c
      call xybvdh(xs,ys,zs,gjxtu,gjytu,gjztu,x,y,z)
c
      r=sqrt(x*x+y*y)
c
      if(modj.ne.'u' .and. modj.ne.'g')
     &         stop 'CAL_B_TUBE     *** ABORTED ! UNDEFINED MODEL OF J'

      if(modj.eq.'u') then
c
c *** modele uniforme
c
          if(r.gt.grtube) then
                          xj=0.
                          yj=0.
                          zj=0.
c 
                          else
                          xj=0.
                          yj=0.
                          zj=gjtube
                          endif
c
      else
c
c *** modele gaussien
c
      xj=0.
      yj=0.
      zj=gjtube*exp(-(r/grtube)**2)
c
      endif
      
      rsj=sqrt(xj*xj+yj*yj+zj*zj)
c
c *** passage au repere structure
c
      call vdhxyb(xj,yj,zj,gjxtu,gjytu,gjztu,xsj,ysj,zsj)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine carsphe(r,tet,phi,x,y,z)
c
c     ******************************************************************
c     auteur   :p. robert
c     categorie:astronomie et changement de base
c     objet    :conversion de coordonnees spheriques en cartesiennes
c     ******************************************************************
c
      pp=3.1415926/180.
c
      x=r*sin(tet*pp)*cos(phi*pp)
      y=r*sin(tet*pp)*sin(phi*pp)
      z=r*cos(tet*pp)
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine sincosp(x,y,z,st,ct,sp,cp)
c
c     ******************************************************************
c     auteur   :p. robert
c     categorie:astronomie et changement de base
c     objet    :sin et cos des angles polaires d"un vecteur cartesien
c     ******************************************************************
c
      r=sqrt(x*x+y*y+z*z)
c
      ct=z/r
      if(abs(ct).gt.1.) then
      stop 'SINCOSP        *** ABORTED ! ERROR 1'
      endif
      
      st=sqrt(1.-ct*ct)
      if(abs(st).gt.1.) then
      stop 'SINCOSP        *** ABORTED ! ERROR 2'
      endif
c
      sp=0.
      cp=1.
      rp=sqrt(x*x+y*y)
c
      if(rp.lt.1.e-30) return
      sp=y/(rp)
      if(abs(sp).gt.1.) then
      stop 'SINCOSP        *** ABORTED ! ERROR 3'
      endif
      
      cp=x/(rp)
      if(abs(cp).gt.1.) then
      stop 'SINCOSP        *** ABORTED ! ERROR 4'
      endif
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine vdhxyb(x,y,b,bv,bd,bh,v,d,h)
c
c     ******************************************************************
c     auteur   :p. robert
c     categorie:astronomie et changement de base
c     objet    :passage du repere ligne de force au repere vdh
c     ******************************************************************
c
c      input=x,y,b avec axe x dans le plan du meridien magnetique
c      output=v,d,h
c
      call sincosp(bv,bd,bh,st,ct,sp,cp)
c
      stv=sqrt(1.-cp*cp*st*st)
      if(stv.gt.1.e-30) then
                        spv=-sp/stv
                        cpv=cp*ct/stv
c
                        else
                        spv=0.
                        cpv=1.
                        endif
c
      b1=cpv*x -spv*y
      b2=spv*x +cpv*y
      b3=b
c
      v= cp*ct*b1 -sp*b2 +cp*st*b3
      d= sp*ct*b1 +cp*b2 +sp*st*b3
      h=-st*b1 +ct*b3
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
      subroutine xybvdh(v,d,h,bv,bd,bh,x,y,b)
c
c     ******************************************************************
c     auteur   :p. robert
c     categorie:astronomie et changement de base
c     objet    :passage du repere vdh au repere ligne de force
c     ******************************************************************
c
c      input=v,d,h
c      output=x,y,b
c             x est dans le plan du meridien magnetique
c
c
      call sincosp(bv,bd,bh,st,ct,sp,cp)
c
      b1= ct*cp*v +ct*sp*d -st*h
      b2=-sp*v +cp*d
      b3= st*cp*v +st*sp*d +ct*h
c
      stv=sqrt(1.-cp*cp*st*st)
      if(stv.gt.1.e-30 ) then
                         spv=-sp/stv
                         cpv=cp*ct/stv
c
                         else
                         spv=0.
                         cpv=1.
                         endif
c
      x= cpv*b1 +spv*b2
      y=-spv*b1 +cpv*b2
      b=b3
c
      return
      end
c
c     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
c
