-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathstraight_y_stroke_2.m
52 lines (52 loc) · 1.67 KB
/
straight_y_stroke_2.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
function PSI=straight_y_stroke_2(i,imax,jstart,jmax,PSIi0,PSIj0,X1,Y1,U1,V1)
PSIi=zeros(jmax,1);
PSIj=zeros(jmax,1);
PSI=zeros(jmax,1);
if jstart==2
PSIj(1)=PSIj0;
end
for j=jstart:jmax
if i==1
if isnan(PSIi0)==false
if length(PSIi0)==1
PSIi0(2:jmax,1)=PSIi0;
end
PSIi(j)=PSIi0(j) -V1(j,imax-i+1)*(X1(j,imax-i+1-1)-X1(j,imax-i+1));
if j==1
if isnan(PSIj0)==false
PSIj(j)=PSIj0 +U1(j,imax-i+1)*(Y1(j+1,imax-i+1)-Y1(j,imax-i+1));
end
else
PSIj(j)=PSI(j-1)+...
0.5*(U1(j,imax-i+1)+U1(j-1,imax-i+1))*(Y1(j,imax-i+1)-Y1(j-1,imax-i+1));
end
if exist('PSIj(j)','var')==1
PSI(j)=0.5*(PSIi(j)+PSIj(j));
else
PSI(j)=PSIi(j);
end
else
if j==1
PSI(j)=PSIj0 +U1(j,imax-i+1)*(Y1(j+1,imax-i+1)-Y1(j,imax-i+1));
else
PSI(j)=PSI(j-1)+...
0.5*(U1(j,imax-i+1)+U1(j-1,imax-i+1))*(Y1(j,imax-i+1)-Y1(j-1,imax-i+1));
end
end
else
PSIi(j)=PSIi0(j) -0.5*(V1(j,imax-i+1)+V1(j,imax-i+1+1))*(X1(j,imax-i+1)-X1(j,imax-i+1+1));
if j==1
if isnan(PSIj0)==false
PSIj(j)=PSIj0 +U1(j,imax-i+1)*(Y1(j+1,imax-i+1)-Y1(j,imax-i+1));
end
else
PSIj(j)=PSIj(j-1)...
+0.5*(U1(j,imax-i+1)+U1(j-1,imax-i+1))*(Y1(j,imax-i+1)-Y1(j-1,imax-i+1));
end
if exist('PSIj(j)','var')==1
PSI(j)=0.5*(PSIi(j)+PSIj(j));
else
PSI(j)=PSIi(j);
end
end
end