function [ux,uy,uz]=DislocOkada(L,W,D,dip,U,lambda,mu,x,y,z)
% 
% D is depth (positive) to bottom of dislocation (c in Okada, 1992)
% W is width (along down-dip direction)
% L is length (along strike)
% U(1) is strike slip
% U(2) is dip slip
% U(3) is tensile slip
% dip is in degrees
% x,y,z is location where deflection is evaluated
% lambda and mu are the elasticity parameters
%
% if slip vector is not horizontal, vertical or tensile use a combination
% of these Okadas to get the total displacement field
%
% Note that these equations provide deflection at any point inside the
% half-space, and the z-axis is positive upwards.  Thus specify z=0 for
% points on the surface, and z negative for points below the surface.
%

alpha = (lambda + mu)/(lambda + 2*mu);

dip = dip*pi/180;
d = D - z;  % yes, d>D for z<0 (beneath the surface)
p = y*cos(dip) + d*sin(dip);
q = y*sin(dip) - d*cos(dip);

ux=0;
uy=0;
uz=0;

%form: [f]=fOkada(xi,eta,z,q,dip,alpha,typ,term,i);

for ityp = 1:3
    if ityp==1
        typ='strike';
    elseif ityp==2
        typ='dip';
    elseif ityp==3
        typ='tensile';
    end
    
    if U(ityp)~=0
        
        u1A = fOkada(x  ,p  ,z,q,dip,alpha,typ,'A',1) - ...
              fOkada(x  ,p-W,z,q,dip,alpha,typ,'A',1) - ...
              fOkada(x-L,p  ,z,q,dip,alpha,typ,'A',1) + ...
              fOkada(x-L,p-W,z,q,dip,alpha,typ,'A',1);

        u1Ahat = fOkada(x  ,p  ,-z,q,dip,alpha,typ,'A',1) - ...
              fOkada(x  ,p-W,-z,q,dip,alpha,typ,'A',1) - ...
              fOkada(x-L,p  ,-z,q,dip,alpha,typ,'A',1) + ...
              fOkada(x-L,p-W,-z,q,dip,alpha,typ,'A',1);

        u1B = fOkada(x  ,p  ,z,q,dip,alpha,typ,'B',1) - ...
              fOkada(x  ,p-W,z,q,dip,alpha,typ,'B',1) - ...
              fOkada(x-L,p  ,z,q,dip,alpha,typ,'B',1) + ...
              fOkada(x-L,p-W,z,q,dip,alpha,typ,'B',1);

        u1C = fOkada(x  ,p  ,z,q,dip,alpha,typ,'C',1) - ...
              fOkada(x  ,p-W,z,q,dip,alpha,typ,'C',1) - ...
              fOkada(x-L,p  ,z,q,dip,alpha,typ,'C',1) + ...
              fOkada(x-L,p-W,z,q,dip,alpha,typ,'C',1);
        
        u2A = fOkada(x  ,p  ,z,q,dip,alpha,typ,'A',2) - ...
              fOkada(x  ,p-W,z,q,dip,alpha,typ,'A',2) - ...
              fOkada(x-L,p  ,z,q,dip,alpha,typ,'A',2) + ...
              fOkada(x-L,p-W,z,q,dip,alpha,typ,'A',2);            
 
        u2Ahat = fOkada(x  ,p  ,-z,q,dip,alpha,typ,'A',2) - ...
                 fOkada(x  ,p-W,-z,q,dip,alpha,typ,'A',2) - ...
                 fOkada(x-L,p  ,-z,q,dip,alpha,typ,'A',2) + ...
                 fOkada(x-L,p-W,-z,q,dip,alpha,typ,'A',2); 
          
        u2B = fOkada(x  ,p  ,z,q,dip,alpha,typ,'B',2) - ...
              fOkada(x  ,p-W,z,q,dip,alpha,typ,'B',2) - ...
              fOkada(x-L,p  ,z,q,dip,alpha,typ,'B',2) + ...
              fOkada(x-L,p-W,z,q,dip,alpha,typ,'B',2); 
        
        u2C = fOkada(x  ,p  ,z,q,dip,alpha,typ,'C',2) - ...
              fOkada(x  ,p-W,z,q,dip,alpha,typ,'C',2) - ...
              fOkada(x-L,p  ,z,q,dip,alpha,typ,'C',2) + ...
              fOkada(x-L,p-W,z,q,dip,alpha,typ,'C',2); 
        
        u3A = fOkada(x  ,p  ,z,q,dip,alpha,typ,'A',3) - ...
              fOkada(x  ,p-W,z,q,dip,alpha,typ,'A',3) - ...
              fOkada(x-L,p  ,z,q,dip,alpha,typ,'A',3) + ...
              fOkada(x-L,p-W,z,q,dip,alpha,typ,'A',3);
          
        u3Ahat = fOkada(x  ,p  ,-z,q,dip,alpha,typ,'A',3) - ...
                 fOkada(x  ,p-W,-z,q,dip,alpha,typ,'A',3) - ...
                 fOkada(x-L,p  ,-z,q,dip,alpha,typ,'A',3) + ...
                 fOkada(x-L,p-W,-z,q,dip,alpha,typ,'A',3); 
                    
        u3B = fOkada(x  ,p  ,z,q,dip,alpha,typ,'B',3) - ...
              fOkada(x  ,p-W,z,q,dip,alpha,typ,'B',3) - ...
              fOkada(x-L,p  ,z,q,dip,alpha,typ,'B',3) + ...
              fOkada(x-L,p-W,z,q,dip,alpha,typ,'B',3); 
       
        u3C = fOkada(x  ,p  ,z,q,dip,alpha,typ,'C',3) - ...
              fOkada(x  ,p-W,z,q,dip,alpha,typ,'C',3) - ...
              fOkada(x-L,p  ,z,q,dip,alpha,typ,'C',3) + ...
              fOkada(x-L,p-W,z,q,dip,alpha,typ,'C',3);
        
        
        ux = ux + (U(ityp)/(2*pi))*(u1A - u1Ahat + u1B + z*u1C);

        uy = uy + (U(ityp)/(2*pi))*((u2A - u2Ahat + u2B + z*u2C)*cos(dip) - ...
                                     (u3A - u3Ahat + u3B + z*u3C)*sin(dip));
        
        uz = uz + (U(ityp)/(2*pi))*((u2A - u2Ahat + u2B - z*u2C)*sin(dip) + ...
                                     (u3A - u3Ahat + u3B - z*u3C)*cos(dip));

    end;
end;


