Define EPSILON 0.05
Define VAR(T,min,max)   (min + CLAMP(T,0,1)*(max - min))
Define RHO(T)            VAR(T, 1., 1e-3)
Define RADIUS(x,y)      (DIAMETER/2.*(1. + EPSILON*cos (2.*atan2 (y, x))))

1 0 GfsSimulation GfsBox GfsGEdge {} {
    Time { end = 1 }

    Refine LEVEL
 
    VariableTracerVOFHeight T
    VariableFiltered T1 T 1
    InitFraction T ({ x += 0.5; y += 0.5; return x*x + y*y - RADIUS(x,y)*RADIUS(x,y); })

    PhysicalParams { alpha = 1./RHO(T1) }
    VariableCurvature K T
    SourceTension T 1. K

    AdaptFunction { istep = 1 } {
        cmax = 0.01
        maxlevel = LEVEL
    } (T > 0 && T < 1 ? 1. : fabs (Vorticity)*dL)

    OutputScalarSum { istep = 1 } k-LEVEL {
        v = RHO(T1)*Velocity2
    }

    EventScript { start = end } {
        cat <<EOF | gnuplot 2>&1 | awk '{if ($1 == "result:") print LEVEL,$2,$3,$4,$5;}'
           k(t)=a*exp(-b*t)*(1.-cos(c*t))
           a = 3e-4
           b = 1.5
  
           D = DIAMETER
           n = 2.
           sigma = 1.
           rhol = 1.
           rhog = 1./1000.
           r0 = D/2.
           omega0 = sqrt((n**3-n)*sigma/((rhol+rhog)*r0**3))

           c = 2.*omega0
           fit k(x) 'k-LEVEL' u 3:5 via a,b,c
           print "result: ", a, b, c, D
EOF
        rm -f fit.log
    }
}
GfsBox {
    left = Boundary
    bottom = Boundary
}