	function [U1,U2,x] = visc_dip(H,xtip,ztip,L,dip,M,t,tR,N,mu,nu,gravity);

%[U1,U2,x] = visc_dip_new(H,xtip,ztip,L,dip,M,t,tR,N,mu);
%Subroutine called by ElasVisco2D
%subsubroutines Pintegrated.m, roots_right.m, makeV_right

xmax = 20*H;       %range of x
delta = 2*xmax/N;  %sampling interval in x
kmax = 0.5/delta;  %maximum wavenumber
x = linspace(-xmax, xmax,N);
ks = linspace(-kmax,kmax,N)*(2*pi);  
% factor of 2*pi accounts for different definition of transform

% some material constants
lam=-2*nu.*mu./(2*nu-1);
K=lam+2*mu/3;   %bulk modulus
g=lam+2*mu;
B = 1/tR;

u1 = zeros(size(ks));
u2 = zeros(size(ks));

 
for j=1:N

	k = ks(j);
    
	%set up A matrix in equation dv=A*v+f
         A=[0  k 1/mu 0; ...
        -k*lam/g 0 0 1/g; ...
        4*k^2*mu*(lam+mu)/g 0 0 k*lam/g; ...
        0 0 -k 0];
  
  rg=3*10^-3*gravity;
  G=[1 0 0 0;0 1 0 0;0 0 1 0;0 -rg 0 1];
  
  
	%propagator from top of halfspace to surface
	z=0; z0=H;
	Ph = expm(A*(z-z0));
	Ph=G*Ph;
	%propagator from source to surface
   Pzs=G*Pintegrated(xtip,ztip,L,dip,lam,mu,ks(j));

	%body forces
	F = [i*M(1,2)/mu; M(2,2)/g; i*i*k*M(1,1) + k*lam*M(2,2)/g; 0];
   
	%compute roots
	[s1, s2, s3, a] = roots_right(k,mu,K,lam,B,Ph);

	%% RETURN VISCOELASTIC PART ONLY
	%% FIRST root:
	V = makeV_right(k,K,lam,s1,B,H,mu);
	%V = makeV_kaj(k,K,s1,B,H,mu);
   A = Ph(3:4,:)*V;
	Adj(1,1) = A(2,2); Adj(2,2) = A(1,1); 
		Adj(1,2) = -A(1,2); Adj(2,1) = -A(2,1);
	Num = -Ph(1:2,:)*V*Adj*Pzs(3:4,:)*F;
	Den = a*s1*(s1-s2)*(s1-s3);
	term1 = Num*(exp(s1*t) -1)/Den;
   
  
	%% SECOND root:
	V = makeV_right(k,K,lam,s2,B,H,mu);
%	V = makeV_kaj(k,K,s2,B,H,mu);
	A = Ph(3:4,:)*V;
	Adj(1,1) = A(2,2); Adj(2,2) = A(1,1); 
		Adj(1,2) = -A(1,2); Adj(2,1) = -A(2,1);
	Num = -Ph(1:2,:)*V*Adj*Pzs(3:4,:)*F;
	Den = a*s2*(s2-s1)*(s2-s3);
	term2 = Num*(exp(s2*t)-1)/Den;

	%% THIRD root:
   V = makeV_right(k,K,lam,s3,B,H,mu);
  % V = makeV_kaj(k,K,s3,B,H,mu);
	A = Ph(3:4,:)*V;
	Adj(1,1) = A(2,2); Adj(2,2) = A(1,1); 
		Adj(1,2) = -A(1,2); Adj(2,1) = -A(2,1);
	Num = -Ph(1:2,:)*V*Adj*Pzs(3:4,:)*F;
	Den =  a*s3*(s3-s1)*(s3-s2);
	term3 = Num*(exp(s3*t)-1)/Den;

	Fhat = term1 + term2 + term3;
   
u1(j) = -i*Fhat(1);
u2(j) = -Fhat(2);
   
end %j (loop over wave number)


