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

      complex function leaf_scat(phi,v,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_l,pzzvva_l,pxzvva_l
      complex pxxvvb_l,pzzvvb_l,pxzvvb_l
      complex pxxvvc_l,pzzvvc_l,pxzvvc_l
      complex pxxvvd_l,pzzvvd_l,pxzvvd_l

      complex pyyhha_l
      complex pyyhhb_l
      complex pyyhhc_l
      complex pyyhhd_l

      complex pxyhva_l,pyzhva_l
      complex pxyhvb_l,pyzhvb_l
      complex pxyhvc_l,pyzhvc_l
      complex pxyhvd_l,pyzhvd_l

      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 pxxva,pxxvb,pxxvc,pxxvd
      complex pzzva,pzzvb,pzzvc,pzzvd
      complex pyyha,pyyhb,pyyhc,pyyhd
      complex pxxm,pzzm,crss,pxzm
      
      real c1,c2,alpha,ell_0
      real theta_i,phi_i,k_0
      real zero,rho_0,v,phi
      real pi,coef

      external pxxvva_l,pzzvva_l,pxzvva_l
      external pxxvvb_l,pzzvvb_l,pxzvvb_l
      external pxxvvc_l,pzzvvc_l,pxzvvc_l
      external pxxvvd_l,pzzvvd_l,pxzvvd_l

      external pyyhha_l
      external pyyhhb_l
      external pyyhhc_l
      external pyyhhd_l

      external pxyhva_l,pyzhva_l
      external pxyhvb_l,pyzhvb_l
      external pxyhvc_l,pyzhvc_l
      external pxyhvd_l,pyzhvd_l      
      
      common /leaf_struct/ c1,c2,alpha,ell_0
      common /incident/ theta_i,phi_i,k_0

      integer l_tensor_flag,l_tau_flag,l_field_flag(2,4)
      common /leaf_flags/ l_tensor_flag,l_tau_flag,l_field_flag

      include 'error.incld'
      
      zero=0.0
      pi=3.1415927
      coef=k_0**4/(16.0*pi**2)

      call leaf_model(v, c1,c2,rho_0,ell_0)
      
      alpha=phi
      
      if (pol .eq. vvvv) then
        
        if (mech .eq. p) then
          
          if (mech_c .eq. p) then
            
            l_tensor_flag=0
            call integral(pxxvva_l,zero,rho_0,pxxa)

            l_tensor_flag=0
            call integral(pzzvva_l,zero,rho_0,pzza)

            l_tensor_flag=0
            call integral(pxzvva_l,zero,rho_0,pxza)

            leaf_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))

            pxxm=pxxa*conjg(pxxa)
            pxzm=4.0*pxza*conjg(pxza)
            crss=(
     &       pxxa*conjg(pzza)+pzza*conjg(pxxa))
            pzzm=pzza*conjg(pzza)
            
