/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 11.08.0 ] */

/* [wxMaxima: title   start ]
Motion Along A Straight Line, 
Part I: 
Position->Velocity->Acceleration
   [wxMaxima: title   end   ] */

/* [wxMaxima: comment start ]
Chapter 3 Technology Application Project
   [wxMaxima: comment end   ] */

/* [wxMaxima: section start ]
Introduction
   [wxMaxima: section end   ] */

/* [wxMaxima: comment start ]
OBJECTIVE: To apply special Maxima functions to analyze 
position, velocity, and acceleration simultaneously.

Note: There are no executable commands in this Introduction 
section.

Structural engineers must design high-rise buildings to 
withstand earthquakes. When a high-rise building moves back 
and forth during an earthquake, the motion of any single floor 
is also back and forth along a straight, horizontal line. For 
earthquake design, the engineer must be able to completely 
describe the position, velocity, and acceleration of each of 
the floors as functions of time, and understand the 
relationships among these three functions. This provides the 
engineer with some of the information and understanding that 
is needed for the seismic design of the structure in a high-
rise building. (In physics and engineering, the study of 
motion is called kinematics and falls in the realm of 
mechanics.) 

This module includes three specially designed Maxima commands 
that present dramatic animated visualizations of the 
derivative relations among the position, velocity, and 
acceleration. In addition to the seismic vibration of each 
floor in a building, a variety of other motions are 
investigated including constant velocity, constant 
acceleration, harmonic oscillation, and decaying oscillation. 
You can also use the specialized Maxima commands to study 
other motions that are of interest to you. 

To measure position along a line, we select a fixed point on 
the line as a reference point or origin, and we scale the line
in appropriate units (e.g., meters, feet, or miles). A 
distance and a direction specify a position on the line. One 
direction from the reference point is taken as positive and 
the other direction as negative. Distance is measured using 
the units of scale along the line.

Before looking at the vibrations of a building shaken by an 
earthquake, let's consider some other motions. 

This module uses three specially designed functions:

velocity ( fnc, var, interval, periodic, fileprefix ), 
acceleration( fnc, var, time_interval, periodic, fileprefix ),
posvelacc( fnc, var, time_interval, periodic, fileprefix ). 
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
NOTE: The special commands velocity( ), acceleration( ), and 
posvelacc( ) are not built-in Maxima functions and are only 
available in this module. The arguments for all three 
functions are: 

fnc: a list of position functions
var: the variable used in the fnc list functions (usually t)
interval: a list of time intervals of the form [x1,xn] over 
    which the motion occurs (one interval for each function in
    the fnc list)
periodic: a flag to indicate whether the motion is periodic or
    not. 
fileprefix: this is a string that will be placed at the 
    beginning of the filename for the animated GIF file.

In the Parts that follow you will see examples of how these 
functions are used.
   [wxMaxima: comment end   ] */

/* [wxMaxima: subsect start ]
Technology Guidelines
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
NOTE:  If you have just finished a document, restart Maxima 
before executing a new document.  This can be done by choosing
"Restart" from the Maxima menu.

TO OPEN OR CLOSE CELLS
Click on the arrow at the top of the cell bracket. 

TO STOP AN EXECUTION
Click on STOP button from the toolbar. 

ORDER OF EXECUTION
Execute commands in the order given. Do not skip any Maxima
Input lines within a given document. 

Alternatively, you can execute the entire worksheet by 
selecting the "Evaluate All Cells" command from the "Cell" drop
down menu or simply press Ctrl-r. 

SAVING WORKSHEETS
You can save anytime to any directory you choose, and it is 
wise to save often.

EXPERIENCING MAJOR PROBLEMS
Save if appropriate, and then shut down Maxima and start it up 
again. 
   [wxMaxima: comment end   ] */

