show_sigma.m 5.09 KB
Newer Older
1
function show_sigma(meshfile,bathfile,sigmafile,varargin)
Pierre Cazenave's avatar
Pierre Cazenave committed
2 3
% plot a sigma distribution along a user-selected line
%
Pierre Cazenave's avatar
Pierre Cazenave committed
4
% show_sigma(meshfile,bathfile,sigmafile) 
Pierre Cazenave's avatar
Pierre Cazenave committed
5 6 7 8 9 10 11 12
%
% DESCRIPTION:
%    plot a sigma distribution along a user-selected line
%
% INPUT:
%    meshfile:   fvcom grid file
%    bathfile:   fvcom bathymetry file
%    sigmafile:  fvcom sigma distribution file
13
%    npts:       number of points in the user-selected line (optional,
14
%                defaults to 25).
Pierre Cazenave's avatar
Pierre Cazenave committed
15 16 17 18
%
% OUTPUT:
%
% EXAMPLE USAGE
19
%    show_sigma('tst_grd.dat', 'tst_dep.dat', 'sigma.dat', 50)
Pierre Cazenave's avatar
Pierre Cazenave committed
20
%
21
% Author(s):
Pierre Cazenave's avatar
Pierre Cazenave committed
22
%    Geoff Cowles (University of Massachusetts Dartmouth)
Pierre Cazenave's avatar
Pierre Cazenave committed
23
%    Pierre Cazenave (Plymouth Marine Laboratory)
Pierre Cazenave's avatar
Pierre Cazenave committed
24 25
%
% Revision history
26 27 28
%   2013-01-08 Added more resilient reading of the sigma coordinates file
%   (can now handle comments and blank lines). Also add extra (optional)
%   argument to increase the profile sampling frequency.
Pierre Cazenave's avatar
Pierre Cazenave committed
29 30
%   
%==============================================================================
31

32 33
close all

34 35
global ftbverbose

36
if nargin == 4
37
    npts = varargin{1};
38
else
39
    npts = 25;
40
end
Pierre Cazenave's avatar
Pierre Cazenave committed
41

42
% read the mesh
Pierre Cazenave's avatar
Pierre Cazenave committed
43 44
fid = fopen(meshfile,'r');
if(fid  < 0)
45 46
    error(['file: ' meshfile ' does not exist']);
end
Pierre Cazenave's avatar
Pierre Cazenave committed
47 48 49 50 51 52
C = textscan(fid, '%s %s %s %d', 1);
Nverts = C{4};
C = textscan(fid, '%s %s %s %d', 1);
Nelems = C{4};
x = zeros(Nverts,3);
tri = zeros(Nelems,3);
53 54 55 56 57
if ftbverbose
    fprintf('Nverts: %d\n',Nverts)
    fprintf('Nelems: %d\n',Nelems)
end

58 59 60 61 62 63 64 65 66 67 68 69 70
for i=1:Nelems
    C = textscan(fid, '%d %d %d %d %d', 1);
    tri(i,1) = C{2};
    tri(i,2) = C{3};
    tri(i,3) = C{4};
end
for i=1:Nverts
    C = textscan(fid, '%d %f %f %f ', 1);
    x(i,1) = C{2};
    x(i,2) = C{3};
end

% read the bathy
Pierre Cazenave's avatar
Pierre Cazenave committed
71 72
fid = fopen(bathfile,'r');
if(fid  < 0)
73 74
    error(['file: ' bathfile ' does not exist']);
end
Pierre Cazenave's avatar
Pierre Cazenave committed
75 76
C = textscan(fid, '%s %s %s %d', 1);
for i=1:Nverts
77 78 79
    C = textscan(fid, '%f %f %f ', 1);
    x(i,3) = C{3};
end
Pierre Cazenave's avatar
Pierre Cazenave committed
80

81 82 83 84
if ftbverbose
    fprintf('min topography %f\n',min(x(:,3)));
    fprintf('max topography %f\n',max(x(:,3)));
end
Pierre Cazenave's avatar
Pierre Cazenave committed
85 86 87 88

% read the sigma file
fid = fopen(sigmafile,'r');
if(fid  < 0)
89
    error(['file: ' sigmafile ' does not exist']);
90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
end

while ~feof(fid)
    line = fgetl(fid);
    if isempty(line) || strncmp(line, '!', 1) || ~ischar(line)
        continue
    end
    key = lower(line(1:3));
    C = strtrim(regexpi(line, '=', 'split'));
    switch key
        case 'num'
            nlev = str2double(C{2});
        case 'sig'
            sigtype = C{2};
        case 'du '
            du = str2double(C{2});
        case 'dl '
            dl = str2double(C{2});
        case 'min'
            min_constant_depth = str2double(C{2});
        case 'ku '
            ku = str2double(C{2});
        case 'kl '
            kl = str2double(C{2});
        case 'zku'
            s = str2double(regexp(C{2}, ' ', 'split'));
            zku = zeros(ku, 1);
            for i = 1:ku
                zku(i) = s(i);
            end
        case 'zkl'
            s = str2double(regexp(C{2}, ' ', 'split'));
            zkl = zeros(kl, 1);
            for i = 1:kl
                zkl(i) = s(i);
            end
    end
