find_mass.py

#!/usr/bin/env python
import math

m1 = 3.0e30
r1 = 1e12
G = 6.67e-11
P = 100.0 * 364.0 * 24.0 * 60.0 * 60.0


def one(r2):
    return m1 * r1 / r2


def two(r2):
    return math.pow(math.pi, 2) * 4.0 * r1 * math.pow(r1 + r2, 2) / (G * math.pow(P, 2))


def solve():
    r2 = 1e10
    print(f"one:{one(r2)} two:{two(r2)}")
    difference = one(r2) - two(r2)
    prev_diff = difference
    cur_step = 1e19
    prev_move = 0
    cur_move = 0
    while abs(difference) / one(r2) > 0.000000000000000000001:
        if difference > 1:
            print(f"    will increase r2 (diff={difference})")
            prev_move = cur_move
            cur_move = 1
        else:
            print(f"    will decrease r2 (diff={difference})")
            prev_move = cur_move
            cur_move = -1

        if prev_move and cur_move != prev_move:
            print(f"    getting closer, halving the step (step={cur_step})")
            cur_step = cur_step / 2

        r2 = r2 + cur_move * cur_step
        prev_diff = difference
        difference = one(r2) - two(r2)
        print(f"  new diff={difference} with r2={r2}")
        if difference == prev_diff:
            print("We are not advancing, stopping...")
            break

    print(f"Found solution: r2={r2} M2={one(r2)}")


solve()