c     Authors: P. C. Chaumet 
c     Date: 24/11/2009

c     Purpose: These routines compute a Gaussian beam in the case of
c     highly nonparaxial field ($w_0<\lambda$). The computation is done
c     in the vectorial case for a polarized Gaussian beam and is valid
c     whatever the position: far away or close to/in the waist. It is
c     also computed from the the power of the Gaussian beam the
c     intensity and the irradiance of the field at the origin. Notice
c     that the origin of the beam is taken at the center of the
c     waist. The profile of the electric field at the origin is
c     E=E0e^{-rho^2/2/w_0^2}
c     
c     Purpose (in french): Cet ensemble de routines calcul le champ
c     electrique d'un faisceau Gaussian dans le cas general d'un
c     faisceau pouvant etre tres fortement non paraxial
c     ($w_0<\lambda$). Le calcul est fait dans le cas vectoriel a trois
c     dimensions pour un faisceau polarise.  Le calcul du faisceau
c     gaussian est base sur un developpement du champ en spectre d'onde
c     plane. La position d'observation peut etre a priori n'importe ou
c     loin comme proche/dans le waist; A noter que le programme n'a pas
c     ete teste en champ tres lointain ou dans ce cas des approximations
c     analytiques se trouvent dans la litterature. Il existe aussi une
c     routine qui sachant la puissance en watt du faisceau Gaussian
c     donne l'intensite du champ a l'origine ainsi que l'irradiance a
c     cette meme position. A noter que quand je parle d'origine il
c     s'agit toujours du centre du waist du faisceau Gaussian. Le profil
c     du champ electrique dans le waist s'ecrit: E=E0e^{-rho^2/2/w_0^2}


c     Reference: if you use this routine in your research, please
c     reference, as appropriate: P. C. Chaumet, Fully vectorial highly
c     non paraxial beam close to the waist. J. Opt. Soc. Am. A 23, 3197
c     (2006) and for the computation of the irradiance P. C. Chaumet,
c     B. Pouligny, R. Dimova and N. Sojic, Optical tweezers in
c     interaction with an apertureless probe J. Appl. Phys. 102, 024915
c     (2007).

c     license: GNU GPL

c     Notice that the routines to compute the Gaussian beam use needs
c     some routines of netlib. For convenience of the users I have
c     enclosed theses routines in one file: routine.f

c     I cannot guarantee correctness of the programs...

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

c     Computation of the intensity and irradiance of the beam at the
c     origin from the power of the Gaussian beam.

c     P: power (Watt)

c     irra: Irradiance (Watt/m^2)

c     w0: waist (m)

c     k0: wavenumber (m^{-1})

c     Igauss: Intensity of the field. Then the magnitude of the elctric
c     field is E0=dsqrt(Igauss) (V/m)

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

      subroutine gaussianpuissance(P,irra,w0,k0,Igauss)
      implicit none
      double precision P,w0,k0,Igauss,c,pi,u,v,const,const0,const1,eps0
      double complex icomp,zarg,uncomp,zf
      logical flag
      double precision Ip, Ipower,ww0,irra

c     declarationlie a netlib
      integer limit,lenw,last,iwork,neval,ier,key
      parameter (limit=100,lenw=limit*5)
      dimension  iwork(limit)
      double precision epsabs,epsrel,abserr,work,rap,borneinf,bornesup
      dimension work(lenw)
      external Ipower
      common/comgaussian0/ww0
      ww0=w0

      c=299792458.d0
      pi=dacos(-1.d0)
      eps0=1.d0/(c*c*4.d0*pi*1.d-7)
      icomp=(0.d0,1.d0)
      uncomp=(1.d0,0.d0)
      
      const0=pi/4.d0*w0*w0*c*eps0
      const1=k0*w0

      zarg=const1*uncomp
      call WOFZGEN (zarg, U, V, FLAG)
      zf=(u*uncomp+v*icomp)
      if (dimag(zf).ne.0.d0) write(*,*) 'puissance imag',zarg,zf
      
      const=const0*(1.d0+(const1*const1-1)/const1*dsqrt(pi)/2.d0
     $     *dimag(zf))
      Igauss=P/const

