
      
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      


      
 
