from CoaSim import *
from CoaSim.popStructure import Population as P, Sample as S, Merge as M
import sys

if len(sys.argv) != 4:
    print 'Usage:', argv[0], 'beta noSeqs merge-time'
    sys.exit(2)

beta = float(sys.argv[1])
noSeqs = int(sys.argv[2])
mergeTime = float(sys.argv[3])

lengths = []
for i in xrange(100000):
    popSpec = P(1, M(mergeTime, [P(1,S(noSeqs)),P(1,S(noSeqs))]))
    arg = simulate([], popSpec, beta=beta, keepEmptyIntervals=True)
    lengths.append(arg.intervals[0].tree.branchLength)
BL = sum(lengths)/len(lengths)
print 'Branch length:', BL
print 'Scaled rho:', 2.0/BL