function [S,Kappa1] = HayterPenfold(q,diam,volfrac,zz,csalt,Temp,dialec)

% function [S,Kappa1] = HayterPenfold(q,diam,volfrac,zz,csalt,Temp,dialec)
%
% This function calculates the structure factor for charged spherical particles
% using the rescaled mean spherical approximation (in Matlab)
% 
% OUTPUT:
%
% S = structure factor
% Kappa1 = Debye-Huckel screening length (in 1/m)
%
% INPUT: 
%
% q = q range (in A^-1)
% diam = diameter of sphere (in A)
% volfrac = volume fraction of spheres
% zz = charge (unitless)
% csalt = salt concentration (in molarity) (the code has not been tested with csalt>0!)
% Temp = temperature (in Kelvin)
% dialec = dielectric constant of solvent (unitless)
%
% This code is adapted from the Igor Pro SANS package (Kline, S. J. Appl. Cryst. 2006, 39, 895–900.)
% and is provided without any warranty. Use at your own risk.
%
% Created 3.9.2009 Ulla Vainio (ulla.vainio@desy.de)
% Edited: 1.3.2010 UV: Added effective kappa and effective charge formulae
% according to Hunter, page 683 or actually Snook & Hayter 

clear a b c f eta gek ak u v gamk seta sgek sak scal g1 fval evar
global a b c f eta gek ak u v gamk seta sgek sak scal g1 fval evar
% -------- Not defined in the original code, the initial number does not
% affect the result
a = 1;
b = 1;
c = 1;
f = 1;
u = 1;
v = 1;
gamk = 1;
seta = 1;
sgek = 1;
scal = 1;
fval = 1;
evar = 1;
g1 = 1; % Not defined in the original code of the Igor Pro SANS package
% --------
eta = volfrac;
Elcharge = 1.602189*10^-19; % in Coulombs
kB = 1.380662*10^-23; % Boltzman constant in J/K
FrSpPerm = 8.85418782*10^-12; % Permittivity of free space C^2/(N m^2)

QQ = q; % To same notation as in SANS package
Beta = 1/(kB*Temp); % in Joules
Perm = dialec*FrSpPerm; % in C^2/(N m^2)
SIdiam = diam*10^-10; % in m
Vp = 4*pi/3*(SIdiam/2)^3; % in m^3
cs = csalt*6.022*10^23*10^3; % salt molecules / m^3
% Ionic strength (C^2/m^3)
IonSt = 0.5*Elcharge^2*(zz*eta/Vp + 2*cs); % original formula in Igor Pro SANS package
% Debye-Huckel screening length in 1/m
Kappa1 = sqrt(2*Beta*IonSt/Perm);
Kappa = Kappa1 - 2*eta^(1/3)*log(1-eta); % Hunter's book, page 683 or Snook & hayter
charge = zz*Elcharge/(1-eta); % in Coulomb , Snook & Hayter
gek = Beta*charge^2/(pi*Perm*SIdiam*(2+Kappa1*SIdiam)^2); % wave 5

% Dimensionless parameters
Qdiam = QQ*diam;
ak = Kappa*SIdiam; % wave 6
if(ak>6)
    disp(sprintf('Kappa*diameter = %.2f > 6! D = %.1f',ak,diam));
end;
%disp(sprintf('Kappa*Diameter = %.2f (should be less than 6), dialec = %.2f',ak,2*Beta*IonSt/Kappa^2/FrSpPerm))
% eta is wave 4, volume fraction

% Calculate coupling
ss = eta^(1/3);
gamk = 2*ss*gek*exp(ak-ak/ss); % wave 9
%disp(sprintf('Coupling gamk = %.3f',gamk))

% Calculate coefficients, check all is well and calculate S(q*sig)

ierr = 0;
[ierr,u] = sqcoef(ierr);
%disp(sprintf('g1 (%.3e), gamma*exp(-k) (%.3e)',g1,gek))

for(mm = 1:length(Qdiam))
    if(ierr>=0)
     S(mm) = sqhcal(Qdiam(mm));
   else
     S(mm) = NaN;
     disp(sprintf('Error level %d',ierr));
     return
   end;