/* [wxMaxima: subsect start ]
Special Functions
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
The code that follows defines the three special functions that 
are used in this worksheet. You do not have to understand the 
code to do this worksheet, however, you will need to execute 
it in order to define the special functions that are used 
later. You can execute the code by placing the cursor anywhere 
in the input cell and then pressing Enter. After you execute 
the code in the next cell, proceed to Part I of the worksheet.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
load(draw)$
numer:true$
ratprint:false$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
velocity (fnc, var, interval, periodic, fileprefix) := block(
    intlen: length(interval),
    fnclen: length(fnc),
    if intlen#fnclen then block(
        print("The number of functions and number of intervals do not match. Try again."),
        return()),
    t_o: interval[1][1],
    t_f: interval[intlen][2],
    
    d1: diff(fnc, t, 1),
    d2: diff(fnc, t, 2),

    fmax: subst(t_o,var,fnc[1]),
    fmin: subst(t_o,var,fnc[1]),
    d1max: subst(t_o,var,d1[1]),
    d1min: subst(t_o,var,d1[1]),
    d2max: subst(t_o,var,d2[1]),
    d2min: subst(t_o,var,d2[1]),

    k: 1,
    for i: 1 thru 500 do block(
        tcurr: i*(t_f-t_o)/500,
        if tcurr>=interval[k][2] then k: k+1,
        if k<=fnclen then block(
            curr: subst(tcurr,var,fnc[k]),
            if curr>fmax then fmax: curr,
            if curr<fmin then fmin: curr),
        if k<=fnclen then block(
            curr: subst(tcurr,var,d1[k]),          
            if curr>d1max then d1max: curr,
            if curr<d1min then d1min: curr)),
 
    if periodic=1 then tend: t_f-(t_f-t_o)/25 else tend: t_f,

    numframes: 27,
    inc: (t_f-t_o)/(numframes-2),
    marginpercent: 20,
    xmargin: (t_f-t_o)*(marginpercent/100), 
    ymargin: (fmax-fmin)*(marginpercent/100),

    pos_time: [],
    fplot: [],

    xmin: t_o-xmargin,
    xmax: t_f+xmargin,
    ymin: fmin-ymargin,
    ymax: fmax+ymargin,

    for i: 1 thru fnclen do block(
        fplot: append(fplot,[explicit(fnc[i],var,interval[i][1],interval[i][2])])),        

    k: 1,    
    for i: t_o thru tend step inc do block(
        if(k<=fnclen) then block(
           
            tanline:  subst(i,var,d1[k])*(var-i)+subst(i,var,fnc[k]),
            box: rectangle([i, subst(i,var,fnc[k])+subst(i,var,d1[k])], [i+1, subst(i,var,fnc[k])]),
            arrow: vector([i, subst(i,var,fnc[k])],[0,subst(i,var,d1[k])]),
            
            pos_time: append(pos_time, [gr2d(
                nticks=80,
                color=blue,
                line_width=3,
                fplot,
                color=orange,
                line_width=1,
                explicit(tanline,var,xmin,xmax),
                color=black,
                fill_color=white,
                transparent=true,
                box,
                color=red,
                transparent=false,
                head_length=0.5,
                head_angle=15,
                line_width=3,
                arrow,
                xrange=[xmin, xmax],
                yrange=[ymin, ymax],
                title="position vs. time")])),
            
            if ( i+inc >= interval[k][2]) then k: k+1),
      
    draw(
        terminal = animated_gif,
        delay = 1,
        file_name  = sconcat(fileprefix,"v_pos_time"),
        pos_time),

    vel_time: [],
    d1plot: [],
    ymargin: (d1max-d1min)*(marginpercent/100),
    ymin: d1min-ymargin,
    ymax: d1max+ymargin,

    k: 1,
    for i: t_o thru tend step inc do block(
        if(k<=fnclen) then block(
            arrow2: vector([i, 0],[0,subst(i,var,d1[k])]),                   

            for j: 1 thru k do block(
                if interval[j][2]<i then 
                    d1plot: append(d1plot,[explicit(d1[j],var,interval[j][1],interval[j][2])]),
                if (interval[j][1]<i and interval[j][2]>i) then block(
                    d1plot: append(d1plot,[explicit(d1[j],var,interval[j][1],i)]))),     
   
            vel_time: append(vel_time, [gr2d(
                nticks=80,
                color=red,
                line_width=1,
                line_type=dots,
                d1plot,
                color=red,
                head_length=0.5,
                head_angle=15,
                line_width=3,
                line_type=solid,
                arrow2,
                xrange=[xmin, xmax],
                yrange=[ymin, ymax],
                title="velocity vs. time")])),
            
            if ( i+inc >= interval[k][2]) then k: k+1),
      
    draw(
        terminal = animated_gif,
        delay = 1,
        file_name  = sconcat(fileprefix,"v_vel_time"),
        vel_time),

    objmot: [],
    ymargin: (fmax-fmin)*(marginpercent/100),
    ymin: fmin-ymargin,
    ymax: fmax+ymargin,

    k: 1,    
    for i: t_o thru tend step inc do block(
        if(k<=fnclen) then block(      
     
            point1: points([[(t_f-t_o)/2, subst(i,var,fnc[k])]]),
            arrow: vector([(t_f-t_o)/2, subst(i,var,fnc[k])],[0,subst(i,var,d1[k])]),
            
            objmot: append(objmot, [gr2d(
                nticks=80,
                color=red,
                head_length=0.5,
                head_angle=15,
                line_width=3,
                arrow,
                color=black,
                point_type=circle,
                point1,
                xtics=false,
                xrange=[xmin, xmax],
                yrange=[ymin, ymax],
                title="objects motion")])),
            
        if ( i+inc >= interval[k][2]) then k: k+1),
      
    draw(
        terminal = animated_gif,
        delay = 1,
        file_name  = sconcat(fileprefix,"v_obj_mot"),
        objmot))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
acceleration (fnc, var, interval, periodic, fileprefix) := block(
    intlen: length(interval),
    fnclen: length(fnc),
    if intlen#fnclen then block(
        print("The number of functions and number of intervals do not match. Try again."),
        return()),
    t_o: interval[1][1],
    t_f: interval[intlen][2],
    
    d1: diff(fnc, t, 1),
    d2: diff(fnc, t, 2),

    fmax: subst(t_o,var,fnc[1]),
    fmin: subst(t_o,var,fnc[1]),
    d1max: subst(t_o,var,d1[1]),
    d1min: subst(t_o,var,d1[1]),
    d2max: subst(t_o,var,d2[1]),
    d2min: subst(t_o,var,d2[1]),

    k: 1,
    for i: 1 thru 500 do block(
        tcurr: i*(t_f-t_o)/500,
        if tcurr>=interval[k][2] then k: k+1,
        if k<=fnclen then block(
            curr: subst(tcurr,var,fnc[k]),
            if curr>fmax then fmax: curr,
            if curr<fmin then fmin: curr),
        if k<=fnclen then block(
            curr: subst(tcurr,var,d1[k]),          
            if curr>d1max then d1max: curr,
            if curr<d1min then d1min: curr),
        if k<=fnclen then block(
            curr: subst(tcurr,var,d2[k]),
            if curr>d2max then d2max: curr,
            if curr<d2min then d2min: curr)),
 
    if periodic=1 then tend: t_f-(t_f-t_o)/25 else tend: t_f,

    numframes: 27,
    inc: (t_f-t_o)/(numframes-2),
    marginpercent: 20,
    xmargin: (t_f-t_o)*(marginpercent/100), 
    ymargin: (d1max-d1min)*(marginpercent/100),

    vel_time: [],
    d1plot: [],

    xmin: t_o-xmargin,
    xmax: t_f+xmargin,
    ymin: d1min-ymargin,
    ymax: d1max+ymargin,

    for i: 1 thru fnclen do block(
        d1plot: append(d1plot,[explicit(d1[i],var,interval[i][1],interval[i][2])])),        

    k: 1,    
    for i: t_o thru tend step inc do block(
        if(k<=fnclen) then block(
           
            tanline:  subst(i,var,d2[k])*(var-i)+subst(i,var,d1[k]),
            box: rectangle([i, subst(i,var,d1[k])+subst(i,var,d2[k])], [i+1, subst(i,var,d1[k])]),
            arrow: vector([i, subst(i,var,d1[k])],[0,subst(i,var,d2[k])]),
            
            vel_time: append(vel_time, [gr2d(
                nticks=80,
                color=red,
                line_width=3,
                d1plot,
                color=orange,
                line_width=1,
                explicit(tanline,var,xmin,xmax),
                color=black,
                fill_color=white,
                transparent=true,
                box,
                color=green,
                transparent=false,
                head_length=0.5,
                head_angle=15,
                line_width=3,
                arrow,
                xrange=[xmin, xmax],
                yrange=[ymin, ymax],
                title="velocity vs. time")])),
            
            if ( i+inc >= interval[k][2]) then k: k+1),
      
    draw(
        terminal = animated_gif,
        delay = 1,
        file_name  = sconcat(fileprefix,"a_vel_time"),
        vel_time),

    acc_time: [],
    d2plot: [],
    ymargin: (d2max-d2min)*(marginpercent/100),
    ymin: d2min-ymargin,
    ymax: d2max+ymargin,

    k: 1,
    for i: t_o thru tend step inc do block(
        if(k<=fnclen) then block(
            arrow2: vector([i, 0],[0,subst(i,var,d2[k])]),                   

            for j: 1 thru k do block(
                if interval[j][2]<i then 
                    d2plot: append(d2plot,[explicit(d2[j],var,interval[j][1],interval[j][2])]),
                if (interval[j][1]<i and interval[j][2]>i) then block(
                    d2plot: append(d2plot,[explicit(d2[j],var,interval[j][1],i)]))),     
   
            acc_time: append(acc_time, [gr2d(
                nticks=80,
                color=green,
                line_width=1,
                line_type=dots,
                d2plot,
                color=green,
                head_length=0.5,
                head_angle=15,
                line_width=3,
                line_type=solid,
                arrow2,
                xrange=[xmin, xmax],
                yrange=[ymin, ymax],
                title="acceleration vs. time")])),
            
            if ( i+inc >= interval[k][2]) then k: k+1),
      
    draw(
        terminal = animated_gif,
        delay = 1,
        file_name  = sconcat(fileprefix,"a_acc_time"),
        acc_time),

    objmot: [],
    ymargin: (fmax-fmin)*(marginpercent/100),
    ymin: fmin-ymargin,
    ymax: fmax+ymargin,

    k: 1,    
    for i: t_o thru tend step inc do block(
        if(k<=fnclen) then block(      
     
            point1: points([[(t_f-t_o)/2, subst(i,var,fnc[k])]]),
            arrow: vector([(t_f-t_o)/2, subst(i,var,fnc[k])],[0,subst(i,var,d2[k])]),
            
            objmot: append(objmot, [gr2d(
                nticks=80,
                color=green,
                head_length=0.5,
                head_angle=15,
                line_width=3,
                arrow,
                color=black,
                point_type=circle,
                point1,
                xtics=false,
                xrange=[xmin, xmax],
                yrange=[ymin, ymax],
                title="objects motion")])),
            
        if ( i+inc >= interval[k][2]) then k: k+1),
      
    draw(
        terminal = animated_gif,
        delay = 1,
        file_name  = sconcat(fileprefix,"a_obj_mot"),
        objmot))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
posvelacc (fnc, var, interval, periodic, fileprefix) := block(
    intlen: length(interval),
    fnclen: length(fnc),
    if intlen#fnclen then block(
        print("The number of functions and number of intervals do not match. Try again."),
        return()),
    t_o: interval[1][1],
    t_f: interval[intlen][2],
    
    d1: diff(fnc, t, 1),
    d2: diff(fnc, t, 2),

    k: 1,
    for i: 1 thru 500 do block(
        tcurr: i*(t_f-t_o)/500,
        if tcurr>=interval[k][2] then k: k+1,
        if k<=fnclen then block(
            curr: subst(tcurr,var,fnc[k]),
            if curr>fmax then fmax: curr,
            if curr<fmin then fmin: curr),
        if k<=fnclen then block(
            curr: subst(tcurr,var,d1[k]),          
            if curr>d1max then d1max: curr,
            if curr<d1min then d1min: curr),
        if k<=fnclen then block(
            curr: subst(tcurr,var,d2[k]),
            if curr>d2max then d2max: curr,
            if curr<d2min then d2min: curr)),
 
    if periodic=1 then tend: t_f-(t_f-t_o)/25 else tend: t_f,

    numframes: 27,
    inc: (t_f-t_o)/(numframes-2),
    marginpercent: 20,
    xmargin: (t_f-t_o)*(marginpercent/100), 
    ymargin: (fmax-fmin)*(marginpercent/100),

    pva: [],

    xmin: t_o-xmargin,
    xmax: t_f+xmargin,
    ymin: fmin-ymargin,
    ymax: fmax+ymargin,
    t_pos: (t_f-t_o)/2,
    t_vel: t_pos+(t_f-t_o)/4,
    t_acc: t_pos-(t_f-t_o)/4,

    k: 1,    
    for i: t_o thru tend step inc do block(
        if(k<=fnclen) then block(
           
            farrow: vector([t_pos, 0], [0, subst(i,var,fnc[k])]),
            d1arrow: vector([t_vel, subst(i,var,fnc[k])],[0,subst(i,var,d1[k])]),
            d2arrow: vector([t_acc, subst(i,var,fnc[k])],[0,subst(i,var,d2[k])]),
            
            pva: append(pva, [gr2d(
                nticks=80,
                color=blue,
                head_length=0.5,
                head_angle=15,
                line_width=3,
                farrow,
                color=red,
                d1arrow,
                color=green,
                d2arrow,
                xtics=false,
                xrange=[xmin, xmax],
                yrange=[ymin, ymax],
                title="s(blue), v(red), a(green) vectors")])),
            
            if ( i+inc >= interval[k][2]) then k: k+1),
      
    draw(
        terminal = animated_gif,
        delay = 1,
        file_name  = sconcat(fileprefix,"pva"),
        pva))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: section start ]
Part I: Constant Velocity
   [wxMaxima: section end   ] */