%filter out high frequency noise at large k
temp=real(u1);
tempo=mean(abs(temp(end-10:end)));
temp1=zeros(size(temp));temp2=zeros(size(temp));temp3=zeros(size(temp));temp4=zeros(size(temp));temp5=zeros(size(temp));
if tempo>1e-3
   temp=diff(temp,10);
   temp1=temp(1:5:end);temp2=temp(2:5:end);temp3=temp(3:5:end);temp4=temp(4:5:end);temp5=temp(5:5:end);
   temp=temp1(1:N/10)+temp2(1:N/10)+temp3(1:N/10)+temp4(1:N/10)+temp5(1:N/10);
   
   index=abs(temp)<10^-4;
   
  
   
   q=0;
   p=0;
   while q==0
      p=p+1;
      q=index(p);
   end

   
u1(1:5*p)=0;
u1(end-5*p:end)=0;
u2(1:5*p)=0;
u2(end-5*p:end)=0;

end


      %inverse FT into physical space

    %divide by sampling interval to turn analytical transform
    %shift negative frequencies to N/2+1 to N-1

    %inverse FFT
	u1s=ifft(fftshift(u1/delta));
	u2s=ifft(fftshift(u2/delta));
    %shift the solution back
	U1=real(fftshift(u1s));
	U2=real(fftshift(u2s));


  
function P=Pintegrated(xtip,ztip,length,dip,lam,mu,k)
%this calculates the P matrix which has been integrated over the point source position to create
%the solution for a dislocation of finite length -- the algebra is a gawdawful mess

a=dip;
g=lam+2*mu;
L=lam;
R=sqrt(xtip^2+ztip^2);

P(1,1)=1/(4*g*k)*(i*(5*g - 3*L - 4*mu)*cos(a + i*k*R*sin(a)) + (g - 3*L - 4*mu)*(k*R*cos(2*a + i*k*R*sin(a)) - ...
   i*cos(3*a + i*k*R*sin(a))) + k*(-g + 3*L + 4*mu)*R*cosh(k*R*sin(a)))/exp(i*k*R*cos(a));
P(1,2)=1/(4*g*k)*((-7*g + 9*L + 12*mu)*sin(a + i*k*R*sin(a)) + (g - 3*L - 4*mu)*(i*k*R*sin(2*a + i*k*R*sin(a)) + ... 
   sin(3*a + i*k*R*sin(a))) + k*(g - 3*L - 4*mu)*R*sinh(k*R*sin(a)))/exp(i*k*R*cos(a));
P(1,3)=1/(4*g*k^2*mu)*((-7*g + 6*L + 9*mu)*sin(a + i*k*R*sin(a)) + (g - 2*L - 3*mu)*(i*k*R*sin(2*a + i*k*R*sin(a)) + ... 
   sin(3*a + i*k*R*sin(a))) + k*(g - 2*L - 3*mu)*R*sinh(k*R*sin(a)))/exp(i*k*R*cos(a));
P(1,4)=((L + mu)*sin(a)*(k*R*sin(a + i*k*R*sin(a)) - i*sin(2*a + i*k*R*sin(a))))/(exp(i*k*R*cos(a))*2*g*k^2*mu);
P(2,1)=1/(4*g^2*k)*((7*g*L - 3*(L + 2*mu)*(3*L + 2*mu))*sin(a + i*k*R*sin(a)) + (-g*L + (L + 2*mu)*(3*L + 2*mu))* ...
     (i*k*R*sin(2*a + i*k*R*sin(a)) + sin(3*a + i*k*R*sin(a))) + k*(-g*L + (L + 2*mu)*(3*L + 2*mu))*R*sinh(k*R*sin(a)))/ ...
     exp(i*k*R*cos(a));
P(2,2)=(i*(5*g + L)*cos(a + i*k*R*sin(a)) + (g + L)*(k*R*cos(2*a + i*k*R*sin(a)) - i*cos(3*a + i*k*R*sin(a))) - ... 
   k*(g + L)*R*cosh(k*R*sin(a)))/(exp(i*k*R*cos(a))*4*g*k);
