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

      subroutine stalk_power(ss_leaf)

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


      complex s,ss_stalk(4,4,4)

      integer i,j,k

      do 10 i=1,4

        do 20 j=1,4

          do 30 k=1,4

          if (i .eq. 1) then

            if (j .eq. 1) then

              if (k .eq. 1) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else if(k .eq. 2) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else if(k .eq. 3) then

                ss_stalk(i,j,k)=ss_stalk(i,j,2)

              else if(k .eq. 4) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s                

              else

                print*, 'error in stalk_power a'

              endif

            else if(j .eq. 2) then

              if (k .eq. 1) then

                ss_stalk(i,j,k)=conjg(ss_stalk(i,k,j))

              else if(k .eq. 2) then
                
                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else if(k .eq. 3) then

                ss_stalk(i,j,k)=ss_stalk(i,j,2)

              else if(k .eq. 4) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else

                print*, 'error in stalk_power b'

              endif

            else if(j .eq. 3) then

              if(k .eq. 1) then

                ss_stalk(i,j,k)=conjg(ss_stalk(i,k,j))

              else if(k .eq. 2) then

                ss_stalk(i,j,k)=conjg(ss_stalk(i,k,j))

              else if(k .eq. 3) then

                ss_stalk(i,j,k)=ss_stalk(i,2,2)

              else if(k .eq. 4) then

                ss_stalk(i,j,k)=ss_stalk(i,2,k)

              else

                print*, 'error in  stalk_power c'

              endif

            if (j .eq. 4) then

              if (k .eq. 1) then

                ss_stalk(i,j,k)=conjg(ss_stalk(i,k,j))

              else if (k .eq. 2) then

                ss_stalk(i,j,k)=conjg(ss_stalk(i,k,j))

              else if (k .eq. 3) then

               ss_stalk(i,j,k)=conjg(ss_stalk(i,k,j))
               
              else if (k .eq. 4) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else

                print*, 'error in  stalk_power d'
                
              endif
              
            else

              print*, 'error in  stalk_power dd'
                
            endif

          else if(i .eq. 2) then

            if (j .eq. 1) then

              if (k .eq. 1) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else if(k .eq. 2) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else if(k .eq. 3) then

                ss_stalk(i,j,k)=ss_stalk(i,j,2)

              else if(k .eq. 4) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s                

              else

                print*, 'error in stalk_power e'

              endif

            else if(j .eq. 2) then

              if (k .eq. 1) then

                ss_stalk(i,j,k)=conjg(ss_stalk(i,k,j))

              else if(k .eq. 2) then
                
                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else if(k .eq. 3) then

                ss_stalk(i,j,k)=ss_stalk(i,j,2)

              else if(k .eq. 4) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else

                print*, 'error in stalk_power f'

              endif

            else if(j .eq. 3) then

              if(k .eq. 1) then

                ss_stalk(i,j,k)=conjg(ss_stalk(i,k,j))

              else if(k .eq. 2) then

                ss_stalk(i,j,k)=conjg(ss_stalk(i,k,j))

              else if(k .eq. 3) then

                ss_stalk(i,j,k)=ss_stalk(i,2,2)

              else if(k .eq. 4) then

                ss_stalk(i,j,k)=ss_stalk(i,2,k)

              else

                print*, 'error in  stalk_power g'

              endif

            if (j .eq. 4) then

              if (k .eq. 1) then

                ss_stalk(i,j,k)=conjg(ss_stalk(i,k,j))

              else if (k .eq. 2) then

                ss_stalk(i,j,k)=conjg(ss_stalk(i,k,j))

              else if (k .eq. 3) then

               ss_stalk(i,j,k)=conjg(ss_stalk(i,k,j))
               
              else if (k .eq. 4) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else

                print*, 'error in  stalk_power h'
                
              endif

            else

              print*, 'error in  stalk_power hh'

            endif                
            
          else if(i .eq. 3) then

            if (j .eq. 1) then

              if (k .eq. 1) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else if(k .eq. 2) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else if(k .eq. 3) then

                call stalk_state_int(i,j,k,s)

              else if(k .eq. 4) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s                

              else

                print*, 'error in stalk_power i'

              endif

            else if(j .eq. 2) then

              if (k .eq. 1) then

                ss_stalk(i,j,k)=conjg(ss_stalk(i,k,j))

              else if(k .eq. 2) then
                
                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else if(k .eq. 3) then
                
                call stalk_state_int(i,j,k,s)

              else if(k .eq. 4) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else

                print*, 'error in stalk_power j'

              endif

            else if(j .eq. 3) then

              if(k .eq. 1) then

                ss_stalk(i,j,k)=conjg(ss_stalk(i,k,j))

              else if(k .eq. 2) then

                ss_stalk(i,j,k)=conjg(ss_stalk(i,k,j))

              else if(k .eq. 3) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s
                
              else if(k .eq. 4) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else

                print*, 'error in  stalk_power k'

              endif

            if (j .eq. 4) then

              if (k .eq. 1) then

                ss_stalk(i,j,k)=conjg(ss_stalk(i,k,j))

              else if (k .eq. 2) then

                ss_stalk(i,j,k)=conjg(ss_stalk(i,k,j))

              else if (k .eq. 3) then

               ss_stalk(i,j,k)=conjg(ss_stalk(i,k,j))
               
              else if (k .eq. 4) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else

                print*, 'error in  stalk_power l'
                
              endif
              
            else

              print*, 'error in  stalk_power ll'
              
            endif
            
          if (i .eq. 4) then

            if (j .eq. 1) then

              if (k .eq. 1) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else if(k .eq. 2) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else if(k .eq. 3) then

                ss_stalk(i,j,k)=ss_stalk(i,j,2)

              else if(k .eq. 4) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s                

              else

                print*, 'error in stalk_power m'

              endif

            else if(j .eq. 2) then

              if (k .eq. 1) then
                
                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else if(k .eq. 2) then
                
                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else if(k .eq. 3) then

                ss_stalk(i,j,k)=ss_stalk(i,j,2)

              else if(k .eq. 4) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else

                print*, 'error in stalk_power n'

              endif

            else if(j .eq. 3) then

              if(k .eq. 1) then

                ss_stalk(i,j,k)=ss_stalk(i,2,k)

              else if(k .eq. 2) then

                ss_stalk(i,j,k)=ss_stalk(i,2,k)

              else if(k .eq. 3) then

                ss_stalk(i,j,k)=ss_stalk(i,2,2)

              else if(k .eq. 4) then

                ss_stalk(i,j,k)=ss_stalk(i,2,k)

              else

                print*, 'error in  stalk_power o'

              endif

            if (j .eq. 4) then

              if (k .eq. 1) then
                
                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else if (k .eq. 2) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s
                
              else if (k .eq. 3) then

               ss_stalk(i,j,k)=ss_stalk(i,j,2)
               
              else if (k .eq. 4) then

                call stalk_state_int(i,j,k,s)

                ss_stalk(i,j,k)=s

              else

                print*, 'error in  stalk_power p'
                
              endif

            else

              print*, 'error in  stalk_power q'
              
            endif

          else

            print*, 'error in  stalk_power r'

          endif

            
   30     continue  

   20   continue

   10 continue

      return

      end

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

      subroutine stalk_state_int(i,j,k, p)

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

      complex h_stalk,p
      real pi,mpi
      integer i,j,k

      common /index/ i,j,k

      pi=4.0*atan(1.0)
      mpi=-pi
      

      external h_stalk
      call qgausx(h_stalk,mpi,pi,p)

      return

      end
      

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

      function h_stalk(phii)

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

      complex h_stalk,s,g_stalk

      real amin,amax,z0min,z0max
      real phi_s,a_s,z0
      real phii
      
      external g_stalk
      
      common /stalk_limits/ amin,amax,z0min,z0max
      common /stalk_state/ phi_s,a_s,z0

      phi_s=phii

      call qgausy(g_stalk,amin,amax,s)

      h_stalk=s

      return

      end

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

      function g_stalk(aa)

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


      complex g_stalk,i_stalk,s

      real amin,amax,z0min,z0max
      real phi_s,a_s,z0
      real aa

      external i_stalk
      
      common /stalk_limits/ amin,amax,z0min,z0max
      common /stalk_state/ phi_s,a_s,z0

      a_s=aa

      call qgausy(i_stalk,z0min,z0max,s)

      g_stalk=s

      return

      end

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

      function i_stalk(zz0)

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


      complex i_stalk,stalk_state_ave

      real phi_s,a_s,z0
      real zz0
      
      common /stalk_state/ phi_s,a_s,z0
      common /index/ i,j,k
      
      z0=zz0
      
      i_stalk=stalk_state_ave(phi_s,a_s,z0,i,j,k)

      return

      end

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

      function stalk_state_ave(alpha_s,a,z0,i,j,k)

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

      complex p11vva_s,p33vva_s,p13vva_s
      complex p11vvb_s,p33vvb_s,p13vvb_s
      complex p11vvc_s,p33vvc_s,p13vvc_s
      complex p11vvd_s,p33vvd_s,p13vvd_s

      complex p22hha_s
      complex p22hhb_s
      complex p22hhc_s
      complex p22hhd_s

      complex p12hva_s,p23hva_s
      complex p12hvb_s,p23hvb_s
      complex p12hvc_s,p23hvc_s
      complex p12hvd_s,p23hvd_s

      complex p11a,p12a,p13a,p22a,p23a,p33a
      complex p11b,p12b,p13b,p22b,p23b,p33b
      complex p11c,p12c,p13c,p22c,p23c,p33c
      complex p11d,p12d,p13d,p22d,p23d,p33d

      complex stalk_state_ave
      
      real dm,z0
      real a,alpha_s
      real theta_i,phi_i,k_0
      real zero
      real stalk_alpha_pdf,stalk_a_pdf,stalk_z0_pdf
      real dm,del_leaf,zlc
      
      integer i,j,k

      external p11vva_s,p33vva_s,p13vva_s
      external p11vvb_s,p33vvb_s,p13vvb_s
      external p11vvc_s,p33vvc_s,p13vvc_s
      external p11vvd_s,p33vvd_s,p13vvd_s

      external p22hha_s
      external p22hhb_s
      external p22hhc_s
      external p22hhd_s

      external p12hva_s,p23hva_s
      external p12hvb_s,p23hvb_s
      external p12hvc_s,p23hvc_s
      external p12hvd_s,p23hvd_s      
      
      common /stalk_struct/ a,alpha_s
      common /incident/ theta_i,phi_i,k_0
      common /layer_struct/ dm,del_leaf,z_lc

      zero=0.0

      if (i .eq. 1) then

        if (j .eq. 1) then

          if (k .eq. 1) then

            call qgaus(p11vva_s,dm,z0,p11a)

            call qgaus(p33vva_s,dm,z0,p33a)

            call qgaus(p13vva_s,dm,z0,p13a)

            stalk_state_ave=stalk_alpha_pdf(alpha)
             *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)**4*p11a*conjg(p11a)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (4.0*p13a*conjg(p13a)
     &       +p11a*conjg(p33a)+p33a*conjg(p11a))
     &       +sin(theta_i)**4*p33a*conjg(p33a))


          else if (k .eq. 2) then

            call qgaus(p11vva_s,dm,z0,p11a)

            call qgaus(p33vva_s,dm,z0,p33a)

            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)*stalk_a_pdf(a)
             *stalk_z0_pdf(z0)*
     &       (-cos(theta_i)**4*p11a*conjg(p11b)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (p11a*conjg(p33b)-p33a*conjg(p11b))
     &       +sin(theta_i)**4*p33a*conjg(p33b))
     
          else if (k .eq. 3) then

            call qgaus(p11vva_s,dm,z0,p11a)

            call qgaus(p33vva_s,dm,z0,p33a)

            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)


            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &      *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &      (-cos(theta_i)**4*p11a*conjg(p11c)
     &      +cos(theta_i)**2*sin(theta_i)**2*
     &      (p11a*conjg(p33c)-p33a*conjg(p11c))
     &      +sin(theta_i)**4*p33a*conjg(p33c))

          else if (k .eq. 4) then

            call qgaus(p11vva_s,dm,z0,p11a)

            call qgaus(p33vva_s,dm,z0,p33a)

            call qgaus(p13vva_s,dm,z0,p13a)

            call qgaus(p11vvd_s,dm,z0,p11d)

            call qgaus(p33vvd_s,dm,z0,p33d)

            call qgaus(p13vvd_s,dm,z0,p13d)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)**4*p11a*conjg(p11d)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (-4.0*p13a*conjg(p13d)+
     &       p11a*conjg(p33d)+p33a*conjg(p11d))
     &       +sin(theta_i)**4*p33a*conjg(p33d))
                      
          else

            print*, 'error in stalk_state_ave a'

          endif

        else if (j .eq. 2) then
          
          if (k .eq. 1) then


            call qgaus(p11vva_s,dm,z0,p11a)

            call qgaus(p33vva_s,dm,z0,p33a)

            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)


            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)**4*p11b*conjg(p11a)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (p33b*conjg(p11a)-p11b*conjg(p33a))
     &       +sin(theta_i)**4*p33b*conjg(p33a))            

          else if (k .eq. 2) then


            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)


            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)**4*p11b*conjg(p11b)
     &       -cos(theta_i)**2*sin(theta_i)**2*
     &       (p33b*conjg(p11b)+p11b*conjg(p33b))
     &       +sin(theta_i)**4*p33b*conjg(p33b))

          else if (k .eq. 3)

            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)

            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)
            
            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)**4*p11b*conjg(p11c)
     &       -cos(theta_i)**2*sin(theta_i)**2*
     &       (p33b*conjg(p11c)+p11b*conjg(p33c))
     &       +sin(theta_i)**4*p33b*conjg(p33c))            

          else if (k .eq. 4)

            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)

            call qgaus(p11vvd_s,dm,z0,p11d)

            call qgaus(p33vvd_s,dm,z0,p33d)
            
            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)**4*p11b*conjg(p11d)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (p33b*conjg(p11d)-p11b*conjg(p33d))
     &       +sin(theta_i)**4*p33b*conjg(p33d))

          else

            print*, 'error in stalk_state_ave b'

          endif

        else if(j .eq. 3) then

          if (k .eq. 1) then


            call qgaus(p11vva_s,dm,z0,p11a)

            call qgaus(p33vva_s,dm,z0,p33a)

            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)


            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)**4*p11c*conjg(p11a)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (p33c*conjg(p11a)-p11c*conjg(p33a))
     &       +sin(theta_i)**4*p33c*conjg(p33a))            

          else if (k .eq. 2) then


            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)

            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)
            
            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)**4*p11c*conjg(p11b)
     &       -cos(theta_i)**2*sin(theta_i)**2*
     &       (p33c*conjg(p11b)+p11c*conjg(p33b))
     &       +sin(theta_i)**4*p33c*conjg(p33b))

          else if (k .eq. 3)

            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)
            
            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)**4*p11c*conjg(p11c)
     &       -cos(theta_i)**2*sin(theta_i)**2*
     &       (p33c*conjg(p11c)+p11c*conjg(p33c))
     &       +sin(theta_i)**4*p33c*conjg(p33c))            

          else if (k .eq. 4)

            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)

            call qgaus(p11vvd_s,dm,z0,p11d)

            call qgaus(p33vvd_s,dm,z0,p33d)
            
            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)**4*p11c*conjg(p11d)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (p33c*conjg(p11d)-p11c*conjg(p33d))
     &       +sin(theta_i)**4*p33c*conjg(p33d))

          else

            print*, 'error in stalk_state_ave c'

          endif

        else if (j .eq. 4) then

          if (k .eq. 1) then

            call qgaus(p11vva_s,dm,z0,p11a)

            call qgaus(p33vva_s,dm,z0,p33a)

            call qgaus(p13vva_s,dm,z0,p13a)
            
            call qgaus(p11vvd_s,dm,z0,p11d)

            call qgaus(p33vvd_s,dm,z0,p33d)

            call qgaus(p13vvd_s,dm,z0,p13d)

            stalk_state_ave=stalk_alpha_pdf*stalk_a_pdf*stalk_z0_pdf*
     &       (cos(theta_i)**4*p11d*conjg(p11a)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (-4.0*p13d*conjg(p13a)
     &       +p11d*conjg(p33a)+p33d*conjg(p11a))
     &       +sin(theta_i)**4*p33d*conjg(p33a))

          else if (k .eq. 2) then

            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)

            call qgaus(p11vvd_s,dm,z0,p11d)

            call qgaus(p33vvd_s,dm,z0,p33d)

            stalk_state_ave=stalk_alpha_pdf*stalk_a_pdf*stalk_z0_pdf*
     &       (-cos(theta_i)**4*p11d*conjg(p11b)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (p11d*conjg(p33b)-p33d*conjg(p11b))
     &       +sin(theta_i)**4*p33d*conjg(p33b))
     
          else if (k .eq. 3) then

            call qgaus(p11vvd_s,dm,z0,p11d)

            call qgaus(p33vvd_s,dm,z0,p33d)

            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)


            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)**4*p11d*conjg(p11c)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (p11d*conjg(p33c)-p33d*conjg(p11c))
     &       +sin(theta_i)**4*p33d*conjg(p33c))

          else if (k .eq. 4) then

            call qgaus(p11vvd_s,dm,z0,p11d)

            call qgaus(p33vvd_s,dm,z0,p33d)

            call qgaus(p13vvd_s,dm,z0,p13d)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)**4*p11d*conjg(p11d)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (4.0*p13d*conjg(p13d)+
     &       p11d*conjg(p33d)+p33d*conjg(p11d))
     &       +sin(theta_i)**4*p33d*conjg(p33d))
                      
          else

            print*, 'error in stalk_state_ave d'

          endif
          
        else

          print*, 'error in stalk_state_ave dd'
            
        endif

      else if (i .eq. 2) then

        if (j .eq. 1) then

          if (k .eq. 1) then

            call qgaus(p22hha_s,dm,z0,p22a)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22a*conjg(p22a)

          else if (k .eq. 2) then

            call qgaus(p22hha_s,dm,z0,p22a)

            call qgaus(p22hhb_s,dm,z0,p22b)            

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22a*conjg(p22b)
         
         else if (k .eq.3) then

            call qgaus(p22hha_s,dm,z0,p22a)

            call qgaus(p22hhc_s,dm,z0,p22c)            

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22a*conjg(p22c)
                     
         else if (k .eq.4) then

            call qgaus(p22hha_s,dm,z0,p22a)

            call qgaus(p22hhd_s,dm,z0,p22d)            

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22a*conjg(p22d)            
        
         else

           print*, 'error in stalk_state_ave e'

         endif

       else if (j .eq.2)

         if (k .eq. 1) then

            call qgaus(p22hhb_s,dm,z0,p22b)            

            call qgaus(p22hha_s,dm,z0,p22a)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22b*conjg(p22a)

          else if (k .eq. 2) then

            call qgaus(p22hhb_s,dm,z0,p22b)            

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22b*conjg(p22b)
         
         else if (k .eq.3) then

            call qgaus(p22hha_s,dm,z0,p22b)

            call qgaus(p22hhc_s,dm,z0,p22c)            

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22b*conjg(p22c)
                     
         else if (k .eq.4) then

            call qgaus(p22hha_s,dm,z0,p22b)

            call qgaus(p22hhd_s,dm,z0,p22d)            

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22b*conjg(p22d)            
        
         else

           print*, 'error in stalk_state_ave f'

         endif

       else if (j .eq. 3)

         if (k .eq. 1) then

            call qgaus(p22hhc_s,dm,z0,p22c)            

            call qgaus(p22hha_s,dm,z0,p22a)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22c*conjg(p22a)

          else if (k .eq. 2) then

            call qgaus(p22hhc_s,dm,z0,p22c)

            call qgaus(p22hhb_s,dm,z0,p22b)            

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22c*conjg(p22b)
         
         else if (k .eq.3) then

            call qgaus(p22hhc_s,dm,z0,p22c)            

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22c*conjg(p22c)
                     
         else if (k .eq.4) then

            call qgaus(p22hhc_s,dm,z0,p22c)

            call qgaus(p22hhd_s,dm,z0,p22d)            

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22c*conjg(p22d)            
        
         else

           print*, 'error in stalk_state_ave g'

         endif

       else if (j .eq. 4)

         if (k .eq. 1) then

            call qgaus(p22hhd_s,dm,z0,p22d)            

            call qgaus(p22hha_s,dm,z0,p22a)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22d*conjg(p22a)

          else if (k .eq. 2) then

            call qgaus(p22hhd_s,dm,z0,p22d)

            call qgaus(p22hhb_s,dm,z0,p22b)            

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22d*conjg(p22b)
         
         else if (k .eq. 3) then

           call qgaus(p22hhd_s,dm,z0,p22d)            

           call qgaus(p22hhc_s,dm,z0,p22c)            

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22d*conjg(p22c)
                     
         else if (k .eq. 4) then

            call qgaus(p22hhd_s,dm,z0,p22d)            

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22d*conjg(p22d)            
        
           else

             print*, 'error in stalk_state_ave h'

           endif

             else

               print*, 'error in stalk_state_ave hh'

             endif
             

       else if (i .eq. 3) then

         if (j .eq. 1) then

           if (k .eq. 1) then

            call qgaus(p12hva_s,dm,z0,p12a)

            call qgaus(p23hva_s,dm,z0,p23a)
            
            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)**2*p12a*conjg(p12a)
     &       +sin(theta_i)**2*p23a*conjg(p23a))

            
           else if (k .eq. 2) then

            call qgaus(p12hva_s,dm,z0,p12a)

            call qgaus(p23hva_s,dm,z0,p23a)
            
            call qgaus(p12hvb_s,dm,z0,p12b)

            call qgaus(p23hvb_s,dm,z0,p23b)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)**2*p12a*conjg(p12b)
     &       +sin(theta_i)**2*p23a*conjg(p23b))

           else if (k .eq. 3) then

            call qgaus(p12hva_s,dm,z0,p12a)

            call qgaus(p23hva_s,dm,z0,p23a)
            
            call qgaus(p12hvc_s,dm,z0,p12c)

            call qgaus(p23hvc_s,dm,z0,p23c)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)**2*p12a*conjg(p12c)
     &       +sin(theta_i)**2*p23a*conjg(p23c))

           else if (k .eq. 4) then

            call qgaus(p12hva_s,dm,z0,p12a)

            call qgaus(p23hva_s,dm,z0,p23a)
            
            call qgaus(p12hvd_s,dm,z0,p12d)

            call qgaus(p23hvd_s,dm,z0,p23d)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)**2*p12a*conjg(p12d)
     &       +sin(theta_i)**2*p23a*conjg(p23d))

          else

            print*, 'error in stalk_state_ave i'

          endif
          
        else

          print*, 'error in stalk_state_ave ii'

        endif

         if (j .eq. 2) then

           if (k .eq. 1) then

            call qgaus(p12hvb_s,dm,z0,p12b)

            call qgaus(p23hvb_s,dm,z0,p23b)

            call qgaus(p12hva_s,dm,z0,p12a)

            call qgaus(p23hva_s,dm,z0,p23a)
            
            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)**2*p12b*conjg(p12a)
     &       +sin(theta_i)**2*p23b*conjg(p23a))

            
           else if (k .eq. 2) then

            call qgaus(p12hvb_s,dm,z0,p12b)

            call qgaus(p23hvb_s,dm,z0,p23b)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)**2*p12b*conjg(p12b)
     &       +sin(theta_i)**2*p23b*conjg(p23b))

           else if (k .eq. 3) then

            call qgaus(p12hvb_s,dm,z0,p12b)

            call qgaus(p23hvb_s,dm,z0,p23b)
            
            call qgaus(p12hvc_s,dm,z0,p12c)

            call qgaus(p23hvc_s,dm,z0,p23c)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)**2*p12b*conjg(p12c)
     &       +sin(theta_i)**2*p23b*conjg(p23c))

           else if (k .eq. 4) then

            call qgaus(p12hvb_s,dm,z0,p12b)

            call qgaus(p23hvb_s,dm,z0,p23b)
            
            call qgaus(p12hvd_s,dm,z0,p12d)

            call qgaus(p23hvd_s,dm,z0,p23d)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (+cos(theta_i)**2*p12b*conjg(p12d)
     &       +sin(theta_i)**2*p23b*conjg(p23d))

          else

            print*, 'error in stalk_state_ave j'

          endif

         else if (j .eq. 3) then

           if (k .eq. 1) then

            call qgaus(p12hvc_s,dm,z0,p12c)

            call qgaus(p23hvc_s,dm,z0,p23c)

            call qgaus(p12hva_s,dm,z0,p12a)

            call qgaus(p23hva_s,dm,z0,p23a)
            
            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)**2*p12c*conjg(p12a)
     &       +sin(theta_i)**2*p23c*conjg(p23a))

            
           else if (k .eq. 2) then

            call qgaus(p12hvc_s,dm,z0,p12c)

            call qgaus(p23hvc_s,dm,z0,p23c)
            
            call qgaus(p12hvb_s,dm,z0,p12b)

            call qgaus(p23hvb_s,dm,z0,p23b)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
            *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
            (-cos(theta_i)**2*p12c*conjg(p12b)
            +sin(theta_i)**2*p23c*conjg(p23b))

           else if (k .eq. 3) then

            call qgaus(p12hvc_s,dm,z0,p12c)

            call qgaus(p23hvc_s,dm,z0,p23c)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)**2*p12c*conjg(p12c)
     &       +sin(theta_i)**2*p23c*conjg(p23c))

           else if (k .eq. 4) then

            call qgaus(p12hvc_s,dm,z0,p12c)

            call qgaus(p23hvc_s,dm,z0,p23c)

            call qgaus(p12hvd_s,dm,z0,p12d)

            call qgaus(p23hvd_s,dm,z0,p23d)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)**2*p12c*conjg(p12d)
     &       +sin(theta_i)**2*p23c*conjg(p23d))

          else

            print*, 'error in stalk_state_ave k'

          endif


         else if (j .eq. 4) then

           if (k .eq. 1) then

            call qgaus(p12hvd_s,dm,z0,p12d)

            call qgaus(p23hvd_s,dm,z0,p23d)

            call qgaus(p12hva_s,dm,z0,p12a)

            call qgaus(p23hva_s,dm,z0,p23a)
            
            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)**2*p12d*conjg(p12a)
     &       +sin(theta_i)**2*p23d*conjg(p23a))

            
           else if (k .eq. 2) then

            call qgaus(p12hvd_s,dm,z0,p12d)

            call qgaus(p23hvd_s,dm,z0,p23d)
            
            call qgaus(p12hvb_s,dm,z0,p12b)

            call qgaus(p23hvb_s,dm,z0,p23b)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)**2*p12d*conjg(p12b)
     &       +sin(theta_i)**2*p23d*conjg(p23b))

           else if (k .eq. 3) then

            call qgaus(p12hvd_s,dm,z0,p12d)

            call qgaus(p23hvd_s,dm,z0,p23d)
            
            call qgaus(p12hvc_s,dm,z0,p12c)

            call qgaus(p23hvc_s,dm,z0,p23c)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)**2*p12d*conjg(p12c)
     &       +sin(theta_i)**2*p23d*conjg(p23c))

           else if (k .eq. 4) then

            call qgaus(p12hvd_s,dm,z0,p12d)

            call qgaus(p23hvd_s,dm,z0,p23d)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)**2*p12d*conjg(p12d)
     &       +sin(theta_i)**2*p23d*conjg(p23d))

          else

            print*, 'error in stalk_state_ave l'

          endif

        else

          print*, 'error in stalk_state_ave m'

        endif

      else if (i .eq. 4) then

        if (j .eq. 1) then

          if (k .eq. 1) then

            call qgaus(p22hha_s,dm,z0,p22ha)

            call qgaus(p11vva_s,dm,z0,p11va)

            call qgaus(p33vva_s,dm,z0,p33va)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)**2*p22ha*conjg(p11va)
     &       +sin(theta_i)**2*p22ha*conjg(p33va))

          else if (k .eq. 2) then

            call qgaus(p22hha_s,dm,z0,p22ha)

            call qgaus(p11vvb_s,dm,z0,p11vb)

            call qgaus(p33vvb_s,dm,z0,p33vb)
            
            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)**2*p22ha*conjg(p11vb)
     &       +sin(theta_i)**2*p22ha*conjg(p33vb))

          else if (k .eq. 3) then

            call qgaus(p22hha_s,dm,z0,p22ha)

            call qgaus(p11vvc_s,dm,z0,p11vc)

            call qgaus(p33vvc_s,dm,z0,p33vc)
            
            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)**2*p22ha*conjg(p11vc)
     &       +sin(theta_i)**2*p22ha*conjg(p33vc))

          else if (k .eq. 4) then

            call qgaus(p22hha_s,dm,z0,p22ha)

            call qgaus(p11vvd_s,dm,z0,p11vd)

            call qgaus(p33vvd_s,dm,z0,p33vd)
            
            stalk_state_ave=stalk_alpha_pdf(alpha_s)
      &      *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
      &      (cos(theta_i)**2*p22ha*conjg(p11vd)
      &      +sin(theta_i)**2*p22ha*conjg(p33vd))

          else

            print*, 'error in stalk_state_ave n'

          endif

        else if (j .eq. 2) then

          if (k .eq. 1) then

            call qgaus(p22hhb_s,dm,z0,p22hb)

            call qgaus(p11vva_s,dm,z0,p11va)

            call qgaus(p33vva_s,dm,z0,p33va)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)**2*p22hb*conjg(p11va)
     &       +sin(theta_i)**2*p22hb*conjg(p33va))

          else if (k .eq. 2) then

            call qgaus(p22hhb_s,dm,z0,p22hb)

            call qgaus(p11vvb_s,dm,z0,p11vb)

            call qgaus(p33vvb_s,dm,z0,p33vb)
            
            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)**2*p22hb*conjg(p11vb)
     &       +sin(theta_i)**2*p22hb*conjg(p33vb))

          else if (k .eq. 3) then

            call qgaus(p22hhb_s,dm,z0,p22hb)

            call qgaus(p11vvc_s,dm,z0,p11vc)

            call qgaus(p33vvc_s,dm,z0,p33vc)
            
            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)**2*p22hb*conjg(p11vc)
     &       +sin(theta_i)**2*p22hb*conjg(p33vc))

          else if (k .eq. 4) then

            call qgaus(p22hhb_s,dm,z0,p22hb)

            call qgaus(p11vvd_s,dm,z0,p11vd)

            call qgaus(p33vvd_s,dm,z0,p33vd)
            
            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)**2*p22hb*conjg(p11vd)
     &       +sin(theta_i)**2*p22hb*conjg(p33vd))

          else

            print*, 'error in stalk_state_ave o'

          endif


        else if (j .eq. 3) then

          if (k .eq. 1) then

            call qgaus(p22hhc_s,dm,z0,p22hc)

            call qgaus(p11vva_s,dm,z0,p11va)

            call qgaus(p33vva_s,dm,z0,p33va)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)**2*p22hc*conjg(p11va)
     &       +sin(theta_i)**2*p22hc*conjg(p33va))

          else if (k .eq. 2) then

            call qgaus(p22hhc_s,dm,z0,p22hc)

            call qgaus(p11vvb_s,dm,z0,p11vb)

            call qgaus(p33vvb_s,dm,z0,p33vb)
            
            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)**2*p22hc*conjg(p11vb)
     &       +sin(theta_i)**2*p22hc*conjg(p33vb))

          else if (k .eq. 3) then

            call qgaus(p22hhc_s,dm,z0,p22hc)

            call qgaus(p11vvc_s,dm,z0,p11vc)

            call qgaus(p33vvc_s,dm,z0,p33vc)
            
            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)**2*p22hc*conjg(p11vc)
     &       +sin(theta_i)**2*p22hc*conjg(p33vc))

          else if (k .eq. 4) then

            call qgaus(p22hhc_s,dm,z0,p22hc)

            call qgaus(p11vvd_s,dm,z0,p11vd)

            call qgaus(p33vvd_s,dm,z0,p33vd)
            
            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)**2*p22hc*conjg(p11vd)
     &       +sin(theta_i)**2*p22hc*conjg(p33vd))

          else

            print*, 'error in stalk_state_ave p'

          endif

        else if (j .eq. 4) then

          if (k .eq. 1) then

            call qgaus(p22hhd_s,dm,z0,p22hd)

            call qgaus(p11vva_s,dm,z0,p11va)

            call qgaus(p33vva_s,dm,z0,p33va)

            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)**2*p22hd*conjg(p11va)
     &       +sin(theta_i)**2*p22hd*conjg(p33va))

          else if (k .eq. 2) then

            call qgaus(p22hhd_s,dm,z0,p22hd)

            call qgaus(p11vvb_s,dm,z0,p11vb)

            call qgaus(p33vvb_s,dm,z0,p33vb)
            
            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)**2*p22hd*conjg(p11vb)
     &       +sin(theta_i)**2*p22hd*conjg(p33vb))

          else if (k .eq. 3) then

            call qgaus(p22hhd_s,dm,z0,p22hd)

            call qgaus(p11vvc_s,dm,z0,p11vc)

            call qgaus(p33vvc_s,dm,z0,p33vc)
            
            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)**2*p22hd*conjg(p11vc)
     &       +sin(theta_i)**2*p22hd*conjg(p33vc))

          else if (k .eq. 4) then

            call qgaus(p22hhd_s,dm,z0,p22hd)

            call qgaus(p11vvd_s,dm,z0,p11vd)

            call qgaus(p33vvd_s,dm,z0,p33vd)
            
            stalk_state_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)**2*p22hd*conjg(p11vd)
     &       +sin(theta_i)**2*p22hd*conjg(p33vd))

          else

            print*, 'error in stalk_state_ave q'

          endif

        else

          print*, 'error in stalk_state_ave r'

        endif

      else

        print*, 'error in stalk_state_ave s'

      endif

      return

      end

