#include <math.h>

double swirl(double x, double y, int m, double s)
{
    double d;
    
    if (x == 0.0 && y == 0.0) {
        return 0.0;
    }
    
    d = sqrt(x*x + y*y);

    return cos(m*(atan2(y, x) + s*d))*0.5 + 0.5;
}

double sink(double x, double y, double r)
{
    double d;
    
    d = sqrt(x*x + y*y);
    
    d = d/r - 1.0;
    if (d < 0.0) {
        d = 0.0;
    }

    return 1.0 - 1.0/exp(d);
}

double ripple(double x, double y, double r)
{
    double d;
    
    d = sqrt(x*x + y*y);

    return cos(d/r*M_PI)*0.5 + 0.5;
}

void bulge(double *x, double *y, double cx, double cy, 
    double r, double e)
{
    double vx, vy;
    double d;
    
    vx = *x - cx;
    vy = *y - cy;
    d = sqrt(vx*vx + vy*vy);
    if (d == 0.0) {
        return;
    }
    
    vx /= d;
    vy /= d;
    
    d = pow(d/r, e)*r;
    *x = cx + vx*d;
    *y = cy + vy*d;
}

double fn(double x, double y)
{
    double f;
    double a;
    double dx, dy;
    double d;
    
    y = 1.0 - y;
    
    dx = x;
    dy = y;
    
    bulge(&x, &y, 0.5, 0.7, 0.6, 3.0);
    
    a = 0.1;
    f = ripple(dx - 0.4, dy - 0.4, a)
        + ripple(dx - 0.8, dy - 0.1, a)
        + ripple(dx - 0.9, dy - 0.3, a)
        + ripple(dx - 0.3, dy - 0.8, a)
        + ripple(dx - 0.7, dy - 0.1, a)
        + ripple(dx - 0.4, dy - 0.5, a);
    f /= 6.0;
    
    a *= sink(dx - 0.5, dy - 0.7, 0.4);
    if (a > 0.0) {
        x += cos(f*2.0*M_PI*2.0)*a;
        y += sin(f*2.0*M_PI*2.0)*a;
    }
    
    a = 0.01;
    f = ripple(x - 0.2, y - 0.4, a)
        + ripple(x - 0.4, y - 0.2, a)
        + ripple(x - 0.6, y - 0.8, a)
        + ripple(x - 0.1, y - 0.2, a)
        + ripple(x - 0.9, y - 0.2, a)
        + ripple(x - 0.2, y - 0.3, a);
    f /= 6.0;
    
    x += cos(f*2.0*M_PI*2.0)*a;
    y += sin(f*2.0*M_PI*2.0)*a;
    
    f = pow(swirl(x - 0.5, y - 0.7, 3, 12.0), 0.7);
    f *= sink(x - 0.5, y - 0.7, 0.005);

    return f;
}