function [res]=assessplan(k,m,n,b,sp2,so2,sd2,sm2)
%assess the plan SP(k,m,n) AP A (b/m,m)
%give results for both gamma and total measurement variation, i.e. so2+sd2+sm2
%if sd2=0 this simplifies to the no interaction model case

%b equals the total number of parts in baseline!

%two options for calling function: provide sp2,so2,sd2,sm2 in that order
%OR provide gamma,soratio,smratio in that order, where all must be between
%0 and 1 and we assume sp2+so2+sd2+sm2=1
%and gamma is sqrt((so2+sd2+sm2)/(sp2+so2+sd2+sm2))

if nargin < 8,   %provided gamma,soratio,smratio instead of sp2,so2,sd2,sm2
    gamma=sp2;  soratio=so2;  smratio=sd2;
    sm2=smratio*gamma^2;
    so2=soratio*(1-smratio)*gamma^2;
    sd2=(1-soratio)*(1-smratio)*gamma^2;   %interaction
    sp2=1-so2-sd2-sm2;   %assume sum is 1
end;

res=[];  %save the results [k,m,n,u se(gamma)]
%k = # parts
%m = # of fixed operators
%n = # measurements on each part
%b = # parts measured once in AP part of plan, b must be a multiple of m
if b/m~=floor(b/m),  disp('choose b as a multiple of number of operators'), end;

N=k*m*n+b;   %total number of measurements

%find feasible values of the mus that match so
%assume close to balance for the relative biases among the operators
muj=zeros(m,1);
if m>1,  %i.e. there is more than a single operator
    if (m/2)==round(m/2),  %m even
        mu=sqrt(so2);
    else
        mu=sqrt(m*so2/(m-1));
    end;
    muj(1:floor(m/2))=mu;
    muj(ceil(m/2)+1:end)=-mu;
end;

if m==1 && so2>0, disp('WARNING, when there is only one operator the so2 value is ignored, it is set to zero'), end;

%find change of variables matrix for gamma
D=cmatrix(muj,sp2,sd2,sm2);
%find change of variables matrix for sigma_meas = so2+sd2+so2
D2=cmatrix2(muj,sp2,sd2,sm2);

