
  program plot_field_lines_XZ
      
  use dim_grid 
  
  character(len=255) :: psfile,flpath,flname,bid
  character(len=1) :: pm,BJ
  character(len=3) :: ckp
  character(len=5) :: valdip

  real (kind=4), dimension(1072) :: bowshox,bowshoz
  real(kind=4) :: romagx(221),romagy(221)

  DIMENSION XXT(2000,5000),YYT(2000,5000),ZZT(2000,5000), LT(5000), iplot(5000)
  DIMENSION BBX(2000,5000),BBY(2000,5000),BBZ(2000,5000), B0(2000,5000)
  
  nlmax=5000
  iplot=1


  print*,'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx'
  print*
  print*, 'FL file name? ex: Field_Lines/FL_BXZ_160160160_p00_AKp.dat'
 read (*,'(a)') flpath
  print*, flpath

  i=index(flpath,'/',back=.true.)
  if(i>0) then
      flname=flpath(i+1:)
          else
      flname=flpath
  endif
  
  print*, 'flname=',trim(flname)
  read(flname,'(3x,a1,13x,a1,i2,1x,a3)') BJ,pm,idta,ckp
  dta=float(idta)
  if(pm == 'm') dta=-dta

  print*, '------------------------------------------------------------------'
  print*,'flname= ',trim(flname),'  dta=',dta,' kp=',ckp
  print*, '------------------------------------------------------------------'

  i=len(trim(flpath))
  psfile=flpath(1:i-4)//'.ps'
  write(valdip,'(f5.1)') dta
  print*,'valdip=',valdip
  
! ============================================================== 

! figure, cadre, legende

  x1=-20.
  x2= 20.
  y1=-20. 
  y2= 20.
  z1=-20.
  z2= 20.

  orix=1.7
  oriy=4.5
  sizx=18.
  sizy=18.

  tc=0.4

  call dopegra_(17,psfile)
  call dfigsiz_(sizx,sizy)
  call dfigori_(orix,oriy)
  call dfiglim_(x1,x2,y1,y2)
  call dlabsiz_(0.4,0.4)
  call dstipos_('oi','oi')

  call ppagcha_(orix+0.75*sizx/2.,oriy+sizy+1.,0,tc,tc,0.,psfile)
 
! ----------------------------------------

  print*,'magnetopause PR'
  
! remplissage magnetosphere

  call comagpau_pr(romagx,romagy,nmagpau)
  call dlinrgb_(0.9,0.9,0.9)
  call dfilzon_
  call pfigcur_(romagx,romagy,nmagpau)
  call pfigpmo_(x1,y2)
  call pfigpmo_(x1,y1)
  call pfigpmo_(romagx(1),y1)
  call pfilzon_

! remplissage vent solaire

  call dlinrgb_(1.,1.,0.9)
  call dfilzon_
  call pfigcur_(romagx,romagy,nmagpau)
  call pfigpmo_(x2,y2)
  call pfigpmo_(x2,y1)
  call pfigpmo_(romagx(1),y1)
  call pfilzon_
  
! plot  magnetopause

  call dlinhsb_(4./36.,1.,1.)
  call dlinwid_(9)
  call pfigcur_(romagx,romagy,nmagpau)
  call dlincol_('n')
  
  print*,'bow shock'

  call cobowsho2(bowshox,bowshoz,nbow)
  where(bowshoz < y1)
     bowshoz=y1
     bowshox=2.9
  endwhere
  
  where(bowshoz > y2)
     bowshoz=y2
     bowshox=2.9
  endwhere
  
  call dlincol_('g')
  call dlinwid_(9)
  call pfigcur_(bowshox,bowshoz,nbow)
  
  call dlincol_('n')
  call dlinwid_(3)
  call pfiggrax(0.,5.,1.,'(f4.0)')
  call pfiggray(0.,5.,1.,'(f4.0)')
  
  call pfiglabx('X GSM')
  call pfiglaby('Z GSM')
  call pfigfra_
  
  bid=' theta='//valdip
  print*,trim(bid)
  call pfigcha_(-19.,17.,-1,0.5,0.5,0.,trim(bid))
  
!-----------------------------------------------------
  
