      program rocot_check_mva 
!                                                                       
!     ***************************************************************0**
! *   Rocot_Check_MVA: program de test des modules varmin de le rocotlib
!                                                                       
!     calcul la base (vecteurs propres) de la matrice de variance       
!     d'un signal bx,by,bz                                              
!                                                                       
!     P. Robert, aout 2001                                              
!     ***************************************************************0**
!                                                                       
      parameter (nbp= 1000) 
!                                                                       
      dimension bx(nbp),by(nbp),bz(nbp) 
      dimension var(3,3),valp(3),vecp(3,3) 
!                                                                       
      print 200 
      print 200, '=======================' 
      print 200, 'rocot_check_mva program' 
      print 200, '=======================' 
      print 200 
!                                                                       
! *** chargement d'un signal test                                       
!                                                                       
      call sigtest(bx,by,bz,nbp) 
!                                                                       
! *** ouverture du fichier resu et chargement du signal d'entree        
!                                                                       
      open(1,file='rocot_check_mva.resu') 
!                                                                       
      call writ_para_tube(1) 
!                                                                       
      write(1,'(a)') 'Input signal:' 
      write(1,300) nbp 
      write(1,100) (bx(i),by(i),bz(i),i=1,nbp) 
!                                                                       
! *** calcul de variance du signal:                                     
!     calcul de la matrice de variance, des  valeurs propres et         
!     des vecteurs propres ; verification de l'orthogonalite            
!     et du sens de la base trouvee, et passage des donnees             
!     dans cette nouvelle base si icrep NE 0.                           
!                                                                       
      icrep=1 
!                                                                       
      call mat_cp_varmin(6,bx,by,bz,nbp,icrep,var,valp,vecp) 
!                                                                       
! *** ecriture du signal dans la base de la variance                    
!                                                                       
      write(1,'(a)') '-------------------------------------' 
      write(1,'(a)') 'MVA signal:' 
      write(1,100) (bx(i),by(i),bz(i),i=1,nbp) 
!                                                                       
! *** fermeture du fichier resultat                                     
!                                                                       
      close(1) 
!                                                                       
      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, '=======================' 
!                                                                       
  100 format(3(1Pe14.5)) 
  200 format(a) 
  300 format(i8) 
!                                                                       
      stop 'rocot_check_mva.exe: normal termination' 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine sigtest(bx,by,bz,nbp) 
!                                                                       
      real bx(nbp),by(nbp),bz(nbp) 
!                                                                       
      common /tube/ grtube,pari,gjtube,tettu,phitu,gjxtu,gjytu,gjztu 
!                                                                       
      call read_para_tube 
!                                                                       
      dx=grtube/60. 
!                                                                       
      print 200 
      print 200,'Computation of B vector along the trajectory' 
      print 200,'--------------------------------------------' 
      print 200 
!                                                                       
      print 200, 'nb. of points=',nbp 
!                                                                       
      call writ_para_tube(6) 
!                                                                       
      do i=1,nbp 
!                                                                       
         sx= float(i-1)*dx -10.*pari 
         sy= -pari 
         sz= 0. 
!                                                                       
         call cal_b_tube(sx,sy,sz,bx(i),by(i),bz(i),bmoy) 
!                                                                       
      enddo 
                                                                        
  200 format(a,i9) 
                                                                        
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
!     magsructlib: "magnetic structure libarary"                        
!                  bibliotheque de modules permettant le calcul         
!                  du champ et de la densite de courant produit         
!                  par une structure de courant telle un tube de courant
!                  a densite uniforme ou gaussienne.                    
!                                                                       
!                  Biblioteque developpee pour les simulation CLUSTER   
!                  Patrick ROBERT, CRPE, 1990-1995                      
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine read_para_tube 
!                                                                       
      common /tube/ grtube,pari,gjtube,tettu,phitu,gjxtu,gjytu,gjztu 
      common /repaj/ modj 
!                                                                       
      character*1 modj 
!                                                                       
!     ***********************************************************       
!     lecture des parametres du tube                                    
!     ***********************************************************       
!                                                                       
      print 200, ' ' 
      print 200, 'Current tube characteristics:' 
      print 200, '-----------------------------' 
!                                                                       
      print 200, ' ' 
      print 200, 'Radius of the current tube, in km' 
      read *,   grtube 
      print 300,grtube 
!                                                                       
      print 200, ' ' 
      print 200, 'impact parameter (min dist to the curr. tube), in km' 
      read *,pari 
      print 300,pari 