c     Compute irradiance
      
      zarg=const1*uncomp/2.d0
      call WOFZGEN (zarg, U, V, FLAG)
      zf=(u*uncomp+v*icomp)
      if (dimag(zf).ne.0.d0) write(*,*) 'irradiance imag',zarg,zf
      irra=c*eps0*(1.d0+(const1*const1-1)/const1*dsqrt(pi/2.d0)*dimag(zf
     $     ))*Igauss

      end

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

c     subroutine to call to compute the Gaussian beam.
c     The origin is defined at the center of the waist

c     w0: waist (m)

c     k0: wavenumber (m^{-1})

c     E0x,E0y,E0z, magnitude of the field at the origin, one can uses
c     the routine gaussianpuissance to know this field if only the power
c     of the Gaussian beam is known. (V/m)

c     x,y,z location where the electric field is wanted (m).

c     x0,y0,z0 location of the center of the waist (origin) (m).

c     Ex,Ey,Ez electric field at the position x,y,z (V/m)

c     vu(3) unit vector which gives the direction of propagation of the
c     Gaussian beam. It should be perpendicular to the field at the
c     origin.

c     tol: tolerance used when the integration is performed

c     nloin: in the angular spectrum representation one can uses only
c     the propagating waves (0), the evanescent waves (1) or both (2).

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

      subroutine gaussianchamp(x,y,z,x0,y0,z0,vu,w0,k0,E0x,E0y,E0z,Ex,Ey
     $     ,Ez,tol,nloin)
      implicit none
      integer nloin,ii
      double precision x,y,z,xx,yy,zz,x0,y0,z0,a,k0,w0,sinthe,costhe
     $     ,sin2the,cos2the,tol,rpos(3),vu(3),Enorm,theta,phi,vn,eps
     $     ,mat(3,3)
      double complex E0x,E0y,E0z,E00x,E00y,E00z,Ex,Ey,Ez,Epx,Epy,Epz
     $     ,Ixxp,Ixxe,Izzp,Izze,Ixx,Izz,icomp,ctmp
      
c     precision machine
      eps=1.d-12

c     check if the unit vector is perpendicular to the incident field.
      Enorm=dsqrt(cdabs(E0x)*cdabs(E0x)+cdabs(E0y)*cdabs(E0y)+cdabs(E0z)
     $     *cdabs(E0z))
      ctmp=E0x*vu(1)+E0y*vu(2)+E0z*vu(3)
      if (cdabs(ctmp)/Enorm.ge.eps) then
         write(*,*)
     $        'unit vector not perpendicular to the incident field'
         stop
      endif
c     check if vv is a unit vector
      vn=vu(1)*vu(1)+vu(2)*vu(2)+vu(3)*vu(3)
      if (vn.le.1.d0-eps.or.vn.ge.1.d0+eps) then
         write(*,*) 'vu is not a unit vector'
         stop
      endif

c     translation
      rpos(1)=x-x0
      rpos(2)=y-y0
      rpos(3)=z-z0
c     compute the matrix of rotation from the unit vector
      vn=vu(1)*vu(1)+vu(2)*vu(2)
      if (vn.le.eps) then
         mat(1,1)=1.d0
         mat(1,2)=0.d0
         mat(1,3)=0.d0
         mat(2,1)=0.d0
         mat(2,2)=1.d0
         mat(2,3)=0.d0
         mat(3,1)=0.d0
         mat(3,2)=0.d0
         mat(3,3)=1.d0
      else
         mat(1,1)=vu(3)+(1.d0-vu(3))*vu(2)*vu(2)/vn
         mat(1,2)=-(1.d0-vu(3))*vu(1)*vu(2)/vn
         mat(1,3)=-vu(1)
         mat(2,1)=mat(1,2)
         mat(2,2)=vu(3)+(1.d0-vu(3))*vu(1)*vu(1)/vn
         mat(2,3)=-vu(2)
         mat(3,1)=vu(1)
         mat(3,2)=vu(2)
         mat(3,3)=vu(3)
      endif