/* [wxMaxima: comment start ]
Section 3.3

Suppose that you drive your car along a straight road at a 
constant speed of 60 mph for 4 hours and then turn around and 
come back at the same speed. Let's take the reference point to
be your starting position, and the direction your car travels 
during the first four hours as the positive direction. The 
position of your car during the trip is given by the piecewise
function s(t) = 60*t for t<=4 and 480-60*t for t>4.  Note that 
s is typically used to represent the position function.

Now we plot this piecewise function over the duration of the 
trip.  
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
s: [60*t, 480-60*t];
sinterval: [[0,4],[4,8]];
slen: length(s)$
t_o: sinterval[1][1]$
t_f: sinterval[slen][2]$
splot: []$

smax: subst(t_o,t,s[1])$
smin: subst(t_o,t,s[1])$

k: 1$
for i: 1 thru 500 do block(
    tcurr: i*(t_f-t_o)/500,
    if (tcurr>=sinterval[k][2] and k+1<=slen) then k: k+1,
    if k<=slen then block(
        curr: subst(tcurr,t,s[k]),
        if curr>smax then smax: curr,
        if curr<smin then smin: curr))$

marginpercent: 10$
xmargin: (t_f-t_o)*(marginpercent/100)$
ymargin: (smax-smin)*(marginpercent/100)$

for i: 1 thru slen do block(
        splot: append(splot,[explicit(s[i],t,sinterval[i][1],sinterval[i][2])]))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=80,
    color=blue,
    splot,
    xrange=[t_o-xmargin, t_f+xmargin],
    yrange=[smin-ymargin, smax+ymargin],
    xlabel="t",
    ylabel="s")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
The velocity is defined to be the rate of change of 
position over time. That is v = ds/dt. Let's find the velocity
function and plot it.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
v: diff(s,t,1);
vinterval: [[0,4],[4,8]];
vlen: length(v)$
vplot: []$

vmax: subst(t_o,t,v[1])$
vmin: subst(t_o,t,v[1])$

k: 1$
for i: 1 thru 500 do block(
    tcurr: i*(t_f-t_o)/500,
    if (tcurr>=vinterval[k][2] and k+1<=vlen) then k: k+1,
    if k<=vlen then block(
        curr: subst(tcurr,t,v[k]),
        if curr>vmax then vmax: curr,
        if curr<vmin then vmin: curr))$

marginpercent: 10$
xmargin: (t_f-t_o)*(marginpercent/100)$
ymargin: (vmax-vmin)*(marginpercent/100)$

for i: 1 thru vlen do block(
        vplot: append(vplot,[explicit(v[i],t,vinterval[i][1],vinterval[i][2])]))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=80,
    color=red,
    vplot,
    xrange=[t_o-xmargin, t_f+xmargin],
    yrange=[vmin-ymargin, vmax+ymargin],
    xlabel="t",
    ylabel="s")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
From the previous graph, we see that the velocity, which is 
the derivative of the position function, is undefined at t = 4
because the two one-sided derivatives of s at t = 4 do not 
match. 

Like position, velocity has a size or magnitude and a 
direction. The magnitude or absolute value of an object's 
velocity is called its speed. If the position function, s(t), 
is increasing with time, then the velocity is positive, and if
it is decreasing then the velocity is negative. Physical 
quantities that have a magnitude and a direction are called 
vectors, and the position and velocity of a moving object are 
vector quantities.

Recall that the derivative of a function at time t is the 
slope of the line tangent to the graph of the function at that
point. Therefore, the velocity at time t is the slope of the 
position function at that time. The velocity( ) command in the
next executable section below illustrates this idea by 
producing a sequence of graphs that show:  

1) the position function, s(t)  
2) the velocity function, v(t)
3) the object as it moves along a straight line. 

The velocity vector is shown on all three graphs. In the 
position-versus-time graph, the width of the small black 
rectangle is always one so that the length of the arrow is the
slope of the orange tangent line. Note that the value for the 
last argument is zero since the motion is not periodic. In the
first two animation graphs, the horizontal axis represents 
time.

After the following function is complete, the animated graphs
can be found in the wxMaxima folder on your computer.  They
will be named 301v_pos_time.gif, 301v_vel_time.gif, and
301v_obj_mot.gif.

The next command may take a little while to execute, depending
on the memory size and speed of your computer.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
velocity ([60*t,480-60*t], t, [[0,4],[4,8]], 0, "301")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
A note about modeling: In the preceding model, the velocity is
undefined at t = 4 hours. In reality, however, your car would 
have to slow down, turn around, and accelerate back to 
cruising speed at t = 4 hours. This maneuver would occur over 
a very short interval of time (say a few minutes) when 
compared to the 4-hour duration of the trip. The sharp corner 
at the turnaround time would in actuality be rounded, and the 
sudden jump from 60 mph to a mph would have a finite negative 
slope to it. On the scale of eight hours, however, the graphs 
would look much like those depicted above, even if the more 
precise model where used. With these observations in mind, the
above model is reasonable for describing your trip.
   [wxMaxima: comment end   ] */

