{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A worksheet to estimate spring and damping constants from experimental data.\n",
    "\n",
    "#The Data: Load in the data, 1460 data points, the position (meters) of the\n",
    "#mass every 1/50 of a second, starting at time t = 0.82 seconds.\n",
    "load('spring_mass_data_clean.sage');\n",
    "N = len(udat0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The data is in the array \"udat0\". A plot:\n",
    "scatter_plot(udat0,facecolor='red',markersize=10).show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Center the data vertically and rescale so time starts at t = 0. First\n",
    "#compute average of data (might be better to use integral number\n",
    "#of cycles):\n",
    "uave = add(udat0[i][1] for i in range(N))/N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now subtract uave, make data start at t = 0:\n",
    "udat = [[udat0[i][0]-0.82,udat0[i][1]-uave] for i in range(N)];\n",
    "Tmax = udat[N-1][0]; #Time of last sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A plot of the recentered data:\n",
    "plt1 = scatter_plot(udat,facecolor='red',markersize=10);\n",
    "show(plt1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Estimating the spring and damping constants: We will fit a\n",
    "#function y(t) of the form\n",
    "var('t d1 alpha omega');\n",
    "y(t,d1,alpha,omega) = d1*exp(-alpha*t)*cos(omega*t);\n",
    "#to the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Based on the plot above we can guess that the initial amplitude\n",
    "#d1 is about 0.05. We can estimate alpha by the rate of decay of the\n",
    "#amplitude of the oscillations. For example, at time t = 25 the amplitude\n",
    "#is down to about 0.035, so 0.05*exp(-alpha*25) = 0.035, which leads to\n",
    "#alpha equal to about 0.014 (solve 0.05*exp(-alpha*25) = 0.035 for alpha).\n",
    "\n",
    "#We can estimate omega by estimating the period of the motion, e.g.,\n",
    "#count how many complete oscillations the mass undergoes during the\n",
    "#approximate 29 second data set. We can then plot y(t) with these\n",
    "#estimated values, as\n",
    "plt2 = plot(y(t,0.05,0.014,7.0),(t,0,Tmax),color='red');\n",
    "show(plt1+plt2);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#It may help to plot on a smaller time range, at first:\n",
    "show(plt1+plt2,xmin=0,xmax=5,ymin=-0.05,ymax=0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Obviously some adjustment in omega and perhaps the other parameters\n",
    "#is in order.  Adjust d1, alpha, and omega to obtain the best (visual)\n",
    "#fit possible, then use the formulas in Modeling Exercise 6.3.5 to\n",
    "#estimate the spring and damping constant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.2",
   "language": "sage",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}