Commit eb4d132a897e8136d70c3c776f288bfc5de24f0f

Authored by Jean-Luc Shaw
1 parent 48d36084
Exists in master

Added m_map, and libjl dependencies.

lib/insertAt.m 0 → 100644
... ... @@ -0,0 +1,25 @@
  1 +function output = insertAt(input,index,value)
  2 +%
  3 +% output = insertAt(input,index,value)
  4 +%
  5 +% this function returns the 'input' vector with 'value' inserted
  6 +% at position 'index' . Picture it as 'bumping up' the values at
  7 +% each of the index positions.
  8 +%
  9 +
  10 +N = numel(input) ;
  11 +Ni = numel(index) ;
  12 +Nv = numel(value) ;
  13 +
  14 +index = sort(index,'ascend') ;
  15 +c = false(N+Nv*Ni,1);
  16 +indexvec = [] ;
  17 +for ii = 1:Ni
  18 + indexvec = [indexvec index(ii)+(ii-1)*Nv:1:index(ii)+(ii-1)*Nv+Nv-1] ;
  19 +end
  20 +c(indexvec) = true;
  21 +output = nan(size(c));
  22 +output(~c) = input;
  23 +output(c) = value ;
  24 +
  25 +end
... ...
lib/llt2spd.m 0 → 100644
... ... @@ -0,0 +1,49 @@
  1 +function [u,v,spd,head] = llt2spd(lon,lat,time)
  2 +
  3 +% Synthax : [u,v,spd,head] = llt2spd(lon,lat,time)
  4 +%
  5 +% Returns the speed 'spd' and it's eastward/northward components 'u' and
  6 +% 'v' in (m/s), and the heading 'head' in degrees with 0 as east increasing
  7 +% towards the north (90 degrees).
  8 +%
  9 +% Values are calculated from a latitude/longitude time series with
  10 +% coordinates in decimal degrees and time in days.
  11 +%
  12 +% Values are estimated between longitude/latitude coordinates and then
  13 +% interpolated linearly to the input 'time' vector. Begining and end
  14 +% values are extrapolated using the "pchip" method.
  15 +
  16 +%% make everything column vectors
  17 +if ~iscolumn(lon) ; lon = lon' ; end
  18 +if ~iscolumn(lat) ; lat = lat' ; end
  19 +if ~iscolumn(time) ; time = time' ; end
  20 +
  21 +%% make differential vectors
  22 +dlon = diff(lon) ;
  23 +dlat = diff(lat) ;
  24 +dt = diff(time) ;
  25 +dl = m_lldist(lon,lat)*1000 ; % (m)
  26 +
  27 +%% calculate norm of speed
  28 +spd = dl./(dt*24*60*60) ;
  29 +
  30 +%% calculate heading
  31 +head = range_pm180_2_360(atan2d(dlat,dlon)) ;
  32 +
  33 +%% decompose into eastward/northward
  34 +u = spd.*cosd(head) ;
  35 +v = spd.*sind(head) ;
  36 +
  37 +%% make the time vector
  38 +t = time(1:end-1)+0.5*dt(1);
  39 +
  40 +%% interpolate to the input time grid
  41 +u = interp1(t,u,time,'linear','extrap') ;
  42 +v = interp1(t,v,time,'linear','extrap') ;
  43 +spd = interp1(t,spd,time) ;
  44 +spd(1) = sqrt(u(1).^2+v(1).^2);
  45 +spd(end) = sqrt(u(end).^2+v(end).^2);
  46 +head = interp1(t,head,time) ;
  47 +
  48 +
  49 +end
... ...
lib/m_lldist.m 0 → 100644
... ... @@ -0,0 +1,65 @@
  1 +function [dist,lons,lats] = m_lldist(long,lat,N)
  2 +% M_LLDIST Spherical earth distance between points in long/lat coordinates.
  3 +% RANGE=M_LLDIST(LONG,LAT) gives the distance in kilometers between
  4 +% successive points in the vectors LONG and LAT, computed
  5 +% using the Haversine formula on a spherical earth of radius
  6 +% 6378.137km. Distances are probably good to better than 1% of the
  7 +% "true" distance on the ellipsoidal earth
  8 +%
  9 +% [RANGE,LONGS,LATS]=M_LLDIST(LONG,LAT,N) computes the N-point geodesics
  10 +% between successive points. Each geodesic is returned on its
  11 +% own row of length N+1.
  12 +%
  13 +% See also M_XYDIST
  14 +
  15 +% Rich Pawlowicz (rich@ocgy.ubc.ca) 6/Nov/00
  16 +% This software is provided "as is" without warranty of any kind. But
  17 +% it's mine, so you can't sell it.
  18 +%
  19 +% 30/Dec/2005 - added n-point geodesic computations, based on an algorithm
  20 +% coded by Jeff Barton at Johns Hopkins APL in an m-file
  21 +% I looked at at mathworks.com.
  22 +
  23 +
  24 +pi180=pi/180;
  25 +earth_radius=6378.137;
  26 +
  27 +m=length(long)-1;
  28 +
  29 +long1=reshape(long(1:end-1),m,1)*pi180;
  30 +long2=reshape(long(2:end) ,m,1)*pi180;
  31 +lat1= reshape(lat(1:end-1) ,m,1)*pi180;
  32 +lat2= reshape(lat(2:end) ,m,1)*pi180;
  33 +
  34 +dlon = long2 - long1;
  35 +dlat = lat2 - lat1;
  36 +a = (sin(dlat/2)).^2 + cos(lat1) .* cos(lat2) .* (sin(dlon/2)).^2;
  37 +angles = 2 * atan2( sqrt(a), sqrt(1-a) );
  38 +dist = earth_radius * angles;
  39 +
  40 +
  41 +if nargin==3 & nargout>1, % Compute geodesics.
  42 +
  43 + % Cartesian unit vectors in rows of v1,v2
  44 + v1=[cos(long1).*cos(lat1) sin(long1).*cos(lat1) sin(lat1) ];
  45 + v2=[cos(long2).*cos(lat2) sin(long2).*cos(lat2) sin(lat2) ];
  46 +
  47 + % We want to get a unit vector tangent to the great circle.
  48 + n1=cross(v1,v2,2);
  49 + t1=cross(n1,v1,2);
  50 + t1=t1./repmat(sqrt(sum(t1.^2,2)),1,3);
  51 +
  52 + lons=zeros(m,N+1);
  53 + lats=zeros(m,N+1);
  54 + for k=1:m,
  55 +
  56 + % Radials for all points
  57 + p1=v1(k,:)'*cos(angles(k)*[0:N]/N) + t1(k,:)'*sin(angles(k)*[0:N]/N);
  58 +
  59 + lons(k,:)=atan2(p1(2,:),p1(1,:))/pi180;
  60 + lats(k,:)=asin(p1(3,:))/pi180;
  61 +
  62 + end;
  63 +
  64 +end;
  65 +
... ...
lib/range_pm180_2_360.m 0 → 100644
... ... @@ -0,0 +1,8 @@
  1 +function output = range_pm180_2_360(input)
  2 +
  3 +I = find(input < 0) ;
  4 +input(I) = input(I) + 360 ;
  5 +
  6 +output = input ;
  7 +
  8 +end
0 9 \ No newline at end of file
... ...