! fleches dipole et x du SM

  call dlincol_('b')
  bodyl=11.
  bodyw=0.3
  headl=1.5
  headw=1.5
  call dfilzon_
  call pfigarr_(0.,0.,-1,bodyl,bodyw,headl,headw,90.-dta)
  call pfilzon_

  call dlincol_('n')
  call dlinwid_(1)
  call pfigarr_(0.,0.,-1,bodyl,bodyw,headl,headw,90-dta)

  bodyl=11.
  call dlincol_('r')
  call dfilzon_
  call pfigarr_(0.,0.,-1,bodyl,bodyw,headl,headw,-dta)
  call pfilzon_

  call dlincol_('n')
  call dlinwid_(1)
  call pfigarr_(0.,0.,-1,bodyl,bodyw,headl,headw,-dta)  
  
  print*, '---------'
  print*, 'read field lines', dta
  print*, '---------'

  nl=0
  
  open(1,file=flpath)
  
  10 continue
  nl=nl+1
  if(nl > nlmax) stop '***pb'
  read(1,*,end=20) LT(nl)
  
  do l=1,LT(nl)
     read(1,*) XXT(l,nl),YYT(l,nl),ZZT(l,nl),BBX(l,nl),BBY(l,nl),BBZ(l,nl)
     B0(l,nl)=sqrt(BBX(l,nl)**2 +BBY(l,nl)**2 +BBZ(l,nl)**2)
  end do
  go to 10
  
  20 continue
  nl=nl-1

  
  print*, 'nb lus    =',nl
  if(nl > nlmax) stop '***pb'
  
  print*, '---------'
  print*, 'select field lines'
  print*, '---------'        

  seuil=0.6
  nbeli=0
  nbplo=0
  
  do n=1,nl
     iplot(n)=1
     call sablier(nl,n)
     
     do nn=1,n-1
     if(iplot(nn) == 1) then
        call compare_line(XXT(1:LT(nn),nn),ZZT(1:LT(nn),nn),LT(nn),&
                          XXT(1:LT(n),n),  ZZT(1:LT(n),n),LT(n), &
                          seuil,idiff)
        if(idiff == 2) cycle
        if(idiff == 0) then
            iplot(n)=0
            nbeli=nbeli+1
            exit
        endif
     endif
     end do
  end do
  
  print*, 'nb elimine=',nbeli
  
  print*, '---------'
  print*, 'plot field lines'
  print*, '---------'        

! n: numero de la ligne contenant LT(n) points  

  call dlinwid_(1)
  
! correction d'amplitude si J
  if(BJ == 'J') B0=B0*50.
  
  do n=1,nl
      call sablier(nl,n)
     if(iplot(n) == 1) then
        nbplo=nbplo+1
        call plot_line(XXT(1:LT(n),n),ZZT(1:LT(n),n),B0(1:LT(n),n),LT(n))
     endif
     
  end do

  print*
  print*, 'nb lus    =',nl
  print*, 'nb elimine=',nbeli
  print*, 'nb plot   =',nbplo
 

! plot de la Terre

  call dlincol_('c')
  call dlinwid_(3)
  call dfilzon_
  call pfigcir_(0.,0.,1.,90.,270.,5.)
  call pfilzon_

  call dlincol_('y')
  call dfilzon_
  call pfigcir_(0.,0.,1.,-90.,90.,5.)
  call pfilzon_

  call dlincol_('n')
  call pfigcir_(0.,0.,1.,0.,360.,5.)
  call pfiglin_(0.,-1.,0.,1.)
  call dlincol_('n')
  call pfigfra_
! ----------------------------------------

  call dclogra_

  stop
  end

!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

  subroutine plot_line(XX,ZZ,B0,k)

  DIMENSION XX(k),ZZ(k),B0(k)
  

  data huep /1./  
  
  call dfigppo_(XX(1),ZZ(1))
  
  do i=2,k
     call cocoul(B0(i),hue)
     if(hue /= huep) then
                      call dfigppo_(XX(i-1),ZZ(i-1))
                      call dlinhsb_(hue,1.,1.)
     endif
     call pfigpmo_(XX(i),ZZ(i)) 
     huep=hue
  end do
  
  return
  end
  
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

  subroutine cocoul(B,hue)

  
  hue=1.
  
  if(B < 500.) hue=34./36.
  if(B < 250.) hue=32./36.
  if(B < 125.) hue=30./36.
  if(B <  60.) hue=28./36.
  if(B <  30.) hue=26./36.
  if(B <  10.) hue=24./36.

! old values for J

!  if(AJ0 < 15.) hue=35./36.
!  if(AJ0 < 10.) hue=33./36.
!  if(AJ0 <  5.) hue=30./36.
!  if(AJ0 <  1.) hue=20./36.
!  if(AJ0 <  0.2) hue=18./36.
!  if(AJ0 <  0.1) hue=16./36.

  return
  end
  
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

   subroutine aera(x,y,nbp,dx,dy,surf)
   
   ! surface d'un polygone ;semble ne pas marcher... pas utile.
   
   real x(4000),y(4000),dx(4000),dy(4000)
   
   do i=1,nbp
      dx(i)= (x(MOD(nbp,i+1)+1) -x(MOD(i+nbp-1,nbp)+1))/2.
      dy(i)= (y(MOD(nbp,i+1)+1) -y(MOD(i+nbp-1,nbp)+1))/2.
   end do
   
   A=0.
   do i=1,nbp
      A = A + (X(i)*DY(i) - Y(i)*DX(i))
  end do
  
  surf=abs(A/2.)
  return
  end

!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

  subroutine surf_poly(x,y,nbp,surf)
  
  real x(4000),y(4000)
  
!  surface d'un polygone ; decomposition du polygone en triangle
! et somme des surfaces

! ne marche pas si polygone convexe

  xo=sum(x(1:nbp))/float(nbp)
  yo=sum(y(1:nbp))/float(nbp)
  
  surf=0.
  
  do i=1,nbp-1
     a=sqrt( (x(i+1)-x(i))**2 +(y(i+1)-y(i))**2)
     b=sqrt( (x(i+1)-xo  )**2 +(y(i+1)-yo  )**2)
     c=sqrt( (x(i  )-xo  )**2 +(y(i  )-yo  )**2)
     
     call surftri(a,b,c,st)  
     surf=surf+st 
  end do
   
  return
  end