c************************************************************************
      
      function p11vva_s(z)

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

       complex ii,phia
       complex p11vva_s
       complex p11_s,p22_s,p33_s
       complex tau_v

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)
       
       x=rho*cos(alpha_s)

       phia=cexp(ii*2.0*k_0*(sin(theta_i)*x+cos(theta_i)*z)
     & + 2.0*tau_v(z)/cos(theta_i))

       drho=a

       dl=sqrt(1.0+drho**2)

       cosb=1.0/dl

       sinb=drho/dl

       p11vva_s=phia*(cos(alpha_s)**2*cosb**2*p11_s(z)
     &  +sin(alpha_s)**2*p22_s(z)
     &  +cos(alpha_s)**2*sinb**2*p33_s(z))*dl

       return

       end
      

c************************************************************************
      
      function p11vvb_s(z)

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


       complex ii,phib
       complex p11vvb_s
       complex p11_s,p22_s,p33_s

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       rho=a(d+z)

       x=rho*cos(alpha_s)

       phib=cexp(ii*2.0*k_0*(sin(theta_i)*x))

       drho=a

       dl=sqrt(1.0+drho**2)

       cosb=1.0/dl

       sinb=drho/dl

       p11vvb_s=phib*(cos(alpha_s)**2*cosb**2*p11_s(z)
     &  +sin(alpha_s)**2*p22_s(z)
     &  +cos(alpha_s)**2*sinb**2*p33_s(z))*dl

       return

       end       

