double lambda_RT = 0.0;
#define L (lambda_RT)  // size of the box
//#define T_END (0.528)
#define T_END (0.4)

#include "../src/properties/film_boiling_case.h"
#define T_inf (5. + T_sat)
#define grav_mag 9.81

#include "../src/01_vaporization/evap_include.h"

T_V[bottom] = dirichlet(T_inf);
T_L[top] = dirichlet(T_sat);

fE_V[bottom] = dirichlet(rhocp_V * T_inf);
fE_L[top] = dirichlet(rhocp_L * T_sat);

u.n[top] = neumann(0.);
p[top] = dirichlet(0.);

u.n[bottom] = dirichlet(0.);
p[bottom] = neumann(0.);

int main(int argc, char* argv[]) {

    lambda_RT = 2.0 * pi * sqrt(3. * sigma_lv / (grav_mag * (rho_L - rho_V)));
    fprintf(stdout, "lambda_RT %g\n", lambda_RT);

    MIN_LEVEL = 5;
    LEVEL = 8;
    MAX_LEVEL = 8;

    setup_evap();

    size(L);
    origin(-0.5*L0, 0.0);
    N = 1 << LEVEL;
    init_grid(N);

    G.y = -grav_mag;

    run();
}

#define RT_y_func(x) ((4. + cos(2. * pi * x / lambda_RT)) * lambda_RT / 128.)
#define RT_func(x, y) (y - RT_y_func(x))

double T_vap_func(double x, double y){
    double y_max = RT_y_func(x);
    double T_y = T_inf - y/y_max*(T_inf-T_sat);
    T_y = clamp(T_y, T_sat, T_inf);
    return T_y;
}

event init(i = 0) {
    if (!restore (file = "restart")) {
        CFL = 0.2;
#if TREE
        refine(level < MAX_LEVEL && RT_func(x, y - dR_refine) < 0. &&
               RT_func(x, y + dR_refine) > 0.);
#endif
        fraction_LS(f, RT_func(x, y));

        foreach () {
            foreach_dimension() {
                u.x[] = 0.;
            }

            double T_vap = T_vap_func(x, y);
            
            T_V[] = (1.0 - f[]) * T_vap + f[] * T_sat;
            T_L[] = T_sat;
        }
        foreach_face() uf.x[] = 0.;
    }

    init_Energy();
}

#define REFINE_VAPOR_REGION
#define femax 1.0e-5
#define Temax 1.0e-2
#define uemax 1.0e-2
#include "../src/01_vaporization/adapt_evap.h"

#include "../src/outputs/progress.h"
event progress(i++) { progress_check(i, t, T_END); }

#include "view.h"
#include "../src/post_processing/movie_maker.h"
event movie_output(t += DELTA_T) {
  scalar temperature[];
  foreach () {
    temperature[] = f[] * T_L[] + (1. - f[]) * T_V[];
  }
  TX = 0.05;
  MESH_ON=0;
  BOX_ON=1;
  movie_maker_i("temperature", i, t);
}