


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

      
      subroutine stalk_leaf(ss_sl,ss_ls)


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


      do 10 i=1,4

        do 20 j=1,4

          do 30 k=1,5
            
           if(i .eq. 1) then
            
            if (j .eq. 1) then

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

                call leaf_field_int(i,j,k, t)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))

              else if(k .eq. 2) then

                call stalk_lphi_int(i,j,k, s)

                call leaf_field_int(i,j,k, t)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))

              else if(k .eq. 3) then

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

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

                call stalk_lphi_int(i,j,k, s)

                call leaf_field_int(i,j,k, t)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))
                
              else

                print*, 'error in stalk_leaf a'

              endif

            else if(j .eq. 2) then

              if (k .eq. 1) then

                call stalk_lphi_int(i,j,k, s)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))
                
              else if(k .eq. 2) then
                
                call stalk_lphi_int(i,j,k, s)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))
                
              else if(k .eq. 3) then

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

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

              else if(k .eq. 4) then

                call stalk_lphi_int(i,j,k, s)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))

              else

                print*, 'error in stalk_leaf b'

              endif

            else if(j .eq. 3) then

              if(k .eq. 1) then

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
                ss_ls(i,k,j)=ss_ls(i,k,2)

              else if(k .eq. 2) then

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
                ss_ls(i,k,j)=ss_ls(i,k,2)

                
              else if(k .eq. 3) then

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
                ss_ls(i,k,j)=ss_ls(i,k,2)

                
              else if(k .eq. 4) then

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
                ss_ls(i,k,j)=ss_ls(i,k,2)
                
              else

                print*, 'error in  stalk_leaf c'

              endif

            else if (j .eq. 4) then

              if (k .eq. 1) then

                call stalk_lphi_int(i,j,k, s)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))


              else if (k .eq. 2) then

                call stalk_lphi_int(i,j,k, s)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))

              else if (k .eq. 3) then

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

                ss_ls(i,k,j)=ss_ls(i,2,j)
               
              else if (k .eq. 4) then
                
                call stalk_lphi_int(i,j,k, s)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))


              else

                print*, 'error in  stalk_leaf d'
                
              endif

            else
              
              print*, 'error in  stalk_leaf e'
                
            endif
            
            
           else if(i .eq. 2) then
            
            if (j .eq. 1) then

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

                call leaf_field_int(i,j,k, t)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))

              else if(k .eq. 2) then

                call stalk_lphi_int(i,j,k, s)

                call leaf_field_int(i,j,k, t)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))

              else if(k .eq. 3) then

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

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

                call stalk_lphi_int(i,j,k, s)

                call leaf_field_int(i,j,k, t)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))
                
              else

                print*, 'error in stalk_leaf f'

              endif

            else if(j .eq. 2) then

              if (k .eq. 1) then

                call stalk_lphi_int(i,j,k, s)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))
                
              else if(k .eq. 2) then
                
                call stalk_lphi_int(i,j,k, s)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))
                
              else if(k .eq. 3) then

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

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

              else if(k .eq. 4) then

                call stalk_lphi_int(i,j,k, s)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))

              else

                print*, 'error in stalk_leaf g'

              endif

            else if(j .eq. 3) then

              if(k .eq. 1) then

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
                ss_ls(i,k,j)=ss_ls(i,k,2)

              else if(k .eq. 2) then

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
                ss_ls(i,k,j)=ss_ls(i,k,2)

                
              else if(k .eq. 3) then

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
                ss_ls(i,k,j)=ss_ls(i,k,2)

                
              else if(k .eq. 4) then

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
                ss_ls(i,k,j)=ss_ls(i,k,2)
                
              else

                print*, 'error in  stalk_leaf h'

              endif

            else if (j .eq. 4) then

              if (k .eq. 1) then

                call stalk_lphi_int(i,j,k, s)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))


              else if (k .eq. 2) then

                call stalk_lphi_int(i,j,k, s)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))

              else if (k .eq. 3) then

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

                ss_ls(i,k,j)=ss_ls(i,2,j)
               
              else if (k .eq. 4) then
                
                call stalk_lphi_int(i,j,k, s)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))


              else

                print*, 'error in  stalk_leaf i'
                
              endif

            else
              
              print*, 'error in  stalk_leaf j'
                
            endif

          else if(i .eq. 3)
            
            ss_sl(i,j,k)=cmplx(0.0,0.0)

            ss_ls(i,k,j)=cmplx(0.0,0.0)
                
            endif
            

          else if(i .eq. 4) then
            
            if (j .eq. 1) then

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

                ss_sl(i,j,k)=s(i,j,k)*t(1,k)

              else if(k .eq. 2) then

                call stalk_lphi_int(i,j,k, s)

                ss_sl(i,j,k)=s(i,j,k)*t(1,k)

              else if(k .eq. 3) then

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

              else if(k .eq. 4) then

                call stalk_lphi_int(i,j,k, s)

                ss_sl(i,j,k)=s(i,j,k)*t(1,k)
                
              else

                print*, 'error in stalk_leaf p'

              endif

            else if(j .eq. 2) then

              if (k .eq. 1) then

                call stalk_lphi_int(i,j,k, s)

                ss_sl(i,j,k)=s(i,j,k)*t(1,k)

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

                ss_sl(i,j,k)=s(i,j,k)*t(1,k)

              else if(k .eq. 3) then

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

              else if(k .eq. 4) then

                call stalk_lphi_int(i,j,k, s)

                ss_sl(i,j,k)=s(i,j,k)*t(1,k)

              else

                print*, 'error in stalk_leaf q'

              endif

            else if(j .eq. 3) then

              if(k .eq. 1) then

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

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
              else if(k .eq. 3) then

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

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
              else

                print*, 'error in  stalk_leaf r'

              endif

            else if (j .eq. 4) then

              if (k .eq. 1) then

                call stalk_lphi_int(i,j,k, s)

                ss_sl(i,j,k)=s(i,j,k)*t(1,k)

              else if (k .eq. 2) then

                call stalk_lphi_int(i,j,k, s)

                ss_sl(i,j,k)=s(i,j,k)*t(1,k)

              else if (k .eq. 3) then

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

              else if (k .eq. 4) then
                
                call stalk_lphi_int(i,j,k, s)

                ss_sl(i,j,k)=s(i,j,k)*t(1,k)

              else

                print*, 'error in  stalk_leaf s'
                
              endif

            else
              
              print*, 'error in  stalk_leaf t'
                
            endif

          else if(i .eq. 5) then
            
            if (j .eq. 1) then

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

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(2,k))

              else if(k .eq. 2) then

                call stalk_lphi_int(i,j,k, s)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(2,k))

              else if(k .eq. 3) then

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

                call stalk_lphi_int(i,j,k, s)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(2,k))
                
              else

                print*, 'error in stalk_leaf u'

              endif

            else if(j .eq. 2) then

              if (k .eq. 1) then

                call stalk_lphi_int(i,j,k, s)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(2,k))
                
              else if(k .eq. 2) then
                
                call stalk_lphi_int(i,j,k, s)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(2,k))
                
              else if(k .eq. 3) then

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

              else if(k .eq. 4) then

                call stalk_lphi_int(i,j,k, s)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(2,k))

              else

                print*, 'error in stalk_leaf v'

              endif
              
            else if(j .eq. 3) then

              if(k .eq. 1) then

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

              else if(k .eq. 2) then

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

              else if(k .eq. 3) then

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

              else if(k .eq. 4) then

                ss_ls(i,k,j)=ss_ls(i,k,2)
                
              else

                print*, 'error in  stalk_leaf w'

              endif

            else if (j .eq. 4) then

              if (k .eq. 1) then

                call stalk_lphi_int(i,j,k, s)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(2,k))

              else if (k .eq. 2) then

                call stalk_lphi_int(i,j,k, s)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(2,k))

              else if (k .eq. 3) then

                ss_ls(i,k,j)=ss_ls(i,2,j)
               
              else if (k .eq. 4) then
                
                call stalk_lphi_int(i,j,k, s)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(2,k))

              else

                print*, 'error in  stalk_leaf x'
                
              endif

            else
              
              print*, 'error in  stalk_leaf y'
                
            endif

          else
            
              print*, 'error in  stalk_leaf y'

          endif

          
   30     continue  

   20   continue

   10 continue

      return

      end

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

      subroutine stalk_lphi_int(m,n,p, s)

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

      common /mn/ m,n,p


      include 'constants.include'

      external h_sl
      call qgausx(h_l,pi,mpi,r)

      return

      end
      

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

      function h_sl(phii)

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

      external g_sl
      common /stalk_stat/ amin,amax
      common /phiv/ phi,a,z0

      phi=phii

      call qgausy(g_sl,amin,amax,s)

      h_sl=s(m,n,p)

      return

      end

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

      function g_sl(aa)

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


      external i_sl
      common /stalk_stat/ amin,amax
      common /phiv/ phi,a,z0

      a=aa

      call qgausy(i_sl,z0min,z0max,s)

      g_sl=s(

      return

      end

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

      function i_sl(zz0)

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

      external stalk_phil_ave

      common /phiv/ phi,a,z0
      common /mn/  m,n,p

      z0=zz0
      
      i_sl=stalk_phil_ave(phi,a,z0,m,n,p)

      return

      end

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

      function stalk_phil_ave(alpha,s,i,j,k)

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

      external p11vva,p33vva,p13vva
      
      common /incident/ theta_i,phi_i
      

      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)

            call phivva_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
             *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (cos(theta_i)**2*p11a
     &       -2.0*cos(theta_i)*sin(theta_i)*p13a
     &       +sin(theta_i)**2*p33a)


          else if (k .eq. 2) then

            call qgaus(p11vva_s,dm,z0,p11a)

            call qgaus(p33vva_s,dm,z0,p33a)

            call qgaus(p13vva_s,dm,z0,p13a)

            call phivvb_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
             *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (cos(theta_i)**2*p11a
     &       -2.0*cos(theta_i)*sin(theta_i)*p13a
     &       +sin(theta_i)**2*p33a)
             
          else if (k .eq. 3) then

            call qgaus(p11vva_s,dm,z0,p11a)

            call qgaus(p33vva_s,dm,z0,p33a)

            call qgaus(p13vva_s,dm,z0,p13a)

            call phivvc_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
             *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (cos(theta_i)**2*p11a
     &       -2.0*cos(theta_i)*sin(theta_i)*p13a
     &       +sin(theta_i)**2*p33a)

          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 phivvd_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
             *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (cos(theta_i)**2*p11a
     &       -2.0*cos(theta_i)*sin(theta_i)*p13a
     &       +sin(theta_i)**2*p33a)
                      
          else

            print*, 'error in stalk_phil_ave a'

          endif

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


            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)

            call phivva_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (-cos(theta_i)**2*p11b*sin(theta_i)**2*p33b)
            
          else if (k .eq. 2) then


            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)

            call phivvb_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (-cos(theta_i)**2*p11b*sin(theta_i)**2*p33b)

          else if (k .eq. 3)


            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)

            call phivvc_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (-cos(theta_i)**2*p11b*sin(theta_i)**2*p33b)

          else if (k .eq. 4)

            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)

            call phivvd_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (-cos(theta_i)**2*p11b*sin(theta_i)**2*p33b)

          else

            print*, 'error in stalk_phil_ave b'

          endif

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


            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)

            call phivva_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (-cos(theta_i)**2*p11c*sin(theta_i)**2*p33c)
            
          else if (k .eq. 2) then


            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)

            call phivvc_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (-cos(theta_i)**2*p11c*sin(theta_i)**2*p33c)

          else if (k .eq. 3)


            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)

            call phivvc_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (-cos(theta_i)**2*p11c*sin(theta_i)**2*p33c)

          else if (k .eq. 4)

            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)

            call phivvd_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (-cos(theta_i)**2*p11c*sin(theta_i)**2*p33c)

          else

            print*, 'error in stalk_phil_ave c'

          endif

        else if (j .eq. 4) then

          if (k .eq. 1) then

            call qgaus(p11vvd_s,dm,z0,p11d)

            call qgaus(p33vvd_s,dm,z0,p33d)

            call qgaus(p13vvd_s,dm,z0,p13d)

            call phivvd_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
             *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (cos(theta_i)**2*p11d
     &       +2.0*cos(theta_i)*sin(theta_i)*p13d
     &       +sin(theta_i)**2*p33d)


          else if (k .eq. 2) then

            call qgaus(p11vvd_s,dm,z0,p11d)

            call qgaus(p33vvd_s,dm,z0,p33d)

            call qgaus(p13vvd_s,dm,z0,p13d)

            call phivvb_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
             *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (cos(theta_i)**2*p11d
     &       +2.0*cos(theta_i)*sin(theta_i)*p13d
     &       +sin(theta_i)**2*p33d)
             
          else if (k .eq. 3) then

            call qgaus(p11vvd_s,dm,z0,p11d)

            call qgaus(p33vvd_s,dm,z0,p33d)

            call qgaus(p13vvd_s,dm,z0,p13d)

            call phivvc_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
             *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (cos(theta_i)**2*p11d
     &       +2.0*cos(theta_i)*sin(theta_i)*p13d
     &       +sin(theta_i)**2*p33d)

          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)

            call phivvd_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
             *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (cos(theta_i)**2*p11d
     &       +2.0*cos(theta_i)*sin(theta_i)*p13d
     &       +sin(theta_i)**2*p33d)
                      
          else

            print*, 'error in stalk_phil_ave d'

          endif

        end if

      else if (i .eq. 2) then

        if (j .eq. 1) then

          if (k .eq. 1) then

            call qgaus(p22hha_s,dm,z0,p22a)

            call phihha_l(phi_l)            

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

          else if (k .eq. 2) then

            call qgaus(p22hha_s,dm,z0,p22a)

            call phihhb_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22a*conjg(phi_l)
         
         else if (k .eq. 3) then

            call qgaus(p22hha_s,dm,z0,p22a)

            call phihhc_l(phi_l)  

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

            call qgaus(p22hha_s,dm,z0,p22a)

            call phihhd_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22a*conjg(phi_l)            
        
         else

           print*, 'error in stalk_phil_ave e'

         endif

       else if (j .eq. 2)

          if (k .eq. 1) then

            call qgaus(p22hhb_s,dm,z0,p22b)

            call phihha_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22b*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p22hhb_s,dm,z0,p22b)

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

            call qgaus(p22hhb_s,dm,z0,p22b)

            call phihhc_l(phi_l)  

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

            call qgaus(p22hhb_s,dm,z0,p22b)

            call phihhd_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22b*conjg(phi_l)            
        
         else

           print*, 'error in stalk_phil_ave f'

         endif

       else if (j .eq. 3)

          if (k .eq. 1) then

            call qgaus(p22hhc_s,dm,z0,p22c)

            call phihha_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22c*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p22hhc_s,dm,z0,p22c)

            call phihhb_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22c*conjg(phi_l)
         
         else if (k .eq. 3) then

            call qgaus(p22hhc_s,dm,z0,p22c)

            call phihhc_l(phi_l)  

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

            call qgaus(p22hhc_s,dm,z0,p22c)

            call phihhd_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22c*conjg(phi_l)            
        
         else

           print*, 'error in stalk_phil_ave g'

         endif

       else if (j .eq. 4)

          if (k .eq. 1) then

            call qgaus(p22hhd_s,dm,z0,p22d)

            call phihha_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22d*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p22hhd_s,dm,z0,p22d)

            call phihhb_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22d*conjg(phi_l)
         
         else if (k .eq. 3) then

            call qgaus(p22hhd_s,dm,z0,p22d)

            call phihhc_l(phi_l)  

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

            call qgaus(p22hhd_s,dm,z0,p22d)

            call phihhd_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22d*conjg(phi_l)            
        
         else

           print*, 'error in stalk_phil_ave h'

         endif

       else if (i .eq. 3) then

        if (j .eq. 1) then

          if (k .eq. 1) then

            call qgaus(p21hva_s,dm,z0,p21a)

            call qgaus(p23hva_s,dm,z0,p23a)

            call phihva_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)*p21a-sin(theta_i)*p23a)*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p21hva_s,dm,z0,p21a)

            call qgaus(p23hva_s,dm,z0,p23a)

            call phihvb_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)*p21a-sin(theta_i)*p23a)*conjg(phi_l)

         
         else if (k .eq. 3) then

            call qgaus(p21hva_s,dm,z0,p21a)

            call qgaus(p23hva_s,dm,z0,p23a)

            call phihvc_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)*p21a-sin(theta_i)*p23a)*conjg(phi_l)
                     
         else if (k .eq. 4) then

            call qgaus(p21hva_s,dm,z0,p21a)

            call qgaus(p23hva_s,dm,z0,p23a)

            call phihvd_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)*p21a-sin(theta_i)*p23a)*conjg(phi_l)

         else

           print*, 'error in stalk_phil_ave e'

         endif

       else if (j .eq. 2)

          if (k .eq. 1) then

            call qgaus(p21hvb_s,dm,z0,p21b)

            call qgaus(p23hvb_s,dm,z0,p23b)

            call phihva_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)*p21b-sin(theta_i)*p23b)*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p21hvb_s,dm,z0,p21b)

            call qgaus(p23hvb_s,dm,z0,p23b)

            call phihvb_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)*p21b-sin(theta_i)*p23b)*conjg(phi_l)

         
         else if (k .eq. 3) then

            call qgaus(p21hvb_s,dm,z0,p21b)

            call qgaus(p23hvb_s,dm,z0,p23b)

            call phihvc_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)*p21b-sin(theta_i)*p23b)*conjg(phi_l)
                     
         else if (k .eq. 4) then

            call qgaus(p21hvb_s,dm,z0,p21b)

            call qgaus(p23hvb_s,dm,z0,p23b)

            call phihvd_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)*p21b-sin(theta_i)*p23b)*conjg(phi_l)

         else

           print*, 'error in stalk_phil_ave f'

         endif

       else if (j .eq. 3)

          if (k .eq. 1) then

            call qgaus(p21hvc_s,dm,z0,p21c)

            call qgaus(p23hvc_s,dm,z0,p23c)

            call phihva_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)*p21c-sin(theta_i)*p23c)*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p21hvc_s,dm,z0,p21c)

            call qgaus(p23hvc_s,dm,z0,p23c)

            call phihvc_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)*p21c-sin(theta_i)*p23c)*conjg(phi_l)

         
         else if (k .eq. 3) then

            call qgaus(p21hvc_s,dm,z0,p21c)

            call qgaus(p23hvc_s,dm,z0,p23c)

            call phihvc_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)*p21c-sin(theta_i)*p23c)*conjg(phi_l)
                     
         else if (k .eq. 4) then

            call qgaus(p21hvc_s,dm,z0,p21c)

            call qgaus(p23hvc_s,dm,z0,p23c)

            call phihvd_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)*p21c-sin(theta_i)*p23c)*conjg(phi_l)

         else
           

           print*, 'error in stalk_phil_ave g'

         endif

       else if (j .eq. 4)

          if (k .eq. 1) then

            call qgaus(p21hvd_s,dm,z0,p21d)

            call qgaus(p23hvd_s,dm,z0,p23d)

            call phihva_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)*p21d-sin(theta_i)*p23d)*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p21hvd_s,dm,z0,p21d)

            call qgaus(p23hvd_s,dm,z0,p23d)

            call phihvd_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)*p21d-sin(theta_i)*p23d)*conjg(phi_l)

         
         else if (k .eq. 3) then

            call qgaus(p21hvd_s,dm,z0,p21d)

            call qgaus(p23hvd_s,dm,z0,p23d)

            call phihvc_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)*p21d-sin(theta_i)*p23d)*conjg(phi_l)
                     
         else if (k .eq. 4) then

            call qgaus(p21hvd_s,dm,z0,p21d)

            call qgaus(p23hvd_s,dm,z0,p23d)

            call phihvd_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)*p21d-sin(theta_i)*p23d)*conjg(phi_l)

         else

           print*, 'error in stalk_phil_ave h'

         endif

      else if (i .eq. 4) then

        if (j .eq. 1) then

          if (k .eq. 1) then

            call qgaus(p22hha_s,dm,z0,p22a)

            call phivva_l(phi_l)            

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

          else if (k .eq. 2) then

            call qgaus(p22hha_s,dm,z0,p22a)

            call phivvb_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22a*conjg(phi_l)
         
         else if (k .eq. 3) then

            call qgaus(p22hha_s,dm,z0,p22a)

            call phivvc_l(phi_l)  

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

            call qgaus(p22hha_s,dm,z0,p22a)

            call phivvd_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22a*conjg(phi_l)            
        
         else

           print*, 'error in stalk_phil_ave e'

         endif

       else if (j .eq. 2)

          if (k .eq. 1) then

            call qgaus(p22hhb_s,dm,z0,p22b)

            call phivva_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22b*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p22hhb_s,dm,z0,p22b)

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

            call qgaus(p22hhb_s,dm,z0,p22b)

            call phivvc_l(phi_l)  

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

            call qgaus(p22hhb_s,dm,z0,p22b)

            call phivvd_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22b*conjg(phi_l)            
        
         else

           print*, 'error in stalk_phil_ave f'

         endif

       else if (j .eq. 3)

          if (k .eq. 1) then

            call qgaus(p22hhc_s,dm,z0,p22c)

            call phivva_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22c*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p22hhc_s,dm,z0,p22c)

            call phivvb_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22c*conjg(phi_l)
         
         else if (k .eq. 3) then

            call qgaus(p22hhc_s,dm,z0,p22c)

            call phivvc_l(phi_l)  

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

            call qgaus(p22hhc_s,dm,z0,p22c)

            call phivvd_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22c*conjg(phi_l)            
        
         else

           print*, 'error in stalk_phil_ave g'

         endif

       else if (j .eq. 4)

          if (k .eq. 1) then

            call qgaus(p22hhd_s,dm,z0,p22d)

            call phivva_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22d*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p22hhd_s,dm,z0,p22d)

            call phivvb_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22d*conjg(phi_l)
         
         else if (k .eq. 3) then

            call qgaus(p22hhd_s,dm,z0,p22d)

            call phivvc_l(phi_l)  

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

            call qgaus(p22hhd_s,dm,z0,p22d)

            call phivvd_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22d*conjg(phi_l)            
        
         else

           print*, 'error in stalk_phil_ave h'

         endif

      if (i .eq. 5) 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)

            call phihha_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
             *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (cos(theta_i)**2*p11a
     &       -2.0*cos(theta_i)*sin(theta_i)*p13a
     &       +sin(theta_i)**2*p33a)


          else if (k .eq. 2) then

            call qgaus(p11vva_s,dm,z0,p11a)

            call qgaus(p33vva_s,dm,z0,p33a)

            call qgaus(p13vva_s,dm,z0,p13a)

            call phihhb_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
             *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (cos(theta_i)**2*p11a
     &       -2.0*cos(theta_i)*sin(theta_i)*p13a
     &       +sin(theta_i)**2*p33a)
             
          else if (k .eq. 3) then

            call qgaus(p11vva_s,dm,z0,p11a)

            call qgaus(p33vva_s,dm,z0,p33a)

            call qgaus(p13vva_s,dm,z0,p13a)

            call phihhc_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
             *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (cos(theta_i)**2*p11a
     &       -2.0*cos(theta_i)*sin(theta_i)*p13a
     &       +sin(theta_i)**2*p33a)

          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 phihhd_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
             *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (cos(theta_i)**2*p11a
     &       -2.0*cos(theta_i)*sin(theta_i)*p13a
     &       +sin(theta_i)**2*p33a)
                      
          else

            print*, 'error in stalk_phil_ave a'

          endif

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


            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)

            call phihha_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (-cos(theta_i)**2*p11b*sin(theta_i)**2*p33b)
            
          else if (k .eq. 2) then


            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)

            call phihhb_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (-cos(theta_i)**2*p11b*sin(theta_i)**2*p33b)

          else if (k .eq. 3)


            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)

            call phihhc_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (-cos(theta_i)**2*p11b*sin(theta_i)**2*p33b)

          else if (k .eq. 4)

            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)

            call phihhd_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (-cos(theta_i)**2*p11b*sin(theta_i)**2*p33b)

          else

            print*, 'error in stalk_phil_ave b'

          endif

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


            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)

            call phihha_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (-cos(theta_i)**2*p11c*sin(theta_i)**2*p33c)
            
          else if (k .eq. 2) then


            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)

            call phihhc_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (-cos(theta_i)**2*p11c*sin(theta_i)**2*p33c)

          else if (k .eq. 3)


            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)

            call phihhc_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (-cos(theta_i)**2*p11c*sin(theta_i)**2*p33c)

          else if (k .eq. 4)

            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)

            call phihhd_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (-cos(theta_i)**2*p11c*sin(theta_i)**2*p33c)

          else

            print*, 'error in stalk_phil_ave c'

          endif

        else if (j .eq. 4) then

          if (k .eq. 1) then

            call qgaus(p11vvd_s,dm,z0,p11d)

            call qgaus(p33vvd_s,dm,z0,p33d)

            call qgaus(p13vvd_s,dm,z0,p13d)

            call phihhd_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
             *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (cos(theta_i)**2*p11d
     &       +2.0*cos(theta_i)*sin(theta_i)*p13d
     &       +sin(theta_i)**2*p33d)


          else if (k .eq. 2) then

            call qgaus(p11vvd_s,dm,z0,p11d)

            call qgaus(p33vvd_s,dm,z0,p33d)

            call qgaus(p13vvd_s,dm,z0,p13d)

            call phihhb_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
             *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (cos(theta_i)**2*p11d
     &       +2.0*cos(theta_i)*sin(theta_i)*p13d
     &       +sin(theta_i)**2*p33d)
             
          else if (k .eq. 3) then

            call qgaus(p11vvd_s,dm,z0,p11d)

            call qgaus(p33vvd_s,dm,z0,p33d)

            call qgaus(p13vvd_s,dm,z0,p13d)

            call phihhc_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
             *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (cos(theta_i)**2*p11d
     &       +2.0*cos(theta_i)*sin(theta_i)*p13d
     &       +sin(theta_i)**2*p33d)

          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)

            call phihhd_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
             *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)
     &       (cos(theta_i)**2*p11d
     &       +2.0*cos(theta_i)*sin(theta_i)*p13d
     &       +sin(theta_i)**2*p33d)
                      
          else

            print*, 'error in stalk_phil_ave d'

          endif

        else

          print*, 'error in stalk_phil_ave d'

        endif

      else

        print*, 'error in stalk_phil_ave s'

      endif

      return

      end
             

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


      subroutine phivva_l(r)


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

       rho=a(d+z_lc)
       
       x=rho*cos(alpha_s)

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

       b=ii*2.0*k_0*(sin(theta_i)*cos(alpha_s)*a+cos(theta_i))
       +2.0*mv_leaf/cos(theta_i)
        
       r=phia*(cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b


       return

       end

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


      subroutine phihha_l(r)


c************************************************************************
      
            
       rho=a(d+z_lc)
       
       x=rho*cos(alpha_s)

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

       b=ii*2.0*k_0*(sin(theta_i)*cos(alpha_s)*a+cos(theta_i))+
       2.0*mh_leaf/cos(theta_i)
        
       r=phia*(cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b


       return

       end

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


      subroutine phihva_l(r)


c************************************************************************
      
            
       rho=a(d+z_lc)
       
       x=rho*cos(alpha_s)

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

       b=ii*2.0*k_0*(sin(theta_i)*cos(alpha_s)*a+cos(theta_i))+
       (mv_leaf+mh_leaf)/cos(theta_i)
        
       r=phia*(cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b


       return

       end      

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


      subroutine phivvb_l(r)


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

       rho=a(d+z_lc)
       
       x=rho*cos(alpha_s)

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

       b=ii*2.0*k_0*sin(theta_i)*cos(alpha_s)*a
        
       r=r_v(theta_i)*phib*
       (cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b

       return

       end

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


      subroutine phihhb_l(r)


c************************************************************************
      
       rho=a(d+z_lc)
       
       x=rho*cos(alpha_s)

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

       b=ii*2.0*k_0*sin(theta_i)*cos(alpha_s)*a
        
       r=r_h(theta_i)*phib*
       (cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b

       return

       end
            
       return

       end

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


      subroutine phihvb_l(r)


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


       rho=a(d+z_lc)
       
       x=rho*cos(alpha_s)

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

       b=ii*2.0*k_0*sin(theta_i)*cos(alpha_s)*a
       +(mh_leaf-mv_leaf)/cos(theta_i)
        
       r=r_v(theta_i)*phib*
       (cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b

       return

       end
            
       return

       end

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


      subroutine phivvc_l(r)


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

       rho=a(d+z_lc)
       
       x=rho*cos(alpha_s)

       phic=cexp(ii*2.0*k_0*(sin(theta_i)*x))
       *cexp(-ii*2.0*k_0*cos(theta_i)*d+2.0*tau_v(-d)/cos(theta_i))

       b=ii*2.0*k_0*sin(theta_i)*cos(alpha_s)*a
        
       r=r_v(theta_i)*phic*
       (cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b

       return

       end

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


      subroutine phihhc_l(r)


c************************************************************************
      
       rho=a(d+z_lc)
       
       x=rho*cos(alpha_s)

       phic=cexp(ii*2.0*k_0*(sin(theta_i)*x))
       *cexp(-ii*2.0*k_0*cos(theta_i)*d+2.0*tau_h(-d)/cos(theta_i))

       b=ii*2.0*k_0*sin(theta_i)*cos(alpha_s)*a
        
       r=r_h(theta_i)*phic*
       (cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b

       return

       end
            
       return

       end

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


      subroutine phihvc_l(r)


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


       rho=a(d+z_lc)
       
       x=rho*cos(alpha_s)

       phic=cexp(ii*2.0*k_0*(sin(theta_i)*x)+(tau_v(z_lc)-tau_h(z_lc))
       /cos(theta_i))
       *cexp(-ii*2.0*k_0*cos(theta_i)*d+2.0*tau_h(-d)/cos(theta_i))

       b=ii*2.0*k_0*sin(theta_i)*cos(alpha_s)*a
       +(mh_leaf-mv_leaf)/cos(theta_i)
        
       r=r_h(theta_i)*phic*
       (cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b

       return

       end
            

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


      subroutine phivvd_l(r)


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

       rho=a(d+z_lc)
       
       x=rho*cos(alpha_s)

       phid=cexp(ii*2.0*k_0*(sin(theta_i)*x-cos(theta_i)*z_lc)
     & - 2.0*tau_v(z_lc)/cos(theta_i))
       *cexp(-ii*4.0*k_0*cos(theta_i)*d+4.0*tau_v(-d)/cos(theta_i))
       
       b=ii*2.0*k_0*(sin(theta_i)*cos(alpha_s)*a-cos(theta_i))
       -2.0*mv_leaf/cos(theta_i)
        
       r=phid*(cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b
       *r_v(theta_i)**2

       return

       end

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


      subroutine phihhd_l(r)


c************************************************************************
      
            
       rho=a(d+z_lc)
       
       x=rho*cos(alpha_s)

       phid=cexp(ii*2.0*k_0*(sin(theta_i)*x-cos(theta_i)*z_lc)
     & - 2.0*tau_h(z_lc)/cos(theta_i))
       *cexp(-ii*4.0*k_0*cos(theta_i)*d+4.0*tau_h(-d)/cos(theta_i))
       
       b=ii*2.0*k_0*(sin(theta_i)*cos(alpha_s)*a-cos(theta_i))-
       2.0*mh_leaf/cos(theta_i)
        
       r=phid*(cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b
       *r_h(theta_i)**2

       return

       end

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


      subroutine phihvd_l(r)


c************************************************************************
      
            
       rho=a(d+z_lc)
       
       x=rho*cos(alpha_s)

       phid=cexp(ii*2.0*k_0*(sin(theta_i)*x-cos(theta_i)*z_lc)
     & - (tau_v(z_lc)+tau_h(z_lc))/cos(theta_i))
       *cexp(-ii*4.0*k_0*cos(theta_i)*d+2.0*(tau_v(-d)+tau_h(-d))
       /cos(theta_i))
       
       b=ii*2.0*k_0*(sin(theta_i)*cos(alpha_s)*a-cos(theta_i))-
       (mv_leaf+mh_leaf)/cos(theta_i)
        
       r=phid*(cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b
       *r_v(theta_i)*r_h(theta_i)

       return

       end




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

      subroutine leaf_field_int(m,n,p, s)

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

      common /mn/ m,n,p


      include 'constants.include'

      external h_l
      call qgausx(h_l,pi,mpi,r)

      return

      end
      

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

      function h_l(phii)

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

      external g_l
      common /stalk_stat/ amin,amax
      common /phiv/ phi,a,z0

      phi=phii

      call qgausy(g_l,vmin,vmax,s)

      h_l=s(m,n,p)

      return

      end

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

      function g_l(vv)

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



      external stalk_phil_ave

      common /phiv/ phi,v
      common /mn/  m,n,p

      v=vv
      
      g_l=leaf_field_ave(phi,v,m,n)

      return

      end

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

      function leaf_field_ave(alpha,s,i,j)

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

      external p11vva,p33vva,p13vva
      
      common /incident/ theta_i,phi_i
      

      if (i .eq. 1) then

        if (j .eq. 1) then


            call qgaus(p11vva_l,zero,p0,p11a)

            call qgaus(p33vva_l,zero,p0,p33a)


            leaf_field_ave=leaf_alpha_pdf(alpha)
             *leaf_v_pdf(v)*
     &       (cos(theta_i)**2*p11a
     &       +sin(theta_i)**2*p33a)


        else if (j .eq. 2) then
          

            call qgaus(p11vvb_l,zero,p0,p11b)

            call qgaus(p33vvb_l,zero,p0,p33b)


            leaf_field_ave=leaf_alpha_pdf(alpha_l)
     &       *leaf_v_pdf(v)* 
     &       (-cos(theta_i)**2*p11b+sin(theta_i)**2*p33b)
            

        else if (j .eq. 3) then
          

            call qgaus(p11vvc_l,zero,p0,p11c)

            call qgaus(p33vvc_l,zero,p0,p33c)


            leaf_field_ave=leaf_alpha_pdf(alpha_l)
     &       *leaf_v_pdf(v)* 
     &       (-cos(theta_i)**2*p11c+sin(theta_i)**2*p33c)
            

        else if (j .eq. 4) then


            call qgaus(p11vvd_l,zero,p0,p11d)

            call qgaus(p33vvd_l,zero,p0,p33d)



            leaf_field_ave=leaf_alpha_pdf(alpha)
             *leaf_v_pdf(v)* 
     &       (cos(theta_i)**2*p11d
     &       +sin(theta_i)**2*p33d)


         else
           
            print*, 'error in leaf_field_ave d'

        end if

      else if (i .eq. 2) then

        if (j .eq. 1) then


            call qgaus(p22hha_l,zero,p0,p22a)


            leaf_field_ave=leaf_alpha_pdf(alpha_l)
     &       *leaf_v_pdf(v)*p22a 

       else if (j .eq. 2)


            call qgaus(p22hhb_l,zero,p0,p22b)


            leaf_field_ave=leaf_alpha_pdf(alpha_l)
     &       *leaf_v_pdf(v)*p22b 


       else if (j .eq. 3)


            call qgaus(p22hhc_l,zero,p0,p22c)

            leaf_field_ave=leaf_alpha_pdf(alpha_l)
     &       *leaf_v_pdf(v)*p22c 


       else if (j .eq. 4)


            call qgaus(p22hhd_l,zero,p0,p22d)

            leaf_field_ave=leaf_alpha_pdf(alpha_l)
     &       *leaf_v_pdf(v)*p22d 

         else

           print*, 'error in leaf_field_ave h'

         endif

       else if (i .eq. 3) then

         leaf_field_ave=cmplx(0.0,0.0)


      else if (i .eq. 4) then

        if (j .eq. 1) then


            call qgaus(p22hha_l,zero,p0,p22a)


            leaf_field_ave=leaf_alpha_pdf(alpha_l)
     &       *leaf_v_pdf(v)*p22a 

       else if (j .eq. 2)


            call qgaus(p22hhb_l,zero,p0,p22b)


            leaf_field_ave=leaf_alpha_pdf(alpha_l)
     &       *leaf_v_pdf(v)*p22b 


       else if (j .eq. 3)


            call qgaus(p22hhc_l,zero,p0,p22c)

            leaf_field_ave=leaf_alpha_pdf(alpha_l)
     &       *leaf_v_pdf(v)*p22c 


       else if (j .eq. 4)


            call qgaus(p22hhd_l,zero,p0,p22d)

            leaf_field_ave=leaf_alpha_pdf(alpha_l)
     &       *leaf_v_pdf(v)*p22d 

         else

           print*, 'error in leaf_field_ave h'

         endif


      if (i .eq. 5) then

        if (j .eq. 1) then


            call qgaus(p11vva_l,zero,p0,p11a)

            call qgaus(p33vva_l,zero,p0,p33a)

            leaf_field_ave=leaf_alpha_pdf(alpha)
             *leaf_v_pdf(v)*
     &       (cos(theta_i)**2*p11a
     &       +sin(theta_i)**2*p33a)


        else if (j .eq. 2) then
          

            call qgaus(p11vvb_l,zero,p0,p11b)

            call qgaus(p33vvb_l,zero,p0,p33b)


            leaf_field_ave=leaf_alpha_pdf(alpha_l)
     &       *leaf_v_pdf(v)* 
     &       (-cos(theta_i)**2*p11b+sin(theta_i)**2*p33b)
            

        else if (j .eq. 3) then
          

            call qgaus(p11vvc_l,zero,p0,p11c)

            call qgaus(p33vvc_l,zero,p0,p33c)


            leaf_field_ave=leaf_alpha_pdf(alpha_l)
     &       *leaf_v_pdf(v)* 
     &       (-cos(theta_i)**2*p11c+sin(theta_i)**2*p33c)
            

        else if (j .eq. 4) then


            call qgaus(p11vvd_l,zero,p0,p11d)

            call qgaus(p33vvd_l,zero,p0,p33d)



            leaf_field_ave=leaf_alpha_pdf(alpha)
             *leaf_v_pdf(v)* 
     &       (cos(theta_i)**2*p11d
     &       +sin(theta_i)**2*p33d)


         else
           
            print*, 'error in leaf_field_ave d'

        end if

      else

        print*, 'error in leaf_field_ave s'

      endif

      return

      end      


      
 
