c $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/raytrace.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
	subroutine raytrace(z1,z2,a1,re,z,rfndx,nz,th2,a2,
     &losdist,rpl,mupl,nplpts,x,layindx,rdebug)
c trace the ray from altitude z1 to z2 in the direction t1 at z1
c given by the sine of t1 (st1).
c
c input:
c z1 - initial alt.
c z2 - final alt.
c a1 - direction at z1
c
c output:
c th2 - angular distance traveled [radians]
c a2 - direction at z2
c
c losdist - segment distance
c rpl  -  radial distance from the center of the earth of the pts along the path
c mupl -  cosine of the local angle at pts along the path
c nplpts - # of path points
c x    -  angle between local zenith at A and local zenith at pts along the path
c         (radians)
c layindx - index of the layer for each segment
c
c
      implicit none
      include 'consts.cmn'
c the rcs revision #
      character*20 rcsrev
      parameter (rcsrev='$Revision: 2.22 $')
c passed
      integer nz
      real*8 z1,z2,a1,re,rfndx(nz),th2,z(nz),a2
      real*8 losdist(2*nz),rpl(2*nz+1),mupl(2*nz+1)
      real*8 x(nz*2+1)
      integer nplpts,layindx(2*nz)
      logical rdebug
c local
      integer klo1,khi1,klo2,khi2,i,klomin,khimin,j,nplpts2
      real*8 repz1,c,cc,nz1,Hrf1,n01,Hrf2,n02,dtht,coem(3,3),xp(3)
      real*8 st1,st2,atmp,thtmp,zmin,Hrfmin,n0min,amin,tmp1,xmax
      logical up,l2space,debug
c external
      real*8 rfintp
      external rfintp
c
      data l2space/.TRUE./
      save l2space
      data debug/.FALSE./
      save debug
      if(a1.le.pi_2)then
c go in the upward direction
        if(debug)write(6,*)'raytrace. going up'
	call traceupr(z1,z2,a1,re,z,rfndx,nz,th2,a2,
     &losdist,rpl,mupl,nplpts,x,layindx)
      else
        if(debug)write(6,*)'raytrace. going down'
c go downward
c
c compute the sine of the inital direction
c
        st1=sin(a1)
c
c find the index of refraction at the 2 altitudes
c find the layer of first index of refraction
c
c get n(r1)
        repz1=re+z1
c
c find the 2 pts that straddle z1
c
        call findex(z,nz,z1,klo1,khi1)
        if(debug)write(6,*)'raytrace. klo1 khi z',klo1,khi1,
     &z(klo1),z1,z(khi1)
c
c calc the constant of the motion
c
        if(z1.ge.z(1))then
          nz1=1.0d0
          call rfscal(rfndx(klo1),rfndx(khi1),z(klo1),z(khi1),Hrf1,
     &n01)
        else
          call rfscal(rfndx(klo1),rfndx(khi1),z(klo1),z(khi1),Hrf1,
     &n01)
          nz1=1.0d0+rfintp(Hrf1,n01,z1)
        endif
c
c constant of the motion
        c=st1*nz1*repz1
c
c find the 2 pts that straddle z2
c
        call findex(z,nz,z2,klo2,khi2)
        if(debug)write(6,*)'raytrace. klo2 khi z',klo2,khi2,
     &z(klo2),z2,z(khi2)
        call rfscal(rfndx(klo2),rfndx(khi2),z(klo2),z(khi2),Hrf2,n02)
c To calc. the angle a2, use the formula
c (re+z1)*nz1*st1=(re+z2)*nz2*st2 to solve for st2 (st2=sine(a2))
        a2=(c/(re+z2)/(1.0d0+rfintp(Hrf2,n02,z2)))
        if(a2.lt.1.0d0)then
          a2=asin(a2)
        else
          a2=pi_2
        endif
