% Conformance probability pc
pc_v = [0.5];
%pc_v = [0.95];
n_pc = length(pc_v);
% MPE
MPE = 0.2;
% Measurement Capability Cm as defined in JCGM 106:2012
Cm_v = [1;2;4;6;10];
n_Cm = length(Cm_v);
% two FOR cycles one for pc (k index) and the other for Cm (q index)
figure
for q = 1 : n_Cm
    Cm = Cm_v(q,1);
    for k = 1 : n_pc
        pc = pc_v(k,1);
        % test standard uncertanty ut
        ut = MPE/(2*Cm);
        % definition of small quantity epsilon eps
        eps = 10^-4;
        % difference between pc computed and the pc value set
        % delta_pc = pc_calc - pc
        delta_pc = 10;
        % until the difference of calculated pc (pc_calc)and the set pc (delta_pc) is greater than eps
        % keep moving the point y -> see WHILE cylce
        pd = makedist('Normal');
        y = 0;
        step_y = MPE/(100*ut);
        if step_y > 10^-4
            step_y = 10^-4;
        end
        while delta_pc > eps
            % compute cumulative probability
            p1 = cdf(pd,(MPE-y)/ut);
            p2 = cdf(pd,(-MPE-y)/ut);
            pc_calc = p1 - p2;
            delta_pc = pc_calc - pc;
            z = (MPE-y)/ut;
            y = y + step_y;
            % record pc_calc values in a matrix PC_calc (k x q)
            PC_calc(k,q)= pc_calc;
            % record z values in a matrix
            Z_calc(k,q)= z;
        end
        z
        % in order to define the x axis, we consider the limits equal to MPE+3*ut and -MPE-3*ut
        % dividing the range [-MPE-3*ut, MPE+3*ut] by the number n1
        % for each point in the range we consider a normal distribution centered in x
        n1 = 1000;
        x = linspace(-MPE-3*ut, MPE+3*ut, n1);
        % we define the vector of mean values mu of the normal distributions. The limits of this vector are the limits of the acceptance
        % zone and depend on the pc value set
        n2 = 1000;
        mu = linspace(-MPE+z*ut, MPE-z*ut, n2+1);
        % we define the guard band value g
        g = MPE - z*ut
        % we determine the normal distribution for each  value "mu" (n2+1)
        for w = 1 : n2+1
            p(w,:) = normpdf(x, mu(1,w),ut);
            X(w,:) = x;
        end
        % sum the normal distributions
        PDF = sum(p);
        % normalization of the resulting distribution PDF (sum of normal
        % distributions)
        norm = sum(PDF);
        norm_1 = norm * (x(n1)-x(1))/n1;
        PDF_norm = PDF/norm_1;                 
        plot(x,PDF_norm)
        hold on  
         
    end
    Ut(q,1) = ut
end