end
Pierre Cazenave's avatar
Pierre Cazenave committed
128

129 130 131 132 133 134 135 136 137 138 139
if ftbverbose
    fprintf('nlev %d\n',nlev)
    fprintf('sigtype %s\n',sigtype)
    fprintf('du %d\n',du)
    fprintf('dl %d\n',dl)
    fprintf('min_constant_depth %f\n',min_constant_depth)
    fprintf('ku %d\n',ku)
    fprintf('kl %d\n',kl)
    fprintf('zku %d\n',zku)
    fprintf('zkl %d\n',zkl)
end
Pierre Cazenave's avatar
Pierre Cazenave committed
140 141

% generate the sigma coordinates
142 143 144 145 146 147

% fix "java.lang.IllegalArgumentException: adding a container to a container
% on a different GraphicsDevice" error. See
% http://www.mathworks.co.uk/matlabcentral/newsreader/view_thread/169024.
set(0,'DefaultFigureRenderer','OpenGL')

148
figure(1)
Pierre Cazenave's avatar
Pierre Cazenave committed
149
patch('Vertices',[x(:,1),x(:,2)],'Faces',tri,...
150
    'Cdata',x(:,3),'edgecolor','interp','facecolor','interp');
Pierre Cazenave's avatar
Pierre Cazenave committed
151 152 153
axis equal

% plot to get a line
154
fprintf('select two end points of a transect with your mouse... ');
Pierre Cazenave's avatar
Pierre Cazenave committed
155 156
[xt,yt] = ginput(2);
hold on
157
fprintf('done.\n')
Pierre Cazenave's avatar
Pierre Cazenave committed
158 159 160 161 162

ds = (xt(2)-xt(1))/(npts-1);
xline = xt(1):ds:xt(2);
ds = (yt(2)-yt(1))/(npts-1);
yline = yt(1):ds:yt(2);
Pierre Cazenave's avatar
Pierre Cazenave committed
163
plot(xline, yline, 'w+')
164
sline = zeros(1, npts);
Pierre Cazenave's avatar
Pierre Cazenave committed
165
for i=2:npts
166
    sline(i) = sline(i-1) + sqrt((xline(i)-xline(i-1))^2 + (yline(i)-yline(i-1))^2);
167
end
Pierre Cazenave's avatar
Pierre Cazenave committed
168 169

% interpolate the bathymetry along the line
170
zline = griddata(x(:,1), x(:,2), x(:,3), xline, yline);
Pierre Cazenave's avatar
Pierre Cazenave committed
171

172 173
figure(2)
plot(sline, -zline)
Pierre Cazenave's avatar
Pierre Cazenave committed
174 175 176 177

% generate the sigma coordinates along the line
xslice = zeros(npts,nlev);
yslice = zeros(npts,nlev);
178
z = nan(npts, nlev);
Pierre Cazenave's avatar
Pierre Cazenave committed
179 180

% calculate the sigma distributions along the transect
181 182 183 184 185 186 187 188 189 190
switch lower(sigtype)
    case 'generalized'
        for i=1:npts
            z(i,1:nlev) = sigma_gen(nlev,dl,du,kl,ku,zkl,zku,zline(i),min_constant_depth);
        end
    case 'uniform'
        for i=1:npts
            z(i,1:nlev) = 0:-1/double(nlev-1):-1;
        end
    otherwise
191
        error('Can''t do that sigtype')
192
end
Pierre Cazenave's avatar
Pierre Cazenave committed
193 194

for i=1:npts
195
    xslice(i,1:nlev) = sline(i);
196
    yslice(i,1:nlev) = z(i,1:nlev) * zline(i);
197
end
Pierre Cazenave's avatar
Pierre Cazenave committed
198 199 200 201


% plot the mesh along the transect
for k=1:nlev
202 203
    plot(xslice(:,k),yslice(:,k),'k'); hold on;
end
Pierre Cazenave's avatar
Pierre Cazenave committed
204
for k=1:npts
205 206
    plot(xslice(k,:),yslice(k,:),'k'); hold on;
end
Pierre Cazenave's avatar
Pierre Cazenave committed
207 208
axis([xslice(1,1),xslice(end,1),min(yslice(:,end)),5])
title('sigma distribution along the transect');