function [q,I,varargout]=perfectpattern_spheres(q,varargin)

%this function calculates the scattering pattern of a distributino of
%spheres. It does so by random generation of spheres until the changes in
%the scattering pattern are smaller than the convergence criterium.
%Implementation inspired by the work of Mulato and Chambouleyron, j. appl.
%cryst. 1996, 29, 29-36. Used sphere scattering function is not the
%function with the Bessel function of the first kind as suggested by
%Mulato, but instead the quicker Rayleigh function. Regarding the size
%distribution function, negative values are rejected by default, not folded into the
%positive region to ensure shape preservation of the distribution.
%Distributions are radius distributions.
%
%first input argument is the q (n-value vector) or q range (2-value vector)
%   of interest. later arguments are parameter-value pairs, with possible:
%'dtype', 'norm' by default, can be anything of the function "random" in
%   the statistics toolbox, or 'gaussian'(default), 'lognormal', 'Schultz'* which are
%   separately implemented (and do not reuquire said toolbox).
%'dparam', a vector with parameters defining the distribution
%*'negfolded', 'true' or 'false'(default)
%'nperblock', number of spheres rendered at once, default 100
%'convcrit', convergence criterium, default 1e-3
%'timeout',if it takes longer than n seconds, quit. The program lets you
%   know the convergence to indicate how close it was.
%'nq', if q is a range, nq indicates the division into n points. 100 def.
%'vwdist', if set to 'true' It is assumed that you have taken the volume weighting into
%   account when calculating the distribution parameters (possible for some).
%   Volume is no longer taken into account in the integration.
%'param', saved parameters from a previous simulation
%Output consists of q, the intensity, and optional arguments:
%r: sizes of particles used in the simulation
%it: number of steps before the convergence criteria was achieved
%param: saved settings of the distribution so it can be easily recreated
%(*not implemented yet)


%set parameters
param.shape='sphere';
param.nq=100;
param.timeout=30; %seconds
param.convcrit=1e-5;
param.negfolded='false'; %not implemented yet.. Not sure if I will..
param.nperblock=100;
param.vwdist='false'; %eventually, I would like to change this into "internal" and "external" and "false", with the second being identical to the current "true", but the first using internal conversion of the fitting parameters.
param.RD.dtype='gaussian'; %RD because later on we want to be able to model other shapes as well
param.RD.dparam=[10 1];%mu sigma

param=paramproc(param,varargin); %process input parameters

kansas=size(q);%something to return to

%process q
if length(q)==2 %flesh it out
    q=min(q):(max(q)-min(q))/(param.nq-1):max(q);
    %get kansas right
    if size(q,1)>size(q,2) %original was in nx1 shape, so let us output that way as well
        kansas=[length(q) 1];
    else
        kansas=[1 length(q)];
    end
    if min(q)==0
        warning('input:badzero','If you let q start at zero, the convergence criteria may not be reached')
    end
end
%reshape q so it is a 1xn vector, i.e. one row
q=reshape(q,[1 size(q,1)*size(q,2)]);

it=1; %first iteration
conv=inf; %convergence criteria
r=[];
tic;
elapsed=0;
Vwsqrtot=0;
while conv>param.convcrit && elapsed<param.timeout;
    R=rsphgen(param); % R should be a mx1 vector, i.e. one column
    r(length(r)+(1:length(R)))=R;
    if it==1
        Iint=zeros(size(q));
    end

    %qR=repmat(q,[length(R) 1]).*repmat(R,[1 length(q)]);
    qR=R*q;
    Fsph=sphereintensity_rayleigh(qR);
    if strcmp(param.vwdist,'false')
        V=(4/3*pi).*R.^3;
    else
        V=ones(size(R));
    end
    %Vwsqrtot=Vwsqrtot+sum(V.^2);
    Iint=(Iint.*(it-1)+sum(repmat(V,[1 length(q)]).^2.*abs(Fsph).^2,1)./sum(V.^2))./it;


    %compute convergence if number of iterations larger than two
    if it==1
        Iold=Iint;
    elseif it>1
        %calculate the convergence criterium
        conv=max(abs(1-Iint/Iold));
        Iold=Iint; %new is the new old
    end
    it=it+1; %number of iterations
    elapsed=toc;
