
clear

N = 10^(-3);
f = 10^(-4);

theta0 = 0.05;
thetaprime = pi/2-theta0;

nx0 = 1;
nz0 = 25;
kx0 = 8*10^(-4);
k0 = kx0/sin(theta0);
kz0 = kx0/tan(theta0);
omega0 = sqrt((N*kx0)^2+(f*kz0)^2)/k0;

%--------
Lx = 2*pi/(k0*sin(theta0))*nx0;
Hz = pi/(k0*cos(theta0))*nz0;

Nx = 1000;
Nz = 500;

dx0 = Lx/Nx;
dz0 = Hz/Nz;

dx = dx0 * ones(Nx,1);
dz = dz0 * ones(Nz,1);
%--------

%--------
x  = (dx0/2):dx0:(Lx-dx0/2);
z  = (-Hz+dz0/2):dz0:(-dz0/2);

[x,z] = meshgrid(x,z);
%--------

%--------
%n = 10; %nb de bosses
xbot = x(1,:);
%Hbot = 110*(sin(n*2*pi.*xbot./(60000))+1)- Hz;
%Hbot = 110*(sin(2*pi*0.5/(6*10^3).*xbot(1:end))+1)- Hz;
%load hbot.mat
Hbot = 0.*xbot - Hz;
%--------

%--------
%Conditions initiales

Fr = 0.001;
alpha = 2*10^(-4);
g = 9.81;
Psi0 = N*Fr/k0^2;

U0 = cos(theta0)*N*Fr/k0;
V0 = (f/omega0)*cos(theta0)*N*Fr/k0;
W0 = (f/omega0)*sin(theta0)*N*Fr/k0;
T0 = -(N^2/(alpha*g*omega0))*(f/omega0)*sin(theta0)*N*Fr/k0;

Lx = 2*pi/(k0*sin(theta0))*nx0;
Lz = pi/(k0*cos(theta0))*nz0;

deltakx = 2*pi/Lx;
deltakz = pi/Lz;

nx1 = 56;
nz1 = 77;
kx1 = nx1*deltakx;
kz1 = -nz1*deltakz;
k1 = sqrt(kx1^2+kz1^2);
omega1 = sqrt((N*kx1/sqrt(kx1^2+kz1^2))^2+(f*kz1/sqrt(kx1^2+kz1^2))^2);

kx2 = kx0-kx1;
kz2 = kz0-kz1;
k2 = sqrt(kx2^2+kz2^2);
omega2 = -sqrt((N*kx2/sqrt(kx2^2+kz2^2))^2+(f*kz2/sqrt(kx2^2+kz2^2))^2);

epsilon = 0.1;
ropt = 1.3265;
A = (k1/k2)^2*(kz2/kz1)*ropt;
B = (k1/k2)^2*(kz2/kz1)*(omega1/omega2)*ropt;
C = (k1/k2)^2*(kx2/kx1)*(omega1/omega2)*ropt;

U = U0.*cos(kz0.*z).*cos(kx0.*x);%epsilon*U0*(cos(kz1.*z).*cos(kx1.*x)+A*cos(kz2.*z).*cos(kx2.*x));
V = V0.*cos(kz0.*z).*sin(kx0.*x);%epsilon*V0*(cos(kz1.*z).*sin(kx1.*x)+B*cos(kz2.*z).*sin(kx2.*x));
T = T0.*sin(kz0.*z).*cos(kx0.*x);%epsilon*T0*(sin(kz1.*z).*cos(kx1.*x)+C*sin(kz2.*z).*cos(kx2.*x));
TRefrep = (N^2/(alpha*g))*(z+Hz);
R = TRefrep(:,1);

%Paramètres de stabilité onde de base
dt = 30;
Ah = 10^(-2);
Av = Ah;
cgx = (((N^2-f^2)/(omega0*k0))*cos(theta0)*sin(theta0))*sin(theta0);
Slaplacianh = norm(4*Ah*dt/(dx(1)^2)); %<0.3
Slaplacianv = norm(4*Av*dt/(dz(1)^2)); %<0.3
Sfoscillation = norm((f*dt)^2); %<0.5
Scfladvection = norm(U0*dt/dx(1));%<0.5
Sigwpropag = norm(cgx*dt/dx(1));%<0.5

%--------

%--------
U  = flipud(U)';
V  = flipud(V)';
T  = flipud(T)';
R  = flipud(R);
dx = flipud(dx);
dz = flipud(dz);

fid = fopen('U.init','w','b'); 
fwrite(fid,U,'real*8'); 
fclose(fid);

fid = fopen('V.init','w','b'); 
fwrite(fid,V,'real*8'); 
fclose(fid);

fid = fopen('T.init','w','b'); 
fwrite(fid,T,'real*8'); 
fclose(fid);

save('R.init','R','-ascii');

fid = fopen('topog.init','w','b'); 
fwrite(fid,Hbot(1,:)','real*8'); 
fclose(fid);

fid=fopen('delX.init','w','b'); 
fwrite(fid,dx,'real*8'); 
fclose(fid);

fid=fopen('delZ.init','w','b'); 
fwrite(fid,dz,'real*8'); 
fclose(fid);

U  = flipud(U');
V  = flipud(V');
T  = flipud(T');
R  = flipud(R);
dx = flipud(dx);
dz = flipud(dz);
%--------
