function CaF2PeakIntensity

Lambda = 1.54056;

aCa = [ 8.627 7.387 1.590 1.021];
bCa = [10.442 0.660 85.748 178.437];
cCa = 1.375;

aF = [ 3.539 2.641 1.517 1.024];
bF = [ 10.283 4.294 0.262 26.148];
cF = 0.278;

xx = 0:0.001:0.6;
fF = 0.0;
fCa = 0.0;
for ii=1:4
    fF = fF + aF(ii)*exp(-bF(ii)*xx.^2);
    fCa = fCa + aCa(ii)*exp(-bCa(ii)*xx.^2);
end
figure(1);
plot(xx,fF);
ylabel('Electron Scattering Factor F');
xlabel('sin\theta/\lambda, angstroms^-^1');
title('Atomic Scattering Factors for Fluorine');

figure(2);
plot(xx,fCa);
ylabel('Electron Scattering Factor Ca');
xlabel('sin\theta/\lambda, angstroms^-^1');
title('Atomic Scattering Factors for Calcium');

a = 5.4626;
a = 5.451;
fprintf(1,'Lambda=%f Angstroms a=%f Angstroms\n',Lambda,a);
fprintf(1,'(h,k,l) M d(Ang)  T(Deg) fF   fCa    nfF   nfCa  S(hkl) LP(T) Intensity\n');
InData = [[1 1 1]; [2 2 0]; [2 0 0]; [1 0 0]; [1 1 0]; [2 1 0]];
charge = 0.0;

for iplane  = 1:6

h=InData(iplane,1);
k=InData(iplane,2);
l=InData(iplane,3);

if (h==k || h==0 | k==0)
    Multiplicity = 4;
else
    Multiplicity = 8;
end
if (l ~= 0)
    Multiplicity = Multiplicity*2;
end
img = sqrt(-1);
d = 1.0/sqrt( (h^2)/(a^2) + (k^2)/(a^2) + (l^2)/(a^2) );

uCa = [0.0 0.0 0.5 0.5];
vCa = [0.0 0.5 0.0 0.5];
wCa = [0.0 0.5 0.5 0.0];

uF = [ 0.25 -0.25  0.25 -0.25  0.75 -0.75  0.75 -0.75];
vF = [ 0.75 -0.75 -0.75  0.75  0.25 -0.25 -0.25  0.25];
wF = [-0.25 -0.25  0.25  0.25  0.25  0.25 -0.25 -0.25];

nfCa = 0.0;
for iatom=1:4
    nfCa = nfCa + exp(-img*2*pi*( h*uCa(iatom) + k*vCa(iatom) + l*wCa(iatom)));
end
nfF = 0.0;
for iatom=1:8
    nfF = nfF + exp(-img*2*pi*( h*uF(iatom) + k*vF(iatom) + l*wF(iatom)));
end

Theta = asin(Lambda/(2*d));
xx = sin(Theta)/Lambda;
fF = 0.0;
fCa = 0.0;
for ii=1:4
    fF = fF + aF(ii)*exp(-bF(ii)*xx^2);
    fCa = fCa + aCa(ii)*exp(-bCa(ii)*xx^2);
end
%fCa = (20-2*charge)*(1.0-sin(Theta)/Lambda);
%fF = (9+charge)*(1.0-sin(Theta)/Lambda);
S = nfF*fF + nfCa*fCa;
LP = (1.0  + (cos(2*Theta))^2)/(8*cos(Theta)*(sin(Theta))^2);
TF = 1.0;
if (h==1 & k==1 & l==1)
maxpeak=S*conj(S)*Multiplicity*LP;
end
fprintf(1,'(%d %d %d) %d %.5f', h, k, l, Multiplicity, d);
fprintf(1,' %5.2f', 360.0*Theta/pi);
fprintf(1,' %5.2f %5.2f %5.2f %5.2f',fF,fCa,nfF,nfCa);
fprintf(1,' %7.2f',S);
fprintf(1,' %5.2f',LP);
fprintf(1,' %9.1f',S*conj(S)*Multiplicity*LP);
fprintf(1,' %6.1f\n',100.0*S*conj(S)*Multiplicity*LP/maxpeak);

end