end;
% Checking up by calculating the radial distribution function
Qdiam = Qdiam(:);
S = S(:);

%counter = 1;
%r = [0.1:0.1:500];
%for(k = 1:length(r))
%   g(counter) = 1 + 1./(12*pi*eta*(r(k)/diam)).*trapz(Qdiam,(S-1).*Qdiam.*sin(Qdiam.*r(k)/diam));
%   counter = counter + 1;
%end;
%subplot(2,1,1);
%plot(Qdiam,S,'b')
%subplot(2,1,2);
%plot(r/diam,g,'g')
%disp(sprintf('%.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f',a, b, c, f, eta, gek, ak, u, v, gamk, seta, sgek, sak, scal, g1,fval,evar))

%--------------------------------- subfunction sqcoef
% calculates rescaled volume fraction and corresponding coefficients
function [ir,u] = sqcoef(ir)
global a b c f eta gek ak u v gamk seta sgek sak scal g1 fval evar

itm = 200; acc = 5*10^-6; del = 10; % itm = 40;
ig = 1;

if(ak>=(eta)) % What is this criterion for? original 1+8*eta
    ig = 0;
    fval = g1;
    evar = eta;
    ix = 1;
    [ir,u] = sqfun(ix,ir);
    g1 = fval;
    eta = evar;
    if(ir<0 | g1 >= 0)
        return;
    end;
end;

seta = min(eta,0.2);
f1 = fval; % Was not defined in the original code.
e1 = evar; % Was not defined in the original code.
f2 = fval; % Was not defined in the original code.
e2 = evar; % Was not defined in the original code.

% If potential is too high, rescale volume fraction and others

if(ig~=1 | gamk >=0.15) 
    ii = 0;
    while(del>acc) % Continue until reach accuracy
        if(ii > itm) % stop if maximum number of iterations is reached
            ir = -1;
            return;
        end;
        ii = ii + 1;
        if(seta<=0)
            seta = eta/ii;
        end;
        if(seta>0.6)
            seta = 0.35/ii;
        end;
        e1 = seta;
        fval = f1;
        evar = e1;
        ix = 2;
        ir = sqfun(ix,ir);
        f1 = fval;
        e1 = evar;
        e2 = seta*1.01;
        fval = f2;
        evar = e2;
        ix = 2;
        ir = sqfun(ix,ir);
        f2 = fval;
        e2 = evar;
        if(f2-f1 == 0)
            e2 = e1-(e2-e1)*f1/0.0000001;
        else
            e2 = e1-(e2-e1)*f1/(f2-f1);
        end;
        seta = e2;
        del = abs((e2-e1)/e1);
    end;
    fval = g1;
    evar = e2;
    ix = 4;
    ir = sqfun(ix,ir);
    g1 = fval;
    e2 = evar;
    ir = ii;
    if(ig~=1 & seta>=eta)
        return
    end;
end;
ix = 3;
[ir,u] = sqfun(ix,ir);
g1 = fval;
eta = evar;
if(ir>=0 & g1 < 0)
    ir = -3;
end;

%--------------------------------- subfunction sqfun
function [ir,u] = sqfun(ix,ir)
global a b c f eta gek ak u v gamk seta sgek sak scal g1 fval evar

acc = 10^-6; itm = 2000; % itm = 40;
reta = evar;
eta2 = reta^2;
eta3 = eta2*reta;
e12 = 12*reta;
e24 = 2*e12;
scal = (eta/evar)^(1/3);
sak = ak/scal;
ibig = 0;
if(sak > 15 & ix == 1)
    ibig = 1;
end;
sgek = gek*scal*exp(ak-sak);
rgek = sgek;
rak = sak;
ak2 = rak^2;
ak1 = 1+rak;
dak2 = 1/ak2;
dak4 = dak2^2;
d = 1-reta;
d2 = d^2;
dak = d/rak;
if(d2 == 0)
   dd2 = 1/0.000001;
else
   dd2 = 1/d2;
