% xTCSim.m                                          by:  Edward T Peltzer, MBARI
%                                                   revised:  2008 Feb 25.
%
% M-file to simulate a tidal current record and to plot a timeseries'
%   a power-spectrum and a histogram plot for comparison w S4AD data.
%
%   Days = set this parameter to change the number of days plotted
%   Amp  = set this parameter to increase/decrease tidal current
%   Luna = set this parameter for the length of the lunar month


% Clear environment space and set orientation

clear all
close all


% Set parameters for calculation

Days = 120;
Amp = 11;
Luna = 29;


% Generate data and assign data vectors

X = [0:0.25:24*Days];
TC = Amp .* ((sin(X/2) + 0.5 .* sin(X/4)) .* (1 + 0.25 .* sin(X/(2*Luna))));


% FFT of current velocity data for power spectrum

V = fft(TC);
Nv = length(V)+1;
V(1) = [];
Vpwr = abs(V(1:Nv/2)).^2;
Vfreq = (1:Nv/2)/(Nv/2).*2;
Vper = 1./Vfreq;


% Bin data for histogram plot

maxI = length(X);

for I = 1:30
    K(I) = length(find(floor(TC)==I-1))*100/maxI;
    L(I) = length(find(abs(TC) < I))*10/maxI;
end


% Define plotting region and make figure 1

figure(1)
orient tall

left = 1.25/8.5;
bottom = 2.5/11;
width = 6/8.5;
height = 6/11;

subplot('position',[left bottom width height])


% Make plot

plot(X/24,TC,'r-')

hold on


% Define axis and set tick-label-size

axis('xy');
axis([0 120 -20 20]);

h=gca;
set(h,'FontSize',16);
set(h,'FontWeight','normal');
set(h,'LineWidth',1);


% Manually define and control X-ticks and X-tick labels

xtl = [0:5:Days];

N = length(xtl);

for I = 1:N
    if rem(I-1,4) == 0
        if xtl(I) == 0
            XTL(I,:) = [' 0 '];
        elseif xtl(I) < 100
            XTL(I,:) = [' ' num2str(xtl(I))];
        else
            XTL(I,:) = [num2str(xtl(I))];
        end
    else
        XTL(I,:) = ['   '];
    end
end

set(h,'XTick',xtl);
set(h,'XTickLabel',XTL);


% Manually define and control Y-ticks and Y-tick labels

ytl = [-20:5:20];

N = length(ytl);

for I = 1:N
    if rem(I-1,2) == 0
        if ytl(I) == 0
            YTL(I,:) = ['   0 '];
        elseif ytl(I) < 0
        	YTL(I,:) = [' ' num2str(ytl(I)) ' '];
        else
            YTL(I,:) = ['  ' num2str(ytl(I)) ' '];
        end
    else
        YTL(I,:) = ['     '];
    end    
end

set(h,'YTick',ytl);
set(h,'YTickLabel',YTL);


% Label plot -- X-Axis

xlabel(['Time (days)']);
h=get(gca,'xlabel');
set(h,'FontSize',20);
set(h,'FontWeight','normal');


% Label plot -- Y-Axis

ylabel('Current velocity (cm/s)');
h=get(gca,'ylabel');
set(h,'FontSize',20);
set(h,'FontWeight','normal');


% Label plot -- Title

title(['Tidal Current Simulation']);
h=get(gca,'title');
set(h,'FontSize',20);
set(h,'FontWeight','normal');


% Label plot -- Legends

leg1 = sprintf(['Mean abs(Vel) = %4.1f ' setstr(177) ' %3.1f cm/s'],mean(abs(TC)),std(abs(TC)));
h = text(25,15,leg1);
set(h,'Color',[0,0,0]);
set(h,'FontSize',16);
set(h,'FontWeight','normal');

hold off


% Define plotting region and make figure 2

figure(2)
orient tall

subplot('position',[left bottom width height])


% Make plot