c     compute the new incident field and position
      xx=mat(1,1)*rpos(1)+mat(1,2)*rpos(2)+mat(1,3)*rpos(3)
      yy=mat(2,1)*rpos(1)+mat(2,2)*rpos(2)+mat(2,3)*rpos(3)
      zz=mat(3,1)*rpos(1)+mat(3,2)*rpos(2)+mat(3,3)*rpos(3)
      
c      write(*,*) 'E',E0x,E0y,E0z

      E00x=mat(1,1)*E0x+mat(1,2)*E0y+mat(1,3)*E0z
      E00y=mat(2,1)*E0x+mat(2,2)*E0y+mat(2,3)*E0z
      E00z=mat(3,1)*E0x+mat(3,2)*E0y+mat(3,3)*E0z

c      write(*,*) 'E',E00x,E00y,E00z

      icomp=(0.d0,1.d0)
      a=dsqrt(xx*xx+yy*yy)
      call gaussian(a,zz,w0,k0,Ixxp,Ixxe,Izzp,Izze,tol,nloin)
      
      Ixx=Ixxp+Ixxe
      Izz=Izzp+Izze

      sinthe=xx/a
      costhe=yy/a
      if (a.eq.0.d0) then
         sinthe=0.d0
         costhe=0.d0
      endif
      sin2the=2.d0*sinthe*costhe
      cos2the=costhe*costhe-sinthe*sinthe
      
      Epx=Ixx*E00x
      Epy=Ixx*E00y
      Epz=-icomp*Izz*(E00x*sinthe+E00y*costhe)

c     compute the final field with the inverse rotation matrix
      Ex=mat(1,1)*Epx+mat(2,1)*Epy+mat(3,1)*Epz
      Ey=mat(1,2)*Epx+mat(2,2)*Epy+mat(3,2)*Epz
      Ez=mat(1,3)*Epx+mat(2,3)*Epy+mat(3,3)*Epz

      end
c*****************************************************************
c*****************************************************************
c*****************************************************************
      subroutine gaussian(a,zz,w0,k0,Ixxp,Ixxe,Izzp,Izze,tol,nloin)
      implicit none
      integer nloin
      double precision a,z,w0,k0,zz,tol
      double complex Ixxp,Ixxe,Izzp,Izze
      double precision expe,spi,w0k0,zp,u,v
      double complex expz,icomp,uncomp,zzz
      logical FLAG

      if (nloin.gt.2.or.nloin.lt.0) then
         write(*,*) 'The value of nloin is not correct',nloin
         stop
      endif

      if (zz.lt.0.d0) then
         z=-zz
      else
         z=zz
      endif

      if (a.eq.0.d0) then
         Izzp=0.d0
         Izze=0.d0

         uncomp=(1.d0,0.d0)
         icomp=(0.d0,1.d0)
         expz=cdexp(icomp*k0*z)
         expe=dexp(-k0*k0*w0*w0/2.d0)
         zp=z/w0/dsqrt(2.d0)
         spi=dsqrt(dacos(-1.d0))
         w0k0=w0*k0/dsqrt(2.d0)
         
         zzz=icomp*zp
         call WOFZ (zzz,u,v,FLAG)

         Ixxe=expe*(1.d0-spi*zp*(u*uncomp+v*icomp))
         zzz=icomp*zp+w0k0
         call WOFZ (zzz,u,v,FLAG)
         Ixxp=expz*(1.d0-spi*zp*(u*uncomp+v*icomp))-Ixxe

      else
         call  gaussianzz(a,z,k0,w0,tol,nloin,Izzp,Izze)      
         call  gaussianxx(a,z,k0,w0,tol,nloin,Ixxp,Ixxe)      
      endif
      
      if (zz.lt.0.d0) then
         Ixxp=dconjg(Ixxp)
         Izzp=dconjg(Izzp)
         Ixxe=dconjg(Ixxe)
         Izze=dconjg(Izze)
      endif


      end