end
if conv>param.convcrit %meaning we ran out of time..
    warning('iterative:convnotreached','Convergence criteria not reached, timeout reached instead.\n Number of iterations %i, Convergence value %f',it,conv);
end

%return the right shape
I=reshape(Iint,kansas);
switch nargout
    case 3
        varargout{1}=r;
    case 4
        varargout{1}=r;
        varargout{2}=it-1; %remove the last it that gets added
    case 5
        varargout{1}=r;
        varargout{2}=it-1;
        varargout{3}=param;
end

end


function Fsph=sphereintensity_rayleigh(qR)

%this is the intensity calcultion at q for a sphere with radius R.
%as given in the paper by J.S. Pedersen, 1997, (eqn. 55)

Fsph=(qR>=0.05).*3.*(sin(qR)-qR.*cos(qR))./qR.^3+(qR<0.05).*(1-0.1.*qR.^2);
end


function R=rsphgen(param)
R=-1.*ones(param.nperblock,1);
while min(R)<=0
    %find out how many values we need to generate:
    R=sort(R);
    maxi=find(R<=0,1,'last'); %this is how many we need to regenerate, because they fall below zero
    
    switch param.RD.dtype
        case 'strict'
            Rn=param.RD.dparam(1).*ones(maxi,1); %the easiest option of all
        case 'none'
            Rn=rand(maxi,1); %not sure if this works, may not produce nx1 matrix
            warning('randomgen:untested','functioning of generator not tested');

        case 'gaussian'
            Rn=normrnd(param.RD.dparam(1),param.RD.dparam(2),maxi,1);
        case 'lognormal'
            R1Mean=param.RD.dparam(1); %mean
            R1Variance=param.RD.dparam(2); %variance
            mu=log(R1Mean^2/(sqrt(R1Variance+R1Mean^2))); %mu
            sigma=sqrt(log(R1Variance/(R1Mean^2)+1)); %sigma
            Rn=lognrnd(mu,sigma,maxi,1);
        otherwise
            switch length(param.RD.dparam)
                case 1
                    p1=param.RD.dparam(1);
                    Rn=random(param.RD.dtype,p1,maxi,1);
                case 2
                    p1=param.RD.dparam(1);
                    p2=param.RD.dparam(2);
                    Rn=random(param.RD.dtype,p1,p2,maxi,1);
                case 3
                    p1=param.RD.dparam(1);
                    p2=param.RD.dparam(2);
                    p3=param.RD.dparam(3);
                    Rn=random(param.RD.dtype,p1,p2,p3,maxi,1);
            end
    end
    
    %now fill in the random generated values wherever we need to 
    R(1:maxi)=Rn;
end
%just in case
R=reshape(R,[length(R) 1]);
end

function param=paramproc(param,vargi)
vi=0;
npv=length(vargi);
while vi<npv
    vi=vi+1;
    switch vargi{vi}
        case 'dtype'
            vi=vi+1;
            param.RD.dtype=vargi{vi};
        case 'dparam'
            vi=vi+1;
            param.RD.dparam=vargi{vi};
        case 'negfolded'
            vi=vi+1;
            param.negfolded=vargi{vi};
            warning('inputprocessing:notimplementedyet','negfolded not implemented yet')
        case 'nperblock'
            vi=vi+1;
            param.nperblock=vargi{vi};
        case 'convcrit'
            vi=vi+1;
            param.convcrit=vargi{vi};
        case 'timeout'
            vi=vi+1;
            param.timeout=vargi{vi};
        case 'nq'
            vi=vi+1;
            param.nq=vargi{vi};
        case 'vwdist'
            vi=vi+1;
            param.vwdist=vargi{vi};
        case 'param'
            vi=vi+1;
            param=vargi{vi};
            warning('inputprocessing:usagewarning','only parameter-value pairs after the param settings will be adjusted!')


        otherwise
            vi=vi+1;
            warning('inputprocessing:invalidparameter','parameter %s is not one of the acceptable parameters. Check your usage',vargi{vi-1});
    end
end

end
