%==========================================================================
% Script for illustrating structure contours of folds
% using simplified Fourier series
% written and copyright by:
%
% Martin Schöpfer
% Fault Analysis Group
% School of Geological Sciences
% University College Dublin
% Belfield, Dublin 4, Ireland
% email: martin@fag.ucd.ie
% Tel.: +353 - 1 - 716 2611
% Fax: +353 - 1 - 716 2607
% Web: http://www.fault-analysis-group.ucd.ie/
%
%==========================================================================
% start of script structure_contours_folds.m
clear

%==========================================================================
% PART 1: Model definition
%==========================================================================

% map dimensions
L_x = 1000;
L_y = 1000;

% height of model
L_z = 1500;

% contour interval, c, and contour vector, v
c = 100;
v = 0:c:L_z;

% grid spacing, s
s = 10;

% define matrices: X- and Y-coordinates, Z (topgraphy), G (geological structure) and O (outrcop)
x = 0:s:L_x;
y = 0:s:L_y;
[X,Y] = ndgrid(x,y);
Z = zeros(size(X));
G = zeros(size(X));
O = zeros(size(X));

% mean elevation of topography and geologic stucture
E = 850;

% amplitude, phase shift, wavelength of topgraphy in x-direction
A_x = 150;
psi_x = 0;
W_x = 1000;

% amplitude, phase shift, wavelength of topgraphy in y-direction
A_y = 150;
psi_y = 1000;
W_y = 2000;

% amplitude and wavelength of folds
A_f = 100;
W_f = 500;

%==========================================================================
% PART 2: User's input
%==========================================================================

% ask user about topography
disp(' ');
disp('Please choose topopgraphy');
q=input('valley (enter 1) or ridge (enter 2) or both (enter 3)? ');
if q == 1
    psi_x = 0;
elseif q == 2
    psi_x = 500;
elseif q == 3
    W_x = 500;
    psi_x = 0;
end

% ask user about geological structure
disp(' ');
disp('Please choose geological structure');
q=input('Box- and Chevron folds (enter 1) or parastitc folds (enter 2)? ');
disp(' ');
if q == 1
    disp('Please enter harmonic coefficients, b1 and b3');
    b_1 = input('b1: ');
    b_3 = input('b3: ');
elseif q == 2
    disp('Please enter harmonic coefficients, b1 and b5');
    b_1 = input('b1: ');
    b_5 = input('b5: ');
end
disp(' ');
phi = input('plunge of fold: ');

%==========================================================================
% PART 3: Calculations of topgraphy (Z), geology (G) and outcrop pattern (O)
%==========================================================================

% calculate topography
Z = A_x*cos(2*pi.*(X-psi_x)./W_x) + A_y*cos(2*pi.*(Y-psi_y)./W_y) + E;

% option 1, Box- and Chevron folds
if q == 1
    G = A_f*b_1*sin(2*pi.*Y./W_f) + A_f*b_3*sin(6*pi.*Y./W_f)...
        + tan(phi*pi/180).*(X - L_x*0.5) + E;
    % option 2, parastitc folds
elseif q == 2
    G = A_f*b_1*sin(2*pi.*Y./W_f) + A_f*b_5*sin(10*pi.*Y./W_f)...
        + tan(phi*pi/180).*(X - L_x*0.5) + E;
end

% outcrop pattern
O = Z - G;

%==========================================================================
% PART 4: Plotting
%==========================================================================

% plot topographic map at sea level and as a 3D object
pcolor(X,Y,Z);shading interp;colormap('summer'); hold on
[T2,t2] = contour(X,Y,Z,v); hold on
surf(X,Y,Z,'EdgeColor','none'); hold on
[T3,t3] = contour3(X,Y,Z,v); hold on

% plot geological structure and structure contours
surface(X,Y,G,'EdgeColor','none','FaceColor','b','FaceAlpha',0.33); hold on
[G3,g3] = contour3(X,Y,G,v); hold on
[G2,g2] = contour(X,Y,G,v); hold on

% plot outcrop pattern at sea level
[O2,o2] = contour(X,Y,O,[0 0]); hold off

% change line styles using object handles t2, t3, g3, g2 and o2
set(t2,'LineWidth',1.5,'Edgecolor','k');
set(t3,'LineWidth',1.5,'Edgecolor','k');
set(g3,'LineWidth',1.5,'Edgecolor','w');
set(g2,'LineWidth',1.5,'Edgecolor','w');
set(o2,'LineWidth',2.5,'Edgecolor',[1 0 .5]);

% set background color
set(gcf,'Color',[1 1 1]);

% set axes properties
set(gca,'Xlim',[0 L_x],'Ylim',[0 L_y],'Zlim',[0 L_z],'LineWidth',1.0,'TickDir','out','FontSize',15);
daspect([1 1 1]);
caxis([(E-A_x) (E+A_x)]);

% 3D persepective, azimuth and elevation, in degrees
view(50,20);

% Use the following line for exporting the figure as a TIFF
% at a resolution given by -r
%print('-f1', 'structure_contour_model.tif','-dtiff','-r500');

% end of script structure_contours_folds.m
%==========================================================================