timour_matlab/mscripts/DIAGWnew.m

 t0=input('Enter temperature : ')
% t0=27;

% Interpolate Dw on T=T0 surface

h=zeros(Nx,Ny);
hh=zeros(Nx,Ny);
hhh=zeros(Nx,Ny);
Dh=zeros(Nx,Ny);

meantheta(:,:,Nz)=0.;
for i=1:Nx
for j=1:Ny
kk=find(meantheta(i,j,:)<t0);

if kk(1)>1
h(i,j)=(kk(1)-1)*dz+dz*(meantheta(i,j,kk(1)-1)-t0)/(meantheta(i,j,kk(1)-1)-meantheta(i,j,kk(1)));

Dh(i,j)=Dw(i,j,kk(1)-1) ...
+(Dw(i,j,kk(1))-Dw(i,j,kk(1)-1)) ...
*(meantheta(i,j,kk(1)-1)-t0)/(meantheta(i,j,kk(1)-1)-meantheta(i,j,kk(1)));
else
h(i,j)=0;
Dh(i,j)=0;
end

end
end

Nz1=4;
Nx1=4;
Nx2=Nx-3;
Ny1=4;
Ny2=Ny-3;

he=dz*(Nz1-1);
hmax=max(max(h))


close all


we=0;
ws=0;
WE=zeros(Nx,Ny);
WS=zeros(Nx,Ny);
hh=zeros(Nx,Ny);
WE1=zeros(Nx,Ny,3);
WS1=zeros(Nx,Ny,3);
hh1=zeros(Nx,Ny,3);

for i=Nx1:Nx2
for j=Ny1:Ny2

WE(i,j)=meanw(i,j,Nz1);
WS(i,j)=Dh(i,j);

if h(i,j)>he
we=meanw(i,j,Nz1)*dx*dx+we;
ws=Dh(i,j)*dx*dx+ws;
end


end
end

we
ws

for i=Nx1:Nx2
for k=Nz1:Nz

z=dz*(k-1);
if z>he 
if z<h(i,Ny1)
we=we-meanv(i,Ny1,k)*dx*dz;
end
end

if z>he
if z<h(i,Ny2)
we=we+meanv(i,Ny2,k)*dx*dz;
end
end

end
end

for j=Ny1:Ny2
for k=Nz1:Nz

z=dz*(k-1);
if z>he
if z<h(Nx1,j)
we=we-meanu(Nx1,j,k)*dy*dz;
end
end

if z>he
if z<h(Nx2,j)
we=we+meanu(Nx2,j,k)*dy*dz;
end
end

end
end

we

'ratio', ws/we


hh=h;

for k=1:3
WE1(:,:,k)=WE(:,:);
WS1(:,:,k)=WS(:,:);
hh1(:,:,k)=hh(:,:);
end


for k=1:8
WE1=smooth3(WE1);
WS1=smooth3(WS1);
hh1=smooth3(hh1);
end

% hhh(2:Nx-1,2:Ny-1)=10000*(hh1(1:Nx-2,2:Ny-1,1)+hh1(3:Nx,2:Ny-1,1)+...
% hh1(2:Nx-1,1:Ny-2,1)+hh1(2:Nx-1,3:Ny,1)-4*hh1(2:Nx-1,2:Ny-1,1));


i=find(h<he);
h(i)=0;
% i=find(abs(hhh)>0.001);
% hhh(i)=0;

for i=1:Nx
for j=1:Ny

if h(i,j)<he
WE1(i,j,1)=0;
WS1(i,j,1)=NaN;
hh(i,j)=NaN;
hhh(i,j)=NaN;
end

end
end


hhh(1:Nx,:)=-10*hhh(1:Nx,:);
% hhh(Nx/2:Nx,:)=100*hhh(Nx/2:Nx,:);
a=max(max(hhh(1:Nx,:)));

