#!/usr/bin/python3
# ----06-06-2024  5:37:06----

# Space nuke effects calculator

# kg m s J K

import math
from collections import namedtuple

def MTtoJ(MT):
    return 4.184e15 * MT

# radiant exposure, J/m^2 https://en.wikipedia.org/wiki/Radiant_exposure
# energyMT is in units of MT TNT equivalent
def radiantExposure(energyMT, radius):
    return MTtoJ(energyMT) / (radius**2 *4*math.pi)

# kg/m^3, J/kg, K, J/kg, m, kg
Armor = namedtuple("Armor", "name density energyToVaporize boilingPoint energyToMelt halvingDepth atomMass")
iron = Armor("iron", 7874, 7.687e6, 3143, 937562, 0.0254, 9.27e-26)
k_B = 1.380649e-23 # Boltzmann's constant, J/K
armorList = [iron]

# J/kg at depth x is A 2^(-x/halvingDepth)
# this function returns A (units of J/kg)
def energyDistributionConstant(exposure, armor):
    # integral_0^inf A 2^(-x/halvingDepth) density dx = exposure
    # A*density*halvingDepth/ln(2) = exposure
    return exposure * math.log(2) / (armor.density*armor.halvingDepth)

def vaporizationDepth(exposure, armor, A):
    # A 2^(-x/halvingDepth) = energyToVaporize
    # x = -halvingDepth * ln_2(energyToVaporize/A)
    return armor.halvingDepth * math.log(A/armor.energyToVaporize, 2)

def meltingDepth(exposure, armor, A):
    return armor.halvingDepth * math.log(A/armor.energyToMelt, 2)

# velocity of gas at depth x
def gasVelocity(exposure, armor, x, A):
    # excess energy per atom
    E = (A * 2**(-x/armor.halvingDepth) - armor.energyToVaporize) * armor.atomMass
    return math.sqrt((3*k_B*armor.boilingPoint + 2*E) / armor.atomMass)

def momentumImparted(energyMT, radius, armor, increment=0.0001):
    exposure = radiantExposure(energyMT, radius)
    A = energyDistributionConstant(exposure, armor)
    d = vaporizationDepth(exposure, armor, A)
    tot = 0
    for i in range(int(d/increment) + 1):
        x = increment*i
        plateMass = 1 * 1 * increment * armor.density # 1x1xincrement cuboid (in meters)
        momentum = gasVelocity(exposure, armor, x, A) * plateMass
        tot += momentum
    return tot

# https://www.icrp.org/docs/Hans%20Menzel%20Doses%20From%20Radiation%20Exposure.pdf
# 450 pSv cm2 "per fluence" for 14 MeV neutrons
# "per fluence" means per neutron, which has 2.243e-12 J
# Converting to Sv / (J/m^2)
SievertPerJm2 = 0.02

# humans are basically water, halving thickness of water is 18 cm
# so we might estimate that 50% of the incident radiation is captured by the body
# a 1m x 1m x 18cm plug weighs .18 * (density of water)
# and takes .5*exposure energy
GrayPerJm2 = .5 / (.18 * 997)

# There's a factor of 10 difference between SvPerJm2 and GrayPerJm2
# Which is a problem since Gray and Sievert should be about equally lethal
# SvPerJm2 was calculated from a medical datasheet for 14 MeV neutrons
# GrayPerJm2 was back of the envelope
# on the other hand the radiation is not entirely 14 MeV neutrons

ld50Gray = 5
ld1Gray = 1
ld50Sievert = 5
ld1Sievert = 1

useGray = True

def lethalUnshieldedRadius(energyMT):
    energy = MTtoJ(energyMT)
    # GrayPerJm2 * energy / (radius**2 *4*math.pi) = ld50Gray
    if useGray:
        factor = GrayPerJm2 / ld50Gray
    else:
        factor = SievertPerJm2 / ld50Sievert
    return math.sqrt(energy * factor / (4 * math.pi))

def shieldingNeededForRadiation(energyMT, radius, armor):
    energy = MTtoJ(energyMT)
    exposure = radiantExposure(energyMT, radius)
    if useGray:
        factor = GrayPerJm2 / ld1Gray
    else:
        factor = SievertPerJm2 / ld1Sievert
    # exposure * 2^(-x / halvingDepth) * GrayPerJm2 = ld1Gray
    # x = halvingDepth * log(GrayPerJm2 * exposure / ld1Gray, 2)
    return armor.halvingDepth * math.log(factor * exposure, 2)

# J/kg
TNTequivalent = 4.184e6
TNTdensity = 1654 # kg/m^3

def summary(energyMT, radius, armor, armorThickness, shipDepth, shipDensity):
    exposure = radiantExposure(energyMT, radius)
    energy = MTtoJ(energyMT)
    A = energyDistributionConstant(exposure, armor)
    print("Exposure: {:e} J/m^2".format(exposure))
    print("TNT equivalent: {:.0f} kg TNT per m^2".format(exposure/TNTequivalent))
    print(" (like covering the hull with {:.2f} m of TNT)".format(exposure/(TNTequivalent * TNTdensity)))
    vapDepthCm = vaporizationDepth(exposure, armor, A)*100
    meltDepthCm = meltingDepth(exposure, armor, A)*100
    print("Hull is vaporized to a depth of {:.2f} cm".format(vapDepthCm))
    print(" and melted at depths between {:.2f} and {:.2f} cm".format(vapDepthCm, meltDepthCm))
    momentum = momentumImparted(energyMT, radius, armor)
    print("Each square meter of hull recoils with {:e} kg m/s of momentum".format(momentum))
    print("The {:.2f} m of outer armor recoils at {:.2f} m/s".format(armorThickness, momentum/(armorThickness * armor.density)))
    print("The ship as a whole recoils at {:.2f} m/s".format(momentum/(shipDepth * shipDensity)))
    print("To protect humans from acute irradiation, the {} hull should be at least {:.2f} m thick".format(armor.name, shieldingNeededForRadiation(energyMT, radius, armor)))
    print("An unshielded human would have 50% risk of dying at a radius of {:e} km".format(lethalUnshieldedRadius(energyMT) / 1000))

if __name__ == "__main__":
    try:
        energyMT = float(input("Nuke energy in MT (default 100 MT): "))
    except:
        energyMT = 100
        print("Using default of 100 MT")
    try:
        radius = float(input("Distance from nuke to spaceship in meters (default 1000 m): "))
    except:
        radius = 1000
        print("Using default of 1000 m")
    if len(armorList) > 1:
        armor = "iron"
        try:
            armorChoices = " ".join([a.name for a in armorList])
            armor = input("Armor material (choices = " + armorChoices + "): ")
        except:
            pass
        if not (armor in [a.name for a in armorList]):
            print("Unrecognized armor; using iron")
            armor = "iron"
        armor = [a for a in armorList if a.name == armor][0]
    else:
        armor = iron
    try:
        armorThickness = float(input("Armor thickness in meters (default 1m): "))
    except:
        armorThickness = 1
        print("Using default of 1 m")
    try:
        shipDepth = float(input("Depth of ship from blast side to opposite side in meters (default 100 m):"))
    except:
        shipDepth = 100
        print("Using default of 100 m")
    try:
        shipDensity = float(input("Density of ship in kg/m^3 (default 500 kg/m^3):"))
    except:
        shipDensity = 500
        print("Using default of 500 kg/m^3")

    summary(energyMT, radius, armor, armorThickness, shipDepth, shipDensity)