/* [wxMaxima: subsect start ]
You Try It: Constant Rate of Change and Extreme
Values
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Over each 4-hour subinterval of the motion in Part I, the 
velocity of your car is constant at ± 60 mph. Each subinterval
provides an example of a situation where the rate of change or
the derivative of a function is constant. Write a brief 
response to each of the following questions.

1.  On the interval where v = ds/dt is positive, what can you 
say about the value of s? 

2.  What can you say about the value of s on the interval 
where v = ds/dt is negative?

3.  When v = ds/dt is constant on an interval, what can you 
say about the relationship between v = ds/dt (the slope of the
tangents in the interval) and delta(s)/delta(t) (the slope of 
any secant line taken inside the interval)?

4.  At a hours, the function s reaches its maximum value of 
240 miles. What happens to the derivative, v = ds/dt, at this 
point? What are the values of the two one-sided derivatives of
s at t = 4? What are the values of v = ds/dt slightly to the 
left of t = 4 and slightly to the right of t = 4?

5.  Over the 8-hour interval, what are the largest and 
smallest values of s and at what times to they occur?

6.  What is the value of v = (d^2 s)/(dt^2), the acceleration,
during each of the 4-hour subintervals? In the mathematical 
model, what happens to the acceleration at t = 4 hours? 
Describe what the acceleration of your car would actually be 
like when you turn around to come back. 
   [wxMaxima: comment end   ] */

/* [wxMaxima: section start ]
Part II: Constant Acceleration
   [wxMaxima: section end   ] */

