Source Code:

%-------------------------------------------------------------------------------
% Name:
%   .dbase/web_programs/plotpro/matlab_scripts/mountain.m
% Description:
%   Creates the monthly spectral energy plot (mountains).
%   A few lines are pre-pended to this script.
% Called by:
%   .f90/plot_utils.make_plot
%-------------------------------------------------------------------------------

path(path,'/project/dbase/web_programs/plotpro/matlab_scripts');

wkg = wkg_path;


% LOAD DATA, GIVE SHORT NAMES FOR LOADED VARIABLES
% NOTE: x = frequency data, y = energy density data

eval(['freq=load(''',wkg,'freq.',pid,'01'');']);
eval(['ener=load(''',wkg,'ener.',pid,'01'');']);

x = freq(1,:);
y = ener;

dt = read_date_file([pid,'01'],wkg);
tim = (dt.day-1)*24 + (dt.hour + dt.min/60); % in hours

yearchar = num2str(dt.year(1));
monthnum = dt.month(1);

clear freq ener date


% CHECK TO SEE IF THERE IS ANY DATA

if ~ ( any(x>0) & any(y>0) )
  clear all;
  quit;
end


% SET FIGURE WINDOW

figure('units','normal','paperunits','normal','position',...
        [.1,.1,.5,.8],'paperposition',[.05,.05,.9,.93],'visible','off');

sz = size(y);
map = [ 0 0 0; 1 0 0; 0 1 0; 0 0 1; .5 .5 .5];
colormap(map);
cd = colormap;


% TAG FIRST RUN OF EACH DAY

pday = 0;
for i=1:sz(1)
  day = floor(tim(i)/24) + 1;
  if day ~= pday
    nday(i) = 1;
    pday = day;
  else
    nday(i) = 0;
  end
end


% SCALE Y VALUES

yscale = 7;
ym = max(max(y));
% Fixed scale: set ym=50.;
%ym=50.;
y = yscale*24*y./ym;

logx = log(x);


% PLOT DATA

eps = 5;
prvtime = 790;
j = 0;
for i=sz(1):-1:1
  if ( (prvtime - tim(i) > eps) | nday(i) == 1) 
    patch('xdata',logx,'ydata',tim(i)+y(i,:),'edgecolor','none','facecolor',cd(1+rem(j,5),:));
    j = j + 1;
    prvtime = tim(i);
  end
end

yaxisheight = .76;
set(gca,'position',[0.1 0.1 0.7 yaxisheight]);


% GET STATION INFO

info = read_info_file([pid,'01'],wkg);


% IDENTIFY IF SURGE BASIN OR STANDARD PROCESSING

surge = 0;
basin = 0;
standard = 0;

if ~isempty( findstr(info.stream_name,'SURGE') )
  surge = 1;
elseif ~isempty( findstr(info.stream_name,'BASIN') )
  basin = 1;
else
  standard = 1;
end

% WRITE STATION INFO ON PLOT

if ( surge == 1 )
  xtext=log(1/128);
elseif ( basin == 1 )
  xtext=log(1/22);
else
  xtext=log(.10);
end

ytext= 860;
txtspace = 24;

text(xtext,ytext+3*txtspace,[char(info.stnid),' ',char(info.stname)],...
  'horizontalalignment','center','fontsize',14);

stream_name = ['(',char(info.stream_name),')'];
if ( length(stream_name) > 5 )
  text(xtext, ytext+2*txtspace, stream_name, 'horizontalalignment',...
  'center','fontsize',12);
end

mname = get_monthname(monthnum,'long');
text(xtext,ytext + -0.3*txtspace,[char(mname),' ',yearchar],'horizontalalignment',...
'center','fontsize',14);

text(xtext,ytext+0.7*txtspace,'ENERGY SPECTRUM','horizontalalignment','center',...
     'fontsize',14);

% MAKE Y AXIS LABELS

ylim = 824;
set(gca,'ylim',[0,ylim],'ytick',[0:24:720],'yticklabel',num2str([1:31]'));
ylabel('Day of Month (UTC)','fontsize',12);


% MAKE X AXIS LABELS

numticks = 5;

if ( surge == 1 )
  lxlim = 1/1024; % (Hz)
  uxlim = 1/16;   % (Hz)
  f_sig_digs = 2;
  p_sig_digs = 4;
  f_labels = 0.20;
elseif ( basin == 1 )
  lxlim = 1/128; % (Hz)
  uxlim = 1/4;   % (Hz)
  f_sig_digs = 2;
  p_sig_digs = 3;
  f_labels = 0.20;
else
  lxlim = 1/25; % (Hz) 
  uxlim = 1/4; % (Hz)
  f_sig_digs = 2;
  p_sig_digs = 2;
  f_labels = 0.07;
end

% note: xtl = "x tick label"

diff = (log(uxlim) - log(lxlim))/(numticks-1);
xt = log(lxlim) + (0:(numticks-1))*diff;
xt(numticks) = log(uxlim);
fxtl = num2str(exp(xt'),f_sig_digs);
pxtl = num2str(1./exp(xt'),p_sig_digs);
set(gca,'xlim',[log(lxlim),log(uxlim)],'xtick',xt,'xticklabel',fxtl);
hold on;
plot([log(lxlim) log(lxlim)],[0 840],'k-');
box on; 
text(xt,-35*ones(1,numticks),pxtl,'horizontalalignment','center');
text(log(uxlim)+f_labels,-15,'Frequency (Hz)','horizontalalignment','left',...
     'fontsize',10);
text(log(uxlim)+f_labels,-35,'Period (s)','horizontalalignment','left',...
     'fontsize',10);


% PUT SCALE ON PLOT

if ym < .1
  rndval = .01;
elseif ym < 1
  rndval = .1;
elseif ym < 10
  rndval = 1;
else
  rndval = 10;
end

rnd = rndval*round(ym/rndval)
yscaleheight = (yaxisheight/ylim) * (rnd/ym)*yscale*24;
sh=subplot('position',[.87,.60,.01,yscaleheight]);

hold on;

%draw the scale line (3 units long)

ytick = [0,rnd/4,rnd/2,3*rnd/4,rnd];
yticklabs = num2str(ytick',4);
plot([0,0],[0,rnd]);
hold on;
set(gca,'ylim',[0,rnd],'ytick',ytick,'yticklabel',yticklabs,'xticklabel',[]);
text(+9,rnd/2,'Energy Density (m^2/Hz)','rotation',90,'horizontalalignment','center','fontsize',12);


% PUT BLUE BOX AROUND PLOT AND ADD WEB ADDRESS

blue_box(0.01);

[pr_axes,ax_orig]=page_relative_axes;

axes(pr_axes);
text(.1,.035,'(Colors delineate data records)','FontSize',11);
axes(ax_orig);

% PRINT THE PLOT

set(gcf, 'visible','off');

png_name = [wkg,'plot_',pid,'.png'];
print('-dpng','-r75',png_name);


quit