show_sigma.m 4.84 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 14
%    spoints:    number of points in the user-selected line (optional,
%                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 35 36 37 38 39 40
% meshfile = 'tst_grd.dat';
% bathfile = 'tst_dep.dat';
% sigmafile = 'sigma.dat';
close all

if nargin == 4
    spoints = varargin{1}
else
    spoints = 25;
end
Pierre Cazenave's avatar
Pierre Cazenave committed
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84

% read the mesh 
fid = fopen(meshfile,'r');
if(fid  < 0)
  error(['file: ' meshfile ' does not exist']);
end;
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)
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 
fid = fopen(bathfile,'r');
if(fid  < 0)
  error(['file: ' bathfile ' does not exist']);
end;
C = textscan(fid, '%s %s %s %d', 1);
for i=1:Nverts
  C = textscan(fid, '%f %f %f ', 1);
  x(i,3) = C{3};
end;

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)
  error(['file: ' sigmafile ' does not exist']);
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 121 122
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
123 124 125 126 127 128 129 130 131 132 133 134 135 136

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
fprintf('select two endpoints of a transect')
137
figure(1)
Pierre Cazenave's avatar
Pierre Cazenave committed
138 139 140 141 142 143 144 145 146 147 148 149 150 151 152
patch('Vertices',[x(:,1),x(:,2)],'Faces',tri,...
       'Cdata',x(:,3),'edgecolor','interp','facecolor','interp');
axis equal

% plot to get a line
%fprintf('select two end points of a transect with your mouse');
[xt,yt] = ginput(2);
hold on

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

% interpolate the bathymetry along the line
159
zline = griddata(x(:,1),x(:,2),x(:,3),xline,yline);
Pierre Cazenave's avatar
Pierre Cazenave committed
160 161 162 163 164 165 166 167 168

figure
plot(sline,-zline)

% generate the sigma coordinates along the line
xslice = zeros(npts,nlev);
yslice = zeros(npts,nlev);

% calculate the sigma distributions along the transect
169 170 171 172 173 174 175 176 177 178 179 180
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
181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198

for i=1:npts
  xslice(i,1:nlev) = sline(i);
  yslice(i,1:nlev) = z(i,1:nlev)*zline(i);
end;


% plot the mesh along the transect
for k=1:nlev
  plot(xslice(:,k),yslice(:,k),'k'); hold on;
end;
for k=1:npts
  plot(xslice(k,:),yslice(k,:),'k'); hold on;
end;
axis([xslice(1,1),xslice(end,1),min(yslice(:,end)),5])
title('sigma distribution along the transect');