/* [wxMaxima: comment start ]
Section 3.3

If you throw an object straight up from the ground with an 
initial speed of 128 ft/s, its height above the ground is 
given by the following function of time.

NOTE: The commands in this worksheet generate a lot of 
graphics, which may fill up your computer's memory. Before 
proceeding with this Part, you should pull down the Cell menu 
at the top of the screen and select "Remove All Output" (From 
Worksheet). Then you will have to return to the place in the 
worksheet where you left off.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
s: [128*t-16*t^2];
sinterval: [[0,8]];
slen: length(s)$
t_o: sinterval[1][1]$
t_f: sinterval[slen][2]$
splot: []$

smax: subst(t_o,t,s[1])$
smin: subst(t_o,t,s[1])$

k: 1$
for i: 1 thru 500 do block(
    tcurr: i*(t_f-t_o)/500,
    if (tcurr>=sinterval[k][2] and k+1<=slen) then k: k+1,
    if k<=slen then block(
        curr: subst(tcurr,t,s[k]),
        if curr>smax then smax: curr,
        if curr<smin then smin: curr))$

marginpercent: 10$
xmargin: (t_f-t_o)*(marginpercent/100)$
ymargin: (smax-smin)*(marginpercent/100)$

for i: 1 thru slen do block(
        splot: append(splot,[explicit(s[i],t,sinterval[i][1],sinterval[i][2])]))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=80,
    color=blue,
    splot,
    xrange=[t_o-xmargin, t_f+xmargin],
    yrange=[smin-ymargin, smax+ymargin],
    xlabel="t",
    ylabel="s");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Note that the origin is taken at the ground surface and the 
positive direction is upward.

Now determine the velocity function and graph it.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
v: diff(s,t,1);
vinterval: [[0,8]];
vlen: length(v)$
vplot: []$

vmax: subst(t_o,t,v[1])$
vmin: subst(t_o,t,v[1])$

k: 1$
for i: 1 thru 500 do block(
    tcurr: i*(t_f-t_o)/500,
    if (tcurr>=vinterval[k][2] and k+1<=vlen) then k: k+1,
    if k<=vlen then block(
        curr: subst(tcurr,t,v[k]),
        if curr>vmax then vmax: curr,
        if curr<vmin then vmin: curr))$

marginpercent: 10$
xmargin: (t_f-t_o)*(marginpercent/100)$
ymargin: (vmax-vmin)*(marginpercent/100)$

for i: 1 thru vlen do block(
        vplot: append(vplot,[explicit(v[i],t,vinterval[i][1],vinterval[i][2])]))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=80,
    color=red,
    vplot,
    xrange=[t_o-xmargin, t_f+xmargin],
    yrange=[vmin-ymargin, vmax+ymargin],
    xlabel="t",
    ylabel="s");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Let's use velocity( ) to depict the motion. Note that the 
value for the last argument is 0 since the motion is not 
periodic.

After the following function is complete, the animated graphs
can be found in the wxMaxima folder on your computer.  They
will be named 302v_pos_time.gif, 302v_vel_time.gif, and
302v_obj_mot.gif.

The next command may take a little while to execute, depending
on the memory size and speed of your computer.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
velocity ([128*t-16*t^2], t, [[0,8]], 0, "302")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
In this example, the rate of change or the slope of the 
velocity function (i.e., the acceleration) is constant at 
-32 ft/s^2.

Acceleration is also a vector quantity. If the acceleration is
in the negative direction, the velocity (which can be positive
or negative) is decreasing with time.  And, if the acceleration
is in the positive direction, the velocity (which can be 
positive or negative) is increasing with time. If the velocity
vector and the acceleration vector are in the same direction, 
then the moving object is speeding up, and if they are in 
opposite directions, then the object is slowing down. 

The acceleration( ) command produces a sequence of graphs that
show: 

1) the velocity function, v(t) 
2) the acceleration function, a(t)
3) the object as it moves along a straight line 

The acceleration vector is shown on all three graphs. The 
width of the small black rectangle on the v(t) graph is always
one so that the length of the arrow is the slope of the orange
tangent line.

After the following function is complete, the animated graphs
can be found in the wxMaxima folder on your computer.  They
will be named 302a_vel_time.gif, 302a_acc_time.gif, and
302a_obj_mot.gif.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
acceleration ([128*t-16*t^2], t, [[0,8]], 0, "302")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
The posvelacc( ) animation shows the motion of the object, 
together with its (blue) position, (red) velocity, and 
(green) acceleration vectors.

After the following function is complete, the animated graph
can be found in the wxMaxima folder on your computer.  It
will be named 302pva.gif.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
posvelacc ([128*t-16*t^2], t, [[0,8]], 0, "302")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Note that on the way up, the velocity and acceleration are in 
opposite directions (the velocity is positive and the 
acceleration is negative) and the object is slowing down, 
whereas, on the way down the velocity and acceleration are in 
the same direction (both negative) and the object is speeding 
up.

Before we continue, there is one last thing we need to point 
out pertaining to the physics and mathematics of motion. While
the purpose of this module is to illustrate the derivative 
relations between position and velocity and between velocity 
and acceleration, real world applications usually work the 
other way around. That is, instead of starting with the 
position function and differentiating it to find the velocity 
and acceleration, we usually start with the acceleration 
function and use it to construct the velocity and position 
functions. The reason for this is that Isaac Newton taught us,
among other things, that the best way to study real world 
phenomena is to look at the forces that bring about change. 
In the context of motion problems, this philosophy is embodied
in Newton's three laws of motion. The second of Newton's three
laws states that if the mass of an object is constant then the
total of all the forces acting on the object is equal to its 
mass times its acceleration, that is, F = m a. This gives us a
way to determine the acceleration function for an object by 
knowing the forces that act on it. Once we know the 
acceleration we use it to construct the velocity function, and
then we us the velocity to construct the position function. 
The problem then is this: if we know the derivative of a 
function, how do we determine what the function is? This 
problem is dealt with in the second big part of calculus - 
antidifferentiation and integration. You can explore these 
ideas in a later module (Motion Along a Straight Line, 
Part II: Acceleration -> Velocity -> Position, which treats 
some of the same motion problems that are in this module, but 
the other way around.
   [wxMaxima: comment end   ] */

/* [wxMaxima: subsect start ]
You Try It - Linear Rate of Change and Extreme Values
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Section 3.3, 4.4

The motion in Part II provides an example of a situation where
the derivative of a function, i.e., the velocity, is a linear
function. Write a brief response to each of the following 
questions.

1.  On what interval is v = ds/dt positive? On what interval 
is it negative?

2.  On the interval where v = ds/dt is positive, what can you 
say about the value of s? 

3.  What can you say about the value of s on the interval 
where v = ds/dt is negative?

4.  At t=4 seconds, the function s reaches its maximum value 
of 256 feet. What happens to the derivative, v = ds/dt, at 
this point? 

5.  What does the value of v = ds/dt do as t increases? How is
this related to the value of a = dv/dt = (d^2 s)/dt^2 and the 
concavity of the graph of s over the interval of motion?

6.  Over the 8-second interval, what are the largest and 
smallest values of s and at what times to they occur?
   [wxMaxima: comment end   ] */

/* [wxMaxima: section start ]
Part III: Oscillations
   [wxMaxima: section end   ] */