!                                                                       
      print 200, ' ' 
      print 200, 'Current density, in A/km2' 
      read *,   gjtube 
      print 300,gjtube 
!                                                                       
      print 200, ' ' 
      print 200, 'Theta angle of the current direction/ given frame (d)' 
      print 200, 'teta ?' 
      read *,   tettu 
      print 300,tettu 
!                                                                       
      print 200, ' ' 
      print 200, 'Phi   angle of the current direction/ given frame (d)' 
      read *,   phitu 
      print 300,phitu 
!                                                                       
      call carsphe(gjtube,tettu,phitu,gjxtu,gjytu,gjztu) 
!                                                                       
!                                                                       
      print 200, ' ' 
      print 200, 'Current density model ("u" or "g"/ uniform or gauss.)' 
      read (5,100) modj 
      print   100, modj 
      print 200 
!                                                                       
  100 format(a1) 
  200 format(a) 
  300 format(f10.3) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine writ_para_tube(ifc) 
!                                                                       
      common /tube/ grtube,pari,gjtube,tettu,phitu,gjxtu,gjytu,gjztu 
      common /repaj/ modj 
!                                                                       
      character*1 modj 
!                                                                       
!     ***********************************************************       
!     ecriture des parametres du tube                                   
!     ***********************************************************       
!                                                                       
      write(ifc,100)  grtube, pari, gjtube, tettu, phitu, modj 
!                                                                       
  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)                       
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cal_b_tube(xs,ys,zs,bxs,bys,bzs,bmod) 
!                                                                       
      common /tube/ grtube,pari,gjtube,tettu,phitu,gjxtu,gjytu,gjztu 
      common /repaj/ modj 
!                                                                       
      character*1 modj 
!                                                                       
!     ***********************************************************       
!     calcul de B dans le repere structure                              
!     repere structure=repere ou le tube est oblique                    
!     common /tube/=axe et module de J dans le repere structure         
!     ***********************************************************       
!                                                                       
      pi=acos(-1.) 
!     valeur de Mu0 en nT*km/A                                          
!     rmu=4.*pi*(1.e-7)*(1.e6)                                          
!     c=Mu0*J/2                                                         
      c=0.2*pi*gjtube 
!                                                                       
!                                                                       
! *** passage au repere tube incline par rapport au repere structure    
!                                                                       
      call xybvdh(xs,ys,zs,gjxtu,gjytu,gjztu,x,y,z) 
!                                                                       
      r2=x*x+y*y 
      gr2=grtube*grtube 
!                                                                       
!                                                                       
      if(modj.ne.'u' .and. modj.ne.'g')                                 &
     &         stop 'CAL_B_TUBE     *** ABORTED ! UNDEFINED MODEL OF J' 
                                                                        
      if(modj.eq.'u') then 
!                                                                       
! *** densite uniforme                                                  
!                                                                       
          if(r2.ge.gr2) then 
                        bx=-c*gr2*y/r2 
                        by= c*gr2*x/r2 
                        bz=0. 
!                                                                       
                        else 
                        bx=-c*y 
                        by= c*x 
                        bz= 0. 
          endif 
!                                                                       
      else 
!                                                                       
! *** densite gaussienne                                                
!                                                                       
      rsr2=r2/gr2 
!                                                                       
      bx=-c*(1./rsr2)*(1.-exp(-rsr2))*y 
      by= c*(1./rsr2)*(1.-exp(-rsr2))*x 
      bz= 0. 
!                                                                       
!     um=1.06                                                           
!     rm=sqrt(um)*grtube                                                
!     bmax=c*rm*(1.-exp(-um))/um                                        
!                                                                       
      endif 
!                                                                       
      bmod=sqrt(bx*bx+by*by+bz*bz) 
!                                                                       
! *** passage au repere structure                                       
!                                                                       
      call vdhxyb(bx,by,bz,gjxtu,gjytu,gjztu,bxs,bys,bzs) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine cal_j_tube(xs,ys,zs,xsj,ysj,zsj,rsj) 
!                                                                       
      common /tube/ grtube,pari,gjtube,tettu,phitu,gjxtu,gjytu,gjztu 
      common /repaj/ modj 
!                                                                       
      character*1 modj 
!                                                                       
!     ***********************************************************       
!     calcul de J dans le repere structure                              
!     repere structure=repere ou le tube est oblique                    
!     common /tube/=axe et module de J dans le repere structure         
!     ***********************************************************       
!                                                                       
!                                                                       
! *** passage au repere tube                                            
!                                                                       
      call xybvdh(xs,ys,zs,gjxtu,gjytu,gjztu,x,y,z) 
