#!/usr/bin/env bx use b cstr usage[] = "m r0 p [steps=10000]", "mass polar-radius spin-period", "5.9736e24 6356800 86164.09 # Earth", "5.6869e26 54364000 37050.56 # Saturn", NULL num G = 6.67428e-11 # m^3/kg/s Main() num M, r0, p, s, a, b, R2, r2, r, x, y, yp, xp, R, dR, r1, steps steps = 10000 Getargs(num, M, r0, p) getargs(num, steps) warn("M=%g r0=%f p=%f", M, r0, p) s = 2*pi/p paper(600, 400) int max_wh = imax(w, h) zoom((max_wh * 0.6) / (r0*2)) blue() circle(0, 0, r0) black() a = 2*G*M/(s*s) b = 1/r0 x = 0, y = r0, yp = y, xp = 0 dR = r0 / steps for R=r0; R R2 + tinynum: break r2 = clamp(r2, 0, R2) r = sqrt(r2) x = r y = sqrt(R2 - r2) if sd(x) > w_2: break point(x, y) point(-x, y) point(x, -y) point(-x, -y) Paint() if xp && hypot(yp-y, xp-x) > r0 / h / 2: R-=dR ; dR /= 2 ; R+=dR ; again eif xp && hypot(yp-y, xp-x) < r0 / h / 8: R-=dR ; dR *= 2 ; R+=dR ; again yp = y ; xp = x r1 = hypot(x, y) warn("r1 ~= %f", r1) warn("equatorial bulge = 2*(r1-r0) = %f km", 2*(r1-r0)/1000)