/* [wxMaxima: comment start ]
Section 3.4, Simple Harmonic Motion

If you hang a mass from a spring, it will eventually come to 
rest in an equilibrium position, where the weight of the mass 
is balanced by the tension in the spring. Now pull the mass 
down slightly and let it go. It will oscillate (move up and 
down) about the equilibrium position. If there were no air 
resistance and/or friction in the mass and spring it would 
oscillate forever with the same amplitude and frequency. Since
the motion is periodic, it is not surprising that the periodic
trig functions can be used to describe such a motion. For 
example, if the amplitude of the vibration is 2 centimeters 
and the period of the oscillations is 4*%pi/3 seconds, then 
the motion is described by the following function.

NOTE: The commands in this worksheet generate a lot of 
graphics, which may fill up your computer's memory. Before 
proceeding with this Part, you should pull down the Cell menu 
at the top of the screen and select "Remove All Output" (From 
Worksheet). Then you will have to return to the place in the 
worksheet where you left off.

   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
s: [-2*cos(3*t/2)];
sinterval: [[0,8*%pi/3]];
slen: length(s)$
t_o: sinterval[1][1]$
t_f: sinterval[slen][2]$
splot: []$

smax: subst(t_o,t,s[1])$
smin: subst(t_o,t,s[1])$

k: 1$
for i: 1 thru 500 do block(
    tcurr: i*(t_f-t_o)/500,
    if (tcurr>=sinterval[k][2] and k+1<=slen) then k: k+1,
    if k<=slen then block(
        curr: subst(tcurr,t,s[k]),
        if curr>smax then smax: curr,
        if curr<smin then smin: curr))$

marginpercent: 10$
xmargin: (t_f-t_o)*(marginpercent/100)$
ymargin: (smax-smin)*(marginpercent/100)$

for i: 1 thru slen do block(
        splot: append(splot,[explicit(s[i],t,sinterval[i][1],sinterval[i][2])]))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=80,
    color=blue,
    splot,
    xrange=[t_o-xmargin, t_f+xmargin],
    yrange=[smin-ymargin, smax+ymargin],
    xlabel="t",
    ylabel="s");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Now we can calculate the velocity of the mass as it moves up 
and down and graph it.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
v: diff(s,t,1);
vinterval: [[0,8*%pi/3]];
vlen: length(v)$
vplot: []$

vmax: subst(t_o,t,v[1])$
vmin: subst(t_o,t,v[1])$

k: 1$
for i: 1 thru 500 do block(
    tcurr: i*(t_f-t_o)/500,
    if (tcurr>=vinterval[k][2] and k+1<=vlen) then k: k+1,
    if k<=vlen then block(
        curr: subst(tcurr,t,v[k]),
        if curr>vmax then vmax: curr,
        if curr<vmin then vmin: curr))$

marginpercent: 10$
xmargin: (t_f-t_o)*(marginpercent/100)$
ymargin: (vmax-vmin)*(marginpercent/100)$

for i: 1 thru vlen do block(
        vplot: append(vplot,[explicit(v[i],t,vinterval[i][1],vinterval[i][2])]))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=80,
    color=red,
    vplot,
    xrange=[t_o-xmargin, t_f+xmargin],
    yrange=[vmin-ymargin, vmax+ymargin],
    xlabel="t",
    ylabel="s");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
The velocity( ) animation can be used to depict the motion. 
Since the motion is periodic, we set the periodic argument
of velocity( ) to 1 instead of 0. 

The animated gif filenames created will be 303v_pos_time.gif, 
303v_vel_time.gif, and 303v_obj_mot.gif.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
velocity ([2*cos(3*t/2)], t, [[0,8*%pi/3]], 1, "303")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
The acceleration of the mass is given by the derivative of 
the velocity.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
a: diff(v,t,1);
ainterval: [[0,8*%pi/3]];
alen: length(a)$
aplot: []$

amax: subst(t_o,t,a[1])$
amin: subst(t_o,t,a[1])$

k: 1$
for i: 1 thru 500 do block(
    tcurr: i*(t_f-t_o)/500,
    if (tcurr>=ainterval[k][2] and k+1<=alen) then k: k+1,
    if k<=alen then block(
        curr: subst(tcurr,t,a[k]),
        if curr>amax then amax: curr,
        if curr<amin then amin: curr))$

marginpercent: 10$
xmargin: (t_f-t_o)*(marginpercent/100)$
ymargin: (amax-amin)*(marginpercent/100)$

for i: 1 thru vlen do block(
        aplot: append(aplot,[explicit(a[i],t,ainterval[i][1],ainterval[i][2])]))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=80,
    color=green,
    aplot,
    xrange=[t_o-xmargin, t_f+xmargin],
    yrange=[amin-ymargin, amax+ymargin],
    xlabel="t",
    ylabel="s");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
The acceleration( ) animation illustrates the relationship 
between the velocity and acceleration as the mass oscillates.

The animated gif filenames created will be 303a_vel_time.gif, 
303a_acc_time.gif, and 303a_obj_mot.gif.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
acceleration ([-2*cos(3*t/2)], t, [[0,8*%pi/3]], 1, "303")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
The posvelacc( ) animation shows the motion of the object, 
position vector, its velocity, and its acceleration.

The animated gif filename created will be 303pva.gif
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
posvelacc ([-2*cos(3*t/2)], t, [[0,8*%pi/3]], 1, "303")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
There are some of observations to make about this motion. The 
first is that the acceleration is in the direction opposite 
the position. When the mass is below the equilibrium position,
the stretch in the spring increases thus exerting a net upward
force on the mass, accelerating in that direction. And, when 
the mass is above the equilibrium position, the stretch in the
spring is reduced resulting in a net downward force on the 
mass, accelerating it in that direction. The second 
observation is that the acceleration is proportional to the 
position. The first two observations can be demonstrated 
mathematically by noting that if s(t) = -A*cos(w*t) then 
a(t) = (d^2 s)/dt^2 = A*w^2*cos(w*t) = w^2*A*cos(w*t) = 
-w^2*s(t).

Another thing to note is that the velocity and acceleration 
are 90 degrees or %pi/2 radians out of phase with one another.
During the first quarter-cycle of the motion, the velocity and
acceleration are both positive and the mass is speeding up, 
moving upward. During the second-quarter cycle, the velocity 
is positive and the acceleration is negative, and the mass is 
still moving upward but now it is slowing down. During the 
third quarter-cycle, the velocity and acceleration are both 
negative, and the mass is speeding up moving downward. During 
the fourth quarter-cycle, the velocity is negative and the 
acceleration is positive, and the mass is slowing down as it 
moves downward. Overlaying the graphs of the velocity and 
acceleration functions shows these relationships.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=80,
    color=blue,
    splot,
    color=red,
    vplot,
    color=green,
    aplot,
    xrange=[t_o-xmargin, t_f+xmargin],
    yrange=[amin-ymargin, amax+ymargin],
    xlabel="t",
    ylabel="s");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
The velocity is v(t) = -A*w*sin(w*t) and the acceleration is 
a(t) = -A*w^2*cos(w*t). From basic trigonometry we know that 
the sine and cosine functions are 90 degrees or %pi/2 radians 
out of phase with one another.

The final observation is that at the instant when the mass 
passes the equilibrium position, the velocity reaches its 
extreme values (positive on the way up and negative on the way
down), its speed is at its maximum value, and its acceleration
is zero.
   [wxMaxima: comment end   ] */