!                                                                       
      r=sqrt(x*x+y*y) 
!                                                                       
      if(modj.ne.'u' .and. modj.ne.'g')                                 &
     &         stop 'CAL_B_TUBE     *** ABORTED ! UNDEFINED MODEL OF J' 
                                                                        
      if(modj.eq.'u') then 
!                                                                       
! *** modele uniforme                                                   
!                                                                       
          if(r.gt.grtube) then 
                          xj=0. 
                          yj=0. 
                          zj=0. 
!                                                                       
                          else 
                          xj=0. 
                          yj=0. 
                          zj=gjtube 
                          endif 
!                                                                       
      else 
!                                                                       
! *** modele gaussien                                                   
!                                                                       
      xj=0. 
      yj=0. 
      zj=gjtube*exp(-(r/grtube)**2) 
!                                                                       
      endif 
                                                                        
      rsj=sqrt(xj*xj+yj*yj+zj*zj) 
!                                                                       
! *** passage au repere structure                                       
!                                                                       
      call vdhxyb(xj,yj,zj,gjxtu,gjytu,gjztu,xsj,ysj,zsj) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine carsphe(r,tet,phi,x,y,z) 
!                                                                       
!     ******************************************************************
!     auteur   :p. robert                                               
!     categorie:astronomie et changement de base                        
!     objet    :conversion de coordonnees spheriques en cartesiennes    
!     ******************************************************************
!                                                                       
      pp=3.1415926/180. 
!                                                                       
      x=r*sin(tet*pp)*cos(phi*pp) 
      y=r*sin(tet*pp)*sin(phi*pp) 
      z=r*cos(tet*pp) 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine sincosp(x,y,z,st,ct,sp,cp) 
!                                                                       
!     ******************************************************************
!     auteur   :p. robert                                               
!     categorie:astronomie et changement de base                        
!     objet    :sin et cos des angles polaires d"un vecteur cartesien   
!     ******************************************************************
!                                                                       
      r=sqrt(x*x+y*y+z*z) 
!                                                                       
      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 
!                                                                       
      sp=0. 
      cp=1. 
      rp=sqrt(x*x+y*y) 
!                                                                       
      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 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine vdhxyb(x,y,b,bv,bd,bh,v,d,h) 
!                                                                       
!     ******************************************************************
!     auteur   :p. robert                                               
!     categorie:astronomie et changement de base                        
!     objet    :passage du repere ligne de force au repere vdh          
!     ******************************************************************
!                                                                       
!      input=x,y,b avec axe x dans le plan du meridien magnetique       
!      output=v,d,h                                                     
!                                                                       
      call sincosp(bv,bd,bh,st,ct,sp,cp) 
!                                                                       
      stv=sqrt(1.-cp*cp*st*st) 
      if(stv.gt.1.e-30) then 
                        spv=-sp/stv 
                        cpv=cp*ct/stv 
!                                                                       
                        else 
                        spv=0. 
                        cpv=1. 
                        endif 
!                                                                       
      b1=cpv*x -spv*y 
      b2=spv*x +cpv*y 
      b3=b 
!                                                                       
      v= cp*ct*b1 -sp*b2 +cp*st*b3 
      d= sp*ct*b1 +cp*b2 +sp*st*b3 
      h=-st*b1 +ct*b3 
!                                                                       
      return 
      END                                           
!                                                                       
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!                                                                       
      subroutine xybvdh(v,d,h,bv,bd,bh,x,y,b) 
!                                                                       
!     ******************************************************************
!     auteur   :p. robert                                               
!     categorie:astronomie et changement de base                        
!     objet    :passage du repere vdh au repere ligne de force          
!     ******************************************************************
!                                                                       
!      input=v,d,h                                                      
!      output=x,y,b                                                     
!             x est dans le plan du meridien magnetique                 
!                                                                       
!                                                                       
      call sincosp(bv,bd,bh,st,ct,sp,cp) 
!                                                                       
      b1= ct*cp*v +ct*sp*d -st*h 
      b2=-sp*v +cp*d 
      b3= st*cp*v +st*sp*d +ct*h 
!                                                                       
      stv=sqrt(1.-cp*cp*st*st) 
      if(stv.gt.1.e-30 ) then 
                         spv=-sp/stv 
                         cpv=cp*ct/stv 
!                                                                       
                         else 
                         spv=0. 
                         cpv=1. 
                         endif 
!                                                                       
      x= cpv*b1 +spv*b2 
      y=-spv*b1 +cpv*b2 
      b=b3 
!                                                                       
      return 
      END                                           
