c $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/soltht.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
        subroutine soltht(z1,z2,re,Hrf1,n01,c,tht,xp,coem,st3,ds,dsN)
c
c calc the refraction integral (angular distance traveled)
c input:
c z1 - first altitude
c z2 - final altitude
c re - radius of earth
c Hrf - refractive scale height
c n01 - constant for index of refraction
c c - optical constant of the motion
c st3 - sine of local zenith angle at point3
c
c output:
c tht - angle traveled (radians)
c xp - reduced corridinate (can be used in other refraction integrals
c coem - coes for quad. fit to xp
c ds - geometric path length
c dsN - index of refraction weighted geometric path length
c
        implicit none
        include 'consts.cmn'
c the rcs revision #
      character*20 rcsrev
      parameter (rcsrev='$Revision: 2.22 $')
c passed
        real*8 z1,z2,re,Hrf1,n01,c,tht,xp(3),coem(3,3),st3
c local
        integer i,j
        real*8 zmid,np(3),rp(3),fpt(3),cc,coe(3),ds,dsN,fpd(3),fpm(3)
        real*8 rfintp,x1,x2,x3,Nm1(3),Nbar,st1,p1
        external rfintp
        logical debug
        data debug /.FALSE./
        save debug
c
c check to see that z1 and z2 differ by amount larger than smallx
      if((z2-z1).gt.smallx)then
        cc=c*c
c
        rp(1)=re+z1
        zmid=0.7*z1+0.3*z2
        rp(2)=re+zmid
        rp(3)=re+z2
c
c
        Nm1(1)=rfintp(Hrf1,n01,z1)
        Nm1(2)=rfintp(Hrf1,n01,zmid)
        Nm1(3)=rfintp(Hrf1,n01,z2)
        np(1)=1.0d0+Nm1(1)
        np(2)=1.0d0+Nm1(2)
        np(3)=1.0d0+Nm1(3)
c
c the sine of the local zenith angle at point 3
c
        st3=c/rp(3)/np(3)
c
        do i=1,3
           xp(i)=(rp(i)*np(i))**2 -cc
           if(xp(i).gt.0.0)then
             xp(i)=sqrt(xp(i))
           else
             xp(i)=0.0d0
           endif
        enddo
        x1=(xp(3)-xp(1))
        x2=(xp(3)**2-xp(1)**2)/2.0d0
        x3=(xp(3)**3-xp(1)**3)/3.0d0
        if(debug)write(6,*)'soltht. Hrf',Hrf1
        if(debug)write(6,*)'soltht. rp np xp',rp,np,xp
c get coes. for quad fit
        call quadcoe(xp,coem)
        if(debug)write(6,*)'soltht. coem',coem
        if(Hrf1.gt.-smallx)then
c scale height is >0
          do i=1,3
             fpt(i)=1.d0/(np(i)*rp(i)*np(i)*rp(i))
             fpd(i)=1.d0/(np(i))
             fpm(i)=Nm1(i)*fpd(i)
          enddo
          write(6,*)'soltht. c ',c
          write(6,*)'soltht. xp ',xp
        else
          do i=1,3
             fpt(i)=1.d0/(1.d0-rp(i)*(1.d0-np(i))/Hrf1)/
     &                 (np(i)*rp(i)*np(i)*rp(i))
             fpd(i)=1.d0/(np(i)-rp(i)*(1.d0-np(i))/Hrf1)
             fpm(i)=Nm1(i)*fpd(i)
          enddo
        endif
c angular distance integral
        if(debug)write(6,*)'soltht. fpt',fpt
        do i=1,3
           coe(i)=0.0d0
           do j=1,3
              coe(i)=coe(i)+fpt(j)*coem(j,i)
           enddo
        enddo
        if(debug)write(6,*)'soltht. coe',coe
        tht=(coe(1)*x1+coe(2)*x2+coe(3)*x3)*c
c geometeric distance
        do i=1,3
           coe(i)=0.0d0
           do j=1,3
              coe(i)=coe(i)+fpd(j)*coem(j,i)
           enddo
        enddo
        ds=(coe(1)*x1+coe(2)*x2+coe(3)*x3)
        st1=c/rp(1)/np(1)
        p1=rp(1)*st1
c use alternate distance formula
        ds=sqrt(rp(3)*rp(3)-p1*p1)-sqrt(rp(1)*rp(1)-p1*p1)
c index of refraction weighted geometeric path length
        do i=1,3
           coe(i)=0.0d0
           do j=1,3
              coe(i)=coe(i)+fpm(j)*coem(j,i)
           enddo
        enddo
        nbar=(Nm1(1)+Nm1(3))*0.5d0
        if(nbar.gt.0.0)then
          dsN=(coe(1)*x1+coe(2)*x2+coe(3)*x3)/nbar
        else
          dsN=0.0d0
        endif
c        if(nbar.gt.smallx)then
c          dsN=(coe(1)*x1+coe(2)*x2+coe(3)*x3)/nbar
c        else
c          dsN=0.0d0
c        endif
      else
        tht=0.0d0
        ds=0.0d0
        dsN=0.0d0
        do i=1,3
           xp(i)=0.0d0
           do j=1,3
              coem(j,i)=0.0d0
           enddo
        enddo
        rp(1)=re+z1
c
        np(1)=1.0d0+rfintp(Hrf1,n01,z1)
c
c the sine of the local zenith angle at point 3 (which is the same as pt1).
c
        st3=c/rp(1)/np(1)
      endif
        if(debug)write(6,*)'soltht. tht ',tht*180.0/acos(-1.0d0)
c        if(debug)stop
        return
        end