c            print*,'alpha,pxxm,pzzm,pxzm,crss,leaf_scat'
c            print*,alpha,pxxm,pzzm,pxzm,crss,leaf_scat
            
          else if (mech_c .eq. gp) then

            l_tensor_flag=0
            call integral(pxxvva_l,zero,rho_0,pxxa)

            l_tensor_flag=0
            call integral(pzzvva_l,zero,rho_0,pzza)

            l_tensor_flag=0
            call integral(pxxvvb_l,zero,rho_0,pxxb)

            l_tensor_flag=0
            call integral(pzzvvb_l,zero,rho_0,pzzb)

            leaf_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

            l_tensor_flag=0
            call integral(pxxvva_l,zero,rho_0,pxxa)

            l_tensor_flag=0
            call integral(pzzvva_l,zero,rho_0,pzza)

            l_tensor_flag=0
            call integral(pxxvvc_l,zero,rho_0,pxxc)

            l_tensor_flag=0
            call integral(pzzvvc_l,zero,rho_0,pzzc)

            
            leaf_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

            l_tensor_flag=0
            call integral(pxxvva_l,zero,rho_0,pxxa)

            l_tensor_flag=0
            call integral(pzzvva_l,zero,rho_0,pzza)

            l_tensor_flag=0
            call integral(pxzvva_l,zero,rho_0,pxza)

            l_tensor_flag=0
            call integral(pxxvvd_l,zero,rho_0,pxxd)

            l_tensor_flag=0
            call integral(pzzvvd_l,zero,rho_0,pzzd)

            l_tensor_flag=0
            call integral(pxzvvd_l,zero,rho_0,pxzd)

            leaf_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=34
            call the_death_code(die,error,code)

          endif

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


            l_tensor_flag=0
            call integral(pxxvva_l,zero,rho_0,pxxa)

            l_tensor_flag=0
            call integral(pzzvva_l,zero,rho_0,pzza)

            l_tensor_flag=0
            call integral(pxxvvb_l,zero,rho_0,pxxb)

            l_tensor_flag=0
            call integral(pzzvvb_l,zero,rho_0,pzzb)


            leaf_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


            l_tensor_flag=0
            call integral(pxxvvb_l,zero,rho_0,pxxb)

            l_tensor_flag=0
            call integral(pzzvvb_l,zero,rho_0,pzzb)


            leaf_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

            l_tensor_flag=0
            call integral(pxxvvb_l,zero,rho_0,pxxb)

            l_tensor_flag=0
            call integral(pzzvvb_l,zero,rho_0,pzzb)

            l_tensor_flag=0
            call integral(pxxvvc_l,zero,rho_0,pxxc)

            l_tensor_flag=0
            call integral(pzzvvc_l,zero,rho_0,pzzc)
            
            leaf_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

            l_tensor_flag=0
            call integral(pxxvvb_l,zero,rho_0,pxxb)

            l_tensor_flag=0
            call integral(pzzvvb_l,zero,rho_0,pzzb)

            l_tensor_flag=0
            call integral(pxxvvd_l,zero,rho_0,pxxd)

            l_tensor_flag=0
            call integral(pzzvvd_l,zero,rho_0,pzzd)
            
            leaf_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=34
            call the_death_code(die,error,code)

          endif

        else if(mech .eq. pg) then

          if (mech_c .eq. p) then


            l_tensor_flag=0
            call integral(pxxvva_l,zero,rho_0,pxxa)

            l_tensor_flag=0
            call integral(pzzvva_l,zero,rho_0,pzza)

            l_tensor_flag=0
            call integral(pxxvvc_l,zero,rho_0,pxxc)

            l_tensor_flag=0
            call integral(pzzvvc_l,zero,rho_0,pzzc)


            leaf_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


            l_tensor_flag=0
            call integral(pxxvvb_l,zero,rho_0,pxxb)

            l_tensor_flag=0
            call integral(pzzvvb_l,zero,rho_0,pzzb)

            l_tensor_flag=0
            call integral(pxxvvc_l,zero,rho_0,pxxc)

            l_tensor_flag=0
            call integral(pzzvvc_l,zero,rho_0,pzzc)
            
            leaf_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

            l_tensor_flag=0
            call integral(pxxvvc_l,zero,rho_0,pxxc)

            l_tensor_flag=0
            call integral(pzzvvc_l,zero,rho_0,pzzc)
            
            leaf_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

            l_tensor_flag=0
            call integral(pxxvvc_l,zero,rho_0,pxxc)

            l_tensor_flag=0
            call integral(pzzvvc_l,zero,rho_0,pzzc)

            l_tensor_flag=0
            call integral(pxxvvd_l,zero,rho_0,pxxd)

            l_tensor_flag=0
            call integral(pzzvvd_l,zero,rho_0,pzzd)
            
            leaf_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=34
            call the_death_code(die,error,code)

          endif

        else if (mech .eq. gpg) then

          if (mech_c .eq. p) then

            l_tensor_flag=0
            call integral(pxxvva_l,zero,rho_0,pxxa)

            l_tensor_flag=0
            call integral(pzzvva_l,zero,rho_0,pzza)

            l_tensor_flag=0
            call integral(pxzvva_l,zero,rho_0,pxza)
            
            l_tensor_flag=0
            call integral(pxxvvd_l,zero,rho_0,pxxd)

            l_tensor_flag=0
            call integral(pzzvvd_l,zero,rho_0,pzzd)

            l_tensor_flag=0
            call integral(pxzvvd_l,zero,rho_0,pxzd)

            leaf_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

            l_tensor_flag=0
            call integral(pxxvvb_l,zero,rho_0,pxxb)

            l_tensor_flag=0
            call integral(pzzvvb_l,zero,rho_0,pzzb)

            l_tensor_flag=0
            call integral(pxxvvd_l,zero,rho_0,pxxd)

            l_tensor_flag=0
            call integral(pzzvvd_l,zero,rho_0,pzzd)

            leaf_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

            l_tensor_flag=0
            call integral(pxxvvd_l,zero,rho_0,pxxd)

            l_tensor_flag=0
            call integral(pzzvvd_l,zero,rho_0,pzzd)

            l_tensor_flag=0
            call integral(pxxvvc_l,zero,rho_0,pxxc)
            
            l_tensor_flag=0
            call integral(pzzvvc_l,zero,rho_0,pzzc)

            leaf_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

            l_tensor_flag=0
            call integral(pxxvvd_l,zero,rho_0,pxxd)

            l_tensor_flag=0
            call integral(pzzvvd_l,zero,rho_0,pzzd)

            l_tensor_flag=0
            call integral(pxzvvd_l,zero,rho_0,pxzd)

            leaf_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=34
            call the_death_code(die,error,code)
            
          endif

        else
            

            write(error,*)'mech = ',mech,' not recognized in',
     &           ' func leaf_scat'
            code=35
            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

            l_tensor_flag=0
            call integral(pyyhha_l,zero,rho_0,pyya)

            leaf_scat=coef*
     &       pyya*conjg(pyya)

          else if (mech_c .eq. gp) then

            l_tensor_flag=0
            call integral(pyyhha_l,zero,rho_0,pyya)

            l_tensor_flag=0
            call integral(pyyhhb_l,zero,rho_0,pyyb)            

            leaf_scat=coef*
     &       pyya*conjg(pyyb)
         
         else if (mech_c .eq. pg) then

            l_tensor_flag=0
            call integral(pyyhha_l,zero,rho_0,pyya)

            l_tensor_flag=0
            call integral(pyyhhc_l,zero,rho_0,pyyc)            

            leaf_scat=coef*
     &       pyya*conjg(pyyc)
                     
         else if (mech_c .eq. gpg) then

            l_tensor_flag=0
            call integral(pyyhha_l,zero,rho_0,pyya)

            l_tensor_flag=0
            call integral(pyyhhd_l,zero,rho_0,pyyd)            

            leaf_scat=coef*
     &       pyya*conjg(pyyd)            
        
         else


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

          endif


       else if (mech .eq. gp) then

         if (mech_c .eq. p) then

            l_tensor_flag=0
            call integral(pyyhhb_l,zero,rho_0,pyyb)            

            l_tensor_flag=0
            call integral(pyyhha_l,zero,rho_0,pyya)

            leaf_scat=coef*
     &       pyyb*conjg(pyya)

          else if (mech_c .eq. gp) then

            l_tensor_flag=0
            call integral(pyyhhb_l,zero,rho_0,pyyb)            

            leaf_scat=coef*
     &       pyyb*conjg(pyyb)

         
         else if (mech_c .eq. pg) then

            l_tensor_flag=0
            call integral(pyyhhb_l,zero,rho_0,pyyb)

            l_tensor_flag=0
            call integral(pyyhhc_l,zero,rho_0,pyyc)            

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

            l_tensor_flag=0
            call integral(pyyhhb_l,zero,rho_0,pyyb)

            l_tensor_flag=0
            call integral(pyyhhd_l,zero,rho_0,pyyd)            

            leaf_scat=coef*
     &       pyyb*conjg(pyyd)            
        
         else

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

         endif

       else if (mech .eq. pg) then

         if (mech_c .eq. p) then

            l_tensor_flag=0
            call integral(pyyhhc_l,zero,rho_0,pyyc)            

            l_tensor_flag=0
            call integral(pyyhha_l,zero,rho_0,pyya)

            leaf_scat=coef*
     &       pyyc*conjg(pyya)

          else if (mech_c .eq. gp) then

            l_tensor_flag=0
            call integral(pyyhhc_l,zero,rho_0,pyyc)

            l_tensor_flag=0
            call integral(pyyhhb_l,zero,rho_0,pyyb)            

            leaf_scat=coef*
     &       pyyc*conjg(pyyb)
         
         else if (mech_c .eq. pg) then

            l_tensor_flag=0
            call integral(pyyhhc_l,zero,rho_0,pyyc)            

            leaf_scat=coef*
     &       pyyc*conjg(pyyc)
                     
         else if (mech_c .eq. gpg) then

            l_tensor_flag=0
            call integral(pyyhhc_l,zero,rho_0,pyyc)

            l_tensor_flag=0
            call integral(pyyhhd_l,zero,rho_0,pyyd)            

            leaf_scat=coef*
     &       pyyc*conjg(pyyd)            
        
         else

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

         endif

       else if (mech .eq. gpg) then

         if (mech_c .eq. p) then

            l_tensor_flag=0
            call integral(pyyhhd_l,zero,rho_0,pyyd)            

            l_tensor_flag=0
            call integral(pyyhha_l,zero,rho_0,pyya)

            leaf_scat=coef*
     &       pyyd*conjg(pyya)

          else if (mech_c .eq. gp) then

            l_tensor_flag=0
            call integral(pyyhhd_l,zero,rho_0,pyyd)

            l_tensor_flag=0
            call integral(pyyhhb_l,zero,rho_0,pyyb)            

            leaf_scat=coef*
     &       pyyd*conjg(pyyb)
         
         else if (mech_c .eq. pg) then

           l_tensor_flag=0
           call integral(pyyhhd_l,zero,rho_0,pyyd)            

           l_tensor_flag=0
           call integral(pyyhhc_l,zero,rho_0,pyyc)            

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

            l_tensor_flag=0
            call integral(pyyhhd_l,zero,rho_0,pyyd)            

            leaf_scat=coef*
     &       pyyd*conjg(pyyd)            
        
         else

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

          endif

        else
          
          write(error,*)'mech = ',mech,' not recognized in',
     &         ' func leaf_scat'
          code=35
          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

            l_tensor_flag=0
            call integral(pxyhva_l,zero,rho_0,pxya)

            l_tensor_flag=0
            call integral(pyzhva_l,zero,rho_0,pyza)
            
            leaf_scat=coef*
     &       (cos(theta_i)**2*pxya*conjg(pxya)
     &       +sin(theta_i)**2*pyza*conjg(pyza))


           else if (mech_c .eq. gp) then

            l_tensor_flag=0
            call integral(pxyhva_l,zero,rho_0,pxya)

            l_tensor_flag=0
            call integral(pyzhva_l,zero,rho_0,pyza)
            
            l_tensor_flag=0
            call integral(pxyhvb_l,zero,rho_0,pxyb)

            l_tensor_flag=0
            call integral(pyzhvb_l,zero,rho_0,pyzb)

            leaf_scat=coef*
     &       (-cos(theta_i)**2*pxya*conjg(pxyb)
     &       +sin(theta_i)**2*pyza*conjg(pyzb))


          else if (mech_c .eq. pg) then

            l_tensor_flag=0
            call integral(pxyhva_l,zero,rho_0,pxya)

            l_tensor_flag=0
            call integral(pyzhva_l,zero,rho_0,pyza)
            
            l_tensor_flag=0
            call integral(pxyhvc_l,zero,rho_0,pxyc)

            l_tensor_flag=0
            call integral(pyzhvc_l,zero,rho_0,pyzc)

            leaf_scat=coef*
     &       (cos(theta_i)**2*pxya*conjg(pxyc)
     &       +sin(theta_i)**2*pyza*conjg(pyzc))


          else if (mech_c .eq. gpg) then

            l_tensor_flag=0
            call integral(pxyhva_l,zero,rho_0,pxya)

            l_tensor_flag=0
            call integral(pyzhva_l,zero,rho_0,pyza)
            
            l_tensor_flag=0
            call integral(pxyhvd_l,zero,rho_0,pxyd)

            l_tensor_flag=0
            call integral(pyzhvd_l,zero,rho_0,pyzd)

            leaf_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=34
            call the_death_code(die,error,code)

          endif

        else if (mech .eq. gp) then

          if (mech_c .eq. p) then

            l_tensor_flag=0
            call integral(pxyhvb_l,zero,rho_0,pxyb)

            l_tensor_flag=0
            call integral(pyzhvb_l,zero,rho_0,pyzb)

            l_tensor_flag=0
            call integral(pxyhva_l,zero,rho_0,pxya)

            l_tensor_flag=0
            call integral(pyzhva_l,zero,rho_0,pyza)
            
            leaf_scat=coef*
     &       (-cos(theta_i)**2*pxyb*conjg(pxya)
     &       +sin(theta_i)**2*pyzb*conjg(pyza))

            

          else if (mech_c .eq. gp) then

            l_tensor_flag=0
            call integral(pxyhvb_l,zero,rho_0,pxyb)

            l_tensor_flag=0
            call integral(pyzhvb_l,zero,rho_0,pyzb)

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


           else if (mech_c .eq. pg) then

            l_tensor_flag=0
            call integral(pxyhvb_l,zero,rho_0,pxyb)

            l_tensor_flag=0
            call integral(pyzhvb_l,zero,rho_0,pyzb)
            
            l_tensor_flag=0
            call integral(pxyhvc_l,zero,rho_0,pxyc)

            l_tensor_flag=0
            call integral(pyzhvc_l,zero,rho_0,pyzc)

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


           else if (mech_c .eq. gpg) then

            l_tensor_flag=0
            call integral(pxyhvb_l,zero,rho_0,pxyb)

            l_tensor_flag=0
            call integral(pyzhvb_l,zero,rho_0,pyzb)
            
            l_tensor_flag=0
            call integral(pxyhvd_l,zero,rho_0,pxyd)

            l_tensor_flag=0
            call integral(pyzhvd_l,zero,rho_0,pyzd)

            leaf_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=34
            call the_death_code(die,error,code)

          endif

         else if (mech .eq. pg) then

           if (mech_c .eq. p) then

            l_tensor_flag=0
            call integral(pxyhvc_l,zero,rho_0,pxyc)

            l_tensor_flag=0
            call integral(pyzhvc_l,zero,rho_0,pyzc)

            l_tensor_flag=0
            call integral(pxyhva_l,zero,rho_0,pxya)

            l_tensor_flag=0
            call integral(pyzhva_l,zero,rho_0,pyza)
            
            leaf_scat=coef*
     &       (cos(theta_i)**2*pxyc*conjg(pxya)
     &       +sin(theta_i)**2*pyzc*conjg(pyza))

            
           else if (mech_c .eq. gp) then

            l_tensor_flag=0
            call integral(pxyhvc_l,zero,rho_0,pxyc)

            l_tensor_flag=0
            call integral(pyzhvc_l,zero,rho_0,pyzc)
            
            l_tensor_flag=0
            call integral(pxyhvb_l,zero,rho_0,pxyb)

            l_tensor_flag=0
            call integral(pyzhvb_l,zero,rho_0,pyzb)

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

           else if (mech_c .eq. pg) then

            l_tensor_flag=0
            call integral(pxyhvc_l,zero,rho_0,pxyc)

            l_tensor_flag=0
            call integral(pyzhvc_l,zero,rho_0,pyzc)

            leaf_scat=coef*
     &       (cos(theta_i)**2*pxyc*conjg(pxyc)
     &       +sin(theta_i)**2*pyzc*conjg(pyzc))

           else if (mech_c .eq. gpg) then

            l_tensor_flag=0
            call integral(pxyhvc_l,zero,rho_0,pxyc)

            l_tensor_flag=0
            call integral(pyzhvc_l,zero,rho_0,pyzc)
            
            l_tensor_flag=0
            call integral(pxyhvd_l,zero,rho_0,pxyd)

            l_tensor_flag=0
            call integral(pyzhvd_l,zero,rho_0,pyzd)

            leaf_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=34
            call the_death_code(die,error,code)

          endif


         else if (mech .eq. gpg) then

           if (mech_c .eq. p) then

            l_tensor_flag=0
            call integral(pxyhvd_l,zero,rho_0,pxyd)

            l_tensor_flag=0
            call integral(pyzhvd_l,zero,rho_0,pyzd)

            l_tensor_flag=0
            call integral(pxyhva_l,zero,rho_0,pxya)

            l_tensor_flag=0
            call integral(pyzhva_l,zero,rho_0,pyza)
            
            leaf_scat=coef*
     &       (-cos(theta_i)**2*pxyd*conjg(pxya)
     &       +sin(theta_i)**2*pyzd*conjg(pyza))

            
           else if (mech_c .eq. gp) then

            l_tensor_flag=0
            call integral(pxyhvd_l,zero,rho_0,pxyd)

            l_tensor_flag=0
            call integral(pyzhvd_l,zero,rho_0,pyzd)
            
            l_tensor_flag=0
            call integral(pxyhvb_l,zero,rho_0,pxyb)

            l_tensor_flag=0
            call integral(pyzhvb_l,zero,rho_0,pyzb)

            leaf_scat=coef*
     &           (cos(theta_i)**2*pxyd*conjg(pxyb)
     &           +sin(theta_i)**2*pyzd*conjg(pyzb))

           else if (mech_c .eq. pg) then

            l_tensor_flag=0
            call integral(pxyhvd_l,zero,rho_0,pxyd)

            l_tensor_flag=0
            call integral(pyzhvd_l,zero,rho_0,pyzd)
            
            l_tensor_flag=0
            call integral(pxyhvc_l,zero,rho_0,pxyc)

            l_tensor_flag=0
            call integral(pyzhvc_l,zero,rho_0,pyzc)

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

           else if (mech_c .eq. gpg) then

            l_tensor_flag=0
            call integral(pxyhvd_l,zero,rho_0,pxyd)

            l_tensor_flag=0
            call integral(pyzhvd_l,zero,rho_0,pyzd)

            leaf_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=34
            call the_death_code(die,error,code)

          endif

        else

          write(error,*)'mech = ',mech,' not recognized in',
     &         ' func leaf_scat'
          code=35
          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

            l_tensor_flag=0
            call integral(pyyhha_l,zero,rho_0,pyyha)

            l_tensor_flag=0
            call integral(pxxvva_l,zero,rho_0,pxxva)

            l_tensor_flag=0
            call integral(pzzvva_l,zero,rho_0,pzzva)

            leaf_scat=coef*
     &       (cos(theta_i)**2*pyyha*conjg(pxxva)
     &       +sin(theta_i)**2*pyyha*conjg(pzzva))

          else if (mech_c .eq. gp) then

            l_tensor_flag=0
            call integral(pyyhha_l,zero,rho_0,pyyha)

            l_tensor_flag=0
            call integral(pxxvvb_l,zero,rho_0,pxxvb)

            l_tensor_flag=0
            call integral(pzzvvb_l,zero,rho_0,pzzvb)
            
            leaf_scat=coef*
     &           (cos(theta_i)**2*pyyha*conjg(pxxvb)
     &           +sin(theta_i)**2*pyyha*conjg(pzzvb))

          else if (mech_c .eq. pg) then

            l_tensor_flag=0
            call integral(pyyhha_l,zero,rho_0,pyyha)

            l_tensor_flag=0
            call integral(pxxvvc_l,zero,rho_0,pxxvc)

            l_tensor_flag=0
            call integral(pzzvvc_l,zero,rho_0,pzzvc)
            
            leaf_scat=coef*
     &           (cos(theta_i)**2*pyyha*conjg(pxxvc)
     &           +sin(theta_i)**2*pyyha*conjg(pzzvc))

          else if (mech_c .eq. gpg) then

            l_tensor_flag=0
            call integral(pyyhha_l,zero,rho_0,pyyha)

            l_tensor_flag=0
            call integral(pxxvvd_l,zero,rho_0,pxxvd)

            l_tensor_flag=0
            call integral(pzzvvd_l,zero,rho_0,pzzvd)
            
            leaf_scat=coef*
     &           (cos(theta_i)**2*pyyha*conjg(pxxvd)
     &           +sin(theta_i)**2*pyyha*conjg(pzzvd))

          else

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

          endif

        else if (mech .eq. gp) then

          if (mech_c .eq. p) then

            l_tensor_flag=0
            call integral(pyyhhb_l,zero,rho_0,pyyhb)

            l_tensor_flag=0
            call integral(pxxvva_l,zero,rho_0,pxxva)

            l_tensor_flag=0
            call integral(pzzvva_l,zero,rho_0,pzzva)

            leaf_scat=coef*
     &       (-cos(theta_i)**2*pyyhb*conjg(pxxva)
     &       +sin(theta_i)**2*pyyhb*conjg(pzzva))

          else if (mech_c .eq. gp) then

            l_tensor_flag=0
            call integral(pyyhhb_l,zero,rho_0,pyyhb)

            l_tensor_flag=0
            call integral(pxxvvb_l,zero,rho_0,pxxvb)

            l_tensor_flag=0
            call integral(pzzvvb_l,zero,rho_0,pzzvb)
            
            leaf_scat=coef*
     &       (-cos(theta_i)**2*pyyhb*conjg(pxxvb)
     &       +sin(theta_i)**2*pyyhb*conjg(pzzvb))

          else if (mech_c .eq. pg) then

            l_tensor_flag=0
            call integral(pyyhhb_l,zero,rho_0,pyyhb)

            l_tensor_flag=0
            call integral(pxxvvc_l,zero,rho_0,pxxvc)

            l_tensor_flag=0
            call integral(pzzvvc_l,zero,rho_0,pzzvc)
            
            leaf_scat=coef*
     &       (-cos(theta_i)**2*pyyhb*conjg(pxxvc)
     &       +sin(theta_i)**2*pyyhb*conjg(pzzvc))

          else if (mech_c .eq. gpg) then

            l_tensor_flag=0
            call integral(pyyhhb_l,zero,rho_0,pyyhb)

            l_tensor_flag=0
            call integral(pxxvvd_l,zero,rho_0,pxxvd)

            l_tensor_flag=0
            call integral(pzzvvd_l,zero,rho_0,pzzvd)
            
            leaf_scat=coef*
     &       (-cos(theta_i)**2*pyyhb*conjg(pxxvd)
     &       +sin(theta_i)**2*pyyhb*conjg(pzzvd))

          else

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

          endif


        else if (mech .eq. pg) then

          if (mech_c .eq. p) then

            l_tensor_flag=0
            call integral(pyyhhc_l,zero,rho_0,pyyhc)

            l_tensor_flag=0
            call integral(pxxvva_l,zero,rho_0,pxxva)

            l_tensor_flag=0
            call integral(pzzvva_l,zero,rho_0,pzzva)

            leaf_scat=coef*
     &       (-cos(theta_i)**2*pyyhc*conjg(pxxva)
     &       +sin(theta_i)**2*pyyhc*conjg(pzzva))

          else if (mech_c .eq. gp) then

            l_tensor_flag=0
            call integral(pyyhhc_l,zero,rho_0,pyyhc)

            l_tensor_flag=0
            call integral(pxxvvb_l,zero,rho_0,pxxvb)

            l_tensor_flag=0
            call integral(pzzvvb_l,zero,rho_0,pzzvb)
            
            leaf_scat=coef*
     &       (-cos(theta_i)**2*pyyhc*conjg(pxxvb)
     &       +sin(theta_i)**2*pyyhc*conjg(pzzvb))

          else if (mech_c .eq. pg) then

            l_tensor_flag=0
            call integral(pyyhhc_l,zero,rho_0,pyyhc)

            l_tensor_flag=0
            call integral(pxxvvc_l,zero,rho_0,pxxvc)

            l_tensor_flag=0
            call integral(pzzvvc_l,zero,rho_0,pzzvc)
            
            leaf_scat=coef*
     &       (-cos(theta_i)**2*pyyhc*conjg(pxxvc)
     &       +sin(theta_i)**2*pyyhc*conjg(pzzvc))

          else if (mech_c .eq. gpg) then

            l_tensor_flag=0
            call integral(pyyhhc_l,zero,rho_0,pyyhc)

            l_tensor_flag=0
            call integral(pxxvvd_l,zero,rho_0,pxxvd)

            l_tensor_flag=0
            call integral(pzzvvd_l,zero,rho_0,pzzvd)
            
            leaf_scat=coef*
     &       (-cos(theta_i)**2*pyyhc*conjg(pxxvd)
     &       +sin(theta_i)**2*pyyhc*conjg(pzzvd))

          else

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

          endif

        else if (mech .eq. gpg) then

          if (mech_c .eq. p) then

            l_tensor_flag=0
            call integral(pyyhhd_l,zero,rho_0,pyyhd)

            l_tensor_flag=0
            call integral(pxxvva_l,zero,rho_0,pxxva)

            l_tensor_flag=0
            call integral(pzzvva_l,zero,rho_0,pzzva)

            leaf_scat=coef*
     &       (cos(theta_i)**2*pyyhd*conjg(pxxva)
     &       +sin(theta_i)**2*pyyhd*conjg(pzzva))

          else if (mech_c .eq. gp) then

            l_tensor_flag=0
            call integral(pyyhhd_l,zero,rho_0,pyyhd)

            l_tensor_flag=0
            call integral(pxxvvb_l,zero,rho_0,pxxvb)

            l_tensor_flag=0
            call integral(pzzvvb_l,zero,rho_0,pzzvb)
            
            leaf_scat=coef*
     &       (cos(theta_i)**2*pyyhd*conjg(pxxvb)
     &       +sin(theta_i)**2*pyyhd*conjg(pzzvb))

          else if (mech_c .eq. pg) then

            l_tensor_flag=0
            call integral(pyyhhd_l,zero,rho_0,pyyhd)

            l_tensor_flag=0
            call integral(pxxvvc_l,zero,rho_0,pxxvc)

            l_tensor_flag=0
            call integral(pzzvvc_l,zero,rho_0,pzzvc)
            
            leaf_scat=coef*
     &       (cos(theta_i)**2*pyyhd*conjg(pxxvc)
     &       +sin(theta_i)**2*pyyhd*conjg(pzzvc))

          else if (mech_c .eq. gpg) then

            l_tensor_flag=0
            call integral(pyyhhd_l,zero,rho_0,pyyhd)

            l_tensor_flag=0
            call integral(pxxvvd_l,zero,rho_0,pxxvd)

            l_tensor_flag=0
            call integral(pzzvvd_l,zero,rho_0,pzzvd)
            
            leaf_scat=coef*
     &       (cos(theta_i)**2*pyyhd*conjg(pxxvd)
     &       +sin(theta_i)**2*pyyhd*conjg(pzzvd))

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

          endif

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

        endif

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

      endif

      return
      end
      