plot(Vper,Vpwr/(2*10^8),'LineWidth',1,'Color',[1 0 0],'LineStyle','-')

hold on


% Define axis and set tick-label-size

axis('xy');
axis([0 40 -1 21]);

h=gca;
set(h,'FontSize',16);
set(h,'FontWeight','normal');
set(h,'LineWidth',1);


% Manually define and control X-ticks and X-tick labels

xtl2 = [0:5:40];

N = length(xtl2);

for I = 1:N
    if rem(I-1,2) == 0
        if xtl2(I) == 0
            XTL2(I,:) = [' 0'];
        elseif xtl2(I) < 10
            XTL2(I,:) = [' ' num2str(xtl2(I))];
        else
            XTL2(I,:) = [num2str(xtl2(I))];
        end
    else
        XTL2(I,:) = ['  '];
    end
end

set(h,'XTick',xtl2);
set(h,'XTickLabel',XTL2);


% Manually define and control Y-ticks and Y-tick labels

ytl2 = [0:1:20];

N = length(ytl2);

for I = 1:N
    if rem(I-1,5) == 0
        if ytl2(I) == 0
            YTL2(I,:) = ['  0 '];
        elseif ytl2(I) < 10
        	YTL2(I,:) = ['  ' num2str(ytl2(I)) ' '];
        else
            YTL2(I,:) = [' ' num2str(ytl2(I)) ' '];
        end
    else
        YTL2(I,:) = ['    '];
    end    
end

set(h,'YTick',ytl2);
set(h,'YTickLabel',YTL2);


% Label plot -- X-Axis

xlabel(['Period (hours/cycle)']);
h=get(gca,'xlabel');
set(h,'FontSize',20);
set(h,'FontWeight','normal');


% Label plot -- Y-Axis

ylabel('Power (arbitrary)');
h=get(gca,'ylabel');
set(h,'FontSize',20);
set(h,'FontWeight','normal');


% Label plot -- Title

title(['TC-Sim: Power Spectrum']);
h=get(gca,'title');
set(h,'FontSize',20);
set(h,'FontWeight','normal');


% Label plot -- Legends

leg1 = ['-- Velocity'];
h = text(24,13,leg1);
set(h,'Color',[1,0,0]);
set(h,'FontSize',16);
set(h,'FontWeight','normal');

hold off


% Define plotting region and make figure 3

figure(3)
orient tall

left = 1.25/8.5;
bottom = 2.5/11;
width = 6/8.5;
height = 6/11;

subplot('position',[left bottom width height])


% Make plot

bar([0.5:1:29.5],K,'b')

hold on

plot([0:30],[0 L],'LineWidth',1,'Color',[1 0 0],'LineStyle','-')


% Define axis and set tick-label-size

axis('xy');
axis([0 30 -1 11]);

h=gca;
set(h,'FontSize',16);
set(h,'FontWeight','normal');
set(h,'LineWidth',1);


% Manually define and control X-ticks and X-tick labels

xtl3 = [0:5:30];

N = length(xtl3);

for I = 1:N
    if rem(I-1,2) == 0
        if xtl3(I) == 0
            XTL3(I,:) = [' 0'];
        elseif xtl3(I) < 10
            XTL3(I,:) = [' ' num2str(xtl3(I))];
        else
            XTL3(I,:) = [num2str(xtl3(I))];
        end
    else
        XTL3(I,:) = ['  '];
    end
end

set(h,'XTick',xtl3);
set(h,'XTickLabel',XTL3);


% Manually define and control Y-ticks and Y-tick labels

ytl3 = [0:1:10];

N = length(ytl3);

for I = 1:N
    if rem(I-1,2) == 0
        if ytl3(I) == 0
            YTL3(I,:) = ['  0 '];
        elseif ytl3(I) < 10
        	YTL3(I,:) = ['  ' num2str(ytl3(I)) ' '];
        else
            YTL3(I,:) = [' ' num2str(ytl3(I)) ' '];
        end
    else
        YTL3(I,:) = ['    '];
    end    
