show_sigma.m 4.97 KB
Newer Older
1
function show_sigma(meshfile,bathfile,sigmafile,varargin)
Pierre Cazenave's avatar
Pierre Cazenave committed
2 3 4 5 6 7 8 9 10 11 12
% plot a sigma distribution along a user-selected line
%
% function [] = show_sigma(meshfile,bathfile,sigmafile) 
%
% 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)
23
%    Pierre Cazenave (Plymouth Marine Laboratory, pica@pml.ac.uk)
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 34
close all

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

40
% read the mesh
Pierre Cazenave's avatar
Pierre Cazenave committed
41 42
fid = fopen(meshfile,'r');
if(fid  < 0)
43 44
    error(['file: ' meshfile ' does not exist']);
end
Pierre Cazenave's avatar
Pierre Cazenave committed
45 46 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);
fprintf('Nverts: %d\n',Nverts)
fprintf('Nelems: %d\n',Nelems)
53 54 55 56 57 58 59 60 61 62 63 64 65
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
66 67
fid = fopen(bathfile,'r');
if(fid  < 0)
68 69
    error(['file: ' bathfile ' does not exist']);
end
Pierre Cazenave's avatar
Pierre Cazenave committed
70 71
C = textscan(fid, '%s %s %s %d', 1);
for i=1:Nverts
72 73 74
    C = textscan(fid, '%f %f %f ', 1);
    x(i,3) = C{3};
end
Pierre Cazenave's avatar
Pierre Cazenave committed
75 76 77 78 79 80 81

fprintf('min topography %f\n',min(x(:,3)));
fprintf('max topography %f\n',max(x(:,3)));

% read the sigma file
fid = fopen(sigmafile,'r');
if(fid  < 0)
82
    error(['file: ' sigmafile ' does not exist']);
83 84 85 86 87 88 89 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
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
121 122 123 124 125 126 127 128 129 130 131 132 133

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)


% generate the sigma coordinates
134 135 136 137 138 139

% 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')

140
figure(1)
Pierre Cazenave's avatar
Pierre Cazenave committed
141
patch('Vertices',[x(:,1),x(:,2)],'Faces',tri,...
142
    'Cdata',x(:,3),'edgecolor','interp','facecolor','interp');
Pierre Cazenave's avatar
Pierre Cazenave committed
143 144 145
axis equal

% plot to get a line
146
fprintf('select two end points of a transect with your mouse');
Pierre Cazenave's avatar
Pierre Cazenave committed
147 148 149 150 151 152 153
[xt,yt] = ginput(2);
hold on

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);
154
plot(xline, yline, 'r+')
155
sline = zeros(1, npts);
Pierre Cazenave's avatar
Pierre Cazenave committed
156
for i=2:npts
157
    sline(i) = sline(i-1) + sqrt((xline(i)-xline(i-1))^2 + (yline(i)-yline(i-1))^2);
158
end
Pierre Cazenave's avatar
Pierre Cazenave committed
159 160

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

163 164
figure(2)
plot(sline, -zline)
Pierre Cazenave's avatar
Pierre Cazenave committed
165 166 167 168

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

% calculate the sigma distributions along the transect
172 173 174 175 176 177 178 179 180 181 182 183
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
        error('cant do that sigtype')
end
Pierre Cazenave's avatar
Pierre Cazenave committed
184 185

for i=1:npts
186 187 188
    xslice(i,1:nlev) = sline(i);
    yslice(i,1:nlev) = z(i,1:nlev)*zline(i);
end
Pierre Cazenave's avatar
Pierre Cazenave committed
189 190 191 192


% plot the mesh along the transect
for k=1:nlev
193 194
    plot(xslice(:,k),yslice(:,k),'k'); hold on;
end
Pierre Cazenave's avatar
Pierre Cazenave committed
195
for k=1:npts
196 197
    plot(xslice(k,:),yslice(k,:),'k'); hold on;
end
Pierre Cazenave's avatar
Pierre Cazenave committed
198 199
axis([xslice(1,1),xslice(end,1),min(yslice(:,end)),5])
title('sigma distribution along the transect');