P(2,3)=(i*(L + mu)*sin(a)*(i*k*R*sin(a + i*k*R*sin(a)) + sin(2*a + i*k*R*sin(a))))/(exp(i*k*R*cos(a))*2*g*k^2*mu);
P(2,4)=1/(4*g^2*k^2*mu)*(-(7*g*mu + 3*L*(L + 2*mu))*sin(a + i*k*R*sin(a)) + (L^2 + g*mu + 2*L*mu)* ...
   (i*k*R*sin(2*a + i*k*R*sin(a)) + sin(3*a + i*k*R*sin(a))) + k*(L^2 + g*mu + 2*L*mu)*R*sinh(k*R*sin(a)))/exp(i*k*R*cos(a));
P(3,1)=1/g^2*(mu*(L + mu)*((-7*g + 6*L + 12*mu)*sin(a + i*k*R*sin(a)) + (g - 2*L - 4*mu)*(i*k*R*sin(2*a + i*k*R*sin(a)) + ...
   sin(3*a + i*k*R*sin(a))) + k*(g - 2*(L + 2*mu))*R*sinh(k*R*sin(a))))/exp(i*k*R*cos(a));
P(3,2)=(2*mu*(L + mu)*sin(a)*(k*R*sin(a + i*k*R*sin(a)) - i*sin(2*a + i*k*R*sin(a))))/(exp(i*k*R*cos(a))*g);
P(3,3)=1/(4*g*k)*(i*(5*g - 3*L - 4*mu)*cos(a + i*k*R*sin(a)) + (g - 3*L - 4*mu)*(k*R*cos(2*a + i*k*R*sin(a)) - ... 
   i*cos(3*a + i*k*R*sin(a))) + k*(-g + 3*L + 4*mu)*R*cosh(k*R*sin(a)))/exp(i*k*R*cos(a));
P(3,4)=1/(4*g^2*k)*((-7*g*L + 3*(L + 2*mu)*(3*L + 2*mu))*sin(a + i*k*R*sin(a)) + (g*L - (L + 2*mu)*(3*L + 2*mu))* ...
     (i*k*R*sin(2*a + i*k*R*sin(a)) + sin(3*a + i*k*R*sin(a))) + k*(g*L - (L + 2*mu)*(3*L + 2*mu))*R*sinh(k*R*sin(a)))/ ...
     exp(i*k*R*cos(a));
P(4,1)=(2*i*mu*(L + mu)*sin(a)*(i*k*R*sin(a + i*k*R*sin(a)) + sin(2*a + i*k*R*sin(a))))/(exp(i*k*R*cos(a))*g);
P(4,2)=(mu*(L + mu)*(-3*sin(a + i*k*R*sin(a)) + i*k*R*sin(2*a + i*k*R*sin(a)) + sin(3*a + i*k*R*sin(a)) + k*R*sinh(k*R*sin(a))))/ ...
   (exp(i*k*R*cos(a))*g);
P(4,3)=1/(4*g*k)*((7*g - 3*(3*L + 4*mu))*sin(a + i*k*R*sin(a)) + (-g + 3*L + 4*mu)*(i*k*R*sin(2*a + i*k*R*sin(a)) + ... 
   sin(3*a + i*k*R*sin(a))) + k*(-g + 3*L + 4*mu)*R*sinh(k*R*sin(a)))/exp(i*k*R*cos(a));
P(4,4)=(i*(5*g + L)*cos(a + i*k*R*sin(a)) + (g + L)*(k*R*cos(2*a + i*k*R*sin(a)) - i*cos(3*a + i*k*R*sin(a))) - ... 
   k*(g + L)*R*cosh(k*R*sin(a)))/(exp(i*k*R*cos(a))*4*g*k);



