function [C, optics, appear] = votf(optics, appear)

% votfoct  Calculate the projected vectorial optical transfer function
% [C, optics, appear] = votfoct(optics, appear)
%
%   Inputs:
%       optics describes the physics of the situation being modeled
%       appear has the details of how to present the results
%
%   Outputs:
%       C = complex vectorial OTF = Cx + Cy + Cz
%       optics = physical paramaters 
%                (so you can see what the defaults were set to)
%       appear = parameters used for display plus extra things added during run
%
%   To run with the defaults, just use votf(0,0);
%
% MRA USyd June 2002
% $Id: votf.m,v 1.4 2003-06-24 11:11:54+10 matthewa Exp matthewa $

% Copyright (C) 2003 Matthew Arnison
% This program is copyleft under the GNU General Public License
% See http://www.gnu.org/copyleft/gpl.html for details

% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.

% todo: truncate size of output C to max bandwidth, as determined by the aperture

system('date;');
disp('$Id: votf.m,v 1.4 2003-06-24 11:11:54+10 matthewa Exp matthewa $')
% bypass this bit in MATLAB for now
if exist('OCTAVE_VERSION') == 102
	rcsdiff('$RCSfile: votf.m,v $');
	fflush(stdout);
end

% save input parameters
opticsin = optics;
appearin = appear;

% defaults
optics.alpha = pi/3; 			% aperture half angle
optics.Apv = 0;					% strength of cubic phase mask, peak to valley, in waves
optics.apodisation = 'sine';	% can be sine or herschel
optics.kmax = 4;				% maximum frequency in k space @ x0 +/- N/2
optics.lambda0 = 0.53;			% wavelength of light in vacuum
optics.n1 = 1.518;				% refractive index in focal region
optics.model = 'vectorial'; 	% accuracy: vectorial, scalar, or paraxial
optics.N = 512;					% accuracy: number of samples of pupil
optics.pupilshape = 'square';	% shape of the aperture: square or circle
optics.zdefocus = 0;			% microns of defocus
optics.z45 = 'clear';			% CJRS aperture filtering: clear, opaque, flip
								% see MRA logbook 13aug02 p. 117
								% DIC options, set dicshear to zero to
								% disable DIC
optics.dicshear = 0;			% DIC shear distance 2 delta x
optics.dicshearxy = 'x';		% DIC shear along x or y axis
optics.dicbias = 0;				% DIC bias 2 theta

appear.result = 'otf';			% what result in C: pupil, psf, ipsf, otf
appear.figstart = 1;			% where to start figure numbers for new plots
appear.logscale = 1;			% whether to use a log10 scale for plots
appear.realabs = 'abs';			% whether plots should show real() or abs()
appear.savemem = 1;				% whether to save memory by clearing temp matrices
appear.showpupil = 0;			% whether to display images of the pupil
appear.showplots = 1;			% whether to display plots
appear.showphase = 1;			% whether to display phase image of result
appear.saveabns = 0;			% whether to save the aberration functions 
								% (in optics.defocus & optics.cubic)
appear.savephase = 0;			% whether to save the phase of the pupil
								% in order to check Nyquist is satisfied
appear.debug = 0;				% debugging flag: set to 1 for extra output

% set values for this run

keys = fieldnames(opticsin);
sizekeys = size(keys); numkeys = sizekeys(1);

for keyi = 1:numkeys
	eval(['optics.' keys{keyi} ' = opticsin.' keys{keyi} ';']);
end

keys = fieldnames(appearin);
sizekeys = size(keys); numkeys = sizekeys(1);

for keyi = 1:numkeys
	eval(['appear.' keys{keyi} ' = appearin.' keys{keyi} ';']);
end

%%% end of config options

% store version of code used to do calculations inside optics output structure
% this gets overwritten each time the code is executed
optics.versionid = '$Id: votf.m,v 1.4 2003-06-24 11:11:54+10 matthewa Exp matthewa $';

% scaling factor to fit over a circular pupil 
% (or a square fitting inside that pupil)
% from logbook 13/2/2 p. 166
% not quite right, but within 5%, close enough for now
Acpv = sqrt(2)*optics.Apv/(2*sin(optics.alpha)^3);

% for max aperture NA=1, alpha = pi/2, radius is 1
% use aperture angle alpha to scale down from there
Rp = sin(optics.alpha);

% find centre
x0 = optics.N/2 + 1;
y0 = x0;

optics.lambda1 = optics.lambda0/optics.n1;
% so k here is actually k1, i.e. wave number in focal region
k = 2*pi/optics.lambda1;

% find width and edges of pupil
appear.pwidth = optics.N/optics.kmax*Rp;
appear.pledge = round(x0 - appear.pwidth/2);
appear.predge = round(x0 + appear.pwidth/2);
appear.prange = appear.pledge:appear.predge;

