c************************************************************************

      complex function stalk_forward_scat(phi,aa,zz0,pol,z)

c************************************************************************

      integer vv,hh,hv
      parameter(vv=1)
      parameter(hh=2)
      parameter(hv=3)
      
      complex pxxfsv_s,pzzfsv_s,pxzfsv_s
      complex pyyfsh_s
      complex pxx,pyy,pzz,pxz
      
      real z0
      real a,alpha_s,ell_0
      real theta_i,phi_i,k_0
      real z
      real pi,coef
      real phi,aa,zz0
      real dl,drho
      
      integer pol
      
      external pxxfsv_s,pzzfsv_s,pxzfsv_s
      external pyyfsh_s
      
      real d,del_leaf,z_lc,hole
      common /layer_struct/ d,del_leaf,z_lc,hole
      
      integer s_tensor_flag,s_tau_flag,z0_flag,s_field_flag(2,4)
      common /stalk_flags/ s_tensor_flag,s_tau_flag
     &     ,z0_flag,s_field_flag

      common /stalk_struct/ a,alpha_s,ell_0
      common /incident/ theta_i,phi_i,k_0

      include 'error.incld'

      alpha_s=phi
      a=aa
      z0=zz0
      
      drho=a
      dl=sqrt(1.0+drho**2)

c      ell_0=dl*(z0+d) <- old ell_0
c     the new ell_0 is

      ell_0=dl*d

      pi=3.141592654      
      coef=k_0**2/(4.0*pi)

c      print*,'in stalk_forward_scat'

c      print*,z
      
      s_tensor_flag=0
          
      if (pol .eq. vv) then

        if (z .ge. z0) then

          stalk_forward_scat=cmplx(0.0,0.0)

        else
          
          call integral(pxxfsv_s,z0,z,pxx)
          
          call integral(pzzfsv_s,z0,z,pzz)

          call integral(pxzfsv_s,z0,z,pxz)

          stalk_forward_scat=coef*
     &    (cos(theta_i)**2*pxx
     &     -2.0*cos(theta_i)*sin(theta_i)*pxz
     &     +sin(theta_i)**2*pzz)

c        this is a test line only!!!
C          if(z .lt.-0.95*d)then
C          call integral(pyyfsh_s,z0,z,pyy)
C          print*,'z,pxx,pyy,pzz',z,pxx,pyy,pzz
C          print*,'stalk_forward_scat
C          endif
          

        endif

      else if (pol .eq. hh) then

        if (z .gt. z0) then 

          stalk_forward_scat=cmplx(0.0,0.0)

        else

          call integral(pyyfsh_s,z0,z,pyy)

c          print*,pyy,z0,z,d

          stalk_forward_scat=coef*pyy

        endif

      else if (pol .eq. hv) then

        stalk_forward_scat=cmplx(0.0,0.0)
        
        write(error,*)'pol = hv in func leaf_forward_scat was',
     &       'encountered but not expected'
        code=26
        call the_death_code(warn,error,code)
       
      else 

        write(error,*)'value of pol = ',pol,
     &   ' not 1,2 or 3 in fun stalk_forward_scat'
        code=27
        call the_death_code(die,error,code)

      endif

      return
      end

      
c************************************************************************
      
      complex function pxxfsv_s(z)

c************************************************************************

       complex p11,p22,p33

       real z,ell
       real dl,cosb,sinb,drho,sina,cosa

       real a,alpha_s,ell_0
       common /stalk_struct/ a,alpha_s,ell_0
       
       real d,del_leaf,z_lc,hole
       common /layer_struct/ d,del_leaf,z_lc,hole

c       print*,'in pxxfsv'

       drho=a
       dl=sqrt(1.0+drho**2)

       ell=dl*(z+d)

       cosb=1.0/dl
       sinb=drho/dl

       cosa=cos(alpha_s)
       sina=sin(alpha_s)

       call stalk_tensor(ell,ell_0, p11,p22,p33)

       pxxfsv_s=(cosa**2*cosb**2*p11
     &  +sina**2*p22+cosa**2*sinb**2*p33)*dl

       return
       end
      

c************************************************************************
      
      complex function pzzfsv_s(z)

c************************************************************************

       complex p11,p22,p33
       real z,ell
       real dl,cosb,sinb,drho,sina,cosa

       real a,alpha_s,ell_0
       common /stalk_struct/ a,alpha_s,ell_0
       
       real d,del_leaf,z_lc,hole
       common /layer_struct/ d,del_leaf,z_lc,hole
       
c       print*,'in pzzfsv'

       drho=a
       dl=sqrt(1.0+drho**2)

       ell=dl*(z+d)

       cosb=1.0/dl
       sinb=drho/dl

       cosa=cos(alpha_s)
       sina=sin(alpha_s)

       call stalk_tensor(ell,ell_0, p11,p22,p33)

       pzzfsv_s=(sinb**2*p11+cosb**2*p33)*dl

       return
       end


c************************************************************************
      
      complex function pxzfsv_s(z)

c************************************************************************

       complex p11,p22,p33
       real z,ell
       real dl,cosb,sinb,drho,sina,cosa

       real a,alpha_s,ell_0
       common /stalk_struct/ a,alpha_s,ell_0
       
       real d,del_leaf,z_lc,hole
       common /layer_struct/ d,del_leaf,z_lc,hole
       
c       print*,'in pzzfsv'

       drho=a
       dl=sqrt(1.0+drho**2)

       ell=dl*(z+d)

       cosb=1.0/dl
       sinb=drho/dl

       cosa=cos(alpha_s)
       sina=sin(alpha_s)

       call stalk_tensor(ell,ell_0, p11,p22,p33)

       pxzfsv_s=(cosa*cosb*sinb*(p33-p11))*dl

       return
       end
       
c************************************************************************
      
      complex function pyyfsh_s(z)

c************************************************************************

       complex p11,p22,p33
       real z,ell
       real dl,cosb,sinb,drho,sina,cosa

       real a,alpha_s,ell_0
       common /stalk_struct/ a,alpha_s,ell_0
       
       real d,del_leaf,z_lc,hole
       common /layer_struct/ d,del_leaf,z_lc,hole

c       print*,'in pyyfsv'
       
       drho=a
       dl=sqrt(1.0+drho**2)

       ell=dl*(z+d)

       cosb=1.0/dl
       sinb=drho/dl

       cosa=cos(alpha_s)
       sina=sin(alpha_s)

       call stalk_tensor(ell,ell_0, p11,p22,p33)

       pyyfsh_s=(sina**2*cosb**2*p11
     &  +cosa**2*p22+sina**2*sinb**2*p33)*dl

       return
       end