Basic support for experimental scattering intensities determination

Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>

zUtil_Cryst/gtCrystComputeExperimentalFamilyIntensities.m

+function fams_ints = gtCrystComputeExperimentalFamilyIntensities(phase_id, p)
+    if (~exist('phase_id', 'var'))
+        phase_id = 1;
+    end
+    if (~exist('p', 'var'))
+        p = gtLoadParameters();
+    end
+
+    sample = GtSample.loadFromFile();
+
+    num_grains = sample.phases{phase_id}.getNumberOfGrains();
+    num_grains = max(num_grains - 100, round(num_grains * 0.9));
+
+    fams_ints = zeros(size(p.cryst(phase_id).int'));
+    fams_counts = zeros(size(p.cryst(phase_id).int'));
+
+    fprintf('Loading info: ')
+    c = tic();
+    for ii_g = 1:num_grains
+        num_chars = fprintf('%03d/%03d', ii_g, num_grains);
+
+        gr = gtLoadGrain(phase_id, ii_g);
+        % Working around wrong formula from Henning's book
+        gr = gtCalculateGrain(gr, p);
+
+%         gr_rec = gtLoadGrain(phase_id, ii_g);
+
+        if (isfield(gr.proj, 'ondet'))
+            included = gr.proj.ondet(gr.proj.included);
+        else
+            included = gr.ondet(gr.included);
+        end
+%         try
+%             num_voxels = numel(gr_rec.SEG.seg);
+%         catch
+%             num_voxels = numel(gr.seg);
+%         end
+
+        not_nan = ~isnan(gr.intensity);
+        inc_good = included(not_nan);
+        L_fac = gr.allblobs.lorentzfac(inc_good);
+        theta_types = gr.allblobs.thetatype(inc_good);
+
+        ints = gr.intensity(not_nan) ./ L_fac;
+%         ints = gr.intensity(not_nan) ./ L_fac / num_voxels;
+        fams_ints = fams_ints + accumarray(theta_types, ints, size(fams_ints));
+        fams_counts = fams_counts + accumarray(theta_types, ones(size(theta_types)), size(fams_ints));
+
+        fprintf(repmat('\b', [1 num_chars]));
+    end
+    fprintf('Done in %f seconds.\n', toc(c))
+
+    usable = fams_counts > 0;
+    fams_ints = fams_ints(usable) ./ fams_counts(usable);
+    fams_ints = fams_ints' / norm(fams_ints);
+
+    fprintf('Theoretical Intensities:\n');
+    disp(p.cryst(phase_id).int(usable) / norm(p.cryst(phase_id).int(usable)));
+    fprintf('Calculated Intensities:\n');
+    disp(fams_ints);
+end