subroutine fminth(rfndx,z,z1,a1,re,nz,zmin,debug)
c
c
c find the minimum tangent height for a ray given the inital altitude,
c and local zenith angle of the ray and the end altitude of the ray
c
c input:
c rfndx - array of refractive index as function of radial distance from
c         center of the earth
c z     - altitude above spherical earth
c z1    - inital altitude distance
c a1    - inital local zenith angle of the ray propagation direction [rad]
c re    - radius of earth
c
c output:
c zmin
c 
c steps:
c 1. compute the optical constant of the motion (c=n(r1)*r1*sin(a1))
c 2. check to see if c is less than n(rsfc)*rsfc
c 3. if does not strike the earth, create array of n(r)*r values
c 4. using c interpolate r as function of n(r)*r to find r_tangent 
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      implicit none
c the rcs revision #
      character*20 rcsrev
      parameter (rcsrev='$Revision: 2.22 $')
c
c passed
      integer nz
      real*8 rfndx(nz),z(nz),z1,a1,re,zmin
      logical debug
c local
      integer klo,khi,i
      real*8 repz1,fi,fip1,zi,zip1,zip2,nz1,Hrf,n0,c,csfc
      external rfintp
      real*8 rfintp
c      data debug/.FALSE./
c      save debug
c
c get n(r1)
      repz1=re+z1
c
c find the 2 pts that straddle z1
c
      call findex(z,nz,z1,klo,khi)
c
c calc the constant of the motion
c
      if(z1.ge.z(1))then
        nz1=0.0d0
      else
        call rfscal(rfndx(klo),rfndx(khi),z(klo),z(khi),Hrf,n0)
        nz1=rfintp(Hrf,n0,z1)
      endif
      if(debug)write(6,*)'fminth. z ',z(klo),z1,z(khi),nz1,z(1)
      if(debug)write(6,*)'fminth. rfndx',rfndx(klo),rfndx(khi),rfndx(nz)
      c=(1.0d0+nz1)*repz1*sin(a1)
c
c check to see if it is less than the value for a tangent beam at the sfc
c
      csfc=(1.0d0+rfndx(nz))*re
      if(debug)write(6,*)'fminth. c ',c,csfc
      if(c.lt.csfc)then
        zmin=0.0d0
      else
c
c search for the minimum altitude starting from the top of atmo
c
        khi=2
        do while((c.lt.(1.0d0+rfndx(khi))*(re+z(khi))).and.(khi.lt.nz))
           khi=khi+1
        enddo
        klo=khi-1
        call rfscal(rfndx(klo),rfndx(khi),z(klo),z(khi),Hrf,n0)
c
c solve for zero using Netwon method
c
        fi=(re+z(klo))*(1.0d0+rfndx(klo))
        fip1=(re+z(khi))*(1.0d0+rfndx(khi))
c
        zi=z(klo)
        zip1=z(khi)
        do i=1,3
           if(fi.ne.fip1)then
             zip2=zip1+(c-fip1)/((fip1-fi)/(zip1-zi))
           else
             zip2=zip1
             goto 10
           endif
           if(i.eq.1)then
             if(abs(zip2-zip1).lt.abs(zip2-zi))then
               zi=zip1
               fi=fip1
             endif
           else
             zi=zip1
             fi=fip1
           endif
           zip1=zip2
           fip1=(re+zip1)*(1.0d0+rfintp(Hrf,n0,zip1))
        enddo
10      zmin=zip2
      endif
      end
c
      subroutine rfscal(f1,f2,z1,z2,Hrf,n0)
c
c compute the scale height and constant given the 2 pts and the altitude of the pts.
c
c input:
c f1 - function at z1
c f2 - function at z2
c z1 - altitude z1
c z2 - altitude
c 
c output:
c Hrf - scale height
c n0 - constant
c
c steps:
c 1. assume that f(z)=n0*exp(z/hrf)
c 2. solve for Hrf=(z2-z1)/log(f1/f2) if (f1 ne f2)
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
c
c passed
      real*8 f1,f2,z1,z2,Hrf,n0
c local
c
      if((f1.ne.f2).and.(f1.ne.0.0).and.(f2.ne.0.0))then
        Hrf=(z2-z1)/log(f2/f1)
        n0=f1*exp(-z1/Hrf)
      else
        Hrf=0.0d0
        n0=f1
      endif
      return
      end
c
c
      real*8 function rfintp(Hrf,n0,z)
c
c compute f(z)=n0*exp(z/Hrf)
c
c input:
c Hrf - scale height
c n0 - constant
c z  - altitude
c
c output:
c rfintp - f(z)
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      include 'consts.cmn'
c passed
      real*8 Hrf,n0,z
c local
      real*8 expon
c
      if(Hrf.gt.-smallx)then
        rfintp=n0
      else
c
        expon=z/Hrf
        if(z/Hrf.lt.-300.0)then
          rfintp=0.0d0
        else if(z/Hrf.gt.300.0)then
c        write(6,*)'rfintp. exponent overflow. Will set to e(300)'
          rfintp=n0*exp(300.0)
        else
          rfintp=n0*exp(expon)
        endif
      endif
      return
      end