% find width and edges of OTF
sqwidth = optics.N/2*Rp/sqrt(2) - 2;
cwidth = optics.N/2*Rp - 2;
width = sqwidth;
ledge = round(x0 - width/2);
redge = round(x0 + width/2);
cledge = round(x0 - cwidth/2);
credge = round(x0 + cwidth/2);

if strcmp(appear.result, 'otf')
	kcircrange = linspace(-2*Rp, 2*Rp, credge - cledge + 1);
	kdiagrange = linspace(-2*Rp, 2*Rp, redge - ledge + 1);
	ksqrange = linspace(-2*Rp/sqrt(2), 2*Rp/sqrt(2), redge - ledge + 1);
elseif strcmp(appear.result, 'ipsf') || strcmp(appear.result, 'psf') ...
	|| strcmp(appear.result, 'pupil')
	
	% scaling for focal region
	% see logbook 23/9/02 p 135
	% physical units the same as lambda0

	% xmax is half the width of PSF as x=0 is in the centre
	% this gives the distance from the middle of the PSF to
	% the horizontal edge in the same units as optics.lambda0
	appear.psfxmax = optics.N*optics.lambda1/(4*optics.kmax);
	
	% xstep is the physical distance between x steps
	appear.psfxstep = appear.psfxmax/(optics.N/2);
	
	% just match up ends of k*range with OTF subN range for now, 
	% i.e. ledge to	redge
	ksqrange = linspace(-appear.psfxmax*Rp/2/sqrt(2), ...
						appear.psfxmax*Rp/2/sqrt(2), redge - ledge + 1);
	kdiagrange = linspace(-appear.psfxmax*Rp/2, appear.psfxmax*Rp/2, ...
							redge - ledge + 1);
	kcircrange = linspace(-appear.psfxmax*Rp/2, appear.psfxmax*Rp/2, ...
							credge - cledge + 1);
end

% save useful paramaters for return to caller
appear.kcircrange = kcircrange;
appear.kdiagrange = kdiagrange;
appear.ksqrange = ksqrange;
appear.ledge = ledge;
appear.redge = redge;
appear.cledge = cledge;
appear.credge = credge;
appear.x0 = x0;
appear.y0 = y0;

