% Conformance probability pc
pc_v = [0.99; 0.975; 0.95; 0.90; 0.85; 0.80; 0.7; 0.5];
n_pc = length(pc_v);
% MPE
MPE = 0.2;
% Measurement Capability Cm as defined in JCGM 106:2012
Cm_v = [0.98;1;1.1;1.2;1.3;1.4;1.5;2;2.1;2.2;2.3;2.4;2.5;2.6;2.7;2.8;2.9;3;4;5;6;10];
n_Cm = length(Cm_v);
% two FOR cycles one for pc (k index) and the other for Cm (q index)
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);
        Freq = PDF/norm; % relative frequency
        % compute mean and sigma of the PDF
        mean = sum(Freq.*x);
        u = sqrt(sum(Freq.*(x.^2)));
        u_v(q,k) = u;
               
    end
    Ut(q,1) = ut
end
u_MPE = MPE/sqrt(3);
y_MPE = repmat(u_MPE,n_Cm,n_pc);
u_v_MPE = u_v./y_MPE;
figure
x_Cm = Cm_v';

% check the PC_calc values vs pc_v deleting points whith PC_calc < pc_v
for k = 1 : n_pc
    c = pc_v(k)-(2*10^-3);
    a = find(PC_calc(k,:) < c);
    u_v_MPE(a,k)= NaN
end

figure
for r = 1 : n_pc
    y = u_v_MPE(:,r)';
    hold on
    plot(x_Cm,y)
end

 
