write_vpv

 WRITE_VPV:  Writes velocity profile data in ASCII VPV format.
  VPV is the freely-available "USGS Velocity Profile Viewer,"
    an Open-GL application for viewing velocity profiles from 
    ADCPs, models.  See <http://wwwdcascr.wr.usgs.gov/projects/vpv>.

  Usage:  vpv_format(VPV_FILE, T, U, V) 
   
  Inputs:  
     VPV_FILE = string containing name for output VPV file.
                (If VPV_FILE is empty or contains a '*' wildcard, 
                 the "uiputfile" dialog is invoked. )

     T = time vector in Matlab's "serial date format" (output of DATENUM)

     U = u velocity component matrix (cm/s): each row is a time record,
            and each column is a time series for each bin (or z level).
            The first row is the bottom bin. 

     V = v velocity component matrix (cm/s): same format as u.

0001 function write_vpv(theFilename, t, u, v)
0002 % WRITE_VPV:  Writes velocity profile data in ASCII VPV format.
0003 %  VPV is the freely-available "USGS Velocity Profile Viewer,"
0004 %    an Open-GL application for viewing velocity profiles from
0005 %    ADCPs, models.  See <http://wwwdcascr.wr.usgs.gov/projects/vpv>.
0006 %
0007 %  Usage:  vpv_format(VPV_FILE, T, U, V)
0008 %
0009 %  Inputs:
0010 %     VPV_FILE = string containing name for output VPV file.
0011 %                (If VPV_FILE is empty or contains a '*' wildcard,
0012 %                 the "uiputfile" dialog is invoked. )
0013 %
0014 %     T = time vector in Matlab's "serial date format" (output of DATENUM)
0015 %
0016 %     U = u velocity component matrix (cm/s): each row is a time record,
0017 %            and each column is a time series for each bin (or z level).
0018 %            The first row is the bottom bin.
0019 %
0020 %     V = v velocity component matrix (cm/s): same format as u.
0021  
0022 % Copyright (C) 1998 Dr. Charles R. Denham, ZYDECO.
0023 %  All Rights Reserved.
0024 %   Disclosure without explicit written consent from the
0025 %    copyright owner does not constitute publication.
0026  
0027 % Version of 16-Nov-1998 14:24:25.
0028 
0029 fclose('all');
0030 
0031 if nargin < 1, help(mfilename), return, end
0032 
0033 if isempty(theFilename), theFilename = '*'; end
0034 
0035 if any(theFilename == '*')
0036         theSuggested = 'vpv.out';
0037         [theFile, thePath] = uiputfile(theSuggested, 'Save VSV Data');
0038         if ~any(theFile), return, end
0039         if thePath(end) ~= filesep, thePath = [thePath filesep]; end
0040         theFilename = [thePath theFile];
0041 end
0042 
0043 if nargin == 3
0044         if size(u, 2) == 2   % Two columns as [u v].
0045                 v = u(:, 2);
0046                 u = u(:, 1);
0047         else   % Complex u + i*v.
0048                 v = imag(u);
0049                 u = real(u);
0050         end
0051 end
0052 
0053 u(u == -999999) = nan;   % WHFC NetCDF ADCP data.
0054 v(v == -999999) = nan;
0055 
0056 u(isnan(u)) = 9999.99;
0057 v(isnan(v)) = 9999.99;
0058 
0059 startDate = t(1);
0060 stopDate = t(end);
0061 
0062 [mEnsembles, nBins] = size(u);
0063 
0064 fmt1 = '%10.5f%10.5f%5i%5i%5i';
0065 fmt2 = ' %s\n';
0066 fmt3 = '%6i%10.2f%10.2f\n';
0067 
0068 fp = fopen(theFilename, 'w');
0069 if fp < 0, close(f), return, end
0070 
0071 fprintf(fp, '%s\n', theFilename);
0072 
0073 yd0 = yearday(startDate);
0074 ydn = yearday(stopDate);
0075 
0076 fprintf(fp, '%03i%4i%03i%4i\n', fix([yd0(2) yd0(1) ydn(2) ydn(1)]));
0077 
0078 lat = [0 0 0];
0079 lon = [0 0 0];
0080 
0081 fprintf(fp, '%2i%2i%2i%2i%2i%2i\n', [lat lon]);
0082 
0083 dt = t(2) - t(1);
0084 binHeight = 1;
0085 
0086 fprintf(fp, '%5i%10i%15.5f%10.5f\n', [1 nBins dt binHeight]);
0087 
0088 d0 = rem(datevec(startDate), 100);
0089 dn = rem(datevec(stopDate), 100);
0090 
0091 fprintf(fp, '%5i%5i%5i%10i%5i%5i\n', [d0(1:3) dn(1:3)]);
0092 
0093 tic
0094 try
0095         for i = 1:mEnsembles
0096                 remaining = mEnsembles-i+1;
0097                 if i == 2 | (i > 1 & rem(remaining, 100) == 0)
0098                         done = mEnsembles - remaining;
0099                         est = (remaining * toc / done) / 86400;
0100                         disp([' ## Remaining: ' int2str(remaining) ' records; ' ... 
0101    dhms(est)]);
0102                 end
0103                 d = datevec(t(i));
0104                 d(1) = rem(d(1), 100);
0105                 hm = int2str(100*d(4)+round(d(5)*100/60));
0106                 while length(hm) < 4, hm = ['0' hm]; end
0107                 dt = t(i)-t(1);
0108                 pressure = 0;
0109                 ensemble_header = [dt pressure d(1:3)];
0110                 fprintf(fp, fmt1, ensemble_header);
0111                 fprintf(fp, fmt2, hm);
0112                 for j = 1:nBins
0113                         ensemble_data = [j u(i, j) v(i, j)];
0114                         fprintf(fp, fmt3, ensemble_data);
0115                 end
0116         end
0117 catch
0118         disp([' ## Error: ' lasterr])
0119 end
0120 
0121 fclose(fp);
0122 
0123 
0124 function theResult = yearday(theYear, theDay)
0125 
0126 % yearday -- Convert date to year and day-of-year.
0127 %  yearday(theDate) returns [year day], where the day
0128 %   is a decimal date-number.  The given date can be
0129 %   a Matlab datenum, datestr, or datevec.
0130 %  yearday([theYear theDay]) returns the Matlab datenum
0131 %   corresponding to theYear and theDay.
0132 %  yearday(theYear, theDay) same as above.
0133 %  yearday (no arguments) demonstrates itself by showing
0134 %   a round-trip, using "now".
0135  
0136 % Copyright (C) 1998 Dr. Charles R. Denham, ZYDECO.
0137 %  All Rights Reserved.
0138 %   Disclosure without explicit written consent from the
0139 %    copyright owner does not constitute publication.
0140  
0141 % Version of 09-Nov-1998 08:17:43.
0142 
0143 if nargin < 1, theYear = 'demo'; end
0144 
0145 if isequal(theYear, 'demo')
0146         help(mfilename)
0147         a = now;
0148         disp(datestr(a))
0149         b = yearday(yearday(a));
0150         disp(datestr(b))
0151         year_day_round_trip_error = b-a
0152         return
0153 end
0154 
0155 % Two arguments: year and day ==> datenum.
0156 
0157 if nargin == 2, theYear = [theYear theDay]; end
0158 
0159 if length(theYear) == 2
0160         theDay = theYear(2);
0161         theYear = theYear(1);
0162         d = [theYear 1 1 0 0 0];   % January 1, midnight.
0163         for i = 1:6
0164                 v{i} = d(i);
0165         end
0166         result = datenum(v{:}) + theDay - 1;
0167         if nargout > 0
0168                 theResult = result;
0169         else
0170                 disp(result)
0171         end
0172         return
0173 end
0174 
0175 % One argument: date ==> year and day.
0176 
0177 theDate = theYear;
0178 
0179 if ischar(theDate)
0180         theDate = datenum(theDate);
0181 elseif length(theDate) > 1
0182         for i = 1:length(theDate)
0183                 v{i} = theDate(i);
0184         end
0185         theDate = datenum(v{:});
0186 end
0187 
0188 d = datevec(theDate);
0189 d(2:6) = [1 1 0 0 0];   % January 1, midnight.
0190 for i = 1:length(d)
0191         v{i} = d(i);
0192 end
0193 newYearsDay = datenum(v{:});
0194 
0195 delta = (theDate-newYearsDay);
0196 
0197 result = [d(1) (1+delta)];
0198 
0199 if nargout > 0
0200         theResult = result;
0201 else
0202         disp(result)
0203 end
0204 
0205 
0206 function theResult = dhms(theArg)
0207 
0208 % dhms -- Convert to/from dhms time-string format.
0209 %  dhms('theTimeString') converts 'theTimeString' of
0210 %   the form '1d2h3m4s' to decimal days.  Multiply
0211 %   the result by 86400 to get decimal seconds.
0212 %  dhms(theDecimalDays) converts theDecimalDays into a
0213 %   time-string of the format '1d2h3m4s', rounded to the
0214 %   nearest second.
0215 %  dhms (no argument) demonstrates itself.
0216  
0217 % Copyright (C) 1998 Dr. Charles R. Denham, ZYDECO.
0218 %  All Rights Reserved.
0219 %   Disclosure without explicit written consent from the
0220 %    copyright owner does not constitute publication.
0221  
0222 % Version of 16-Oct-1998 18:09:35.
0223 
0224 if nargin < 1, theArg = 'demo'; end
0225 
0226 if isequal(theArg, 'demo')
0227         help(mfilename)
0228         format long
0229         theTimeString = '9d23h59m59s'
0230         theDecimalDays = dhms2d(theTimeString)
0231         theTimeString = d2dhms(theDecimalDays)
0232         format short
0233         return
0234 end
0235 
0236 switch class(theArg)
0237 case 'double'
0238         result = d2dhms(theArg);
0239 case 'char'
0240         result = dhms2d(theArg);
0241 otherwise
0242         help(mfilename)
0243         warning([' ## Invalid argument type: ' class(theArg)])
0244         result = [];
0245 end
0246 
0247 if nargout > 0
0248         theResult = result;
0249 else
0250         disp(result)
0251 end
0252 
0253 
0254 function theResult = d2dhms(theDecimalDays)
0255 
0256 % d2dhms -- Convert decimal days to '1d2h3m4s' time-string.
0257 %  d2dhms(theDecimalDays) converts theDecimalDays into a
0258 %   time-string of the format '1d2h3m4s', rounded to the
0259 %   nearest second.
0260 %  d2dhms (no argument) demonstrates itself.
0261 %
0262 % Also see: dhms2d.
0263  
0264 % Copyright (C) 1998 Dr. Charles R. Denham, ZYDECO.
0265 %  All Rights Reserved.
0266 %   Disclosure without explicit written consent from the
0267 %    copyright owner does not constitute publication.
0268  
0269 % Version of 07-Oct-1998 08:16:05.
0270 
0271 if nargin < 1
0272         help(mfilename)
0273         theDecimalDays = 'demo';
0274 end
0275 
0276 if isequal(theDecimalDays, 'demo')
0277         theDecimalDays = dhms2d('9d23h59m59s')
0278         theTimeString = d2dhms(theDecimalDays)
0279         return
0280 end
0281 
0282 t = theDecimalDays;
0283 
0284 d = fix(t);
0285 t = rem(t, 1) * 24;
0286 h = fix(t);
0287 t = rem(t, 1) * 60;
0288 m = fix(t);
0289 t = rem(t, 1) * 60;
0290 s = round(t);
0291 
0292 result = '';
0293 if any(d), result = [result int2str(d) 'd']; end
0294 if any(h), result = [result int2str(h) 'h']; end
0295 if any(m), result = [result int2str(m) 'm']; end
0296 if any(s), result = [result int2str(s) 's']; end
0297 
0298 if isempty(result), result = '0s'; end
0299 
0300 if nargout > 0
0301         theResult = result;
0302 else
0303         disp(result)
0304 end
0305 
0306 
0307 function theResult = dhms2d(theTimeString)
0308 
0309 % dhms2d -- Convert dhms time-string to decimal days.
0310 %  dhms('theTimeString') converts 'theTimeString' of
0311 %   the form '1d2h3m4s' to decimal days.  Multiply
0312 %   the result by 86400 to get decimal seconds.
0313 %  dhms2d (no argument) demonstrates itself.
0314 %
0315 % Also see: d2dhms.
0316  
0317 % Copyright (C) 1998 Dr. Charles R. Denham, ZYDECO.
0318 %  All Rights Reserved.
0319 %   Disclosure without explicit written consent from the
0320 %    copyright owner does not constitute publication.
0321  
0322 % Version of 06-Oct-1998 21:19:25.
0323 
0324 if nargin < 1
0325         help(mfilename)
0326         theTimeString = 'demo';
0327 end
0328 
0329 if isequal(theTimeString, 'demo')
0330         theTimeString = '9d23h59m60s'
0331         theDecimalDays = dhms2d(theTimeString)
0332         return
0333 end
0334 
0335 t = lower(theTimeString);
0336 
0337 flag = 0;
0338 
0339 if any(t == 's')
0340         t = strrep(t, 's', '*1');
0341         flag = 1;
0342 end
0343 
0344 if any(t == 'm')
0345         if flag
0346                 t = strrep(t, 'm', '*60+');
0347         else
0348                 t = strrep(t, 'm', '*60');
0349         end
0350         flag = 1;
0351 end
0352 
0353 if any(t == 'h')
0354         if flag
0355                 t = strrep(t, 'h', '*3600+');
0356         else
0357                 t = strrep(t, 'h', '*3600');
0358         end
0359         flag = 1;
0360 end
0361 
0362 if any(t == 'd')
0363         if flag
0364                 t = strrep(t, 'd', '*86400+');
0365         else
0366                 t = strrep(t, 'd', '*86400');
0367         end
0368 end
0369 
0370 result = eval(t) / 86400;
0371 
0372 if nargout > 0
0373         theResult = result;
0374 else
0375         disp(result)
0376 end