c there are 2 possibilities: 1. pt2 is between pt1 and tangent
c 2. tangent is between pt1 and pt2.
c We will choose the solution based upon the switch l2space
c if l2space is .TRUE., then case 2 will be used.
c
        if(l2space)then
          if(debug)write(6,*)'raytrace. will go 2 space'
c
c will go tangent
c
c first from the point z1 to the tangent point
c find the tangent height for this direction
          call fminth(rfndx,z,z1,a1,re,nz,zmin,debug)
          if(debug)write(6,*)'raytrace. z1 zmin',z1,zmin,a1/dr
c check to see what the angle is at altitude zmin
c
c find the 2 pts that straddle zmin
c
          call findex(z,nz,zmin,klomin,khimin)
          if(debug)write(6,*)'raytrace. klomin khi z',klomin,khimin,
     &z(klomin),zmin,z(khimin)
          call rfscal(rfndx(klomin),rfndx(khimin),z(klomin),z(khimin),
     &Hrfmin,n0min)
          amin=(c/(re+zmin)/(1.0d0+rfintp(Hrfmin,n0min,zmin)))
          if(debug)write(6,*)c,
     &(re+zmin)*(1.0d0+rfintp(Hrfmin,n0min,zmin))
          if((zmin-z(nz)).le.smallx)then
            if(amin.lt.1.0)then
              amin=asin(amin)
            else
              amin=pi_2
            endif
c it is not tangent, but strikes the earth.  So, will go from the ground up to z1.
	    call traceupr(zmin,z1,amin,re,z,rfndx,nz,th2,atmp,
     &losdist,rpl,mupl,nplpts,x,layindx)
c must reverse the arrays since we expect to go from z1 to z2 or zmin
            amin=pi-amin
            atmp=pi-atmp
            if(debug)
     &write(6,*)'this beam strikes the earth with a local ',
     &'zenith angle of ',real(amin/dr),real((a1-atmp)/dr)
            a2=amin
            j=nplpts-1
            call revrarr(losdist,j)
            call reviarr(layindx,j)
            j=nplpts
            call revrarr(rpl,j)
            call revrarr(mupl,j)
            call revrarr(x,j)
          else
c from the tangent point to z2
            amin=pi_2
	    call traceupr(zmin,z2,amin,re,z,rfndx,nz,th2,a2,
     &losdist,rpl,mupl,nplpts,x,layindx)
            do i=1,nplpts
               if(debug)write(6,50)i,rpl(i)-re,x(i),losdist(i)
50     format('ray1. ',i4,1P3e12.4)
            enddo
c store th2 while we calc. the other portion
            thtmp=th2
c also store the values of losdist,rpl,mupl,layindx,nplpts
            nplpts2=nplpts
c shift the values to the last part of the array
            do i=nplpts2,1,-1
               j=2*nz+1-(nplpts2-i)
               rpl(j)=rpl(i)
               x(j)=x(i)
               mupl(j)=mupl(i)
            enddo
            do i=nplpts2-1,1,-1
               j=2*nz-(nplpts2-1-i)
               losdist(j)=losdist(i)
               layindx(j)=layindx(i)
            enddo
            do i=1,nplpts*2+1
               if(debug)write(6,55)i,rpl(i)-re,x(i),losdist(i)
55     format('ray1b. ',i4,1P3e12.4)
            enddo
c now from the tangent point to z1
	    call traceupr(zmin,z1,amin,re,z,rfndx,nz,th2,atmp,
     &losdist,rpl,mupl,nplpts,x,layindx)
            atmp=pi-atmp
            do i=1,nplpts
               if(debug)write(6,60)i,rpl(i)-re,x(i),losdist(i)
60     format('ray2. ',i4,1P3e12.4)
            enddo
c now switch the order of the new stuff (zmin->z1 to z1->zmin)
            if(nplpts.gt.1)then
              xmax=x(nplpts)
              j=nplpts-1
              call revrarr(losdist,j)
              call reviarr(layindx,j)
              j=nplpts
              call revrarr(rpl,j)
              call revrarr(mupl,j)
              call revrarr(x,j)
