%this function will determine the best fitting tilt direction and amount by %minimizing energy on the vertical component. %input are X(east), Y(north) and Z(vertical). %output are Znew (tilt corrected Z), Xnew (rotated X), Ynew (rotated Y), %angy (azimuthal rotation in degrees) and vrt (tilt angle from vertical in %degrees) %Note this operates on the second half of the time series-can be modified %as need to operate on whole time series, etc by modifying inline function %ff below. %requires optimization toolbox function [Znew,angy,vrt,Xnew,Ynew] = tilt_corr_rot(X,Y,Z); %% first do polarization analysis on the data to get a guess of the tilt direction. S=cov([Z(end-length(Z)/2:end) X(end-length(Z)/2:end) Y(end-length(Z)/2:end)]); [V,D]=eig(S); linny=1-(D(2,2)+D(1,1))/2/D(3,3); planny=1-2*D(1,1)/(D(2,2)+D(3,3)); if linny>0.1 %this is a relatively linear phase recorded %calculate the rotation angle angy=atan(V(3,3)*sign(V(1,3))/(V(2,3)))*180/pi; end %% search for the best fitting azimuth (angy) and tilt (vrt) [~,idd]=max(abs(V(1,:))); vrt=acosd(V(1,idd)); ff=@(x)Zrot(Z(end-length(Z)/2:end),X(end-length(Z)/2:end),Y(end-length(Z)/2:end),x); options=optimset('Display','iter'); warning off [vrt,zvr2]=fminsearch(ff,[vrt;angy]); angy=vrt(2); vrt=vrt(1); [Znew,Ynew,Xnew]=xyzrot(X,Y,Z,[vrt;angy]); end function zvr = Zrot(Z,X,Y,vrt) [Znew,ynew,xnew]=xyzrot(X,Y,Z,vrt); zvr=sqrt(Znew(:)'*Znew(:)); end function [znew,ynew,xnew]=xyzrot(x,y,z,tst) vrt=tst(1);angy=tst(2); dipZ=-90+vrt;%sin change on vrt due to how it is determined dipY=vrt*cosd(angy); dipX=vrt*sind(angy); %make vectors from dip and azimuth zv=[-sind(dipZ);cosd(dipZ).*cosd(angy);sind(angy).*cosd(dipZ)]; yv=[-sind(dipY);cosd(dipY).*cosd(0);sind(0).*cosd(dipY)]; xv=[-sind(dipX);cosd(dipX).*cosd(90);sind(90).*cosd(dipX)]; C=[zv yv xv]; jk=inv(C)*[z y x]'; znew=jk(1,:)'; ynew=jk(2,:)'; xnew=jk(3,:)'; end