c $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/refracray.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
        subroutine refracray(z1,tht0,re,rfndx,z,nz,thtap,iapp,
     &nplpts,layindx,losdist,rpl,mupl,x,ishadow,debug)
c 
c compute the refracted path from alt. z1 in the direction
c tht0
c
c input:
c z1 - initial altitude [km]
c tht0 - inital local zenith angle [rad]
c re - radius of earth [km]
c rfndx - refractive modulus (n-1) profile
c z - altitudes of the profile [km]
c nz - # of points in alt. profile
c iapp - switch for using tht0 as the true or apparent angle
c        [iapp=0 for true and iapp=1 for apparent]
c 
c output:
c thtap - apparent local zenith angle [rad]
c nplpts - # of points in path
c layindx(2*nz) - layer index for each segment
c losdist(2*nz) - distance of each segment [km]
c rpl(2*nz+1) - radial distance of each point [km]
c mupl(2*nz+1) - cosine of local zenith angle [rad]
c x(nz*2+1) - angular distance from the observer [rad]
c ishadow - flag for path that strikes the earth [1=strikes the earth]
c
        implicit none
        include 'consts.cmn'
c the rcs revision #
      character*20 rcsrev
      parameter (rcsrev='$Revision: 2.22 $')
c passed
        integer nz,nplpts,layindx(2*nz),ishadow,iapp
        real*8 rfndx(nz),z1,tht0,z(nz),re,thtap
        real*8 losdist(2*nz),rpl(2*nz+1),mupl(2*nz+1),x(nz*2+1)
c local
        integer i,j,itmax
        real*8 pc,pcpc,noref_dist,tht_shadow,thtap_shadow
        real*8 a1,psi,a2,diff,zst,zfn,dbendmin,dbendmax,dbend
        real*8 fip1mfi,xi(3),fi(3),li(3),xi_interp
        logical debug,usenewton
        integer klo1,khi1
        real*8 Hrf1,n01,nzfn,thtap_shadow_rfndx
c external
        real*8 rfintp
        external rfintp
c set the threshold at 1.0e-03 deg
c
        data dbendmin/1.7453292d-05/
        save dbendmin
c        data debug/.FALSE./
c        save debug
        data itmax/5/
        data usenewton/.FALSE./
        save usenewton
c
      if(iapp.eq.0)then
c check for the zenith angle allowed at this altitude (find the shadow).
c
c 1. calc. the angular dist. from the sfc to the toa along a ray
c    that is tangent at the sfc.  Then add the angle a2 to this to get the
c    true solar zenith angle at the beginning of the shadow at the sfc.
c
c 2. calc. the angular dist. from the sfc to the point z1 along a tangent
c    ray from the sfc.
c
c from sfc to toa
        a1=pi_2
        zst=z(nz)
        zfn=z(1)
	call raytrace(zst,zfn,a1,re,z,rfndx,nz,psi,a2,
     &losdist,rpl,mupl,nplpts,x,layindx,debug)
        tht_shadow=a2+psi
        if(debug)write(6,*)'for z=',zst,'the shadow begins at',
     &tht_shadow/dr
c from sfc to z1
        a1=pi_2
        zst=z(nz)
        zfn=z1
	call raytrace(zst,zfn,a1,re,z,rfndx,nz,psi,a2,
     &losdist,rpl,mupl,nplpts,x,layindx,debug)
        tht_shadow=psi+tht_shadow
        thtap_shadow=pi-a2
        if(debug)write(6,*)'ref. for z=',z1,
     &'the shadow begins at',tht_shadow/dr,
     &'and has an apparent sza of ',thtap_shadow/dr
        if(debug)write(6,*)'refracray. tht0-tht_shadow',
     &tht0-tht_shadow
c
        if((tht0-tht_shadow).gt.-smallx)then
c go from the observer to the earth
          zfn=z(nz)
          ishadow=1
        else
c slant path to space
          zfn=z(1)
          ishadow=0
        endif
        zst=z1
        a1=tht0
	call raytrace(zst,zfn,a1,re,z,rfndx,nz,psi,a2,
     &losdist,rpl,mupl,nplpts,x,layindx,debug)
      else if(iapp.eq.1)then