c************************************************************************
      
      function p11vvc_s(z)

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

       complex ii,phic
       complex p11vvc_s
       complex p11_s,p22_s,p33_s

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       x=rho*cos(alpha_s)

       phic=cexp(ii*2.0*k_0*(sin(theta_i)*x))

       drho=a

       dl=sqrt(1.0+drho**2)

       cosb=1.0/dl

       sinb=drho/dl

       p11vvc_s=phic*(cos(alpha_s)**2*cosb**2*p11_s(z)
     &  +sin(alpha_s)**2*p22_s(z)
     &  +cos(alpha_s)**2*sinb**2*p33_s(z))*dl

       return

       end              

c************************************************************************
      
      function p11vvd_s(z)

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

       complex ii,phid
       complex p11vvd_s
       complex p11_s,p22_s,p33_s
       complex tau_v

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       x=rho*cos(alpha_s)

       phid=cexp(ii*2.0*k_0*(sin(theta_i)*x-cos(theta_i)*z)
     &  -2.0*tau_v(z)/cos(theta_i))

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

       cosb=1.0/dl

       sinb=drho/dl
       
       p11vvd_s=phid*(cos(alpha_s)**2*cosb**2*p11_s(z)
     &  +sin(alpha_s)**2*p22_s(z)
     &  +cos(alpha_s)**2*sinb**2*p33_s(z))*dl

       return

       end              

c************************************************************************
      
      function p33vva_s(z)

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

       complex ii,phia
       complex p33vva_s
       complex p11_s,p22_s,p33_s
       complex tau_v

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       x=rho*cos(alpha_s)

       phia=cexp(ii*2.0*k_0*(sin(theta_i)*x+cos(theta_i)*z)
     & + 2.0*tau_v(z)/cos(theta_i))

       drho=a

       dl=sqrt(1.0+drho**2)

       cosb=1.0/dl

       sinb=drho/dl
       

       p33vva_s=phia*(sinb**2*p11_s(z)
     &  +cosb**2*p33_s(z))*dl

       return

       end
      

c************************************************************************
      
      function p33vvb_s(z)

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


       complex ii,phib
       complex p33vva_s
       complex p11_s,p22_s,p33_s

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       rho=a(d+z)

       x=rho*cos(alpha_s)

       phib=cexp(ii*2.0*k_0*(sin(theta_i)*x))

       drho=a

       cosb=dz/sqrt(1.0+dz**2)

       sinb=1.0/sqrt(1.0+dz**2)

       p33vvb_s=phib*(sinb**2*p11_s(z)
     &  +cosb**2*p33_s(z))*dl

       return

       end       

c************************************************************************
      
      function p33vvc_s(z)

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

       complex ii,phic
       complex p33vvc_s
       complex p11_s,p22_s,p33_s

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       x=rho*cos(alpha_s)

       phic=cexp(ii*2.0*k_0*(sin(theta_i)*x))

       drho=a

       dl=sqrt(1.0+drho**2)

       cosb=1.0/dl

       sinb=drho/dl       

       p33vvc_s=phic*(sinb**2*p11_s(z)
     &  +cosb**2*p33_s(z))*dl

       return

       end              

c************************************************************************
      
      function p33vvd_s(z)

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

       complex ii,phid
       complex p33vvd_s
       complex p11_s,p22_s,p33_s
       complex tau_v

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       x=rho*cos(alpha_s)

       phid=cexp(ii*2.0*k_0*(sin(theta_i)*x-cos(theta_i)*z)
     &  -2.0*tau_v(z)/cos(theta_i))

       drho=a

       dl=sqrt(1.0+drho**2)

       cosb=1.0/dl

       sinb=drho/dl       

       p33vvd_s=phid*(sinb**2*p11_s(z)
     &  +cosb**2*p33_s(z))*dl

       return

       end              

c************************************************************************
      
      function p13vva_s(z)

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

       complex ii,phia
       complex p13vva_s
       complex p11_s,p22_s,p33_s
       complex tau_v

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       x=rho*cos(alpha_s)

       phia=cexp(ii*2.0*k_0*(sin(theta_i)*x+cos(theta_i)*z)
     & + 2.0*tau_v(z)/cos(theta_i))

       drho=a

       dl=sqrt(1.0+drho**2)

       cosb=1.0/dl

       sinb=drho/dl
       
       p13vva_s=phia*(cos(alpha_s)*cosb*sinb*(p33_s(z)-p11_s(z)))

       return

       end
      

