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

      complex function stalk_phig_scat(phi,aa,z0,pol,mech,mech_c)

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

      integer pol,mech,mech_c
      
      integer hhhh,vvvv,hvhv,hhvv,vvhh
      parameter(vvvv=1)
      parameter(hhhh=2)
      parameter(hvhv=3)
      parameter(hhvv=4)
      parameter(vvhh=5)
      
      integer p,gp,pg,gpg
      parameter(p=1)
      parameter(gp=2)
      parameter(pg=3)
      parameter(gpg=4)
      
      complex pxxvva_s,pzzvva_s,pxzvva_s
      complex pxxvvb_s,pzzvvb_s,pxzvvb_s
      complex pxxvvc_s,pzzvvc_s,pxzvvc_s
      complex pxxvvd_s,pzzvvd_s,pxzvvd_s

      complex pyyhha_s
      complex pyyhhb_s
      complex pyyhhc_s
      complex pyyhhd_s

      complex pxyhva_s,pyzhva_s
      complex pxyhvb_s,pyzhvb_s
      complex pxyhvc_s,pyzhvc_s
      complex pxyhvd_s,pyzhvd_s

      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 pxxa_mn,pxya_mn,pxza_mn,pyya_mn,pyza_mn,pzza_mn
      complex pxxb_mn,pxyb_mn,pxzb_mn,pyyb_mn,pyzb_mn,pzzb_mn
      complex pxxc_mn,pxyc_mn,pxzc_mn,pyyc_mn,pyzc_mn,pzzc_mn
      complex pxxd_mn,pxyd_mn,pxzd_mn,pyyd_mn,pyzd_mn,pzzd_mn
      
      complex phi_g
      
      real z0
      real a,alpha_s,ell_0
      real theta_i,phi_i,k_0
      real phi,aa
      real d,del_leaf,z_lc,hole
      real pi,coef
      real dl,drho
      real z0_mn,z0_sig

      external pxxvva_s,pzzvva_s,pxzvva_s
      external pxxvvb_s,pzzvvb_s,pxzvvb_s
      external pxxvvc_s,pzzvvc_s,pxzvvc_s
      external pxxvvd_s,pzzvvd_s,pxzvvd_s

      external pyyhha_s
      external pyyhhb_s
      external pyyhhc_s
      external pyyhhd_s

      external pxyhva_s,pyzhva_s
      external pxyhvb_s,pyzhvb_s
      external pxyhvc_s,pyzhvc_s
      external pxyhvd_s,pyzhvd_s
      
      integer s_tensor_flag,s_tau_flag,z0_flag,s_field_flag(2,4)
      common /stalk_flags/ s_tensor_flag,s_tau_flag
     &     ,z0_flag,s_field_flag      
      
      common /stalk_struct/ a,alpha_s,ell_0
      common /incident/ theta_i,phi_i,k_0
      common /layer_struct/ d,del_leaf,z_lc,hole
      common /z0_stat/ z0_mn,z0_sig
      
      common /p_store/
     &     pxxa_mn,pxya_mn,pxza_mn,pyya_mn,pyza_mn,pzza_mn,     
     &     pxxb_mn,pxyb_mn,pxzb_mn,pyyb_mn,pyzb_mn,pzzb_mn,
     &     pxxc_mn,pxyc_mn,pxzc_mn,pyyc_mn,pyzc_mn,pzzc_mn,
     &     pxxd_mn,pxyd_mn,pxzd_mn,pyyd_mn,pyzd_mn,pzzd_mn

      save /p_store/
      
      include 'error.incld'
      
      alpha_s=phi
      a=aa

      drho=a
      dl=sqrt(1.0+drho**2)

c     this has been changed to allow for efficient integration
c      from -d to z0, ell_0 is now not the length of the stalk,
c      but is instead is greater than or equal to the MAXIMUM
c      length of the stalk.
c     
c      ell_0=dl*(z+d) <- old ell_0
c
c     new ell_0 follows:

      ell_0= dl*d

      pi=3.1415927
      coef=k_0**2/(4.0*pi)
      s_tensor_flag=0

      if (pol .eq. vvvv) then

        if (mech .eq. p) then

          if (mech_c .eq. p) then
            
            if(z0_flag .eq. 0) then

              call integral(pxxvva_s,-d,z0_mn,pxxa_mn)
              call integral(pzzvva_s,-d,z0_mn,pzza_mn)
              call integral(pxzvva_s,-d,z0_mn,pxza_mn)

              z0_flag=1

            else
              
              continue
              
            endif
            
            call integral(pxxvva_s,z0_mn,z0,pxxa)
            call integral(pzzvva_s,z0_mn,z0,pzza)
            call integral(pxzvva_s,z0_mn,z0,pxza)

            pxxa=pxxa+pxxa_mn
            pzza=pzza+pzza_mn
            pxza=pxza+pxza_mn

            call phivva_g(z0, phi_g)