/* [wxMaxima: subsect start ]
You Try It - The Second Derivative, Concavity, and 
Inflection Points
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Section 3.3, 4.4

After you have viewed the velocity( ), acceleration( ), and 
posvelacc( ) animations, write a brief response to each of the
following questions.

1. During intervals where the tangent line is rotating 
clockwise the graph of s is said to be "concave down". Over 
what intervals is the graph of s concave down? What is the 
value of v = ds/dt doing in these intervals? What is the sign 
of a = (d^2 s)/dt^2 in these intervals?

2. During the intervals where the tangent line is rotating 
counterclockwise, the graph of s is said to be "concave up". 
Over what intervals is the graph of s concave up? What is the 
value of v = ds/dt doing in these intervals? What is the sign 
of a = (d^2 s)/dt^2 in these intervals?

3. An inflection point is a point on the graph of s where the 
concavity changes from up to down, or visa versa. Describe the
motion of the green tangent line as it passes over the 
inflection points on the graph of s. What happens to the value
of v = ds/dt as the point on the graph passes over the 
inflection points? What happens to the value of 
a = dv/dt = (d^2 s)/dt^2 as the point on the graph passes 
over the inflection points?

4. At what points on the graph of s is the slope the steepest?
At what points on the graph of s does the value of v = ds/dt
reach its extreme (maximum and minimum) values?

5. Describe the motion of an object that moves up and down 
along a straight line for each of the following conditions: 
a) v > 0 and a > 0
b) v > 0 and a < 0
c) v < 0 and a < 0 
d) v < 0 and a > 0
With up taken as the positive direction, specify the direction
the object is moving and whether it is speeding up or slowing 
down.
   [wxMaxima: comment end   ] */

/* [wxMaxima: section start ]
Part IV: Decaying Oscillations
   [wxMaxima: section end   ] */

/* [wxMaxima: comment start ]
Chapter 3, Section 4

If you wish, you can quit here. The next two Parts, IV and V, 
are just for fun. 

A more realistic model for oscillations of real objects is to 
account for the loss of energy, which results in a decay in 
the amplitude of the oscillations. In mechanical systems, this
energy loss is usually due to friction (objects rubbing against
one another in some way). For the motion of the mass hanging 
from a spring, a more realistic position function has an 
amplitude that decays exponentially.

NOTE: The commands in this worksheet generate a lot of 
graphics, which may fill up your computer's memory. Before 
proceeding with this Part, you should pull down the Cell menu 
at the top of the screen and select "Remove All Output" (From 
Worksheet). Then you will have to return to the place in the 
worksheet where you left off.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
s: [-2*%e^(-t/5)*cos(t)];
sinterval: [[0,6*%pi]];
slen: length(s)$
t_o: sinterval[1][1]$
t_f: sinterval[slen][2]$
splot: []$

smax: subst(t_o,t,s[1])$
smin: subst(t_o,t,s[1])$

k: 1$
for i: 1 thru 500 do block(
    tcurr: i*(t_f-t_o)/500,
    if (tcurr>=sinterval[k][2] and k+1<=slen) then k: k+1,
    if k<=slen then block(
        curr: subst(tcurr,t,s[k]),
        if curr>smax then smax: curr,
        if curr<smin then smin: curr))$

marginpercent: 10$
xmargin: (t_f-t_o)*(marginpercent/100)$
ymargin: (smax-smin)*(marginpercent/100)$

for i: 1 thru slen do block(
        splot: append(splot,[explicit(s[i],t,sinterval[i][1],sinterval[i][2])]))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=80,
    color=blue,
    splot,
    xrange=[t_o-xmargin, t_f+xmargin],
    yrange=[smin-ymargin, smax+ymargin],
    xlabel="t",
    ylabel="s");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
v: diff(s,t,1);
vinterval: [[0,6*%pi]];
vlen: length(v)$
vplot: []$

vmax: subst(t_o,t,v[1])$
vmin: subst(t_o,t,v[1])$

k: 1$
for i: 1 thru 500 do block(
    tcurr: i*(t_f-t_o)/500,
    if (tcurr>=vinterval[k][2] and k+1<=vlen) then k: k+1,
    if k<=vlen then block(
        curr: subst(tcurr,t,v[k]),
        if curr>vmax then vmax: curr,
        if curr<vmin then vmin: curr))$

marginpercent: 10$
xmargin: (t_f-t_o)*(marginpercent/100)$
ymargin: (vmax-vmin)*(marginpercent/100)$

for i: 1 thru vlen do block(
        vplot: append(vplot,[explicit(v[i],t,vinterval[i][1],vinterval[i][2])]))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=80,
    color=red,
    vplot,
    xrange=[t_o-xmargin, t_f+xmargin],
    yrange=[vmin-ymargin, vmax+ymargin],
    xlabel="t",
    ylabel="s");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
a: diff(v,t,1);
ainterval: [[0,6*%pi]];
alen: length(a)$
aplot: []$

amax: subst(t_o,t,a[1])$
amin: subst(t_o,t,a[1])$

k: 1$
for i: 1 thru 500 do block(
    tcurr: i*(t_f-t_o)/500,
    if (tcurr>=ainterval[k][2] and k+1<=alen) then k: k+1,
    if k<=alen then block(
        curr: subst(tcurr,t,a[k]),
        if curr>amax then amax: curr,
        if curr<amin then amin: curr))$

marginpercent: 10$
xmargin: (t_f-t_o)*(marginpercent/100)$
ymargin: (amax-amin)*(marginpercent/100)$

for i: 1 thru vlen do block(
        aplot: append(aplot,[explicit(a[i],t,ainterval[i][1],ainterval[i][2])]))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=80,
    color=green,
    aplot,
    xrange=[t_o-xmargin, t_f+xmargin],
    yrange=[amin-ymargin, amax+ymargin],
    xlabel="t",
    ylabel="s");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
It is apparent that the amplitudes of the velocity and 
acceleration also decay as time progress.

The motion can be analyzed using the velocity( ), 
acceleration( ), and posvelacc( ) functions.

The following functions will create gif filenames that begin 
with "304".
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
velocity ([-2*%e^(-t/5)*cos(t)], t, [[0,6*%pi]], 1, "304")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
acceleration ([-2*%e^(-t/5)*cos(t)], t, [[0,6*%pi]], 1, "304")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
posvelacc ([-2*%e^(-t/5)*cos(t)], t, [[0,6*%pi]], 1, "304")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: section start ]
Part V: Earthquake
   [wxMaxima: section end   ] */

