! -----------------------------------------------------------------------------
! Frame   : Procedures and program used in JGR data report 
!           "Use of twenty years CLUSTER/FGM data to observe the mean behavior 
!           of the magnetic field and current density of Earth's magnetosphere" 
!           By P. Robert and M. Dunlop, October 2021
!
! Library : mag_bow.f90
! Object  : routine for comuting bow shock and MP model
! Author  : P. Robert, LPP-ScientiDev, 2020-Oct.2021
! Mail    : patrick.robert@lpp.polytechnique.fr
! -----------------------------------------------------------------------------

! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

  subroutine comagpau_pr(romagx,romagy,nmagpau)
  
  real(kind=4) :: romagx(221),romagy(221)
  
  p=18.
  s=12.
  e=p/s -1.
  
  do i=-110,110
    t=float(i)*3.14159/180.
    r=p/(1.+e*cos(t))
    romagx(i+111)=r*cos(t)
    romagy(i+111)=r*sin(t)
  enddo
  
  nmagpau=221
  
  return
  end
  
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine cobowsho2(bowshox,bowshoz,nbow)

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

!     ******************************************************************
!     calcul du bowshock  en GSE ou GSM (idem car symetrie autour de x)
!     ******************************************************************

      x1=-40.
      x2= 13.5
      dx=  0.1
      nx= int((x2-x1+0.0001)/dx) +1
      nbow=2*nx

      if(nbow.gt.1072) then
                       print*, 'nbow=',nbow
                       print*, 'GT.1072 *** ERROR***'
                       stop 'error cobowsho'
                       endif

      do 10 ix=1,nx

      x=x1+float(ix-1)*dx
      if(ix.eq.1) x=-75.
      if(x .gt. x2) then
            nx=ix-1
            nbow=2*nx
            go to 11
      endif

      call co_zbowsho(x,z,ierr)

      if(ierr.ne.0) go to 10

      bowshox(ix)= x
      bowshoz(ix)= z

   10 continue
   11 continue


! *** remplissage pour la partie des z negatifs

      do 20 ix=nx+1,nbow
      bowshox(ix)= bowshox(nbow-ix+1)
      bowshoz(ix)=-bowshoz(nbow-ix+1)
   20 continue

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine co_zbowsho(x,z,ierr)

!     ******************************************************************
!     calcul du bowshock  en GSE ou GSM (idem car symetrie autour de x)
!     d'apres CLUSTER: Mission payload and supporting activities,
!     ESA SP-1159, page 259
!@INPROCEEDINGS{1993ESASP1159..259R,
!       author = {{Rodriguez-Canabal}, J. and {Warhaut}, M. and {Schmidt}, R. and {Bello-Mora}, M.},
!        title = "{The Cluster orbit and mission scenario}",
!    booktitle = {Cluster: Mission, Payload and Supporting Activities},
!         year = 1993,
!       editor = {{Burke}, W.~R.},
!       series = {ESA Special Publication},
!       volume = {1159},
!        month = mar,
!        pages = {259},
!       adsurl = {https://ui.adsabs.harvard.edu/abs/1993ESASP1159..259R},
!      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
!}
!     ******************************************************************

      c=  56.7
      d= 177.7
      s=  13.5

      p1= (d+s-x)**2
      p2= d*d
      p3=p1/p2-1.

      r2=p3*c*c + x*x
      ierr=0

      if(r2.lt.0.) then
                   ierr=1
                   return
                   endif

      z2=p3*c*c

      if(z2.lt.0.) then
                   ierr=2
                   return
                   endif

      z=sqrt(z2)

      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

