Subversion Repositories lagranto.um

Rev

Rev 3 | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 3 Rev 12
1 
function out = get_index (vert,position,grid);
 1 
function out = get_index (vert,position,grid);
2 
 
 2 
 
3 
% General index calculation for data grid
 3 
% General index calculation for data grid
4 
% Usage:   out = int_index(vert,coords,grid)
4 
% Usage:   out = int_index(vert,coords,grid)
5 
% Example: grid.xmin = -180;
 5 
% Example: grid.xmin = -180;
6 
%          grid.dx   =    1;
 6 
%          grid.dx   =    1;
7 
%          grid.nx   =  361;
 7 
%          grid.nx   =  361;
8 
%          grid.ymin =    0;
 8 
%          grid.ymin =    0;
9 
%          grid.dy   =    1;
 9 
%          grid.dy   =    1;
10 
%          grid.ny   =   91;
 10 
%          grid.ny   =   91;
11 
%          out get_index(inp.P.data,[ 500 45 78; 250 56.7 78.4],grid)
11 
%          out get_index(inp.P.data,[ 500 45 78; 250 56.7 78.4],grid)
12 
 
 12 
 
13 
% Get the position arrays
 13 
% Get the position arrays
14 
xpos = position(:,3)';
 14 
xpos = position(:,3)';
15 
ypos = position(:,2)';
 15 
ypos = position(:,2)';
16 
zpos = position(:,1)';
 16 
zpos = position(:,1)';
17 
 
 17 
 
18 
% Get the grid dimension
18 
% Get the grid dimension
19 
nx = size(vert,3);
 19 
nx = size(vert,3);
20 
ny = size(vert,2);
 20 
ny = size(vert,2);
21 
nz = size(vert,1);
 21 
nz = size(vert,1);
22 
 
 22 
 
23 
% Get the horiozontal grid indices
 23 
% Get the horiozontal grid indices
24 
rindx = (xpos-grid.xmin)/grid.dx+1;
 24 
rindx = (xpos-grid.xmin)/grid.dx+1;
25 
rindy = (ypos-grid.ymin)/grid.dy+1;
 25 
rindy = (ypos-grid.ymin)/grid.dy+1;
26 
 
 26 
 
27 
% Special treatment if the the <vert> field is two-dimensional
 27 
% Special treatment if the the <vert> field is two-dimensional
28 
if ( nz == 1 )
 28 
