z=NewtonRingsI(params,r)

%---------------------------------------------------------------
% PROGRAM: NRI
% AUTHOR:  Andri M. Gretarsson
% DATE:    8/5/03
%
%
% SYNTAX: z=NwetonRingsI([w1,w2,R1,R2,P1,P2,lambda,deltaphi],r);
%
% Returns the intensity of a set of Newton's rings as a function 
% of the radial coordinate R. (See e.g. Fundamentals of Physical Optics 
% by Jenkins and White, 1st ed. Section 4.4).  The Newton's rings are
% generated by the interference of two coaxial beams with different radii
% of curvature and/or width.
%
% w1       = 1/e field amplitude radius of the 1st of the interfering beams
% w2       = 1/e field amplitude radius of the 2nd of the interfering beams
% R1       = Radius of curvature of the phasefront of the 1st of the interfering beams
% R2       = Radius of curvature of the phasefront of the 2nd of the interfering beams
% P1       = Radius of curvature of the phasefront of the 2nd of the interfering beams
% P2       = Radius of curvature of the phasefront of the 2nd of the interfering beams
% lambda   = Radius of curvature of the phasefront of the 2nd of the interfering beams
% deltaphi = Radius of curvature of the phasefront of the 2nd of the interfering beams
% z(k)     = Intensity of the interference pattern at radiusr(k).
%
% If x and y are not vectors but matrixes generated by
% the matlab function meshgrid, then the output variable z
% is a matrix rather than a vector.  The matrix form allows
% the function to have a plane as it's domain rather than
% a curve. To generate a plane domain one uses meshgrid.
%
% EXAMPLE:  thetaoffset=20*pi/180;
%           rseed=[0:0.01:sqrt(0.3)].^2;
%           thetaseed=[0:1:360]*pi/180;
%           [theta,r]=meshgrid(thetaseed,rseed);
%           [x,y]=pol2cart(theta,r);
%           h=pcolor(x,y,NewtonRingsI([0.01,100,2e3,1e99,1,1,1.024e-6,thetaoffset],r));
%           colormap('bone');
%           set(h,'EdgeColor','none');
%           set(h,'FaceColor','interp');
%           set(gca,'Visible','off');
%           set(gcf,'Color','black');
%           axis square
%           shg;
%
% Last updated: 8/5/03 by AMG
%
%---------------------------------------------------------------
%% SYNTAX: z=NewtonRingsI([w1,w2,R1,R2,P1,P2,lambda,deltaphi],r);
%---------------------------------------------------------------

function z=NewtonRingsI(params,r)

w1=params(1);
w2=params(2);
R1=params(3);
R2=params(4);
P1=params(5);
P2=params(6);
lambda=params(7);
deltaphi=params(8);
a=1/2/R1-1/2/R2;

z = exp(-r.^2/w1).*exp(-r.^2/w2).*...
    ( sin( (a/2/lambda)* r.^2 - deltaphi) ).^2 + abs(P1-P2);