Commit 768c5e8b authored by Ricardo Torres's avatar Ricardo Torres 💬
Browse files

Added coloured vector plots, ERSEM OBC write out is more verbouse and input...

Added coloured vector plots, ERSEM OBC write out is more verbouse and input data is interpolated to same frequency for all obc variables. Removed double precision of X,Y in nesting file to enable single precision runs
parent 894fc35a
......@@ -175,6 +175,7 @@ nemo.o = (nemo.o / 16) *1000 ./ tmp; % Nemo oxygen concentrations are for O rath
nemo.p = (nemo.p / 35.5)*1000 ./ tmp;%g/s to mmol/m3
nemo.sio3 = (nemo.sio3 / 28) *1000./ tmp;%g/2 to mmol/m3
nemo.bioalk= nemo.bioalk./ tmp / 1000; % bioalk is in umol/s need umol/kg
nemo.dic= nemo.dic./12./ tmp *1000; % dic is in gC/s need mmol/m3
% total alkalinity is already in umol/Kg as expected by ERSEM.
clear tmp
......
......@@ -193,45 +193,45 @@ Times_varid = netcdf.defVar(nc, 'Times' ,'NC_CHAR', ...
[datestrlen_dimid, time_dimid]);
netcdf.putAtt(nc, Times_varid, 'time_zone', 'UTC');
x_varid = netcdf.defVar(nc, 'x', 'double', ...
x_varid = netcdf.defVar(nc, 'x', 'NC_FLOAT', ...
node_dimid);
netcdf.putAtt(nc, x_varid, 'units', 'meters');
netcdf.putAtt(nc, x_varid, 'long_name', 'nodal x-coordinate');
y_varid = netcdf.defVar(nc, 'y', 'double', ...
y_varid = netcdf.defVar(nc, 'y', 'NC_FLOAT', ...
node_dimid);
netcdf.putAtt(nc, y_varid, 'units', 'meters');
netcdf.putAtt(nc, y_varid, 'long_name', 'nodal y-coordinate');
xc_varid = netcdf.defVar(nc, 'xc', 'double', ...
xc_varid = netcdf.defVar(nc, 'xc', 'NC_FLOAT', ...
elem_dimid);
netcdf.putAtt(nc, xc_varid, 'units', 'meters');
netcdf.putAtt(nc, xc_varid, 'long_name', 'zonal x-coordinate');
yc_varid = netcdf.defVar(nc, 'yc', 'double', ...
yc_varid = netcdf.defVar(nc, 'yc', 'NC_FLOAT', ...
elem_dimid);
netcdf.putAtt(nc, yc_varid, 'units', 'meters');
netcdf.putAtt(nc, yc_varid, 'long_name', 'zonal y-coordinate');
lon_varid = netcdf.defVar(nc, 'lon', 'double', ...
lon_varid = netcdf.defVar(nc, 'lon', 'NC_FLOAT', ...
node_dimid);
netcdf.putAtt(nc, lon_varid, 'units', 'degrees_east');
netcdf.putAtt(nc, lon_varid, 'standard_name', 'longitude');
netcdf.putAtt(nc, lon_varid, 'long_name', 'nodal longitude');
lat_varid = netcdf.defVar(nc, 'lat', 'double', ...
lat_varid = netcdf.defVar(nc, 'lat', 'NC_FLOAT', ...
node_dimid);
netcdf.putAtt(nc, lat_varid, 'units', 'degrees_north');
netcdf.putAtt(nc, lat_varid, 'standard_name', 'latitude');
netcdf.putAtt(nc, lat_varid, 'long_name', 'nodal latitude');
lonc_varid = netcdf.defVar(nc, 'lonc', 'double', ...
lonc_varid = netcdf.defVar(nc, 'lonc', 'NC_FLOAT', ...
elem_dimid);
netcdf.putAtt(nc, lonc_varid, 'units', 'degrees_east');
netcdf.putAtt(nc, lonc_varid, 'standard_name', 'longitude');
netcdf.putAtt(nc, lonc_varid, 'long_name', 'zonal longitude');
latc_varid = netcdf.defVar(nc, 'latc', 'double', ...
latc_varid = netcdf.defVar(nc, 'latc', 'NC_FLOAT', ...
elem_dimid);
netcdf.putAtt(nc, latc_varid, 'units', 'degrees_north');
netcdf.putAtt(nc, latc_varid, 'standard_name', 'latitude');
......
......@@ -130,7 +130,7 @@ end
%--------------------------------------------------------------------------
% Set netCDF variables and dump to file
%--------------------------------------------------------------------------
if FileExist
if FileExists
% open boundary forcing
nc = netcdf.open(tsOBCFile, 'WRITE');
% read dimensions from the
......@@ -268,27 +268,30 @@ if ~FileExists
end
for nuts=1:NNuts
eval(['netcdf.putVar(nc,',varidN{nuts},',Mobj.(ERSEMdata(nuts).name));'])
end
disp(['Finished with variable, ',ERSEMdata(nuts).name])
end
else % file exist and time could be different... check and interpolate if necessary
if length(time) < length(file_timem)
for nuts=1:NNuts
data= Mobj.(ERSEMdata(nuts).name);
dataint = nans(size(data, 1),size(data, 2),length(file_timem));
dataint = nan(size(data, 1),size(data, 2),length(file_timem));
for nn=1:size(data, 1)
[X1,Y1]=meshgrid(file_timem,siglay(nn,:));
[X1,Y1]=meshgrid(file_timem- 678942,siglay(nn,:));
[X,Y]=meshgrid(time,siglay(nn,:));
% interpolate ERSEMdata...
dataint(nn,:,:) = interp2(X,Y,squeeze(data(nn,:,:),X1,Y1));
dataint(nn,:,:) = interp2(X,Y,squeeze(data(nn,:,:)),X1,Y1);
end
netcdf.putVar(nc,varidN{nuts},dataint);
Nut_id= netcdf.inqVarID(nc,ERSEMdata(nuts).name);
netcdf.putVar(nc,Nut_id,dataint);
disp(['Finished with variable, ',ERSEMdata(nuts).name])
end
else
% everything is in the same time frequency
for nuts=1:NNuts
eval(['netcdf.putVar(nc,',varidN{nuts},',Mobj.(ERSEMdata(nuts).name));'])
end
disp(['Finished with variable, ',ERSEMdata(nuts).name])
end
end
end
......
function [Plots]=do_vector_plot_MatlabMapC(plotOPTS,FVCOM,tt)
%
% Function to display vector maps of FVCOM currents (i.e. U,V)
%
% [Plots]=do_vector_plot(plotOPTS,FVCOM)
%
% DESCRIPTION:
% Generates vector maps of currents using m_map toolbox
%
% INPUT:
% plotOPTS = structure array with predefined options for generating the
% maps
% FVCOM = data extracted from FVCOM NC file. See read_netCDF_FVCOM for
% details
%
% plotOPTS.range_lat: [50.1000 50.4000]
% plotOPTS.range_lon: [-4.5000 -3.8500]
% plotOPTS.fig_name: 'co2_S5_slowleak'
% plotOPTS.mesh: [1x1 struct]
% plotOPTS.coastline_file: '../mat/tamar3_0coast.mat'
% plotOPTS.zone: 30
% plotOPTS.ell: 'grs80'
% plotOPTS.var_plot: 'PH'
% plotOPTS.clims: [6 8]
% plotOPTS.do_mesh: 0
% plotOPTS.nz_plot: 1
% plotOPTS.figure: 1
% plotOPTS.Time_record: 7.3271e+05
% plotOPTS.nz_plot_vec: 1
% plotOPTS.data_dec: 5
% plotOPTS.vel_sca: 5
% plotOPTS.pause: 0.5000
%
% OUTPUT:
% Plots = structure array with figure handles
%
% EXAMPLE USAGE
% [Plots]=do_vector_plot(plotOPTS,FVCOM)
%
% Author(s):
% Ricardo Torres and Pierre Cazenave (Plymouth Marine Laboratory)
%
% Revision history
%
%==============================================================================
% adds m_map to matlab paths. file is in utilities directory.
% ammend according to your m_map installation paths
figure(plotOPTS.figure);
% generate figure with correct projection lat and lon range ellipsoid and
% zone.
axesm('mercator','MapLatLimit',plotOPTS.range_lat,'MapLonLimit',[plotOPTS.range_lon],'MeridianLabel','on',...
'ParallelLabel','on','MLineLocation',[0.1],'PLineLocation',[0.1],'LabelUnits','dm')
% add coastline if present
if (isfield(plotOPTS,'coastline_file') && ~isempty(plotOPTS.coastline_file) )
coast=load(plotOPTS.coastline_file);
geoshow([coast.ncst(:,2)],[coast.ncst(:,1)],'Color','black')
end
%-----------------------------------------------------------------------
% Convert element positions from FVCOM to lat and lon using m_ll2ll.m from
% utilities directory. This accesses m_map functions.
% In my case it needs adding 6 because of discrepancies between proj and m_map.
% Proj automatically determines a
% reference long in strides of 6deg while m_map doesn't
%------------------------------------------------------------------------
% only plot vectors inside lat and lon range and ...
if ~isfield(plotOPTS.mesh,'latc')
try
plotOPTS.mesh.latc = nodes2elems(plotOPTS.mesh.lat,plotOPTS.mesh);
plotOPTS.mesh.lonc = nodes2elems(plotOPTS.mesh.lon,plotOPTS.mesh);
catch
disp(['We have no access to nodes2elems in the fvcom_prepro directory'])
end
end
igood = find ( plotOPTS.mesh.lonc < plotOPTS.range_lon(2) & plotOPTS.mesh.lonc > plotOPTS.range_lon(1) &...
plotOPTS.mesh.latc < plotOPTS.range_lat(2) & plotOPTS.mesh.latc > plotOPTS.range_lat(1));
% decimate positions. Plot every plotOPTS.data_dec position.
igood=igood(1:plotOPTS.data_dec:end);
%------------------------------------------------------------------------
% Select how many layers to plot
%------------------------------------------------------------------------
ND=ndims(FVCOM.(plotOPTS.var_plotu))
switch ND
case 2
nLayers=1;
colourSpec=[0 0 0];
case 3
if isfield(plotOPTS,'nz_plot_vec')
nLayers=size(plotOPTS.nz_plot_vec,2);
nLayersRange=plotOPTS.nz_plot_vec;
else
nLayers=size(plotOPTS.nz_plot,2);
nLayersRange=plotOPTS.nz_plot;
end
% choose colors for vectors
if nLayers==1
colourSpec=[0 0 0];
else
colourSpec=colormap(hsv(nLayers));
end
end
% Preallocate outputs
u=nan(size(igood,1),nLayers,1);
v=nan(size(igood,1),nLayers,1);
% Check if we're running
% for aa=1:length(plotOPTS.Time_record)
aa=tt;
fprintf('Time step %i of %i\n',aa,length(FVCOM.Time_record));
% Mesh goes underneath vectors.
% if plotOPTS.do_mesh
% % plot vertices
% Plots(plotOPTS.figure).mesh=patch('Vertices',[double(plotOPTS.mesh.lat),double(plotOPTS.mesh.lon)],'Faces',plotOPTS.mesh.tri,...
% 'EdgeColor',[0.6 0.6 0.6],'FaceColor','none');hold on
% end
switch ND
case 2
u(:,1)=squeeze(FVCOM.(plotOPTS.var_plotu)(igood,(aa)));
v(:,1)=squeeze(FVCOM.(plotOPTS.var_plotv)(igood,(aa)));
case 3
for ii=1:nLayers
u(:,ii)=squeeze(FVCOM.(plotOPTS.var_plotu)(igood,nLayersRange(ii),(aa)));
v(:,ii)=squeeze(FVCOM.(plotOPTS.var_plotv)(igood,nLayersRange(ii),(aa)));
end
end
for jj=1:nLayers
% [Plots(plotOPTS.figure).handles{jj}]=quiver(double(plotOPTS.mesh.lonc(igood)),double(plotOPTS.mesh.latc(igood)),u,v,plotOPTS.vel_sca,'k','filled');
[Plots(plotOPTS.figure).handles{jj}]=quiverwcolorbar(double(plotOPTS.mesh.lonc(igood)),double(plotOPTS.mesh.latc(igood)),u,v,...
plotOPTS.vel_sca,'bounds',plotOPTS.clims)
% [Plots(plotOPTS.figure).handles{jj}]=quiverm(double(FVCOM.latc(igood)),double(FVCOM.lonc(igood)),u./100,v./100,'k',0,'filled'),
% [Plots(plotOPTS.figure).handles(jj),~]=m_vec(plotOPTS.vel_sca,x(igood),y(igood),...
% u(:,jj),v(:,jj),colourSpec(jj,:),...
% 'shaftwidth',1,'headwidth',2);
% hold on
end
pause(plotOPTS.pause)
% if aa~=length(plotOPTS.Time_record)
% delete(Plots(plotOPTS.figure).handles(:))
% end
% end
function hh = quiverwcolorbar(varargin)
% Melissa Day 5/24/2013
% Upgrade of Andrew Diamond's quiverc2wcmap to generate a quiver plot
% with arrows colored according to vector magnitude.
% Functional differences from quiverc2wcmap:
% 1) Allows user to specify colormap boundaries using 'bound': changes
% colorbar axes AND corresponding vector coloring
% (much more useful for intercomparison of datasets)
% 2) Improved fidelity to magnitude colors in large datasets
% (at a cost of increased computation time...)
% Uses some improvements from DS's vfield_color
% 3) Small clarifications added
% Example:
% x = rand(1,50).*100;
% y = rand(1,50).*100;
% u = rand(1,50) .* 10;
% v = rand(1,50) .* 10;
% scale = 0;
% figure; quiverwcolorbar(x',y',u',v',scale); %compare to:
% figure; quiverwcolorbar(x',y',u',v',scale,'bounds',[0 10]);
%----------------
% function hh = quiverc2wcmap(varargin)
% Andrew Diamond 3/17/2005
% This is based off of Tobias Höffken which was based off of Bertrand Dano
% keeping with their main intent to show a color w/r to magnitude quiver plot
% while maintaining a large amount of quiver API backwards compatability.
% Functional differences from quiverc2:
% 1) This works under 6.5.1
% 2) It handles NaNs
% 3) It draws a colormap that is w/r to the quiver magnitudes (hard coded to
% 20 segments of the colormap as per quiverc2 - seems fine for a quiver).
% 4) Various bug fixes (I think)
% In order to do this I needed some small hacks on 6.5.1's quiver to take a
% color triplet. I have included as part of this file in a subfunction below.
%----------------
% Comments from quiverc2
% changed Tobias Höffken 3-14-05
% totally downstripped version of the former
% split input field into n segments and do a quiver qhich each of them
% Modified version of Quiver to plots velocity vectors as arrows
% with components (u,v) at the points (x,y) using the current colormap
% Bertrand Dano 3-3-03
% Copyright 1984-2002 The MathWorks, Inc.
% changed T. Höffken 14.03.05, for high data resolution
% using fixed color "spacing" of 20
%QUIVERC Quiver color plot.
% QUIVERC(X,Y,U,V) plots velocity vectors as arrows with components (u,v)
% at the points (x,y). The matrices X,Y,U,V must all be the same size
% and contain corresponding position and velocity components (X and Y
% can also be vectors to specify a uniform grid). QUIVER automatically
% scales the arrows to fit within the grid.
%
% QUIVERC(U,V) plots velocity vectors at equally spaced points in
% the x-y plane.
%
% QUIVERC(U,V,S) or QUIVER(X,Y,U,V,S) automatically scales the
% arrows to fit within the grid and then stretches them by S. Use
% S=0 to plot the arrows without the automatic scaling.
%
% QUIVERC(...,LINESPEC) uses the plot linestyle specified for
% the velocity vectors. Any marker in LINESPEC is drawn at the base
% instead of an arrow on the tip. Use a marker of '.' to specify
% no marker at all. See PLOT for other possibilities.
%
% QUIVERC(...,'filled') fills any markers specified.
%
% H = QUIVERC(...) returns a vector of line handles.
%
% Example:
% [x,y] = meshgrid(-2:.2:2,-1:.15:1);
% z = x .* exp(-x.^2 - y.^2); [px,py] = gradient(z,.2,.15);
% contour(x,y,z), hold on
% quiverc(x,y,px,py), hold off, axis image
%
% See also FEATHER, QUIVER3, PLOT.
% Clay M. Thompson 3-3-94
% Copyright 1984-2002 The MathWorks, Inc.
% $Revision: 5.21 $ $Date: 2002/06/05 20:05:16 $
%-------------------------------------------------------------
nin = nargin; %number of inputs
% error(nargchk(2,5,nin));
error(nargchk(2,7,nin)); %added +2 to maxargs to account for 'bounds' add
% Check numeric input arguments
if nin<4, % quiver(u,v) or quiver(u,v,s)
[msg,x,y,u,v] = xyzchk(varargin{1:2});
else % quiver(x,y,u,v) and beyond
[msg,x,y,u,v] = xyzchk(varargin{1:4});
end
if ~isempty(msg), error(msg); end
scale=1; % This is the default I think.
if(nin == 3) % quiver(u,v,s)
if(isscalar(varargin{nin}))
scale = varargin{nin};
end
elseif(nin >= 5) % quiver(x,y,u,v,s) or quiver(x,y,u,v,s,'bounds',[start end])
if(isscalar(varargin{5}))
scale = varargin{5};
end
end
%-------------Define matrix of vector magnitudes-------------
% Define matrix of vector magnitudes
vr = sqrt(u.^2+v.^2);
% From quiverc2wcmap:
% if data has Nan, don't let it contaminate the computations that segment the
% data; I could just do this with vr and then get clever with the indices but
% this make for easy debugging and as this is a graphics routine the computation
% time is completely irrelevant.
nonNaNind = find(~isnan(vr(:)));
xyuvvrNN = [x(nonNaNind),y(nonNaNind),u(nonNaNind),v(nonNaNind),vr(nonNaNind)];
[xyuvvrNNs, xyuvvrNNsi] = sortrows(xyuvvrNN,5);
% From quiverc2wcmap, no longer necessary
% n = 20; %number of colors
% CC = colormap;
% iCCs=round(linspace(1,size(CC,1),n));
% iData = round(linspace(0,size(xyuvvrNNs,1),n+1));
% figure;
%-------------Generate colorbar-------------
% Includes a clever way to generate (and subsequently hide) an image that
% is required for colorbar without running out of memory for large datasets.
%
% Condensed comments from quiverc2wcmap:
% In 6.5.1 if you ever want a colorbar tick marks to reflect real data ranges
% (versus just the indices of a colormap) then there apparently has to be an
% image somewhere in the figure. Of course, I don't want an image but I figured
% I just make it invisible and then draw the quiver plot over it.
% Unfortunately, it seems that colorbar uses caxis to retrive the data range in
% the image and for invisible images it always seems to be 0 UNLESS you
% explictly reset the caxis.
% This will work but then the axis will be oversized to accomodate the image
% because images have their first and last vitual "pixels" CENTERED around the
% implicit or explict xmin,xmax,ymin,ymax (as per imagesc documentation) but the
% start and end of each of these "pixels" is +/- half a unit where the unit
% corresponds to subdividing the limits by the number of pixels (-1). Given
% that formula and given my invisible 2x2 image for which it is desired to
% diplay in an axis that ISN'T oversized (i.e. min(x), max(x), min(y),max(y)) it
% is possible to solve for the limits (i.e. an artifically reduced limit)
% that need to be specified for imagesc to make its final oversized axis to be the
% non oversized axis that we really want.
% xa,xb,ya,yb compenstates for the axis extention given by imagesc to make
% pixels centered at the limit extents (etc.). Note, this "easy"
% formula would only work for 2x2 pixel images.
xs = min(x); %x(1);
xf = max(x); %x(end)
xa = (3 * xs + xf)/4;
xb = (3 * xf + xs)/4;
ys = min(y); %y(1);
yf = max(y); %y(end)
ya = (3 * ys + yf)/4;
yb = (3 * yf + ys)/4;
% Determine magnitude min/max (which is reflected in colorbar)
colormin = min(xyuvvrNNs(:,5)); %column 5 is NaN-cleared vr
colormax = max(xyuvvrNNs(:,5));
% Allow user to edit bounds using "'bounds',[colormin colormax]" input
for k=1:nin
if (k~=nin) && (length(varargin{k})==6) && strcmp(varargin{k},'bounds')
bounds = varargin{k+1};
if isempty(bounds), error('Specify colormap boundaries'); end
colormin = bounds(1);
colormax = bounds(2);
end
end
% mapbounds = reshape(xyuvvrNNs([1,end;1,end],5),2,2); %from quiverc2wcmap
mapbounds = [colormin colormax; colormin colormax];
h=imagesc([xa,xb],[ya,yb],mapbounds);
set(h,'visible','off');
% Prep colorbar
rang = (colormax-colormin)/colormax;
ticknum = 6; %if you want to toggle number of ticks on colorbar
incr = rang./(ticknum-1);
B = [colormin/colormax:incr:1];
B = B.*colormax;
C = sprintf(['%4.2e',repmat([' \n%4.2e'], 1, ticknum)],B);
C = str2num(C);
caxis([colormin colormax])
colorbar('EastOutside','ytick',B,'yticklabel',C,...
'ticklength',[0.04 0.1],'YLim',[B(1) B(ticknum)],'FontSize',20)
% In quiverc2wcmap this loop plotted for each color level (n=20) and was very
% fast, but I found it did not plot some large data sets with enough color
% accuracy. Switched the loop to examine each data point individually.
% Takes longer but I'm more confident that the colors are correct. Note:
% much of this overhaul was inspired by DS's vfield_color.
hold on;
cmap = jet(64); %toggle type of colormap
CC = colormap(cmap);
cm_stepsize = (colormax-colormin)/length(CC);
for it=1:size(xyuvvrNNs,1) %takes ~13 seconds for a 8730-point dataset
% Some quiverc2wcmap fragments for reference:
% for it=1:n %10x faster, but may not be accurate
% c = CC(iCCs(it),:); %colormap([1:64](it),:) %ie "This RGB color row ="
% si = iData(it)+1; %[1:size(data)](it)+1; %ie "This start row"
% ei = iData(it+1); %[1:size(data)](it+1); %ie "This end row"
% hh=quiver(xyuvvrNNs(si:ei,1),xyuvvrNNs(si:ei,2),...
% xyuvvrNNs(si:ei,3),xyuvvrNNs(si:ei,4),scale*it/n,'Color',c)
cm_index = floor( (xyuvvrNNs(it,5) - colormin) / ( cm_stepsize ) ) + 1;
if cm_index == 1 %in case colormin is zero
c = CC(cm_index,:);
elseif cm_index > length(CC) %in case max(xyuvvrNNs) > colormax
cm_index = length(CC);
c = CC(cm_index,:);
elseif cm_index <= 0 %in case min(xyuvvrNNs) < colormin
cm_index = 1;
c = CC(cm_index,:);
else
cm_index = cm_index-1;
c = CC(cm_index,:);
end
hh=quiver(xyuvvrNNs(it,1),xyuvvrNNs(it,2),...
xyuvvrNNs(it,3),xyuvvrNNs(it,4),scale,'Color',c);
end
%----------Rest of document is from quiverc2wcmap---------------
% This is Matlab's 6.5.1 quiver. I figure that ensures a fair amouint of backward
% compatibility but I needed to hack it to allow a 'Color' property. Obviously
% a person could do more.
function hh = quiver(varargin)
%QUIVER Quiver plot.
% QUIVER(X,Y,U,V) plots velocity vectors as arrows with components (u,v)
% at the points (x,y). The matrices X,Y,U,V must all be the same size
% and contain corresponding position and velocity components (X and Y
% can also be vectors to specify a uniform grid). QUIVER automatically
% scales the arrows to fit within the grid.
%
% QUIVER(U,V) plots velocity vectors at equally spaced points in
% the x-y plane.
%
% QUIVER(U,V,S) or QUIVER(X,Y,U,V,S) automatically scales the
% arrows to fit within the grid and then stretches them by S. Use
% S=0 to plot the arrows without the automatic scaling.
%
% QUIVER(...,LINESPEC) uses the plot linestyle specified for
% the velocity vectors. Any marker in LINESPEC is drawn at the base
% instead of an arrow on the tip. Use a marker of '.' to specify
% no marker at all. See PLOT for other possibilities.
%
% QUIVER(...,'filled') fills any markers specified.
%
% H = QUIVER(...) returns a vector of line handles.
%
% Example:
% [x,y] = meshgrid(-2:.2:2,-1:.15:1);
% z = x .* exp(-x.^2 - y.^2); [px,py] = gradient(z,.2,.15);
% contour(x,y,z), hold on
% quiver(x,y,px,py), hold off, axis image
%
% See also FEATHER, QUIVER3, PLOT.
% Clay M. Thompson 3-3-94
% Copyright 1984-2002 The MathWorks, Inc.
% $Revision: 5.21 $ $Date: 2002/06/05 20:05:16 $
% Arrow head parameters
alpha = 0.33; % Size of arrow head relative to the length of the vector
beta = 0.33; % Width of the base of the arrow head relative to the length
autoscale = 1; % Autoscale if ~= 0 then scale by this.
plotarrows = 1; % Plot arrows
sym = '';
filled = 0;
ls = '-';
ms = '';
col = 'b';
nin = nargin;
ColorSpecInd = find(strcmpi(varargin, 'color'));
if(length(ColorSpecInd)==1 & nin > ColorSpecInd)
col = varargin{ColorSpecInd+1};
varargin = varargin([1:ColorSpecInd-1,ColorSpecInd+2:nin]);
nin = nin-2;
end
% Parse the string inputs
while isstr(varargin{nin}),
vv = varargin{nin};
if ~isempty(vv) & strcmp(lower(vv(1)),'f')
filled = 1;
nin = nin-1;
else
[l,c,m,msg] = colstyle(vv);
if ~isempty(msg),
error(sprintf('Unknown option "%s".',vv));
end
if ~isempty(l), ls = l; end
if ~isempty(c), col = c; end
if ~isempty(m), ms = m; plotarrows = 0; end
if isequal(m,'.'), ms = ''; end % Don't plot '.'
nin = nin-1;
end
end
error(nargchk(2,5,nin));
% Check numeric input arguments
if nin<4, % quiver(u,v) or quiver(u,v,s)
[msg,x,y,u,v] = xyzchk(varargin{1:2});
else
[msg,x,y,u,v] = xyzchk(varargin{1:4});
end
if ~isempty(msg), error(msg); end
if nin==3 | nin==5, % quiver(u,v,s) or quiver(x,y,u,v,s)
autoscale = varargin{nin};
end
% Scalar expand u,v
if prod(size(u))==1, u = u(ones(size(x))); end
if prod(size(v))==1, v = v(ones(size(u))); end
if autoscale,
% Base autoscale value on average spacing in the x and y
% directions. Estimate number of points in each direction as
% either the size of the input arrays or the effective square
% spacing if x and y are vectors.
if min(size(x))==1, n=sqrt(prod(size(x))); m=n; else [m,n]=size(x); end
delx = diff([min(x(:)) max(x(:))])/n;
dely = diff([min(y(:)) max(y(:))])/m;
del = delx.^2 + dely.^2;
if del>0
len = sqrt((u.^2 + v.^2)/del);
maxlen = max(len(:));
else
maxlen = 0;
end