%
% --------------------------------------------------------------------
% function field = FT_ulp_cos (lambda,w0, r, phi, z, p, l)
% (required LaguerrePoly.m)
%
% Octave/Matlab function to compute the field of a sinusoidal 2D 
% Laguerre-Gauss function u_lp(x,y,z), i.e. uses cos(m phi) instead of 
% exp( i m phi) as in FT_ulp.m, see
% http://relativity.livingreviews.org/Articles/lrr-2010-1/
% 
% lambda = wavelength
% w0 (real):  beam widths [m]
% r (real):   distance to optical axis [m]
% phi (real): position around optical axis [rad]
% z (real):   distance to waist [m]
% l (int):    order of mode (azimuthal index, p>=|l|>=0)
% p (int):    order of mode (radial index, p>=0)
%
% field (complex): field amplitude at (x,y,z)
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
% Andreas Freise 25.03.2007
% --------------------------------------------------------------------
%

% Description: Computes the field amplitude of a Laguerre-Gauss mode (sinusoidal)
% Keywords: Laguerre, mode, Gauss, field, amplitude

function field = FT_ulp_cos(lambda,w0,r,phi,z,p,l)

  sl=l;
  l=abs(l);
  
  sigl=sign(sl+0.1);
  
  field=0;
   
  k=2*pi/lambda;
  zr=pi*w0^2/lambda;
  wz=w0*sqrt(1+(z/zr)^2);
  sr=sqrt(2)*r/wz;
  qz=z+i*zr;
  psi=atan(z/zr);
  
  t1=sqrt(4*factorial(p)/((1+FT_kd(0,l))*pi*factorial(l+p)))/wz;
  t2=exp(i *(2*p+l+1)*psi);
  t3=sr.^l.*FT_LaguerrePol(p,l,sr.^2);
  t4=exp(-i * k*r.^2./(2*qz)).*((1+sigl)/2.*cos(sl*phi)+(1-sigl)/2.*sin(sl*phi));
  
  field=t1*t2*t3.*t4;

  return