end

set(h,'YTick',ytl3);
set(h,'YTickLabel',YTL3);


% Label plot -- Title

title(['TC-Sim: Histogram']);
h=get(gca,'title');
set(h,'FontSize',20);
set(h,'FontWeight','normal');


% Label plot -- X-Axis

xlabel(['Velocity (cm/s)']);
h=get(gca,'xlabel');
set(h,'Color',[0,0,0]);
set(h,'FontSize',20);
set(h,'FontWeight','normal');


% Label plot -- Y-Axis

ylabel('Frequency (%)');
h=get(gca,'ylabel');
set(h,'Color',[0,0,1]);
set(h,'FontSize',20);
set(h,'FontWeight','normal');


% Label plot -- Y-Axis 2

ylab2 = ['Cumulative frequency (%)'];
h = text(33.5,8.3,ylab2);
set(h,'Color',[1,0,0]);
set(h,'FontSize',20);
set(h,'FontWeight','normal');
set(h,'VerticalAlignment','middle','HorizontalAlignment','left')
set(h,'Rotation',-90);

leg10 = ['100'];
h = text(30.2,10,leg10);
set(h,'Color',[0,0,0]);
set(h,'FontSize',16);
set(h,'FontWeight','normal');

leg08 = [' 80'];
h = text(30.2,8,leg08);
set(h,'Color',[0,0,0]);
set(h,'FontSize',16);
set(h,'FontWeight','normal');

leg06 = [' 60'];
h = text(30.2,6,leg06);
set(h,'Color',[0,0,0]);
set(h,'FontSize',16);
set(h,'FontWeight','normal');

leg04 = [' 40'];
h = text(30.2,4,leg04);
set(h,'Color',[0,0,0]);
set(h,'FontSize',16);
set(h,'FontWeight','normal');

leg02 = [' 20'];
h = text(30.2,2,leg02);
set(h,'Color',[0,0,0]);
set(h,'FontSize',16);
set(h,'FontWeight','normal');

leg00 = [' 0'];
h = text(30.2,0,leg00);
set(h,'Color',[0,0,0]);
set(h,'FontSize',16);
set(h,'FontWeight','normal');


% Label plot -- Legends

leg1 = sprintf(['Mean = %4.1f ' setstr(177) ' %3.1f cm/s'],mean(abs(TC)),std(abs(TC)));
h = text(15,8.5,leg1);
set(h,'Color',[0,0,0]);
set(h,'FontSize',16);
set(h,'FontWeight','normal');

leg2 = sprintf(['Median = %4.1f cm/s'],median(abs(TC)));
h = text(15,7.5,leg2);
set(h,'Color',[0,0,0]);
set(h,'FontSize',16);
set(h,'FontWeight','normal');

leg3 = sprintf(['Max velocity = %4.1f cm/s'],max(abs(TC)));
h = text(15,6.5,leg3);
set(h,'Color',[0,0,0]);
set(h,'FontSize',16);
set(h,'FontWeight','normal');

leg4 = sprintf(['Min velocity = %4.2f cm/s'],min(abs(TC)));
h = text(15,5.5,leg4);
set(h,'Color',[0,0,0]);
set(h,'FontSize',16);
set(h,'FontWeight','normal');

leg5 = sprintf(['Freq < 10 cm/s = %4.1f %%'],length(find(abs(TC)<10))*100/maxI);
h = text(15,3.5,leg5);
set(h,'Color',[0,0,0]);
set(h,'FontSize',16);
set(h,'FontWeight','normal');

leg6 = sprintf(['Freq  <  1 cm/s  =  %3.1f %%'],length(find(abs(TC)<1))*100/maxI);
h = text(15,2.5,leg6);
set(h,'Color',[0,0,0]);
set(h,'FontSize',16);
set(h,'FontWeight','normal');

hold off
