{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# ces_01\n",
    "%reset -f\n",
    "\n",
    "from math import exp, pi, sin, cos, tan, atan, sqrt, floor, log10, nan, log   # math functions\n",
    "import numpy as np\n",
    "import scipy.linalg as la             # linear algebra functions\n",
    "from scipy.integrate import solve_ivp\n",
    "import matplotlib\n",
    "from matplotlib import pyplot as plt  # general plotting \n",
    "from matplotlib.colors import ListedColormap\n",
    "import time;\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": 47,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 25880 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Caption: Initial conditions, C0\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 21776 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Caption: Source, S (assumed time-indpendent)\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 20872 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Caption: C at time t = 0.0000\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 21760 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Caption: C at time t = 8.0000\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 22292 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Caption: C at time t = 16.0000\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 22500 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Caption: C at time t = 24.0000\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 23024 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Caption: C at time t = 32.0000\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 23224 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Caption: C at time t = 40.0000\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 23816 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Caption: C at time t = 48.0000\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 23832 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Caption: C at time t = 56.0000\n",
      " \n",
      "done\n"
     ]
    }
   ],
   "source": [
    "# Makes figures for lecture\n",
    "# not part of group assignments\n",
    "\n",
    "kappa = 1.0; # diffusivity\n",
    "PLOTSTEP = 8; # plot every PLOTSTEP time slice\n",
    "\n",
    "# standard set-up of position x-axis \n",
    "Nx=64;\n",
    "Dx=1.0;\n",
    "x=gda_cvec(Dx*np.linspace(0,Nx-1,Nx));\n",
    "xmin = 0.0;\n",
    "xmax = Dx*(Nx-1);\n",
    "ixmid = int(floor(Nx/2));\n",
    "xmid=x[ixmid,0];\n",
    "\n",
    "# kx-axis\n",
    "Nkx=int(Nx/2+1);\n",
    "kxmax=2.0*pi/(2.0*Dx);\n",
    "Dkx=kxmax/(Nx/2);\n",
    "kxp = gda_cvec( Dkx*np.linspace(0,Nkx-1,Nkx));\n",
    "kxn = -np.flipud(kxp[1:Nkx-1,0:1]);\n",
    "kx = gda_cvec(kxp,kxn);\n",
    "\n",
    "# standard set-up of position y and wavenumber ky-axis\n",
    "Ny=64;\n",
    "Dy=1.0;\n",
    "y=gda_cvec(Dy*np.linspace(0,Ny-1,Ny));\n",
    "ymax = Dy*(Ny-1);\n",
    "iymid = int(floor(Ny/2));\n",
    "ymid=y[iymid,0]\n",
    "\n",
    "# ky-axis\n",
    "Nky=int(Ny/2+1);\n",
    "kymax=2.0*pi/(2.0*Dy);\n",
    "Dky=kymax/(Ny/2);\n",
    "kyp = gda_cvec( Dky*np.linspace(0,Nky-1,Nky));\n",
    "kyn = -np.flipud(kyp[1:Nky-1,0:1]);\n",
    "ky = gda_cvec(kyp,kyn);\n",
    "\n",
    "# time axis\n",
    "Nt = 64;\n",
    "Dt = 1.0;   \n",
    "t = Dt*np.linspace(0.0, Nt-1, Ny);  # note is an (Nt,) array\n",
    "t0 = t[0];\n",
    "tmax = t[Nt-1];\n",
    "\n",
    "# order of variables in arrays: x, y, t\n",
    "C = np.zeros((Nx,Ny,Nt));\n",
    "Ct = np.zeros((Nx,Ny,Nt),dtype=complex);\n",
    "\n",
    "# initial condition\n",
    "C0 = np.zeros((Nx,Ny));\n",
    "Csx = 5.0;\n",
    "Csy = 10.0;\n",
    "# CA = 0.0;\n",
    "CA = 1.0;\n",
    "for ix in range(Nx):\n",
    "    for iy in range(Ny):\n",
    "        r2 = ((x[ix,0]-xmid)/Csx)**2 + ((y[iy,0]-ymid)/Csy)**2;\n",
    "        C0[ix,iy] = CA*exp(-0.5*r2);\n",
    "\n",
    "fig2 = plt.figure(figsize=(12,5));   # smallish figure\n",
    "ax1 = plt.subplot(1, 1, 1);\n",
    "if( 0 ): # check prientation of plot\n",
    "    C0=np.zeros((Nx,Ny));\n",
    "    ix=2; iy=2; C0[ix,iy]=1.0\n",
    "    ix=38; iy=2; C0[ix,iy]=2.0;\n",
    "    ix=2; iy=38; C0[ix,iy]=3.0; # to check orientation\n",
    "plt.axis([0,xmax,0,ymax]);\n",
    "mycmap=matplotlib.colormaps['jet'];\n",
    "dmax = np.max(np.abs(C0));\n",
    "left=0; right=xmax;\n",
    "top=ymax; bottom=0;\n",
    "plt.imshow(np.flipud(C0.T), cmap=mycmap, vmin=-dmax, vmax=dmax, extent=(left,right,bottom,top) );\n",
    "plt.colorbar();\n",
    "plt.xlabel(\"x\");\n",
    "plt.ylabel(\"y\");\n",
    "plt.show();\n",
    "print('Caption: Initial conditions, C0');\n",
    "print(\" \");\n",
    "\n",
    "# source (assumed time-independent)\n",
    "S = np.zeros((Nx,Ny));\n",
    "Ssx = 10.0;\n",
    "Ssy = 5.0;\n",
    "SA = 0.05;\n",
    "# SA = 0.0;\n",
    "for ix in range(Nx):\n",
    "    for iy in range(Ny):\n",
    "        r2 = ((x[ix,0]-xmid)/Ssx)**2 + ((y[iy,0]-ymid)/Ssy)**2;\n",
    "        S[ix,iy] = SA*exp(-0.5*r2);\n",
    "\n",
    "fig1 = plt.figure(figsize=(12,5));   # smallish figure\n",
    "ax1 = plt.subplot(1, 1, 1);\n",
    "if( 0 ): # check prientation of plot\n",
    "    C0=np.zeros((Nx,Ny));\n",
    "    ix=2; iy=2; S[ix,iy]=1.0\n",
    "    ix=38; iy=2; S[ix,iy]=2.0;\n",
    "    ix=2; iy=38; S[ix,iy]=3.0; # to check orientation\n",
    "plt.axis([0,xmax,0,ymax]);\n",
    "mycmap=matplotlib.colormaps['jet'];\n",
    "dmax = np.max(np.abs(S));\n",
    "left=0; right=xmax;\n",
    "top=ymax; bottom=0;\n",
    "plt.imshow(np.flipud(S.T), cmap=mycmap, vmin=-dmax, vmax=dmax, extent=(left,right,bottom,top) );\n",
    "plt.colorbar();\n",
    "plt.xlabel(\"x\");\n",
    "plt.ylabel(\"y\");\n",
    "plt.show();\n",
    "print('Caption: Source, S (assumed time-indpendent)');\n",
    "print(\" \");\n",
    "\n",
    "# Fourier Transforms\n",
    "C0t = np.fft.fft2(C0);\n",
    "St = np.fft.fft2(S);\n",
    "\n",
    "# function for differential equation du/dt = f(u,t)\n",
    "# here u=C(kx,ky,t) and f=-k(kx**2+ky**2)+S\n",
    "def myfun( t, u ):\n",
    "    f = np.zeros((1,),dtype=complex);\n",
    "    f[0] =  -b * u[0] + St[ikx,iky];\n",
    "    return f;\n",
    "\n",
    "for ikx in range(Nx):\n",
    "    kx0 = kx[ikx,0];\n",
    "    for iky in range(Ny):\n",
    "        ky0 = ky[iky,0];\n",
    "        b = kappa*(kx0**2 + ky0**2);\n",
    "        u0 = np.zeros((1,),dtype=complex);\n",
    "        u0[0] = C0t[ikx,iky];\n",
    "        sol=solve_ivp(myfun, [t0, tmax], u0, method = 'RK45', t_eval=t);\n",
    "        Et = np.sum( np.power((t-sol.t)/tmax,2) );\n",
    "        Ct[ikx,iky,0:Nt]=sol.y[0];\n",
    "        \n",
    "for i in range(Nt):\n",
    "    C[0:Nx,0:Ny,i] = np.real(np.fft.ifft2( Ct[0:Nx,0:Ny,i] ));\n",
    "    \n",
    "dmax = np.max(np.abs(C));    \n",
    "for i in range(0,Nt,PLOTSTEP):\n",
    "    Ci = np.squeeze( C[0:Nx,0:Ny,i] );\n",
    "    figi = plt.figure(figsize=(12,5));   # smallish figure\n",
    "    ax1 = plt.subplot(1, 1, 1);\n",
    "    if( 0 ): # check prientation of plot\n",
    "        C0=np.zeros((Nx,Ny));\n",
    "        ix=2; iy=2; Ci[ix,iy]=1.0\n",
    "        ix=38; iy=2; Ci[ix,iy]=2.0;\n",
    "        ix=2; iy=38; Ci[ix,iy]=3.0; # to check orientation\n",
    "    plt.axis([0,xmax,0,ymax]);\n",
    "    mycmap=matplotlib.colormaps['jet'];\n",
    "    left=0; right=xmax;\n",
    "    top=ymax; bottom=0;\n",
    "    plt.imshow(np.flipud(Ci.T), cmap=mycmap, vmin=-dmax, vmax=dmax, extent=(left,right,bottom,top) );\n",
    "    plt.colorbar();\n",
    "    plt.xlabel(\"x\");\n",
    "    plt.ylabel(\"y\");\n",
    "    plt.show();\n",
    "    print(\"Caption: C at time t = %.4f\" % (t[i]) );\n",
    "    print(\" \");\n",
    "    \n",
    "print('done');\n",
    "\n",
    "      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 40264 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "saved file: d.txt\n"
     ]
    }
   ],
   "source": [
    "# Part 1\n",
    "# make uncorrelated Normally-distrumeted random noise\n",
    "# and save in file\n",
    "\n",
    "dfilename='d.txt';\n",
    "\n",
    "# standard set-up of position x-axis \n",
    "Nx=64;\n",
    "Dx=1.0;\n",
    "x=gda_cvec(Dx*np.linspace(0,Nx-1,Nx));\n",
    "xmin = 0.0;\n",
    "xmax = Dx*(Nx-1);\n",
    "ixmid = int(floor(Nx/2));\n",
    "xmid=x[ixmid,0];\n",
    "\n",
    "# standard set-up of position y and wavenumber ky-axis\n",
    "Ny=64;\n",
    "Dy=1.0;\n",
    "y=gda_cvec(Dy*np.linspace(0,Ny-1,Ny));\n",
    "ymin=0.0;\n",
    "ymax = Dy*(Ny-1);\n",
    "iymid = int(floor(Ny/2));\n",
    "ymid=y[iymid,0]\n",
    "\n",
    "d = np.random.normal(loc=0.0,scale=1.0,size=(Nx,Ny));     # random noise\n",
    "\n",
    "if( 0 ): # check prientation of plot\n",
    "    d=np.zeros((Nx,Ny));\n",
    "    ix=2; iy=2; d[ix,iy]=1.0\n",
    "    ix=Nx-5; iy=2; d[ix,iy]=2.0;\n",
    "    ix=2; iy=Ny-5; d[ix,iy]=3.0; # to check orientation\n",
    "\n",
    "fig1 = plt.figure(1,figsize=(12,5));   # smallish figure\n",
    "ax1 = plt.subplot(1, 1, 1);\n",
    "plt.axis([xmin,xmax,xmin,ymax]);\n",
    "mycmap=matplotlib.colormaps['jet'];\n",
    "dmax = np.max(np.abs(d));\n",
    "left=xmin; right=xmax;\n",
    "top=ymax; bottom=ymin;\n",
    "plt.imshow(np.flipud(d.T), cmap=mycmap, vmin=0, vmax=dmax, extent=(left,right,bottom,top) );\n",
    "plt.colorbar();\n",
    "plt.xlabel(\"x\");\n",
    "plt.ylabel(\"y\");\n",
    "plt.title(\"2D field of uncorrelated noise\");\n",
    "plt.show();\n",
    "\n",
    "# save matrix to file of tab-separated values\n",
    "np.savetxt(dfilename,d,delimiter=\"\\t\");\n",
    "print(\"saved file: %s\" % (dfilename) );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 7744 characters]",
      "text/plain": [
       "<Figure size 1600x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Caption: Exponential autocovariance, (left) target, (right) actual.\n",
      " \n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 52324 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "saved file: C0.txt\n"
     ]
    }
   ],
   "source": [
    "# Part2: Make C0.txt from d.txt\n",
    "dfilename='d.txt'\n",
    "C0filename='C0.txt';\n",
    "A = 1.0;  # 0 < C0 < A\n",
    "\n",
    "# TYPE=1 Gaussian\n",
    "# TYPE=2 Exponential\n",
    "# TYPE=3 Gaussian, Oscillatory\n",
    "# TYPE=4 Exponential, Oscillatory\n",
    "TYPE = 2;\n",
    "\n",
    "# EXACT=1: exactly satisfy c(x,y) on only approximately EXACT=0\n",
    "EXACT = 1;\n",
    "\n",
    "# field should have variance gamma2\n",
    "sigmad=1.0;\n",
    "sigmad2=sigmad**2;\n",
    " \n",
    "# field should have autocorrelation with scale factor, s\n",
    "sx = 0.1;   # scale factor\n",
    "sy = 0.1;   # scale factor\n",
    "sx2 = sx**2;\n",
    "sy2 = sy**2;\n",
    "\n",
    "# autocorrelation oscilation frequency (only some types)\n",
    "# Lx = wavelength in x, Ly = wavelength in y;\n",
    "Lx = 0.5;\n",
    "Ly = 0.5;\n",
    "kx0 = 2.0*pi/Lx;\n",
    "ky0 = 2.0*pi/Ly;\n",
    " \n",
    "# standard set-up of position x-axis \n",
    "Nx=64;\n",
    "xmin = -1.0;\n",
    "xmax = 1.0;\n",
    "x=gda_cvec(np.linspace(xmin,xmax,Nx));\n",
    "Dx = x[1,0]-x[0,0];\n",
    "ixmid = int(floor(Nx/2));\n",
    "xmid=x[ixmid,0];\n",
    "\n",
    "# kx-axis\n",
    "Nkx=int(Nx/2+1);\n",
    "kxmax=2.0*pi/(2.0*Dx);\n",
    "Dkx=kxmax/(Nx/2);\n",
    "kxp = gda_cvec( Dkx*np.linspace(0,Nkx-1,Nkx));\n",
    "kxn = -np.flipud(kxp[1:Nkx-1,0:1]);\n",
    "kx = gda_cvec(kxp,kxn);\n",
    "\n",
    "# standard set-up of position y and wavenumber ky-axis\n",
    "Ny=64;\n",
    "ymin = -1.0;\n",
    "ymax = 1.0;\n",
    "y=gda_cvec(np.linspace(ymin,ymax,Ny));\n",
    "Dy = y[1,0]-y[0,0]\n",
    "iymid = int(floor(Ny/2));\n",
    "ymid=y[iymid,0]\n",
    "\n",
    "# ky-axis\n",
    "Nky=int(Ny/2+1);\n",
    "kymax=2.0*pi/(2.0*Dy);\n",
    "Dky=kymax/(Ny/2);\n",
    "kyp = gda_cvec( Dky*np.linspace(0,Nky-1,Nky));\n",
    "kyn = -np.flipud(kyp[1:Nky-1,0:1]);\n",
    "ky = gda_cvec(kyp,kyn);\n",
    "\n",
    "c = np.zeros((Nx,Ny));\n",
    "if( TYPE==1 ):\n",
    "    covname = 'Gaussian';\n",
    "    for ix in range(Nx):\n",
    "        xi = x[ix,0];\n",
    "        for iy in range(Ny):\n",
    "            yi = y[iy,0];\n",
    "            r2 = (xi/sx)**2 + (yi/sy)**2;\n",
    "            c[ix,iy] = exp(-0.5*r2);\n",
    "elif( TYPE==2 ):\n",
    "    covname = 'Exponential';\n",
    "    for ix in range(Nx):\n",
    "        xi = x[ix,0];\n",
    "        for iy in range(Ny):\n",
    "            yi = y[iy,0];\n",
    "            r = abs(xi/sx) + abs(yi/sy);\n",
    "            c[ix,iy] = exp(-r);\n",
    "elif( TYPE==3 ):\n",
    "    covname = 'Gaussian-Oscillatory';\n",
    "    for ix in range(Nx):\n",
    "        xi = x[ix,0];\n",
    "        for iy in range(Ny):\n",
    "            yi = y[iy,0];\n",
    "            r2 = (xi/sx)**2 + (yi/sy)**2;\n",
    "            c[ix,iy] = exp(-0.5*r2)*cos(kx0*xi)*cos(ky0*yi);\n",
    "elif( TYPE==4 ):\n",
    "    covname = 'Exponential-Oscillatory';\n",
    "    for ix in range(Nx):\n",
    "        xi = x[ix,0];\n",
    "        for iy in range(Ny):\n",
    "            yi = y[iy,0];\n",
    "            r = abs(xi/sx) + abs(yi/sy);\n",
    "            c[ix,iy] = exp(-r)*cos(kx0*xi)*cos(ky0*yi);\n",
    "        \n",
    "# make m(x,y) with desired autocovariance, c using fact that\n",
    "# amplitude spectrum of autocorrelation is power spectrum of data\n",
    "\n",
    "casd = np.sqrt(np.abs(np.fft.fft2(c)));  # target amplitude spectral denisty\n",
    "d = np.genfromtxt(dfilename, delimiter=\"\\t\");\n",
    "dt = np.fft.fft2(d);\n",
    "dasd = np.abs(dt);\n",
    "\n",
    "if( EXACT ):\n",
    "    dasd = np.abs(dt);\n",
    "    dt = np.multiply(casd, np.divide(dt,dasd));\n",
    "else:\n",
    "    dt = np.multiply(casd, dt);\n",
    "    \n",
    "d = np.real(np.fft.ifft2(dt)); # shape spectrum on n into d\n",
    "sigmad_est = np.std(d);    # sqrt variance\n",
    "d = (sigmad/sigmad_est)*d; # normalize to correct variance\n",
    " \n",
    "cest = np.fft.ifft2(np.power(np.abs(np.fft.fft2(d)),2)); # autocorrelation of m\n",
    "cest = np.real(cest);\n",
    "\n",
    "# put into more reasonable order, with origin in the center\n",
    "cest2=np.zeros((Nx,Ny));\n",
    "cest2[Nx-Nkx:Nx,Ny-Nky:Ny] = cest[0:Nkx,0:Nky];\n",
    "cest2[0:Nkx-2,Ny-Nky:Ny] = cest[Nkx:Nx,0:Nky];\n",
    "cest2[0:Nkx-2,0:Nky-2] = cest[Nkx:Nx,Nky:Ny];\n",
    "cest2[Nkx-2:Nx,0:Nky-2] = cest[0:Nkx,Nky+0:Ny];\n",
    " \n",
    "gda_draw(\"title true\",c.T,\"title est\",cest2.T,);\n",
    "print(\"Caption: %s autocovariance, (left) target, (right) actual.\" % (covname) );\n",
    "print(' ');\n",
    "print(' ');\n",
    "\n",
    "# rescale to go from zero to A\n",
    "dmin=np.min(d);\n",
    "dmax=np.max(d);\n",
    "Dd = dmax-dmin;\n",
    "d = A*(d - dmin)/Dd;\n",
    "    \n",
    "if( 0 ): # check prientation of plot\n",
    "    d=np.zeros((Nx,Ny));\n",
    "    ix=2; iy=2; d[ix,iy]=1.0\n",
    "    ix=Nx-5; iy=2; d[ix,iy]=2.0;\n",
    "    ix=2; iy=Ny-5; d[ix,iy]=3.0; # to check orientation\n",
    "\n",
    "fig1 = plt.figure(1,figsize=(12,5));   # smallish figure\n",
    "ax1 = plt.subplot(1, 1, 1);\n",
    "plt.axis([xmin,xmax,xmin,ymax]);\n",
    "mycmap=matplotlib.colormaps['jet'];\n",
    "dmax = np.max(np.abs(d));\n",
    "left=xmin; right=xmax;\n",
    "top=ymax; bottom=ymin;\n",
    "plt.imshow(np.flipud(d.T), cmap=mycmap, vmin=0, vmax=dmax, extent=(left,right,bottom,top) );\n",
    "plt.colorbar();\n",
    "plt.xlabel(\"x\");\n",
    "plt.ylabel(\"y\");\n",
    "plt.title(\"2D field with %s autocovariance\" % (covname));\n",
    "plt.show();\n",
    "\n",
    "# save matrix to file of tab-separated values\n",
    "np.savetxt(C0filename,d,delimiter=\"\\t\");\n",
    "print(\"saved file: %s\" % (C0filename) );\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 7328 characters]",
      "text/plain": [
       "<Figure size 1600x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Caption: Gaussian autocovariance, (left) target, (right) actual.\n",
      " \n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 49368 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "saved file: S.txt\n"
     ]
    }
   ],
   "source": [
    "# Part 3: Make S.txt from d.txt\n",
    "\n",
    "dfilename='d.txt'\n",
    "Sfilename='S.txt';\n",
    "A = 0.05;  # 0 < S < A\n",
    "\n",
    "# TYPE=1 Gaussian\n",
    "# TYPE=2 Exponential\n",
    "# TYPE=3 Gaussian, Oscillatory\n",
    "# TYPE=4 Exponential, Oscillatory\n",
    "TYPE = 1;\n",
    "\n",
    "# EXACT=1: exactly satisfy c(x,y) on only approximately EXACT=0\n",
    "EXACT = 1;\n",
    "\n",
    "# field should have variance gamma2\n",
    "sigmad=1.0;\n",
    "sigmad2=sigmad**2;\n",
    " \n",
    "# field should have autocorrelation with scale factor, s\n",
    "sx = 0.1;   # scale factor\n",
    "sy = 0.1;   # scale factor\n",
    "sx2 = sx**2;\n",
    "sy2 = sy**2;\n",
    "\n",
    "# autocorrelation oscilation frequency (only some types)\n",
    "# Lx = wavelength in x, Ly = wavelength in y;\n",
    "Lx = 0.5;\n",
    "Ly = 0.5;\n",
    "kx0 = 2.0*pi/Lx;\n",
    "ky0 = 2.0*pi/Ly;\n",
    " \n",
    "# standard set-up of position x-axis \n",
    "Nx=64;\n",
    "xmin = -1.0;\n",
    "xmax = 1.0;\n",
    "x=gda_cvec(np.linspace(xmin,xmax,Nx));\n",
    "Dx = x[1,0]-x[0,0];\n",
    "ixmid = int(floor(Nx/2));\n",
    "xmid=x[ixmid,0];\n",
    "\n",
    "# kx-axis\n",
    "Nkx=int(Nx/2+1);\n",
    "kxmax=2.0*pi/(2.0*Dx);\n",
    "Dkx=kxmax/(Nx/2);\n",
    "kxp = gda_cvec( Dkx*np.linspace(0,Nkx-1,Nkx));\n",
    "kxn = -np.flipud(kxp[1:Nkx-1,0:1]);\n",
    "kx = gda_cvec(kxp,kxn);\n",
    "\n",
    "# standard set-up of position y and wavenumber ky-axis\n",
    "Ny=64;\n",
    "ymin = -1.0;\n",
    "ymax = 1.0;\n",
    "y=gda_cvec(np.linspace(ymin,ymax,Ny));\n",
    "Dy = y[1,0]-y[0,0]\n",
    "iymid = int(floor(Ny/2));\n",
    "ymid=y[iymid,0]\n",
    "\n",
    "# ky-axis\n",
    "Nky=int(Ny/2+1);\n",
    "kymax=2.0*pi/(2.0*Dy);\n",
    "Dky=kymax/(Ny/2);\n",
    "kyp = gda_cvec( Dky*np.linspace(0,Nky-1,Nky));\n",
    "kyn = -np.flipud(kyp[1:Nky-1,0:1]);\n",
    "ky = gda_cvec(kyp,kyn);\n",
    "\n",
    "c = np.zeros((Nx,Ny));\n",
    "if( TYPE==1 ):\n",
    "    covname = 'Gaussian';\n",
    "    for ix in range(Nx):\n",
    "        xi = x[ix,0];\n",
    "        for iy in range(Ny):\n",
    "            yi = y[iy,0];\n",
    "            r2 = (xi/sx)**2 + (yi/sy)**2;\n",
    "            c[ix,iy] = exp(-0.5*r2);\n",
    "elif( TYPE==2 ):\n",
    "    covname = 'Exponential';\n",
    "    for ix in range(Nx):\n",
    "        xi = x[ix,0];\n",
    "        for iy in range(Ny):\n",
    "            yi = y[iy,0];\n",
    "            r = abs(xi/sx) + abs(yi/sy);\n",
    "            c[ix,iy] = exp(-r);\n",
    "elif( TYPE==3 ):\n",
    "    covname = 'Gaussian-Oscillatory';\n",
    "    for ix in range(Nx):\n",
    "        xi = x[ix,0];\n",
    "        for iy in range(Ny):\n",
    "            yi = y[iy,0];\n",
    "            r2 = (xi/sx)**2 + (yi/sy)**2;\n",
    "            c[ix,iy] = exp(-0.5*r2)*cos(kx0*xi)*cos(ky0*yi);\n",
    "elif( TYPE==4 ):\n",
    "    covname = 'Exponential-Oscillatory';\n",
    "    for ix in range(Nx):\n",
    "        xi = x[ix,0];\n",
    "        for iy in range(Ny):\n",
    "            yi = y[iy,0];\n",
    "            r = abs(xi/sx) + abs(yi/sy);\n",
    "            c[ix,iy] = exp(-r)*cos(kx0*xi)*cos(ky0*yi);\n",
    "        \n",
    "# make m(x,y) with desired autocovariance, c using fact that\n",
    "# amplitude spectrum of autocorrelation is power spectrum of data\n",
    "\n",
    "casd = np.sqrt(np.abs(np.fft.fft2(c)));  # target amplitude spectral denisty\n",
    "d = np.genfromtxt(dfilename, delimiter=\"\\t\");\n",
    "dt = np.fft.fft2(d);\n",
    "dasd = np.abs(dt);\n",
    "\n",
    "if( EXACT ):\n",
    "    dasd = np.abs(dt);\n",
    "    dt = np.multiply(casd, np.divide(dt,dasd));\n",
    "else:\n",
    "    dt = np.multiply(casd, dt);\n",
    "    \n",
    "d = np.real(np.fft.ifft2(dt)); # shape spectrum on n into d\n",
    "sigmad_est = np.std(d);    # sqrt variance\n",
    "d = (sigmad/sigmad_est)*d; # normalize to correct variance\n",
    " \n",
    "cest = np.fft.ifft2(np.power(np.abs(np.fft.fft2(d)),2)); # autocorrelation of m\n",
    "cest = np.real(cest);\n",
    "\n",
    "# put into more reasonable order, with origin in the center\n",
    "cest2=np.zeros((Nx,Ny));\n",
    "cest2[Nx-Nkx:Nx,Ny-Nky:Ny] = cest[0:Nkx,0:Nky];\n",
    "cest2[0:Nkx-2,Ny-Nky:Ny] = cest[Nkx:Nx,0:Nky];\n",
    "cest2[0:Nkx-2,0:Nky-2] = cest[Nkx:Nx,Nky:Ny];\n",
    "cest2[Nkx-2:Nx,0:Nky-2] = cest[0:Nkx,Nky+0:Ny];\n",
    " \n",
    "gda_draw(\"title true\",c.T,\"title est\",cest2.T,);\n",
    "print(\"Caption: %s autocovariance, (left) target, (right) actual.\" % (covname) );\n",
    "print(' ');\n",
    "print(' ');\n",
    "\n",
    "# rescale to go from zero to A\n",
    "dmin=np.min(d);\n",
    "dmax=np.max(d);\n",
    "Dd = dmax-dmin;\n",
    "d = A*(d - dmin)/Dd;\n",
    "    \n",
    "if( 0 ): # check prientation of plot\n",
    "    d=np.zeros((Nx,Ny));\n",
    "    ix=2; iy=2; d[ix,iy]=1.0\n",
    "    ix=Nx-5; iy=2; d[ix,iy]=2.0;\n",
    "    ix=2; iy=Ny-5; d[ix,iy]=3.0; # to check orientation\n",
    "\n",
    "fig1 = plt.figure(1,figsize=(12,5));   # smallish figure\n",
    "ax1 = plt.subplot(1, 1, 1);\n",
    "plt.axis([xmin,xmax,xmin,ymax]);\n",
    "mycmap=matplotlib.colormaps['jet'];\n",
    "dmax = np.max(np.abs(d));\n",
    "left=xmin; right=xmax;\n",
    "top=ymax; bottom=ymin;\n",
    "plt.imshow(np.flipud(d.T), cmap=mycmap, vmin=0.0, vmax=dmax, extent=(left,right,bottom,top) );\n",
    "plt.colorbar();\n",
    "plt.xlabel(\"x\");\n",
    "plt.ylabel(\"y\");\n",
    "plt.title(\"2D field with %s autocovariance\" % (covname));\n",
    "plt.show();\n",
    "\n",
    "# save matrix to file of tab-separated values\n",
    "np.savetxt(Sfilename,d,delimiter=\"\\t\");\n",
    "print(\"saved file: %s\" % (Sfilename) );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 38900 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Caption: Initial conditions, C0\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 36716 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Caption: Source, S (assumed time-indpendent)\n",
      " \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 32696 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 25808 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 25288 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 25152 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 25356 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 24900 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 24616 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 24016 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 23696 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 23552 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 23196 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 23360 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 23332 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 23304 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 23372 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 23380 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 23404 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 23424 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 23536 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 23692 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 23352 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 23496 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 23496 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 24036 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 24444 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 24768 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 24952 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 25112 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 25016 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 24932 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 24436 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 24392 characters]",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "done\n"
     ]
    }
   ],
   "source": [
    "# PART 4: Solve equation to obtain C(t)\n",
    "\n",
    "C0filename = 'C0.txt';\n",
    "Sfilename = 'S.txt';\n",
    "\n",
    "kappa = 1.0; # diffusivity\n",
    "PLOTSTEP = 2; # plot every PLOTSTEP time slice\n",
    "\n",
    "# standard set-up of position x-axis \n",
    "Nx=64;\n",
    "Dx=1.0;\n",
    "x=gda_cvec(Dx*np.linspace(0,Nx-1,Nx));\n",
    "xmin = 0.0;\n",
    "xmax = Dx*(Nx-1);\n",
    "ixmid = int(floor(Nx/2));\n",
    "xmid=x[ixmid,0];\n",
    "\n",
    "# kx-axis\n",
    "Nkx=int(Nx/2+1);\n",
    "kxmax=2.0*pi/(2.0*Dx);\n",
    "Dkx=kxmax/(Nx/2);\n",
    "kxp = gda_cvec( Dkx*np.linspace(0,Nkx-1,Nkx));\n",
    "kxn = -np.flipud(kxp[1:Nkx-1,0:1]);\n",
    "kx = gda_cvec(kxp,kxn);\n",
    "\n",
    "# standard set-up of position y and wavenumber ky-axis\n",
    "Ny=64;\n",
    "Dy=1.0;\n",
    "y=gda_cvec(Dy*np.linspace(0,Ny-1,Ny));\n",
    "ymax = Dy*(Ny-1);\n",
    "iymid = int(floor(Ny/2));\n",
    "ymid=y[iymid,0]\n",
    "\n",
    "# ky-axis\n",
    "Nky=int(Ny/2+1);\n",
    "kymax=2.0*pi/(2.0*Dy);\n",
    "Dky=kymax/(Ny/2);\n",
    "kyp = gda_cvec( Dky*np.linspace(0,Nky-1,Nky));\n",
    "kyn = -np.flipud(kyp[1:Nky-1,0:1]);\n",
    "ky = gda_cvec(kyp,kyn);\n",
    "\n",
    "# time axis\n",
    "Nt = 64;\n",
    "Dt = 1.0;   \n",
    "t = Dt*np.linspace(0.0, Nt-1, Ny);  # note is an (Nt,) array\n",
    "t0 = t[0];\n",
    "tmax = t[Nt-1];\n",
    "\n",
    "# order of variables in arrays: x, y, t\n",
    "C = np.zeros((Nx,Ny,Nt));\n",
    "Ct = np.zeros((Nx,Ny,Nt),dtype=complex);\n",
    "\n",
    "# initial condition\n",
    "# read into C0\n",
    "C0 = np.genfromtxt(C0filename, delimiter=\"\\t\");\n",
    "\n",
    "fig1 = plt.figure(figsize=(12,5));   # smallish figure\n",
    "ax1 = plt.subplot(1, 1, 1);\n",
    "if( 0 ): # check prientation of plot\n",
    "    C0=np.zeros((Nx,Ny));\n",
    "    ix=2; iy=2; C0[ix,iy]=1.0\n",
    "    ix=38; iy=2; C0[ix,iy]=2.0;\n",
    "    ix=2; iy=38; C0[ix,iy]=3.0; # to check orientation\n",
    "plt.axis([0,xmax,0,ymax]);\n",
    "mycmap=matplotlib.colormaps['jet'];\n",
    "dmax = np.max(np.abs(C0));\n",
    "left=0; right=xmax;\n",
    "top=ymax; bottom=0;\n",
    "plt.imshow(np.flipud(C0.T), cmap=mycmap, vmin=0.0, vmax=dmax, extent=(left,right,bottom,top) );\n",
    "plt.colorbar();\n",
    "plt.xlabel(\"x\");\n",
    "plt.ylabel(\"y\");\n",
    "plt.show();\n",
    "print('Caption: Initial conditions, C0');\n",
    "print(\" \");\n",
    "\n",
    "# source (assumed time-independent)\n",
    "# read into S\n",
    "S = np.genfromtxt(Sfilename, delimiter=\"\\t\");\n",
    "\n",
    "fig2 = plt.figure(figsize=(12,5));   # smallish figure\n",
    "ax1 = plt.subplot(1, 1, 1);\n",
    "if( 0 ): # check prientation of plot\n",
    "    C0=np.zeros((Nx,Ny));\n",
    "    ix=2; iy=2; S[ix,iy]=1.0\n",
    "    ix=38; iy=2; S[ix,iy]=2.0;\n",
    "    ix=2; iy=38; S[ix,iy]=3.0; # to check orientation\n",
    "plt.axis([0,xmax,0,ymax]);\n",
    "mycmap=matplotlib.colormaps['jet'];\n",
    "dmax = np.max(np.abs(S));\n",
    "left=0; right=xmax;\n",
    "top=ymax; bottom=0;\n",
    "plt.imshow(np.flipud(S.T), cmap=mycmap, vmin=0.0, vmax=dmax, extent=(left,right,bottom,top) );\n",
    "plt.colorbar();\n",
    "plt.xlabel(\"x\");\n",
    "plt.ylabel(\"y\");\n",
    "plt.show();\n",
    "print('Caption: Source, S (assumed time-indpendent)');\n",
    "print(\" \");\n",
    "\n",
    "# Fourier Transforms\n",
    "C0t = np.fft.fft2(C0);\n",
    "St = np.fft.fft2(S);\n",
    "\n",
    "# function for differential equation du/dt = f(u,t)\n",
    "# here u=C(kx,ky,t) and f=-k(kx**2+ky**2)+S\n",
    "def myfun( t, u ):\n",
    "    f = np.zeros((1,),dtype=complex);\n",
    "    f[0] =  -b * u[0] + St[ikx,iky];\n",
    "    return f;\n",
    "\n",
    "for ikx in range(Nx):\n",
    "    kx0 = kx[ikx,0];\n",
    "    for iky in range(Ny):\n",
    "        ky0 = ky[iky,0];\n",
    "        b = kappa*(kx0**2 + ky0**2);\n",
    "        u0 = np.zeros((1,),dtype=complex);\n",
    "        u0[0] = C0t[ikx,iky];\n",
    "        sol=solve_ivp(myfun, [t0, tmax], u0, method = 'RK45', t_eval=t);\n",
    "        Et = np.sum( np.power((t-sol.t)/tmax,2) );\n",
    "        Ct[ikx,iky,0:Nt]=sol.y[0];\n",
    "        \n",
    "for i in range(Nt):\n",
    "    C[0:Nx,0:Ny,i] = np.real(np.fft.ifft2( Ct[0:Nx,0:Ny,i] ));\n",
    "    \n",
    "dmax = np.max(np.abs(C));    \n",
    "for i in range(0,Nt,PLOTSTEP):\n",
    "    Ci = np.squeeze( C[0:Nx,0:Ny,i] );\n",
    "    figi = plt.figure(figsize=(12,5));   # smallish figure\n",
    "    ax1 = plt.subplot(1, 1, 1);\n",
    "    if( 0 ): # check prientation of plot\n",
    "        C0=np.zeros((Nx,Ny));\n",
    "        ix=2; iy=2; Ci[ix,iy]=1.0\n",
    "        ix=38; iy=2; Ci[ix,iy]=2.0;\n",
    "        ix=2; iy=38; Ci[ix,iy]=3.0; # to check orientation\n",
    "    plt.axis([0,xmax,0,ymax]);\n",
    "    mycmap=matplotlib.colormaps['jet'];\n",
    "    left=0; right=xmax;\n",
    "    top=ymax; bottom=0;\n",
    "    plt.imshow(np.flipud(Ci.T), cmap=mycmap, vmin=-dmax, vmax=dmax, extent=(left,right,bottom,top) );\n",
    "    plt.colorbar();\n",
    "    plt.xlabel(\"x\");\n",
    "    plt.ylabel(\"y\");\n",
    "    plt.title(\"C at time t = %.4f\" % (t[i]) );\n",
    "    plt.show();\n",
    "\n",
    "print('done');"
   ]
  },
  {
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