%type A AP
Jmat = J_kmn(k,m,n,sp2,sd2,sm2) + J_typeA(m,b/m,sp2,sd2,sm2);  %add the two information matrices       
varmat=inv(Jmat);  %diagonals give variance for sp2, sd2 and sm2
SEsp=sqrt(varmat(m+1,m+1))/(2*sqrt(sp2));      %determine the SE for sp  (not sp2)
temp=D*inv(Jmat)*D';
temp2=inv(D2*Jmat*D2');   %for total measurement variance
 %temp=inv(Jfinal);   %diagonals give SE for sp2, sd2 and sm2
if sd2~=0,   res=[res;k,m,n,b,SEsp,sqrt(temp(2,2)),sqrt(temp(3,3)),sqrt(temp(4,4)),sqrt(temp(1,1)),sqrt(temp2(1,1))];
else
    if m>=2,  res=[res;k,m,n,b,SEsp,sqrt(temp(2,2)),sqrt(temp(3,3)),sqrt(temp(1,1)),sqrt(temp2(1,1))]; 
    else  res=[res;k,m,n,b,SEsp,sqrt(temp(2,2)),sqrt(temp(1,1)),sqrt(temp2(1,1))]; end;
end;
         


[temp,ind]=sort(res(:,end));
res=res(ind,:);

gam=sqrt((so2+sd2+sm2)/(sp2+so2+sd2+sm2));

dispres=1;
if dispres==1,
    %display results
    disp(['number of operators = ',num2str(m),'  total number of measurements = ',num2str(N),'  gamma = ',num2str(gam)])
    disp(['sigmap2 = ',num2str(sp2),'  sigmao2 = ',num2str(so2),'  sigmad2 = ',num2str(sd2),'  sigmam2 = ',num2str(sm2)])
    if sd2~=0,   disp('    #parts     #oper     #meas #baseline      SEsp      SEso      SEsd      SEsm   SEgamma  SE_Mtotal')
    else
        if m>=2,  disp('    #parts     #oper     #meas #baseline      SEsp      SEso      SEsm   SEgamma  SE_Mtotal') 
        else  disp('    #parts     #oper     #meas #baseline      SEsp      SEsm   SEgamma   SE_Mtotal'),   end;
    end
    disp(res)
end

function [mat]=cmatrix(mu,sp2,sd2,sm2)
%returns the change of variables matrix based on the partial derivatives
%go from muis sp2, sd2 and sm2 (note variance scale) to gamma, so, sd, sm in that order
%skip row and column for sd if we assume the no interaction model
%also remove the so column if there is only one operat, i.e. m=1
%for fixed operator case!

m=length(mu);  %number of operators
if sd2~=0,   %interaction model, this doesn't make sense if there is only one operator
    mat=zeros(4,m+3);   %four rows
    if m==1, disp('WARNING, using the interaction model when there is only one operator does not make sense, select sd2=0'), end;
    mat(3,m+2)=1/(2*sqrt(sd2));  %only element in third row for sd
    mat(4,m+3)=1/(2*sqrt(sm2));  %only element in fourth row for sm
else    %no interaction model
    if m>=2
       mat=zeros(3,m+2);   %three rows
       mat(3,m+2)=1/(2*sqrt(sm2));  %only element in third row for sm
    else  %only one operator case
       mat=zeros(2,m+2);   %two rows, only for gamma and sm
       mat(2,m+2)=1/(2*sqrt(sm2));  %only element in second row for sm
    end;
end;

muiminusmubar=mu-mean(mu);
so=sqrt(sum(muiminusmubar.^2)/m);
if m>=2,  mat(2,1:m)=muiminusmubar/(m*so);   end; %for so

%can't include sp in change of variable matrix, otherwise there is a redundancy
%mat(2,m+1)=1/(2*sqrt(sp2));   %only element in second row for sp

%first row for gamma
gam=sqrt((so^2+sd2+sm2)/(sp2+so^2+sd2+sm2));
mat(1,1:m)=muiminusmubar*(1-gam^2)/gam/(sp2+so^2+sd2+sm2)/m;
so2=so^2;  sp=sqrt(sp2); sd=sqrt(sd2);  sm=sqrt(sm2);

%change of variables on the variance scale
mat(1,m+1) = -((so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2)) ^ (-0.1e1 / 0.2e1) * (so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2) ^ 2 / 0.2e1;
if sd2~=0, 
   mat(1,m+2) = ((so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2)) ^ (-0.1e1 / 0.2e1) * (0.1e1 / (sp2 + so2 + sd2 + sm2) - (so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2) ^ 2) / 0.2e1;
   mat(1,m+3) = ((so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2)) ^ (-0.1e1 / 0.2e1) * (0.1e1 / (sp2 + so2 + sd2 + sm2) - (so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2) ^ 2) / 0.2e1;
else
   mat(1,m+2) = ((so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2)) ^ (-0.1e1 / 0.2e1) * (0.1e1 / (sp2 + so2 + sd2 + sm2) - (so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2) ^ 2) / 0.2e1;
end;

function [mat]=cmatrix2(mu,sp2,sd2,sm2)
%returns the change of variables matrix based on the partial derivatives
%go from muis sp2, sd2 and sm2 (note variance scale) to sigma_meas
%skip row and column for sd if we assume the no interaction model
%also remove the so column if there is only one operat, i.e. m=1

m=length(mu);  %number of operators
if sd2~=0,   %interaction model, this doesn't make sense if there is only one operator
    mat=zeros(1,m+3);   
    if m==1, disp('WARNING, using the interaction model when there is only one operator does not make sense, select sd2=0'), end;
else    %no interaction model
    mat=zeros(1,m+2);   
end;

muiminusmubar=mu-mean(mu);
mat(1,1:m)=2*muiminusmubar/m;  %derivatives wrt muj

%change of variables on the variance scale
mat(1,m+1) = 0;   %wrt sp2
if sd2~=0, 
   mat(1,m+2) = 1;
   mat(1,m+3) = 1;
else
   mat(1,m+2) = 1;
end;

function [Jmat]=J_kmn(k,m,n,sp2,sd2,sm2)
%find the information matrix for the SP given k (parts), m (fixed operators), n (repeats)

if sd2 ~=0,  Jmat=zeros(m+3,m+3);   else   Jmat=zeros(m+2,m+2);   end;

%Note that the information matrix Jmat is always block diagonal with terms
%related to the mujs and another block for the rest

b1 = 1/sm2;
b2 = -sd2 / (sm2 * (sm2 + n * sd2));
b3 = -sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2);

%setup top bit related to the mus
Jmat(1:m,1:m)=-k*n*n*b3;
for i=1:m,   %replace the diagonal
    Jmat(i,i)=-k*n*(b1+n*(b2+b3));
end;

%expected vlaues of the various sums of squares
A=k*m*n*(sp2+sd2+sm2);
B=k*m*n*(n*sp2+n*sd2+sm2);
C=k*m*n*(m*n*sp2+n*sd2+sm2);