c
c let's try to find the apparent tht0 value
c
c check for the maximum solar zenith angle allowed at this altitude (find the shadow).
c
c 1. calc. the angular dist. from the sfc to the toa along a ray
c    that is tangent at the sfc.  Then add the angle a2 to this to get the
c    true solar zenith angle at the beginning of the shadow at the sfc.
c
c 2. calc. the angular dist. from the sfc to the point z1 along a tangent
c    ray from the sfc.
c
c from sfc to toa
        a1=pi_2
        zst=z(nz)
        zfn=z(1)
	call raytrace(zst,zfn,a1,re,z,rfndx,nz,psi,a2,
     &losdist,rpl,mupl,nplpts,x,layindx,debug)
        tht_shadow=a2+psi
        if(debug)write(6,*)'for z=',zst,'the shadow begins at',
     &tht_shadow/dr
c from sfc to z1
        a1=pi_2
        zst=z(nz)
        zfn=z1
	call raytrace(zst,zfn,a1,re,z,rfndx,nz,psi,a2,
     &losdist,rpl,mupl,nplpts,x,layindx,debug)
        tht_shadow=psi+tht_shadow
c the angle at the observer
        thtap_shadow=pi-a2
        if(debug)write(6,*)'ref. for z=',z1,
     &'the shadow begins at',tht_shadow/dr,
     &'and has an apparent sza of ',thtap_shadow/dr
        if(debug)write(6,*)'refracray. tht0-tht_shadow',
     &(tht0-tht_shadow)/dr
c now what does the tabulated values of rfndx give for the shadow
        call findex(z,nz,zfn,klo1,khi1)
        call rfscal(rfndx(klo1),rfndx(khi1),z(klo1),z(khi1),Hrf1,
     &n01)
        nzfn=1.0d0+rfintp(Hrf1,n01,zfn)
        thtap_shadow_rfndx=(pi-asin((re*(1.+rfndx(nz))/
     &((re+zfn)*nzfn))))
        if(debug)write(6,*)'refracray. the apparent sza for the shadow ',
     &'from rfndx is ',thtap_shadow_rfndx/dr
        if(dbendmin.lt.abs(thtap_shadow-thtap_shadow_rfndx))
     &dbendmin=abs(thtap_shadow-thtap_shadow_rfndx)
        if(debug)write(6,*)'have set dbendmin to ',dbendmin/dr
c
        if((tht0-tht_shadow).gt.-dbendmin)then
c in the shadow
          ishadow=1
          thtap=tht0
          if(debug)write(6,*)'refracray.  This point ',
     &'is in the shadow'
        else
c the amount of bending for the shadow ray at this altitude
          dbendmax=tht_shadow-thtap_shadow
          ishadow=0
c first guess is the lesser of tht0 and thtap_shadow-smallx
c (tht0 will be the answer if the index to refraction is 1.0)
          a1=tht0
          if(tht0.gt.(thtap_shadow-dbendmin))
     &a1=thtap_shadow-dbendmin
          zst=z1
          zfn=z(1)
          if(debug)write(6,*)'refracray. calling raytrace for 1st ',
     &'guess',a1/dr,zst,zfn
	  call raytrace(zst,zfn,a1,re,z,rfndx,nz,psi,a2,
     &losdist,rpl,mupl,nplpts,x,layindx,debug)
          if(debug)write(6,*)'refracray. rpl ',rpl(1)-re,rpl(nplpts)-re
          if(debug)write(6,*)'zst zfn ',zst,zfn,a1/dr,tht0/dr
c wish to find root of psi+a2-tht0=0
          fi(1)=(psi+a2)
          xi(1)=a1
          if(debug)write(6,1020)'xi fi ',xi(1)/dr,
     &fi(1)/dr,(fi(1)-tht0)/dr
1020  format(a11,1P3E16.8)
          if(abs(fi(1)-tht0).gt.dbendmin)then