c************************************************************************
       
      complex function pxxvva_l(rho)

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

       complex ii,phia
       complex p11,p22,p33
       complex mv_leaf,mh_leaf

       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh

       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /leaf_prop/ mv_leaf, mh_leaf

       external asinh

       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)
       
       x=rho*cosa

       z=c1*rho+c2*rho**2

c       phia=cexp(ii*2.0*k_0*(sin(theta_i)*x+cos(theta_i)*z)
c     & + 2.0*mv_leaf*z/cos(theta_i))
       
       phia=cexp(ii*2.0*k_0*(sin(theta_i)*x+cos(theta_i)*z))
       
       dz=c1+2.0*c2*rho

       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl
       
      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif


       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pxxvva_l=phia*(cosa**2*cosb**2*p11
     & +sina**2*p22
     & +cosa**2*sinb**2*p33)*dl

c       print*,'p11,p22,p33,cosa,cosb,sina,sinb,dl'
c       print*,p11,p22,p33,cosa,cosb,sina,sinb,dl
       
       return

       end
      

c************************************************************************
      
      complex function pxxvvb_l(rho)

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

       complex ii,phib
       complex p11,p22,p33

       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0

       external asinh
      
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa

       z=c1*rho+c2*rho**2

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

       dz=c1+2.0*c2*rho

       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl

      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pxxvvb_l=phib*(cosa**2*cosb**2*p11
     &  +sina**2*p22
     &  +cosa**2*sinb**2*p33)*dl

       return

       end       