%subplot(2,1,1)
%pcolor(hhh(:,:).');caxis([0 a]);shading flat; axis square
%set(gca,'DataAspectRatio',[2,2,2],'XtickLabel','|','YtickLabel','|')
%title('-\nabla^2 h','FontSize',15);
%xlabel('X','FontSize',15);
%ylabel('Y','Rotation',[0],'FontSize',13);

figure
WS1(1:Nx,:,1)=-WS1(1:Nx,:,1)*16667;
a=max(max(WS1(1:Nx,:,1)));
pcolor(squeeze(WS1(:,:,1)).');caxis([-0.3*a a]);shading flat; axis square;colorbar;
set(gca,'DataAspectRatio',[2,2,2],'XtickLabel','|','YtickLabel','|')
title('w^*','FontSize',15);
xlabel('X','FontSize',15);
ylabel('Y','Rotation',[0],'FontSize',13);


%WS2=squeeze(WS1(:,:,1));
%i=find(~isnan(hhh));
%WS2=WS2(i);
%WS=-WS*16667;
%hhh=hhh(i);
%sqrt(mean(WS2.*WS2)/ mean(hhh.*hhh))
%mean(WS2)/mean(hhh)
%max(WS2)/max(hhh)

1	t0=input('Enter temperature : ')
2	% t0=27;
3
4	% Interpolate Dw on T=T0 surface
5
6	h=zeros(Nx,Ny);
7	hh=zeros(Nx,Ny);
8	hhh=zeros(Nx,Ny);
9	Dh=zeros(Nx,Ny);
10
11	meantheta(:,:,Nz)=0.;
12	for i=1:Nx
13	for j=1:Ny
14	kk=find(meantheta(i,j,:)<t0);
15
16	if kk(1)>1
17	h(i,j)=(kk(1)-1)dz+dz(meantheta(i,j,kk(1)-1)-t0)/(meantheta(i,j,kk(1)-1)-meantheta(i,j,kk(1)));
18
19	Dh(i,j)=Dw(i,j,kk(1)-1) ...
20	+(Dw(i,j,kk(1))-Dw(i,j,kk(1)-1)) ...
21	*(meantheta(i,j,kk(1)-1)-t0)/(meantheta(i,j,kk(1)-1)-meantheta(i,j,kk(1)));
22	else
23	h(i,j)=0;
24	Dh(i,j)=0;
25	end
26
27	end
28	end
29
30	Nz1=4;
31	Nx1=4;
32	Nx2=Nx-3;
33	Ny1=4;
34	Ny2=Ny-3;
35
36	he=dz*(Nz1-1);
37	hmax=max(max(h))
38
39
40	close all
41
42
43	we=0;
44	ws=0;
45	WE=zeros(Nx,Ny);
46	WS=zeros(Nx,Ny);
47	hh=zeros(Nx,Ny);
48	WE1=zeros(Nx,Ny,3);
49	WS1=zeros(Nx,Ny,3);
50	hh1=zeros(Nx,Ny,3);
51
52	for i=Nx1:Nx2
53	for j=Ny1:Ny2
54
55	WE(i,j)=meanw(i,j,Nz1);
56	WS(i,j)=Dh(i,j);
57
58	if h(i,j)>he
59	we=meanw(i,j,Nz1)dxdx+we;
60	ws=Dh(i,j)dxdx+ws;
61	end
62
63
64
65	end
66	end
67
68	we
69	ws
70
71	for i=Nx1:Nx2
72	for k=Nz1:Nz
73
74	z=dz*(k-1);
75	if z>he
76	if z<h(i,Ny1)
77	we=we-meanv(i,Ny1,k)dxdz;
78	end
79	end
80
81	if z>he
82	if z<h(i,Ny2)
83	we=we+meanv(i,Ny2,k)dxdz;
84	end
85	end
86
87	end
88	end
89
90	for j=Ny1:Ny2
91	for k=Nz1:Nz
92
93	z=dz*(k-1);
94	if z>he
95	if z<h(Nx1,j)
96	we=we-meanu(Nx1,j,k)dydz;
97	end
98	end
99
100	if z>he
101	if z<h(Nx2,j)
102	we=we+meanu(Nx2,j,k)dydz;
103	end
104	end
105
106	end
107	end
108
109	we
110
111	'ratio', ws/we
112
113
114	hh=h;
115
116	for k=1:3
117	WE1(:,:,k)=WE(:,:);
118	WS1(:,:,k)=WS(:,:);
119	hh1(:,:,k)=hh(:,:);
120	end
121
122
123	for k=1:8
124	WE1=smooth3(WE1);
125	WS1=smooth3(WS1);
126	hh1=smooth3(hh1);
127	end
128
129	% hhh(2:Nx-1,2:Ny-1)=10000*(hh1(1:Nx-2,2:Ny-1,1)+hh1(3:Nx,2:Ny-1,1)+...
130	% hh1(2:Nx-1,1:Ny-2,1)+hh1(2:Nx-1,3:Ny,1)-4*hh1(2:Nx-1,2:Ny-1,1));
131
132
133	i=find(h<he);
134	h(i)=0;
135	% i=find(abs(hhh)>0.001);
136	% hhh(i)=0;
137
138	for i=1:Nx
139	for j=1:Ny
140
141	if h(i,j)<he
142	WE1(i,j,1)=0;
143	WS1(i,j,1)=NaN;
144	hh(i,j)=NaN;
145	hhh(i,j)=NaN;
146	end
147
148	end
149	end
150
151
152	hhh(1:Nx,:)=-10*hhh(1:Nx,:);
153	% hhh(Nx/2:Nx,:)=100*hhh(Nx/2:Nx,:);
154	a=max(max(hhh(1:Nx,:)));
155
156	%subplot(2,1,1)
157	%pcolor(hhh(:,:).');caxis([0 a]);shading flat; axis square
158	%set(gca,'DataAspectRatio',[2,2,2],'XtickLabel','\|','YtickLabel','\|')
159	%title('-\nabla^2 h','FontSize',15);
160	%xlabel('X','FontSize',15);
161	%ylabel('Y','Rotation',[0],'FontSize',13);
162
163	figure
164	WS1(1:Nx,:,1)=-WS1(1:Nx,:,1)*16667;
165	a=max(max(WS1(1:Nx,:,1)));
166	pcolor(squeeze(WS1(:,:,1)).');caxis([-0.3*a a]);shading flat; axis square;colorbar;
167	set(gca,'DataAspectRatio',[2,2,2],'XtickLabel','\|','YtickLabel','\|')
168	title('w^*','FontSize',15);
169	xlabel('X','FontSize',15);
170	ylabel('Y','Rotation',[0],'FontSize',13);
171
172
173	%WS2=squeeze(WS1(:,:,1));
174	%i=find(~isnan(hhh));
175	%WS2=WS2(i);
176	%WS=-WS*16667;
177	%hhh=hhh(i);
178	%sqrt(mean(WS2.WS2)/ mean(hhh.hhh))
179	%mean(WS2)/mean(hhh)
180	%max(WS2)/max(hhh)
181