c second guess will be tht0- amount bent for tht0
            dbend=fi(1)-tht0
            a1=xi(1)-dbend
            if(abs(zst-z(1)).lt.smallx)then
              if(a1-pi_2.lt.smallx)a1=pi_2+smallx
            endif
	    call raytrace(zst,zfn,a1,re,z,rfndx,nz,psi,a2,
     &losdist,rpl,mupl,nplpts,x,layindx,debug)
            fi(2)=(psi+a2)
            xi(2)=a1
            xi(3)=xi(2)
            if(debug)write(6,1020)'xi fi ',xi(2)/dr,
     &fi(2)/dr,(fi(2)-tht0)/dr
            if(usenewton)then
              if(debug)write(6,*)'refracray. will use newton for thtap'
c now calc. the new position (will use 3 trys)
              i=1
              diff=abs(dbend)
              do while((i.le.itmax).and.(diff.gt.dbendmin))
                 fip1mfi=fi(2)-fi(1)
                 if(abs(fip1mfi).gt.smallx)then
                   xi(3)=xi(2)+(tht0-fi(2))/((fip1mfi)/(xi(2)-xi(1)))
                 else
                   xi(3)=xi(2)
                   goto 10
                 endif
                 if(i.eq.1)then
                   if(abs(xi(3)-xi(2)).lt.abs(xi(3)-xi(1)))then
                     xi(1)=xi(2)
                     fi(1)=fi(2)
                   endif
                 else
                   xi(1)=xi(2)
                   fi(1)=fi(2)
                 endif
                 xi(2)=xi(3)
                 if(abs(zst-z(1)).lt.smallx)then
                   if(xi(2)-pi_2.lt.smallx)xi(2)=pi_2+smallx
                 endif
	         call raytrace(zst,zfn,xi(2),re,z,rfndx,nz,psi,a2,
     &losdist,rpl,mupl,nplpts,x,layindx,debug)
                 fi(2)=(psi+a2)
                 diff=abs(fi(2)-tht0)
                 if(debug)write(6,1020)'xi(2) fi(2) ',xi(2)/dr,
     &fi(2)/dr,(fi(2)-tht0)/dr
                 i=i+1
              enddo
            else
c use a quad fit
c solve for point in the middle
              fip1mfi=fi(2)-fi(1)
              if(abs(fip1mfi).gt.smallx)then
                xi(3)=xi(2)+(tht0-fi(2))/((fip1mfi)/(xi(2)-xi(1)))
              else
                xi(3)=xi(2)
                goto 10
              endif
	      call raytrace(zst,zfn,xi(3),re,z,rfndx,nz,psi,a2,
     &losdist,rpl,mupl,nplpts,x,layindx,debug)
              fi(3)=(psi+a2)
              if(debug)write(6,1020)'xi fi ',xi(3)/dr,
     &fi(3)/dr,(fi(3)-tht0)/dr
c now do quad interpolate
              xi_interp=0.0d0
              do i=1,3
                 li(i)=1.0d0
                 do j=1,3
                    if(i.ne.j)li(i)=li(i)*(tht0-fi(j))/(fi(i)-fi(j))
                 enddo
                 xi_interp=xi_interp+li(i)*xi(i)
                 if(debug)write(6,*)'refracray. li ',li(i),xi_interp/dr
              enddo
              xi(3)=xi_interp
            endif
10          xi(1)=xi(3)
          endif
          thtap=xi(1)
c now let's check it
	  call raytrace(zst,zfn,thtap,re,z,rfndx,nz,psi,a2,
     &losdist,rpl,mupl,nplpts,x,layindx,debug)
          if(debug)write(6,*)'the apparent sza is ',thtap/dr,
     &'for tht0 ',tht0/dr,(psi+a2)/dr
c          i=1
c          if(debug)write(6,1010)i,rpl(i)-re,acos(mupl(i))/dr,
c     &x(i)/dr,zst
        endif
      endif
        return
1000    format(1P3e12.4)
1010    format(i5,1P4e12.4,i5)
        end