c************************************************************************
      
      complex function pxxvvc_l(rho)

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

       complex ii,phic
       complex p11,p22,p33

       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0

       external asinh
      
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa

       z=c1*rho+c2*rho**2

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

       dz=c1+2.0*c2*rho

       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl

      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pxxvvc_l=phic*(cosa**2*cosb**2*p11
     &  +sina**2*p22
     &  +cosa**2*sinb**2*p33)*dl

       return

       end              

c************************************************************************
      
      complex function pxxvvd_l(rho)

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

       complex ii,phid
       complex p11,p22,p33
       complex mv_leaf,mh_leaf
       
       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /leaf_prop/ mv_leaf, mh_leaf

       external asinh
       
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa

       z=c1*rho+c2*rho**2

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

       dz=c1+2.0*c2*rho
       
       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl
       
      if (c2 .eq. 0.0) then

         ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pxxvvd_l=phid*(cosa**2*cosb**2*p11
     &  +sina**2*p22
     &  +cosa**2*sinb**2*p33)*dl

       return

       end              

c************************************************************************
      
      complex function pzzvva_l(rho)

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

       complex ii,phia
       complex p11,p22,p33
       complex mv_leaf,mh_leaf
       
       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /leaf_prop/ mv_leaf, mh_leaf

       external asinh
       
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa
       
       z=c1*rho+c2*rho**2