/* [wxMaxima: comment start ]
Section 3.4

During an earthquake, the floors of building act like masses, 
and the columns between the floors act like elastic springs. 
Consequently, the motion of a single floor can be modeled with
a mass that is equal to the mass of the floor, and a spring 
with a stiffness equal to that provided by the columns above 
and below the floor. In part, the motion will consist of 
oscillations along a straight line. Initially, when the 
earthquake begins, the building is at rest and the ground 
motion puts energy into the building, resulting in back and 
forth oscillations of increasing amplitude. If the earthquake 
is long enough in duration, a steady state condition is 
reached where the amount of energy that goes into the building
is equal to the amount dissipated by friction forces as parts 
of the building rub against one another. During this phase of 
the earthquake the building sways with oscillations of 
constant amplitude. Then, when the earthquake stops the 
oscillations will decay and the building will eventually come 
to rest (hopefully with little or no damage). 

The following function can be used to describe the position of
the floor in a building during an earthquake. To see where we 
got this function see the note at the end of this Part.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
s: [2*%e^(-t/5)*cos(0.9798*t)-2*cos(t)+0.40825*%e^(-t/5)*sin(0.9798*t),
    (-698.8)*%e^(-t/5)*cos(0.9798*t) + 478.3*%e^(-t/5)*sin(0.9798*t)];
sinterval: [[0,30],[30,70]];
slen: length(s)$
t_o: sinterval[1][1]$
t_f: sinterval[slen][2]$
splot: []$

smax: subst(t_o,t,s[1])$
smin: subst(t_o,t,s[1])$

k: 1$
for i: 1 thru 500 do block(
    tcurr: i*(t_f-t_o)/500,
    if (tcurr>=sinterval[k][2] and k+1<=slen) then k: k+1,
    if k<=slen then block(
        curr: subst(tcurr,t,s[k]),
        if curr>smax then smax: curr,
        if curr<smin then smin: curr))$

marginpercent: 10$
xmargin: (t_f-t_o)*(marginpercent/100)$
ymargin: (smax-smin)*(marginpercent/100)$

for i: 1 thru slen do block(
        splot: append(splot,[explicit(s[i],t,sinterval[i][1],sinterval[i][2])]))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=80,
    color=blue,
    splot,
    xrange=[t_o-xmargin, t_f+xmargin],
    yrange=[smin-ymargin, smax+ymargin],
    xlabel="t",
    ylabel="s");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
And, the graphs of the velocity and acceleration look like 
this.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
v: diff(s,t,1);
vinterval: [[0,30],[30,70]];
vlen: length(v)$
vplot: []$

vmax: subst(t_o,t,v[1])$
vmin: subst(t_o,t,v[1])$

k: 1$
for i: 1 thru 500 do block(
    tcurr: i*(t_f-t_o)/500,
    if (tcurr>=vinterval[k][2] and k+1<=vlen) then k: k+1,
    if k<=vlen then block(
        curr: subst(tcurr,t,v[k]),
        if curr>vmax then vmax: curr,
        if curr<vmin then vmin: curr))$

marginpercent: 10$
xmargin: (t_f-t_o)*(marginpercent/100)$
ymargin: (vmax-vmin)*(marginpercent/100)$

for i: 1 thru vlen do block(
        vplot: append(vplot,[explicit(v[i],t,vinterval[i][1],vinterval[i][2])]))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=80,
    color=red,
    vplot,
    xrange=[t_o-xmargin, t_f+xmargin],
    yrange=[vmin-ymargin, vmax+ymargin],
    xlabel="t",
    ylabel="s");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
a: diff(v,t,1);
ainterval: [[0,30],[30,70]];
alen: length(a)$
aplot: []$

amax: subst(t_o,t,a[1])$
amin: subst(t_o,t,a[1])$

k: 1$
for i: 1 thru 500 do block(
    tcurr: i*(t_f-t_o)/500,
    if (tcurr>=ainterval[k][2] and k+1<=alen) then k: k+1,
    if k<=alen then block(
        curr: subst(tcurr,t,a[k]),
        if curr>amax then amax: curr,
        if curr<amin then amin: curr))$

marginpercent: 10$
xmargin: (t_f-t_o)*(marginpercent/100)$
ymargin: (amax-amin)*(marginpercent/100)$

for i: 1 thru vlen do block(
        aplot: append(aplot,[explicit(a[i],t,ainterval[i][1],ainterval[i][2])]))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
wxdraw2d(
    nticks=80,
    color=green,
    aplot,
    xrange=[t_o-xmargin, t_f+xmargin],
    yrange=[amin-ymargin, amax+ymargin],
    xlabel="t",
    ylabel="s");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
The velocity( ), acceleration( ), and posvelacc( ) functions 
depict the motion.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
velocity ([2*%e^(-t/5)*cos(0.9798*t)-2*cos(t)+0.40825*%e^(-t/5)*sin(0.9798*t),
    (-698.8)*%e^(-t/5)*cos(0.9798*t) + 478.3*%e^(-t/5)*sin(0.9798*t)], 
    t, [[0,30],[30,70]], 1, "305")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
acceleration ([2*%e^(-t/5)*cos(0.9798*t)-2*cos(t)+0.40825*%e^(-t/5)*sin(0.9798*t),
    (-698.8)*%e^(-t/5)*cos(0.9798*t) + 478.3*%e^(-t/5)*sin(0.9798*t)], 
    t, [[0,30],[30,70]], 1, "305")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
posvelacc ([2*%e^(-t/5)*cos(0.9798*t)-2*cos(t)+0.40825*%e^(-t/5)*sin(0.9798*t),
    (-698.8)*%e^(-t/5)*cos(0.9798*t) + 478.3*%e^(-t/5)*sin(0.9798*t)], 
    t, [[0,30],[30,70]], 1, "305")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Note: The position function used for the earthquake motion 
above was obtained by solving the following initial value 
problem. The ground motion excitation occurs between t = 0 and
t = 30, and after t = 30 the motion is free or unforced 
vibration.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
kill(s)$
atemp: 'diff(s,t,2) + .4*'diff(s,t)+ s = .8*sin(t)$
atemp2: ode2(atemp, s, t)$
atemp3: ic2(atemp2, t=0, s=0, 'diff(s,t)=0);
btemp: 'diff(s,t,2) + .4*'diff(s,t)+ s = 0$
btemp2: ode2(btemp, s, t)$
btemp3: ic2(btemp2, t=0, s=0, 'diff(s,t)=0);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
f2: (-698.8)*%e^(-t/5)*cos(0.9798*t) + 478.3*%e^(-t/5)*sin(0.9798*t)$
wxdraw2d(
    nticks=80,
    color=blue,
    explicit(rhs(atemp3),t,0,30),
    explicit(f2,t,30,50));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: section start ]
Part VI: Some Help with Your Homework
   [wxMaxima: section end   ] */

/* [wxMaxima: comment start ]
Section 3.3, Exercise 3

The special commands in this section may be helpful with some 
of the homework Exercises from Chapter 3, Section 3. For 
example, we use them to visualize some of the motion for 
Exercise 3.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
s:-t^3+3*t^2-3*t;
v: diff(s,t,1);
a: diff(v,t,1);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
velocity ([t^3+3*t^2-3*t], t, [[0,3]], 0, "306")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
acceleration ([t^3+3*t^2-3*t], t, [[0,3]], 0, "306")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
posvelacc ([t^3+3*t^2-3*t], t, [[0,3]], 0, "306")$
/* [wxMaxima: input   end   ] */

/* Maxima can't load/batch files which end with a comment! */
"Created with wxMaxima"$