Newer
Older
leprob001
committed
function [velocityTotal, velocityLatitude, velocityLongitude] = horizontalVelocity(pLatitude, pLongitude, pInterval, pMethod, pNop)
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
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
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
% 29 October 2003, HHFurey. This routine does not account for
% discontinuous segments of track being passed, and therefore returns
% erroneous results. Not true - data passed in is bad - some segments that
% are discontinuous are passed through with NaNs as filler, but others are not!
% Hmmm ... working on r023 as a test case, this has three discontinous segments
% of track.
if nargin < 4
pMethod='average';
end
if nargin < 5
pNop=1;
end
%save debug_02
% nautic miles in kilometers
nm_in_km=60*1.852;
velocityTotal=-9.99*ones(size(pLatitude))*NaN;
velocityLatitude=-9.99*ones(size(pLatitude))*NaN;
velocityLongitude=-9.99*ones(size(pLatitude))*NaN;
if strcmp(pMethod,'average'); % centered mean between point ix+nop and ix-nop, centered at ix
for ix=1:length(pLatitude)
if ~((ix-pNop < 1) | (ix+pNop > length(pLatitude)))
if ~any(isnan(pLatitude([ix-pNop:ix+pNop])))
dkm = artoa.data.calculateGeodist([pLatitude(ix-pNop) pLongitude(ix-pNop)] , [pLatitude(ix+pNop) pLongitude(ix+pNop)]); %[m/s]
velocityTotal(ix)=(dkm*1000)/(pInterval*2*pNop);
velocityLatitude(ix)=((pLatitude(ix+pNop)-pLatitude(ix-pNop))*nm_in_km*1000)/(pInterval*2*pNop);
velocityLongitude(ix)=((pLongitude(ix+pNop)-pLongitude(ix-pNop))*nm_in_km*1000*cos(pLatitude(ix)/180*pi))/(pInterval*2*pNop);
end
end
end
elseif strcmp(pMethod,'backward'); % shifted means between points ix and ix-nop
for ix=1:length(pLatitude)
if ~(ix-pNop < 1)
if ~any(isnan(pLatitude([ix-pNop,ix])))
dkm=artoa.data.calculateGeodist([pLatitude(ix-pNop) pLongitude(ix-pNop)], [pLatitude(ix) pLongitude(ix)]);
velocityTotal(ix)=(dkm*1000)/(pInterval*pNop);
velocityLatitude(ix)=((pLatitude(ix)-pLatitude(ix-pNop))*nm_in_km*1000)/(pInterval*pNop);
velocityLongitude(ix)=((pLongitude(ix)-pLongitude(ix-pNop))*nm_in_km*1000*cos(pLatitude(ix)/180*pi))/(pInterval*pNop);
end
end
end
elseif strcmp(pMethod,'forward'); % shifted means between points ix and ix+nop
for ix=1:length(pLatitude)
if ~(ix+pNop > length(pLatitude))
if ~any(isnan(pLatitude([ix,ix+pNop])))
dkm=artoa.data.calculateGeodist([pLatitude(ix+pNop) pLongitude(ix+pNop)], [pLatitude(ix) pLongitude(ix)]);
velocityTotal(ix)=(dkm*1000)/(pInterval*pNop);
velocityLatitude(ix)=((pLatitude(ix+pNop)-pLatitude(ix))*nm_in_km*1000)/(pInterval*pNop);
velocityLongitude(ix)=((pLongitude(ix+pNop)-pLongitude(ix))*nm_in_km*1000*cos(pLatitude(ix)/180*pi))/(pInterval*pNop);
end
end
end
elseif strcmp(pMethod,'spline'); % spline at 2hour intervals
% spline latitude and longitudes and calculate velocities
% dummies must be excluded for the spline
resolution = 2; % this is the resolution at which spline will be performed
step = pInterval/3600/24; % interval lenght between successive data points [days]
time = 0:step:length(pLatitude)*step; % creates virtural time basis for raw data
time = time(:);
index = find(pLatitude<999); % find dummies, might be different value within artoa
pLatitude = pLatitude(index); pLongitude=pLongitude(index); time = time(index); % filter out dummies
% new timebase: every 2 hours including first and last time
if isempty(time)
dummy = 3;
end
time_min = time(1);
time_max = time(size(time,1)); % save for later
%time_spl = linspace(time_min, time_max,(time_max-time_min)*24/resolution+1 ); % time basis for spline [days]
%Corrected by T. Reynaud: 22/06/2005
time_spl = time_min:resolution/24:time_max;
time_spl = time_spl(:);
lat_spl_f = csape(time,pLatitude,[2 1]);
lat_spl = fnval(lat_spl_f,time_spl);
v_spl_f = fnder(lat_spl_f);
v_spl = fnval(v_spl_f,time_spl) *60*1.852*1000/(3600*24) ; % [m/s]
lon_spl_f = csape(time,pLongitude,[2 1]);
lon_spl = fnval(lon_spl_f,time_spl);
u_spl_f = fnder(lon_spl_f);
u_spl = fnval(u_spl_f,time_spl) .*cos(lat_spl/180*pi)*60*1.852*1000/(3600*24) ; % [m/s]
% subsample to original rate:
%save debug_tr
for i = 1:length(time)
% indi(i) = find(time_spl == time(i));
%Modified by T.Reynaud: 22/03/2005
indi(i) = find(abs(time_spl-time(i)) <= 0.000000000001);
end
%figure
%title('control plot for spline function')
%plot(lon_spl,lat_spl); hold on
%plot(lon,lat,'rx');
%zoom on;
velocityLongitude(index) = u_spl(indi);
velocityLatitude(index) = v_spl(indi);
velocityTotal(index) = sqrt(u_spl(indi).^2 + v_spl(indi).^2);
end
%%%%%%%%%%%%% idea for faster calculation of 'cetered mean'
%v = NaN*ones(size(data(:,7)))
%v(2:length(data)-1,1) = data(1:length(data)-2,6)
%v(2:length(data)-1,2) = data(3:length(data),6)
%v = diff(v')'
%v = v*60*1.852*1000*100/3600/24