function f = funcao_galaxia(x,y)
R0 = 0.04;
R1 = 0.45;
F = 8;
r = sqrt(R0^2 + x^2 + y^2);
s = atan2(y,x);
t = s + F*R0/r;
u = (x^2+y^2)/R1^2;
if u < 1
    b = (r^2)/R0;
    a = b/hypot(1,b)*(1 - u)^2; 
else 
    a = 0;
end
f = a*sin(t);


