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

      complex function leaf_forward_scat(phi,v,pol)

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

      integer vv,hh,hv
      parameter(vv=1)
      parameter(hh=2)
      parameter(hv=3)

      complex pxxfsv_l,pzzfsv_l,pxzfsv_l
      complex pyyfsh_l
      complex pxx,pyy,pzz,pxz
      
      real c1,c2,alpha,ell_0
      real theta_i,phi_i,k_0
      real zero,rho_0,v,phi
      real pi,coef
      real asinh
      
      integer pol

      external pxxfsv_l,pzzfsv_l,pxzfsv_l
      external pyyfsh_l
      external asinh
      
      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**2/(4.0*pi)

      call leaf_model(v, c1,c2,rho_0,ell_0)
      
      alpha=phi
      
      if (pol .eq. vv) then

        l_tensor_flag=0
        call integral(pxxfsv_l,zero,rho_0,pxx)
          
        l_tensor_flag=0
        call integral(pzzfsv_l,zero,rho_0,pzz)

        l_tensor_flag=0
        call integral(pxzfsv_l,zero,rho_0,pxz)
          
c     write(6,*) 'phi = ',alpha
c     print*,'pxx = ',pxx
c     print*,'pzz = ',pzz
c     print*,'coef = ',coef

        leaf_forward_scat=coef*
     &       (cos(theta_i)**2*pxx 
     &       -2.0*cos(theta_i)*sin(theta_i)*pxz
     &       +sin(theta_i)**2*pzz)


      else if (pol .eq. hh) then
        
        l_tensor_flag=0
        call integral(pyyfsh_l,zero,rho_0,pyy)

c     write(6,*)'pyy = ',pyy

        leaf_forward_scat=coef*pyy

c        print*,'leaf_fs',leaf_forward_scat,coef,rho_0,pyy

      else if (pol .eq. hv) then

        leaf_forward_scat=cmplx(0.0,0.0)

        write(error,*)'pol = hv in func leaf_fs was encountered',
     &       ' but not expected'
        code=23
        call the_death_code(warn,error,code)

      else
          
        write(error,*)'pol = ',pol,' in func leaf_fs is not a',
     &       ' valid integer (1,2,or 3)'
        code=24
        call the_death_code(die,error,code)

      endif

      return
      end      


c************************************************************************
      
       complex function pxxfsv_l(rho)

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

       complex p11,p22,p33

       real c1,c2,alpha
       real rho
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0
       real asinh
      
       external asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       

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

       cosa=cos(alpha)
       sina=sin(alpha)
       
       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
c        print*, ell
       call leaf_tensor(ell,ell_0, p11,p22,p33)

       pxxfsv_l=(cosa**2*cosb**2*p11
     & + sina**2*p22 + cosa**2*sinb**2*p33)*dl

c      print*,'p11fsv_l'
c       print*,'p11 = ',p11
c       print*,'p22 = ',p22
c       print*,'p33 = ',p33

       return
       end
      
c************************************************************************
      
       complex function pzzfsv_l(rho)

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

       complex p11,p22,p33
       
       real c1,c2,alpha
       real rho
       real dz,dl,cosb,sinb
       real ell,ell_0
       real asinh
       
       external asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0

       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
c        print*, ell
       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pzzfsv_l=(sinb**2*p11+cosb**2*p33)*dl
       
C       print*,'p33fsv_l'
C       print*,'p11 = ',p11
C       print*,'p22 = ',p22
C       print*,'p33 = ',p33

       return
       end

      
c************************************************************************
      
       complex function pxzfsv_l(rho)

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

       complex p11,p22,p33
       
       real c1,c2,alpha
       real rho
       real dz,dl,cosb,sinb,cosa,sina
       real ell,ell_0
       real asinh
       
       external asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0

       dz=c1+2.0*c2*rho
       dl=sqrt(1.0+dz**2)
       
       cosa=cos(alpha)
       sina=sin(alpha)
       
       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
c        print*, ell
       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pxzfsv_l=(cosa*cosb*sinb*(p33-p11))*dl
       
C       print*,'p33fsv_l'
C       print*,'p11 = ',p11
C       print*,'p22 = ',p22
C       print*,'p33 = ',p33

       return
       end      

c************************************************************************
      
       complex function pyyfsh_l(rho)

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

       complex p11,p22,p33

       real c1,c2,alpha
       real rho
       real dz,dl,cosb,sinb,cosa,sina
       real asinh
       real ell,ell_0
       
       external asinh
       
       common /leaf_struct/ c1,c2,alpha,ell_0
       
       dz=c1+2.0*c2*rho
       dl=sqrt(1.0+dz**2)
       
       cosa=cos(alpha)
       sina=sin(alpha)
       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
c        print*, ell
       call leaf_tensor(ell,ell_0, p11,p22,p33)
       
       pyyfsh_l=(sina**2*cosb**2*p11+cosa**2*p22
     &      +sina**2*sinb**2*p33)*dl
       
C       print*,'p22fsv_l'
C       print*,'p11 = ',p11
C       print*,'p22 = ',p22
C       print*,'p33 = ',p33
       
       return
       end