end;
dd4 = dd2^2;
dd45 = dd4*0.2;
eta3d = 3*reta;
eta6d = eta3d*2;
eta32 = eta3*2;
eta2d = reta + 2;
eta2d2 = eta2d^2;
eta21 = 2*reta + 1;
eta22 = eta21^2;

% Alpha(I)

al1 = -eta21*dak;
al2 = (14*eta2-4*reta-1)*dak2;
al3 = 36*eta2*dak4;

% Beta(I)

be1 = -(eta2+7*reta + 1)*dak;
be2 = 9*reta*(eta2+4*reta-2)*dak2;
be3 = 12*reta*(2*eta2+8*reta-1)*dak4;

% Nu(I)

vu1 = -(eta3+3*eta2+45*reta+5)*dak;
vu2 = (eta32+3*eta2+42*reta-20)*dak2;
vu3 = (eta32+30*reta-5)*dak4;
vu4 = vu1+e24*rak*vu3;
vu5 = eta6d*(vu2+4*vu3);

% Phi(I)

ph1 = eta6d/rak;
ph2 = d-e12*dak2;

% Tau(I)

ta1 = (reta+5)/(5*rak);
ta2 = eta2d*dak2;
ta3 = -e12*rgek*(ta1+ta2);
ta4 = eta3d*ak2*(ta1^2-ta2^2);
ta5 = eta3d*(reta+8)*0.1 - 2*eta22*dak2;

% Sinh(k), cosh(k)

ex1 = exp(rak);
ex2 = 0;
if(sak < 20)
    ex2 = exp(-rak);
end;
sk = 0.5*(ex1-ex2);
ck = 0.5*(ex1+ex2);
ckma = ck - 1 - rak*sk;
skma = sk - rak*ck;

% a(I)

a1 = (e24*rgek*(al1+al2+ak1*al3)-eta22)*dd4;
if(ibig == 0)
    a2 = e24*(al3*skma + al2*sk - al1*ck)*dd4;
    a3 = e24*(eta22*dak2-0.5*d2+al3*ckma-al1*sk+al2*ck)*dd4;
end;

% b(I)

b1 = (1.5*reta*eta2d2 - e12*rgek*(be1+be2+ak1*be3))*dd4;
if(ibig == 0)
    b2 = e12*(-be3*skma-be2*sk+be1*ck)*dd4;
    b3 = e12*(0.5*d2*eta2d-eta3d*eta2d2*dak2-be3*ckma+be1*sk-be2*ck)*dd4;
end;

% V(I)

v1 = (eta21*(eta2-2*reta+10)*0.25 - rgek*(vu4+vu5))*dd45;
if(ibig == 0)
    v2 = (vu4*ck-vu5*sk)*dd45;
    v3 = ((eta3-6*eta2+5)*d-eta6d*(2*eta3-3*eta2+18*reta+10)*dak2+e24*vu3+vu4*sk-vu5*ck)*dd45;
end;

%sprintf('%.4f %.4f %.4f %.4f %.4f %.4f %.4f %.4f %.4f',a1,a2,a3,b1,b2,b3,v1,v2,v3)


% P(I)

pp1 = ph1^2;
pp2 = ph2^2;
pp = pp1 + pp2;
p1p2 = ph1*ph2*2;
p1 = (rgek*(pp1+pp2-p1p2)-0.5*eta2d)*dd2;
if(ibig == 0)
    p2 = (pp*sk + p1p2*ck)*dd2;
    p3 = (pp*ck + p1p2*sk + pp1 - pp2)*dd2;
end;

% T(I)

t1 = ta3+ta4*a1+ta5*b1;

