<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Modeling and Simulation in Python
Case study.
Copyright 2017 Allen Downey
License
Step1:  Unrolling
Let's simulate a kitten unrolling toilet paper.  As reference material, see this video.
The interactions of the kitten and the paper roll are complex.  To keep things simple, let's assume that the kitten pulls down on the free end of the roll with constant force.  Also, we will neglect the friction between the roll and the axle.  
This figure shows the paper roll with $r$, $F$, and $\tau$.  As a vector quantity, the direction of $\tau$ is into the page, but we only care about its magnitude for now.
We'll start by loading the units we need.
Step2:  And a few more parameters in the Params object.
Step4:  make_system computes rho_h, which we'll need to compute moment of inertia, and k, which we'll use to compute r.
Step5:  Testing make_system
Step7:  Here's how we compute I as a function of r
Step8:  When r is Rmin, I is small.
Step9:  As r increases, so does I.
Step11:  Exercises
Write a slope function we can use to simulate this system.  Here are some suggestions and hints
Step12:  Test slope_func with the initial conditions.
Step13:  Run the simulation.
Step14:  And look at the results.
Step15:  Check the results to see if they seem plausible
Step16:  Plot omega
Step17:  Plot y
<ASSISTANT_TASK:>
Python Code:
# Configure Jupyter so figures appear in the notebook
%matplotlib inline
# Configure Jupyter to display the assigned value after an assignment
%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'
# import functions from the modsim.py module
from modsim import *
Explanation: Modeling and Simulation in Python
Case study.
Copyright 2017 Allen Downey
License: Creative Commons Attribution 4.0 International
End of explanation
radian = UNITS.radian
m = UNITS.meter
s = UNITS.second
kg = UNITS.kilogram
N = UNITS.newton
Explanation: Unrolling
Let's simulate a kitten unrolling toilet paper.  As reference material, see this video.
The interactions of the kitten and the paper roll are complex.  To keep things simple, let's assume that the kitten pulls down on the free end of the roll with constant force.  Also, we will neglect the friction between the roll and the axle.  
This figure shows the paper roll with $r$, $F$, and $\tau$.  As a vector quantity, the direction of $\tau$ is into the page, but we only care about its magnitude for now.
We'll start by loading the units we need.
End of explanation
params = Params(Rmin = 0.02 * m,
                Rmax = 0.055 * m,
                Mcore = 15e-3 * kg,
                Mroll = 215e-3 * kg,
                L = 47 * m,
                tension = 2e-4 * N,
                t_end = 120 * s)
Explanation: And a few more parameters in the Params object.
End of explanation
def make_system(params):
    Make a system object.
    
    params: Params with Rmin, Rmax, Mcore, Mroll,
                              L, tension, and t_end
    
    returns: System with init, k, rho_h, Rmin, Rmax,
                         Mcore, Mroll, ts
    
    L, Rmax, Rmin = params.L, params.Rmax, params.Rmin
    Mroll = params.Mroll
    
    init = State(theta = 0 * radian,
                 omega = 0 * radian/s,
                 y = L)
    
    area = pi * (Rmax**2 - Rmin**2)
    rho_h = Mroll / area
    k = (Rmax**2 - Rmin**2) / 2 / L / radian    
    
    return System(params, init=init, area=area, rho_h=rho_h, k=k)
Explanation: make_system computes rho_h, which we'll need to compute moment of inertia, and k, which we'll use to compute r.
End of explanation
system = make_system(params)
system.init
Explanation: Testing make_system
End of explanation
def moment_of_inertia(r, system):
    Moment of inertia for a roll of toilet paper.
    
    r: current radius of roll in meters
    system: System object with Mcore, rho, Rmin, Rmax
    
    returns: moment of inertia in kg m**2
    
    Mcore, Rmin, rho_h = system.Mcore, system.Rmin, system.rho_h
    
    Icore = Mcore * Rmin**2   
    Iroll = pi * rho_h / 2 * (r**4 - Rmin**4)
    return Icore + Iroll
Explanation: Here's how we compute I as a function of r:
End of explanation
moment_of_inertia(system.Rmin, system)
Explanation: When r is Rmin, I is small.
End of explanation
moment_of_inertia(system.Rmax, system)
Explanation: As r increases, so does I.
End of explanation
# Solution
def slope_func(state, t, system):
    Computes the derivatives of the state variables.
    
    state: State object with theta, omega, y
    t: time
    system: System object with Rmin, k, Mcore, rho_h, tension
    
    returns: sequence of derivatives
    
    theta, omega, y = state
    k, Rmin, tension = system.k, system.Rmin, system.tension
    
    r = sqrt(2*k*y + Rmin**2)
    I = moment_of_inertia(r, system)
    tau = r * tension
    alpha = tau / I
    dydt = -r * omega
    
    return omega, alpha, dydt      
Explanation: Exercises
Write a slope function we can use to simulate this system.  Here are some suggestions and hints:
r is no longer part of the State object.  Instead, we compute r at each time step, based on the current value of y, using
$y = \frac{1}{2k} (r^2 - R_{min}^2)$
Angular velocity, omega, is no longer constant.  Instead, we compute torque, tau, and angular acceleration, alpha, at each time step.
I changed the definition of theta so positive values correspond to clockwise rotation, so dydt = -r * omega; that is, positive values of omega yield decreasing values of y, the amount of paper still on the roll.
Your slope function should return omega, alpha, and dydt, which are the derivatives of theta, omega, and y, respectively.
Because r changes over time, we have to compute moment of inertia, I, at each time step.
That last point might be more of a problem than I have made it seem.  In the same way that $F = m a$ only applies when $m$ is constant, $\tau = I \alpha$ only applies when $I$ is constant.  When $I$ varies, we usually have to use a more general version of Newton's law.  However, I believe that in this example, mass and moment of inertia vary together in a way that makes the simple approach work out.  Not all of my collegues are convinced.
End of explanation
# Solution
slope_func(system.init, 0*s, system)
Explanation: Test slope_func with the initial conditions.
End of explanation
# Solution
results, details = run_ode_solver(system, slope_func)
details
Explanation: Run the simulation.
End of explanation
results.tail()
Explanation: And look at the results.
End of explanation
def plot_theta(results):
    plot(results.theta, color='C0', label='theta')
    decorate(xlabel='Time (s)',
             ylabel='Angle (rad)')
    
plot_theta(results)
Explanation: Check the results to see if they seem plausible:
The final value of theta should be about 220 radians.
The final value of omega should be near 4 radians/second, which is less one revolution per second, so that seems plausible.
The final value of y should be about 35 meters of paper left on the roll, which means the kitten pulls off 12 meters in two minutes.  That doesn't seem impossible, although it is based on a level of consistency and focus that is unlikely in a kitten.
Angular velocity, omega, should increase almost linearly at first, as constant force yields almost constant torque.  Then, as the radius decreases, the lever arm decreases, yielding lower torque, but moment of inertia decreases even more, yielding higher angular acceleration.
Plot theta
End of explanation
def plot_omega(results):
    plot(results.omega, color='C2', label='omega')
    decorate(xlabel='Time (s)',
             ylabel='Angular velocity (rad/s)')
    
plot_omega(results)
Explanation: Plot omega
End of explanation
def plot_y(results):
    plot(results.y, color='C1', label='y')
    decorate(xlabel='Time (s)',
             ylabel='Length (m)')
    
plot_y(results)
Explanation: Plot y
End of explanation
<END_TASK>