c            print*, z0,real(pxxa),aimag(pxxa),
c     &           real(pzza),aimag(pzza),real(pxza)
            
            stalk_phig_scat=coef*
     &           conjg(phi_g)*
     &           (cos(theta_i)**2*pxxa
     &           -2.0*cos(theta_i)*sin(theta_i)*pxza
     &           +sin(theta_i)**2*pzza)
            
            
          else if (mech_c .eq. gp) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvva_s,-d,z0_mn,pxxa_mn)
              call integral(pzzvva_s,-d,z0_mn,pzza_mn)
              call integral(pxzvva_s,-d,z0_mn,pxza_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxxvva_s,z0_mn,z0,pxxa)
            call integral(pzzvva_s,z0_mn,z0,pzza)
            call integral(pxzvva_s,z0_mn,z0,pxza)
            
            pxxa=pxxa+pxxa_mn
            pzza=pzza+pzza_mn
            pxza=pxza+pxza_mn
            
            call phivvb_g(z0, phi_g)             
            stalk_phig_scat=coef*
     &           conjg(phi_g)*
     &           (cos(theta_i)**2*pxxa
     &           -2.0*cos(theta_i)*sin(theta_i)*pxza
     &           +sin(theta_i)**2*pzza)
            
          else if (mech_c .eq. pg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvva_s,-d,z0_mn,pxxa_mn)
              call integral(pzzvva_s,-d,z0_mn,pzza_mn)
              call integral(pxzvva_s,-d,z0_mn,pxza_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxxvva_s,z0_mn,z0,pxxa)
            call integral(pzzvva_s,z0_mn,z0,pzza)
            call integral(pxzvva_s,z0_mn,z0,pxza)
            
            pxxa=pxxa+pxxa_mn
            pzza=pzza+pzza_mn
            pxza=pxza+pxza_mn
            
            call phivvc_g(z0, phi_g)
            
            stalk_phig_scat=coef*
     &           conjg(phi_g)*
     &           (cos(theta_i)**2*pxxa
     &           -2.0*cos(theta_i)*sin(theta_i)*pxza
     &           +sin(theta_i)**2*pzza)
            
          else if (mech_c .eq. gpg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvva_s,-d,z0_mn,pxxa_mn)
              call integral(pzzvva_s,-d,z0_mn,pzza_mn)
              call integral(pxzvva_s,-d,z0_mn,pxza_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxxvva_s,z0_mn,z0,pxxa)
            call integral(pzzvva_s,z0_mn,z0,pzza)
            call integral(pxzvva_s,z0_mn,z0,pxza)
            
            pxxa=pxxa+pxxa_mn
            pzza=pzza+pzza_mn
            pxza=pxza+pxza_mn
            
            call phivvd_g(z0, phi_g)
             
            stalk_phig_scat=coef*
     &           conjg(phi_g)*
     &           (cos(theta_i)**2*pxxa
     &           -2.0*cos(theta_i)*sin(theta_i)*pxza
     &           +sin(theta_i)**2*pzza)
            
          else
            
            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func stalk_phig_scat'
            code=49
            call the_death_code(die,error,code)
            
          endif
          
        else if (mech .eq. gp) then
          
          if (mech_c .eq. p) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvvb_s,-d,z0_mn,pxxb_mn)
              call integral(pzzvvb_s,-d,z0_mn,pzzb_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxxvvb_s,z0_mn,z0,pxxb)
            call integral(pzzvvb_s,z0_mn,z0,pzzb)
            
            pxxb=pxxb+pxxb_mn
            pzzb=pzzb+pzzb_mn
            
            call phivva_g(z0, phi_g)
             
                        
            stalk_phig_scat=coef*
     &           conjg(phi_g)*
     &           (-cos(theta_i)**2*pxxb+sin(theta_i)**2*pzzb)
            
          else if (mech_c .eq. gp) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvvb_s,-d,z0_mn,pxxb_mn)
              call integral(pzzvvb_s,-d,z0_mn,pzzb_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxxvvb_s,z0_mn,z0,pxxb)
            call integral(pzzvvb_s,z0_mn,z0,pzzb)
            
            pxxb=pxxb+pxxb_mn
            pzzb=pzzb+pzzb_mn
            
            call phivvb_g(z0, phi_g)
             
                        
            stalk_phig_scat=coef*
     &           conjg(phi_g)*
     &           (-cos(theta_i)**2*pxxb+sin(theta_i)**2*pzzb)
            
          else if (mech_c .eq. pg)then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvvb_s,-d,z0_mn,pxxb_mn)
              call integral(pzzvvb_s,-d,z0_mn,pzzb_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxxvvb_s,z0_mn,z0,pxxb)
            call integral(pzzvvb_s,z0_mn,z0,pzzb)
            
            pxxb=pxxb+pxxb_mn
            pzzb=pzzb+pzzb_mn
            
            call phivvc_g(z0, phi_g)
             
                        
            stalk_phig_scat=coef*
     &           conjg(phi_g)*
     &           (-cos(theta_i)**2*pxxb+sin(theta_i)**2*pzzb)
            
          else if (mech_c .eq. gpg)then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvvb_s,-d,z0_mn,pxxb_mn)
              call integral(pzzvvb_s,-d,z0_mn,pzzb_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxxvvb_s,z0_mn,z0,pxxb)
            call integral(pzzvvb_s,z0_mn,z0,pzzb)
            
            pxxb=pxxb+pxxb_mn
            pzzb=pzzb+pzzb_mn
            
            call phivvd_g(z0, phi_g)
             
                        
            stalk_phig_scat=coef*
     &           conjg(phi_g)*
     &           (-cos(theta_i)**2*pxxb+sin(theta_i)**2*pzzb)
            
          else
            
            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func stalk_phig_scat'
            code=49
            call the_death_code(die,error,code)
            
          endif
          
        else if (mech .eq. pg) then
          
          if (mech_c .eq. p) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvvc_s,-d,z0_mn,pxxc_mn)
              call integral(pzzvvc_s,-d,z0_mn,pzzc_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxxvvc_s,z0_mn,z0,pxxc)
            call integral(pzzvvc_s,z0_mn,z0,pzzc)
            
            pxxc=pxxc+pxxc_mn
            pzzc=pzzc+pzzc_mn
            
            call phivva_g(z0, phi_g)
             
                        
            stalk_phig_scat=coef*
     &           conjg(phi_g)*
     &           (-cos(theta_i)**2*pxxc+sin(theta_i)**2*pzzc)
            
          else if (mech_c .eq. gp) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvvc_s,-d,z0_mn,pxxc_mn)
              call integral(pzzvvc_s,-d,z0_mn,pzzc_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxxvvc_s,z0_mn,z0,pxxc)
            call integral(pzzvvc_s,z0_mn,z0,pzzc)
            
            pxxc=pxxc+pxxc_mn
            pzzc=pzzc+pzzc_mn
            
            call phivvb_g(z0, phi_g)
             
                        
            stalk_phig_scat=coef*
     &           conjg(phi_g)*
     &           (-cos(theta_i)**2*pxxc+sin(theta_i)**2*pzzc)
            
          else if (mech_c .eq. pg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvvc_s,-d,z0_mn,pxxc_mn)
              call integral(pzzvvc_s,-d,z0_mn,pzzc_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxxvvc_s,z0_mn,z0,pxxc)
            call integral(pzzvvc_s,z0_mn,z0,pzzc)
            
            pxxc=pxxc+pxxc_mn
            pzzc=pzzc+pzzc_mn
            
            call phivvc_g(z0, phi_g)
                        
            stalk_phig_scat=coef*
     &           conjg(phi_g)*
     &           (-cos(theta_i)**2*pxxc+sin(theta_i)**2*pzzc)
            
            else if (mech_c .eq. gpg) then
                
              if(z0_flag .eq. 0) then
                
                call integral(pxxvvc_s,-d,z0_mn,pxxc_mn)
                call integral(pzzvvc_s,-d,z0_mn,pzzc_mn)
                
                z0_flag=1
                
              else
                
                continue
                
              endif
              
              call integral(pxxvvc_s,z0_mn,z0,pxxc)
              call integral(pzzvvc_s,z0_mn,z0,pzzc)
              
              pxxc=pxxc+pxxc_mn
              pzzc=pzzc+pzzc_mn
              
              call phivvd_g(z0, phi_g)
                          
              stalk_phig_scat=coef*
     &             conjg(phi_g)*
     &             (-cos(theta_i)**2*pxxc+sin(theta_i)**2*pzzc)
              
            else
              
              write(error,*)'mech_c = ',mech_c,' not recognized',
     &             ' in func stalk_phig_scat'
              code=49
              call the_death_code(die,error,code)
              
            endif
            
          else if (mech .eq. gpg) then
            
            if (mech_c .eq. p) then
              
              if(z0_flag .eq. 0) then
                
                call integral(pxxvvd_s,-d,z0_mn,pxxd_mn)
                call integral(pzzvvd_s,-d,z0_mn,pzzd_mn)
                call integral(pxzvvd_s,-d,z0_mn,pxzd_mn)
                
                z0_flag=1
                
              else
                
                continue
                
              endif
              
              call integral(pxxvvd_s,z0_mn,z0,pxxd)
              call integral(pzzvvd_s,z0_mn,z0,pzzd)
              call integral(pxzvvd_s,z0_mn,z0,pxzd)
              
              pxxd=pxxd+pxxd_mn
              pzzd=pzzd+pzzd_mn
              pxzd=pxzd+pxzd_mn

              call phivva_g(z0, phi_g)
                          
              stalk_phig_scat=coef*
     &             conjg(phi_g)*
     &             (cos(theta_i)**2*pxxd
     &             +2.0*cos(theta_i)*sin(theta_i)*pxzd
     &             +sin(theta_i)**2*pzzd)
              
              
            else if (mech_c .eq. gp) then
              
              if(z0_flag .eq. 0) then
                
                call integral(pxxvvd_s,-d,z0_mn,pxxd_mn)
                call integral(pzzvvd_s,-d,z0_mn,pzzd_mn)
                call integral(pxzvvd_s,-d,z0_mn,pxzd_mn)
                
                z0_flag=1
                
              else
                
                continue
                
              endif
              
              call integral(pxxvvd_s,z0_mn,z0,pxxd)
              call integral(pzzvvd_s,z0_mn,z0,pzzd)
              call integral(pxzvvd_s,z0_mn,z0,pxzd)
              
              pxxd=pxxd+pxxd_mn
              pzzd=pzzd+pzzd_mn
              pxzd=pxzd+pxzd_mn
              
              call phivvb_g(z0, phi_g)
             
                          
              stalk_phig_scat=coef*
     &             conjg(phi_g)*
     &             (cos(theta_i)**2*pxxd
     &             +2.0*cos(theta_i)*sin(theta_i)*pxzd
     &             +sin(theta_i)**2*pzzd)
              
            else if (mech_c .eq. pg) then
              
              if(z0_flag .eq. 0) then
                
                call integral(pxxvvd_s,-d,z0_mn,pxxd_mn)
                call integral(pzzvvd_s,-d,z0_mn,pzzd_mn)
                call integral(pxzvvd_s,-d,z0_mn,pxzd_mn)
                
                z0_flag=1
                
              else
                
                continue
                
              endif
              
              call integral(pxxvvd_s,z0_mn,z0,pxxd)
              call integral(pzzvvd_s,z0_mn,z0,pzzd)
              call integral(pxzvvd_s,z0_mn,z0,pxzd)
              
              pxxd=pxxd+pxxd_mn
              pzzd=pzzd+pzzd_mn
              pxzd=pxzd+pxzd_mn
              
              call phivvc_g(z0, phi_g)
             
                          
              stalk_phig_scat=coef*
     &             conjg(phi_g)*
     &             (cos(theta_i)**2*pxxd
     &             +2.0*cos(theta_i)*sin(theta_i)*pxzd
     &             +sin(theta_i)**2*pzzd)
              
            else if (mech_c .eq. gpg) then
              
              if(z0_flag .eq. 0) then
                
                call integral(pxxvvd_s,-d,z0_mn,pxxd_mn)
                call integral(pzzvvd_s,-d,z0_mn,pzzd_mn)
                call integral(pxzvvd_s,-d,z0_mn,pxzd_mn)
                
                z0_flag=1
                
              else
                
                continue
                
              endif
              
              call integral(pxxvvd_s,z0_mn,z0,pxxd)
              call integral(pzzvvd_s,z0_mn,z0,pzzd)
              call integral(pxzvvd_s,z0_mn,z0,pxzd)
              
              pxxd=pxxd+pxxd_mn
              pzzd=pzzd+pzzd_mn
              pxzd=pxzd+pxzd_mn
              
              call phivvd_g(z0, phi_g)
             
                          
              stalk_phig_scat=coef*
     &             conjg(phi_g)*
     &             (cos(theta_i)**2*pxxd
     &             +2.0*cos(theta_i)*sin(theta_i)*pxzd
     &             +sin(theta_i)**2*pzzd)
              
            else
              
              write(error,*)'mech_c = ',mech_c,' not recognized',
     &             ' in func stalk_phig_scat'
              code=49
              call the_death_code(die,error,code)
              
            endif
            
          else
            
            write(error,*)'mech = ',mech,' not recognized',
     &           ' in func stalk_phig_scat'
            code=50
            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
              
              if(z0_flag .eq. 0) then
                
                call integral(pyyhha_s,-d,z0_mn,pyya_mn)
                
                z0_flag=1
                
              else
                
                continue
                
              endif
              
              call integral(pyyhha_s,z0_mn,z0,pyya)
              
              pyya=pyya+pyya_mn
              
              call phihha_g(z0, phi_g)            
             
                          
              stalk_phig_scat=coef*
     &             pyya*conjg(phi_g)
              
            else if (mech_c .eq. gp) then
              
              if(z0_flag .eq. 0) then
                
                call integral(pyyhha_s,-d,z0_mn,pyya_mn)
                
                z0_flag=1
                
              else
                
                continue
                
              endif
              
              call integral(pyyhha_s,z0_mn,z0,pyya)
              
              pyya=pyya+pyya_mn
              
              call phihhb_g(z0, phi_g)
             
                          
              stalk_phig_scat=coef*
     &             pyya*conjg(phi_g)
              
            else if (mech_c .eq. pg) then
              
              if(z0_flag .eq. 0) then
                
                call integral(pyyhha_s,-d,z0_mn,pyya_mn)
                
                z0_flag=1
                
              else
                
                continue
                
              endif
              
              call integral(pyyhha_s,z0_mn,z0,pyya)
              
              pyya=pyya+pyya_mn
              
              call phihhc_g(z0, phi_g)  
             
                          
              stalk_phig_scat=coef*
     &             pyya*conjg(phi_g)
              
            else if (mech_c .eq. gpg) then
              
              if(z0_flag .eq. 0) then
                
                call integral(pyyhha_s,-d,z0_mn,pyya_mn)
                
                z0_flag=1
                
              else
                
                continue
                
              endif
              
              call integral(pyyhha_s,z0_mn,z0,pyya)
              
              pyya=pyya+pyya_mn
              
              call phihhd_g(z0, phi_g)
             
                          
              stalk_phig_scat=coef*
     &             pyya*conjg(phi_g)            
              
            else
              
              write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func stalk_phig_scat'
              code=49
              call the_death_code(die,error,code)
              
            endif
            
          else if (mech .eq. gp) then
            
            if (mech_c .eq. p) then
              
              if(z0_flag .eq. 0) then
                
                call integral(pyyhhb_s,-d,z0_mn,pyyb_mn)
                
                z0_flag=1
                
              else
                
                continue
                
              endif
              
              call integral(pyyhhb_s,z0_mn,z0,pyyb)
              
              pyyb=pyyb+pyyb_mn
              
              call phihha_g(z0, phi_g)            
             
                          
              stalk_phig_scat=coef*
     &             pyyb*conjg(phi_g)
              
            else if (mech_c .eq. gp) then
              
              if(z0_flag .eq. 0) then
                
                call integral(pyyhhb_s,-d,z0_mn,pyyb_mn)
                
                z0_flag=1
                
              else
                
                continue
                
              endif
              
              call integral(pyyhhb_s,z0_mn,z0,pyyb)
              
              pyyb=pyyb+pyyb_mn
                          
              call phihhb_g(z0, phi_g)
             
                          
              stalk_phig_scat=coef*
     &             pyyb*conjg(phi_g)
              
            else if (mech_c .eq. pg) then
              
              if(z0_flag .eq. 0) then
                
                call integral(pyyhhb_s,-d,z0_mn,pyyb_mn)
                
                z0_flag=1
                
              else
                
                continue
                
              endif
              
              call integral(pyyhhb_s,z0_mn,z0,pyyb)
              
              pyyb=pyyb+pyyb_mn
              
              call phihhc_g(z0, phi_g)  
             
            
              stalk_phig_scat=coef*
     &             pyyb*conjg(phi_g)
              
            else if (mech_c .eq. gpg) then
              
              if(z0_flag .eq. 0) then
                
                call integral(pyyhhb_s,-d,z0_mn,pyyb_mn)
                
                z0_flag=1
                
              else
                
                continue
                
              endif
              
              call integral(pyyhhb_s,z0_mn,z0,pyyb)
              
              pyyb=pyyb+pyyb_mn
              
              s_tensor_flag=0
              call integral(pyyhhb_s,-d,z0,pyyb)
              
              call phihhd_g(z0, phi_g)
             
                          
              stalk_phig_scat=coef*
     &             pyyb*conjg(phi_g)            
              
            else
              
              write(error,*)'mech_c = ',mech_c,' not recognized',
     &             ' in func stalk_phig_scat'
              code=49
              call the_death_code(die,error,code)
              
            endif
            
          else if (mech .eq. pg) then
            
            if (mech_c .eq. p) then
              
              if(z0_flag .eq. 0) then
                
                call integral(pyyhhc_s,-d,z0_mn,pyyc_mn)
                
                z0_flag=1
                
              else
                
                continue
                
              endif
              
              call integral(pyyhhc_s,z0_mn,z0,pyyc)
              
              pyyc=pyyc+pyyc_mn
              
              call phihha_g(z0, phi_g)            
             
                          
            stalk_phig_scat=coef*
     &             pyyc*conjg(phi_g)
            
          else if (mech_c .eq. gp) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pyyhhc_s,-d,z0_mn,pyyc_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhhc_s,z0_mn,z0,pyyc)
            
            pyyc=pyyc+pyyc_mn
            
            call phihhb_g(z0, phi_g)
             
                        
            stalk_phig_scat=coef*
     &           pyyc*conjg(phi_g)
            
          else if (mech_c .eq. pg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pyyhhc_s,-d,z0_mn,pyyc_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhhc_s,z0_mn,z0,pyyc)
            
            pyyc=pyyc+pyyc_mn
            
            call phihhc_g(z0, phi_g)  
             
                        
            stalk_phig_scat=coef*
     &           pyyc*conjg(phi_g)
            
          else if (mech_c .eq. gpg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pyyhhc_s,-d,z0_mn,pyyc_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhhc_s,z0_mn,z0,pyyc)
            
            pyyc=pyyc+pyyc_mn
            
            call phihhd_g(z0, phi_g)
             
                        
            stalk_phig_scat=coef*
     &           pyyc*conjg(phi_g)            
            
          else
            
            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func stalk_phig_scat'
            code=49
            call the_death_code(die,error,code)
            
          endif
          
        else if (mech .eq. gpg) then
          
          if (mech_c .eq. p) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pyyhhd_s,-d,z0_mn,pyyd_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhhd_s,z0_mn,z0,pyyd)
            
            pyyd=pyyd+pyyd_mn
            
            call phihha_g(z0, phi_g)            
             
                        
            stalk_phig_scat=coef*
     &           pyyd*conjg(phi_g)
            
          else if (mech_c .eq. gp) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pyyhhd_s,-d,z0_mn,pyyd_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhhd_s,z0_mn,z0,pyyd)
            
            pyyd=pyyd+pyyd_mn
            
            call phihhb_g(z0, phi_g)
             
                        
            stalk_phig_scat=coef*
     &           pyyd*conjg(phi_g)
            
          else if (mech_c .eq. pg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pyyhhd_s,-d,z0_mn,pyyd_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhhd_s,z0_mn,z0,pyyd)
            
            pyyd=pyyd+pyyd_mn
            
            call phihhc_g(z0, phi_g)  
             
                        
            stalk_phig_scat=coef*
     &           pyyd*conjg(phi_g)
            
          else if (mech_c .eq. gpg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pyyhhd_s,-d,z0_mn,pyyd_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhhd_s,z0_mn,z0,pyyd)
            
            pyyd=pyyd+pyyd_mn
            
            call phihhd_g(z0, phi_g)
             
                        
            stalk_phig_scat=coef*
     &           pyyd*conjg(phi_g)            
            
          else
            
            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func stalk_phig_scat'
            code=49
            call the_death_code(die,error,code)
            
          endif
          
        else
          
          write(error,*)'mech = ',mech,' not recognized',
     &         ' in func stalk_phig_scat'
          code=50
          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
            
            if(z0_flag .eq. 0) then
              
              call integral(pxyhva_s,-d,z0_mn,pxya_mn)
              call integral(pyzhva_s,-d,z0_mn,pyza_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxyhva_s,z0_mn,z0,pxya)
            call integral(pyzhva_s,z0_mn,z0,pyza)
            
            pxya=pxya+pxya_mn
            pyza=pyza+pyza_mn
            
            call phihva_g(z0, phi_g)            
             
                        
            stalk_phig_scat=coef*
     &           (cos(theta_i)*pxya-sin(theta_i)*pyza)*conjg(phi_g)
            
          else if (mech_c .eq. gp) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxyhva_s,-d,z0_mn,pxya_mn)
              call integral(pyzhva_s,-d,z0_mn,pyza_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxyhva_s,z0_mn,z0,pxya)
            call integral(pyzhva_s,z0_mn,z0,pyza)
            
            pxya=pxya+pxya_mn
            pyza=pyza+pyza_mn
            
            call phihvb_g(z0, phi_g)            
             
                        
            stalk_phig_scat=coef*
     &           (cos(theta_i)*pxya-sin(theta_i)*pyza)*conjg(phi_g)
            
            
          else if (mech_c .eq. pg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxyhva_s,-d,z0_mn,pxya_mn)
              call integral(pyzhva_s,-d,z0_mn,pyza_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxyhva_s,z0_mn,z0,pxya)
            call integral(pyzhva_s,z0_mn,z0,pyza)
            
            pxya=pxya+pxya_mn
            pyza=pyza+pyza_mn
            
            call phihvc_g(z0, phi_g)            
             
                        
            stalk_phig_scat=coef*
     &           (cos(theta_i)*pxya-sin(theta_i)*pyza)*conjg(phi_g)
            
          else if (mech_c .eq. gpg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxyhva_s,-d,z0_mn,pxya_mn)
              call integral(pyzhva_s,-d,z0_mn,pyza_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxyhva_s,z0_mn,z0,pxya)
            call integral(pyzhva_s,z0_mn,z0,pyza)
            
            pxya=pxya+pxya_mn
            pyza=pyza+pyza_mn
            
            call phihvd_g(z0, phi_g)            
             
                        
            stalk_phig_scat=coef*
     &           (cos(theta_i)*pxya-sin(theta_i)*pyza)*conjg(phi_g)
            
          else
            
            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func stalk_phig_scat'
            code=49
           call the_death_code(die,error,code)
           
         endif
         
       else if (mech .eq. gp) then
         
         if (mech_c .eq. p) then
           
           if(z0_flag .eq. 0) then
             
             call integral(pxyhvb_s,-d,z0_mn,pxyb_mn)
             call integral(pyzhvb_s,-d,z0_mn,pyzb_mn)
             
             z0_flag=1
             
           else
             
             continue
             
           endif
           
           call integral(pxyhvb_s,z0_mn,z0,pxyb)
           call integral(pyzhvb_s,z0_mn,z0,pyzb)
           
           pxyb=pxyb+pxyb_mn
           pyzb=pyzb+pyzb_mn
           
           call phihva_g(z0, phi_g)            
             
                       
           stalk_phig_scat=coef*
     &          (-cos(theta_i)*pxyb-sin(theta_i)*pyzb)*conjg(phi_g)
           
         else if (mech_c .eq. gp) then
           
           if(z0_flag .eq. 0) then
             
             call integral(pxyhvb_s,-d,z0_mn,pxyb_mn)
             call integral(pyzhvb_s,-d,z0_mn,pyzb_mn)
             
             z0_flag=1
             
           else
             
             continue
             
           endif
           
            call integral(pxyhvb_s,z0_mn,z0,pxyb)
            call integral(pyzhvb_s,z0_mn,z0,pyzb)
            
            pxyb=pxyb+pxyb_mn
            pyzb=pyzb+pyzb_mn
            
            call phihvb_g(z0, phi_g)            
             
                        
            stalk_phig_scat=coef*
     &           (-cos(theta_i)*pxyb-sin(theta_i)*pyzb)*conjg(phi_g)
            
            
          else if (mech_c .eq. pg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxyhvb_s,-d,z0_mn,pxyb_mn)
              call integral(pyzhvb_s,-d,z0_mn,pyzb_mn)
              
              z0_flag=1
              
            else
              
              continue
               
            endif
            
            call integral(pxyhvb_s,z0_mn,z0,pxyb)
            call integral(pyzhvb_s,z0_mn,z0,pyzb)
            
            pxyb=pxyb+pxyb_mn
            pyzb=pyzb+pyzb_mn
            
            call phihvc_g(z0, phi_g)            
             
                        
            stalk_phig_scat=coef*
     &          (-cos(theta_i)*pxyb-sin(theta_i)*pyzb)*conjg(phi_g)
            
          else if (mech_c .eq. gpg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxyhvb_s,-d,z0_mn,pxyb_mn)
              call integral(pyzhvb_s,-d,z0_mn,pyzb_mn)
              
              z0_flag=1
             
            else
              
              continue
              
            endif
            
            call integral(pxyhvb_s,z0_mn,z0,pxyb)
            call integral(pyzhvb_s,z0_mn,z0,pyzb)
            
            pxyb=pxyb+pxyb_mn
            pyzb=pyzb+pyzb_mn
            
            call phihvd_g(z0, phi_g)            
             
                        
            stalk_phig_scat=coef*
     &           (-cos(theta_i)*pxyb-sin(theta_i)*pyzb)*conjg(phi_g)
            
          else
            
            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func stalk_phig_scat'
            code=49
            call the_death_code(die,error,code)
            
          endif
          
        else if (mech .eq. pg) then
          
          if (mech_c .eq. p) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxyhvc_s,-d,z0_mn,pxyc_mn)
              call integral(pyzhvc_s,-d,z0_mn,pyzc_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxyhvc_s,z0_mn,z0,pxyc)
            call integral(pyzhvc_s,z0_mn,z0,pyzc)
            
            pxyc=pxyc+pxyc_mn
            pyzc=pyzc+pyzc_mn
            
            call phihva_g(z0, phi_g)            
             
                                    
            stalk_phig_scat=coef*
     &           (cos(theta_i)*pxyc-sin(theta_i)*pyzc)*conjg(phi_g)
            
          else if (mech_c .eq. gp) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxyhvc_s,-d,z0_mn,pxyc_mn)
              call integral(pyzhvc_s,-d,z0_mn,pyzc_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxyhvc_s,z0_mn,z0,pxyc)
            call integral(pyzhvc_s,z0_mn,z0,pyzc)
            
            pxyc=pxyc+pxyc_mn
            pyzc=pyzc+pyzc_mn
            
            call phihvb_g(z0, phi_g)            
             
                                    
            stalk_phig_scat=coef*
     &           (cos(theta_i)*pxyc-sin(theta_i)*pyzc)*conjg(phi_g)
            
            
          else if (mech_c .eq. pg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxyhvc_s,-d,z0_mn,pxyc_mn)
              call integral(pyzhvc_s,-d,z0_mn,pyzc_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxyhvc_s,z0_mn,z0,pxyc)
            call integral(pyzhvc_s,z0_mn,z0,pyzc)
            
            pxyc=pxyc+pxyc_mn
            pyzc=pyzc+pyzc_mn
            
            call phihvc_g(z0, phi_g)            
             
                                    
            stalk_phig_scat=coef*
     &           (cos(theta_i)*pxyc-sin(theta_i)*pyzc)*conjg(phi_g)
            
          else if (mech_c .eq. gpg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxyhvc_s,-d,z0_mn,pxyc_mn)
              call integral(pyzhvc_s,-d,z0_mn,pyzc_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxyhvc_s,z0_mn,z0,pxyc)
            call integral(pyzhvc_s,z0_mn,z0,pyzc)
            
            pxyc=pxyc+pxyc_mn
            pyzc=pyzc+pyzc_mn
            
            call phihvd_g(z0, phi_g)            
             
                                    
            stalk_phig_scat=coef*
     &           (cos(theta_i)*pxyc-sin(theta_i)*pyzc)*conjg(phi_g)
            
          else
            
            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func stalk_phig_scat'
            code=49
            call the_death_code(die,error,code)
            
          endif
          
        else if (mech .eq. gpg) then
          
          if (mech_c .eq. p) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxyhvd_s,-d,z0_mn,pxyd_mn)
              call integral(pyzhvd_s,-d,z0_mn,pyzd_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxyhvd_s,z0_mn,z0,pxyd)
            call integral(pyzhvd_s,z0_mn,z0,pyzd)
            
            pxyd=pxyd+pxyd_mn
            pyzd=pyzd+pyzd_mn
            
            call phihva_g(z0, phi_g)            
             
                                    
            stalk_phig_scat=coef*
     &           (-cos(theta_i)*pxyd-sin(theta_i)*pyzd)*conjg(phi_g)
            
          else if (mech_c .eq. gp) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxyhvd_s,-d,z0_mn,pxyd_mn)
              call integral(pyzhvd_s,-d,z0_mn,pyzd_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxyhvd_s,z0_mn,z0,pxyd)
            call integral(pyzhvd_s,z0_mn,z0,pyzd)
            
            pxyd=pxyd+pxyd_mn
            pyzd=pyzd+pyzd_mn
            
            call phihvb_g(z0, phi_g)            
             
                                    
            stalk_phig_scat=coef*
     &           (-cos(theta_i)*pxyd-sin(theta_i)*pyzd)*conjg(phi_g)
            
            
          else if (mech_c .eq. pg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxyhvd_s,-d,z0_mn,pxyd_mn)
              call integral(pyzhvd_s,-d,z0_mn,pyzd_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxyhvd_s,z0_mn,z0,pxyd)
            call integral(pyzhvd_s,z0_mn,z0,pyzd)
            
            pxyd=pxyd+pxyd_mn
            pyzd=pyzd+pyzd_mn
            
            call phihvc_g(z0, phi_g)            
             
                                    
            stalk_phig_scat=coef*
     &           (-cos(theta_i)*pxyd-sin(theta_i)*pyzd)*conjg(phi_g)
            
          else if (mech_c .eq. gpg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxyhvd_s,-d,z0_mn,pxyd_mn)
              call integral(pyzhvd_s,-d,z0_mn,pyzd_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxyhvd_s,z0_mn,z0,pxyd)
            call integral(pyzhvd_s,z0_mn,z0,pyzd)
            
            pxyd=pxyd+pxyd_mn
            pyzd=pyzd+pyzd_mn
            
            call phihvd_g(z0, phi_g)            
             
                                    
            stalk_phig_scat=coef*
     &           (-cos(theta_i)*pxyd-sin(theta_i)*pyzd)*conjg(phi_g)
            
          else
            
            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func stalk_phig_scat'
            code=49
            call the_death_code(die,error,code)
            
          endif
          
        else
          
          write(error,*)'mech = ',mech,' not recognized',
     &         ' in func stalk_phig_scat'
          code=50
          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
            
            if(z0_flag .eq. 0) then
              
              call integral(pyyhha_s,-d,z0_mn,pyya_mn)               
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhha_s,z0_mn,z0,pyya)               
            
            pyya=pyya+pyya_mn
            
            call phivva_g(z0, phi_g)            
             
                                    
            stalk_phig_scat=coef*
     &           pyya*conjg(phi_g)
            
          else if (mech_c .eq. gp) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pyyhha_s,-d,z0_mn,pyya_mn)               
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhha_s,z0_mn,z0,pyya)               
            
            pyya=pyya+pyya_mn
            
            call phivvb_g(z0, phi_g)
             
                                    
            stalk_phig_scat=coef*
     &           pyya*conjg(phi_g)
            
          else if (mech_c .eq. pg) then
           
            if(z0_flag .eq. 0) then
              
              call integral(pyyhha_s,-d,z0_mn,pyya_mn)               
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhha_s,z0_mn,z0,pyya)               
            
            pyya=pyya+pyya_mn
            
            call phivvc_g(z0, phi_g)  
             
                                    
            stalk_phig_scat=coef*
     &           pyya*conjg(phi_g)
            
          else if (mech_c .eq. gpg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pyyhha_s,-d,z0_mn,pyya_mn)               
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            pxxa=pxxa+pxxa_mn
            
            call phivvd_g(z0, phi_g)
             
                                    
            stalk_phig_scat=coef*
     &           pyya*conjg(phi_g)            
            
          else
            
            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func stalk_phig_scat'
            code=49
            call the_death_code(die,error,code)
            
          endif
          
        else if (mech .eq. gp) then
          
          if (mech_c .eq. p) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pyyhhb_s,-d,z0_mn,pyyb_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhhb_s,z0_mn,z0,pyyb)
            
            pyyb=pyyb+pyyb_mn
            
            call phivva_g(z0, phi_g)            
             
                        
            stalk_phig_scat=coef*
     &           pyyb*conjg(phi_g)
            
          else if (mech_c .eq. gp) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pyyhhb_s,-d,z0_mn,pyyb_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhhb_s,z0_mn,z0,pyyb)
            
            pyyb=pyyb+pyyb_mn
            
            call phivvb_g(z0, phi_g)
             
                                    
            stalk_phig_scat=coef*
     &           pyyb*conjg(phi_g)
            
          else if (mech_c .eq. pg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pyyhhb_s,-d,z0_mn,pyyb_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhhb_s,z0_mn,z0,pyyb)
            
            pyyb=pyyb+pyyb_mn
            
            call phivvc_g(z0, phi_g)  
             
                                    
            stalk_phig_scat=coef*
     &           pyyb*conjg(phi_g)
            
          else if (mech_c .eq. gpg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pyyhhb_s,-d,z0_mn,pyyb_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhhb_s,z0_mn,z0,pyyb)
            
            pyyb=pyyb+pyyb_mn
            
            call phivvd_g(z0, phi_g)
             
                                    
            stalk_phig_scat=coef*
     &           pyyb*conjg(phi_g)            
            
          else
            
            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func stalk_phig_scat'
            code=49
            call the_death_code(die,error,code)
            
          endif
          
        else if (mech .eq. pg) then
          
          if (mech_c .eq. p) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pyyhhc_s,-d,z0_mn,pyyc_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhhc_s,z0_mn,z0,pyyc)
            
            pyyc=pyyc+pyyc_mn
            
            call phivva_g(z0, phi_g)            
             
                                    
            stalk_phig_scat=coef*
     &           pyyc*conjg(phi_g)
            
          else if (mech_c .eq. gp) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pyyhhc_s,-d,z0_mn,pyyc_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhhc_s,z0_mn,z0,pyyc)
            
            pyyc=pyyc+pyyc_mn
            
            call phivvb_g(z0, phi_g)
             
                                    
            stalk_phig_scat=coef*
     &           pyyc*conjg(phi_g)
            
          else if (mech_c .eq. pg) then
           
            if(z0_flag .eq. 0) then
              
              call integral(pyyhhc_s,-d,z0_mn,pyyc_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhhc_s,z0_mn,z0,pyyc)
            
            pyyc=pyyc+pyyc_mn
            
            call phivvc_g(z0, phi_g)  
             
                                    
            stalk_phig_scat=coef*
     &           pyyc*conjg(phi_g)
            
          else if (mech_c .eq. gpg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pyyhhc_s,-d,z0_mn,pyyc_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhhc_s,z0_mn,z0,pyyc)
            
            pyyc=pyyc+pyyc_mn
            
            call phivvd_g(z0, phi_g)
             
                                    
            stalk_phig_scat=coef*
     &           pyyc*conjg(phi_g)            
            
          else
            
            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func stalk_phig_scat'
            code=49
            call the_death_code(die,error,code)
            
          endif
          
        else if (mech .eq. gpg) then
          
          if (mech_c .eq. p) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pyyhhd_s,-d,z0_mn,pyyd_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhhd_s,z0_mn,z0,pyyd)
            
            pyyd=pyyd+pyyd_mn
            
            call phivva_g(z0, phi_g)
             
            
            stalk_phig_scat=coef*
     &           pyyd*conjg(phi_g)
            
          else if (mech_c .eq. gp) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pyyhhd_s,-d,z0_mn,pyyd_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhhd_s,z0_mn,z0,pyyd)
            
            pyyd=pyyd+pyyd_mn
            
            call phivvb_g(z0, phi_g)
             
                                    
            stalk_phig_scat=coef*
     &           pyyd*conjg(phi_g)
            
          else if (mech_c .eq. pg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pyyhhd_s,-d,z0_mn,pyyd_mn)
              
              z0_flag=1

            else
              
              continue
              
            endif
            
            call integral(pyyhhd_s,z0_mn,z0,pyyd)
            
            pyyd=pyyd+pyyd_mn
            
            call phivvc_g(z0, phi_g)  
             
                                    
            stalk_phig_scat=coef*
     &           pyyd*conjg(phi_g)
            
          else if (mech_c .eq. gpg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pyyhhd_s,-d,z0_mn,pyyd_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pyyhhd_s,z0_mn,z0,pyyd)
            
            pyyd=pyyd+pyyd_mn
            
            call phivvd_g(z0, phi_g)
             
                                    
            stalk_phig_scat=coef*
     &           pyyd*conjg(phi_g)            
            
          else
            
            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func stalk_phig_scat'
            code=49
            call the_death_code(die,error,code)
            
          endif
          
        else
          
          write(error,*)'mech = ',mech,' not recognized',
     &         ' in func stalk_phig_scat'
          code=50
          call the_death_code(die,error,code)
          
        endif
        
        
      else if (pol .eq. vvhh) then
        
        if (mech .eq. p) then
          
          if (mech_c .eq. p) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvva_s,-d,z0_mn,pxxa_mn)
              call integral(pzzvva_s,-d,z0_mn,pzza_mn)
              call integral(pxzvva_s,-d,z0_mn,pxza_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxxvva_s,z0_mn,z0,pxxa)
            call integral(pzzvva_s,z0_mn,z0,pzza)
            call integral(pxzvva_s,z0_mn,z0,pxza)
            
            pxxa=pxxa+pxxa_mn
            pzza=pzza+pzza_mn
            pxza=pxza+pxza_mn
            
            call phihha_g(z0, phi_g)
             
                                    
            stalk_phig_scat=coef*
     &           conjg(phi_g)*
     &           (cos(theta_i)**2*pxxa
     &           -2.0*cos(theta_i)*sin(theta_i)*pxza
     &           +sin(theta_i)**2*pzza)
            
            
          else if (mech_c .eq. gp) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvva_s,-d,z0_mn,pxxa_mn)
              call integral(pzzvva_s,-d,z0_mn,pzza_mn)
              call integral(pxzvva_s,-d,z0_mn,pxza_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxxvva_s,z0_mn,z0,pxxa)
            call integral(pzzvva_s,z0_mn,z0,pzza)
            call integral(pxzvva_s,z0_mn,z0,pxza)
            
            pxxa=pxxa+pxxa_mn
            pzza=pzza+pzza_mn
            pxza=pxza+pxza_mn
            
            call phihhb_g(z0, phi_g)
             
                                    
            stalk_phig_scat=coef*
     &           conjg(phi_g)*
     &           (cos(theta_i)**2*pxxa
     &           -2.0*cos(theta_i)*sin(theta_i)*pxza
     &           +sin(theta_i)**2*pzza)
            
          else if (mech_c .eq. pg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvva_s,-d,z0_mn,pxxa_mn)
              call integral(pzzvva_s,-d,z0_mn,pzza_mn)
              call integral(pxzvva_s,-d,z0_mn,pxza_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxxvva_s,z0_mn,z0,pxxa)
            call integral(pzzvva_s,z0_mn,z0,pzza)
            call integral(pxzvva_s,z0_mn,z0,pxza)
            
            pxxa=pxxa+pxxa_mn
            pzza=pzza+pzza_mn
            pxza=pxza+pxza_mn
            
            call phihhc_g(z0, phi_g)
             
                                    
            stalk_phig_scat=coef*
     &           conjg(phi_g)*
     &           (cos(theta_i)**2*pxxa
     &           -2.0*cos(theta_i)*sin(theta_i)*pxza
     &           +sin(theta_i)**2*pzza)
            
          else if (mech_c .eq. gpg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvva_s,-d,z0_mn,pxxa_mn)
              call integral(pzzvva_s,-d,z0_mn,pzza_mn)
              call integral(pxzvva_s,-d,z0_mn,pxza_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxxvva_s,z0_mn,z0,pxxa)
            call integral(pzzvva_s,z0_mn,z0,pzza)
            call integral(pxzvva_s,z0_mn,z0,pxza)
            
            pxxa=pxxa+pxxa_mn
            pzza=pzza+pzza_mn
            pxza=pxza+pxza_mn
            
            call phihhd_g(z0, phi_g)
             
                                    
            stalk_phig_scat=coef*
     &           conjg(phi_g)*
     &           (cos(theta_i)**2*pxxa
     &           -2.0*cos(theta_i)*sin(theta_i)*pxza
     &           +sin(theta_i)**2*pzza)
            
          else
            
            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func stalk_phig_scat'
            code=49
            call the_death_code(die,error,code)
            
          endif
          
        else if (mech .eq. gp) then
          
          if (mech_c .eq. p) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvvb_s,-d,z0_mn,pxxb_mn)
              call integral(pzzvvb_s,-d,z0_mn,pzzb_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxxvvb_s,z0_mn,z0,pxxb)
            call integral(pzzvvb_s,z0_mn,z0,pzzb)
            
            pxxb=pxxb+pxxb_mn
            pzzb=pzzb+pzzb_mn
            
            call phihha_g(z0, phi_g)
             
                                    
            stalk_phig_scat=coef*
     &           conjg(phi_g)*
     &           (-cos(theta_i)**2*pxxb+sin(theta_i)**2*pzzb)
            
          else if (mech_c .eq. gp) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvvb_s,-d,z0_mn,pxxb_mn)
              call integral(pzzvvb_s,-d,z0_mn,pzzb_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxxvvb_s,z0_mn,z0,pxxb)
            call integral(pzzvvb_s,z0_mn,z0,pzzb)
            
            pxxb=pxxb+pxxb_mn
            pzzb=pzzb+pzzb_mn
            
            call phihhb_g(z0, phi_g)
             
                                    
            stalk_phig_scat=coef*
     &           conjg(phi_g)*
     &           (-cos(theta_i)**2*pxxb+sin(theta_i)**2*pzzb)
            
          else if (mech_c .eq. pg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvvb_s,-d,z0_mn,pxxb_mn)
              call integral(pzzvvb_s,-d,z0_mn,pzzb_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxxvvb_s,z0_mn,z0,pxxb)
            call integral(pzzvvb_s,z0_mn,z0,pzzb)
            
            pxxb=pxxb+pxxb_mn
            pzzb=pzzb+pzzb_mn
            
            call phihhc_g(z0, phi_g)
             
                                    
            stalk_phig_scat=coef*
     &           conjg(phi_g)*
     &           (-cos(theta_i)**2*pxxb+sin(theta_i)**2*pzzb)
            
          else if (mech_c .eq. gpg) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvvb_s,-d,z0_mn,pxxb_mn)
              call integral(pzzvvb_s,-d,z0_mn,pzzb_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxxvvb_s,z0_mn,z0,pxxb)
            call integral(pzzvvb_s,z0_mn,z0,pzzb)
            
            pxxb=pxxb+pxxb_mn
            pzzb=pzzb+pzzb_mn
            
            call phihhd_g(z0, phi_g)
             
                                    
            stalk_phig_scat=coef*
     &           conjg(phi_g)*
     &           (-cos(theta_i)**2*pxxb+sin(theta_i)**2*pzzb)
            
          else
            
            write(error,*)'mech_c = ',mech_c,' not recognized',
     &           ' in func stalk_phig_scat'
            code=49
            call the_death_code(die,error,code)
            
            
          endif
          
        else if (mech .eq. pg) then
          
          if (mech_c .eq. p) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvvc_s,-d,z0_mn,pxxc_mn)
              call integral(pzzvvc_s,-d,z0_mn,pzzc_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
            call integral(pxxvvc_s,z0_mn,z0,pxxc)
            call integral(pzzvvc_s,z0_mn,z0,pzzc)
            
            pxxc=pxxc+pxxc_mn
            pzzc=pzzc+pzzc_mn
            
            call phihha_g(z0, phi_g)
             
                                    
            stalk_phig_scat=coef*
     &           conjg(phi_g)*
     &           (-cos(theta_i)**2*pxxc+sin(theta_i)**2*pzzc)
            
          else if (mech_c .eq. gp) then
            
            if(z0_flag .eq. 0) then
              
              call integral(pxxvvc_s,-d,z0_mn,pxxc_mn)
              call integral(pzzvvc_s,-d,z0_mn,pzzc_mn)
              
              z0_flag=1
              
          else
            
            continue
            
          endif
          
          call integral(pxxvvc_s,z0_mn,z0,pxxc)
          call integral(pzzvvc_s,z0_mn,z0,pzzc)
          
          pxxc=pxxc+pxxc_mn
          pzzc=pzzc+pzzc_mn
          
          call phihhb_g(z0, phi_g)
             
                                  
          stalk_phig_scat=coef*
     &         conjg(phi_g)*
     &         (-cos(theta_i)**2*pxxc+sin(theta_i)**2*pzzc)
          
        else if (mech_c .eq. pg) then
          if(z0_flag .eq. 0) then
            
            call integral(pxxvvc_s,-d,z0_mn,pxxc_mn)
            call integral(pzzvvc_s,-d,z0_mn,pzzc_mn)
            
            z0_flag=1
            
          else
            
            continue
            
          endif
          
          call integral(pxxvvc_s,z0_mn,z0,pxxc)
          call integral(pzzvvc_s,z0_mn,z0,pzzc)
          
          pxxc=pxxc+pxxc_mn
          pzzc=pzzc+pzzc_mn
          
          call phihhc_g(z0, phi_g)
             
                                  
            stalk_phig_scat=coef*
     &         conjg(phi_g)*
     &         (-cos(theta_i)**2*pxxc+sin(theta_i)**2*pzzc)
            
          else if (mech_c .eq. gpg) then
            if(z0_flag .eq. 0) then
              
              call integral(pxxvvc_s,-d,z0_mn,pxxc_mn)
              call integral(pzzvvc_s,-d,z0_mn,pzzc_mn)
              
              z0_flag=1
              
            else
              
              continue
              
            endif
            
              call integral(pxxvvc_s,z0_mn,z0,pxxc)
              call integral(pzzvvc_s,z0_mn,z0,pzzc)
              
              pxxc=pxxc+pxxc_mn
              pzzc=pzzc+pzzc_mn
              
              call phihhd_g(z0, phi_g)
             
                                      
              stalk_phig_scat=coef*
     &             conjg(phi_g)*
     &             (-cos(theta_i)**2*pxxc+sin(theta_i)**2*pzzc)
              
            else
              
              write(error,*)'mech_c = ',mech_c,' not recognized',
     &             ' in func stalk_phig_scat'
              code=49
              call the_death_code(die,error,code)
              
              
            endif
            
          else if (mech .eq. gpg) then
            
            if (mech_c .eq. p) then
              
              if(z0_flag .eq. 0) then
                
                call integral(pxxvvd_s,-d,z0_mn,pxxd_mn)
                call integral(pzzvvd_s,-d,z0_mn,pzzd_mn)
                call integral(pxzvvd_s,-d,z0_mn,pxzd_mn)
                
                z0_flag=1
                
              else
                
                continue
                
              endif
              
              call integral(pxxvvd_s,z0_mn,z0,pxxd)
              call integral(pzzvvd_s,z0_mn,z0,pzzd)
              call integral(pxzvvd_s,z0_mn,z0,pxzd)
              
              pxxd=pxxd+pxxd_mn
              pzzd=pzzd+pzzd_mn
              pxzd=pxzd+pxzd_mn
              
              call phihha_g(z0, phi_g)
              
                                      
              stalk_phig_scat=coef*
     &             conjg(phi_g)*
     &             (cos(theta_i)**2*pxxd
     &             +2.0*cos(theta_i)*sin(theta_i)*pxzd
     &             +sin(theta_i)**2*pzzd)
              
              
            else if (mech_c .eq. gp) then
              
              if(z0_flag .eq. 0) then
                
                call integral(pxxvvd_s,-d,z0_mn,pxxd_mn)
                call integral(pzzvvd_s,-d,z0_mn,pzzd_mn)
                call integral(pxzvvd_s,-d,z0_mn,pxzd_mn)
                
                z0_flag=1
                
              else
                
                continue
                
              endif
              
              call integral(pxxvvd_s,z0_mn,z0,pxxd)
              call integral(pzzvvd_s,z0_mn,z0,pzzd)
              call integral(pxzvvd_s,z0_mn,z0,pxzd)
              
              pxxd=pxxd+pxxd_mn
              pzzd=pzzd+pzzd_mn
              pxzd=pxzd+pxzd_mn
              
              call phihhb_g(z0, phi_g)
              
                                     
              stalk_phig_scat=coef*
     &             conjg(phi_g)*
     &             (cos(theta_i)**2*pxxd
     &             +2.0*cos(theta_i)*sin(theta_i)*pxzd
     &             +sin(theta_i)**2*pzzd)
              
            else if (mech_c .eq. pg) then
              
              if(z0_flag .eq. 0) then
                
                call integral(pxxvvd_s,-d,z0_mn,pxxd_mn)
                call integral(pzzvvd_s,-d,z0_mn,pzzd_mn)
                call integral(pxzvvd_s,-d,z0_mn,pxzd_mn)
                
                z0_flag=1
                
              else
                
                continue
                
              endif
              
              call integral(pxxvvd_s,z0_mn,z0,pxxd)
              call integral(pzzvvd_s,z0_mn,z0,pzzd)
              call integral(pxzvvd_s,z0_mn,z0,pxzd)
              
              pxxd=pxxd+pxxd_mn
              pzzd=pzzd+pzzd_mn
              pxzd=pxzd+pxzd_mn
              
              call phihhc_g(z0, phi_g)
             
                                      
              stalk_phig_scat=coef*
     &             conjg(phi_g)*
     &             (cos(theta_i)**2*pxxd
     &             +2.0*cos(theta_i)*sin(theta_i)*pxzd
     &             +sin(theta_i)**2*pzzd)
              
            else if (mech_c .eq. gpg) then
              
              if(z0_flag .eq. 0) then
                
                call integral(pxxvvd_s,-d,z0_mn,pxxd_mn)
                call integral(pzzvvd_s,-d,z0_mn,pzzd_mn)
                call integral(pxzvvd_s,-d,z0_mn,pxzd_mn)
                
                z0_flag=1
                
              else
                
                continue
                
              endif
              
              call integral(pxxvvd_s,z0_mn,z0,pxxd)
              call integral(pzzvvd_s,z0_mn,z0,pzzd)
              call integral(pxzvvd_s,z0_mn,z0,pxzd)
              
              pxxd=pxxd+pxxd_mn
              pzzd=pzzd+pzzd_mn
              pxzd=pxzd+pxzd_mn
              
              call phihhd_g(z0, phi_g)
             
                                      
              stalk_phig_scat=coef*
     &             conjg(phi_g)*
     &             (cos(theta_i)**2*pxxd
     &             +2.0*cos(theta_i)*sin(theta_i)*pxzd
     &             +sin(theta_i)**2*pzzd)
              
            else
              
              write(error,*)'mech_c = ',mech_c,' not recognized',
     &             ' in func stalk_phig_scat'
              code=49
              call the_death_code(die,error,code)
              
              
            endif
            
          else
            
            write(error,*)'mech = ',mech,' not recognized',
     &           ' in func stalk_phig_scat'
            code=50
            call the_death_code(die,error,code)
            
            
          endif
          
        else
          
          write(error,*)'pol = ',pol,' not recognized',
     &         ' in func stalk_phig_scat'
          code=51
          call the_death_code(die,error,code)
          
          
        endif
        
        return
        end

      
