function [z1, z2, z3, a] = roots_right(k,u,K,L,B,Ph);
                                     
P31=Ph(3,1);P32=Ph(3,2);P33=Ph(3,3);P34=Ph(3,4);
P41=Ph(4,1);P42=Ph(4,2);P43=Ph(4,3);P44=Ph(4,4);

%positive wave numbers k
aposk3 = (-L*P32*P41 + L*P31*P42 - 3*P32*P41*u + ...
    2*k*L*P34*P41*u + 3*P31*P42*u - 2*k*L*P33*P42*u + 2*k*L*P32*P43*u - ...
    2*k*L*P31*P44*u + 2*k*P33*P41*u^2 + 4*k*P34*P41*u^2 - 4*k*P33*P42*u^2 - ...
    2*k*P34*P42*u^2 - 2*k*P31*P43*u^2 + 4*k*P32*P43*u^2 - 4*k^2*L*P34*P43*u^2 - ...
    4*k*P31*P44*u^2 + 2*k*P32*P44*u^2 + 4*k^2*L*P33*P44*u^2 - ...
    4*k^2*P34*P43*u^3 + 4*k^2*P33*P44*u^3);

aposk2 = (-2*B*K*P32*P41 - 4*B*L*P32*P41 + 2*B*K*P31*P42 + 4*B*L*P31*P42 - ...
    12*B*P32*P41*u + 4*B*k*K*P34*P41*u + 4*B*k*L*P34*P41*u + 12*B*P31*P42*u - ...
    4*B*k*K*P33*P42*u - 4*B*k*L*P33*P42*u + 4*B*k*K*P32*P43*u + ...
    4*B*k*L*P32*P43*u - 4*B*k*K*P31*P44*u - 4*B*k*L*P31*P44*u + ...
    4*B*k*P33*P41*u^2 + 8*B*k*P34*P41*u^2 - 8*B*k*P33*P42*u^2 - ...
    4*B*k*P34*P42*u^2 - 4*B*k*P31*P43*u^2 + 8*B*k*P32*P43*u^2 - ...
    8*B*k^2*K*P34*P43*u^2 - 8*B*k*P31*P44*u^2 + 4*B*k*P32*P44*u^2 + ...
    8*B*k^2*K*P33*P44*u^2);

aposk1 = (-8*B^2*K*P32*P41 - 4*B^2*L*P32*P41 + ...
    8*B^2*K*P31*P42 + 4*B^2*L*P31*P42 - 12*B^2*P32*P41*u + ....
    8*B^2*k*K*P34*P41*u + 12*B^2*P31*P42*u - 8*B^2*k*K*P33*P42*u + ...
    8*B^2*k*K*P32*P43*u - 8*B^2*k*K*P31*P44*u);

aposk0 = -8*B^3*K*P32*P41 + 8*B^3*K*P31*P42;

b2 = aposk2/aposk3;
b1 = aposk1/aposk3;
b0 = aposk0/aposk3;
a = aposk3;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Solve CUBIC


%% associatd cubic

	q = b1/3 - b2^2/9;
	r = (b1*b2-3*b0)/6 - b2^3/27;
	rt = q^3 + r^2;
	s1 = (r + sqrt(rt))^(1/3);
	s2 = (r - sqrt(rt))^(1/3);
	z1 = s1 + s2 -b2/3;
	z2 = -(s1 + s2)/2 - b2/3 +i*sqrt(3)*(s1 - s2)/2;
	z3 = -(s1 + s2)/2 - b2/3 -i*sqrt(3)*(s1 - s2)/2;

