      program coanim_magsphe_cluorbit
!
! ----------------------------------------------------------------------
!
!     trace d'une carte meridienne des lignes de force
!     du modele de tsyganenko., pour chaque intervalle de temps,
!     ainsi que la trajectoire et les positions des 4 sat. CLUSTER,
!     et production des utilitaires pour l'animation finale.
!
!     P. Robert, CETP, 1996
!     Revision Mars 2003 pour version 2001 de la comagnetolib.
!     Revision Juin 2003 pour generalisation.
!     Revision Janv 2012 pour compatibilite coplot_magsphe et Windows
!
!                     ---------------
!
!     Ce programme est derive de la souche comagneto.f, mais avec
!     une boucle sur le temps pour l'animation.
!     Il genere une serie de .ps pour  les differents pas de temps.
!     La sortie est en mode portrait pour pouvoir convertir
!     les .ps en image apres filtre de position/selection.
!     le repere de sortie est gsm ou gse
!
!     Le soleil est a gauche (isun=-1) ou a droite (isun=+1),
!     On a la Terre et le bowshock.
!     les couleurs de remplissage sont choisies pour mettre en valeur
!     les zones d'interet
!
!     On trace une carte meridienne des lignes de force
!     du modele de tsyganenko.
!
!     Calcul selon ecrements axe X le long de SM pour les lignes fermees
!     (et donc le long d'un axe X perpendiculaire au dipole Z et
!     avec la direction du soleil dans le plan XZ)
!     et le long d'un axe // a Z GSM pour les lignes ouvertes
!     (et donc perpendiculaire a la direction du soleil)
!     Rappel: SM et GSM ont Y en commun.
!     Le depart des lignes de force est dans le plan XZ.
!
!     Utilise les bibliotheques suivantes:
!          cluster_orbitlib.o
!          comagnetolib.o
!          geopacklib.o
!          maglinelib.o
!          rocotlib.o
!          rogralib_V9p6.o
!          tsauxlib.o
!          tsygalib.o
!          worldlib.o
!     ------------------------------------------------------------------
!
! 
! coanim_magsphe_cluorbit.sh
! --------------------------
! 
! 1) detruit les eventuels directories resu_ps, resu_png et resu_x
! 
! 2) recrer ces directories vides
!
!    resu_ps  pour les amag_XXX.ps
!    resu_png pour les amag_XXX.png
! 
! 3) lance coanim_magsphe_cluorbit.exe <inimagneto.in >coanim_magsphe_cluorbit.out
! 
! 
! coanim_magsphe_cluorbit.exe:
! ----------------------------
! 
! 1) production en chaine des amag_XXX.ps
! 
! 2) produit le conv_all_ps_to_png.sh de conversion des .ps en .png:
! 
!
!                     -----------------------------------
!
! *** declaration pour le co_plot_magneto
!
 
      real(kind=4), dimension(10) :: parmod

      character*6 ex
      character*3 rep
      character*5 sunpo
      character*2 stipox,stipoy,ckp
      character*80 orbidir, bid, cmon, ctime

      character(len=1) :: cocli,cooli,codip,corot,coarr

!     real (kind=4), dimension(3,1200,20) :: linjou,linnui,linnor,linsud
      real (kind=4), dimension(3,2000,20) :: linjou,linnui,linnor,linsud
      integer,       dimension(20)        :: npljou,nplnui,nplnor,nplsud

!     real (kind=4), dimension(3,3600) :: magpau,bowsho ! 3600=3*1200
      real (kind=4), dimension(3,6000) :: magpau,bowsho ! 6000=3*2000

      integer, dimension(180,360) :: icowor
      integer, dimension(101,101) :: imawor
!
! *** declarations propres a coanim_magsphe_cluorbit
!
      character*3 ch3
      character*21 fileps,filpng
      character*64 com

      common /sunsid/ isun
      common /dimlin/ Nlimax,Nptmax