c       phia=cexp(ii*2.0*k_0*(sin(theta_i)*x+cos(theta_i)*z)
c     & + 2.0*mv_leaf*z/cos(theta_i))
       
       phia=cexp(ii*2.0*k_0*(sin(theta_i)*x+cos(theta_i)*z))
       
       dz=c1+2.0*c2*rho

       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl
       

      if (c2 .eq. 0.0) then

         ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pzzvva_l=phia*(sinb**2*p11
     &  +cosb**2*p33)*dl
       
       return

       end
      

c************************************************************************
      
      complex function pzzvvb_l(rho)

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

       complex ii,phib
       complex p11,p22,p33

       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0

       external asinh
      
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       
       x=rho*cosa

       z=c1*rho+c2*rho**2

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

       dz=c1+2.0*c2*rho
       dl=sqrt(1.0+dz**2)

       cosb=dz/dl
       sinb=1.0/dl

      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pzzvvb_l=phib*(sinb**2*p11
     &  +cosb**2*p33)*dl

       return

       end       

************************************************************************
      
      complex function pzzvvc_l(rho)

************************************************************************

       complex ii,phic
       complex p11,p22,p33

       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0

       external asinh
      
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa

       z=c1*rho+c2*rho**2

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

       dz=c1+2.0*c2*rho

       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl       

      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pzzvvc_l=phic*(sinb**2*p11
     &  +cosb**2*p33)*dl

       return

       end              

