inverted_barometer/input/gendata.m

395b093796 Mart*0001 ieee='b';
                0002 accuracy='real*8';
                0003
                0004 Ho=2000;
                0005 nx=60;
                0006 ny=60;
                0007
                0008 % Flat bottom at z=-Ho
                0009 h=-Ho*ones(nx,ny);
                0010 % Walls
                0011 h(end,:)=0;
                0012 h(:,end)=0;
                0013
                0014 % center point
                0015 xc = nx/2;
                0016 yc = ny/2;
                0017 [x y]  = meshgrid([1:nx]',[1:ny]');
                0018 r0 = nx/2;
                0019 r = sqrt((x-xc).^2+(y-yc).^2);
                0020 h(find(r.^2>r0.^2)) = 0;
                0021 % coordinate transformation
                0022 phi = acos((x-xc)./r);
                0023 phi(find(r==0))=0;
                0024
                0025 kr = 2*pi/r0;
                0026
96d647b91a Jean*0027 %mask = change(h,'~=',0,1);
                0028 mask=ones(nx,ny); mask(find(h==0))=0;
395b093796 Mart*0029
                0030 fid=fopen('topog.box','w',ieee); fwrite(fid,h,accuracy); fclose(fid);
                0031
                0032 % Atmospheric pressure in (Pa)
                0033 pMean=0e5;
                0034 pMax = 1e4;
                0035 nu = 2;
                0036 slp = pMean + pMax*besselj(nu,kr*r).*cos(nu*phi).*mask;
                0037
                0038 % remove mean
                0039 slp = (slp-sum(slp(:))/sum(mask(:))).*mask;
                0040 fid=fopen('pLoad.bin','w',ieee); fwrite(fid,slp,accuracy); fclose(fid);