1 |
function R=timeribbon(icq4liststruct,mesh,xi,yi,t,interpint) |
---|
2 |
|
---|
3 |
% R=timeribbon(icq4liststruct,varname,mesh,[x y t],interpint) |
---|
4 |
% |
---|
5 |
% R = x: [nnv x npts] |
---|
6 |
% y: [nnv x npts] |
---|
7 |
% z: [nnv x npts] |
---|
8 |
% d: [nnv x npts] |
---|
9 |
% v: [nnv x npts] |
---|
10 |
% t: [nnv x npts] |
---|
11 |
% |
---|
12 |
% Q and Z are expected to be [nn x nnv] |
---|
13 |
|
---|
14 |
flds{1}='ZMID'; |
---|
15 |
flds{2}='UZMID'; |
---|
16 |
flds{3}='VZMID'; |
---|
17 |
flds{4}='WZMID'; |
---|
18 |
flds{5}='Q2MID'; |
---|
19 |
flds{6}='Q2LMID'; |
---|
20 |
flds{7}='TMPMID'; |
---|
21 |
flds{8}='SALMID'; |
---|
22 |
|
---|
23 |
if length(xi)==1 |
---|
24 |
xi=xi*ones(size(t)); |
---|
25 |
yi=yi*ones(size(t)); |
---|
26 |
end |
---|
27 |
if length(t)==1 |
---|
28 |
t=t*ones(size(xi)); |
---|
29 |
else |
---|
30 |
t=t(:); |
---|
31 |
end |
---|
32 |
|
---|
33 |
% OPNML Compliant 3-18-99 CVL |
---|
34 |
|
---|
35 |
mesh=belint(mesh); |
---|
36 |
mesh=el_areas(mesh); |
---|
37 |
|
---|
38 |
% INTERPOLATE X & Y ALONG TRACK IF REQUESTED |
---|
39 |
|
---|
40 |
d=[0; cumsum(((diff(xi(:))).^2+(diff(yi(:))).^2).^0.5)]; |
---|
41 |
xy=[xi(:),yi(:)]; |
---|
42 |
if exist('interpint') |
---|
43 |
if ~isnan(interpint) |
---|
44 |
di=d; |
---|
45 |
d=sort([d(:);[interpint:interpint:max(d)]']); |
---|
46 |
xy=interp1(di,xy,d); |
---|
47 |
t=interp1(di,t,d); |
---|
48 |
end |
---|
49 |
end |
---|
50 |
|
---|
51 |
% BUILD SPATIAL BASIS FN LIST |
---|
52 |
|
---|
53 |
ll=findelem(mesh,xy); |
---|
54 |
ind=find(~isnan(ll)); |
---|
55 |
if isempty(ind) |
---|
56 |
error('section off grid!') |
---|
57 |
end |
---|
58 |
p=basis2d(mesh,xy); |
---|
59 |
|
---|
60 |
% INTERPOLATE Z & Q AT EACH SIGMA LEVEL |
---|
61 |
|
---|
62 |
D=icq4liststruct; |
---|
63 |
NFILE=length(D); |
---|
64 |
for in=1:NFILE |
---|
65 |
f{1}=read_icq4(D(in).name,0); |
---|
66 |
jday(in)=julian([f{1}.year,f{1}.month,f{1}.day,... |
---|
67 |
0,0,f{1}.curr_seconds]); |
---|
68 |
end |
---|
69 |
[jday,Dsi]=sort(jday); |
---|
70 |
D=D(Dsi); |
---|
71 |
fileind=floor(interp1(jday,1:NFILE,t)); |
---|
72 |
fileind=[fileind,fileind+1]'; |
---|
73 |
if t(1)<jday(1) |
---|
74 |
fileind(1:2,find(t<jday(1)))=1; |
---|
75 |
end |
---|
76 |
if isnan(fileind(1,length(t))) |
---|
77 |
fileind(1:2,find(t>jday(NFILE)))=NFILE; |
---|
78 |
end |
---|
79 |
|
---|
80 |
% FILL STRUCTURE TO NPTS x NNV |
---|
81 |
|
---|
82 |
R.x=meshgrid(xy(:,1),1:f{1}.nnv)'; |
---|
83 |
R.y=meshgrid(xy(:,2),1:f{1}.nnv)'; |
---|
84 |
R.d=d*ones(1,f{1}.nnv); |
---|
85 |
R.t=meshgrid(t,1:f{1}.nnv)'; |
---|
86 |
for if=1:8 |
---|
87 |
R=setfield(R,flds{if},nan*R.x); |
---|
88 |
end |
---|
89 |
|
---|
90 |
% READ FIRST TWO FILES |
---|
91 |
f{1}=read_icq4(D(fileind(1,1)).name); |
---|
92 |
f{2}=read_icq4(D(fileind(2,1)).name); |
---|
93 |
|
---|
94 |
% FOR LOOP THROUGH TIME POINTS |
---|
95 |
|
---|
96 |
for it=ind(:)' |
---|
97 |
tf=jday(fileind(1:2,it)); |
---|
98 |
tb=min([max([(t(it)-tf(1))/(eps+diff(tf)) 0]) 1]); |
---|
99 |
R.tb(it,:)=[tb 1-tb]; |
---|
100 |
tb=[1-tb tb]; |
---|
101 |
|
---|
102 |
for if=1:8 |
---|
103 |
Q=tb(1)*getfield(f{1},flds{if})+tb(2)*getfield(f{2},flds{if}); |
---|
104 |
for iv=1:f{1}.nnv |
---|
105 |
R=setfield(R,flds{if},{it,iv},sum( ( p(it,:)' .* Q( mesh.e(ll(it),:),iv ) )' )); |
---|
106 |
end |
---|
107 |
end |
---|
108 |
|
---|
109 |
nit=min(length(t), it+1); |
---|
110 |
if fileind(1,it)~=fileind(1,nit) |
---|
111 |
if fileind(1,nit)==fileind(2,it) |
---|
112 |
f{1}=f{2}; |
---|
113 |
elseif fileind(1,nit)~=fileind(1,it) |
---|
114 |
disp(['reading 'D(fileind(1,nit)).name]) |
---|
115 |
f{1}=read_icq4(D(fileind(1,nit)).name); |
---|
116 |
end |
---|
117 |
end |
---|
118 |
if fileind(2,it)~=fileind(2,nit) |
---|
119 |
if fileind(2,nit)==fileind(1,nit) |
---|
120 |
f{2}=f{1}; |
---|
121 |
else |
---|
122 |
disp(['reading 'D(fileind(2,nit)).name]) |
---|
123 |
f{2}=read_icq4(D(fileind(2,nit)).name); |
---|
124 |
end |
---|
125 |
end |
---|
126 |
end |
---|
127 |
|
---|
128 |
|
---|