R=sqrt((xtip+length*cos(a))^2+(ztip+length*sin(a))^2);
P2(1,1)=1/(4*g*k)*(i*(5*g - 3*L - 4*mu)*cos(a + i*k*R*sin(a)) + (g - 3*L - 4*mu)*(k*R*cos(2*a + i*k*R*sin(a)) - ...
   i*cos(3*a + i*k*R*sin(a))) + k*(-g + 3*L + 4*mu)*R*cosh(k*R*sin(a)))/exp(i*k*R*cos(a));
P2(1,2)=1/(4*g*k)*((-7*g + 9*L + 12*mu)*sin(a + i*k*R*sin(a)) + (g - 3*L - 4*mu)*(i*k*R*sin(2*a + i*k*R*sin(a)) + ... 
   sin(3*a + i*k*R*sin(a))) + k*(g - 3*L - 4*mu)*R*sinh(k*R*sin(a)))/exp(i*k*R*cos(a));
P2(1,3)=1/(4*g*k^2*mu)*((-7*g + 6*L + 9*mu)*sin(a + i*k*R*sin(a)) + (g - 2*L - 3*mu)*(i*k*R*sin(2*a + i*k*R*sin(a)) + ... 
   sin(3*a + i*k*R*sin(a))) + k*(g - 2*L - 3*mu)*R*sinh(k*R*sin(a)))/exp(i*k*R*cos(a));
P2(1,4)=((L + mu)*sin(a)*(k*R*sin(a + i*k*R*sin(a)) - i*sin(2*a + i*k*R*sin(a))))/(exp(i*k*R*cos(a))*2*g*k^2*mu);
P2(2,1)=1/(4*g^2*k)*((7*g*L - 3*(L + 2*mu)*(3*L + 2*mu))*sin(a + i*k*R*sin(a)) + (-g*L + (L + 2*mu)*(3*L + 2*mu))* ...
     (i*k*R*sin(2*a + i*k*R*sin(a)) + sin(3*a + i*k*R*sin(a))) + k*(-g*L + (L + 2*mu)*(3*L + 2*mu))*R*sinh(k*R*sin(a)))/ ...
     exp(i*k*R*cos(a));
P2(2,2)=(i*(5*g + L)*cos(a + i*k*R*sin(a)) + (g + L)*(k*R*cos(2*a + i*k*R*sin(a)) - i*cos(3*a + i*k*R*sin(a))) - ... 
   k*(g + L)*R*cosh(k*R*sin(a)))/(exp(i*k*R*cos(a))*4*g*k);
P2(2,3)=(i*(L + mu)*sin(a)*(i*k*R*sin(a + i*k*R*sin(a)) + sin(2*a + i*k*R*sin(a))))/(exp(i*k*R*cos(a))*2*g*k^2*mu);
P2(2,4)=1/(4*g^2*k^2*mu)*(-(7*g*mu + 3*L*(L + 2*mu))*sin(a + i*k*R*sin(a)) + (L^2 + g*mu + 2*L*mu)* ...
   (i*k*R*sin(2*a + i*k*R*sin(a)) + sin(3*a + i*k*R*sin(a))) + k*(L^2 + g*mu + 2*L*mu)*R*sinh(k*R*sin(a)))/exp(i*k*R*cos(a));
P2(3,1)=1/g^2*(mu*(L + mu)*((-7*g + 6*L + 12*mu)*sin(a + i*k*R*sin(a)) + (g - 2*L - 4*mu)*(i*k*R*sin(2*a + i*k*R*sin(a)) + ...
   sin(3*a + i*k*R*sin(a))) + k*(g - 2*(L + 2*mu))*R*sinh(k*R*sin(a))))/exp(i*k*R*cos(a));
P2(3,2)=(2*mu*(L + mu)*sin(a)*(k*R*sin(a + i*k*R*sin(a)) - i*sin(2*a + i*k*R*sin(a))))/(exp(i*k*R*cos(a))*g);
P2(3,3)=1/(4*g*k)*(i*(5*g - 3*L - 4*mu)*cos(a + i*k*R*sin(a)) + (g - 3*L - 4*mu)*(k*R*cos(2*a + i*k*R*sin(a)) - ... 
   i*cos(3*a + i*k*R*sin(a))) + k*(-g + 3*L + 4*mu)*R*cosh(k*R*sin(a)))/exp(i*k*R*cos(a));
