#!/usr/bin/python
# -*- coding: utf-8 -*-
"""What happens when a sledgehammer hits a rock?

Because I’m stupid, I’m pretending that the average force is the force
at half of the deformation, and I’m solving this iteratively instead
of in closed form.

"""
from __future__ import division

if __name__ == '__main__':
    energy_J = 500
    modulus_GPa = 200
    mass_kg = 4
    stop_time_ms = 10
    cross_section_mm_2 = 2800
    head_length_mm = 180

    # E = ½mv²; 2E = mv²; 2E/m = v²; v = √(2E/m)
    speed_m_per_s = (2 * energy_J / mass_kg)**0.5
    print "A hammer massing %.3g kg is going %.3g m/s for an energy of %s J." % (mass_kg, speed_m_per_s, mass_kg * speed_m_per_s**2 / 2)
    print "It’s made from tool steel with a %.3g GPa Young’s modulus." % modulus_GPa
    print "Its head is a cylinder %.3g mm long and %.3g mm² in cross section." % (head_length_mm, cross_section_mm_2)
    print "(That works out to a density of %.3gg/cc.  Steel is 7.8.)" % (mass_kg / (head_length_mm / 100) / (cross_section_mm_2 / 10000))
    print "It hits a rock.  How long does it take to stop?  Let’s guess."

    for ii in range(1024):
        print
        print "So, let’s try with a %.3g ms stop time." % stop_time_ms
        deceleration_m_per_s_2 = 1000 * speed_m_per_s / stop_time_ms
        print "To stop that fast, we must decelerate at an average of %.3g m/s²." % deceleration_m_per_s_2
        motion_during_stop_mm = speed_m_per_s * stop_time_ms / 2
        print "As it decelerates, the hammer deforms by %.3g mm." % motion_during_stop_mm
        average_deformation_mm = motion_during_stop_mm / 2
        average_strain = average_deformation_mm / head_length_mm
        average_stress_GPa = modulus_GPa * average_strain
        print "This gives rise to an average deformation of %.3g mm." % average_deformation_mm
        print "On a %.3g-mm-long hammer head, this amounts to a strain of %.3g%%." % (head_length_mm, average_strain * 100)
        if average_strain > 0.09:
            print "(This is enough to shatter the hammer.)"
        print "At a %.3g GPa Young’s modulus, that creates a stress of %.3g GPa." % (modulus_GPa, average_stress_GPa)
        force_N = 1000 * average_stress_GPa * cross_section_mm_2
        print "Over a %.3g mm² cross section, that creates a force of %.3g N." % (cross_section_mm_2, force_N)
        deceleration_from_force_m_per_s_2 = force_N / mass_kg
        print "That decelerates the hammer at %.3g m/s²." % deceleration_from_force_m_per_s_2
        print "That’s %s deceleration than we need." % ("more" if deceleration_from_force_m_per_s_2 > deceleration_m_per_s_2 else "less")
        ratio = deceleration_from_force_m_per_s_2 / deceleration_m_per_s_2
        print "The ratio is actually %.20g." % ratio
        stop_time_ms = (stop_time_ms + stop_time_ms / ratio) / 2
        if abs(ratio - 1) < 1e-5:
            break

    print
    print "That’d better be good enough."
    print "While stopping in %.3g ms, the hammer dissipates an average of %.3g kW." % (stop_time_ms, energy_J / stop_time_ms)