!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

  subroutine surftri(a,b,c,surf)
   
! surface d'un triangle par formule de heron:

  p=(a+b+c)/2.
  qq=p*(p-a)*(p-b)*(p-c)

  if(qq < 0.) then
     surf=0.
     return
  endif
  
  surf=sqrt(qq)
  
  return
  end
  
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

  subroutine compare_line(x1,z1,n1,x2,z2,n2,seuil,idiff)
  
  real x1(2000),z1(2000),x2(2000),z2(2000)
! real p1(2000),p2(2000),q1(2000),q2(2000)
 ! real x(4000),y(4000)
  
  if(n1 < 3 .or. n2 < 3) then
              idiff=2
              return
  endif
  
! ----------------------------------------------  
! ! test par surface entre les deux lignes
! ! marche pas bien; pas utile

! ! fabrication du polynome

  ! do i=1,n1
    ! x(i)=x1(i)
    ! y(i)=z1(i)
  ! end do
  
  ! do i=1,n2
     ! ii=i+n1
     ! kk=n2-i+1
     ! x(ii)=x2(kk)
     ! y(ii)=z2(kk)
  ! end do
  
  ! nbp=ii+1
  ! x(nbp)=x1(1)
  ! y(nbp)=z1(1)
  
! ! calcul de la surface

  ! call surf_poly(x,y,nbp,surf)
  
! ! perimetre 

  ! pp=0.
  ! do i=2,nbp
    ! a=sqrt( (x(i-1)-x(i))**2 +(y(i-1)-y(i))**2)
    ! pp=pp+a
  ! end do
  
  ! test=2.*surf/pp
  
! ---------------------------------------------- 
  
! test utile qui marche
  
  k1=2*n1/3
  k2=2*n2/3
  
  dist=sqrt((x1(k1)-x2(k2))**2 + (z1(k1)-z2(k2))**2)
  test=dist
  
  xmax=(x1(k1)+x2(k2))/2.
  radi=( sqrt(x1(k1)**2 +z1(k1)**2) + sqrt(x2(k2)**2 +z2(k2)**2) )/2.
  
  if(radi > 11.5 .and. xmax > 3.) then
       test=test*2.
  endif
! ---------------------------------------------- 

! ! test par surface pour ligne ~horizontales

  ! xmia=min(x1(n1),x2(n2))  
  ! xmib=min( x1(1),x2(1)) 
  ! xmin=min(xmia,xmib)
  ! xmax=max(xmia,xmib)
  
  ! k1=0
  ! k2=0
  
! ! tableau des x communs
  
  ! do i=1,n1
    ! if(x1(i) > xmin .and. x1(i) < xmax) then
       ! k1=k1+1
       ! p1(k1)=x1(i)
       ! q1(k1)=z1(i)
    ! endif
  ! enddo

  ! do i=1,n2
    ! if(x2(i) > xmin .and. x2(i) < xmax) then
       ! k2=k2+1
       ! p2(k2)=x2(i)
       ! q2(k2)=z2(i)
    ! endif
  ! enddo  
  
  ! ! if(k1 <2 .or. k2 <2) then
           ! ! idiff=2
           ! ! return
  ! ! endif
  
 ! ! surface
 
  ! s1=0.
  ! s2=0.
  
  ! do i=1,k1-1
    ! s1=s1+0.5*abs(q1(i+1)-q1(i))*abs(p1(i+1)-p1(i))
  ! enddo
  
  ! do i=1,k2-1
    ! s2=s2+0.5*abs(q2(i+1)-q2(i))*abs(p2(i+1)-p2(i))
  ! enddo  
  
 ! ! distance moyenne

 ! ! test=abs(s2-s1)/abs(xmax-xmin) 
  
  ! !print*,'n1,n2,k1,k2,xmin,xmax,s1,s2,test=',n1,n2,k1,k2,xmin,xmax,s1,s2,test

! ---------------------------------------------- 
  
! et on fait le test  
  if(test > seuil) then
           idiff=1
                else
           idiff=0
  endif
  
  return
  end

!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
      
      subroutine sablier(ntot,n)
 
      if(n == ntot/20)   print*, '  5 % n=',n
      if(n == ntot/10)   print*, ' 10 % n=',n
      if(n == 2*ntot/10) print*, ' 20 % n=',n
      if(n == 3*ntot/10) print*, ' 30 % n=',n
      if(n == 4*ntot/10) print*, ' 40 % n=',n
      if(n == 5*ntot/10) print*, ' 50 % n=',n
      if(n == 6*ntot/10) print*, ' 60 % n=',n
      if(n == 7*ntot/10) print*, ' 70 % n=',n
      if(n == 8*ntot/10) print*, ' 80 % n=',n
      if(n == 9*ntot/10) print*, ' 90 % n=',n
      if(n ==   ntot   ) print*, '100 % n=',n

      return
      end

!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX



