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

      complex function stalk_scat(phi,aa,z0,pol,mech,mech_c)

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

      integer pol,mech,mech_c
      

      integer hhhh,vvvv,hvhv,hhvv
      parameter(vvvv=1)
      parameter(hhhh=2)
      parameter(hvhv=3)
      parameter(hhvv=4)

      integer p,gp,pg,gpg
      parameter(p=1)
      parameter(gp=2)
      parameter(pg=3)
      parameter(gpg=4)
      
      complex pxxvva_s,pzzvva_s,pxzvva_s
      complex pxxvvb_s,pzzvvb_s,pxzvvb_s
      complex pxxvvc_s,pzzvvc_s,pxzvvc_s
      complex pxxvvd_s,pzzvvd_s,pxzvvd_s

      complex pyyhha_s
      complex pyyhhb_s
      complex pyyhhc_s
      complex pyyhhd_s

      complex pxyhva_s,pyzhva_s
      complex pxyhvb_s,pyzhvb_s
      complex pxyhvc_s,pyzhvc_s
      complex pxyhvd_s,pyzhvd_s

      complex pxxa,pxya,pxza,pyya,pyza,pzza
      complex pxxb,pxyb,pxzb,pyyb,pyzb,pzzb
      complex pxxc,pxyc,pxzc,pyyc,pyzc,pzzc
      complex pxxd,pxyd,pxzd,pyyd,pyzd,pzzd

      complex pxxa_mn,pxya_mn,pxza_mn,pyya_mn,pyza_mn,pzza_mn
      complex pxxb_mn,pxyb_mn,pxzb_mn,pyyb_mn,pyzb_mn,pzzb_mn
      complex pxxc_mn,pxyc_mn,pxzc_mn,pyyc_mn,pyzc_mn,pzzc_mn
      complex pxxd_mn,pxyd_mn,pxzd_mn,pyyd_mn,pyzd_mn,pzzd_mn

      real z0
      real a,alpha_s,ell_0
      real theta_i,phi_i,k_0
      real phi,aa
      real d,del_leaf,z_lc,hole
      real pi,coef
      real dl,drho
      real z0_mn,z0_sig
      
      external pxxvva_s,pzzvva_s,pxzvva_s
      external pxxvvb_s,pzzvvb_s,pxzvvb_s
      external pxxvvc_s,pzzvvc_s,pxzvvc_s
      external pxxvvd_s,pzzvvd_s,pxzvvd_s

      external pyyhha_s
      external pyyhhb_s
      external pyyhhc_s
      external pyyhhd_s

      external pxyhva_s,pyzhva_s
      external pxyhvb_s,pyzhvb_s
      external pxyhvc_s,pyzhvc_s
      external pxyhvd_s,pyzhvd_s
      
      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
      common /layer_struct/ d,del_leaf,z_lc,hole
      common /z0_stat/ z0_mn,z0_sig

      common /p_store/
     &     pxxa_mn,pxya_mn,pxza_mn,pyya_mn,pyza_mn,pzza_mn,     
     &     pxxb_mn,pxyb_mn,pxzb_mn,pyyb_mn,pyzb_mn,pzzb_mn,
     &     pxxc_mn,pxyc_mn,pxzc_mn,pyyc_mn,pyzc_mn,pzzc_mn,
     &     pxxd_mn,pxyd_mn,pxzd_mn,pyyd_mn,pyzd_mn,pzzd_mn

      save /p_store/

      include 'error.incld'
      
      alpha_s=phi
      a=aa

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

c     this has been changed to allow for efficient integration
c      from -d to z0, ell_0 is now not the length of the stalk,
c      but is instead is greater than or equal to the MAXIMUM
c      length of the stalk.
c     
c      ell_0=dl*(z+d) <- old ell_0
c
c     new ell_0 follows:

      ell_0= dl*d
      
      pi=3.1415927
      coef=k_0**4/(16.0*pi**2)
      s_tensor_flag=0

      
      if (pol .eq. vvvv) then

        if (mech .eq. p) then

          if (mech_c .eq. p) then

            if(z0_flag .eq. 0) then

              call integral(pxxvva_s,-d,z0_mn,pxxa_mn)
              call integral(pzzvva_s,-d,z0_mn,pzza_mn)
              call integral(pxzvva_s,-d,z0_mn,pxza_mn)

              z0_flag=1

          else
              
              continue
              
            endif
            
              call integral(pxxvva_s,z0_mn,z0,pxxa)
              call integral(pzzvva_s,z0_mn,z0,pzza)
              call integral(pxzvva_s,z0_mn,z0,pxza)

              pxxa=pxxa+pxxa_mn
              pzza=pzza+pzza_mn
              pxza=pxza+pxza_mn            
            
            
            stalk_scat=coef* 
     &       (cos(theta_i)**4*pxxa*conjg(pxxa)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (4.0*pxza*conjg(pxza)
     &       +pxxa*conjg(pzza)+pzza*conjg(pxxa))
     &       +sin(theta_i)**4*pzza*conjg(pzza))


          else if (mech_c .eq. gp) then

            if(z0_flag .eq. 0) then

              call integral(pxxvva_s,-d,z0_mn,pxxa_mn)
              call integral(pzzvva_s,-d,z0_mn,pzza_mn)
              call integral(pxxvvb_s,-d,z0_mn,pxxb_mn)
              call integral(pzzvvb_s,-d,z0_mn,pzzb_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
              call integral(pxxvva_s,z0_mn,z0,pxxa)
              call integral(pzzvva_s,z0_mn,z0,pzza)
              call integral(pxxvvb_s,z0_mn,z0,pxxb)
              call integral(pzzvvb_s,z0_mn,z0,pzzb)

              pxxa=pxxa+pxxa_mn
              pzza=pzza+pzza_mn
              pxxb=pxxb+pxxb_mn
              pzzb=pzzb+pzzb_mn

            stalk_scat=coef*
     &       (-cos(theta_i)**4*pxxa*conjg(pxxb)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (pxxa*conjg(pzzb)-pzza*conjg(pxxb))
     &             +sin(theta_i)**4*pzza*conjg(pzzb))
     
          else if (mech_c .eq. pg) then
            
            if(z0_flag .eq. 0) then

              call integral(pxxvva_s,-d,z0_mn,pxxa_mn)
              call integral(pzzvva_s,-d,z0_mn,pzza_mn)
              call integral(pxxvvc_s,-d,z0_mn,pxxc_mn)
              call integral(pzzvvc_s,-d,z0_mn,pzzc_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
              call integral(pxxvva_s,z0_mn,z0,pxxa)
              call integral(pzzvva_s,z0_mn,z0,pzza)
              call integral(pxxvvc_s,z0_mn,z0,pxxc)
              call integral(pzzvvc_s,z0_mn,z0,pzzc)

              pxxa=pxxa+pxxa_mn
              pzza=pzza+pzza_mn
              pxxc=pxxc+pxxc_mn
              pzzc=pzzc+pzzc_mn

            stalk_scat=coef*
     &      (-cos(theta_i)**4*pxxa*conjg(pxxc)
     &      +cos(theta_i)**2*sin(theta_i)**2*
     &      (pxxa*conjg(pzzc)-pzza*conjg(pxxc))
     &      +sin(theta_i)**4*pzza*conjg(pzzc))

          else if (mech_c .eq. gpg) then
            
            if(z0_flag .eq. 0) then

              call integral(pxxvva_s,-d,z0_mn,pxxa_mn)
              call integral(pzzvva_s,-d,z0_mn,pzza_mn)
              call integral(pxzvva_s,-d,z0_mn,pxza_mn)
              call integral(pxxvvd_s,-d,z0_mn,pxxd_mn)
              call integral(pzzvvd_s,-d,z0_mn,pzzd_mn)
              call integral(pxzvvd_s,-d,z0_mn,pxzd_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
              call integral(pxxvva_s,z0_mn,z0,pxxa)
              call integral(pzzvva_s,z0_mn,z0,pzza)
              call integral(pxzvva_s,z0_mn,z0,pxza)
              call integral(pxxvvd_s,z0_mn,z0,pxxd)
              call integral(pzzvvd_s,z0_mn,z0,pzzd)
              call integral(pxzvvd_s,z0_mn,z0,pxzd)
              
              pxxa=pxxa+pxxa_mn
              pzza=pzza+pzza_mn
              pxza=pxza+pxza_mn
              pxxd=pxxd+pxxd_mn
              pzzd=pzzd+pzzd_mn
              pxzd=pxzd+pxzd_mn
              
            stalk_scat=coef*
     &       (cos(theta_i)**4*pxxa*conjg(pxxd)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (-4.0*pxza*conjg(pxzd)+
     &       pxxa*conjg(pzzd)+pzza*conjg(pxxd))
     &       +sin(theta_i)**4*pzza*conjg(pzzd))
                      
          else
            
            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func leaf_scat'
            code=43
            call the_death_code(die,error,code)
            
          endif

        else if (mech .eq. gp) then
          
          if (mech_c .eq. p) then

            if(z0_flag .eq. 0) then
              
              call integral(pxxvva_s,-d,z0_mn,pxxa_mn)
              call integral(pzzvva_s,-d,z0_mn,pzza_mn)
              call integral(pxxvvb_s,-d,z0_mn,pxxb_mn)
              call integral(pzzvvb_s,-d,z0_mn,pzzb_mn)
              
            z0_flag=1
            
          else
              
              continue
              
            endif
            
              call integral(pxxvva_s,z0_mn,z0,pxxa)
              call integral(pzzvva_s,z0_mn,z0,pzza)
              call integral(pxxvvb_s,z0_mn,z0,pxxb)
              call integral(pzzvvb_s,z0_mn,z0,pzzb)

              pxxa=pxxa+pxxa_mn
              pzza=pzza+pzza_mn
              pxxb=pxxb+pxxb_mn
              pzzb=pzzb+pzzb_mn
            

            stalk_scat=coef*
     &       (-cos(theta_i)**4*pxxb*conjg(pxxa)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (pzzb*conjg(pxxa)-pxxb*conjg(pzza))
     &       +sin(theta_i)**4*pzzb*conjg(pzza))            

          else if (mech_c .eq. gp) then

            if(z0_flag .eq. 0) then
              
              call integral(pxxvvb_s,-d,z0_mn,pxxb_mn)
              call integral(pzzvvb_s,-d,z0_mn,pzzb_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
              call integral(pxxvvb_s,z0_mn,z0,pxxb)
              call integral(pzzvvb_s,z0_mn,z0,pzzb)

              pxxb=pxxb+pxxb_mn
              pzzb=pzzb+pzzb_mn

            stalk_scat=coef*
     &       (cos(theta_i)**4*pxxb*conjg(pxxb)
     &       -cos(theta_i)**2*sin(theta_i)**2*
     &       (pzzb*conjg(pxxb)+pxxb*conjg(pzzb))
     &       +sin(theta_i)**4*pzzb*conjg(pzzb))

          else if (mech_c .eq. pg) then

            if(z0_flag .eq. 0) then
              
              call integral(pxxvvb_s,-d,z0_mn,pxxb_mn)
              call integral(pzzvvb_s,-d,z0_mn,pzzb_mn)
              call integral(pxxvvc_s,-d,z0_mn,pxxc_mn)
              call integral(pzzvvc_s,-d,z0_mn,pzzc_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
              call integral(pxxvvb_s,z0_mn,z0,pxxb)
              call integral(pzzvvb_s,z0_mn,z0,pzzb)
              call integral(pxxvvc_s,z0_mn,z0,pxxc)
              call integral(pzzvvc_s,z0_mn,z0,pzzc)

              pxxb=pxxb+pxxb_mn
              pzzb=pzzb+pzzb_mn
              pxxc=pxxc+pxxc_mn
              pzzc=pzzc+pzzc_mn

            stalk_scat=coef*
     &       (cos(theta_i)**4*pxxb*conjg(pxxc)
     &       -cos(theta_i)**2*sin(theta_i)**2*
     &       (pzzb*conjg(pxxc)+pxxb*conjg(pzzc))
     &       +sin(theta_i)**4*pzzb*conjg(pzzc))            

          else if (mech_c .eq. gpg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvvb_s,-d,z0_mn,pxxb_mn)
              call integral(pzzvvb_s,-d,z0_mn,pzzb_mn)
              call integral(pxxvvd_s,-d,z0_mn,pxxd_mn)
              call integral(pzzvvd_s,-d,z0_mn,pzzd_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
              call integral(pxxvvb_s,z0_mn,z0,pxxb)
              call integral(pzzvvb_s,z0_mn,z0,pzzb)
              call integral(pxxvvd_s,z0_mn,z0,pxxd)
              call integral(pzzvvd_s,z0_mn,z0,pzzd)

              pxxb=pxxb+pxxb_mn
              pzzb=pzzb+pzzb_mn
              pxxd=pxxd+pxxd_mn
              pzzd=pzzd+pzzd_mn

            stalk_scat=coef*
     &       (-cos(theta_i)**4*pxxb*conjg(pxxd)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (pzzb*conjg(pxxd)-pxxb*conjg(pzzd))
     &       +sin(theta_i)**4*pzzb*conjg(pzzd))

          else
            
            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func leaf_scat'
            code=43
            call the_death_code(die,error,code)
            
          endif

        else if(mech .eq. pg) then

          if (mech_c .eq. p) then
            
            if(z0_flag .eq. 0) then

              call integral(pxxvva_s,-d,z0_mn,pxxa_mn)
              call integral(pzzvva_s,-d,z0_mn,pzza_mn)
              call integral(pxxvvc_s,-d,z0_mn,pxxc_mn)
              call integral(pzzvvc_s,-d,z0_mn,pzzc_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
              call integral(pxxvva_s,z0_mn,z0,pxxa)
              call integral(pzzvva_s,z0_mn,z0,pzza)
              call integral(pxxvvc_s,z0_mn,z0,pxxc)
              call integral(pzzvvc_s,z0_mn,z0,pzzc)

              pxxa=pxxa+pxxa_mn
              pzza=pzza+pzza_mn
              pxxc=pxxc+pxxc_mn
              pzzc=pzzc+pzzc_mn


            stalk_scat=coef*
     &       (-cos(theta_i)**4*pxxc*conjg(pxxa)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (pzzc*conjg(pxxa)-pxxc*conjg(pzza))
     &       +sin(theta_i)**4*pzzc*conjg(pzza))            

          else if (mech_c .eq. gp) then

            if(z0_flag .eq. 0) then
              
              call integral(pxxvvb_s,-d,z0_mn,pxxb_mn)
              call integral(pzzvvb_s,-d,z0_mn,pzzb_mn)
              call integral(pxxvvc_s,-d,z0_mn,pxxc_mn)
              call integral(pzzvvc_s,-d,z0_mn,pzzc_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
              call integral(pxxvvb_s,z0_mn,z0,pxxb)
              call integral(pzzvvb_s,z0_mn,z0,pzzb)
              call integral(pxxvvc_s,z0_mn,z0,pxxc)
              call integral(pzzvvc_s,z0_mn,z0,pzzc)

              pxxb=pxxb+pxxb_mn
              pzzb=pzzb+pzzb_mn
              pxxc=pxxc+pxxc_mn
              pzzc=pzzc+pzzc_mn

            stalk_scat=coef*
     &       (cos(theta_i)**4*pxxc*conjg(pxxb)
     &       -cos(theta_i)**2*sin(theta_i)**2*
     &       (pzzc*conjg(pxxb)+pxxc*conjg(pzzb))
     &       +sin(theta_i)**4*pzzc*conjg(pzzb))

          else if (mech_c .eq. pg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvvc_s,-d,z0_mn,pxxc_mn)
              call integral(pzzvvc_s,-d,z0_mn,pzzc_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
              call integral(pxxvvc_s,z0_mn,z0,pxxc)
              call integral(pzzvvc_s,z0_mn,z0,pzzc)

              pxxc=pxxc+pxxc_mn
              pzzc=pzzc+pzzc_mn
              

            stalk_scat=coef*
     &       (cos(theta_i)**4*pxxc*conjg(pxxc)
     &       -cos(theta_i)**2*sin(theta_i)**2*
     &       (pzzc*conjg(pxxc)+pxxc*conjg(pzzc))
     &       +sin(theta_i)**4*pzzc*conjg(pzzc))            

          else if (mech_c .eq. gpg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvvd_s,-d,z0_mn,pxxd_mn)
              call integral(pzzvvd_s,-d,z0_mn,pzzd_mn)
              call integral(pxxvvc_s,-d,z0_mn,pxxc_mn)
              call integral(pzzvvc_s,-d,z0_mn,pzzc_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
              call integral(pxxvvd_s,z0_mn,z0,pxxd)
              call integral(pzzvvd_s,z0_mn,z0,pzzd)
              call integral(pxxvvc_s,z0_mn,z0,pxxc)
              call integral(pzzvvc_s,z0_mn,z0,pzzc)

              pxxd=pxxd+pxxd_mn
              pzzd=pzzd+pzzd_mn
              pxxc=pxxc+pxxc_mn
              pzzc=pzzc+pzzc_mn
              
            stalk_scat=coef*
     &       (-cos(theta_i)**4*pxxc*conjg(pxxd)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (pzzc*conjg(pxxd)-pxxc*conjg(pzzd))
     &       +sin(theta_i)**4*pzzc*conjg(pzzd))

          else

            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func leaf_scat'
            code=43
            call the_death_code(die,error,code)

          endif

        else if (mech .eq. gpg) then

          if (mech_c .eq. p) then
            
            if(z0_flag .eq. 0) then

              call integral(pxxvva_s,-d,z0_mn,pxxa_mn)
              call integral(pzzvva_s,-d,z0_mn,pzza_mn)
              call integral(pxzvva_s,-d,z0_mn,pxza_mn)
              call integral(pxxvvd_s,-d,z0_mn,pxxd_mn)
              call integral(pzzvvd_s,-d,z0_mn,pzzd_mn)
              call integral(pxzvvd_s,-d,z0_mn,pxzd_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
              call integral(pxxvva_s,z0_mn,z0,pxxa)
              call integral(pzzvva_s,z0_mn,z0,pzza)
              call integral(pxzvva_s,z0_mn,z0,pxza)
              call integral(pxxvvd_s,z0_mn,z0,pxxd)
              call integral(pzzvvd_s,z0_mn,z0,pzzd)
              call integral(pxzvvd_s,z0_mn,z0,pxzd)
              
              pxxa=pxxa+pxxa_mn
              pzza=pzza+pzza_mn
              pxza=pxza+pxza_mn
              pxxd=pxxd+pxxd_mn
              pzzd=pzzd+pzzd_mn
              pxzd=pxzd+pxzd_mn
              
            stalk_scat=coef*
     &       (cos(theta_i)**4*pxxd*conjg(pxxa)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (-4.0*pxzd*conjg(pxza)
     &       +pxxd*conjg(pzza)+pzzd*conjg(pxxa))
     &       +sin(theta_i)**4*pzzd*conjg(pzza))

          else if (mech_c .eq. gp) then

            if(z0_flag .eq. 0) then
              
              call integral(pxxvvb_s,-d,z0_mn,pxxb_mn)
              call integral(pzzvvb_s,-d,z0_mn,pzzb_mn)
              call integral(pxxvvd_s,-d,z0_mn,pxxd_mn)
              call integral(pzzvvd_s,-d,z0_mn,pzzd_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
              call integral(pxxvvb_s,z0_mn,z0,pxxb)
              call integral(pzzvvb_s,z0_mn,z0,pzzb)
              call integral(pxxvvd_s,z0_mn,z0,pxxd)
              call integral(pzzvvd_s,z0_mn,z0,pzzd)

              pxxb=pxxb+pxxb_mn
              pzzb=pzzb+pzzb_mn
              pxxd=pxxd+pxxd_mn
              pzzd=pzzd+pzzd_mn

            stalk_scat=coef*
     &       (-cos(theta_i)**4*pxxd*conjg(pxxb)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (pxxd*conjg(pzzb)-pzzd*conjg(pxxb))
     &       +sin(theta_i)**4*pzzd*conjg(pzzb))
     
          else if (mech_c .eq. pg) then

            if(z0_flag .eq. 0) then
              
              call integral(pxxvvd_s,-d,z0_mn,pxxd_mn)
              call integral(pzzvvd_s,-d,z0_mn,pzzd_mn)
              call integral(pxxvvc_s,-d,z0_mn,pxxc_mn)
              call integral(pzzvvc_s,-d,z0_mn,pzzc_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
              call integral(pxxvvd_s,z0_mn,z0,pxxd)
              call integral(pzzvvd_s,z0_mn,z0,pzzd)
              call integral(pxxvvc_s,z0_mn,z0,pxxc)
              call integral(pzzvvc_s,z0_mn,z0,pzzc)

              pxxd=pxxd+pxxd_mn
              pzzd=pzzd+pzzd_mn
              pxxc=pxxc+pxxc_mn
              pzzc=pzzc+pzzc_mn

            stalk_scat=coef*
     &       (-cos(theta_i)**4*pxxd*conjg(pxxc)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (pxxd*conjg(pzzc)-pzzd*conjg(pxxc))
     &       +sin(theta_i)**4*pzzd*conjg(pzzc))

          else if (mech_c .eq. gpg) then

            if(z0_flag .eq. 0) then

              call integral(pxxvvd_s,-d,z0_mn,pxxd_mn)
              call integral(pzzvvd_s,-d,z0_mn,pzzd_mn)
              call integral(pxzvvd_s,-d,z0_mn,pxzd_mn)

            z0_flag=1

          else
              
              continue
              
            endif
            
              call integral(pxxvvd_s,z0_mn,z0,pxxd)
              call integral(pzzvvd_s,z0_mn,z0,pzzd)
              call integral(pxzvvd_s,z0_mn,z0,pxzd)

              pxxd=pxxd+pxxd_mn
              pzzd=pzzd+pzzd_mn
              pxzd=pxzd+pxzd_mn

            stalk_scat=coef*
     &       (cos(theta_i)**4*pxxd*conjg(pxxd)
     &       +cos(theta_i)**2*sin(theta_i)**2*
     &       (4.0*pxzd*conjg(pxzd)+
     &       pxxd*conjg(pzzd)+pzzd*conjg(pxxd))
     &       +sin(theta_i)**4*pzzd*conjg(pzzd))
                      
          else

            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func leaf_scat'
            code=43
            call the_death_code(die,error,code)
            
          endif
          
        else

          write(error,*)'mech = ',mech,' not recognized',
     &         ' in func leaf_scat'
          code=44
          call the_death_code(die,error,code)
            
        endif

      else if (pol .eq. hhhh) then

        if (mech .eq. p) then

          if (mech_c .eq. p) then
            
           if(z0_flag .eq. 0) then
               
              call integral(pyyhha_s,-d,z0_mn,pyya_mn)

            z0_flag=1

          else
              
              continue
              
            endif
            
            call integral(pyyhha_s,z0_mn,z0,pyya)
            
            pyya=pyya+pyya_mn

            stalk_scat=coef*pyya*conjg(pyya)
            
          else if (mech_c .eq. gp) then
            
            if(z0_flag .eq. 0) then
               
              call integral(pyyhha_s,-d,z0_mn,pyya_mn)
              call integral(pyyhhb_s,-d,z0_mn,pyyb_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
            call integral(pyyhha_s,z0_mn,z0,pyya)
            call integral(pyyhhb_s,z0_mn,z0,pyyb)
            
            pyya=pyya+pyya_mn
            pyyb=pyyb+pyyb_mn
            
            stalk_scat=coef*pyya*conjg(pyyb)
         
         else if (mech_c .eq.pg) then

            if(z0_flag .eq. 0) then
               
              call integral(pyyhha_s,-d,z0_mn,pyya_mn)
              call integral(pyyhhc_s,-d,z0_mn,pyyc_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
            call integral(pyyhha_s,z0_mn,z0,pyya)
            call integral(pyyhhc_s,z0_mn,z0,pyyc)
            
            pyya=pyya+pyya_mn
            pyyc=pyyc+pyyc_mn

            stalk_scat=coef*pyya*conjg(pyyc)
                     
         else if (mech_c .eq.gpg) then
           
            if(z0_flag .eq. 0) then
               
              call integral(pyyhha_s,-d,z0_mn,pyya_mn)
              call integral(pyyhhd_s,-d,z0_mn,pyyd_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
            call integral(pyyhha_s,z0_mn,z0,pyya)
            call integral(pyyhhd_s,z0_mn,z0,pyyd)
            
            pyya=pyya+pyya_mn
            pyyd=pyyd+pyyd_mn
            

            stalk_scat=coef*pyya*conjg(pyyd)            
        
         else

            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func leaf_scat'
            code=43
            call the_death_code(die,error,code)
            
         endif

       else if (mech .eq.gp) then

         if (mech_c .eq. p) then

           if(z0_flag .eq. 0) then
               
              call integral(pyyhha_s,-d,z0_mn,pyya_mn)
              call integral(pyyhhb_s,-d,z0_mn,pyyb_mn)

              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhha_s,z0_mn,z0,pyya)
            call integral(pyyhhb_s,z0_mn,z0,pyyb)
            
            pyya=pyya+pyya_mn
            pyyb=pyyb+pyyb_mn


            stalk_scat=coef*pyyb*conjg(pyya)

          else if (mech_c .eq. gp) then

           if(z0_flag .eq. 0) then
               
              call integral(pyyhhb_s,-d,z0_mn,pyyb_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
            call integral(pyyhhb_s,z0_mn,z0,pyyb)
            
            pyyb=pyyb+pyyb_mn

            stalk_scat=coef*pyyb*conjg(pyyb)

            
         else if (mech_c .eq.pg) then

           if(z0_flag .eq. 0) then
               
              call integral(pyyhhc_s,-d,z0_mn,pyyc_mn)
              call integral(pyyhhb_s,-d,z0_mn,pyyb_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
            call integral(pyyhhc_s,z0_mn,z0,pyyc)
            call integral(pyyhhb_s,z0_mn,z0,pyyb)
            
            pyyc=pyyc+pyyc_mn
            pyyb=pyyb+pyyb_mn

            stalk_scat=coef*pyyb*conjg(pyyc)
                     
         else if (mech_c .eq.gpg) then

           
           if(z0_flag .eq. 0) then
               
              call integral(pyyhhd_s,-d,z0_mn,pyyd_mn)
              call integral(pyyhhb_s,-d,z0_mn,pyyb_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
            call integral(pyyhhd_s,z0_mn,z0,pyyd)
            call integral(pyyhhb_s,z0_mn,z0,pyyb)
            
            pyyd=pyyd+pyyd_mn
            pyyb=pyyb+pyyb_mn

            stalk_scat=coef*pyyb*conjg(pyyd)            
        
         else

            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func leaf_scat'
            code=43
            call the_death_code(die,error,code)
            
         endif

       else if (mech .eq. pg) then

         if (mech_c .eq. p) then
           
            if(z0_flag .eq. 0) then
               
              call integral(pyyhha_s,-d,z0_mn,pyya_mn)
              call integral(pyyhhc_s,-d,z0_mn,pyyc_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
            call integral(pyyhha_s,z0_mn,z0,pyya)
            call integral(pyyhhc_s,z0_mn,z0,pyyc)
            
            pyya=pyya+pyya_mn
            pyyc=pyyc+pyyc_mn

            stalk_scat=coef*pyyc*conjg(pyya)

          else if (mech_c .eq. gp) then
            
           if(z0_flag .eq. 0) then
               
              call integral(pyyhhc_s,-d,z0_mn,pyyc_mn)
              call integral(pyyhhb_s,-d,z0_mn,pyyb_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
            call integral(pyyhhc_s,z0_mn,z0,pyyc)
            call integral(pyyhhb_s,z0_mn,z0,pyyb)
            
            pyyc=pyyc+pyyc_mn
            pyyb=pyyb+pyyb_mn
            
            stalk_scat=coef*pyyc*conjg(pyyb)
         
         else if (mech_c .eq.pg) then

            if(z0_flag .eq. 0) then
               
              call integral(pyyhhc_s,-d,z0_mn,pyyc_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
            call integral(pyyhhc_s,z0_mn,z0,pyyc)
            
            pyyc=pyyc+pyyc_mn

            stalk_scat=coef*pyyc*conjg(pyyc)

            print*,'pyyc,stalk_scat,coef'
            print*,pyyc,stalk_scat,coef
            
         else if (mech_c .eq.gpg) then

            if(z0_flag .eq. 0) then
               
              call integral(pyyhhc_s,-d,z0_mn,pyyc_mn)
              call integral(pyyhhd_s,-d,z0_mn,pyyd_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
            call integral(pyyhhc_s,z0_mn,z0,pyyc)
            call integral(pyyhhd_s,z0_mn,z0,pyyd)
            
            pyyc=pyyc+pyyc_mn
            pyyd=pyyd+pyyd_mn

            stalk_scat=coef*pyyc*conjg(pyyd)            
        
         else

            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func leaf_scat'
            code=43
            call the_death_code(die,error,code)
            
         endif

       else if (mech .eq. gpg) then

         if (mech_c .eq. p) then

            if(z0_flag .eq. 0) then
               
              call integral(pyyhha_s,-d,z0_mn,pyya_mn)
              call integral(pyyhhd_s,-d,z0_mn,pyyd_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
            call integral(pyyhha_s,z0_mn,z0,pyya)
            call integral(pyyhhd_s,z0_mn,z0,pyyd)
            
            pyya=pyya+pyya_mn
            pyyd=pyyd+pyyd_mn

            stalk_scat=coef*pyyd*conjg(pyya)

          else if (mech_c .eq. gp) then

           if(z0_flag .eq. 0) then
               
              call integral(pyyhhd_s,-d,z0_mn,pyyd_mn)
              call integral(pyyhhb_s,-d,z0_mn,pyyb_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
            call integral(pyyhhd_s,z0_mn,z0,pyyd)
            call integral(pyyhhb_s,z0_mn,z0,pyyb)
            
            pyyd=pyyd+pyyd_mn
            pyyb=pyyb+pyyb_mn
            
            stalk_scat=coef*pyyd*conjg(pyyb)
         
         else if (mech_c .eq. pg) then

            if(z0_flag .eq. 0) then
               
              call integral(pyyhhc_s,-d,z0_mn,pyyc_mn)
              call integral(pyyhhd_s,-d,z0_mn,pyyd_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
            call integral(pyyhhc_s,z0_mn,z0,pyyc)
            call integral(pyyhhd_s,z0_mn,z0,pyyd)
            
            pyyc=pyyc+pyyc_mn
            pyyd=pyyd+pyyd_mn

            stalk_scat=coef*pyyd*conjg(pyyc)
                     
         else if (mech_c .eq. gpg) then

            if(z0_flag .eq. 0) then
               
              call integral(pyyhhd_s,-d,z0_mn,pyyd_mn)
              
            z0_flag=1

          else
              
              continue
              
            endif
            
            call integral(pyyhhd_s,z0_mn,z0,pyyd)
            
            pyyd=pyyd+pyyd_mn
            
            stalk_scat=coef*pyyd*conjg(pyyd)            
        
           else

            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func leaf_scat'
            code=43
            call the_death_code(die,error,code)
            
          endif

        else

          write(error,*)'mech = ',mech,' not recognized',
     &         ' in func leaf_scat'
          code=44
          call the_death_code(die,error,code)
            
        endif
             

       else if (pol .eq. hvhv) then

         if (mech .eq. p) then

           if (mech_c .eq. p) then
             
             if(z0_flag .eq. 0) then
               
               call integral(pxyhva_s,-d,z0_mn,pxya_mn)
               call integral(pyzhva_s,-d,z0_mn,pyza_mn)
             
             z0_flag=1

           else
              
               continue
              
             endif
            
             call integral(pxyhva_s,z0_mn,z0,pxya)
             call integral(pyzhva_s,z0_mn,z0,pyza)
             
             pxya=pxya+pxya_mn
             pyza=pyza+pyza_mn
            
             stalk_scat=coef*
     &            (cos(theta_i)**2*pxya*conjg(pxya)
     &            +sin(theta_i)**2*pyza*conjg(pyza))
            
           else if (mech_c .eq. gp) then
             
             if(z0_flag .eq. 0) then
               
               call integral(pxyhva_s,-d,z0_mn,pxya_mn)
               call integral(pyzhva_s,-d,z0_mn,pyza_mn)
               call integral(pxyhvb_s,-d,z0_mn,pxyb_mn)
               call integral(pyzhvb_s,-d,z0_mn,pyzb_mn)
               
             z0_flag=1

           else
              
               continue
              
             endif
            
             call integral(pxyhva_s,z0_mn,z0,pxya)
             call integral(pyzhva_s,z0_mn,z0,pyza)
             call integral(pxyhvb_s,z0_mn,z0,pxyb)
             call integral(pyzhvb_s,z0_mn,z0,pyzb)
                          
             pxya=pxya+pxya_mn
             pyza=pyza+pyza_mn
             pxyb=pxyb+pxyb_mn
             pyzb=pyzb+pyzb_mn
                        
             stalk_scat=coef*
     &            (-cos(theta_i)**2*pxya*conjg(pxyb)
     &            +sin(theta_i)**2*pyza*conjg(pyzb))

           else if (mech_c .eq. pg) then
             
             if(z0_flag .eq. 0) then
               
               call integral(pxyhva_s,-d,z0_mn,pxya_mn)
               call integral(pyzhva_s,-d,z0_mn,pyza_mn)
               call integral(pxyhvc_s,-d,z0_mn,pxyc_mn)
               call integral(pyzhvc_s,-d,z0_mn,pyzc_mn)
               
             z0_flag=1

           else
              
               continue
              
             endif
            
             call integral(pxyhva_s,z0_mn,z0,pxya)
             call integral(pyzhva_s,z0_mn,z0,pyza)
             call integral(pxyhvc_s,z0_mn,z0,pxyc)
             call integral(pyzhvc_s,z0_mn,z0,pyzc)
                          
             pxya=pxya+pxya_mn
             pyza=pyza+pyza_mn
             pxyc=pxyc+pxyc_mn
             pyzc=pyzc+pyzc_mn
             
            stalk_scat=coef*
     &       (cos(theta_i)**2*pxya*conjg(pxyc)
     &       +sin(theta_i)**2*pyza*conjg(pyzc))

           else if (mech_c .eq. gpg) then
             
             if(z0_flag .eq. 0) then
               
               call integral(pxyhva_s,-d,z0_mn,pxya_mn)
               call integral(pyzhva_s,-d,z0_mn,pyza_mn)
               call integral(pxyhvd_s,-d,z0_mn,pxyd_mn)
               call integral(pyzhvd_s,-d,z0_mn,pyzd_mn)
               
             z0_flag=1

           else
              
               continue
              
             endif
            
             call integral(pxyhva_s,z0_mn,z0,pxya)
             call integral(pyzhva_s,z0_mn,z0,pyza)
             call integral(pxyhvd_s,z0_mn,z0,pxyd)
             call integral(pyzhvd_s,z0_mn,z0,pyzd)
                          
             pxya=pxya+pxya_mn
             pyza=pyza+pyza_mn
             pxyd=pxyd+pxyd_mn
             pyzd=pyzd+pyzd_mn
             
            stalk_scat=coef*
     &           (-cos(theta_i)**2*pxya*conjg(pxyd)
     &           +sin(theta_i)**2*pyza*conjg(pyzd))

          else

            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func leaf_scat'
            code=43
            call the_death_code(die,error,code)
            
          endif
          
        else if (mech .eq. gp) then

          if (mech_c .eq. p) then
            
             if(z0_flag .eq. 0) then
               
               call integral(pxyhva_s,-d,z0_mn,pxya_mn)
               call integral(pyzhva_s,-d,z0_mn,pyza_mn)
               call integral(pxyhvb_s,-d,z0_mn,pxyb_mn)
               call integral(pyzhvb_s,-d,z0_mn,pyzb_mn)
               
             z0_flag=1

           else
              
               continue
              
             endif
            
             call integral(pxyhva_s,z0_mn,z0,pxya)
             call integral(pyzhva_s,z0_mn,z0,pyza)
             call integral(pxyhvb_s,z0_mn,z0,pxyb)
             call integral(pyzhvb_s,z0_mn,z0,pyzb)
                          
             pxya=pxya+pxya_mn
             pyza=pyza+pyza_mn
             pxyb=pxyb+pxyb_mn
             pyzb=pyzb+pyzb_mn
             
            stalk_scat=coef*
     &       (-cos(theta_i)**2*pxyb*conjg(pxya)
     &       +sin(theta_i)**2*pyzb*conjg(pyza))

            
           else if (mech_c .eq. gp) then
             
             if(z0_flag .eq. 0) then
               
               call integral(pxyhvb_s,-d,z0_mn,pxyb_mn)
               call integral(pyzhvb_s,-d,z0_mn,pyzb_mn)
               
             z0_flag=1

           else
              
               continue
              
             endif
            
             call integral(pxyhvb_s,z0_mn,z0,pxyb)
             call integral(pyzhvb_s,z0_mn,z0,pyzb)
                          
             pxyb=pxyb+pxyb_mn
             pyzb=pyzb+pyzb_mn

            stalk_scat=coef*
     &       (cos(theta_i)**2*pxyb*conjg(pxyb)
     &       +sin(theta_i)**2*pyzb*conjg(pyzb))

           else if (mech_c .eq. pg) then

             if(z0_flag .eq. 0) then
               
               call integral(pxyhvc_s,-d,z0_mn,pxyc_mn)
               call integral(pyzhvc_s,-d,z0_mn,pyzc_mn)
               call integral(pxyhvb_s,-d,z0_mn,pxyb_mn)
               call integral(pyzhvb_s,-d,z0_mn,pyzb_mn)
               
             z0_flag=1

           else
              
               continue
              
             endif
            
             call integral(pxyhvc_s,z0_mn,z0,pxyc)
             call integral(pyzhvc_s,z0_mn,z0,pyzc)
             call integral(pxyhvb_s,z0_mn,z0,pxyb)
             call integral(pyzhvb_s,z0_mn,z0,pyzb)
                          
             pxyc=pxyc+pxyc_mn
             pyzc=pyzc+pyzc_mn
             pxyb=pxyb+pxyb_mn
             pyzb=pyzb+pyzb_mn

            stalk_scat=coef*
     &       (-cos(theta_i)**2*pxyb*conjg(pxyc)
     &       +sin(theta_i)**2*pyzb*conjg(pyzc))

           else if (mech_c .eq. gpg) then
             
             if(z0_flag .eq. 0) then
               
               call integral(pxyhvd_s,-d,z0_mn,pxyd_mn)
               call integral(pyzhvd_s,-d,z0_mn,pyzd_mn)
               call integral(pxyhvb_s,-d,z0_mn,pxyb_mn)
               call integral(pyzhvb_s,-d,z0_mn,pyzb_mn)
               
             z0_flag=1

           else
              
               continue
              
             endif
            
             call integral(pxyhvd_s,z0_mn,z0,pxyd)
             call integral(pyzhvd_s,z0_mn,z0,pyzd)
             call integral(pxyhvb_s,z0_mn,z0,pxyb)
             call integral(pyzhvb_s,z0_mn,z0,pyzb)
                          
             pxyd=pxyd+pxyd_mn
             pyzd=pyzd+pyzd_mn
             pxyb=pxyb+pxyb_mn
             pyzb=pyzb+pyzb_mn
             
            stalk_scat=coef*
     &       (+cos(theta_i)**2*pxyb*conjg(pxyd)
     &       +sin(theta_i)**2*pyzb*conjg(pyzd))

          else

             write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func leaf_scat'
            code=43
            call the_death_code(die,error,code)
            
          endif

         else if (mech .eq. pg) then

           if (mech_c .eq. p) then
             
             if(z0_flag .eq. 0) then
               
               call integral(pxyhva_s,-d,z0_mn,pxya_mn)
               call integral(pyzhva_s,-d,z0_mn,pyza_mn)
               call integral(pxyhvc_s,-d,z0_mn,pxyc_mn)
               call integral(pyzhvc_s,-d,z0_mn,pyzc_mn)
               
             z0_flag=1

           else
              
               continue
              
             endif
            
             call integral(pxyhva_s,z0_mn,z0,pxya)
             call integral(pyzhva_s,z0_mn,z0,pyza)
             call integral(pxyhvc_s,z0_mn,z0,pxyc)
             call integral(pyzhvc_s,z0_mn,z0,pyzc)
                          
             pxya=pxya+pxya_mn
             pyza=pyza+pyza_mn
             pxyc=pxyc+pxyc_mn
             pyzc=pyzc+pyzc_mn
             
            stalk_scat=coef*
     &       (cos(theta_i)**2*pxyc*conjg(pxya)
     &       +sin(theta_i)**2*pyzc*conjg(pyza))

            
           else if (mech_c .eq. gp) then
             
             if(z0_flag .eq. 0) then
               
               call integral(pxyhvb_s,-d,z0_mn,pxyb_mn)
               call integral(pyzhvb_s,-d,z0_mn,pyzb_mn)
               call integral(pxyhvc_s,-d,z0_mn,pxyc_mn)
               call integral(pyzhvc_s,-d,z0_mn,pyzc_mn)
               
             z0_flag=1

           else
              
               continue
              
             endif
            
             call integral(pxyhvb_s,z0_mn,z0,pxyb)
             call integral(pyzhvb_s,z0_mn,z0,pyzb)
             call integral(pxyhvc_s,z0_mn,z0,pxyc)
             call integral(pyzhvc_s,z0_mn,z0,pyzc)
                          
             pxyb=pxyb+pxyb_mn
             pyzb=pyzb+pyzb_mn
             pxyc=pxyc+pxyc_mn
             pyzc=pyzc+pyzc_mn
             

            stalk_scat=coef*
     &           (-cos(theta_i)**2*pxyc*conjg(pxyb)
     &           +sin(theta_i)**2*pyzc*conjg(pyzb))

           else if (mech_c .eq. pg) then
             
             if(z0_flag .eq. 0) then
               
               call integral(pxyhvc_s,-d,z0_mn,pxyc_mn)
               call integral(pyzhvc_s,-d,z0_mn,pyzc_mn)
               
             z0_flag=1

           else
              
               continue
              
             endif
            
             call integral(pxyhvc_s,z0_mn,z0,pxyc)
             call integral(pyzhvc_s,z0_mn,z0,pyzc)
                          
             pxyc=pxyc+pxyc_mn
             pyzc=pyzc+pyzc_mn
             
            stalk_scat=coef*
     &       (cos(theta_i)**2*pxyc*conjg(pxyc)
     &       +sin(theta_i)**2*pyzc*conjg(pyzc))

           else if (mech_c .eq. gpg) then
             
             if(z0_flag .eq. 0) then
               
               call integral(pxyhvd_s,-d,z0_mn,pxyd_mn)
               call integral(pyzhvd_s,-d,z0_mn,pyzd_mn)
               call integral(pxyhvc_s,-d,z0_mn,pxyc_mn)
               call integral(pyzhvc_s,-d,z0_mn,pyzc_mn)
               
             z0_flag=1

           else
              
               continue
              
             endif
            
             call integral(pxyhvd_s,z0_mn,z0,pxyd)
             call integral(pyzhvd_s,z0_mn,z0,pyzd)
             call integral(pxyhvc_s,z0_mn,z0,pxyc)
             call integral(pyzhvc_s,z0_mn,z0,pyzc)
                          
             pxyd=pxyd+pxyd_mn
             pyzd=pyzd+pyzd_mn
             pxyc=pxyc+pxyc_mn
             pyzc=pyzc+pyzc_mn
             
            stalk_scat=coef*
     &       (-cos(theta_i)**2*pxyc*conjg(pxyd)
     &       +sin(theta_i)**2*pyzc*conjg(pyzd))

          else

            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func leaf_scat'
            code=43
            call the_death_code(die,error,code)
            
          endif


         else if (mech .eq. gpg) then

           if (mech_c .eq. p) then
             
             if(z0_flag .eq. 0) then
               
               call integral(pxyhvd_s,-d,z0_mn,pxyd_mn)
               call integral(pyzhvd_s,-d,z0_mn,pyzd_mn)
               call integral(pxyhva_s,-d,z0_mn,pxya_mn)
               call integral(pyzhva_s,-d,z0_mn,pyza_mn)
               
             z0_flag=1

           else
              
               continue
              
             endif
            
             call integral(pxyhvd_s,z0_mn,z0,pxyd)
             call integral(pyzhvd_s,z0_mn,z0,pyzd)
             call integral(pxyhva_s,z0_mn,z0,pxya)
             call integral(pyzhva_s,z0_mn,z0,pyza)
                          
             pxyd=pxyd+pxyd_mn
             pyzd=pyzd+pyzd_mn
             pxya=pxya+pxya_mn
             pyza=pyza+pyza_mn
             
            stalk_scat=coef*
     &       (-cos(theta_i)**2*pxyd*conjg(pxya)
     &       +sin(theta_i)**2*pyzd*conjg(pyza))

            
           else if (mech_c .eq. gp) then
             
             if(z0_flag .eq. 0) then
               
               call integral(pxyhvd_s,-d,z0_mn,pxyd_mn)
               call integral(pyzhvd_s,-d,z0_mn,pyzd_mn)
               call integral(pxyhvb_s,-d,z0_mn,pxyb_mn)
               call integral(pyzhvb_s,-d,z0_mn,pyzb_mn)
               
             z0_flag=1

           else
              
               continue
              
             endif
            
             call integral(pxyhvd_s,z0_mn,z0,pxyd)
             call integral(pyzhvd_s,z0_mn,z0,pyzd)
             call integral(pxyhvb_s,z0_mn,z0,pxyb)
             call integral(pyzhvb_s,z0_mn,z0,pyzb)
                          
             pxyd=pxyd+pxyd_mn
             pyzd=pyzd+pyzd_mn
             pxyb=pxyb+pxyb_mn
             pyzb=pyzb+pyzb_mn
             
            stalk_scat=coef*
     &       (cos(theta_i)**2*pxyd*conjg(pxyb)
     &       +sin(theta_i)**2*pyzd*conjg(pyzb))

           else if (mech_c .eq. pg) then

             if(z0_flag .eq. 0) then
               
               call integral(pxyhvd_s,-d,z0_mn,pxyd_mn)
               call integral(pyzhvd_s,-d,z0_mn,pyzd_mn)
               call integral(pxyhvc_s,-d,z0_mn,pxyc_mn)
               call integral(pyzhvc_s,-d,z0_mn,pyzc_mn)
               
             z0_flag=1

           else
              
               continue
              
             endif
            
             call integral(pxyhvd_s,z0_mn,z0,pxyd)
             call integral(pyzhvd_s,z0_mn,z0,pyzd)
             call integral(pxyhvc_s,z0_mn,z0,pxyc)
             call integral(pyzhvc_s,z0_mn,z0,pyzc)
                          
             pxyd=pxyd+pxyd_mn
             pyzd=pyzd+pyzd_mn
             pxyc=pxyc+pxyc_mn
             pyzc=pyzc+pyzc_mn

            stalk_scat=coef*
     &       (-cos(theta_i)**2*pxyd*conjg(pxyc)
     &       +sin(theta_i)**2*pyzd*conjg(pyzc))

           else if (mech_c .eq. gpg) then

             if(z0_flag .eq. 0) then
               
               call integral(pxyhvd_s,-d,z0_mn,pxyd_mn)
               call integral(pyzhvd_s,-d,z0_mn,pyzd_mn)
               
             z0_flag=1

           else
              
               continue
              
             endif
            
             call integral(pxyhvd_s,z0_mn,z0,pxyd)
             call integral(pyzhvd_s,z0_mn,z0,pyzd)
                          
             pxyd=pxyd+pxyd_mn
             pyzd=pyzd+pyzd_mn

            stalk_scat=coef*
     &       (cos(theta_i)**2*pxyd*conjg(pxyd)
     &       +sin(theta_i)**2*pyzd*conjg(pyzd))

          else

            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func leaf_scat'
            code=43
            call the_death_code(die,error,code)
            

          endif

        else

          write(error,*)'mech_c = ',mech_c,' not recognized',
     &         ' in func leaf_scat'
          code=44
          call the_death_code(die,error,code)
            
        endif

      else if (pol .eq. hhvv) then

        if (mech .eq. p) then

          if (mech_c .eq. p) then
            
             if(z0_flag .eq. 0) then
               
               call integral(pyyhha_s,-d,z0_mn,pyya_mn)               
               call integral(pxxvva_s,-d,z0_mn,pxxa_mn)            
               call integral(pzzvva_s,-d,z0_mn,pzza_mn)

             z0_flag=1

           else
              
               continue
              
             endif
            
               call integral(pyyhha_s,z0_mn,z0,pyya)               
               call integral(pxxvva_s,z0_mn,z0,pxxa)            
               call integral(pzzvva_s,z0_mn,z0,pzza)
               
             pxxa=pxxa+pxxa_mn
             pyya=pyya+pyya_mn
             pzza=pzza+pzza_mn             
             
            stalk_scat=coef*
     &       (cos(theta_i)**2*pyya*conjg(pxxa)
     &       +sin(theta_i)**2*pyya*conjg(pzza))

          else if (mech_c .eq. gp) then
            
             if(z0_flag .eq. 0) then
               
               call integral(pyyhha_s,-d,z0_mn,pyya_mn)               
               call integral(pxxvvb_s,-d,z0_mn,pxxb_mn)            
               call integral(pzzvvb_s,-d,z0_mn,pzzb_mn)

             z0_flag=1

           else
              
               continue
              
             endif
            
               call integral(pyyhha_s,z0_mn,z0,pyya)               
               call integral(pxxvvb_s,z0_mn,z0,pxxb)            
               call integral(pzzvvb_s,z0_mn,z0,pzzb)
               
             pxxb=pxxb+pxxb_mn
             pyya=pyya+pyya_mn
             pzzb=pzzb+pzzb_mn
             
            stalk_scat=coef*
     &       (cos(theta_i)**2*pyya*conjg(pxxb)
     &       +sin(theta_i)**2*pyya*conjg(pzzb))

          else if (mech_c .eq. pg) then
            
             if(z0_flag .eq. 0) then
               
               call integral(pyyhha_s,-d,z0_mn,pyya_mn)               
               call integral(pxxvvc_s,-d,z0_mn,pxxc_mn)            
               call integral(pzzvvc_s,-d,z0_mn,pzzc_mn)

             z0_flag=1

           else
              
               continue
              
             endif
            
               call integral(pyyhha_s,z0_mn,z0,pyya)               
               call integral(pxxvvc_s,z0_mn,z0,pxxc)            
               call integral(pzzvvc_s,z0_mn,z0,pzzc)
               
             pxxc=pxxc+pxxc_mn
             pyya=pyya+pyya_mn
             pzzc=pzzc+pzzc_mn
             
            stalk_scat=coef*
     &       (cos(theta_i)**2*pyya*conjg(pxxc)
     &       +sin(theta_i)**2*pyya*conjg(pzzc))

          else if (mech_c .eq. gpg) then

             if(z0_flag .eq. 0) then
               
               call integral(pyyhhd_s,-d,z0_mn,pyyd_mn)               
               call integral(pxxvvd_s,-d,z0_mn,pxxd_mn)            
               call integral(pzzvvd_s,-d,z0_mn,pzzd_mn)

             z0_flag=1

           else
              
               continue
              
             endif
            
               call integral(pyyhhd_s,z0_mn,z0,pyyd)               
               call integral(pxxvvd_s,z0_mn,z0,pxxd)            
               call integral(pzzvvd_s,z0_mn,z0,pzzd)
               
             pxxd=pxxd+pxxd_mn
             pyyd=pyyd+pyyd_mn
             pzzd=pzzd+pzzd_mn

            stalk_scat=coef*
     &           (cos(theta_i)**2*pyya*conjg(pxxd)
     &           +sin(theta_i)**2*pyya*conjg(pzzd))

          else

            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func leaf_scat'
            code=43
            call the_death_code(die,error,code)
            
          endif

        else if (mech .eq. gp) then

          if (mech_c .eq. p) then
            
             if(z0_flag .eq. 0) then
               
               call integral(pyyhhb_s,-d,z0_mn,pyyb_mn)               
               call integral(pxxvva_s,-d,z0_mn,pxxa_mn)            
               call integral(pzzvva_s,-d,z0_mn,pzza_mn)

             z0_flag=1

           else
              
               continue
              
             endif
            
               call integral(pyyhhb_s,z0_mn,z0,pyyb)               
               call integral(pxxvva_s,z0_mn,z0,pxxa)            
               call integral(pzzvva_s,z0_mn,z0,pzza)
               
             pxxa=pxxa+pxxa_mn
             pyyb=pyyb+pyyb_mn
             pzza=pzza+pzza_mn

            stalk_scat=coef*
     &       (-cos(theta_i)**2*pyyb*conjg(pxxa)
     &       +sin(theta_i)**2*pyyb*conjg(pzza))

          else if (mech_c .eq. gp) then
            
             if(z0_flag .eq. 0) then
               
               call integral(pyyhhb_s,-d,z0_mn,pyyb_mn)               
               call integral(pxxvvb_s,-d,z0_mn,pxxb_mn)            
               call integral(pzzvvb_s,-d,z0_mn,pzzb_mn)

             z0_flag=1

           else
              
               continue
              
             endif
            
               call integral(pyyhhb_s,z0_mn,z0,pyyb)               
               call integral(pxxvvb_s,z0_mn,z0,pxxb)            
               call integral(pzzvvb_s,z0_mn,z0,pzzb)
               
             pxxb=pxxb+pxxb_mn
             pyyb=pyyb+pyyb_mn
             pzzb=pzzb+pzzb_mn
             
            stalk_scat=coef*
     &       (-cos(theta_i)**2*pyyb*conjg(pxxb)
     &       +sin(theta_i)**2*pyyb*conjg(pzzb))

          else if (mech_c .eq. pg) then
            
             if(z0_flag .eq. 0) then
               
               call integral(pyyhhb_s,-d,z0_mn,pyyb_mn)               
               call integral(pxxvvc_s,-d,z0_mn,pxxc_mn)            
               call integral(pzzvvc_s,-d,z0_mn,pzzc_mn)

             z0_flag=1

           else
              
               continue
              
             endif
            
               call integral(pyyhhb_s,z0_mn,z0,pyyb)               
               call integral(pxxvvc_s,z0_mn,z0,pxxc)           
               call integral(pzzvvc_s,z0_mn,z0,pzzc)
               
             pxxc=pxxc+pxxc_mn
             pyyb=pyyb+pyyb_mn
             pzzc=pzzc+pzzc_mn
             
            stalk_scat=coef*
     &       (-cos(theta_i)**2*pyyb*conjg(pxxc)
     &       +sin(theta_i)**2*pyyb*conjg(pzzc))

          else if (mech_c .eq. gpg) then
            
             if(z0_flag .eq. 0) then
               
               call integral(pyyhhb_s,-d,z0_mn,pyyb_mn)               
               call integral(pxxvvd_s,-d,z0_mn,pxxd_mn)            
               call integral(pzzvvd_s,-d,z0_mn,pzzd_mn)

             z0_flag=1

           else
              
               continue
              
             endif
            
               call integral(pyyhhb_s,z0_mn,z0,pyyb)               
               call integral(pxxvvd_s,z0_mn,z0,pxxd)           
               call integral(pzzvvd_s,z0_mn,z0,pzzd)
               
             pxxd=pxxd+pxxd_mn
             pyyb=pyyb+pyyb_mn
             pzzd=pzzd+pzzd_mn
             
            stalk_scat=coef*
     &       (-cos(theta_i)**2*pyyb*conjg(pxxd)
     &       +sin(theta_i)**2*pyyb*conjg(pzzd))

          else

            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func leaf_scat'
            code=43
            call the_death_code(die,error,code)
            

          endif


        else if (mech .eq. pg) then

          if (mech_c .eq. p) then
            
             if(z0_flag .eq. 0) then
               
               call integral(pyyhhc_s,-d,z0_mn,pyyc_mn)               
               call integral(pxxvva_s,-d,z0_mn,pxxa_mn)            
               call integral(pzzvva_s,-d,z0_mn,pzza_mn)

             z0_flag=1

           else
              
               continue
              
             endif
            
               call integral(pyyhhc_s,z0_mn,z0,pyyc)               
               call integral(pxxvva_s,z0_mn,z0,pxxa)           
               call integral(pzzvva_s,z0_mn,z0,pzza)
               
             pxxa=pxxa+pxxa_mn
             pyyc=pyyc+pyyc_mn
             pzza=pzza+pzza_mn
             
            stalk_scat=coef*
     &       (-cos(theta_i)**2*pyyc*conjg(pxxa)
     &       +sin(theta_i)**2*pyyc*conjg(pzza))

          else if (mech_c .eq. gp) then
            
             if(z0_flag .eq. 0) then
               
               call integral(pyyhhc_s,-d,z0_mn,pyyc_mn)               
               call integral(pxxvvb_s,-d,z0_mn,pxxb_mn)            
               call integral(pzzvvb_s,-d,z0_mn,pzzb_mn)

             z0_flag=1

           else
              
               continue
              
             endif
            
               call integral(pyyhhc_s,z0_mn,z0,pyyc)               
               call integral(pxxvvb_s,z0_mn,z0,pxxb)           
               call integral(pzzvvb_s,z0_mn,z0,pzzb)
               
             pxxb=pxxb+pxxb_mn
             pyyc=pyyc+pyyc_mn
             pzzb=pzzb+pzzb_mn
             
             stalk_scat=coef*
     &       (-cos(theta_i)**2*pyyc*conjg(pxxb)
     &       +sin(theta_i)**2*pyyc*conjg(pzzb))

          else if (mech_c .eq. pg) then
            
             if(z0_flag .eq. 0) then
               
               call integral(pyyhhc_s,-d,z0_mn,pyyc_mn)               
               call integral(pxxvvc_s,-d,z0_mn,pxxc_mn)            
               call integral(pzzvvc_s,-d,z0_mn,pzzc_mn)

             z0_flag=1

           else
              
               continue
              
             endif
            
               call integral(pyyhhc_s,z0_mn,z0,pyyc)               
               call integral(pxxvvc_s,z0_mn,z0,pxxc)           
               call integral(pzzvvc_s,z0_mn,z0,pzzc)
               
             pxxc=pxxc+pxxc_mn
             pyyc=pyyc+pyyc_mn
             pzzc=pzzc+pzzc_mn
             
            stalk_scat=coef*
     &       (-cos(theta_i)**2*pyyc*conjg(pxxc)
     &       +sin(theta_i)**2*pyyc*conjg(pzzc))

          else if (mech_c .eq. gpg) then
            
             if(z0_flag .eq. 0) then
               
               call integral(pyyhhc_s,-d,z0_mn,pyyc_mn)               
               call integral(pxxvvd_s,-d,z0_mn,pxxd_mn)            
               call integral(pzzvvd_s,-d,z0_mn,pzzd_mn)

             z0_flag=1

           else
              
               continue
              
             endif
            
               call integral(pyyhhc_s,z0_mn,z0,pyyc)               
               call integral(pxxvvd_s,z0_mn,z0,pxxd)           
               call integral(pzzvvd_s,z0_mn,z0,pzzd)
               
             pxxd=pxxd+pxxd_mn
             pyyc=pyyc+pyyc_mn
             pzzd=pzzd+pzzd_mn
             
            stalk_scat=coef*
     &       (-cos(theta_i)**2*pyyc*conjg(pxxd)
     &       +sin(theta_i)**2*pyyc*conjg(pzzd))

          else

            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func leaf_scat'
            code=43
            call the_death_code(die,error,code)
            

          endif

        else if (mech .eq. gpg) then

          if (mech_c .eq. p) then
            
             if(z0_flag .eq. 0) then
               
               call integral(pyyhhd_s,-d,z0_mn,pyyd_mn)               
               call integral(pxxvva_s,-d,z0_mn,pxxa_mn)            
               call integral(pzzvva_s,-d,z0_mn,pzza_mn)

             z0_flag=1

           else
              
               continue
              
             endif
            
               call integral(pyyhhd_s,z0_mn,z0,pyyd)               
               call integral(pxxvva_s,z0_mn,z0,pxxa)           
               call integral(pzzvva_s,z0_mn,z0,pzza)
               
             pxxa=pxxa+pxxa_mn
             pyyd=pyyd+pyyd_mn
             pzza=pzza+pzza_mn
             
            stalk_scat=coef*
     &       (cos(theta_i)**2*pyyd*conjg(pxxa)
     &       +sin(theta_i)**2*pyyd*conjg(pzza))

          else if (mech_c .eq. gp) then
            
             if(z0_flag .eq. 0) then
               
               call integral(pyyhhd_s,-d,z0_mn,pyyd_mn)               
               call integral(pxxvvb_s,-d,z0_mn,pxxb_mn)            
               call integral(pzzvvb_s,-d,z0_mn,pzzb_mn)

             z0_flag=1

           else
              
               continue
              
             endif
            
               call integral(pyyhhd_s,z0_mn,z0,pyyd)               
               call integral(pxxvvb_s,z0_mn,z0,pxxb)           
               call integral(pzzvvb_s,z0_mn,z0,pzzb)
               
             pxxb=pxxb+pxxb_mn
             pyyd=pyyd+pyyd_mn
             pzzb=pzzb+pzzb_mn
             
            stalk_scat=coef*
     &       (cos(theta_i)**2*pyyd*conjg(pxxb)
     &       +sin(theta_i)**2*pyyd*conjg(pzzb))

          else if (mech_c .eq. pg) then

             if(z0_flag .eq. 0) then
               
               call integral(pyyhhd_s,-d,z0_mn,pyyd_mn)               
               call integral(pxxvvc_s,-d,z0_mn,pxxc_mn)            
               call integral(pzzvvc_s,-d,z0_mn,pzzc_mn)

             z0_flag=1

           else
              
               continue
              
             endif
            
               call integral(pyyhhd_s,z0_mn,z0,pyyd)               
               call integral(pxxvvc_s,z0_mn,z0,pxxc)           
               call integral(pzzvvc_s,z0_mn,z0,pzzc)
               
             pxxc=pxxc+pxxc_mn
             pyyd=pyyd+pyyd_mn
             pzzc=pzzc+pzzc_mn

            stalk_scat=coef*
     &       (cos(theta_i)**2*pyyd*conjg(pxxc)
     &       +sin(theta_i)**2*pyyd*conjg(pzzc))

          else if (mech_c .eq. gpg) then
            
             if(z0_flag .eq. 0) then
               
               call integral(pyyhhd_s,-d,z0_mn,pyyd_mn)               
               call integral(pxxvvd_s,-d,z0_mn,pxxd_mn)            
               call integral(pzzvvd_s,-d,z0_mn,pzzd_mn)

             z0_flag=1

           else
              
               continue
              
             endif
            
               call integral(pyyhhd_s,z0_mn,z0,pyyd)               
               call integral(pxxvvd_s,z0_mn,z0,pxxd)           
               call integral(pzzvvd_s,z0_mn,z0,pzzd)
               
             pxxd=pxxd+pxxd_mn
             pyyd=pyyd+pyyd_mn
             pzzd=pzzd+pzzd_mn
             
            stalk_scat=coef*
     &       (cos(theta_i)**2*pyyd*conjg(pxxd)
     &       +sin(theta_i)**2*pyyd*conjg(pzzd))

          else

            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func leaf_scat'
            code=43
            call the_death_code(die,error,code)
            
          endif

        else

          write(error,*)'mech = ',mech,' not recognized',
     &         ' in func leaf_scat'
          code=44
          call the_death_code(die,error,code)

        endif

      else

        write(error,*)'pol = ',pol,' not recognized',
     &       ' in func leaf_scat'
        code=45
        call the_death_code(die,error,code)
            
      endif

      return

      end

c************************************************************************
      
      complex function pxxvva_s(z)

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

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

       real rho,x,z,ell
       real dl,cosb,sinb,drho,sina,cosa
       real a,alpha_s,ell_0
       real theta_i,phi_i,k_0
       real d,del_leaf,z_lc,hole
       
       common /stalk_struct/ a,alpha_s,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ d,del_leaf,z_lc,hole
      
       ii=cmplx(0.0,1.0)

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

       rho=a*(z+d)
       x=rho*cosa

       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

       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

       pxxvva_s=phia*(cosa**2*cosb**2*p11_s
     &  +sina**2*p22_s
     &  +cosa**2*sinb**2*p33_s)*dl

       return

       end
      

c************************************************************************
      
      complex function pxxvvb_s(z)

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


       complex ii,phib
       complex p11_s,p22_s,p33_s

       real rho,x,z,ell
       real dl,cosb,sinb,drho,sina,cosa
       real a,alpha_s,ell_0
       real theta_i,phi_i,k_0
       real d,del_leaf,z_lc,hole
       
       common /stalk_struct/ a,alpha_s,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ d,del_leaf,z_lc,hole
      
       ii=cmplx(0.0,1.0)      

       cosa=cos(alpha_s)
       sina=sin(alpha_s)
       
       rho=a*(z+d)
       x=rho*cosa

       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
       
       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

       pxxvvb_s=phib*(cosa**2*cosb**2*p11_s
     &  +sina**2*p22_s
     &  +cosa**2*sinb**2*p33_s)*dl

       return

       end       

c************************************************************************
      
      complex function pxxvvc_s(z)

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

       complex ii,phic
       complex p11_s,p22_s,p33_s

       real rho,x,z,ell
       real dl,cosb,sinb,drho,sina,cosa
       real a,alpha_s,ell_0
       real theta_i,phi_i,k_0
       real d,del_leaf,z_lc,hole
       
       common /stalk_struct/ a,alpha_s,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ d,del_leaf,z_lc,hole
      
       ii=cmplx(0.0,1.0)      

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

       rho=a*(z+d)
       x=rho*cosa

       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

       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

       pxxvvc_s=phic*(cosa**2*cosb**2*p11_s
     &  +sina**2*p22_s
     &  +cosa**2*sinb**2*p33_s)*dl

       return

       end              

c************************************************************************
      
      complex function pxxvvd_s(z)

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

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

       real rho,x,z,ell
       real dl,cosb,sinb,drho,sina,cosa
       real a,alpha_s,ell_0
       real theta_i,phi_i,k_0
       real d,del_leaf,z_lc,hole
       
       common /stalk_struct/ a,alpha_s,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ d,del_leaf,z_lc,hole
      
       ii=cmplx(0.0,1.0)      

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

       rho=a*(z+d)
       x=rho*cosa

       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

       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

       pxxvvd_s=phid*(cosa**2*cosb**2*p11_s
     &  +sina**2*p22_s
     &  +cosa**2*sinb**2*p33_s)*dl

       return
       end              

c************************************************************************
      
      complex function pzzvva_s(z)

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

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

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

       cosa=cos(alpha_s)

       rho=a*(z+d)
       x=rho*cosa

       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
       
       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

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

       return

       end
      

c************************************************************************
      
      complex function pzzvvb_s(z)

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


       complex ii,phib
       complex p11_s,p22_s,p33_s

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

       cosa=cos(alpha_s)

       rho=a*(z+d)
       x=rho*cosa

       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

       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

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

       return

       end       

c************************************************************************
      
      complex function pzzvvc_s(z)

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

       complex ii,phic
       complex p11_s,p22_s,p33_s

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

       cosa=cos(alpha_s)
       
       rho=a*(z+d)
       x=rho*cosa

       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       

       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

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

       return

       end              

c************************************************************************
      
      complex function pzzvvd_s(z)

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

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

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

       cosa=cos(alpha_s)

       rho=a*(z+d)
       x=rho*cosa

       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       

       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

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

       return

       end              

c************************************************************************
      
      complex function pxzvva_s(z)

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

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

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

       cosa=cos(alpha_s)

       rho=a*(z+d)
       x=rho*cosa

       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
       
       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

       pxzvva_s=phia*(cosa*cosb*sinb*(p33_s-p11_s))

       return

       end
      

c************************************************************************
      
      complex function pxzvvb_s(z)

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


       complex ii,phib
       complex p11_s,p22_s,p33_s

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

       cosa=cos(alpha_s)

       rho=a*(z+d)
       x=rho*cosa

       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

       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

       pxzvvb_s=phib*(cosa*cosb*sinb*(p33_s-p11_s))

       return
       end       

c************************************************************************
      
      complex function pxzvvc_s(z)

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

       complex ii,phic
       complex p11_s,p22_s,p33_s

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

       cosa=cos(alpha_s)

       rho=a*(z+d)
       x=rho*cosa

       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       

       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

       pxzvvc_s=phic*(cosa*cosb*sinb*(p33_s-p11_s))

       return

       end              

c************************************************************************
      
      complex function pxzvvd_s(z)

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

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

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

       cosa=cos(alpha_s)

       rho=a*(z+d)
       x=rho*cosa

       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       

       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

       pxzvvd_s=phid*(cosa*cosb*sinb*(p33_s-p11_s))

       return

       end              

c************************************************************************
      
      complex function pyyhha_s(z)

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

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

       real rho,x,z,ell
       real dl,cosb,sinb,drho,sina,cosa
       real a,alpha_s,ell_0
       real theta_i,phi_i,k_0
       real d,del_leaf,z_lc,hole
       
       common /stalk_struct/ a,alpha_s,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ d,del_leaf,z_lc,hole
      
       ii=cmplx(0.0,1.0)      

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

       rho=a*(z+d)
       x=rho*cosa

       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

       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

       pyyhha_s=phia*(sina**2*cosb**2*p11_s
     &  +cosa**2*p22_s
     &  +sina**2*sinb**2*p33_s)*dl

       return

       end
      

c************************************************************************
      
      complex function pyyhhb_s(z)

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

       complex ii,phib
       complex p11_s,p22_s,p33_s

       real rho,x,z,ell
       real dl,cosb,sinb,drho,sina,cosa
       real a,alpha_s,ell_0
       real theta_i,phi_i,k_0
       real d,del_leaf,z_lc,hole
       
       common /stalk_struct/ a,alpha_s,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ d,del_leaf,z_lc,hole
      
       ii=cmplx(0.0,1.0)      

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

       rho=a*(z+d)
       x=rho*cosa

       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

       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

       pyyhhb_s=phib*(sina**2*cosb**2*p11_s
     &  +cosa**2*p22_s
     &  +sina**2*sinb**2*p33_s)*dl

C       print*,'a,alpha_s,p11_s,p22_s,p33_s,phib,pyyhhb_s'
C       print*, a,alpha_s,p11_s,p22_s,p33_s,phib,pyyhhb_s
       
       return

       end       

c************************************************************************
      
      complex function pyyhhc_s(z)

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

       complex ii,phic
       complex p11_s,p22_s,p33_s

       real rho,x,z,ell
       real dl,cosb,sinb,drho,sina,cosa
       real a,alpha_s,ell_0
       real theta_i,phi_i,k_0
       real d,del_leaf,z_lc,hole
       
       common /stalk_struct/ a,alpha_s,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ d,del_leaf,z_lc,hole
      
       ii=cmplx(0.0,1.0)      

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

       rho=a*(z+d)
       x=rho*cosa

       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

       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

       pyyhhc_s=phic*(sina**2*cosb**2*p11_s
     &  +cosa**2*p22_s
     &  +sina**2*sinb**2*p33_s)*dl

       return

       end              

c************************************************************************
      
      complex function pyyhhd_s(z)

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

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

       real rho,x,z,ell
       real dl,cosb,sinb,drho,sina,cosa
       real a,alpha_s,ell_0
       real theta_i,phi_i,k_0
       real d,del_leaf,z_lc,hole
       
       common /stalk_struct/ a,alpha_s,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ d,del_leaf,z_lc,hole
      
       ii=cmplx(0.0,1.0)      

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

       rho=a*(z+d)
       x=rho*cosa

       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
       
       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

       pyyhhd_s=phid*(sina**2*cosb**2*p11_s
     &  +cosa**2*p22_s
     &  +sina**2*sinb**2*p33_s)*dl
  
       return

       end              

c************************************************************************
      
      complex function pxyhva_s(z)

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

       real rho,x,z,ell
       real dl,cosb,sinb,drho,sina,cosa
       real a,alpha_s,ell_0
       real theta_i,phi_i,k_0
       real d,del_leaf,z_lc,hole
       
       common /stalk_struct/ a,alpha_s,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ d,del_leaf,z_lc,hole
      
       ii=cmplx(0.0,1.0)      

       rho=a*(z+d)

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

       x=rho*cosa

       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

       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

       pxyhva_s=phia*(sina*cosa*cosb**2*p11_s
     &  -cosa*sina*p22_s
     &  +sina*cosa*sinb**2*p33_s)*dl

       return

       end
      

c************************************************************************
      
      complex function pxyhvb_s(z)

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

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

       real rho,x,z,ell
       real dl,cosb,sinb,drho,sina,cosa
       real a,alpha_s,ell_0
       real theta_i,phi_i,k_0
       real d,del_leaf,z_lc,hole
       
       common /stalk_struct/ a,alpha_s,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ d,del_leaf,z_lc,hole
      
       ii=cmplx(0.0,1.0)      

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

       rho=a*(z+d)
       x=rho*cosa

       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

       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

       pxyhvb_s=phib*(sina*cosa*cosb**2*p11_s
     &  -cosa*sina*p22_s
     &  +sina*cosa*sinb**2*p33_s)*dl

       return

       end       

c************************************************************************
      
      complex function pxyhvc_s(z)

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

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

       real rho,x,z,ell
       real dl,cosb,sinb,drho,sina,cosa
       real a,alpha_s,ell_0
       real theta_i,phi_i,k_0
       real d,del_leaf,z_lc,hole
       
       common /stalk_struct/ a,alpha_s,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ d,del_leaf,z_lc,hole
      
       ii=cmplx(0.0,1.0)      

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

       rho=a*(z+d)
       x=rho*cosa

       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

       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

       pxyhvc_s=phic*(sina*cosa*cosb**2*p11_s
     &  -cosa*sina*p22_s
     &  +sina*cosa*sinb**2*p33_s)*dl
       
       return

       end              

c************************************************************************
      
      complex function pxyhvd_s(z)

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


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

       real rho,x,z,ell
       real dl,cosb,sinb,drho,sina,cosa
       real a,alpha_s,ell_0
       real theta_i,phi_i,k_0
       real d,del_leaf,z_lc,hole
       
       common /stalk_struct/ a,alpha_s,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ d,del_leaf,z_lc,hole
      
       ii=cmplx(0.0,1.0)      

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

       rho=a*(z+d)
       x=rho*cosa

       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
       

       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

       pxyhvd_s=phid*(sina*cosa*cosb**2*p11_s
     &  -cosa*sina*p22_s
     & +sina*cosa*sinb**2*p33_s)*dl

       return

       end

      
c************************************************************************
      
      complex function pyzhva_s(z)

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

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

       real rho,x,z,ell
       real dl,cosb,sinb,drho,sina,cosa
       real a,alpha_s,ell_0
       real theta_i,phi_i,k_0
       real d,del_leaf,z_lc,hole
       
       common /stalk_struct/ a,alpha_s,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ d,del_leaf,z_lc,hole
      
       ii=cmplx(0.0,1.0)      

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

       rho=a*(z+d)
       x=rho*cosa

       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

       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

       pyzhva_s=phia*(sina*sinb*cosb*(p33_s-p11_s))*dl

       return

       end
      

c************************************************************************
      
      complex function pyzhvb_s(z)

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

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

       real rho,x,z,ell
       real dl,cosb,sinb,drho,sina,cosa
       real a,alpha_s,ell_0
       real theta_i,phi_i,k_0
       real d,del_leaf,z_lc,hole
       
       common /stalk_struct/ a,alpha_s,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ d,del_leaf,z_lc,hole
      
       ii=cmplx(0.0,1.0)      

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

       rho=a*(z+d)
       x=rho*cosa

       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

       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

       pyzhvb_s=phib*(sina*sinb*cosb*(p33_s-p11_s))*dl

       return

       end       

c************************************************************************
      
      complex function pyzhvc_s(z)

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

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

       real rho,x,z,ell
       real dl,cosb,sinb,drho,sina,cosa
       real a,alpha_s,ell_0
       real theta_i,phi_i,k_0
       real d,del_leaf,z_lc,hole
       
       common /stalk_struct/ a,alpha_s,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ d,del_leaf,z_lc,hole
      
       ii=cmplx(0.0,1.0)      

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

       rho=a*(z+d)
       x=rho*cosa

       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

       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

       pyzhvc_s=phic*(sina*sinb*cosb*(p33_s-p11_s))*dl

       return

       end              

c************************************************************************
      
      complex function pyzhvd_s(z)

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

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

       real rho,x,z,ell
       real dl,cosb,sinb,drho,sina,cosa
       real a,alpha_s,ell_0
       real theta_i,phi_i,k_0
       real d,del_leaf,z_lc,hole
       
       common /stalk_struct/ a,alpha_s,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /layer_struct/ d,del_leaf,z_lc,hole
      
       ii=cmplx(0.0,1.0)      

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

       rho=a*(z+d)
       x=rho*cosa

       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

       ell=dl*(z+d)

       call stalk_tensor(ell,ell_0, p11_s,p22_s,p33_s)

       pyzhvd_s=phid*(sina*sinb*cosb*(p33_s-p11_s))*dl

       return

       end              

