NCCOOS Trac Projects: Top | Web | Platforms | Processing | Viz | Sprints | Sandbox | (Wind)

root/gliderproc/trunk/MATLAB/opnml/FCAST_1.2/timeribbon.m

Revision 495 (checked in by cbc, 12 years ago)

Initial import of Stark code.

Line 
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
Note: See TracBrowser for help on using the browser.