c            do i=1,nplpts/2
c               tmp1=rpl(i)
c               rpl(i)=rpl(nplpts+1-i)
c               rpl(nplpts+1-i)=tmp1
c               tmp1=mupl(i)
c               mupl(i)=mupl(nplpts+1-i)
c               mupl(nplpts+1-i)=tmp1
c               tmp1=x(i)
c               x(i)=x(nplpts+1-i)
c               x(nplpts+1-i)=tmp1
c            enddo
c
c            do i=1,(nplpts-1)/2
c               tmp1=losdist(i)
c               losdist(i)=losdist(nplpts-i)
c               losdist(nplpts-i)=tmp1
c               tmp1=layindx(i)
c               layindx(i)=layindx(nplpts-i)
c               layindx(nplpts-i)=tmp1
c            enddo
c by my definition x(1)=0, so
              do i=1,nplpts
                 x(i)=xmax-x(i)
              enddo
c change the sign of mupl since we are going downward
              do i=1,nplpts-1
                 mupl(i)=-mupl(i)
              enddo
c now shift the old stuff to butt up against the new stuff
              do i=2,nplpts2
                 j=nplpts+i-1
                 rpl(j)=rpl(2*nz+1-(nplpts2-i))
                 mupl(j)=mupl(2*nz+1-(nplpts2-i))
                 x(j)=x(2*nz+1-(nplpts2-i))+x(nplpts)
              enddo
              do i=1,nplpts2-1
                 j=nplpts-1+i
                 losdist(j)=losdist(2*nz-(nplpts2-1-i))
                 layindx(j)=layindx(2*nz-(nplpts2-1-i))
              enddo
              nplpts=nplpts+nplpts2-1
            else
c shift the first stuff back to where it was originally
              do i=nplpts2,1,-1
                 j=2*nz+1-(nplpts2-i)
                 rpl(i)=rpl(j)
                 mupl(i)=mupl(j)
                 x(i)=x(j)
              enddo
              do i=nplpts2-1,1,-1
                 j=2*nz-(nplpts2-1-i)
                 losdist(i)=losdist(j)
                 layindx(i)=layindx(j)
              enddo
              nplpts=nplpts2
            endif
c check to see that atmp is equal to a1
            if(debug)
     &write(6,*)'raytrace. going from zmin to z1',zmin,z1,
     &a1/dr,atmp/dr,(a1-atmp)/dr,th2/dr
            th2=th2+thtmp
            if(debug)write(6,*)'raytrace ',nplpts,nplpts2
            do i=1,nplpts
               if(debug)write(6,70)i,rpl(i)-re,x(i),losdist(i)
70     format('ray3. ',i4,1P3e12.4)
            enddo
            if(debug)write(6,*)'*****'
          endif
        else
c will only go down, so will start at the opposite end (z2) since
c it is lower than z1.  To calc. the angle a2, use the formula
c (re+z1)*nz1*st1=(re+z2)*nz2*st2 to solve for a1 (st2=sine(a2))
	  call traceupr(z2,z1,a2,re,z,rfndx,nz,th2,atmp,
     &losdist,rpl,mupl,nplpts,x,layindx)
c check to see that atmp is equal to a1
          if(debug)
     &write(6,*)'raytrace. doing down a1 atmp',a1,atmp,a1-atmp
        endif
      endif
c zero-out all the other entries
      i=2*nz+1
      j=nplpts+1
      call zeroout1dr(rpl,i,j,i)
      call zeroout1dr(mupl,i,j,i)
      call zeroout1dr(x,i,j,i)
      i=2*nz
      j=nplpts
      call zeroout1dr(losdist,i,j,i)
      call zeroout1di(layindx,i,j,i)
      return
      end