inverted_barometer/input/gendata.m

ieee='b';
accuracy='real*8';

Ho=2000;
nx=60;
ny=60;

% Flat bottom at z=-Ho
h=-Ho*ones(nx,ny);
% Walls
h(end,:)=0;
h(:,end)=0;

% center point
xc = nx/2;
yc = ny/2;
[x y]  = meshgrid([1:nx]',[1:ny]');
r0 = nx/2;
r = sqrt((x-xc).^2+(y-yc).^2);
h(find(r.^2>r0.^2)) = 0;
% coordinate transformation
phi = acos((x-xc)./r);
phi(find(r==0))=0;

kr = 2*pi/r0;

%mask = change(h,'~=',0,1);
mask=ones(nx,ny); mask(find(h==0))=0;

fid=fopen('topog.box','w',ieee); fwrite(fid,h,accuracy); fclose(fid);

% Atmospheric pressure in (Pa)
pMean=0e5;
pMax = 1e4;
nu = 2;
slp = pMean + pMax*besselj(nu,kr*r).*cos(nu*phi).*mask;

% remove mean
slp = (slp-sum(slp(:))/sum(mask(:))).*mask;
fid=fopen('pLoad.bin','w',ieee); fwrite(fid,slp,accuracy); fclose(fid);
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	%mask = change(h,'~=',0,1);
28	mask=ones(nx,ny); mask(find(h==0))=0;
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 + pMaxbesselj(nu,krr).cos(nuphi).*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);