if(ibig~=0)
    % very large screening: asymptotic solution
    a3 = e24*(eta22*dak2-0.5*d2-al3)*dd4; % in original code these (a3 and b3) were later
    b3 = e12*(0.5*d2*eta2d-eta3d*eta2d2*dak2+be3)*dd4;
    v3 = ((eta3-6*eta2+5)*d-eta6d*(2*eta3-3*eta2+18*reta + 10)*dak2 + e24*vu3)*dd45;
    t3 = ta4*a3+ta5*b3+e12*ta2-0.4*reta*(reta+10)-1;
    p3 = (pp1-pp2)*dd2;
    um6 = t3*a3-e12*v3^2;
    um5 = t1*a3+a1*t3-e24*v1*v3;
    um4 = t1*a1-e12*v1^2;
    al6 = e12*p3^2;
    al5 = e24*p1*p3-2*b3-ak2;
    al4 = e12*p1^2-2*b1;
    w56 = um5*al6-al5*um6;
    w46 = um4*al6-al4*um6;
    fa = -w46/w56;
    ca = -fa;
    f = fa;
    c = ca;
    b = b1+b3*fa;
    a = a1+a3*fa;
    v = v1+v3*fa;
    g1 = -(p1+p3*fa);
    fval = g1;
    if(abs(fval)<0.001)
        fval = 0;
    end;
    seta = evar;
%    disp('Very large screening!')
else
    t2 = ta4*a2+ta5*b2+e12*(ta1*ck-ta2*sk);
    t3 = ta4*a3+ta5*b3+e12*(ta1*sk-ta2*(ck-1))-0.4*reta*(reta+10)-1;
   
    % Mu(I)
    
    um1 = t2*a2-e12*v2^2;
    um2 = t1*a2+t2*a1-e24*v1*v2;
    um3 = t2*a3+t3*a2-e24*v2*v3;
    um4 = t1*a1-e12*v1^2;
    um5 = t1*a3+t3*a1-e24*v1*v3;
    um6 = t3*a3-e12*v3^2;
    
    % Gillan condition?
    % Yes ---- G(X = 1+) = 0
    % Coefficients and function value
    if(ix == 1 | ix == 3)
        % No, calculate remaining coefficients
        % Lambda(I)
        al1 = e12*p2^2;
        al2 = e24*p1*p2-2*b2;
        al3 = e24*p2*p3;
        al4 = e12*p1^2-2*b1;
        al5 = e24*p1*p3 - 2*b3-ak2;
        al6 = e12*p3^2;
        
        % Omega(I)
        w16 = um1*al6-al1*um6;
        w15 = um1*al5-al1*um5;
        w14 = um1*al4-al1*um4;
        w13 = um1*al3-al1*um3;
        w12 = um1*al2-al1*um2;
        
        w26 = um2*al6-al2*um6;
        w25 = um2*al5-al2*um5;
        w24 = um2*al4-al2*um4;
        
        w36 = um3*al6-al3*um6;
        w35 = um3*al5-al3*um5;
        w34 = um3*al4-al3*um4;
        w32 = um3*al2-al3*um2;
        
        w46 = um4*al6-al4*um6;
        w56 = um5*al6-al5*um6;
        w3526 = w35+w26;
        w3425 = w34+w25;
        
        % Quadratic coefficients
        
        w4 = w16^2-w13*w36;
        w3 = 2*w16*w15-w13*w3526-w12*w36;
        w2 = w15^2+2*w16*w14-w13*w3425-w12*w3526;
        w1 = 2*w15*w14-w13*w24-w12*w3425;
        w0 = w14^2-w12*w24;
        
        % Estimate starting value of f
        
        if(ix ~= 1)
            % Assume not too far from Gillan condition.
            % if both rgek and rak as small, use P-W estimate.
            g1 = 0.5*eta2d*dd2*exp(-rgek);
            if(sgek<=2 & sgek >= 0 & sak <= 1)
                e24g = e24*rgek*exp(rak);
                pwk = sqrt(e24g);
                qpw = (1-sqrt(1+2*d2*d*pwk/eta22))*eta21/d;
                g1 = -qpw^2/e24+0.5*eta2d*dd2;
            end;
            pg = p1+g1;
            ca = ak2*pg+2*(b3*pg-b1*p3)+e12*g1^2*p3;
            ca = -ca/(ak2*p2+2*(b3*p2-b2*p3));
            fap = -(pg+p2*ca)/p3; 
         else
            % Large k
            if(w12-w15+w35-w26+w56-w32 ~=0)
               fap = (w14-w34-w46)/(w12-w15+w35-w26+w56-w32);
            else
               % Assume not too far from Gillan condition.
             % if both rgek and rak as small, use P-W estimate.
               g1 = 0.5*eta2d*dd2*exp(-rgek);
               if(sgek<=2 & sgek >= 0 & sak <= 1)
                   e24g = e24*rgek*exp(rak);
                   pwk = sqrt(e24g);
                   qpw = (1-sqrt(1+2*d2*d*pwk/eta22))*eta21/d;
                   g1 = -qpw^2/e24+0.5*eta2d*dd2;
               end;
               pg = p1+g1;
               ca = ak2*pg+2*(b3*pg-b1*p3)+e12*g1^2*p3;
               ca = -ca/(ak2*p2+2*(b3*p2-b2*p3));
               fap = -(pg+p2*ca)/p3;           
           end;
        end;

        % And refine it according to Newton
        
        ii = 0;
        del = 10;
        while(del>acc)
            if(ii > itm)
                % Failed to converge in itm iterations
                ir = -2;
                return
            end;
            ii = ii + 1;
            fa = fap;
            fun = w0+(w1+(w2+(w3+w4*fa)*fa)*fa)*fa;
            fund = w1 + (2*w2+(3*w3+4*w4*fa)*fa)*fa;
            if(fund == 0)
               fap = fa - fun/0.000001;
            else
               fap = fa - fun/fund;
            end;
            if(fa == 0)
               del = abs((fap)/0.0000001);                
            else
               del = abs((fap-fa)/fa);
            end;
        end;
        ir = ir + ii;
        fa = fap;
        ca = -(w16*fa^2+w15*fa+w14)/(w13*fa+w12);
        g1 = -(p1+p2*ca+p3*fa);
        fval = g1;
        if(abs(fval)<0.001)
            fval = 0;
        end;
        seta = evar;
    else
        ca = ak2*p1+2*(b3*p1-b1*p3);
        ca = -ca/(ak2*p2+2*(b3*p2-b2*p3));
        fa = -(p1+p2*ca)/p3;
        if(ix == 2)
            fval = um1*ca^2+(um2+um3*fa)*ca+um4+um5*fa+um6*fa^2;
        end;
        if(ix == 4)
            fval = -(p1+p2*ca+p3*fa);
        end;
    end;
    f = fa;
    c = ca;
    b = b1+b2*ca+b3*fa;
    a = a1+a2*ca+a3*fa;
    v = (v1+v2*ca+v3*fa)/a;
