


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

      
      subroutine leaf_leaf(ss_ll)


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


      complex ss_sl(4,4,4)
      complex s(4,4,4),t(4,4)
      complex ss,tt
      
      
      
      do 10 i=1,4

        do 20 j=1,4

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

              if (k .eq. 1) then
 
                call lphi_lphi_int(i,j,k, ss)

                call leaf_field_int(i,k, tt)

                s(i,j,k)=ss

                t(i,k)=tt

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

              else if(k .eq. 2) then

                call lphi_lphi_int(i,j,k, ss)

                call leaf_field_int(i,k, tt)

                s(i,j,k)=ss

                t(i,k)=tt

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

              else if(k .eq. 3) then

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

                
              else if(k .eq. 4) then

                call lphi_lphi_int(i,j,k, ss)

                call leaf_field_int(i,k, tt)

                s(i,j,k)=ss

                t(i,k)=tt

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

              else

                print*, 'error in stalk_leaf a'

              endif

            else if(j .eq. 2) then

              if (k .eq. 1) then

                ss_ll(i,j,k)=conjg(ss_ll(i,k,j))
                
              else if(k .eq. 2) then
                
                call lphi_lphi_int(i,j,k, ss)

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

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


              else if(k .eq. 4) then

                call lphi_lphi_int(i,j,k, ss)

                s(i,j,k)=ss

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

              else

                print*, 'error in stalk_leaf b'

              endif

            else if(j .eq. 3) then

              if(k .eq. 1) then

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

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

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

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

                print*, 'error in  stalk_leaf c'

              endif

            else if (j .eq. 4) then

              if (k .eq. 1) then

                call lphi_lphi_int(i,j,k, ss)

                s(i,j,k)=ss

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

              else if (k .eq. 2) then

                ss_ll(i,j,k)=conjg(ss_ll(i,k,j))

              else if (k .eq. 3) then

                ss_ll(i,j,k)=ss_ll(i,j,2)
               
              else if (k .eq. 4) then
                
                call lphi_lphi_int(i,j,k, ss)

                s(i,j,k)=ss
                
                ss_ll(i,j,k)=s(i,j,k)*t(i,j)*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 lphi_lphi_int(i,j,k, ss)

                call leaf_field_int(i,k, tt)

                s(i,j,k)=ss

                t(i,k)=tt

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

              else if(k .eq. 2) then

                call lphi_lphi_int(i,j,k, ss)

                call leaf_field_int(i,k, tt)

                s(i,j,k)=ss

                t(i,k)=tt

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

              else if(k .eq. 3) then

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

                
              else if(k .eq. 4) then

                call lphi_lphi_int(i,j,k, ss)

                call leaf_field_int(i,k, tt)

                s(i,j,k)=ss

                t(i,k)=tt

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

              else

                print*, 'error in stalk_leaf a'

              endif

            else if(j .eq. 2) then

              if (k .eq. 1) then

                ss_ll(i,j,k)=conjg(ss_ll(i,k,j))
                
              else if(k .eq. 2) then
                
                call lphi_lphi_int(i,j,k, ss)

                s(i,j,k)=ss

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

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


              else if(k .eq. 4) then

                call lphi_lphi_int(i,j,k, ss)

                s(i,j,k)=ss

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

              else

                print*, 'error in stalk_leaf b'

              endif

            else if(j .eq. 3) then

              if(k .eq. 1) then

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

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

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

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

                print*, 'error in  stalk_leaf c'

              endif

            else if (j .eq. 4) then

              if (k .eq. 1) then

                call lphi_lphi_int(i,j,k, ss)

                s(i,j,k)=ss

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

              else if (k .eq. 2) then

                ss_ll(i,j,k)=conjg(ss_ll(i,k,j))

              else if (k .eq. 3) then

                ss_ll(i,j,k)=ss_ll(i,j,2)
               
              else if (k .eq. 4) then
                
                call lphi_lphi_int(i,j,k, ss)

                s(i,j,k)=ss

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

              else

                print*, 'error in  stalk_leaf d'
                
              endif

            else
              
              print*, 'error in  stalk_leaf e'
                
            endif

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

            endif
            

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

              if (k .eq. 1) then
 
                call lphi_lphi_int(i,j,k, ss)

                s(i,j,k)=ss

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

              else if(k .eq. 2) then

                call lphi_lphi_int(i,j,k, ss)

                s(i,j,k)=ss

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

              else if(k .eq. 3) then

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

              else if(k .eq. 4) then

                call lphi_lphi_int(i,j,k, ss)

                s(i,j,k)=ss

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

                print*, 'error in stalk_leaf p'

              endif

            else if(j .eq. 2) then

              if (k .eq. 1) then

                call lphi_lphi_int(i,j,k, ss)

                s(i,j,k)=ss

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


              else if(k .eq. 2) then
                
                call lphi_lphi_int(i,j,k, ss)

                s(i,j,k)=ss

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

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

              else if(k .eq. 4) then
                
                call lphi_lphi_int(i,j,k, ss)

                s(i,j,k)=ss

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

              else

                print*, 'error in stalk_leaf q'

              endif

            else if(j .eq. 3) then

              if(k .eq. 1) then

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

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

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

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

                print*, 'error in  stalk_leaf r'

              endif

            else if (j .eq. 4) then

              if (k .eq. 1) then

                call lphi_lphi_int(i,j,k, ss)

                s(i,j,k)=ss

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

              else if (k .eq. 2) then

                call lphi_lphi_int(i,j,k, ss)

                s(i,j,k)=ss

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

              else if (k .eq. 3) then

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

              else if (k .eq. 4) then
                
                call lphi_lphi_int(i,j,k, ss)

                s(i,j,k)=ss

                ss_ll(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
            
              print*, 'error in  stalk_leaf y'

          endif

          
   30     continue  

   20   continue

   10 continue

      return

      end



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

      subroutine lphi_lphi_int(i,j,k, p)

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

      complex h_ll,p
      real pi,mpi
      integer i,j,k

      common /index/ i,j,k

      pi=4.0*atan(1.0)
      mpi=-pi
      

      external h_ll
      call qgausx(h_ll,mpi,pi,p)

      return

      end
      

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

      function h_ll(phii)

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

      complex h_ll,s,g_ll

      real amin,amax,z0min,z0max
      real phi_s,a_s,z0
      real phii
      
      external g_ll
      
      common /stalk_limits/ amin,amax,z0min,z0max
      common /stalk_state/ phi_s,a_s,z0

      phi_s=phii

      call qgausy(g_ll,amin,amax,s)

      h_ll=s

      return

      end

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

      function g_ll(aa)

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


      complex g_ll,i_ll,s

      real amin,amax,z0min,z0max
      real phi_s,a_s,z0
      real aa

      external i_ll
      
      common /stalk_limits/ amin,amax,z0min,z0max
      common /stalk_state/ phi_s,a_s,z0

      a_s=aa

      call qgausy(i_ll,z0min,z0max,s)

      g_ll=s

      return

      end

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

      function i_ll(zz0)

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


      complex i_ll,lphi_lphi_ave

      real phi_s,a_s,z0
      real zz0
      
      common /stalk_state/ phi_s,a_s,z0
      common /index/ i,j,k
      
      z0=zz0
      
      i_ll=lphi_lphi_ave(phi_s,a_s,z0,i,j,k)

      return

      end

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

      function lphi_lphi_ave(alpha_s,a,z0,i,j,k)

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


      complex lphi_lphi_ave
      complex kxz,kzz
      complex k1_c,k1,k2_c,k2
      complex ii
      complex tau_v,tau_h
      complex mv_leaf,mh_leaf
      complex r
      
      
      real dm,z0
      real a,alpha_s
      real theta_i,phi_i,k_0
      real zero
      real stalk_alpha_pdf,stalk_a_pdf,stalk_z0_pdf
      real dm,del_leaf,z_lc,hole
      
      integer i,j,k

      common /incident/ theta_i,phi_i,k_0
      common /layer_struct/ dm,del_leaf,z_lc,hole
      common /leaf_prop/  mv_leaf,mh_leaf

      zero=0.0
      ii=cmplx(0.0,1.0)


      kxz=sin(theta_i)*cos(alpha)*a

      kzz=cos(theta_i)


      if (i .eq. 1) then

        if (j .eq. 1) then

          if (k .eq. 1) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &      +2.0*tau_v(z_lc)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz+kzz)+2.0*mv_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &      +2.0*tau_v(z_lc)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz+kzz)+2.0*mv_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)*
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if(k .eq. 2) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +2.0*tau_v(z_lc)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz+kzz)+2.0*mv_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &      +2.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_v(theta_i))
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 3) then
     
            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
            +2.0*tau_v(z_lc)/cos(theta_i)
     
            k1=ii*2.0*k_0*(kxz+kzz)+2.0*mv_leaf/cos(theta_i)
     
            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_v(dm)/cos(theta_i))
     
            k2=conjg(ii*2.0*k_0*kxz)
     
            r=cexp(k1_c)*cexp(k2_c)*conjg(r_v(theta_i))
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       (-cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2)

          else if (k .eq. 4) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +2.0*tau_v(z_lc)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz+kzz)+2.0*mv_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -2.0*tau_v(z_lc)/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm+4.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz-kzz)-2.0*mv_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)
     &       *conjg(r_v(theta_i))*conjg(r_v(theta_i))*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
            /(k1**2+k1*k2))

          else

            print*, 'error in leaf_phase_int a'

          endif

        else if (j .eq. 2) then


          if (k .eq. 1) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_v(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +2.0*tau_v(z_lc)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz+kzz)+2.0*mv_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*r_v(theta_i)
     &       *(-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 2) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_v(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*r_v(theta_i)*conjg(r_v(theta_i))*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 3) then

            k1_c=ii*2.0*k_0*(kxz)*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_v(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*(kxz)*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*r_v(theta_i)*conjg(r_v(theta_i))*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       (-cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2)

          else if (k .eq. 4) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_v(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -2.0*tau_v(z_lc)/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm+4.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz-kzz)-2.0*mv_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*r_v(theta_i)*
     &       *conjg(r_v(theta_i))*conjg(r_v(theta_i))*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))

          else

            print*, 'error in leaf_phase_int b'

          endif

          
        else if (j .eq. 3) then

          if (k .eq. 1) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
            +2.0*tau_v(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
            +2.0*tau_v(z_lc)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz+kzz)+2.0*mv_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*r_v(theta_i)
     &       *(-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if(k .eq. 2) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_v(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*r_v(theta_i)*conjg(r_v(theta_i))
     &       *(-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 3) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_v(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*r_v(theta_i)*conjg(r_v(theta_i))
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       (-cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2)

          else if (k .eq. 4) then

            k1_c=ii*2.0*k_0*(kxz)*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_v(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -2.0*tau_v(z_lc)/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm+4.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz-kzz)-2.0*mv_leaf/cos(theta_i))

             r=cexp(k1_c)*cexp(k2_c)*r_v(theta_i)
     &       *conjg(r_v(theta_i))*conjg(r_v(theta_i))*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))

          else

            print*, 'error in leaf_phase_int c'

          endif


        else if (j .eq. 4) then

          if (k .eq. 1) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -2.0*tau_v(z_lc)/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm+4.0*tau_v(dm)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz-kzz)-2.0*mv_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +2.0*tau_v(z_lc)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz+kzz)+2.0*mv_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*r_v(theta_i)**2*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if(k .eq. 2) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -2.0*tau_v(z_lc)/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm+4.0*tau_v(dm)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz-kzz)-2.0*mv_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_v(theta_i))
     &       *r_v(theta_i)**2*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 3) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -2.0*tau_v(z_lc)/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm+4.0*tau_v(dm)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz-kzz)-2.0*mv_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_v(theta_i))*
     &       *r_v(theta_i)**2*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       (-cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2)

          else if (k .eq. 4) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -2.0*tau_v(z_lc)/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm+4.0*tau_v(dm)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz-kzz)-2.0*mv_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -2.0*tau_v(z_lc)/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm+4.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz-kzz)-2.0*mv_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*r_v(theta_i)**2
     &       *conjg(r_v(theta_i))*conjg(r_v(theta_i))*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))

          else

            print*, 'error in leaf_phase_int d'

          endif

        else

          print*, 'error in leaf_phase_int e'

        endif

      else if (i .eq. 2) then

        if (j .eq. 1) then

          if (k .eq. 1) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +2.0*tau_h(z_lc)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz+kzz)+2.0*mh_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +2.0*tau_h(z_lc)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz+kzz)+2.0*mh_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if(k .eq. 2) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +2.0*tau_h(z_lc)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz+kzz)+2.0*mh_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_h(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_h(theta_i))*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 3) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +2.0*tau_h(z_lc)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz+kzz)+2.0*mh_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_h(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_h(theta_i))*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       (-cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2)

          else if (k .eq. 4) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +2.0*tau_h(z_lc)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz+kzz)+2.0*mh_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -2.0*tau_h(z_lc)/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm+4.0*tau_h(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz-kzz)-2.0*mh_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)
     &       *conjg(r_h(theta_i))*conjg(r_h(theta_i))*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))

          else

            print*, 'error in leaf_phase_int f'

          endif

        else if (j .eq. 2) then


          if (k .eq. 1) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +2.0*tau_h(z_lc)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz+kzz)+2.0*mh_leaf/cos(theta_i))

     &       r=cexp(k1_c)*cexp(k2_c)*r_h(theta_i)
     &       *(-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 2) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_h(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_h(theta_i))*r_h(theta_i)
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 3) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
            +2.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
            +2.0*tau_h(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_h(theta_i))*r_h(theta_i)*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       (-cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2)

          else if (k .eq. 4) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -2.0*tau_h(z_lc)/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm+4.0*tau_h(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz-kzz)-2.0*mh_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*r_h(theta_i)
     &       *conjg(r_h(theta_i))*conjg(r_h(theta_i))*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))

          else

            print*, 'error in leaf_phase_int g'

          endif

          
        else if (j .eq. 3) then

          if (k .eq. 1) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +2.0*tau_h(z_lc)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz+kzz)+2.0*mh_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*r_h(theta_i)
     &        (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &        +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &        /(k1**2+k1*k2)
     &        +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &        /(k2**2+k1*k2)
     &        -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &        +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &        /(k2**2+k2*k1)
     &        +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &        /(k1**2+k1*k2))
             
          else if (k .eq. 2) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_h(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_h(theta_i))*r_h(theta_i)*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 3) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_h(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_h(theta_i))*r_h(theta_i)*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       (-cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2)

          else if (k .eq. 4) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -2.0*tau_h(z_lc)/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm+4.0*tau_h(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz-kzz)-2.0*mh_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*r_h(theta_i)
     &       *conjg(r_h(theta_i))*conjg(r_h(theta_i))*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))

          else

            print*, 'error in leaf_phase_int h'

          endif


        else if (j .eq. 4) then

          if (k .eq. 1) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -2.0*tau_h(z_lc)/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm+4.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz-kzz)-2.0*mh_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +2.0*tau_h(z_lc)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz+kzz)+2.0*mh_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*r_h(theta_i)**2*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 2) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -2.0*tau_h(z_lc)/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm+4.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz-kzz)-2.0*mh_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_h(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_h(theta_i))*
     &       r_h(theta_i)**2*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 3) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -2.0*tau_h(z_lc)/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm+4.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz-kzz)-2.0*mh_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_h(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_h(theta_i))*
     &       r_h(theta_i)**2*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       (-cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2)
     
          else if (k .eq. 4) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -2.0*tau_h(z_lc)/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm+4.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz-kzz)-2.0*mh_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -2.0*tau_h(z_lc)/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm+4.0*tau_h(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz-kzz)-2.0*mh_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*r_h(theta_i)**2
     &       *conjg(r_h(theta_i))*conjg(r_h(theta_i))*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))

          else

            print*, 'error in leaf_phase_int i'

          endif

        else

          print*, 'error in leaf_phase_int j'

        endif


      else if (i .eq. 3) then

        if (j .eq. 1) then

          if (k .eq. 1) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +(tau_h(z_lc)+tau_v(z_lc))/cos(theta_i)

            k1=ii*2.0*k_0*(kxz+kzz)+(mh_leaf+mv_leaf)/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +(tau_h(z_lc)+tau_v(z_lc))/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz+kzz)
     &       +(mh_leaf+mv_leaf)/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))

          else if(k .eq. 2) then


            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +(tau_h(z_lc)+tau_v(z_lc))/cos(theta_i)

            k1=ii*2.0*k_0*(kxz+kzz)+(mh_leaf+mv_leaf)/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)
     &       +(tau_h(z_lc)-tau_v(z_lc)) /cos(theta_i)
     &       +ii*2.0*k_0*kzz*dm+2.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz+(mh_leaf-mv_leaf)/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_v(theta_i))*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 3) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
            +(tau_h(z_lc)+tau_v(z_lc))/cos(theta_i)

            k1=ii*2.0*k_0*(kxz+kzz)+(mh_leaf+mv_leaf)/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)
     &       +(tau_v(z_lc)-tau_h(z_lc))/cos(theta_i)
     &       +ii*2.0*k_0*kzz*dm+2.0*tau_h(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz+(mv_leaf-mh_leaf)/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_h(theta_i))
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       (-cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2)

          else if (k .eq. 4) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +(tau_h(z_lc)+tau_v(z_lc))/cos(theta_i)

            k1=ii*2.0*k_0*(kxz+kzz)+(mh_leaf+mv_leaf)/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -(tau_h(z_lc)+tau_v(z_lc))/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm+2.0*(tau_v(dm)+tau_h(dm))/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz-kzz)-(mh_leaf+mv_leaf)/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)
     &       *conjg(r_h(theta_i))*conjg(r_v(theta_i))*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))

          else

            print*, 'error in leaf_phase_int k'

          endif

        else if (j .eq. 2) then

          if (k .eq. 1) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)
     &       +(tau_h(z_lc)-tau_v(z_lc))/cos(theta_i)
     &       +ii*2.0*k_0*kzz*dm+2.0*tau_v(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz+(mh_leaf-mv_leaf)/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +(tau_h(z_lc)+tau_v(z_lc))/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz+kzz)
     &       +(mh_leaf+mv_leaf)/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*r_v(theta_i)
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 2) then


            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)
     &       +(tau_h(z_lc)-tau_v(z_lc))/cos(theta_i)
     &       +ii*2.0*k_0*kzz*dm+2.0*tau_v(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz+(mh_leaf-mv_leaf)/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)
     &       +(tau_h(z_lc)-tau_v(z_lc)) /cos(theta_i)
     &       +ii*2.0*k_0*kzz*dm+2.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz+(mh_leaf-mv_leaf)/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*r_v(theta_i)*conjg(r_v(theta_i))*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 3) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)
     &       +(tau_h(z_lc)-tau_v(z_lc))/cos(theta_i)
     &       +ii*2.0*k_0*kzz*dm+2.0*tau_v(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz+(mh_leaf-mv_leaf)/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)
     &       +(tau_v(z_lc)-tau_h(z_lc))/cos(theta_i)
     &       +ii*2.0*k_0*kzz*dm+2.0*tau_h(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz+(mv_leaf-mh_leaf)/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_h(theta_i))*r_v(theta_i)
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       (-cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2)

          else if (k .eq. 4) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)
     &       +(tau_h(z_lc)-tau_v(z_lc))/cos(theta_i)
     &       +ii*2.0*k_0*kzz*dm+2.0*tau_v(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz+(mh_leaf-mv_leaf)/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -(tau_h(z_lc)+tau_v(z_lc))/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm+2.0*(tau_v(dm)+tau_h(dm))/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz-kzz)-(mh_leaf+mv_leaf)/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)
     &       *conjg(r_h(theta_i))*conjg(r_v(theta_i))*r_v(theta_i)*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))

          else

            print*, 'error in leaf_phase_int l'

          endif

        else if (j .eq. 3) then

          if (k .eq. 1) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)
     &       +(tau_v(z_lc)-tau_h(z_lc))/cos(theta_i)
     &       +ii*2.0*k_0*kzz*dm+2.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz+(mv_leaf-mh_leaf)/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +(tau_h(z_lc)+tau_v(z_lc))/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz+kzz)
     &       +(mh_leaf+mv_leaf)/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*r_h(theta_i)*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 2) then
     

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)
     &       +(tau_v(z_lc)-tau_h(z_lc))/cos(theta_i)
     &       +ii*2.0*k_0*kzz*dm+2.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz+(mv_leaf-mh_leaf)/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)
     &       +(tau_h(z_lc)-tau_v(z_lc))/cos(theta_i)
     &       +ii*2.0*k_0*kzz*dm+2.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz+(mh_leaf-mv_leaf)/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_v(theta_i))*r_h(theta_i)*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 3) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)
     &       +(tau_v(z_lc)-tau_h(z_lc))/cos(theta_i)
     &       +ii*2.0*k_0*kzz*dm+2.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz+(mv_leaf-mh_leaf)/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)
     &       +(tau_v(z_lc)-tau_h(z_lc))/cos(theta_i)
     &       +ii*2.0*k_0*kzz*dm+2.0*tau_h(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz+(mv_leaf-mh_leaf)/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_h(theta_i))*r_h(theta_i)*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       (-cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2)

          else if (k .eq. 4) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)
     &       +(tau_v(z_lc)-tau_h(z_lc))/cos(theta_i)
     &       +ii*2.0*k_0*kzz*dm+2.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz+(mv_leaf-mh_leaf)/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -(tau_h(z_lc)+tau_v(z_lc))/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm+2.0*(tau_v(dm)+tau_h(dm))/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz-kzz)
     &       -(mh_leaf+mv_leaf)/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)
     &       *conjg(r_h(theta_i))*conjg(r_v(theta_i))*r_h(theta_i)*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))

          else

            print*, 'error in leaf_phase_int m'

          endif

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

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -(tau_h(z_lc)+tau_v(z_lc))/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm
     &       +2.0*(tau_v(dm)+tau_h(dm))/cos(theta_i)

            k1=ii*2.0*k_0*(kxz-kzz)-(mh_leaf+mv_leaf)/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +(tau_h(z_lc)+tau_v(z_lc))/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz+kzz)
     &       +(mh_leaf+mv_leaf)/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*r_h(theta_i)*r_v(theta_i)*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 2) then


            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -(tau_h(z_lc)+tau_v(z_lc))/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm
     &       +2.0*(tau_v(dm)+tau_h(dm))/cos(theta_i)

            k1=ii*2.0*k_0*(kxz-kzz)-(mh_leaf+mv_leaf)/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm) 
     &       +(tau_h(z_lc)-tau_v(z_lc))/cos(theta_i)
     &       +ii*2.0*k_0*kzz*dm+2.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz+(mh_leaf-mv_leaf)/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_v(theta_i))
     &       *r_h(theta_i)*r_v(theta_i)*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 3) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -(tau_h(z_lc)+tau_v(z_lc))/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm
     &       +2.0*(tau_v(dm)+tau_h(dm))/cos(theta_i)

            k1=ii*2.0*k_0*(kxz-kzz)-(mh_leaf+mv_leaf)/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)
     &       +(tau_v(z_lc)-tau_h(z_lc))/cos(theta_i)
     &       +ii*2.0*k_0*kzz*dm+2.0*tau_h(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz+(mv_leaf-mh_leaf)/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_h(theta_i))
     &       *r_h(theta_i)*r_v(theta_i)*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       (-cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2)

          else if (k .eq. 4) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -(tau_h(z_lc)+tau_v(z_lc))/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm
     &       +2.0*(tau_v(dm)+tau_h(dm))/cos(theta_i)

            k1=ii*2.0*k_0*(kxz-kzz)-(mh_leaf+mv_leaf)/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -(tau_h(z_lc)+tau_v(z_lc))/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm
     &       +2.0*(tau_v(dm)+tau_h(dm))/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz-kzz)
     &       -(mh_leaf+mv_leaf)/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*r_h(theta_i)*r_v(theta_i)
     &       *conjg(r_h(theta_i))*conjg(r_v(theta_i))*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))

          else

            print*, 'error in leaf_phase_int m'

          endif


          else

            print*, 'error in leaf_phase_int n'

          endif         


      else if (i .eq. 4) then

        if (j .eq. 1) then

          if (k .eq. 1) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +2.0*tau_h(z_lc)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz+kzz)+2.0*mh_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +2.0*tau_v(z_lc)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz+kzz)
     &       +2.0*mv_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 2) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &      +2.0*tau_h(z_lc)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz+kzz)+2.0*mh_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)
     &       +ii*2.0*k_0*kzz*dm+2.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_v(theta_i))*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 3) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +2.0*tau_h(z_lc)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz+kzz)+2.0*mh_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_v(theta_i))*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       (-cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2)

          else if (k .eq. 4) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +2.0*tau_h(z_lc)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz+kzz)+2.0*mh_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -2.0*tau_v(z_lc)/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm+4.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz-kzz)-2.0*mv_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)
     &       *conjg(r_v(theta_i))*conjg(r_v(theta_i))*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))

          else

            print*, 'error in leaf_phase_int o'

          endif

        else if (j .eq. 2) then

          if (k .eq. 1) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +2.0*tau_v(z_lc)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz+kzz)+2.0*mv_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*r_h(theta_i)*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 2) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_v(theta_i))*r_h(theta_i)*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 3) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_v(theta_i))*r_h(theta_i)*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       (-cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2)

          else if (k .eq. 4) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz
            
            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
     &       -2.0*tau_v(z_lc)/cos(theta_i)
     &       +ii*4.0*k_0*kzz*dm+4.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz-kzz)-2.0*mv_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*r_h(theta_i)
     &       *conjg(r_v(theta_i))*conjg(r_v(theta_i))*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))

          else

            print*, 'error in leaf_phase_int p'

          endif

          
        else if (j .eq. 3) then

          if (k .eq. 1) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
     &       +2.0*tau_v(z_lc)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz+kzz)+2.0*mv_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*r_h(theta_i)*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 2) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_v(theta_i))*r_h(theta_i)*
     &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
     &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
     &       /(k1**2+k1*k2)
     &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
     &       /(k2**2+k1*k2)
     &       -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
     &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
     &       /(k2**2+k2*k1)
     &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
     &       /(k1**2+k1*k2))
            
          else if (k .eq. 3) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
     &       +2.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_v(theta_i))*r_h(theta_i)*
      &       (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
      &       +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
      &       /(k1**2+k1*k2)
      &       +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
      &       /(k2**2+k1*k2)
      &       (-cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
      &       +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
      &       /(k2**2+k2*k1)
      &       +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
      &       /(k1**2+k1*k2)

          else if (k .eq. 4) then

            k1_c=ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
      &      +2.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*kxz

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
      &      -2.0*tau_v(z_lc)/cos(theta_i)
      &      +ii*4.0*k_0*kzz*dm+4.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz-kzz)-2.0*mv_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*r_h(theta_i)*
      &      *conjg(r_v(theta_i))*conjg(r_v(theta_i))*
      &      (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
      &      +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
      &      /(k1**2+k1*k2)
      &      +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
      &      /(k2**2+k1*k2)
      &      -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
      &      +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
      &      /(k2**2+k2*k1)
      &      +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
      &      /(k1**2+k1*k2))

          else

            print*, 'error in leaf_phase_int q'

          endif


        else if (j .eq. 4) then

          if (k .eq. 1) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
      &      -2.0*tau_h(z_lc)/cos(theta_i)
      &      +ii*4.0*k_0*kzz*dm+4.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz-kzz)-2.0*mh_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)+kzz*z_lc)
      &      +2.0*tau_v(z_lc)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz+kzz)+2.0*mv_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)*r_h(theta_i)**2*
      &      (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
      &      +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
      &      /(k1**2+k1*k2)
      &      +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
      &      /(k2**2+k1*k2)
      &      -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
      &      +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
      &      /(k2**2+k2*k1)
      &      +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
      &      /(k1**2+k1*k2))
            
          else if (k .eq. 2) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
      &      -2.0*tau_h(z_lc)/cos(theta_i)
      &      +ii*4.0*k_0*kzz*dm+4.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz-kzz)-2.0*mh_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
      &      +2.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_v(theta_i))*
      &      r_h(theta_i)**2*
      &      (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
      &      +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
      &      /(k1**2+k1*k2)
      &      +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
      &      /(k2**2+k1*k2)
      &      -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
      &      +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
      &      /(k2**2+k2*k1)
      &      +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
      &      /(k1**2+k1*k2))
            
          else if (k .eq. 3) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
      &      -2.0*tau_h(z_lc)/cos(theta_i)
      &      +ii*4.0*k_0*kzz*dm+4.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz-kzz)-2.0*mh_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*kxz*(z_lc-dm)+ii*2.0*k_0*kzz*dm
      &      +2.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*kxz)

            r=cexp(k1_c)*cexp(k2_c)*conjg(r_v(theta_i))*
      &      r_h(theta_i)**2*
      &      (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
      &      +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
      &      /(k1**2+k1*k2)
      &      +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
      &      /(k2**2+k1*k2)
      &      (-cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
      &      +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
      &      /(k2**2+k2*k1)
      &      +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
      &      /(k1**2+k1*k2)
      
          else if (k .eq. 4) then

            k1_c=ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
      &      -2.0*tau_h(z_lc)/cos(theta_i)
      &      +ii*4.0*k_0*kzz*dm+4.0*tau_h(dm)/cos(theta_i)

            k1=ii*2.0*k_0*(kxz-kzz)-2.0*mh_leaf/cos(theta_i)

            k2_c=conjg(ii*2.0*k_0*(kxz*(z_lc-dm)-kzz*z_lc)
      &      -2.0*tau_v(z_lc)/cos(theta_i)
      &      +ii*4.0*k_0*kzz*dm+4.0*tau_v(dm)/cos(theta_i))

            k2=conjg(ii*2.0*k_0*(kxz-kzz)-2.0*mv_leaf/cos(theta_i))

            r=cexp(k1_c)*cexp(k2_c)
      &      *conjg(r_v(theta_i))*conjg(r_v(theta_i))*r_h(theta_i)**2*
      &      (-cexp((k1-k2)*del_leaf/2.0)/(k1*k2)
      &      +cexp(-k1*del_leaf/2.0-k2*del_leaf/2.0 + k1*hole)
      &      /(k1**2+k1*k2)
      &      +cexp(k1*del_leaf/2.0+k2*del_leaf/2.0-k2*hole)
      &      /(k2**2+k1*k2)
      &      -cexp((k2-k1)*del_leaf/2.0)/(k2*k1)
      &      +cexp(-k2*del_leaf/2.0-k1*del_leaf/2.0 + k2*hole)
      &      /(k2**2+k2*k1)
      &      +cexp(k2*del_leaf/2.0+k1*del_leaf/2.0-k1*hole)
      &      /(k1**2+k1*k2))

          else

            print*, 'error in leaf_phase_int r'

          endif

        else

          print*, 'error in leaf_phase_int s'

        endif                                  

      else                                    
                                                 
        print*, 'error in leaf_phase_int t'
                                          
      endif                                       

      return

      end

      


      
 
