function [ Itot, M_I_rho, M_A_rho ] = fresnel_occultation_fs(rhoStar, rhoOccultor, vectorOfSamplePoints, impactParameter, spatialRes)
% Monochromatic diffraction pattern by finite source
% Package: AstroUtil.Occultation
% Description: Calculate the diffraction pattern generated by a finite,
%              non-point, monochromatic source.
% Input  : - Radius of the source (e.g., star) in units of the fresnel
%            scale (radius / F), where F=sqrt(\lambda d/2).
%            Default is 0.5.
%          - Radius of the occulter in units of the fresnel scale
%            (radius / F). Default is 0.5.
%          - A one dimensional row vector specifies the x-coordinates of 
%            the occulter trajectory in which to calculate the diffraction.
%            [Fresnel radius units]. Default is (-10:0.1:10).
%          - The impact parameter in units of the Fresnel scale radius.
%            I.e., The minimum distance from the source center to the
%            occulter trajectory. Default is 0.
%          - The integration grid spatial resolution, in units of the
%            stellar radius. Default is 0.1.
% Output : - Row vector of the total intensity, cooresponding to each
%            position along the occulter trajectory.
%          - Matrix of total intensity. Rows correspond to diiferent
%            point along the trajectory while columns are for the
%            various resolution elements of the source.
%          - Amplitude and phase (complex) matrix.
%            Rows correspond to diiferent
%            point along the trajectory while columns are for the
%            various resolution elements of the source.
% Tested : Matlab R2011b
%     By : Lior Rubanenko                  Jun 2009
%    URL : http://weizmann.ac.il/home/eofek/matlab/
% See also: fresnel_occultation_ps.m
% Examples: [sum, M_I, M_A] = AstroUtil.Occultation.fresnel_occultation_fs; plot((-10:0.1:10).',sum)
% Reliable: 2
%--------------------------------------------------------------------------

import AstroUtil.Occultation.*

Def.rhoStar               = 0.5;
Def.rhoOccultor           = 0.5;
Def.vectorOfSamplePoints  = (-10:0.1:10);
Def.impactParameter       = 0;
Def.spatialRes            = 0.1;

if (nargin==0),
    rhoStar               = Def.rhoStar;
    rhoOccultor           = Def.rhoOccultor;
    vectorOfSamplePoints  = Def.vectorOfSamplePoints;
    impactParameter       = Def.impactParameter;
    spatialRes            = Def.spatialRes;
elseif (nargin==1),
    rhoOccultor           = Def.rhoOccultor;
    vectorOfSamplePoints  = Def.vectorOfSamplePoints;
    impactParameter       = Def.impactParameter;
    spatialRes            = Def.spatialRes;
elseif (nargin==2),
    vectorOfSamplePoints  = Def.vectorOfSamplePoints;
    impactParameter       = Def.impactParameter;
    spatialRes            = Def.spatialRes;
elseif (nargin==3),
    impactParameter       = Def.impactParameter;
    spatialRes            = Def.spatialRes;
elseif (nargin==4),
    spatialRes            = Def.spatialRes;
elseif (nargin==5),
    % do nothing
else
    error('Illegal number of input arguments');
end

% Field and KBO trajectory calculation:
starCenterCoordinates = [0 0];
x = -rhoStar:spatialRes:rhoStar;
y = rhoStar:-spatialRes:-rhoStar;

[X,Y] = meshgrid(x,y);

field = zeros(size(X));
starIndices = (X - starCenterCoordinates(1)).^2 + (Y - starCenterCoordinates(2)).^2 <= rhoStar.^2;
%   Put ones where the star is located:
field(starIndices) = 1;  

%%  r(t) calculation:
starCoordinatesIndices = find(field);

%   Positions of occulter:
occultorCoordinatesAtTimeT = [vectorOfSamplePoints; ones(size(vectorOfSamplePoints)) * impactParameter]; 

deltaX = bsxfun(@minus, X(starCoordinatesIndices), occultorCoordinatesAtTimeT(1,:));
deltaY = bsxfun(@minus, Y(starCoordinatesIndices), occultorCoordinatesAtTimeT(2,:));

D = sqrt(deltaX.^2 + deltaY.^2);

%   Calculating the diffraction pattern:
[M_I_rho, M_A_rho] = fresnel_occultation_ps(rhoOccultor, D);

s = size(M_I_rho);
Itot = sum(M_I_rho) / s(1);