c************************************************************************
      
      function p13vvb_s(z)

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


       complex ii,phib
       complex p13vvb_s
       complex p11_s,p22_s,p33_s

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       x=rho*cos(alpha_s)

       phib=cexp(ii*2.0*k_0*(sin(theta_i)*x))

       drho=a

       cosb=dz/sqrt(1.0+dz**2)

       sinb=1.0/sqrt(1.0+dz**2)

       p13vvb_s=phib*(cos(alpha_s)*cosb*sinb*(p33_s(z)-p11_s(z)))

       return

       end       

c************************************************************************
      
      function p13vvc_s(z)

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

       complex ii,phic
       complex p13vvc_s
       complex p11_s,p22_s,p33_s

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       x=rho*cos(alpha_s)

       phic=cexp(ii*2.0*k_0*(sin(theta_i)*x))

       drho=a

       dl=sqrt(1.0+drho**2)

       cosb=1.0/dl

       sinb=drho/dl       

       p13vvc_s=phic*(cos(alpha_s)*cosb*sinb*(p33_s(z)-p11_s(z)))

       return

       end              

c************************************************************************
      
      function p13vvd_s(z)

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

       complex ii,phid
       complex p11vvd_s
       complex p11_s,p22_s,p33_s
       complex tau_v

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       x=rho*cos(alpha_s)

       phid=cexp(ii*2.0*k_0*(sin(theta_i)*x-cos(theta_i)*z)
     &  -2.0*tau_v(z)/cos(theta_i))

       drho=a

       dl=sqrt(1.0+drho**2)

       cosb=1.0/dl

       sinb=drho/dl       

       p13vvd_s=phid*(cos(alpha_s)*cosb*sinb*(p33_s(z)-p11_s(z)))

       return

       end              

