{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np                    # matrices & etc\n",
    "from matplotlib import pyplot as plt  # general plotting \n",
    "from matplotlib import patches        # plotting shapes\n",
    "from datetime import date             # date and time arithmetic\n",
    "from math import exp, pi, sin, cos, tan, sqrt, floor, log10, nan, sqrt   # math functions\n",
    "import scipy.linalg as la             # linear algebra functions\n",
    "import os                             # operating system\n",
    "import matplotlib.cm as cm\n",
    "from matplotlib.colors import ListedColormap\n",
    "import scipy.signal as sg\n",
    "import scipy.stats as st\n",
    "import scipy.sparse.linalg as las\n",
    "import matplotlib\n",
    "\n",
    "# function to make a numpy N-by-1 column vector\n",
    "# c=eda_cvec(v1, v2, ...) from a list of several\n",
    "# array-like entities v1, v2, including a number\n",
    "# a list of numbers, a tuple of numbers, an N-by-0 np array\n",
    "# and a N-by-1 np array. The function\n",
    "# also insures that, if v is an np array, that\n",
    "# c is a copy, as contrasted to a view, of it\n",
    "# It promotes integers to floats, and integers\n",
    "# and floats to complex, by context.\n",
    "# This version concatenates many argments,\n",
    "# whereas c=eda_cvec1(v1) takes just one argiment.\n",
    "# I recommend always using eda_cvec(v1, v2, ...)\n",
    "def gda_cvec(*argv):\n",
    "    t = int;\n",
    "    Nt = 0;\n",
    "    for a in argv:\n",
    "        v = gda_cvec1(a);\n",
    "        N,M = np.shape(v);\n",
    "        Nt = Nt + N;\n",
    "        if( N==0 ):\n",
    "            continue; # skip vector of zero length\n",
    "        if (t==int) and isinstance(v[0,0],float):\n",
    "            t=float;\n",
    "        elif isinstance(v[0,0],complex):\n",
    "            t=complex;\n",
    "    w = np.zeros((Nt,1),dtype=t);\n",
    "    Nt = 0;\n",
    "    for a in argv:\n",
    "        v = gda_cvec1(a);\n",
    "        N,M = np.shape(v);\n",
    "        w[Nt:Nt+N,0:1] = v;  # patch 20230418 was #w[Nt:Nt+N,0] = v[0:N,0];\n",
    "        Nt = Nt + N;\n",
    "    return w;\n",
    "\n",
    "# function to make a numpy N-by-1 column vector\n",
    "# c=gda_cvec1(v) from entity v that is array-like,\n",
    "# including a number, a list of numbers, a tuple\n",
    "# of numbers, an N-by-0 np array and a N-by1 np array.\n",
    "# It promotes integers to floats, and integers\n",
    "# and floats to complex, by context. The function\n",
    "# also insures that, if v is an np array, that\n",
    "# c is a copy, as contrasted to a view, of it.\n",
    "# This version takes just one input argment.\n",
    "# whereas c=gda_cvec(v1,v2,...) concatenates\n",
    "# many argiments.\n",
    "def gda_cvec1(v):\n",
    "    if isinstance(v, int) or isinstance(v, np.int32):\n",
    "        w = np.zeros((1,1),dtype=int);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, float):\n",
    "        w = np.zeros((1,1),dtype=float);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, complex):\n",
    "        w = np.zeros((1,1),dtype=complex);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, np.ndarray):\n",
    "        s = np.shape(v);\n",
    "        if len(s) == 1:\n",
    "            return np.copy(np.reshape(v,(s[0],1)));\n",
    "        else:\n",
    "            [r,c]=s;\n",
    "            if( c==1 ):\n",
    "                return(np.copy(v));\n",
    "            elif(r==1):\n",
    "                return(np.copy(v.T));\n",
    "            else:\n",
    "                raise TypeError(\"gda_cvec: %d by %d ndarray not allowed\" % (r, c));\n",
    "    elif isinstance(v, list):\n",
    "        r = len(v);\n",
    "        t = int;\n",
    "        for vi in v:\n",
    "            if isinstance(vi,int) or isinstance(vi, np.int32): #patch v->vi 20230418\n",
    "                pass;\n",
    "            elif isinstance(vi,float):\n",
    "                t=float;\n",
    "            elif isinstance(vi,complex):\n",
    "                t=complex;\n",
    "                break;\n",
    "            else:\n",
    "                raise TypeError(\"gda_cvec: list contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = np.array(v); # patch v -> np.array(v)\n",
    "        return w;\n",
    "    elif isinstance(v, tuple):\n",
    "        r = len(v);\n",
    "        t = int;\n",
    "        for vi in v:\n",
    "            if isinstance(vi,int) or isinstance(vi, np.int32): #patch v->vi 20230418\n",
    "                pass;\n",
    "            elif isinstance(vi,float):\n",
    "                t=float;\n",
    "            elif isinstance(vi,complex):\n",
    "                t=complex;\n",
    "                break;\n",
    "            else:\n",
    "                raise TypeError(\"gda_cvec: tuple contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = np.array(list(v)); # patch v -> np.array(list(v));\n",
    "        return w;\n",
    "    else:\n",
    "        raise TypeError(\"gda_cvec: %s not supported\" % type(v));\n",
    "\n",
    "        \n",
    "# gda_draw function makes a \"pictorial matrix equation\"\n",
    "# arguments are vectors, matrices and strings\n",
    "# which are plotted in the order that the appear\n",
    "# except that strings starting with 'title ' are plotted\n",
    "# under the subseqeunt matrix or vector\n",
    "# always returns a status of 1\n",
    "def gda_draw(*argv):\n",
    "    DOCOLOR=True;\n",
    "    if( DOCOLOR ):\n",
    "        bwcmap = matplotlib.colormaps['jet'];\n",
    "    else:\n",
    "        bw = np.zeros((256,4));\n",
    "        v = 0.9*(256 - np.linspace( 0, 255, 256 ))/255;\n",
    "        bw[:,0] = v;\n",
    "        bw[:,1] = v;\n",
    "        bw[:,2] = v;\n",
    "        bw[:,3] = np.ones(256);\n",
    "        bwcmap = ListedColormap(bw);\n",
    "    # size of plot\n",
    "    W = 16;\n",
    "    H = 4;\n",
    "    fig1 = plt.figure(1);\n",
    "    # figsize width and height in inches\n",
    "    fig1.set_size_inches(W,H);\n",
    "    ax1 = plt.subplot(1,1,1);\n",
    "    plt.axis([0, W, -H/2, H/2]);\n",
    "    plt.axis('off');\n",
    "    LM = W/6;    # matrix width and heoght\n",
    "    LV = W/40;   # vector width\n",
    "    FS = 0.12;    # character width\n",
    "    TO = 0.4;    # title vertical offset\n",
    "    SP = 0.2;    # space between objects\n",
    "    LS = 0.2;    # leading space\n",
    "    p = LS; # starting x-position\n",
    "    istitle=0; # flags presence of a title\n",
    "    for a in argv:\n",
    "        if isinstance(a,np.ndarray):\n",
    "            sh = np.shape(a);\n",
    "            if len(sh) == 1:  # conversion to nx1 array\n",
    "                n = sh[0];\n",
    "                m = 1;\n",
    "                ap = a;\n",
    "                a = np.zeros((n,1));\n",
    "                a[:,0] = ap;\n",
    "            else:\n",
    "                n = sh[0];\n",
    "                m = sh[1];\n",
    "            if m==1:\n",
    "                pold=p;\n",
    "                left=p;\n",
    "                right=p+LV;\n",
    "                bottom=-LM/2;\n",
    "                top=LM/2;\n",
    "                plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );\n",
    "                p = p+LV;\n",
    "                pm = (p+pold)/2;\n",
    "                if istitle:\n",
    "                    plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');\n",
    "                    istitle=0;\n",
    "                p = p+SP;\n",
    "            else:\n",
    "                pold=p;\n",
    "                left=p;\n",
    "                right=p+LM;\n",
    "                bottom=-LM/2;\n",
    "                top=LM/2;\n",
    "                plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );\n",
    "                p = p+LM;\n",
    "                pm = (p+pold)/2;\n",
    "                if istitle:\n",
    "                    plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');\n",
    "                    istitle=0;\n",
    "                p = p+SP;\n",
    "        elif isinstance(a,str):\n",
    "            ns = len(a);\n",
    "            istitle=0;\n",
    "            if( ns>=6 ):\n",
    "                if 'title ' in a[0:6]:\n",
    "                    istitle=1;\n",
    "                    titlestr=a[6:];\n",
    "            if( istitle != 1):\n",
    "                plt.text(p,0,a);\n",
    "                p = p + ns*FS + SP;\n",
    "    plt.show();\n",
    "    return 1;\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Check normalization of eigenfunctions. Normalization should be close to unity.\n",
      "min norm 0.9999999999999989  max norm  1.0000000000000009\n",
      "Check eigenfrequencies. Errors should be close to zero.\n",
      "min error 4.2131957854972795e-11  max error  0.004160475109148909\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 121980 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Caption: Pressure field (colors) for source at 64.00, with observer (triangle) at 25.50\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 31348 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Caption: Pressure field for observer at 25.50\n"
     ]
    }
   ],
   "source": [
    "# Groups should use NEXT section, not this one\n",
    "# Excitation of an Organ Pipe by a impulsive force\n",
    "\n",
    "ipx = 51; # observer location (as an index into x)\n",
    "ix0 = 128; # source location (as an index into x)\n",
    "sf = 2.0; # width of the source\n",
    "\n",
    "rho = 1.0;\n",
    "lam = 1.0;\n",
    "c2 = lam/rho;\n",
    "c = sqrt(c2);\n",
    "\n",
    "# standard set-up of position x and time t\n",
    "Nx=512;\n",
    "Dx=0.5;\n",
    "x=gda_cvec(Dx*np.linspace(0,Nx-1,Nx));\n",
    "xmin=0.0;\n",
    "xmax = Dx*(Nx-1);\n",
    "\n",
    "Nt=512;\n",
    "Dt=0.5;\n",
    "t=gda_cvec(Dt*np.linspace(0,Nt-1,Nt));\n",
    "tmin=0.0;\n",
    "tmax = Dt*(Nt-1);\n",
    "\n",
    "# eigenfrequencies and eigenfunction of organ pipe\n",
    "Nw = 100;\n",
    "n = gda_cvec( np.linspace(1,Nw,Nw) );\n",
    "wn = n*pi*c/xmax;        # eigenfrequencies w(n)\n",
    "pn = sqrt(2*lam/xmax) * np.sin( (pi/xmax)*np.matmul(x, n.T) );  # eigenfunctions p(ix,n)\n",
    "\n",
    "# 1st derivative operator\n",
    "myr = gda_cvec( 1, -1, np.zeros((Nx-2)) );\n",
    "myc = gda_cvec( 1, np.zeros((Nx-1)) );\n",
    "D1 = la.toeplitz( myr.ravel(), myc.ravel() );\n",
    "D1 = D1 / Dx;\n",
    "\n",
    "# 2nd derivative operator\n",
    "myr = gda_cvec( -2, 1, np.zeros((Nx-2)) );\n",
    "myc = gda_cvec( -2, 1, np.zeros((Nx-2)) );\n",
    "D2 = la.toeplitz( myr.ravel(), myc.ravel() );\n",
    "D2[0,0]=1.0;\n",
    "D2[0,1]=0.0;\n",
    "D2[Nx-1,Nx-1]=1.0;\n",
    "D2[Nx-1,Nx-2]=0.0;\n",
    "D2 = D2 / (Dx**2);\n",
    "\n",
    "# check noralization of eigenfunctions\n",
    "norm = np.zeros((Nw,1));     \n",
    "for iw in range(Nw):\n",
    "    p = pn[0:Nx,iw:iw+1];\n",
    "    norm[iw,0] = Dx * np.sum( np.power(p,2) / lam );\n",
    "print('Check normalization of eigenfunctions. Normalization should be close to unity.');\n",
    "print( 'min norm', np.min(norm), ' max norm ', np.max(norm) );\n",
    "\n",
    "# check eigenfrequencies\n",
    "# d2/dt2 pn = lam d/dx (1/rho) d/dx pn\n",
    "# -wn2 pn = lam d/dx (1/rho) d/dx pn\n",
    "# -wn2 pn = c2 d2/dx2 pn\n",
    "# wn2 pn + c2 d2/dx2 pn = 0\n",
    "error = np.zeros((Nw,1));     \n",
    "for iw in range(Nw):\n",
    "    # error = wn2 pn + c2 d2/dx2 pn\n",
    "    p = pn[0:Nx,iw:iw+1];\n",
    "    error[iw,0] = np.max( np.abs( (wn[iw,0]**2)*p + c2 * np.matmul( D2, p ) ) );\n",
    "print('Check eigenfrequencies. Errors should be close to zero.');\n",
    "print( 'min error', np.min(error), ' max error ', np.max(error) );\n",
    "\n",
    "# assign excitation function for force, f\n",
    "sf2 = sf**2;\n",
    "x0 = x[ix0,0];\n",
    "f = (1.0/(sqrt(2*pi)*sf)) * np.exp( -0.5*np.power(x-x0,2)/sf2 ); # gaussian\n",
    "# the source in the differential equation is actually -c2 df/dx\n",
    "excit = np.zeros((Nw,1));     \n",
    "for iw in range(Nw):\n",
    "    # excit = (1/lam) (integral) pn(x) f(s)\n",
    "    p = pn[0:Nx,iw:iw+1];\n",
    "    g = -c2 * np.matmul(D1,f);\n",
    "    excit[iw,0] = (1.0/lam) * Dx * np.sum( np.multiply(p,g) );\n",
    "\n",
    "# compute pressure by performing mode sum\n",
    "P = np.zeros((Nt,Nx));\n",
    "for iw in range(Nw):\n",
    "    p = pn[0:Nx,iw:iw+1];\n",
    "    tfun = np.sin( wn[iw,0] * t );\n",
    "    P = P + (1.0/wn[iw,0]) * excit[iw,0] * np.matmul(tfun,p.T);\n",
    "    \n",
    "# extract pressure time series for one position\n",
    "px = P[0:Nt,ipx:ipx+1];\n",
    "pxabsmax=np.max(np.abs(px));\n",
    "\n",
    "fig1 = plt.figure(figsize=(12,5));   # smallish figure\n",
    "ax1 = plt.subplot(1, 1, 1);\n",
    "plt.axis( [ 0, xmax, 0, tmax] );\n",
    "mycm = plt.get_cmap('jet');\n",
    "plt.axis( [xmin, xmax, tmin, tmax ] );\n",
    "plt.imshow(P,cmap=mycm,aspect='auto',extent=(xmin,xmax,tmax,tmin));\n",
    "plt.plot(x[ipx,0],0.0,'k^',lw=4);\n",
    "plt.colorbar();\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('t');\n",
    "plt.show();\n",
    "print(\"Caption: Pressure field (colors) for source at %.2f, with observer (triangle) at %.2f\" % (x0,x[ipx,0]) );\n",
    "\n",
    "\n",
    "fig2 = plt.figure(figsize=(12,5));   # smallish figure\n",
    "ax1 = plt.subplot(1, 1, 1);\n",
    "plt.axis( [ 0, tmax, 0, tmax] );\n",
    "mycm = plt.get_cmap('jet');\n",
    "plt.axis( [tmin, tmax, -1.1*pxabsmax, 1.1*pxabsmax ] );\n",
    "plt.plot(t,px,'k-');\n",
    "plt.xlabel('t');\n",
    "plt.ylabel('p(t)');\n",
    "plt.show();\n",
    "print(\"Caption: Pressure field for observer at %.2f\" % (x[ipx,0]) );\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Check normalization of eigenfunctions. Normalization should be close to unity.\n",
      "min norm 0.9999999999999989  max norm  1.0000000000000009\n",
      "Check eigenfrequencies. Errors should be close to zero.\n",
      "min error 4.2131957854972795e-11  max error  0.004160475109148909\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 123736 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Caption: Pressure field (colors) for source at 64.00, with observer (triangle) at 25.50\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 43596 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Caption: Pressure field for observer at 25.50\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 60728 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Caption: ASD of pressure field for observer at 25.50\n"
     ]
    }
   ],
   "source": [
    "# Groups should use this section\n",
    "# Excitation of an Organ Pipe by a impulsive force\n",
    "\n",
    "ipx = 51; # observer location (as an index into x) (was 51)\n",
    "ix0 = 128; # source location (as an index into x) (was 128)\n",
    "Dt = 0.5; # time increment (was 0.5)\n",
    "sf = 2.0; # width of the source (was 2.0)\n",
    "\n",
    "rho = 1.0;  # density\n",
    "lam = 1.0;  # compressibility\n",
    "\n",
    "c2 = lam/rho;\n",
    "c = sqrt(c2);\n",
    "\n",
    "# standard set-up of position x and time t\n",
    "Nx=512;\n",
    "Dx=0.5;\n",
    "x=gda_cvec(Dx*np.linspace(0,Nx-1,Nx));\n",
    "xmin=0.0;\n",
    "xmax = Dx*(Nx-1);\n",
    "\n",
    "Nt=1024*4;\n",
    "t=gda_cvec(Dt*np.linspace(0,Nt-1,Nt));\n",
    "tmin=0.0;\n",
    "tmax = Dt*(Nt-1);\n",
    "\n",
    "# eigenfrequencies and eigenfunction of organ pipe\n",
    "Nw = 100;\n",
    "n = gda_cvec( np.linspace(1,Nw,Nw) );\n",
    "wn = n*pi*c/xmax;        # eigenfrequencies w(n)\n",
    "pn = sqrt(2*lam/xmax) * np.sin( (pi/xmax)*np.matmul(x, n.T) );  # eigenfunctions p(ix,n)\n",
    "\n",
    "# 1st derivative operator\n",
    "myr = gda_cvec( 1, -1, np.zeros((Nx-2)) );\n",
    "myc = gda_cvec( 1, np.zeros((Nx-1)) );\n",
    "D1 = la.toeplitz( myr.ravel(), myc.ravel() );\n",
    "D1 = D1 / Dx;\n",
    "\n",
    "# 2nd derivative operator\n",
    "myr = gda_cvec( -2, 1, np.zeros((Nx-2)) );\n",
    "myc = gda_cvec( -2, 1, np.zeros((Nx-2)) );\n",
    "D2 = la.toeplitz( myr.ravel(), myc.ravel() );\n",
    "D2[0,0]=1.0;\n",
    "D2[0,1]=0.0;\n",
    "D2[Nx-1,Nx-1]=1.0;\n",
    "D2[Nx-1,Nx-2]=0.0;\n",
    "D2 = D2 / (Dx**2);\n",
    "\n",
    "# check noralization of eigenfunctions\n",
    "norm = np.zeros((Nw,1));     \n",
    "for iw in range(Nw):\n",
    "    p = pn[0:Nx,iw:iw+1];\n",
    "    norm[iw,0] = Dx * np.sum( np.power(p,2) / lam );\n",
    "print('Check normalization of eigenfunctions. Normalization should be close to unity.');\n",
    "print( 'min norm', np.min(norm), ' max norm ', np.max(norm) );\n",
    "\n",
    "# check eigenfrequencies\n",
    "# d2/dt2 pn = lam d/dx (1/rho) d/dx pn\n",
    "# -wn2 pn = lam d/dx (1/rho) d/dx pn\n",
    "# -wn2 pn = c2 d2/dx2 pn\n",
    "# wn2 pn + c2 d2/dx2 pn = 0\n",
    "error = np.zeros((Nw,1));     \n",
    "for iw in range(Nw):\n",
    "    # error = wn2 pn + c2 d2/dx2 pn\n",
    "    p = pn[0:Nx,iw:iw+1];\n",
    "    error[iw,0] = np.max( np.abs( (wn[iw,0]**2)*p + c2 * np.matmul( D2, p ) ) );\n",
    "print('Check eigenfrequencies. Errors should be close to zero.');\n",
    "print( 'min error', np.min(error), ' max error ', np.max(error) );\n",
    "\n",
    "# assign excitation function for force, f\n",
    "sf2 = sf**2;\n",
    "x0 = x[ix0,0];\n",
    "f = (1.0/(sqrt(2*pi)*sf)) * np.exp( -0.5*np.power(x-x0,2)/sf2 ); # gaussian\n",
    "# the source in the differential equation is actually -c2 df/dx\n",
    "excit = np.zeros((Nw,1));     \n",
    "for iw in range(Nw):\n",
    "    # excit = (1/lam) (integral) pn(x) f(s)\n",
    "    p = pn[0:Nx,iw:iw+1];\n",
    "    g = -c2 * np.matmul(D1,f);\n",
    "    excit[iw,0] = (1.0/lam) * Dx * np.sum( np.multiply(p,g) );\n",
    "\n",
    "# compute pressure by performing mode sum\n",
    "P = np.zeros((Nt,Nx));\n",
    "for iw in range(Nw):\n",
    "    p = pn[0:Nx,iw:iw+1];\n",
    "    tfun = np.sin( wn[iw,0] * t );\n",
    "    P = P + (1.0/wn[iw,0]) * excit[iw,0] * np.matmul(tfun,p.T);\n",
    "    \n",
    "# extract pressure time series for one position\n",
    "px = P[0:Nt,ipx:ipx+1];\n",
    "pxabsmax=np.max(np.abs(px));\n",
    "\n",
    "fig1 = plt.figure(figsize=(12,5));   # smallish figure\n",
    "ax1 = plt.subplot(1, 1, 1);\n",
    "plt.axis( [ 0, xmax, 0, tmax] );\n",
    "mycm = plt.get_cmap('jet');\n",
    "plt.axis( [xmin, xmax, tmin, tmax ] );\n",
    "plt.imshow(P,cmap=mycm,aspect='auto',extent=(xmin,xmax,tmax,tmin));\n",
    "plt.plot(x[ipx,0],0.0,'k^',lw=4);\n",
    "plt.colorbar();\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('t');\n",
    "plt.show();\n",
    "print(\"Caption: Pressure field (colors) for source at %.2f, with observer (triangle) at %.2f\" % (x0,x[ipx,0]) );\n",
    "\n",
    "\n",
    "fig2 = plt.figure(figsize=(12,5));   # smallish figure\n",
    "ax1 = plt.subplot(1, 1, 1);\n",
    "plt.axis( [ 0, tmax, 0, tmax] );\n",
    "mycm = plt.get_cmap('jet');\n",
    "plt.axis( [tmin, tmax, -1.1*pxabsmax, 1.1*pxabsmax ] );\n",
    "plt.plot(t,px,'k-');\n",
    "plt.xlabel('t');\n",
    "plt.ylabel('p(t)');\n",
    "plt.show();\n",
    "print(\"Caption: Pressure field for observer at %.2f\" % (x[ipx,0]) );\n",
    "\n",
    "# amplitude spectral density\n",
    "Nf = int(Nt/2)+1;\n",
    "fmax = 1.0/(2.0*Dt);\n",
    "f = gda_cvec(np.linspace(0.0,fmax,Nf));\n",
    "Df = f[1,0]-f[0,0];\n",
    "pxt = Dt*np.fft.fft(px,axis=0);\n",
    "asd = sqrt(2.0/(Nt*Dt)) * np.abs(pxt[0:Nf,0:1]);\n",
    "asdmax = np.max(asd);\n",
    "\n",
    "w0 = pi*c2/xmax;\n",
    "f0 = w0/(2.0*pi);\n",
    "\n",
    "fig3 = plt.figure(figsize=(12,5));   # smallish figure\n",
    "ax1 = plt.subplot(1, 1, 1);\n",
    "plt.axis( [ 0, tmax, 0, tmax] );\n",
    "mycm = plt.get_cmap('jet');\n",
    "plt.axis( [0, fmax/20, 0, 1.1*asdmax ] );\n",
    "plt.plot(f,asd,'k-', lw=2);\n",
    "for i in range(1,100):\n",
    "    plt.plot( gda_cvec(i*f0,i*f0), gda_cvec(0.0,asdmax/5.0), 'r-' );\n",
    "plt.xlabel('f (Kz)');\n",
    "plt.ylabel('asd');\n",
    "plt.show();\n",
    "print(\"Caption: ASD of pressure field for observer at %.2f\" % (x[ipx,0]) );\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