!     data Nlimax, Nptmax /20, 600/ ! 600*2=1200 points (// et \//)
      data Nlimax, Nptmax /20,1000/ ! 1000*2=2000 points (// et \//)
      data idiwor /101/
!
! *   taille et positions de la figure (cm)
!     (attention echelles, pour coherence sfx/dx = sfy/dy)
!
      data sfx,sfy /18. ,18./
      data ofx,ofy / 1. , 1./
!
      data x1,x2   /-60., 30./
      data y1,y2   /-40., 50./
!
! *   position des fleches pour le soleil (en RT)
!
      data oarx,oary / 21.,  8.5/
      data tc       /  0.25 /
!
! *   taille et position des graduations
!
      data stipox,stipoy /'oi','oi'/
      data stisix,stisiy /0.30,0.20/
      data grasix,grasiy /0.30,0.30/
!
! *   couleurs de la legende (date, numero et label LPP)
!     en unite ih (1-22), is(1-10), ib(1-10) selon mire rogralib
!
      data icdat1, cdat2 , cdat3   / 3, 1.0 , 0.7/
      data icnum1, cnum2 , cnum3   /20, 1.0 , 1./
      data iccet1, ccet2 , ccet3   / 3, 0.5 , 0.5/
!
! *   files codes: pour le conv_all_ps_to_png.bat
!     Le code 1 est utilise par la wordlib et on utilise le code 7
!     pour la rogralib qui produit les .ps
!
      data icpng /2/
!
! *   valeurs par defaut de la resolution en ppi (pts per inch)
!     pour le grand png, de la taille de l'image finale (en pts),
!     (redefini par la lecture des donnees)
!
      data iresx, iresy /  300, 300/
      data isizx, isizy / 1440, 960/

!     Valeur de l'offset pour le recadrage des .ps avant conversion en image
!     (redefini par la lecture des donnees)

      data tranx, trany /2.45, 2.45/
!
!                     ---------------
!
      print*
      print*, 'lecture des donnees pour le programme co_magneto'
      print*, '************************************************'

      call re_co_magneto(iyear,imon,iday,ih,im,is,ex,ikpl, &
                         parmod,rep,sunpo)

      print*
      print*, 'lecture des donnees pour l''animation'
      print*, '*************************************'

      call rap_dtnbp(deltat,nbp)

!     Valeur de l'offset pour le recadrage des .ps avant conversion en image

      call rap_qualimage(iresx, iresy,isizx, isizy,tranx, trany)

      print*
      print*, 'lecture des donnees pour le programme plot_magneto'
      print*, '**************************************************'

      call re_plot_magneto(ifill,isolwi, &
                           irvs1,irvs2,icbs1,icbs2,irbs1,irbs2, &
                           icmp1,icmp2,irmp1,irmp2,ircl1,ircl2, &
                           irps1,irps2,cocli,cooli,codip,corot,coarr)


      print*, 'lecture du directory des fichiers d''orbite'
      print*, '*******************************************'

      read(*,100) bid
      print*, 'directory des fichiers d''orbite ?'
      print*, 'ex: /NFS/nas-cluster1/CLUSTER/ORBIT/DATA'
      read(*,'(a)') orbidir
      print*, trim(orbidir)

      print*, 'plot of field lines crossing each spacecraft ? (y/n)'
      print*, '****************************************************'

      read(*,*) kline
      print * , kline

!            ----------------------------
!
! *** calcul de l'instant de debut
!
      call cdoyear(iyear,imon,iday,idoty1)
      call chouday(ih,im,is,dech1)
!
      iyear1=iyear
!
! *** ouverture et entete du fichier de conversion ps to png
!
      call wheader_ps_to_png(icpng,'conv_all_ps_to_png.sh')

! *** lecture du fichier monde et calcul de l'image projetee

      close(1)
      call rewormap(1,icowor,3)
!
!
!     -------------------
! *** BOUCLE SUR LE TEMPS
!     -------------------
!
      do 10 itime=1,nbp
!
      print*
      print*, 'itime=',itime
      print*, '------'
      print*
!
! *** calcul du nom generique des fichiers
!
      call calch3(itime-1,ch3)
!
      fileps='resu_ps/amag_'//ch3//'.ps'
      filpng='resu_png/amag_'//ch3//'.png'
!
! *** remplissage du fichier de conversion conv_all_ps_to_png.bat
!
      call fill_ps_to_png(icpng,fileps,filpng,iresx,iresy, &
                                                     isizx,isizy)
!
! *** initialisation des parametres dependants du temps
!
      dech=dech1 + (itime-1)*deltat/60.
!
! *   correction de l'annee et du jour de l'annee si idoty
!     depasse le nb de jour de l'annee courante
!
      call corr_date(dech,iyear1,idoty1,iyear,idoty)
!
! *   calcul de la date et de l'heure courante
!
      call cdatdoy(idoty,iyear,imon,iday)
      call ctimhou(dech,ih,im,is)
!
! *   calcul du ikp pour la date et l'heure consideree
!
      if(ikpl.lt.0) then
                    call cokp(1,iyear,imon,iday,ih,rkp)
                    call convkptsy(ex,rkp,ikp)
                    if(ikp.lt.0) ikp=4 ! valeurs par defaut si pas de kp

                    else
                    ikp=ikpl
                    call convkptsy2(ex,ikp,rkp)

      endif

      call convrkp(rkp,ckp)
!
! *** ouverture du fichier graphique
!
      call dopegra_(7,fileps)
      call dpagfor_('p')
!
! *** attributs de la page et recadrage de la partie interressante
!
      call dpagori_(tranx,trany)
      call dfontyp_('h')
!
! *** calcul de la magnetosphere complete
!
      call co_magneto(iyear,imon,iday,ih,im,is,ex,ikp,parmod,rep, &
                    linjou,linnui,linnor,linsud, &
                    nlijou,nlinui,nlinor,nlisud, &
                    npljou,nplnui,nplnor,nplsud, &
                    magpau,npau,bowsho,nbow,     &
                    dxrep,dyrep,dzrep, rxrep,ryrep,rzrep)

! *** calcul de l'image projetee de la Terre

      call coworpic(icowor,imawor,idiwor,rep,sunpo)

! *** trace de la magnetosphere complete

      call plot_magneto(linjou,linnui,linnor,linsud, &
                        nlijou,nlinui,nlinor,nlisud, &
                        npljou,nplnui,nplnor,nplsud, &
                        magpau,npau,bowsho,nbow,     &
                        imawor,idiwor,               &
                        dxrep,dyrep,dzrep, rxrep,ryrep,rzrep,  &
                        ofx,ofy,sfx,sfy,x1,x2,y1,y2,oarx,oary, &
                        stipox,stipoy,stisil,stisis,grasix,grasiy,&
                        ifill,isolwi, &
                        irvs1,irvs2,icbs1,icbs2,irbs1,irbs2, &
                        icmp1,icmp2,irmp1,irmp2,ircl1,ircl2, &
                        irps1,irps2,cocli,cooli,codip,corot,coarr)
!
! *** calcul et trace du vent solaire
!
      gdx=1.5
      eps=0.5
      ts=0.05
      istep=5
      nprot=800
!
      call co_plot_flow(gdx-eps, 11.0,sunpo,y1,y2,ts,itime,istep,nprot)
      call co_plot_flow(gdx    ,  6.5,sunpo,y1,y2,ts,itime,istep,nprot)
      call co_plot_flow(gdx+eps,  2.0,sunpo,y1,y2,ts,itime,istep,nprot)
!
      call co_plot_flow(gdx+eps, -2.0,sunpo,y1,y2,ts,itime,istep,nprot)
      call co_plot_flow(gdx    , -6.5,sunpo,y1,y2,ts,itime,istep,nprot)
      call co_plot_flow(gdx-eps,-11.0,sunpo,y1,y2,ts,itime,istep,nprot)

! *** calcul et plot des 4 lignes de force passant par les sat'

      if(kline.eq.1) then
              call co_clu_pos(orbidir,iyear,imon,iday,ih,im,is,rep)
              call coplot_linesat(ex,ikp,parmod,rep)
                     endif
!
! *** calcul et trace l'orbite moyenne et la position de 4 sat Cluster
!
      call co_plot_cluorbit(orbidir,iyear,imon,iday,ih,im,is,sunpo,rep)
!
! *** retrace du cadre
!
      call dlincol_('n')
      call dlinwid_(3)
      call pfigfra_
!
! *** impression sur la figure de la date courante, de l'heure
!     du # de frame, du label LPP
!
      call upper_case(rep,3)
      call upper_case(ex,1)
      call cmoncar_(imon,cmon)
      call ctimcha_(ih,im,is,ctime)
      call cchalen_(cmon,ncmon)

!   * date/time

      write(com,200) iyear, trim(cmon),iday,ctime(1:5)
  200 format(i4,1x,a,i3.2,', 'a)

      call dlinihsb(icdat1,cdat2,cdat3)
      call dfontyp_('hb')
      call pfigcha_(22.0, 14.5 ,-1, 0.33,0.33,0.,com)

!   * Tsy, Kp, rep

      write(com,300) trim(ex),ckp,rep
  300 format(a,'  Kp ',a2,'  ',a)
      call dfontyp_('h')
      call pfigcha_(22.0, 11.5 ,-1, 0.30,0.30,0.,com)

!   * numero de frame

      call dfontyp_('hb')
      call dlinihsb(icnum1,cnum2,cnum3)
      call pfigcha_(41.0, 11.5,-1, 0.30,0.30,0.,ch3)

!    * label LPP

      call dlinihsb(iccet1,ccet2,ccet3)
      call dfontyp_('hbobl')
      call pfigcha_( 36.0,-13.7 ,-1, 0.60,0.60,0.,'LPP')
      call dfontyp_('h')

      call lower_case(rep,3)
      call lower_case(ex ,6)
!
! *** fermeture du fichier graphique
!
      call dclogra_
!
   10 continue
!
! *** fermeture des fichiers .sh
!
      close(icpng)
!
  100 format(80a)
!
      stop 'OK coanim_magsphe_cluorbit'
      end
!
!     /////////////////////////////////////////////////////////////////
!
      subroutine co_plot_cluorbit(orbidir,iyear,imon,iday,ih,im,is,sunpo,rep)
!
! *** calcul et trace l'orbite moyenne et la position de 4 sat Cluster
!
      integer ifc(4)
      real   satpos(3,4),satvit(3,4)
      real   satrev(  4)
      real posgx(1000),posgz(1000)
      real posx(4),posy(4),posz(4)
!
      character*(*) orbidir
      character*(*) sunpo,rep
      character*1 coul(4)
!
      data ifc /11,12,13,14/
!
      save posgx,posgz, icall
!
!
! *** ouverture des fichiers d'orbite
!
      call openorbi(orbidir,iyear,imon,iday,ifc)
!
! *** calcul des positions de cluster dans le GEI pour l'instant donne
!
      call cluposvit(iyear,imon,iday,ih,im,is,ifc,satpos,satvit,satrev)
!
! *   passage dans le gse ou le gsm
!
      call ctimpar(iyear,imon,iday,ih,im,is)
!
      if(rep(1:3).eq.'GSE'.or.rep(1:3).eq.'gse') call gei_to_gse(satpos)
      if(rep(1:3).eq.'GSM'.or.rep(1:3).eq.'gsm') call gei_to_gsm(satpos)
!
! *** passage en RT
!
      do j=1,4
      posx(j)=satpos(1,j)/6370.
      posy(j)=satpos(2,j)/6370.
      posz(j)=satpos(3,j)/6370.
      enddo
!
! *** trace des 4 satellites dasn le plan x-z
!
      sx=0.2
      sy=0.2
!
      if(sunpo(1:1).eq.'l') then
                            do j=1,4
                            posx(j)= -posx(j)
                            enddo
                            isens= -1
                            else
                            isens=  1
                            endif
!
! *   ordre du trace
!
      coul(1)='n'
      coul(2)='r'
      coul(3)='g'
      coul(4)='b'
!
      print*
      do j=1,4
      print*, 'sat ',j,coul(j),' xyz= ',posx(j),posy(j),posz(j)
      enddo
!
      call sort_arr_xyzc(posy,posx,posz,coul,4,isens)
!
      print*
      print*, 'ordre du trace: ',coul
!
! *   trace des cercles pleins et des contours en noir
!
      call dlinwid_(1)
!
      call dlincol_(coul(1))
      call pfigsym_(posx(1),posz(1),0,sx,sy,'dode')
      call pfilzon_
      call dlincol_('n')
      call pfigsym_(posx(1),posz(1),0,sx,sy,'dode')
!
      call dlincol_(coul(2))
      call pfigsym_(posx(2),posz(2),0,sx,sy,'dode')
      call pfilzon_
      call dlincol_('n')
      call pfigsym_(posx(2),posz(2),0,sx,sy,'dode')
!
      call dlincol_(coul(3))
      call pfigsym_(posx(3),posz(3),0,sx,sy,'dode')
      call pfilzon_
      call dlincol_('n')
      call pfigsym_(posx(3),posz(3),0,sx,sy,'dode')
!
      call dlincol_(coul(4))
      call pfigsym_(posx(4),posz(4),0,sx,sy,'dode')
      call pfilzon_
      call dlincol_('n')
      call pfigsym_(posx(4),posz(4),0,sx,sy,'dode')
!
! *** position du centre de gravite 
!
      pxg= (posx(1) + posx(2) +posx(3) +posx(4))/4.
      pyg= (posy(1) + posy(2) +posy(3) +posy(4))/4.
      pzg= (posz(1) + posz(2) +posz(3) +posz(4))/4.
!
      print*
      print*, 'position du centre de gravite de Cluster: ', &
     &          pxg,pyg,pzg, ' rep=',rep, ' sunpo=',sunpo
!
! *** trace de l'orbite
!
      icall =icall +1
      if(icall.gt.1000) stop '*** coanim_magsphe_cluorbit: plus de 1000 pos.'
!
      posgx(icall)=pxg
      posgz(icall)=pzg
!
      call dlinihs_(2,9)
      call dlinwid_(4)
      call pfigcur_(posgx,posgz,icall)
      call dlincol_('n')
!
      return
      end
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

       subroutine co_clu_pos(orbidir,iyear,imon,iday,ih,im,is,rep)
 
! *** calcul l'orbite moyenne et la position de 4 sat Cluster
 
      integer ifc(4)

      real(kind=4), dimension(3,4) :: satpos, satvit
      real(kind=4), dimension(4)   :: satrev

      real(kind=4), dimension(4)  :: sposx,sposy,sposz
      
      common /sat_spos/ sposx,sposy,sposz, spgx,spgy,spgz
 
      character*(*) orbidir
      character*(*) rep
 
      data ifc /11,12,13,14/
 
! *** ouverture des fichiers d'orbite
 
      call openorbi(orbidir,iyear,imon,iday,ifc)
 
! *** calcul des positions de cluster dans le GEI pour l'instant donne
 
      call cluposvit(iyear,imon,iday,ih,im,is,ifc,satpos,satvit,satrev)
 
! *   passage dans le gse ou le gsm
 
      call ctimpar(iyear,imon,iday,ih,im,is)
 
      if(rep(1:3).eq.'GSE'.or.rep(1:3).eq.'gse') call gei_to_gse(satpos)
      if(rep(1:3).eq.'GSM'.or.rep(1:3).eq.'gsm') call gei_to_gsm(satpos)
 
! *** passage en RT
 
      do j=1,4
      sposx(j)=satpos(1,j)/6370.
      sposy(j)=satpos(2,j)/6370.
      sposz(j)=satpos(3,j)/6370.
      enddo
 
! *** position du centre de gravite
 
      spgx= (sposx(1) + sposx(2) +sposx(3) +sposx(4))/4.
      spgy= (sposy(1) + sposy(2) +sposy(3) +sposy(4))/4.
      spgz= (sposz(1) + sposz(2) +sposz(3) +sposz(4))/4.

      call closeorbi(ifc)
 
      return
      end
 
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
 