if ( nz == 1 )
29 
   out = [ ones(size(rindy))', rindy', rindx' ];
 29 
   out = [ ones(size(rindy))', rindy', rindx' ];
30 
   return;
 30 
   return;
31 
end
31 
end
32
32
33 
% Direction of the vertical axis (Pressure: descencding, Height: ascending)
 33 
% Direction of the vertical axis (Pressure: descencding, Height: ascending)
34 
dir = vert(2:nz,:,:) - vert(1:nz-1,:,:);
 34 
dir = vert(2:nz,:,:) - vert(1:nz-1,:,:);
35 
if ( mean(dir(:)) > 0 )
35 
if ( mean(dir(:)) > 0 )
36 
     zpos = -zpos;
 36 
     zpos = -zpos;
37 
     vert = -vert;
 37 
     vert = -vert;
38 
end
 38 
end
39 
 
 39 
 
40 
% Get the weights for the individual grid points
 40 
% Get the weights for the individual grid points
41 
fracx = 1-rindx+floor(rindx);
 41 
fracx = 1-rindx+floor(rindx);
42 
fracy = 1-rindy+floor(rindy);
 42 
fracy = 1-rindy+floor(rindy);
43 
 
 43 
 
44 
% Combine the indices i and j to a common index ij
 44 
% Combine the indices i and j to a common index ij
45 
vert1 = reshape(vert,nz,nx*ny);
 45 
vert1 = reshape(vert,nz,nx*ny);
46 
 
 46 
 
47 
% The position array and the vertical index is needed in the same format
 47 
% The position array and the vertical index is needed in the same format
48 
zpos1 = kron(ones(nz,1),zpos);
 48 
zpos1 = kron(ones(nz,1),zpos);
49 
zind1 = kron([1:nz]',ones(size(zpos)));
 49 
zind1 = kron([1:nz]',ones(size(zpos)));
50 
shif1 = (0:length(zpos)-1) * nz;
 50 
shif1 = (0:length(zpos)-1) * nz;
51 
 
 51 
 
52 
% Get the vertical indices at position 00
 52 
% Get the vertical indices at position 00
53 
indx           = floor(rindx);
 53 
indx           = floor(rindx);
54 
indy           = floor(rindy);
 54 
indy           = floor(rindy);
55 
indxy          = indy+(indx-1)*ny;
 55 
indxy          = indy+(indx-1)*ny;
56 
isnok1         = (indx<1) | (indx>nx) | (indy<1) | (indy>ny);
 56 
isnok1         = (indx<1) | (indx>nx) | (indy<1) | (indy>ny);
57 
indxy(isnok1)  = 1;
57 
indxy(isnok1)  = 1;
58 
vert           = vert1(:,indxy);
 58 
vert           = vert1(:,indxy);
59 
indz           = max((vert>=zpos1).*zind1);
 59 
indz           = max((vert>=zpos1).*zind1);
60 
isnok2         = (indz<1) | (indz>=nz);
 60 
isnok2         = (indz<1) | (indz>=nz);
61 
indz(isnok2)   = 1;
 61 
indz(isnok2)   = 1;
62 
indu           = indz     + shif1;
 62 
indu           = indz     + shif1;
63 
indo           = indz + 1 + shif1;
 63 
indo           = indz + 1 + shif1;
64 
diff           = vert(indo)-vert(indu);
 64 
diff           = vert(indo)-vert(indu);
65 
isnok3         = (diff==0);
 65 
isnok3         = (diff==0);
66 
diff(isnok3)   = 1;
 66 
diff(isnok3)   = 1;
67 
rind00         = indz + (zpos-vert(indu)) ./ diff;
 67 
rind00         = indz + (zpos-vert(indu)) ./ diff;
68 
isnok4         = isnok1 | isnok2 | isnok3;
 68 
isnok4         = isnok1 | isnok2 | isnok3;
69 
rind00(isnok4) = NaN;
 69 
rind00(isnok4) = NaN;
70 
 
 70 
 
71 
% Get the vertical indices at positions 01, 10, 11 (if necessary)
 71 
% Get the vertical indices at positions 01, 10, 11 (if necessary)
72 
if ( any(fracx ~= 1) | any(fracy ~=1) )
72 
if ( any(fracx ~= 1) | any(fracy ~=1) )
73
73
74 
    indx           = floor(rindx+1);
 74 
    indx           = floor(rindx+1);
75 
    indy           = floor(rindy);
 75 
    indy           = floor(rindy);
76 
    indxy          = indy+(indx-1)*ny;
 76 
    indxy          = indy+(indx-1)*ny;
77 
    isnok1         = (indx<1) | (indx>nx) | (indy<1) | (indy>ny);
 77 
    isnok1         = (indx<1) | (indx>nx) | (indy<1) | (indy>ny);
78 
    indxy(isnok1)  = 1;
78 
    indxy(isnok1)  = 1;
79 
    vert           = vert1(:,indxy);
 79 
    vert           = vert1(:,indxy);
80 
    indz           = max((vert>zpos1).*zind1);
 80 
    indz           = max((vert>zpos1).*zind1);
81 
    isnok2         = (indz<1) | (indz>=nz);
 81 
    isnok2         = (indz<1) | (indz>=nz);
82 
    indz(isnok2)   = 1;
 82 
    indz(isnok2)   = 1;
83 
    indu           = indz     + shif1;
83 
    indu           = indz     + shif1;
84 
    indo           = indz + 1 + shif1;
 84 
    indo           = indz + 1 + shif1;
85 
    diff           = vert(indo)-vert(indu);
 85 
    diff           = vert(indo)-vert(indu);
86 
    isnok3         = (diff==0);
 86 
    isnok3         = (diff==0);
87 
    diff(isnok3)   = 1;
 87 
    diff(isnok3)   = 1;
88 
    rind01         = indz + (zpos-vert(indu)) ./ diff;
 88 
    rind01         = indz + (zpos-vert(indu)) ./ diff;
89 
    isnok4         = isnok1 | isnok2 | isnok3;
 89 
    isnok4         = isnok1 | isnok2 | isnok3;
90 
    rind01(isnok4) = NaN;
 90 
    rind01(isnok4) = NaN;
91 
 
 91 
 
92 
    indx           = floor(rindx);
 92 
    indx           = floor(rindx);
93 
    indy           = floor(rindy+1);
 93 
    indy           = floor(rindy+1);
94 
    indxy          = indy+(indx-1)*ny;
 94 
    indxy          = indy+(indx-1)*ny;
95 
    isnok1         = (indx<1) | (indx>nx) | (indy<1) | (indy>ny);
 95 
    isnok1         = (indx<1) | (indx>nx) | (indy<1) | (indy>ny);
96 
    indxy(isnok1)  = 1;
96 
    indxy(isnok1)  = 1;
97 
    vert           = vert1(:,indxy);
 97 
    vert           = vert1(:,indxy);
98 
    indz           = max((vert>zpos1).*zind1);
 98 
    indz           = max((vert>zpos1).*zind1);
99 
    isnok2         = (indz<1) | (indz>=nz);
 99 
    isnok2         = (indz<1) | (indz>=nz);
100 
    indz(isnok2)   = 1;
 100 
    indz(isnok2)   = 1;
101 
    indu           = indz     + shif1;
101 
    indu           = indz     + shif1;
102 
    indo           = indz + 1 + shif1;
 102 
    indo           = indz + 1 + shif1;
103 
    diff           = vert(indo)-vert(indu);
 103 
    diff           = vert(indo)-vert(indu);
104 
    isnok3         = (diff==0);
 104 
    isnok3         = (diff==0);
105 
    diff(isnok3)   = 1;
 105 
    diff(isnok3)   = 1;
106 
    rind10         = indz + (zpos-vert(indu)) ./ diff;
 106 
    rind10         = indz + (zpos-vert(indu)) ./ diff;
107 
    isnok4         = isnok1 | isnok2 | isnok3;
 107 
    isnok4         = isnok1 | isnok2 | isnok3;
108 
    rind10(isnok4) = NaN;
 108 
    rind10(isnok4) = NaN;
109 
 
 109 
 
110 
    indx           = floor(rindx+1);
 110 
    indx           = floor(rindx+1);
111 
    indy           = floor(rindy+1);
 111 
    indy           = floor(rindy+1);
112 
    indxy          = indy+(indx-1)*ny;
 112 
    indxy          = indy+(indx-1)*ny;
113 
    isnok1         = (indx<1) | (indx>nx) | (indy<1) | (indy>ny);
 113 
    isnok1         = (indx<1) | (indx>nx) | (indy<1) | (indy>ny);
114 
    indxy(isnok1)  = 1;
114 
    indxy(isnok1)  = 1;
115 
    vert           = vert1(:,indxy);
 115 
    vert           = vert1(:,indxy);
116 
    indz           = max((vert>zpos1).*zind1);
 116 
    indz           = max((vert>zpos1).*zind1);
117 
    isnok2         = (indz<1) | (indz>=nz);
 117 
    isnok2         = (indz<1) | (indz>=nz);
118 
    indz(isnok2)   = 1;
 118 
    indz(isnok2)   = 1;
119 
    indu           = indz     + shif1;
119 
    indu           = indz     + shif1;
120 
    indo           = indz + 1 + shif1;
 120 
    indo           = indz + 1 + shif1;
121 
    diff           = vert(indo)-vert(indu);
 121 
    diff           = vert(indo)-vert(indu);
122 
    isnok3         = (diff==0);
 122 
    isnok3         = (diff==0);
123 
    diff(isnok3)   = 1;
 123 
    diff(isnok3)   = 1;
124 
    rind11         = indz + (zpos-vert(indu)) ./ diff;
 124 
    rind11         = indz + (zpos-vert(indu)) ./ diff;
125 
    isnok4         = isnok1 | isnok2 | isnok3;
 125 
    isnok4         = isnok1 | isnok2 | isnok3;
126 
    rind11(isnok4) = NaN;
 126 
    rind11(isnok4) = NaN;
127 
 
 127 
 
128 
else
 128 
else
129 
 
 129 
 
130 
    rind01 = rind00;
 130 
    rind01 = rind00;
131 
    rind10 = rind00;
 131 
    rind10 = rind00;
132 
    rind11 = rind00;
 132 
    rind11 = rind00;
133
133
134 
end
 134 
end
135
135
136 
% Vectorised interpolation
 136 
% Vectorised interpolation
137 
rindz = fracx     .* fracy     .* rind00 + ...
 137 
rindz = fracx     .* fracy     .* rind00 + ...
138 
        (1-fracx) .* fracy     .* rind01 + ...
 138 
        (1-fracx) .* fracy     .* rind01 + ...
139 
        fracx     .* (1-fracy) .* rind10 + ...
 139 
        fracx     .* (1-fracy) .* rind10 + ...
140 
        (1-fracx) .* (1-fracy) .* rind11;
 140 
        (1-fracx) .* (1-fracy) .* rind11;
141
141
142 
% Set the output of the function
 142 
% Set the output of the function
143 
nok = isnan(rindz') | isnan(rindy') | isnan(rindx');
 143 
nok = isnan(rindz') | isnan(rindy') | isnan(rindx');
144 
out = [ rindz', rindy', rindx' ];
 144 
out = [ rindz', rindy', rindx' ];
145 
 
 145 
 
146 
 
 146 
 
147 
 
 147