%
%---------------------------------------------------------------------  
% function [gp] = FT_init_gauss_param(lambda,nr,varargin)
% 
% Matlab function that creates a structure `gp' storing the
% Gaussian beam parameter 'q' and possible derived parameters.
%
% Usage for example:
%  gp=FT_init_gauss_param(lambda,nr,'q',q)
% or
%  gp=FT_init_gauss_param(lambda,nr,'w0',w0,'z',z)
%
% lambda (real): wavelength of laser light [m]
% nr (real):     index of refraction
%
% The function can compute the beam parameter `q' from a variety of
% input arguments, possible options are:
% q  or  q0, z  or  w0, z  or  w, z  or  zr, z  or  Rc, z  or  Rc, w
%
% If any of the above combinations is provided (and the wavelength and
% an index of refraction) the Gaussian beam paramter q is computed.
% In addition, the following derived paramters are computed and stored
% in the structure `gp':
% q0, w, w0, z, zr, Rc, Psi
%
% You can also specify an astigmatic mode by using combinations such as
% qx, qy  or  q0x, zx, q0y, zy   or  w0x, w0y, zx, zy   ...
%
% If an astigmatic beam is provided the variable gp.is_astigmatic is 
% set to 1 (otherwise 0).
%
% The structure returns:
% lambda (real): wavelength of laser light [m]
% nr (real):     index of refraction
% k (real):      wavenumber (as 2pi/lambda) [1/m] 
%
% (the following are provided as vectors with two elements where the first
% element refers to the value in the x-z plane, the second element for the y-z plane.
% A typical notation would be w0x=w0(1), w0y=w0(2), etc.
% q0 (complex): the Gaussian beam parameter at the waist
% w (real):     beam radius (x-z,y-z plane) [m]
% w0 (real):    waist radius (x-z,y-z plane) [m]
% z (real):     position on z-axis (distance from beam waist at z=0) [m]                    
%               negative values: running towards the waist
%               positive values: running away from the waist
% zr (real):    Rayleigh range [m]
% Rc (real):    radius of curvature of phase front [m]
%               negative values mean, centre of curvature lies in the direction of propagation
% Psi (real):   Gouy phase [rad]
%
% Part of the SimTools package
% Andreas Freise, 19.12.2009 afreise@googlemail.com
%---------------------------------------------------------------------  
%

% Description: Creates a structure "gp" storing the Gaussian beam parameter "q"
% Keywords: Gaussian, parameter, q, gp, create, init

 function [gp] = FT_init_gauss_param(lambda,nr,varargin)

   baseid='init_gauss_param';

   %%% initialise `gp' structure
   gp.lambda=lambda;
   gp.k=2*pi/lambda;
   gp.nr=nr;
   gp.is_astigmatic=0;
   
   gp.q=complex(zeros(2,1,'double'));
   gp.q0=complex(zeros(2,1,'double'));
   gp.w=zeros(2,1,'double');
   gp.w0=zeros(2,1,'double');
   gp.zr=zeros(2,1,'double');
   gp.z=zeros(2,1,'double');
   gp.Rc=zeros(2,1,'double');
   gp.Psi=zeros(2,1,'double');
   
   
   q_set=0;
   q0_set=0;
   w0_set=0;
   w_set=0;
   zr_set=0;
   z_set=0;
   Rc_set=0;
   
   %%% parse input arguments
   if nargin > 1
     % check which calling syntax is used
     actual_n=1;
     while actual_n<=nargin-2
       switch lower(varargin{actual_n})
        case 'q0',
         gp.q0(1:2) = varargin{actual_n+1};
         q0_set=bitor(q0_set,1);
        case 'q0x',
         gp.q0(1) = varargin{actual_n+1};
         q0_set=bitor(q0_set,2);
         gp.is_astigmatic=1;
        case 'q0y',
         gp.q0(2) = varargin{actual_n+1};
         q0_set=bitor(q0_set,4);
         gp.is_astigmatic=1;
        case 'q',
         gp.q(1:2) = varargin{actual_n+1};
         q_set=bitor(q_set,1);
        case 'qx',
         gp.q(1) = varargin{actual_n+1};
         q_set=bitor(q_set,2);
         gp.is_astigmatic=1;
        case 'qy',
         gp.q(2) = varargin{actual_n+1};
         q_set=bitor(q_set,4);
         gp.is_astigmatic=1;
        case 'w0',
         gp.w0(1:2) = varargin{actual_n+1};
         w0_set=bitor(w0_set,1);
        case 'w0x',
         gp.w0(1) = varargin{actual_n+1};
         w0_set=bitor(w0_set,2);
         gp.is_astigmatic=1;
        case 'w0y',
         gp.w0(2) = varargin{actual_n+1};
         w0_set=bitor(w0_set,4);
         gp.is_astigmatic=1;
        case 'w',
         gp.w(1:2) = varargin{actual_n+1};
         w_set=bitor(w_set,1);
        case 'wx',
         gp.w(1) = varargin{actual_n+1};
         w_set=bitor(w_set,2);
         gp.is_astigmatic=1;
        case 'wy',
         gp.w(2) = varargin{actual_n+1};
         w_set=bitor(w_set,4);
         gp.is_astigmatic=1;
        case 'zr',
         gp.zr(1:2) = varargin{actual_n+1};
         zr_set=bitor(zr_set,1);
        case 'zrx',
         gp.zr(1) = varargin{actual_n+1};
         zr_set=bitor(zr_set,2);
         gp.is_astigmatic=1;
        case 'zry',
         gp.zr(2) = varargin{actual_n+1};
         zr_set=bitor(zr_set,4);
         gp.is_astigmatic=1;
        case 'z',
         gp.z(1:2) = varargin{actual_n+1};
         z_set=bitor(z_set,1);
        case 'zx',
         gp.z(1) = varargin{actual_n+1};
         z_set=bitor(z_set,2);
         gp.is_astigmatic=1;
        case 'zy',
         gp.z(2) = varargin{actual_n+1};
         z_set=bitor(z_set,4);
         gp.is_astigmatic=1;
        case 'rc',
         gp.Rc(1:2) = varargin{actual_n+1};
         Rc_set=bitor(Rc_set,1);
        case 'rcx',
         gp.Rc(1) = varargin{actual_n+1};
         Rc_set=bitor(Rc_set,2);
         gp.is_astigmatic=1;
        case 'rcy',
         gp.Rc(2) = varargin{actual_n+1};
         Rc_set=bitor(Rc_set,4);
         gp.is_astigmatic=1;
        otherwise ,
          msgid=[baseid,':checkinarg'];
          result=sprintf(['Unknown option: ', varargin{actual_n}, '.\n']);
          error(msgid,result);
         return
       end;   
       actual_n=actual_n+2;
     end;  
   end
   
   % check for consistency of arguments
   
   any_set=bitor(bitor(q0_set,q_set),bitor(bitor(w0_set,w_set),bitor(z_set,Rc_set)));
   
   if (bitand(any_set,1) && gp.is_astigmatic==1 )
       msgid=[baseid,':astigmatism'];
       result=sprintf('Input parameters inconsistent: do not use q,q0,w,w0,z or Rc for an astigmatic beam.');
       error(msgid,result);
   end
   
   % check possible parameter combinations:
   % - q
   % - q0, z
   % - w0, z
   % - w, z
   % - zr, z
   % - Rc, z
   % - w, Rc
   
   complete_params=0;
   
   if (gp.is_astigmatic)
     for_idx=[2,4];
   else
     for_idx=[1];
   end
   
   for l=for_idx
     z_required=(bitand(q0_set,l)+bitand(w0_set,l)+bitand(zr_set,l)+bitand(w_set,l)+bitand(Rc_set,l))/l;
     if (bitand(q_set,l))
       if (z_required>0 || bitand(z_set,l))
         msgid=[baseid,':input_params'];
         result=sprintf('Input parameters inconsistent: "q" fully characterises the beam.');
         error(msgid,result);
       else
         complete_params=1;
       end
     elseif (z_required>1)
       if ((z_required~=2) || ((bitand(w_set,l)+bitand(Rc_set,l))/l~=2))
       msgid=[baseid,':input_params'];
       result=sprintf('Input parameters inconsistent: multiple parameters given which (indirectly) define the beam size.');
       error(msgid,result);     
       else
         complete_params=1;
       end                
     else
       if (bitand(z_set,l))
         complete_params=1;
       else
         msgid=[baseid,':input_params'];
         switch l
           case 1,
            result=sprintf('Input parameters inconsistent: "z" required but not given.');
           case 2,
            result=sprintf('Input parameters inconsistent: "zx" required but not given.');
           case 4,
            result=sprintf('Input parameters inconsistent: "zy" required but not given.');
         end
         error(msgid,result);     
       end
     end
   end
    
   if (complete_params~=1)
     msgid=[baseid,':input_params'];
     result=sprintf('Input parameters inconsistent.');
     error(msgid,result);     
   end
   
   % converting the following parameters combinations into q   
   % - q0, z
   % - w0, z
   % - w, z
   % - zr, z
   % - Rc, z
   % - w, Rc
   
   for l=[1,2]
     if (gp.is_astigmatic)
       idx=2*l;
     else
       idx=1;
     end
     if (bitand(q0_set,idx))
       gp.q(l)=gp.q0(l)+gp.z(l);
     end
     if (bitand(w0_set,idx))
       zr=pi./(gp.lambda./gp.nr).*gp.w0(l).^2;
       gp.q(l)=1i*zr+gp.z(l);
     end
     if (bitand(w_set,idx))
       tmp=pi/2./(gp.lambda./gp.nr)*gp.w(l).^2;
       zr=tmp+sqrt(tmp.^2-gp.z(l).^2);
       gp.q(l)=1i*zr+gp.z(l);
     end
     if (bitand(zr_set,idx))
       gp.q(l)=1i*gp.zr(l)+gp.z(l);
     end
     if (bitand(Rc_set,idx))
       if (bitand(z_set,idx))
         zr=sqrt(gp.z(l)*(gp.Rc(l)-gp.z(l)));
         gp.q(l)=1i*zr+gp.z(l);
       else
         [w0,z]=FT_Rcw_to_w0z(gp.Rc(l),gp.w(l),gp.lambda,gp.nr);
         zr=pi./(gp.lambda./gp.nr)*w0.^2;
         gp.q(l)=1i*zr+z;
       end
     end
        
   end
    
   % Now we know q and can (re-) compute the rest:
   gp=FT_update_gauss_param(gp,gp.q,gp.nr);
 return