************************************************************************
      
      complex function pzzvvd_l(rho)

************************************************************************

       complex ii,phid
       complex p11,p22,p33
       complex mv_leaf,mh_leaf
       
       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /leaf_prop/ mv_leaf, mh_leaf

       external asinh
       
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa

       z=c1*rho+c2*rho**2

c       phid=cexp(ii*2.0*k_0*(sin(theta_i)*x-cos(theta_i)*z)
c     &  -2.0*mv_leaf*z/cos(theta_i))
       
       phid=cexp(ii*2.0*k_0*(sin(theta_i)*x-cos(theta_i)*z))
       
       dz=c1+2.0*c2*rho

       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl       

      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pzzvvd_l=phid*(sinb**2*p11
     &  +cosb**2*p33)*dl

       return

       end              

************************************************************************
      
      complex function pxzvva_l(rho)

************************************************************************

       complex ii,phia
       complex p11,p22,p33
       complex mv_leaf,mh_leaf
       
       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /leaf_prop/ mv_leaf, mh_leaf

       external asinh
       
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa

       z=c1*rho+c2*rho**2

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

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

       dz=c1+2.0*c2*rho

       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl
       
      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pxzvva_l=phia*(cosa*cosb*sinb*(p33-p11))
       
       return
       end
      

************************************************************************
      
      complex function pxzvvb_l(rho)
 