end;
g24 = e24*rgek*ex1;
u = (rak*ak2*ca-g24)/(ak2*g24);


%--------------------------------- subfunction sphcal
function S = sqhcal(qq)
global a b c f eta gek ak u v gamk seta sgek sak scal g1 fval evar

etaz = seta;
akz = sak;
gekz = sgek;
e24 = 24*etaz;
x1 = exp(akz);
x2 = 0;
if(sak<20)
    x2 = exp(-akz);
end;
ck = 0.5*(x1+x2);
sk = 0.5*(x1-x2);
ak2 = akz^2;
if(qq <= 0)
     S = -1/a;
else
    qk = qq/scal;
    q2k = qk.^2;
    qk2 = 1./q2k;
    qk3 = qk2./qk;
    qqk = 1./(qk.*(q2k+ak2));
    sink = sin(qk);
    cosk = cos(qk);
    asink = akz.*sink;
    qcosk = qk.*cosk;
    aqk = a*(sink-qcosk);
    aqk = aqk + b*((2*qk2-1).*qcosk+2*sink-2./qk);
    inter = 24*qk3+4*(1-6*qk2).*sink;
    aqk = (aqk+0.5*etaz*a.*(inter-(1-12*qk2+24*qk2.^2).*qcosk)).*qk3;
    aqk = aqk + c*(ck.*asink-sk.*qcosk).*qqk;
    aqk = aqk + f*(sk.*asink-qk.*(ck.*cosk-1)).*qqk;
    aqk = aqk + f*(cosk-1).*qk2;
    aqk = aqk - gekz.*(asink+qcosk).*qqk;
    S = 1./(1-e24.*aqk);
end;