function [Dx_phase,C_phase,M_phase,minphase,Dx_maxcorr,C_maxcorr,M_maxcorr,Dx_maxE] = FindDisparityFT(ikx,iky,FTimL,FTimR,FTRFE,FTRFO,xpref,ypref,Dxprefs)
% Estimate the disparity of the stereogram imL,imR at position (xpref,ypref) using the RFs specified.
% Evaluates the convolution of images and RFs in Fourier space
% I have only written it for isotropic Gaussian envelopes

% All cells have the same preferred cyclopean position
% (xpref,ypref) and preferred vertical disparity zero, but differ in their
% preferred horizontal disparity and phase disparity.

% Dxprefs is an array of position disparities first used to find maxima.

nDx = length(Dxprefs);
options=optimset('tolx',1e-6);

for jdx=1:nDx
    [zerophasedispE(jdx),zerophaseC(jdx)] = EnergyModel(Dxprefs(jdx),ikx,iky,FTimL,FTimR,FTRFE,FTRFO,xpref,ypref,0);
end
% Locate turning points:
df=diff(zerophasedispE);
jlocmax = find((df<0&shiftx(df,1)>0) | (df==0 & shiftx(df,-1)<0&shiftx(df,1)>0));
jlocmin = find((df>0&shiftx(df,1)<0) | (df==0 & shiftx(df,-1)>0&shiftx(df,1)<0));
% Now look for these turning points more accurately:
Dxprefmin = [];
for j=1:length(jlocmin)
    [Dxprefmin(j),fval,exitflag] = fminsearch( @EnergyModel,Dxprefs(jlocmin(j)),options, ikx,iky,FTimL,FTimR,FTRFE,FTRFO,xpref,ypref,0 );
    if exitflag<0
        disp('failed when finding minimum E as a function of pos-disp')
    end
end
Dxprefmax=[];
minusEpeak=[];
for j=1:length(jlocmax)
    [Dxprefmax(j),minusEpeak(j),exitflag] = fminsearch( @MinusEnergyModel,Dxprefs(jlocmax(j)),options,ikx,iky,FTimL,FTimR,FTRFE,FTRFO,xpref,ypref,0 );
    if exitflag<0
        disp('failed when finding maximum E as a function of pos-disp')
    end
end
% Reject TPs outside the range it was told to look for:
Dxprefmin = Dxprefmin(Dxprefmin>=min(Dxprefs) & Dxprefmin<=max(Dxprefs)); 
Dxprefmax = Dxprefmax(Dxprefmax>=min(Dxprefs) & Dxprefmax<=max(Dxprefs)); 
minusEpeak = minusEpeak(Dxprefmax>=min(Dxprefs) & Dxprefmax<=max(Dxprefs)); 

% I'm alos going to find accurately the peak in corrtelation , for comparison
% Peak in C is estimated:
[mx,j] = max(zerophaseC);
if j==1 | j==nDx
    % then the maximum in C was the last value in the array, ie is not a true maximum, so ignore:
    Dx_maxcorr = NaN;
    E_maxcorr = NaN;
    C_maxcorr = NaN;
    M_maxcorr = NaN;
else
    Dx_maxcorr = Dxprefs(j);
    [Dx_maxcorr,fval,exitflag] = fminsearch( @MinusCorrelation,Dx_maxcorr,options, ikx,iky,FTimL,FTimR,FTRFE,FTRFO,xpref,ypref,0 );
    % Finally record value of monoc responses at this fitted disparity:
    [E_maxcorr,C_maxcorr,M_maxcorr] = EnergyModel(Dx_maxcorr,ikx,iky,FTimL,FTimR,FTRFE,FTRFO,xpref,ypref,0);
end

% If there was at least 1 local maximum in E, record it for the maxE method:
if length(Dxprefmax)>0
    % Find biggest peak:
    [mx,j] = min(minusEpeak);
    Dx_maxE = Dxprefmax(j);
else
    Dx_maxE = max(zerophasedispE);
end

% Amalgamate all TPs:
DxprefTP = sort([Dxprefmin Dxprefmax]);
nTP = length(Dxprefmin)+length(Dxprefmax);
% If no TPs at all, return NaN:    
if nTP==0
   [Dx_phase,minphase,E_phase,C_phase,M_phase]  = deal(NaN);
    return
end
DxprefTP = sort([Dxprefmin Dxprefmax]);

