c This program estimates the distance to a pulsar given the
c galactic longitude, latitude and dispersion measure, using
c the exp(-x) model.
c On output, a flag=1 means that distance is a lower bound.

      program exponen
      implicit none

      integer i,flag(3)
      real*8 par(6),l,b,DM,D(3),l0,b0
      real*8 dens,dispmeas

      common /param/ par
      common /direc/ l,b

c----> Parameters of the model.
      data par / 0.0203, 0.0071, 1.0684, 0.0533, 30.4410, 1.4851 /

c----> Reads position and DM
      write(*,*) 'l ?'
      read(*,*) l0
      l=l0*3.141592/180.
      write(*,*) 'b ?'
      read(*,*) b0
      b=b0*3.141592/180.
      write(*,*) 'DM [cm^-3 pc] ?'
      read(*,*) DM

c----> Calculates distance and prints result.
      call dist(DM,D,flag)
      write(*,1000) 'l','b','D-','D','D+','DM','flags'
      write(*,1001) l0,b0,(D(i),i=1,3),DM,(flag(i),i=1,3)
 1000 format(7a10)
 1001 format(6f10.2,4x,3i2)

      end

c----------------------------------------------------------------
c Returns the density at a given position.

      real*8 function dens(D)
      implicit none

      common /param/ n0,n1,z0,z1,r0,r1
      common /direc/ l,b

      real*8 D
      real*8 r,z,Rsol,rg,c0,c1
      real*8 l,b,n0,n1,z0,z1,r0,r1

      parameter(Rsol=8.5)

c----> Calculates the distance from the center of the galaxy (r)
c      and the distance from the plane.
      z=abs(D*sin(b))
      rg=D*cos(b)
      r=sqrt(Rsol**2 + rg**2 - 2*Rsol*rg*cos(l))

      c0=0.
      c1=0.

c----> If it is farther than 20 scale lengths in either direction,
c      the density of that component is zero.
      if (z/z0 .le. 20. .and. r/r0 .le. 20.)
     1  c0=n0*exp(-z/z0-r/r0+Rsol/r0)
      if (z/z1 .le. 20. .and. r/r1 .le. 20.)
     1  c1=n1*exp(-z/z1-r/r1+Rsol/r1)

      dens=c0+c1

      end

c----------------------------------------------------------------
c Returns the dispersion measure between the distances D0 and D1
c from the Sun. It uses the Simpson's rule, with a step size
c equal to the smallest of a fifth of the scale lengths or a
c fifth of the distance.

      real*8 function dispmeas(D0,D1)
      implicit none

      common /param/ n0,n1,z0,z1,r0,r1

      real*8 D0,D1
      real*8 n0,n1,z0,z1,r0,r1
      real*8 dens

      real*8 h,extreme,odd,even
      integer i,n

      h=min(z1,r1)/5.

c----> Guesses the number of steps.
      n=1+dabs(D1-D0)/h

      if (n .lt. 4) n=4
c----> Force n to be even.
      if ((n/2)*2 .ne. n) n=n+1

c----> Picks the step size.
      h=(D1-D0)/n

      extreme=dens(D0)+dens(D1)
      odd=0.
      do 10 i=1,n-1,2
        odd=odd+dens(D0+i*h)
 10   continue
      even=0.
      do 11 i=2,n-1,2
        even=even+dens(D0+i*h)
 11   continue

      dispmeas=1e3*(h/3.)*(extreme+4.*odd+2.*even)

      end

c-------------
      subroutine dist(DM,D,flag)
      implicit none

      integer i,flag(3)
      real*8 DM,DMaux,dispmeas,D(3)
      real*8 D0,D1,D2,DM0,DM1,DM2
      real*8 eps,zero
      parameter(eps=1e-4)

      do 40 i=1,3
        if (i .eq. 1) then
          DMaux=DM*0.7
        else if (i .eq. 2) then
          DMaux=DM
        else
          DMaux=DM*1.3
        endif

        D0=0.
        D1=1d4
        DM0=0.
        DM1=dispmeas(D0,D1)

        if (DM1 .lt. DMaux) then
          flag(i)=1
          zero=0.
 20       continue
            D2=D0+D1/2.
            DM2=dispmeas(zero,D2)
            if (D1 .lt. eps) then
              D(i)=D0
            elseif ((DM1-DM2) .le. eps) then
              D1=D2-D0
              goto 20
            else
              D1=D2-D0
              D0=D2
              goto 20
            endif
        else
 30       continue
            D(i)=(D1+D0)/2.
            DM2=DM0+dispmeas(D0,D(i))
            if (dabs(DM2-DMaux)/DMaux .le. eps) then
              goto 40
            elseif (DM2 .gt. DMaux) then
              D1=D(i)
              DM1=DM2
            else
              D0=D(i)
              DM0=DM2
            end if
          goto 30
        endif
 40   continue

      end