c************************************************************************
      
      function p22hha_s(z)

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

       complex ii,phia
       complex p22hha_s
       complex p11_s,p22_s,p33_s
       complex tau_h

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       x=rho*cos(alpha_s)

       phia=cexp(ii*2.0*k_0*(sin(theta_i)*x+cos(theta_i)*z)
     & + 2.0*tau_h(z)/cos(theta_i))

       drho=a

       dl=sqrt(1.0+drho**2)

       cosb=1.0/dl

       sinb=drho/dl

       p22hha_s=phia*(sin(alpha_s)**2*cosb**2*p11_s(z)
     &  +cos(alpha_s)**2*p22_s(z)
     &  +sin(alpha_s)**2*sinb**2*p33_s(z))*dl

       return

       end
      

c************************************************************************
      
      function p22hhb_s(z)

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

       complex ii,phib
       complex p22hhb_s
       complex p11_s,p22_s,p33_s

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       x=rho*cos(alpha_s)

       phib=cexp(ii*2.0*k_0*(sin(theta_i)*x))

       drho=a

       dl=sqrt(1.0+drho**2)

       cosb=1.0/dl

       sinb=drho/dl

       p22hhb_s=phib*(sin(alpha_s)**2*cosb**2*p11_s(z)
     &  +cos(alpha_s)**2*p22_s(z)
     &  +sin(alpha_s)**2*sinb**2*p33_s(z))*dl

       return

       end       

c************************************************************************
      
      function p22hhc_s(z)

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

       complex ii,phic
       complex p22hhc_s
       complex p11_s,p22_s,p33_s

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       x=rho*cos(alpha_s)

       phic=cexp(ii*2.0*k_0*(sin(theta_i)*x))

       drho=a

       dl=sqrt(1.0+drho**2)

       cosb=1.0/dl

       sinb=drho/dl

       p22hhc=phic*(sin(alpha_s)**2*cosb**2*p11_s(z)
     &  +cos(alpha_s)**2*p22_s(z)
     &  +sin(alpha_s)**2*sinb**2*p33_s(z))*dl

       return

       end              

c************************************************************************
      
      function p22hhd_s(z)

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

       complex ii,phid
       complex p22hhd_s
       complex p11_s,p22_s,p33_s
       complex tau_h

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       x=rho*cos(alpha_s)

       phid=cexp(ii*2.0*k_0*(sin(theta_i)*x-cos(theta_i)*z)
     &  -2.0*tau_h(z)/cos(theta_i))

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

       cosb=1.0/dl

       sinb=drho/dl
       
       p11vvd_s=phid*(sin(alpha_s)**2*cosb**2*p11_s(z)
     &  +cos(alpha_s)**2*p22_s(z)
     &  +sin(alpha_s)**2*sinb**2*p33_s(z))*dl
  
       return

       end              