% For each turning-point on the zero-phase-disp line, calculate the response of phase disparity detectors:
for jTP = 1:nTP
    % Find phase disparity at which E peaks, keeping position disparity constant:
    [fittedphase(jTP),fval,exitflag] = fminsearch(@negEasfnofphase,0,options,ikx,iky,FTimL,FTimR,FTRFE,FTRFO,xpref,ypref,DxprefTP(jTP));
    if exitflag<0
        disp('failed when finding maximum E as a function of phase-disp')
    end
end

% Now return the disparity of the TP with the smallest absolute phase:
[minphase,jfit]=min(abs(fittedphase));
Dx_phase = DxprefTP(jfit);
% Finally record value of monoc responses at this fitted place:
[E_phase,C_phase,M_phase] = EnergyModel(Dx_phase,ikx,iky,FTimL,FTimR,FTRFE,FTRFO,xpref,ypref,0);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function E = MinusEnergyModel(Dxpref,ikx,iky,FTimL,FTimR,FTRFE,FTRFO,xpref,ypref,Dphipref)
E = - EnergyModel(Dxpref,ikx,iky,FTimL,FTimR,FTRFE,FTRFO,xpref,ypref,Dphipref);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [E,varargout] = EnergyModel(Dxpref,ikx,iky,FTimL,FTimR,FTRFE,FTRFO,xpref,ypref,Dphipref)
% Energy model as function of postion disparity
xL = xpref - Dxpref/2;
xR = xpref + Dxpref/2;

cs = cos(Dphipref/2);
sn = sin(Dphipref/2);
LO = sum(FTRFO.*exp(ikx.*xL+iky.*ypref).*FTimL);
RO = sum(FTRFO.*exp(ikx.*xR+iky.*ypref).*FTimR);
LE = sum(FTRFE.*exp(ikx.*xL+iky.*ypref).*FTimL);
RE = sum(FTRFE.*exp(ikx.*xR+iky.*ypref).*FTimR);
% RF 1
L1 = real(cs * LE + sn * LO);
R1 = real(cs * RE - sn * RO);
% RF 2
L2 = real(-cs * LO + sn * LE);
R2 = real(-cs * RO - sn * RE);

E = (L1+R1)^2+(L2+R2)^2;
if nargout>1
    varargout{1} = 2*(L1*R1+L2*R2)/(L1^2+L2^2+R1^2+R2^2);
end
if nargout>2
    varargout{2} = (L1^2+L2^2+R1^2+R2^2);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function E = negEasfnofphase(Dphipref,ikx,iky,FTimL,FTimR,FTRFE,FTRFO,xpref,ypref,Dxpref)
% minus Energy model ouput as function of phase disparity
xL = xpref - Dxpref/2;
xR = xpref + Dxpref/2;

cs = cos(Dphipref/2);
sn = sin(Dphipref/2);
LO = sum(FTRFO.*exp(ikx.*xL+iky.*ypref).*FTimL);
RO = sum(FTRFO.*exp(ikx.*xR+iky.*ypref).*FTimR);
LE = sum(FTRFE.*exp(ikx.*xL+iky.*ypref).*FTimL);
RE = sum(FTRFE.*exp(ikx.*xR+iky.*ypref).*FTimR);

% RF 1
L1 = real(cs * LE + sn * LO);
R1 = real(cs * RE - sn * RO);
% RF 2
L2 = real(cs * LO - sn * LE);
R2 = real(cs * RO + sn * RE);

E = -(L1+R1)^2-(L2+R2)^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function mC = MinusCorrelation(Dxpref,ikx,iky,FTimL,FTimR,FTRFE,FTRFO,xpref,ypref,Dphipref)
xL = xpref - Dxpref/2;
xR = xpref + Dxpref/2;

cs = cos(Dphipref/2);
sn = sin(Dphipref/2);
LO = sum(FTRFO.*exp(ikx.*xL+iky.*ypref).*FTimL);
RO = sum(FTRFO.*exp(ikx.*xR+iky.*ypref).*FTimR);
LE = sum(FTRFE.*exp(ikx.*xL+iky.*ypref).*FTimL);
RE = sum(FTRFE.*exp(ikx.*xR+iky.*ypref).*FTimR);

% RF 1
L1 = real(cs * LE + sn * LO);
R1 = real(cs * RE - sn * RO);
% RF 2
L2 = real(cs * LO - sn * LE);
R2 = real(cs * RO + sn * RE);
mC = -2*(L1*R1+L2*R2)/(L1^2+L2^2+R1^2+R2^2);
