Revision 495
(checked in by cbc, 12 years ago)
|
Initial import of Stark code.
|
Line | |
---|
1 |
% |
---|
2 |
% Load data from files |
---|
3 |
% |
---|
4 |
clear all |
---|
5 |
!ls *.m3d *.o3d |
---|
6 |
filename=input('Enter the name of observational file: ','s'); |
---|
7 |
[obs,gridname,year,ncol,gmt]=read_obs(filename); |
---|
8 |
% |
---|
9 |
% Filter data |
---|
10 |
% |
---|
11 |
size(obs); |
---|
12 |
nnv=21; |
---|
13 |
nn=ans(1)/21; |
---|
14 |
ii=0; |
---|
15 |
for i=1:10:nn |
---|
16 |
for j=1:nnv |
---|
17 |
ij=(i-1)*nnv+j; |
---|
18 |
ii=ii+1; |
---|
19 |
obs1(ii,:)=obs(ij,:); |
---|
20 |
end |
---|
21 |
end |
---|
22 |
clear obs; |
---|
23 |
obs=obs1; |
---|
24 |
% |
---|
25 |
% Enter mesh generation info |
---|
26 |
% |
---|
27 |
for i=1:ncol |
---|
28 |
fprintf(1,' first col %4.0f entry is %12.4e \n',i,obs(1,i)) |
---|
29 |
end |
---|
30 |
xcol=input('Enter the column number which corresponds to the x axis: '); |
---|
31 |
ycol=input('Enter the column number which corresponds to the y axis: '); |
---|
32 |
x=obs(:,xcol); |
---|
33 |
y=obs(:,ycol); |
---|
34 |
% |
---|
35 |
% Create mesh and create boundary |
---|
36 |
% |
---|
37 |
nnv=21; |
---|
38 |
size(x); |
---|
39 |
nn=ans(1) |
---|
40 |
nseg=nn/nnv-1 |
---|
41 |
femgen |
---|
42 |
bnd=detbndy(in); |
---|
43 |
% |
---|
44 |
% Begin plotting loop |
---|
45 |
% |
---|
46 |
newplot='y'; |
---|
47 |
while newplot == 'y', |
---|
48 |
% |
---|
49 |
% Plot boundary |
---|
50 |
% |
---|
51 |
figure |
---|
52 |
whitebg('w') |
---|
53 |
hold on |
---|
54 |
bndo=plotbnd(x,y,bnd); |
---|
55 |
set(bndo,'Color','k') |
---|
56 |
axis('off') |
---|
57 |
% |
---|
58 |
% Set contour levels and make plot |
---|
59 |
% |
---|
60 |
newcont='y'; |
---|
61 |
while newcont == 'y', |
---|
62 |
for i=1:ncol |
---|
63 |
fprintf(1,' first col %4.0f entry is %12.4e \n',i,obs(1,i)) |
---|
64 |
end |
---|
65 |
scol=input('Enter the column number which corresponds to the scalar you desire to contour: '); |
---|
66 |
scalar=obs(:,scol); |
---|
67 |
scrange(scalar) |
---|
68 |
cint=input('Enter the contour interval: '); |
---|
69 |
cmin=cint*ceil(min(scalar)/cint); |
---|
70 |
cmax=cint*floor(max(scalar)/cint); |
---|
71 |
clear cval |
---|
72 |
i=1; |
---|
73 |
cval(i)=cmin; |
---|
74 |
while cval(i) < cmax |
---|
75 |
i=i+1; |
---|
76 |
cval(i)=cval(i-1)+cint; |
---|
77 |
end |
---|
78 |
cval |
---|
79 |
hc=lcontour2(in,x,y,scalar,cval); |
---|
80 |
title(filename) |
---|
81 |
zoom on |
---|
82 |
newcont=input('New contour interval? (y/n): ','s'); |
---|
83 |
if newcont == 'y', |
---|
84 |
delete(hc); |
---|
85 |
end |
---|
86 |
end |
---|
87 |
newplot=input('New plot? (y/n): ','s'); |
---|
88 |
end |
---|