c************************************************************************
      
      function p12hva_s(z)

c************************************************************************
       complex ii,phia
       complex p12hva_s
       complex p11_s,p22_s,p33_s
       complex tau_v,tau_h

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       x=rho*cos(alpha_s)

       phia=cexp(ii*2.0*k_0*(sin(theta_i)*x+cos(theta_i)*z)
     & + (tau_v(z)+tau_h(z))/cos(theta_i))

       drho=a

       dl=sqrt(1.0+drho**2)

       cosb=1.0/dl

       sinb=drho/dl

       p12hva_s=phia*(sin(alpha_s)*cos(alpha_s)*cosb**2*p11_s(z)
     &  -cos(alpha_s)*sin(alpha_s)*p22_s(z)
     &  +sin(alpha_s)*cos(alpha_s)*sinb**2*p33_s(z))*dl

       return

       end
      

c************************************************************************
      
      function p12hvb_s(z)

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

       complex ii,phib
       complex p12hvb_s
       complex p11_s,p22_s,p33_s
       complex tau_v,tau_h

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       x=rho*cos(alpha_s)

       phib=cexp(ii*2.0*k_0*(sin(theta_i)*x)
     &  +(tau_h(z)-tau_v(z))/cos(theta_i))

       drho=a

       dl=sqrt(1.0+drho**2)

       cosb=1.0/dl

       sinb=drho/dl

       p12hvb_s=phib*(sin(alpha_s)*cos(alpha_s)*cosb**2*p11_s(z)
     &  -cos(alpha_s)*sin(alpha_s)*p22_s(z)
     &  +sin(alpha_s)*cos(alpha_s)*sinb**2*p33_s(z))*dl

       return

       end       

c************************************************************************
      
      function p12hvc_s(z)

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

       complex ii,phic
       complex p12hvc_s
       complex p11_s,p22_s,p33_s
       complex tau_v,tau_h

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       x=rho*cos(alpha_s)

       phic=cexp(ii*2.0*k_0*(sin(theta_i)*x))
     &  +(tau_v(z)-tau_h(z))/cos(theta_i))
       
       drho=a

       dl=sqrt(1.0+drho**2)

       cosb=1.0/dl

       sinb=drho/dl

       p12hvc_s=phic*(sin(alpha_s)*cos(alpha_s)*cosb**2*p11_s(z)
     &  -cos(alpha_s)*sin(alpha_s)*p22_s(z)
     &  +sin(alpha_s)*cos(alpha_s)*sinb**2*p33_s(z))*dl
       
       return

       end              

c************************************************************************
      
      function p12hvd_s(z)

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


       complex ii,phid
       complex p12hvd_s
       complex p11_s,p22_s,p33_s
       complex tau_v,tau_h

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       x=rho*cos(alpha_s)

       phid=cexp(ii*2.0*k_0*(sin(theta_i)*x-cos(theta_i)*z)
     &  -(tau_v(z)+tau_h(z))/cos(theta_i))

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

       cosb=1.0/dl

       sinb=drho/dl
       

       p12hvd_s=phid*(sin(alpha_s)*cos(alpha_s)*cosb**2*p11_s(z)
     &  -cos(alpha_s)*sin(alpha_s)*p22_s(z)
     & +sin(alpha_s)*cos(alpha_s)*sinb**2*p33_s(z))*dl

       return

       end

      
c************************************************************************
      
      function p23hva_s(z)

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

       complex ii,phia
       complex p23hva_s
       complex p11_s,p22_s,p33_s
       complex tau_v,tau_h

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       x=rho*cos(alpha_s)

       phia=cexp(ii*2.0*k_0*(sin(theta_i)*x+cos(theta_i)*z)
     & + (tau_v(z)+tau_h(z))/cos(theta_i))

       drho=a

       dl=sqrt(1.0+drho**2)

       cosb=1.0/dl

       sinb=drho/dl

       p23hva_s=phia*(sin(alpha_s)*sinb*cosb*(p33_s(z)-p11_s(z)))*dl

       return

       end
      

c************************************************************************
      
      function p23hvb_s(z)

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

       complex ii,phib
       complex p23hvb_s
       complex p11_s,p22_s,p33_s
       complex tau_v,tau_h

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)
       
       x=rho*cos(alpha_s)

       phib=cexp(ii*2.0*k_0*(sin(theta_i)*x)
     &  +(tau_h(z)-tau_v(z))/cos(theta_i))

       drho=a

       dl=sqrt(1.0+drho**2)

       cosb=1.0/dl

       sinb=drho/dl

       p23hvb_s=phib*(sin(alpha_s)*sinb*cosb*(p33_s(z)-p11_s(z)))*dl

       return

       end       

c************************************************************************
      
      function p23hvc_s(z)

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

       complex ii,phic
       complex p23hvc_s
       complex p11_s,p22_s,p33_s
       complex tau_v,tau_h

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       x=rho*cos(alpha_s)

       phic=cexp(ii*2.0*k_0*(sin(theta_i)*x))
     &  +(tau_v(z)-tau_h(z))/cos(theta_i))
       
       drho=a

       dl=sqrt(1.0+drho**2)

       cosb=1.0/dl

       sinb=drho/dl

       p23hvc_s=phic*(sin(alpha_s)*sinb*cosb*(p33_s(z)-p11_s(z)))*dl

       return

       end              

c************************************************************************
      
      function p23hvd_s(z)

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

       complex ii,phid
       complex p23hvd_s
       complex p11_s,p22_s,p33_s
       complex tau_v,tau_h

       real rho,x,z
       real dz,dl,cosb,sinb,drho
       real a,alpha_s
       real theta_i,phi_i,k_0
       real dm,del_leaf,z_lc
       
       common /stalk_struct/ a,alpha_s
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ dm,del_leaf,z_lc
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z-dm)

       x=rho*cos(alpha_s)

       phid=cexp(ii*2.0*k_0*(sin(theta_i)*x-cos(theta_i)*z)
     &  -(tau_v(z)+tau_h(z))/cos(theta_i))

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

       cosb=1.0/dl

       sinb=drho/dl

       p23hvd_s=phid*(sin(alpha_s)*sinb*cosb*(p33_s(z)-p11_s(z)))*dl

       return

       end              

