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	edhill	1.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