************************************************************************

       complex ii,phib
       complex p11,p22,p33

       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0

       external asinh
      
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa

       z=c1*rho+c2*rho**2

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

       dz=c1+2.0*c2*rho
       
       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl

      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pxzvvb_l=phib*(cosa*cosb*sinb*(p33-p11))

       return

       end       

************************************************************************
      
      complex function pxzvvc_l(rho)

************************************************************************

       complex ii,phic
       complex p11,p22,p33

       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0

       external asinh
      
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa

       z=c1*rho+c2*rho**2

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

       dz=c1+2.0*c2*rho

       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl       

      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pxzvvc_l=phic*(cosa*cosb*sinb*(p33-p11))

       return

       end              

************************************************************************
      
      complex function pxzvvd_l(rho)

************************************************************************

       complex ii,phid
       complex p11,p22,p33
       complex mv_leaf,mh_leaf       

       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /leaf_prop/ mv_leaf, mh_leaf

       external asinh
       
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa

       z=c1*rho+c2*rho**2

c       phid=cexp(ii*2.0*k_0*(sin(theta_i)*x-cos(theta_i)*z)
c     &  -2.0*mv_leaf*z/cos(theta_i))
       
       phid=cexp(ii*2.0*k_0*(sin(theta_i)*x-cos(theta_i)*z))
       
       dz=c1+2.0*c2*rho

       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl       

      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pxzvvd_l=phid*(cosa*cosb*sinb*(p33-p11))

       return

       end              

************************************************************************
      
      complex function pyyhha_l(rho)

************************************************************************

       complex ii,phia
       complex p11,p22,p33
       complex mv_leaf,mh_leaf       

       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /leaf_prop/ mv_leaf, mh_leaf

       external asinh
       
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa       
       z=c1*rho+c2*rho**2

c       phia=cexp(ii*2.0*k_0*(sin(theta_i)*x+cos(theta_i)*z)
c     & + 2.0*mh_leaf*z/cos(theta_i))
       
       phia=cexp(ii*2.0*k_0*(sin(theta_i)*x+cos(theta_i)*z))
       
       dz=c1+2.0*c2*rho
       dl=sqrt(1.0+dz**2)

       cosb=dz/dl
       sinb=1.0/dl
       
      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pyyhha_l=phia*(sina**2*cosb**2*p11
     &  +cosa**2*p22
     &  +sina**2*sinb**2*p33)*dl

       return

       end
      

************************************************************************
      
      complex function pyyhhb_l(rho)

************************************************************************

       complex ii,phib
       complex p11,p22,p33

       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0

       external asinh
      
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa

       z=c1*rho+c2*rho**2

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

       dz=c1+2.0*c2*rho

       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl

      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pyyhhb_l=phib*(sina**2*cosb**2*p11
     &  +cosa**2*p22
     &  +sina**2*sinb**2*p33)*dl

C       print*,'alpha,c1,p11,p22,p33,pyyhhb_l'
C       print*, alpha,c1,p11,p22,p33,pyyhhb_l
       
       return

       end       

************************************************************************
      
      complex function pyyhhc_l(rho)

************************************************************************

       complex ii,phic
       complex p11,p22,p33

       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0
       
       external asinh
       
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa

       z=c1*rho+c2*rho**2

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

       dz=c1+2.0*c2*rho

       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl

      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pyyhhc_l=phic*(sina**2*cosb**2*p11
     &  +cosa**2*p22
     &  +sina**2*sinb**2*p33)*dl

       return

       end              

c************************************************************************
      
      complex function pyyhhd_l(rho)

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

       complex ii,phid
       complex p11,p22,p33
       complex mv_leaf,mh_leaf       

       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /leaf_prop/ mv_leaf, mh_leaf

       external asinh
       
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa

       z=c1*rho+c2*rho**2

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

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

       dz=c1+2.0*c2*rho
       
       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl
       
      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pyyhhd_l=phid*(sina**2*cosb**2*p11
     &  +cosa**2*p22
     &  +sina**2*sinb**2*p33)*dl

       return

       end              

c************************************************************************
      
      complex function pxyhva_l(rho)

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

       complex ii,phia
       complex p11,p22,p33
       complex mv_leaf,mh_leaf       

       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /leaf_prop/ mv_leaf, mh_leaf

       external asinh       
      
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa

       z=c1*rho+c2*rho**2

c       phia=cexp(ii*2.0*k_0*(sin(theta_i)*x+cos(theta_i)*z)
c     & + (mv_leaf+mh_leaf)*z/cos(theta_i))
       
       phia=cexp(ii*2.0*k_0*(sin(theta_i)*x+cos(theta_i)*z))
       
       dz=c1+2.0*c2*rho

       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl

      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pxyhva_l=phia*(sina*cosa*cosb**2*p11
     &  -cosa*sina*p22
     &  +sina*cosa*sinb**2*p33)*dl

       return

       end
      

