function [F,J,Eccc,Ecc,c] = FDF(csmall,gamma,m,ell,K,SG)
% [F,J,Eccc,Ecc,c] = FDF(csmall,gamma,m,ell,K,SG)
% computes the F=F(csmall) where csmall is in symmetry reduced form
% as well as the Jacobian J at that point.
% The physical parameter values are gamma, m and ell 
% (see the description of the Ohta-Kawasaki model in the paper).
%
% The integer K determines the size of the finite dimensional projection
% The symmetry group SG can be either 230 or 229.
% The default symmetry group is 230.
%
% The function works both for floating point numbers and for intervals
%
% The extra outputs are the extended convolutions (all nonvanishing terms)
% Eccc for c^3 and Ecc for c^2 as well as c itself
% all in tensor format

if nargin == 5
    SG = 230;
    disp(['Space Group not set, defaulting to ',int2str(SG)]);
else
    disp(['The Space Group is ',int2str(SG)]);
end;

[kx,ky,kz]=ndgrid(-K:K);
% choose fundamental domain
ksym=(kx>=ky & ky>=kz & kz>=0);

% determine sum(alpha_g(k) , g in G)
symmetrizedcount=symmetrize(ksym,SG);
% determine the number of elements in the stabilizer
stabilizercount=abs(symmetrizedcount);
stabilizercount(stabilizercount==0)=1; % prevent division by zero
% tildealpha keeps track of prefactors (see paper)
tildealpha=symmetrizedcount./stabilizercount;

% find k that are not in Z_{triv} 
% within the ball of radius K 
% and exclude k=0
k2=kx.^2+ky.^2+kz.^2;
ksmall=(symmetrizedcount~=0 & k2<=K^2 & k2>0);
% find symmetry reduced indices 
ksymsmall=(ksym & ksmall);

ksymindexsmall=(1:length(find(ksymsmall)))';
ksymindex=zeros([2*K+1,2*K+1,2*K+1]);
ksymindex(ksymsmall)=ksymindexsmall;

% ksymindex keeps track of to which of the group orbits 
% (which are indexed by the symmetry reduced indices in ksymsmall) 
% elements in ksmall belong
ksymindex=abs(symmetrize(ksymindex,SG))./stabilizercount;

c=ILzeros([2*K+1,2*K+1,2*K+1],csmall(1));
c(ksymsmall)=csmall;
% determine tensor representation 
c=symmetrize(c,SG);
c=c(:)./stabilizercount(:);  
% reshape necessary to work with intlab
c=reshape(c,2*K+1,2*K+1,2*K+1); 
% since c_0=m is not a variable we need to put it in explicitly
c(K+1,K+1,K+1)=m; 

lk=k2/ell^2;
lambda=gamma^(-2)*lk.^2-lk+1;

[ccc,Eccc]=powconv3(c,3);
% we explicitly use the fact that we have real Fourier coefficients here
ccc=real(ccc);
F=lambda.*c +lk.*ccc;
F=F(ksymsmall);

[~,Ecc]=powconv3(c,2);
% we explicitly use the fact that we have real Fourier coefficients here
Ecc=real(Ecc);
JM=length(find(ksymsmall));

% determine the derivative of the convolution c^3 
% (factor 3 is taken care of below)
J=ILzeros([JM,JM],csmall(1));
Sc=2*K+1;
Sx0=Sc+kx(ksymsmall);
Sy0=Sc+ky(ksymsmall)-1;
Sz0=Sc+kz(ksymsmall)-1;
My=4*K+1;
Mz=(4*K+1)^2;
for j=find(ksmall)'
  kj=ksymindex(j);
  % uses shift property of convolution property
  Sx=Sx0-kx(j);
  Sy=Sy0-ky(j);
  Sz=Sz0-kz(j);
% Sx+My*Sy+Mz*Sz turns tensor indexing into linear indexing
% so that the entire kj-th column of values corresponding to 
% symmetry reduced indices can be extracted from the tensor Ecc
  J(:,kj)=J(:,kj)+tildealpha(j)*Ecc(Sx+My*Sy+Mz*Sz);
end

% full Jacobian 
J=diag(lambda(ksymsmall))+3*diag(lk(ksymsmall))*J;

end