1 |
mlosch |
1.1 |
ieee='b'; |
2 |
|
|
accuracy='real*8'; |
3 |
|
|
|
4 |
|
|
Ho=2000; |
5 |
|
|
nx=60; |
6 |
|
|
ny=60; |
7 |
|
|
|
8 |
|
|
% Flat bottom at z=-Ho |
9 |
|
|
h=-Ho*ones(nx,ny); |
10 |
|
|
% Walls |
11 |
|
|
h(end,:)=0; |
12 |
|
|
h(:,end)=0; |
13 |
|
|
|
14 |
|
|
% center point |
15 |
|
|
xc = nx/2; |
16 |
|
|
yc = ny/2; |
17 |
|
|
[x y] = meshgrid([1:nx]',[1:ny]'); |
18 |
|
|
r0 = nx/2; |
19 |
|
|
r = sqrt((x-xc).^2+(y-yc).^2); |
20 |
|
|
h(find(r.^2>r0.^2)) = 0; |
21 |
|
|
% coordinate transformation |
22 |
|
|
phi = acos((x-xc)./r); |
23 |
|
|
phi(find(r==0))=0; |
24 |
|
|
|
25 |
|
|
kr = 2*pi/r0; |
26 |
|
|
|
27 |
|
|
|
28 |
|
|
mask = change(h,'~=',0,1); |
29 |
|
|
|
30 |
|
|
fid=fopen('topog.box','w',ieee); fwrite(fid,h,accuracy); fclose(fid); |
31 |
|
|
|
32 |
|
|
% Atmospheric pressure in (Pa) |
33 |
|
|
pMean=0e5; |
34 |
|
|
pMax = 1e4; |
35 |
|
|
nu = 2; |
36 |
|
|
slp = pMean + pMax*besselj(nu,kr*r).*cos(nu*phi).*mask; |
37 |
|
|
|
38 |
|
|
% remove mean |
39 |
|
|
slp = (slp-sum(slp(:))/sum(mask(:))).*mask; |
40 |
|
|
fid=fopen('pLoad.bin','w',ieee); fwrite(fid,slp,accuracy); fclose(fid); |