P2(3,4)=1/(4*g^2*k)*((-7*g*L + 3*(L + 2*mu)*(3*L + 2*mu))*sin(a + i*k*R*sin(a)) + (g*L - (L + 2*mu)*(3*L + 2*mu))* ...
     (i*k*R*sin(2*a + i*k*R*sin(a)) + sin(3*a + i*k*R*sin(a))) + k*(g*L - (L + 2*mu)*(3*L + 2*mu))*R*sinh(k*R*sin(a)))/ ...
     exp(i*k*R*cos(a));
P2(4,1)=(2*i*mu*(L + mu)*sin(a)*(i*k*R*sin(a + i*k*R*sin(a)) + sin(2*a + i*k*R*sin(a))))/(exp(i*k*R*cos(a))*g);
P2(4,2)=(mu*(L + mu)*(-3*sin(a + i*k*R*sin(a)) + i*k*R*sin(2*a + i*k*R*sin(a)) + sin(3*a + i*k*R*sin(a)) + k*R*sinh(k*R*sin(a))))/ ...
   (exp(i*k*R*cos(a))*g);
P2(4,3)=1/(4*g*k)*((7*g - 3*(3*L + 4*mu))*sin(a + i*k*R*sin(a)) + (-g + 3*L + 4*mu)*(i*k*R*sin(2*a + i*k*R*sin(a)) + ... 
   sin(3*a + i*k*R*sin(a))) + k*(-g + 3*L + 4*mu)*R*sinh(k*R*sin(a)))/exp(i*k*R*cos(a));
P2(4,4)=(i*(5*g + L)*cos(a + i*k*R*sin(a)) + (g + L)*(k*R*cos(2*a + i*k*R*sin(a)) - i*cos(3*a + i*k*R*sin(a))) - ... 
   k*(g + L)*R*cosh(k*R*sin(a)))/(exp(i*k*R*cos(a))*4*g*k);


P=P2-P;


	function [V] = makeV_right(k,K,L,s,B,H,u);
	

if k < 0

	V = [ -(s+2*B)  		(-2*B - s)*(4*B*K + 2*L*s + 3*s*u + 4*B*H*k*K*s*u + 2*H*k*L*s^2*u + 2*H*k*s^2*u^2); ...
	       (s+2*B)   		(2*B + s)*(2*B*K + L*s + 4*B*H*k*K*s*u + 2*H*k*L*s^2*u + 2*H*k*s^2*u^2);  ...
	      -2*u*s*k       -4*k*s*u*(2*B*K + L*s + s*u)*(1 + H*k*s*u); ...
	       2*u*s*k        2*k*s*u*(2*B*K + L*s + s*u)*(1 + 2*H*k*s*u)];

else
	V = [-(s+2*B)  -s*(2*B + s)*u*(-1 + 4*B*H*k*K + 2*H*k*L*s + 2*H*k*s*u); ...
	     -(s+2*B)  (-2*B - s)*(2*B*K + L*s + 2*s*u + 4*B*H*k*K*s*u + 2*H*k*L*s^2*u + 2*H*k*s^2*u^2); ...
	     2*u*s*k        4*H*k^2*s^2*u^2*(2*B*K + L*s + s*u); ...
	     2*u*s*k        2*k*s*u*(2*B*K + L*s + s*u)*(1 + 2*H*k*s*u)];

end


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);

if k < 0

%neg wave numbers k

anegk3 =  (-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);

anegk2 = (-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);

anegk1 = (-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);

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

b2 = anegk2/anegk3;
b1 = anegk1/anegk3;
b0 = anegk0/anegk3;
a = anegk3;

else
%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;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%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;