c*****************************************************************
c*****************************************************************
c*****************************************************************
      subroutine gaussianxx(a,z,k0,w0,tol,nloin,Ixxp,Ixxe)
      implicit none
      integer nloin,test1,test2
      double precision a,z,k0,w0,aa,zz,kk0,ww0,tol
      double precision I1xxpr,I1xxpi,I1xxer,I1xxei
      double complex Ixxp,Ixxe,icomp

c     declarationlie a netlib
      integer limit,lenw,last,iwork,neval,ier,key
      parameter (limit=100,lenw=limit*5)
      dimension  iwork(limit),ier(4),neval(4)
      double precision epsabs,epsrel,abserr,work,rap,borneinf,bornesup
      dimension work(lenw),abserr(4),rap(4)

      external Ixxgaussian


      common/comgaussian/aa,zz,kk0,ww0
      common/testgaussian/test1,test2

      icomp=(0.d0,1.d0)

c     reinitialise a cause du common
      aa=a
      zz=z
      kk0=k0*k0
      ww0=w0*w0

      Ixxp=0.d0
      Ixxe=0.d0

      if (nloin.eq.0.or.nloin.eq.1) then
c     PREMIERE INTEGRALE DE 0 A K0
         test2=0
         borneinf=max(0.d0,kk0-100.d0/ww0)
         borneinf=dsqrt(borneinf)
         bornesup=k0
c     PARTIE REELLE
         key=6
         epsrel=tol
         epsabs=0.d0
         test1=0
         
         call dqag(Ixxgaussian,borneinf,bornesup,EPSABS,EPSREL,KEY
     $        ,I1xxpr,ABSERR(1),NEVAL(1),IER(1),LIMIT,LENW,LAST,IWORK
     $        ,WORK)
         
C     PARTIE IMAGINAIRE
         key=6
         epsrel=tol
         epsabs=0.d0
         test1=1
         
         call dqag(Ixxgaussian,borneinf,bornesup,EPSABS,EPSREL,KEY
     $        ,I1xxpi,ABSERR(2),NEVAL(2),IER(2),LIMIT,LENW,LAST,IWORK
     $        ,WORK)
         
C     TOTAL
         Ixxp=I1xxpr+icomp*I1xxpi
      endif
         
      if (nloin.eq.0.or.nloin.eq.2) then
c     DEUXIEME INTEGRALE DE 0 A L INFINI
         test2=1
         borneinf=0.d0
c     bornesup soit avec z soit avec w0
         bornesup=min(300.d0/(z+1.d-10),10.d0/w0)
c     PARTIE IMAGINAIRE
         I1xxei=0.d0
C     PARTIE REELLE
         key=6
         epsrel=tol
         epsabs=I1xxpr*epsrel
         test1=1
         
         call dqag(Ixxgaussian,borneinf,bornesup,EPSABS,EPSREL,KEY
     $        ,I1xxer,ABSERR(3),NEVAL(3),IER(3),LIMIT,LENW,LAST,IWORK
     $        ,WORK)        
         Ixxe=-icomp*I1xxei+I1xxer
      endif

      end
c*******************************************************************
c*******************************************************************
c*******************************************************************
      double precision function Ixxgaussian(kz)
      implicit none
      integer test1,test2
      double precision a,z,k0,w0,kz,kpara,cons,dbesj0
      double complex kzc

      common/comgaussian/a,z,k0,w0
      common/testgaussian/test1,test2

      if (test2.eq.0) kzc=kz*(1.d0,0.d0)
      if (test2.eq.1) kzc=kz*(0.d0,1.d0)

      kpara=dsqrt(k0-dreal(kzc*kzc))
      cons=w0*dexp(-kpara*kpara*w0/2.d0)*dbesj0(a*kpara)

      if(test1.eq.0)Ixxgaussian=cons*dreal(kzc*cdexp((0.d0,1.d0)*kzc*z))
      if(test1.eq.1)Ixxgaussian=cons*dimag(kzc*cdexp((0.d0,1.d0)*kzc*z))

      end