%now look at bottom right 3x3 part related to sp2, sd2 and sm2 in that order!
if sd2~=0, 
    Jmat(m+1,m+1) = -0.1e1 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * m * n + sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C * m ^ 2 * n ^ 2 + k * m ^ 2 * n ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 / 0.2e1;
    Jmat(m+1,m+2) = -0.1e1 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) * C * n / 0.2e1 - 0.1e1 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * n / 0.2e1 + sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * m * n ^ 2 / 0.2e1 + sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C * m * n ^ 2 + k * m * n ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 / 0.2e1;
    Jmat(m+2,m+1) = Jmat(m+1,m+2);
    Jmat(m+1,m+3) = -0.1e1 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) * C / 0.2e1 - 0.1e1 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C / 0.2e1 + sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * m * n / 0.2e1 + sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C * m * n + k * m * n / (sm2 + n * sd2 + m * n * sp2) ^ 2 / 0.2e1;
    Jmat(m+3,m+1) = Jmat(m+1,m+3);
    Jmat(m+2,m+2) = -(1 / sm2 / (sm2 + n * sd2) ^ 2 * B * n) + (sd2 / sm2 / (sm2 + n * sd2) ^ 3 * B * n ^ 2) + (sp2 / (sm2 + n * sd2) ^ 3 / (sm2 + n * sd2 + m * n * sp2) * C * n ^ 2) + (sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * n ^ 2) + (sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C * n ^ 2) - (k * (2 * n ^ 2 * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) ^ 2 * n ^ 2 / (sm2 + n * sd2) ^ 2 * sm2 ^ (m * n - m) - (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) * n ^ 2 / (sm2 + n * sd2) ^ 2 * sm2 ^ (m * n - m)) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m)) / 0.2e1 + (k * (n * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) * n / (sm2 + n * sd2) * sm2 ^ (m * n - m)) / (sm2 + n * sd2 + m * n * sp2) ^ 2 / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m) * n) / 0.2e1 + (k * (n * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) * n / (sm2 + n * sd2) * sm2 ^ (m * n - m)) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m) * (m - 1) * n / (sm2 + n * sd2)) / 0.2e1;
    Jmat(m+2,m+3) = -0.1e1 / sm2 ^ 2 / (sm2 + n * sd2) * B / 0.2e1 - 0.1e1 / sm2 / (sm2 + n * sd2) ^ 2 * B / 0.2e1 + sd2 / sm2 ^ 2 / (sm2 + n * sd2) ^ 2 * B * n / 0.2e1 + sd2 / sm2 / (sm2 + n * sd2) ^ 3 * B * n + sp2 / (sm2 + n * sd2) ^ 3 / (sm2 + n * sd2 + m * n * sp2) * C * n + sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * n + sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C * n - k * (0.2e1 * n * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + n * (sm2 + n * sd2) ^ (m - 0.1e1) * sm2 ^ (m * n - m) * (m * n - m) / sm2 + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) ^ 2 / (sm2 + n * sd2) ^ 2 * n * sm2 ^ (m * n - m) - (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) * n / (sm2 + n * sd2) ^ 2 * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) * n / (sm2 + n * sd2) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 0.1e1) / sm2 ^ (m * n - m) / 0.2e1 + k * (n * (sm2 + n * sd2) ^ (m - 0.1e1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) * n / (sm2 + n * sd2) * sm2 ^ (m * n - m)) / (sm2 + n * sd2 + m * n * sp2) ^ 2 / (sm2 + n * sd2) ^ (m - 0.1e1) / sm2 ^ (m * n - m) / 0.2e1 + k * (n * (sm2 + n * sd2) ^ (m - 0.1e1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) * n / (sm2 + n * sd2) * sm2 ^ (m * n - m)) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 0.1e1) / sm2 ^ (m * n - m) * (m - 0.1e1) / (sm2 + n * sd2) / 0.2e1 + k * (n * (sm2 + n * sd2) ^ (m - 0.1e1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) * n / (sm2 + n * sd2) * sm2 ^ (m * n - m)) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 0.1e1) / sm2 ^ (m * n - m) * (m * n - m) / sm2 / 0.2e1;
    Jmat(m+3,m+2) = Jmat(m+2,m+3);
    Jmat(m+3,m+3) = -(1 / sm2 ^ 3 * A) + (sd2 / sm2 ^ 3 / (sm2 + n * sd2) * B) + (sd2 / sm2 ^ 2 / (sm2 + n * sd2) ^ 2 * B) + (sd2 / sm2 / (sm2 + n * sd2) ^ 3 * B) + (sp2 / (sm2 + n * sd2) ^ 3 / (sm2 + n * sd2 + m * n * sp2) * C) + (sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C) + (sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C) - (k * (2 * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + 2 * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2 + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) ^ 2 / (sm2 + n * sd2) ^ 2 * sm2 ^ (m * n - m) - (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) ^ 2 * sm2 ^ (m * n - m) + 2 * (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) * (m * n - m) / sm2 + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) ^ 2 / sm2 ^ 2 - (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2 ^ 2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m)) / 0.2e1 + (k * ((sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m)) / 0.2e1 + (k * ((sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m) * (m - 1) / (sm2 + n * sd2)) / 0.2e1 + (k * ((sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m) * (m * n - m) / sm2) / 0.2e1;
else  %remove bits related to sd2
    Jmat(m+1,m+1) = -0.1e1 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * m * n + sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C * m ^ 2 * n ^ 2 + k * m ^ 2 * n ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 / 0.2e1;
    Jmat(m+1,m+2) = -0.1e1 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) * C / 0.2e1 - 0.1e1 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C / 0.2e1 + sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * m * n / 0.2e1 + sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C * m * n + k * m * n / (sm2 + n * sd2 + m * n * sp2) ^ 2 / 0.2e1;
    Jmat(m+2,m+1) = Jmat(m+1,m+2);
    Jmat(m+2,m+2) = -(1 / sm2 ^ 3 * A) + (sd2 / sm2 ^ 3 / (sm2 + n * sd2) * B) + (sd2 / sm2 ^ 2 / (sm2 + n * sd2) ^ 2 * B) + (sd2 / sm2 / (sm2 + n * sd2) ^ 3 * B) + (sp2 / (sm2 + n * sd2) ^ 3 / (sm2 + n * sd2 + m * n * sp2) * C) + (sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C) + (sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C) - (k * (2 * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + 2 * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2 + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) ^ 2 / (sm2 + n * sd2) ^ 2 * sm2 ^ (m * n - m) - (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) ^ 2 * sm2 ^ (m * n - m) + 2 * (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) * (m * n - m) / sm2 + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) ^ 2 / sm2 ^ 2 - (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2 ^ 2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m)) / 0.2e1 + (k * ((sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m)) / 0.2e1 + (k * ((sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m) * (m - 1) / (sm2 + n * sd2)) / 0.2e1 + (k * ((sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m) * (m * n - m) / sm2) / 0.2e1;
end;
%change the sign since we want the negative of the expected value instead of the expected value
Jmat=-Jmat;


function [Jmat]=J_typeA(m,u,sp2,sd2,sm2)
%find the special extra information matrix for the type A AP given m and u
%with AP plan A m operators each measure u *different* parts once each
%total number of parts in extra AP bit is m*u

if sd2 ~=0,  Jmat=zeros(m+3,m+3);
else   Jmat=zeros(m+2,m+2);   end;

%setup top bit related to the mus
Jmat(1:m,1:m)=-u/(sp2+sd2+sm2)*eye(m);  

%n=1; k=u;   %special case
%now look at bottom right 3x3 part related to sp2, sd2, sm2
A=m*u*(sp2+sd2+sm2);
ka=m*u;
if sd2~=0, 
    Jmat(m+1:m+3,m+1:m+3) = -A / (sp2 + sd2 + sm2) ^ 3 + ka / (sp2 + sd2 + sm2) ^ 2 / 0.2e1;   %3x3 matrix all the same
    %off diagonal elements also all the same
    Jmat(m+1,m+2)=-A / (sp2 + sd2 + sm2) ^ 3 + ka / (sp2 + sd2 + sm2) ^ 2 / 0.2e1;
    Jmat(m+1,m+3)=Jmat(m+1,m+2);
    Jmat(m+2,m+3)=Jmat(m+1,m+2);
    Jmat(m+2,m+1)=Jmat(m+1,m+2);
    Jmat(m+3,m+2)=Jmat(m+1,m+2);
    Jmat(m+3,m+1)=Jmat(m+1,m+2);
else  %remove bits related to sd2
    Jmat(m+1:m+2,m+1:m+2) = -A / (sp2 + sd2 + sm2) ^ 3 + ka / (sp2 + sd2 + sm2) ^ 2 / 0.2e1;   %3x3 matrix all the same
    %off diagonal elements also all the same
    Jmat(m+1,m+2)=-A / (sp2 + sd2 + sm2) ^ 3 + ka / (sp2 + sd2 + sm2) ^ 2 / 0.2e1;
    Jmat(m+2,m+1)=Jmat(m+1,m+2);
end;
%change the sign since we want the negative of the expected value instead of the expected value
Jmat=-Jmat;