c************************************************************************
      
      complex function pxyhvb_l(rho)

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

       complex ii,phib
       complex p11,p22,p33
       complex mv_leaf,mh_leaf
       
       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /leaf_prop/ mv_leaf, mh_leaf

       external asinh       
      
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa

       z=c1*rho+c2*rho**2

c       phib=cexp(ii*2.0*k_0*(sin(theta_i)*x)
c     &  +(mh_leaf-mv_leaf)*z/cos(theta_i))
       
       phib=cexp(ii*2.0*k_0*(sin(theta_i)*x))
       
       dz=c1+2.0*c2*rho

       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl

      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pxyhvb_l=phib*(sina*cosa*cosb**2*p11
     &  -cosa*sina*p22
     &  +sina*cosa*sinb**2*p33)*dl

       return

       end       

c************************************************************************
      
      complex function pxyhvc_l(rho)

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

       complex ii,phic
       complex p11,p22,p33
       complex mv_leaf,mh_leaf       

       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /leaf_prop/ mv_leaf, mh_leaf

       external asinh
      
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa

       z=c1*rho+c2*rho**2

c       phic=cexp(ii*2.0*k_0*(sin(theta_i)*x)
c     &  +(mv_leaf-mh_leaf)*z/cos(theta_i))
       
       phic=cexp(ii*2.0*k_0*(sin(theta_i)*x))
       
       dz=c1+2.0*c2*rho

       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl

      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pxyhvc_l=phic*(sina*cosa*cosb**2*p11
     &  -cosa*sina*p22
     &  +sina*cosa*sinb**2*p33)*dl
       
       return

       end              

c************************************************************************
      
      complex function pxyhvd_l(rho)

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

       complex ii,phid
       complex p11,p22,p33
       complex mv_leaf,mh_leaf
       
       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /leaf_prop/ mv_leaf, mh_leaf

       external asinh
      
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa

       z=c1*rho+c2*rho**2

c       phid=cexp(ii*2.0*k_0*(sin(theta_i)*x-cos(theta_i)*z)
c     &  -(mv_leaf+mh_leaf)*z/cos(theta_i))
       
       phid=cexp(ii*2.0*k_0*(sin(theta_i)*x-cos(theta_i)*z))

       dz=c1+2.0*c2*rho
       
       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl

      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pxyhvd_l=phid*(sina*cosa*cosb**2*p11
     &  -cosa*sina*p22
     &  +sina*cosa*sinb**2*p33)*dl

       return

       end

      
c************************************************************************
      
      complex function pyzhva_l(rho)

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

       complex ii,phia
       complex p11,p22,p33
       complex mv_leaf,mh_leaf
       
       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /leaf_prop/ mv_leaf, mh_leaf

       external asinh
       
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa

       z=c1*rho+c2*rho**2

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

       dz=c1+2.0*c2*rho

       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl

      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pyzhva_l=phia*(sina*sinb*cosb*(p33-p11))*dl

       return

       end
      

c************************************************************************
      
      complex function pyzhvb_l(rho)

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

       complex ii,phib
       complex p11,p22,p33
       complex mv_leaf,mh_leaf
       
       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /leaf_prop/ mv_leaf, mh_leaf

       external asinh       
      
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa

       z=c1*rho+c2*rho**2

c       phib=cexp(ii*2.0*k_0*(sin(theta_i)*x)
c     &  +(mh_leaf-mv_leaf)*z/cos(theta_i))
       
       phib=cexp(ii*2.0*k_0*(sin(theta_i)*x))
       
       dz=c1+2.0*c2*rho

       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl

      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pyzhvb_l=phib*(sina*sinb*cosb*(p33-p11))*dl

       return

       end       

c************************************************************************
      
      complex function pyzhvc_l(rho)

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

       complex ii,phic
       complex p11,p22,p33
       complex mv_leaf,mh_leaf

       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /leaf_prop/ mv_leaf, mh_leaf

       external asinh       
      
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa

       z=c1*rho+c2*rho**2

c       phic=cexp(ii*2.0*k_0*(sin(theta_i)*x)
c     &  +(mv_leaf-mh_leaf)*z/cos(theta_i))
       
       phic=cexp(ii*2.0*k_0*(sin(theta_i)*x))

       dz=c1+2.0*c2*rho

       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl

      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pyzhvc_l=phic*(sina*sinb*cosb*(p33-p11))*dl

       return

       end              

c************************************************************************
      
      complex function pyzhvd_l(rho)

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

       complex ii,phid
       complex p11,p22,p33
       complex mv_leaf,mh_leaf
       
       real c1,c2,alpha
       real rho,x,z
       real theta_i,phi_i,k_0
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0,asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       common /incident/ theta_i,phi_i,k_0
       common /leaf_prop/ mv_leaf, mh_leaf
       
       external asinh
      
       ii=cmplx(0.0,1.0)

       cosa=cos(alpha)
       sina=sin(alpha)

       x=rho*cosa

       z=c1*rho+c2*rho**2

c       phid=cexp(ii*2.0*k_0*(sin(theta_i)*x-cos(theta_i)*z)
c    &  -(mv_leaf+mh_leaf)*z/cos(theta_i))
       
       phid=cexp(ii*2.0*k_0*(sin(theta_i)*x-cos(theta_i)*z))
       
       dz=c1+2.0*c2*rho
       
       dl=sqrt(1.0+dz**2)

       cosb=dz/dl

       sinb=1.0/dl

      if (c2 .eq. 0.0) then

        ell=dl*rho

      else

        ell=((dz*sqrt(1+dz**2))-c1*sqrt(1+c1**2)-asinh(c1)+asinh(dz))
     &      /(4.0*c2)

      endif

       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pyzhvd_l=phid*(sina*sinb*cosb*(p33-p11))*dl

       return

       end              