% setup 2D array where real part is x co-ord and imag part is y co-ord
% tricky N+1 and cut fiddle to get a zero at (N/2+1,N/2+1) i.e. (x0,y0)
% k1width is width of k=1 aperture (i.e. full aperture if alpha is max)
% round is needed for odd optics.N
appear.k1width = round(optics.N/optics.kmax);
[m n] = meshgrid(linspace(-1, 1, appear.k1width + 1));
m = m(1:appear.k1width, 1:appear.k1width);
n = n(1:appear.k1width, 1:appear.k1width);
% m increases from left to right
% n increases from top to bottom (reverse of usual)
% so when displaying things need to flipud() first
% (I tried to reverse it but couldn't get the spiral stuff to still work)
mn = m + i*n;

if strcmp(optics.pupilshape, 'circle')
	% create a filled circle of radius Rp: 1 inside, 0 outside
	circle = mn/Rp;
	circle(abs(circle) <= 1) = 1;
	circle(abs(circle) > 1) = 0;
	pupil = circle;
else
	% create a filled square that fits inside the circle
	sq = mn;
	sq(abs(m) <= Rp/sqrt(2)) = 1;
	sq(abs(n) <= Rp/sqrt(2)) = 1;
	sq(abs(m) > Rp/sqrt(2)) = 0;
	sq(abs(n) > Rp/sqrt(2)) = 0;
	pupil = sq;
end

if appear.savemem
	clear circle sq
end

if appear.savephase
	optics.pupilphase = zeros(optics.N/optics.kmax, optics.N/optics.kmax);
end

if strcmp(optics.z45, 'clear')
	z45filter = 1;
else
	R45 = sin(pi/4);
	if strcmp(optics.z45, 'opaque')
		z45filter = ones(optics.N/optics.kmax, optics.N/optics.kmax);
		z45filter(abs(m) > R45) = 0;
	else
		% assume optics.z45 = flip
		phase = zeros(optics.N/optics.kmax, optics.N/optics.kmax);
		phase(m > R45) = -pi/2;
		phase(m < R45) = pi/2;
		z45filter = exp(i*phase);
		if appear.savephase
			optics.pupilphase = optics.pupilphase + phase;
		end
		clear phase;
	end
	z45filter = z45filter.*pupil;
end

% debug: save z45filter to test that I am setting it right
if appear.debug
	optics.z45filter = z45filter;
end

l = abs(mn).*pupil;

splus = sqrt(1 - l.^2);

if optics.zdefocus ~= 0
	if strcmp(optics.model, 'paraxial')
		% 1 out the front of l.^2 has no physical effect (phase piston)
		% but makes it easier to compare the unwrapped pupil phase
		phase = k*optics.zdefocus*(1-l.^2/2);
		defocus = exp(i*phase);
	else
		% high NA defocus pupil multiplier
		phase = k*optics.zdefocus*splus;
		defocus = exp(i*phase);
	end
	optics.defocusmax = max(max(phase));
	optics.defocusmin = min(min(phase));
	optics.defocuspv = optics.defocusmax - optics.defocusmin;
	if appear.savephase
		optics.pupilphase = optics.pupilphase + phase;
	end
	clear phase;
else
	defocus = 1;
end

% DIC optics - see logbook 6nov02 pp 163-165
% start with shear in x direction
if optics.dicshear ~= 0
	if strcmp(optics.dicshearxy, 'x')
		dic = -2*i*sin(optics.dicbias/2 + k*optics.dicshear*0.5.*m);
	else
		dic = -2*i*sin(optics.dicbias/2 + k*optics.dicshear*0.5.*n);
	end
	if appear.savephase
		optics.pupilphase = optics.pupilphase + phase;
	end
else
	dic = 1;
end

if appear.savemem
	clear l mn
end

if strcmp(optics.model, 'vectorial')
	ax = (1 - m.^2./(1 + splus)).*pupil;
	ay = (-m.*n./(1 + splus)).*pupil;
	az = -m.*pupil;
else
	ax = 1; ay = 1; az = 1;
end

% create a cubic phase function
if Acpv ~= 0
	% we need to scale the Acpv range to lie between the pupil cutoffs
	% ah well another day
	phase = 2*pi*Acpv*(m.^3 + n.^3);
	cubic = exp(i*phase);
	if appear.savephase
		optics.pupilphase = optics.pupilphase + phase;
	end
	clear phase;
else
	cubic = 1;
end

if strcmp(optics.apodisation, 'sine') && ~strcmp(optics.model, 'paraxial')
	S = sqrt(splus);
else
	S = 1;
end

T = cubic.*z45filter.*dic.*defocus.*S.*pupil;

optics.cubicmax = 2*pi*Acpv*(max(max(m.*pupil))^3 + max(max(m.*pupil))^3);
optics.cubicmin = 2*pi*Acpv*(min(min(m.*pupil))^3 + min(min(m.*pupil))^3);
optics.cubicpv = optics.cubicmax - optics.cubicmin;

if appear.saveabns
	optics.cubic = cubic;
	optics.z45filter = z45filter;
	optics.dic = dic;
	optics.defocus = defocus;
	optics.S = S;
end

if appear.savephase
	optics.pupilphase = optics.pupilphase.*pupil;
	[ppdx, ppdy] = gradient(optics.pupilphase);
	optics.ppgrad = abs(ppdx + i*ppdy).*pupil;
	clear ppdx ppdy

	% crop phase gradient to just over pi
	% 'cos if it goes over pi we're stuffed anyway
	optics.ppgrad(abs(optics.ppgrad) > pi) = 3.5;
end

if appear.savemem
	clear cubic z45filter dic defocus S m n pupil
end

if strcmp(optics.model, 'vectorial')
	Qx = (1./splus).*ax.*T;
	Qy = (1./splus).*ay.*T;
	Qz = (1./splus).*az.*T;
elseif strcmp(optics.model, 'scalar')
	Qx = (1./splus).*T;
	Qy = Qx; Qz = Qx;
else
	% paraxial approximation
	Qx = T; 
	Qy = Qx; Qz = Qx;
end

if appear.savemem
	clear T ax ay az splus
end

if appear.showpupil
	pscale = [-1 1];

	figure(appear.figstart); 

	subplot(2,3,1);
	imagesc(abs(Qx)); 
	%axis square; axis xy; colorbar; colormap(jet);

	subplot(2,3,2);
	imagesc(abs(Qy)); 
	%axis square; axis xy; colorbar; colormap(jet);

	subplot(2,3,3);
	imagesc(abs(Qz)); 
	%axis square; axis xy; colorbar; colormap(jet);

	subplot(2,3,4);
	imagesc(angle(Qz)); 
	%axis square; axis xy; colorbar; colormap(jet);

	subplot(2,3,5);
	imagesc(angle(Qz)); 
	%axis square; axis xy; colorbar; colormap(jet);

	subplot(2,3,6);
	imagesc(angle(Qz)); 
	%axis square; axis xy; colorbar; colormap(jet);
end % appear.showplots

% pad pupils out to full size (optics.N) before doing each FFT chain
% impad() is from octave-forge

% might be worth trying to use the padding built into fft2 instead
% fft2(matrix, N, M) pads to size NxM

appear.padding = round((optics.N - optics.N/optics.kmax)/2);

Qx = impad(Qx, appear.padding, appear.padding);

if strcmp(appear.result, 'otf')
	C = fftshift(ifft2(abs(fft2(fftshift(Qx))).^2));
elseif strcmp(appear.result, 'ipsf')
	% need to add intensities to get vectorial PSF
	C = abs(fftshift(fft2(fftshift(Qx)))).^2;
elseif strcmp(appear.result, 'psf')
	C.x = fftshift(fft2(fftshift(Qx)));
elseif strcmp(appear.result, 'pupil')
	C.x = Qx;
end

if appear.saveabns
	optics.Qx = Qx;
end

if appear.savemem clear Qx; end

Qy = impad(Qy, appear.padding, appear.padding);

if strcmp(appear.result, 'otf')
	C = C + fftshift(ifft2(abs(fft2(fftshift(Qy))).^2));
elseif strcmp(appear.result, 'ipsf')
	C = C + abs(fftshift(fft2(fftshift(Qy)))).^2;
elseif strcmp(appear.result, 'psf')
	C.y = fftshift(fft2(fftshift(Qy)));
elseif strcmp(appear.result, 'pupil')
	C.y = Qy;
end

if appear.saveabns
	optics.Qy = Qy;
end

if appear.savemem clear Qy; end

Qz = impad(Qz, appear.padding, appear.padding);

if strcmp(appear.result, 'otf')
	C = C + fftshift(ifft2(abs(fft2(fftshift(Qz))).^2));
elseif strcmp(appear.result, 'ipsf')
	C = C + abs(fftshift(fft2(fftshift(Qz)))).^2;
elseif strcmp(appear.result, 'psf')
	C.z = fftshift(fft2(fftshift(Qz)));
elseif strcmp(appear.result, 'pupil')
	C.z = Qz;
end

if appear.saveabns
	optics.Qz = Qz;
end

if appear.savemem clear Qz; end

% normalise
if strcmp(appear.result, 'otf')
	C = C/abs(C(x0,y0));
end
 
if strcmp(appear.result, 'psf') || strcmp(appear.result, 'pupil')
	% plotting code not setup to deal with C.x, C.y, C.z, so
	% just force plots off for PSF case for now
	appear.showplots = 0;
end

if appear.showplots
	scale = [-2*Rp/sqrt(2) 2*Rp/sqrt(2)];

	if strcmp(appear.realabs, 'real')
		ptmp = real(C(ledge:redge,ledge:redge));
	else
		if strcmp(appear.result, 'otf')
			ptmp = abs(C(ledge:redge,ledge:redge));
		else
			% show intensity for the IPSF or image case - 
			% C is already intensity - see above
			ptmp = C;
		end
	end

	if appear.logscale
		ptmp = log10(ptmp);
	end

	figure(appear.figstart + 1); %clf reset;
	imagesc(scale, scale, ptmp); 
	%axis square; axis xy; colorbar; colormap(jet);

	figure(appear.figstart + 2); %clf reset;
	if appear.showphase
		imagesc(scale, scale, angle(C(ledge:redge,ledge:redge))); 
	end
	%axis square; axis xy; colorbar; colormap(jet);

	figure(appear.figstart + 3); %clf reset;

	% can plot 1D diagonals using something like this:
	% Ccpm = C;

	% dodgy normalise x axis to match BW99 Fig. 9.13
	% Rp = 1;

	diagC = ones(1,optics.N)*(C.*eye(optics.N));

	dtmp = diagC(ledge:redge);

	if strcmp(optics.pupilshape, 'circle')
		htmp = C(x0, cledge:credge);
		vtmp = C(cledge:credge, y0);
		axesrange = kcircrange;
	else
		htmp = C(x0, ledge:redge);
		vtmp = C(ledge:redge, y0);
		axesrange = ksqrange;
	end

	if strcmp(appear.realabs, 'real')
		dtmp = real(dtmp);
		htmp = real(htmp);
		vtmp = real(vtmp);
	else
		dtmp = abs(dtmp);
		htmp = abs(htmp);
		vtmp = abs(vtmp);
	end
	
	if appear.logscale
		semilogy(kdiagrange, dtmp,'r', axesrange, htmp,'g', axesrange, vtmp,'b');
		legend('diagonal', 'horizontal', 'vertical');
		axis('auto');
	else
		plot(kdiagrange, dtmp,'r', axesrange, htmp,'g', axesrange, vtmp,'b');
		legend('diagonal', 'horizontal', 'vertical');
		axis('auto');
		grid('on');
	end

	if appear.debug
		appear.axesrange = axesrange;
		appear.dtmp = dtmp;
		appear.htmp = htmp;
		appear.vtmp = vtmp;
	end
	
end % appear.showplots