c*****************************************************************
c*****************************************************************
c*****************************************************************
      subroutine gaussianzz(a,z,k0,w0,tol,nloin,Izzp,Izze)
      implicit none
      integer nloin,test1,test2
      double precision a,z,k0,w0,aa,zz,kk0,ww0,tol
      double precision I1zzpr,I1zzpi,I1zzer,I1zzei
      double complex Izzp,Izze,icomp

c     declarationlie a netlib
      integer limit,lenw,last,iwork,neval,ier,key
      parameter (limit=100,lenw=limit*5)
      dimension  iwork(limit),ier(4),neval(4)
      double precision epsabs,epsrel,abserr,work,rap,borneinf,bornesup
      dimension work(lenw),abserr(4),rap(4)

      external Izzgaussian


      common/comgaussian/aa,zz,kk0,ww0
      common/testgaussian/test1,test2

      icomp=(0.d0,1.d0)

c     reinitialise a cause du common
      aa=a
      zz=z
      kk0=k0*k0
      ww0=w0*w0

      Izzp=0.d0
      Izze=0.d0

      if (nloin.eq.0.or.nloin.eq.1) then
c     PREMIERE INTEGRALE DE 0 A K0
         test2=0
         borneinf=max(0.d0,kk0-100.d0/ww0)
         borneinf=dsqrt(borneinf)
         bornesup=k0
c     PARTIE REELLE
         key=6
         epsrel=tol
         epsabs=0.d0
         test1=0
         
         call dqag(Izzgaussian,borneinf,bornesup,EPSABS,EPSREL,KEY
     $        ,I1zzpr,ABSERR(1),NEVAL(1),IER(1),LIMIT,LENW,LAST,IWORK
     $        ,WORK)

C     PARTIE IMAGINAIRE
         key=6
         epsrel=tol
         epsabs=0.d0
         test1=1
         
         call dqag(Izzgaussian,borneinf,bornesup,EPSABS,EPSREL,KEY
     $        ,I1zzpi,ABSERR(2),NEVAL(2),IER(2),LIMIT,LENW,LAST,IWORK
     $        ,WORK)
         
C     TOTAL
         Izzp=I1zzpr+icomp*I1zzpi
      endif
         
      if (nloin.eq.0.or.nloin.eq.2) then
c     DEUXIEME INTEGRALE DE 0 A L INFINI
         test2=1
         borneinf=0.d0
c     bornesup soit avec z soit avec w0
         bornesup=min(300.d0/(z+1.d-10),20.d0/w0)
c     PARTIE IMAGINAIRE
         key=6
         epsrel=tol
         epsabs=I1zzpi*epsrel
         test1=0   
         call dqag(Izzgaussian,borneinf,bornesup,EPSABS,EPSREL,KEY
     $        ,I1zzei,ABSERR(2),NEVAL(2),IER(2),LIMIT,LENW,LAST,IWORK
     $        ,WORK)        
C     PARTIE REELLE
         I1zzer=0.d0
         Izze=-icomp*I1zzei+I1zzer
      endif
      end
c*******************************************************************
c*******************************************************************
c*******************************************************************
      double precision function Izzgaussian(kz)
      implicit none
      integer test1,test2
      double precision a,z,k0,w0,kz,kpara,cons,dbesj1
      double complex kzc

      common/comgaussian/a,z,k0,w0
      common/testgaussian/test1,test2

      if (test2.eq.0) kzc=kz*(1.d0,0.d0)
      if (test2.eq.1) kzc=kz*(0.d0,1.d0)

      kpara=dsqrt(k0-dreal(kzc*kzc))
      cons=w0*dexp(-kpara*kpara*w0/2.d0)*dbesj1(a*kpara)

      if(test1.eq.0)Izzgaussian=cons*dreal(kpara*cdexp((0.d0,1.d0)*kzc*z
     $     ))
      if(test1.eq.1)Izzgaussian=cons*dimag(kpara*cdexp((0.d0,1.d0)*kzc*z
     $     ))

      end
