lab_fft.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Applied Seismology, GEOS 626, University of Alaska Fairbanks\n",
    "\n",
    "- spectral analysis of the motion of the Sun due to planetary positions."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget\n",
    "\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.io\n",
    "import warnings\n",
    "\n",
    "from lib_geos import fftvec\n",
    "from lib_seis import markp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# script settings\n",
    "\n",
    "warnings.filterwarnings('ignore')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# USER PARAMETERS\n",
    "\n",
    "# switch between using long time series and short time series\n",
    "longdata = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# synthetically derived sample points of the motion of our sun due to planet positions\n",
    "\n",
    "solarshort = np.loadtxt( './data/solarshort.dat' ) # Make sure that solarshort.dat is in your working directory\n",
    "\n",
    "ts = solarshort[:,0]\n",
    "xs = solarshort[:,1]\n",
    "ys = solarshort[:,2]\n",
    "#cs = [xs + i*ys]\n",
    "cs = np.vectorize(complex)(xs, ys)\n",
    "\n",
    "solarlong = np.loadtxt( './data/solarlong.dat' ) # Make sure that solarlong.dat is in your working directory\n",
    "tL = solarlong[:,0]\n",
    "xL = solarlong[:,1]\n",
    "yL = solarlong[:,2]\n",
    "#cL = [xL + i*yL]\n",
    "cL = np.vectorize(complex)(xL, yL)\n",
    "numt = len(tL)\n",
    "\n",
    "fig1 = plt.figure()\n",
    "plt.rc('xtick',labelsize=8)\n",
    "plt.rc('ytick',labelsize=8)\n",
    "\n",
    "ax1 = plt.subplot2grid((3, 2), (0, 0), colspan=1)\n",
    "ax1.plot(xs, ys,'r', linewidth=0.5)\n",
    "plt.xlim(-0.0025,0.0025)\n",
    "plt.xlabel('x distance, AU', fontsize=9)\n",
    "plt.ylabel('y distance, AU', fontsize=9)\n",
    "plt.ticklabel_format(axis='both',style='sci',scilimits=(0,0))\n",
    "ax1.set_aspect('equal')\n",
    "\n",
    "ax2 = plt.subplot2grid((3, 2), (0, 1), colspan=1)\n",
    "ax2.plot(xL, yL,'b', linewidth=0.5)\n",
    "plt.xlim(-0.0025,0.0025)\n",
    "plt.xlabel('x distance, AU', fontsize=9)\n",
    "plt.ylabel('y distance, AU', fontsize=9)\n",
    "plt.ticklabel_format(axis='both',style='sci',scilimits=(0,0))\n",
    "ax2.set_aspect('equal')\n",
    "\n",
    "ax3 = plt.subplot2grid((3, 2), (1, 0), colspan=2)\n",
    "ax3.plot(tL, xL,'b', linewidth=0.5)\n",
    "ax3.plot(ts, xs,'r--', linewidth=0.5)\n",
    "plt.xlim(0,max(tL))\n",
    "plt.xlabel('time, years', fontsize=9)\n",
    "plt.ylabel('x distance, AU', fontsize=9)\n",
    "plt.ticklabel_format(axis='y',style='sci',scilimits=(0,0))\n",
    "\n",
    "ax4 = plt.subplot2grid((3, 2), (2, 0), colspan=2)\n",
    "ax4.plot(tL, yL,'b', linewidth=0.5)\n",
    "ax4.plot(ts, ys,'r--', linewidth=0.5)\n",
    "plt.xlim(0,max(tL))\n",
    "plt.xlabel('time, years', fontsize=9)\n",
    "plt.ylabel('y distance, AU', fontsize=9)\n",
    "plt.ticklabel_format(axis='y',style='sci',scilimits=(0,0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if longdata==0:\n",
    "    t = ts; c = cs; ftag = 'short'\n",
    "else:\n",
    "    t = tL; c = cL; ftag = 'long'\n",
    "\n",
    "n = len(t)\n",
    "\n",
    "# data type, length, size\n",
    "print('t:',t.dtype, t.size, t.shape)\n",
    "# data type, length, size\n",
    "print('c:',c.dtype, c.size, c.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sidereal orbit periods (years)\n",
    "# http://ssd.jpl.nasa.gov/?planet_phys_par\n",
    "\n",
    "periods = np.array([0.2408467,0.61519726,1.0000174,1.8808476,11.862615,29.447498,84.016846,164.79132,247.92065])\n",
    "nump = len(periods)\n",
    "freqs = []\n",
    "omegs = []\n",
    "\n",
    "print('frequency for Mercury:',1/periods[0],' 1/years\\n')\n",
    "\n",
    "for i in range(nump):\n",
    "    freq=(1/periods[i])\n",
    "    freqs.append(freq)\n",
    "    omegs.append(2*math.pi*float(freq))\n",
    "\n",
    "stvec = ['Mer','Ven','Ear','Mar','Jup','Sat','Ura','Nep','Plu']\n",
    "\n",
    "print('ALL TIME UNITS IN YEARS')\n",
    "print('%22s %14s %20s ' %('period', 'frequency',  'angular frequency'))\n",
    "\n",
    "for ii in range(nump):\n",
    "    print('%8s %14s %12s %14s ' %(stvec[ii],round(periods[ii],6),round(freqs[ii],6),round(omegs[ii],6)))\n",
    "\n",
    "#          period     frequency  angular frequency\n",
    "# Mer       0.240847    4.152019   26.087903\n",
    "# Ven       0.615197    1.625495   10.213286\n",
    "# Ear       1.000017    0.999983    6.283076\n",
    "# Mar       1.880848    0.531675    3.340614\n",
    "# Jup      11.862615    0.084298    0.529663\n",
    "# Sat      29.447498    0.033959    0.213369\n",
    "# Ura      84.016846    0.011902    0.074785\n",
    "# Nep     164.791320    0.006068    0.038128\n",
    "# Plu     247.920650    0.004034    0.025344"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# what does fft do?\n",
    "np.fft.fft?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#==========================================================================\n",
    "# FFT ANALYSIS HERE\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# example using markp (click on the plot to see values displayed)\n",
    "\n",
    "p = np.arange(2,9.0,0.1)\n",
    "y = np.sin(p)\n",
    "fig = plt.figure(); plt.plot(p,y,'-')\n",
    "cid = fig.canvas.mpl_connect('button_press_event', markp) #displays the value 1/x for a selected point (x,y)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "name": "python3",
   "language": "python",
   "display_name": "Python 3 (ipykernel)"
  },
  "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}