clear all;close all
amp=5;x0=80;Rx=20;dx=1;y0=0;Ry=64/3;dy=1
nx=160;ny=64;
for j=1:ny;
    for i=1:nx;
        c(i,j)=0;
        for id=1:4
        c(i,j)=c(i,j)+amp*exp(-((i*dx-x0)/Rx)^2)*exp(-((j*dy-y0-(id-1)*3*Ry)/Ry)^2);
        end
    end
end
figure(2);
surf(c');
save 'h_init.dat' c -ascii;