function y=pow_slopes_lith(rx_effect,sample_size,ratio,power,visit,vs,sigma2,drop_month,alpha)
%pow_slopes(rx_effect,sample_size,ratio,power,visit,vs,sigma2,drop_month,alpha);
%Finds power for two fixed sample sizes in ratio;
%rx_effect: difference in slopes
%visits: cell array giving first treatment visit schedule then control visit schedule
%ratio: a 1 x 2 matrix which is the randomization ratio [2 1] implies 2 treatment to 1 control
%vs: variance covariance matrix of random effects
%sigma2: Error variance;
%drop_month: If it is a single number than it is the drop out rate per month i.e. 0.02;
%if it is a cell array with two vectors then it is the probability of having exactly  that number of visits; 
%alpha: One sided alpha level of test
%y:  Treatment effect or sample_size or power depending on what has been set to zero.
%EXAMPLE  ALFRS
%vs=[35.4864 0.6848;0.6848 .5540];
%vs=[35.4864 0.7031;0.7031 .5724];
%sigma2: 2.9739;
%pow_slopes(0,40,[1 1],.8,{[0 1 2 4 6],[0 1 2 4 6]},[31.8 .8527;.8527 .6687],2.7085,.02,.05)
%ans=0.7173
vsv=zeros(3,3);
s2b=vs(2,2);
s2e=sigma2;
dflag=isa(drop_month,'numeric');
for i=1:2
    visits=visit{i};
    for j=2:length(visits)
        if j<length(visits)
            %Proportion of patients with only j visits;
            %pp=prob_drop(-1+rx_effect,s2b,s2e,visits(2:(j+1)))-prob_drop(-1+rx_effect,s2b,s2e,visits(2:j));
            if dflag 
                 pp=exp(-visits(j)*drop_month)-exp(-visits(j+1)*drop_month);
            else
                pp=drop_month{i}(j);
            end
        else
            if dflag
                  pp=exp(-visits(j)*drop_month);
            else
                pp=drop_month{i}(j);
            end
            %pp=1-prob_drop(-1+rx_effect,s2b,s2e,visits(2:j));
        end
        vis=visits(1:j);
        sz=ratio(i)*length(vis);
        ssz=length(vis);
        x=[ones(sz,1) (kron(ones(ratio(i),1),vis')*ones(1,2)).*[ones(sz,1) (2-i)*ones(sz,1)]];
        z=kron(diag(ones(ratio(i),1)), [ones(ssz,1) vis']);
        v=x'*inv((sigma2*diag(ones(sz,1))+z*kron(diag(ones(ratio(i),1)),vs)*z'))*x;
        vsv=vsv+v*pp;

    end
end
vsvt=inv(vsv);
sdb=sqrt(vsvt(3,3));
mult=sample_size/sum(ratio);
if power==0
    y=normcdf(sqrt(mult)*rx_effect/sdb-norminv(1-alpha));
else if rx_effect==0
        y=(norminv(power)+norminv(1-alpha))*sdb/sqrt(mult);
    else
        y=sum(ratio)*((norminv(power)+norminv(1-alpha))*sdb/rx_effect)^2;
    end
end
