{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n",
    "%reset -f\n",
    "\n",
    "# import various packages\n",
    "\n",
    "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   # 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",
    "from scipy import sparse\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 = cm.get_cmap('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",
    "# bandpass filter, used in seismological example, but hand\n",
    "# in a variety of settings involving time series\n",
    "def gda_chebyshevfilt(d, Dt, flow, fhigh):\n",
    "    # (dout,u,v)=gda_chebyshevfilt(d, Dt, flow, fhigh);\n",
    "    # chebyshev IIR bandpass filter\n",
    "    # d - input array of data\n",
    "    # Dt - sampling interval\n",
    "    # flow - low pass frequency, Hz\n",
    "    # fhigh - high pass frequency, Hz\n",
    "    # dout - output array of data\n",
    "    # u - the numerator filter\n",
    "    # v - the denominator filter\n",
    "    # these filters can be used again using dout=filter(u,v,din);\n",
    "\n",
    "    # make sure input timeseries is a column vector\n",
    "    s = np.shape(d);\n",
    "    N = s[0];\n",
    "    if(N==1):\n",
    "        dd = np.zeros((N,1));\n",
    "        dd[:,0] = d;\n",
    "    else:\n",
    "        dd=d;\n",
    "        \n",
    "    # sampling rate\n",
    "    rate=1/Dt;\n",
    "\n",
    "    # ripple parameter, set to ten percent\n",
    "    ripple=0.1;  \n",
    "\n",
    "    # normalise frequency\n",
    "    fl=2.0*flow/rate;\n",
    "    fh=2.0*fhigh/rate;\n",
    "\n",
    "    # center frequency \n",
    "    cf = 4 * tan( (fl*pi/2) ) * tan( (fh*pi/2) );\n",
    "\n",
    "    # bandwidth\n",
    "    bw = 2 * ( tan( (fh*pi/2) ) - tan( (fl*pi/2) ) );\n",
    "\n",
    "    # ripple parameter factor\n",
    "    rpf = sqrt((sqrt((1.0+1.0/(ripple*ripple))) + 1.0/ripple));\n",
    "    a = 0.5*(rpf-1.0/rpf);\n",
    "    b = 0.5*(rpf+1.0/rpf);\n",
    "\n",
    "    u=np.zeros((5,1));\n",
    "    v=np.zeros((5,1));\n",
    "    theta = 3*pi/4;\n",
    "    sr = a * cos(theta);\n",
    "    si = b * sin(theta);\n",
    "    es = sqrt(sr*sr+si*si);\n",
    "    tmp= 16.0 - 16.0*bw*sr + 8.0*cf + 4.0*es*es*bw*bw - 4.0*bw*cf*sr + cf*cf;\n",
    "    v[0,0] = 1.0;\n",
    "    v[1,0] = 4.0*(-16.0 + 8.0*bw*sr - 2.0*bw*cf*sr + cf*cf)/tmp;\n",
    "    v[2,0] = (96.0 - 16.0*cf - 8.0*es*es*bw*bw + 6.0*cf*cf)/tmp;\n",
    "    v[3,0] = (-64.0 - 32.0*bw*sr + 8.0*bw*cf*sr + 4.0*cf*cf)/tmp;\n",
    "    v[4,0] = (16.0 + 16.0*bw*sr + 8.0*cf + 4.0*es*es*bw*bw + 4.0*bw*cf*sr + cf*cf)/tmp;\n",
    "    tmp = 4.0*es*es*bw*bw/tmp;\n",
    "    u[0,0] = tmp;\n",
    "    u[1,0] = 0.0;\n",
    "    u[2,0] = -2.0*tmp;\n",
    "    u[3,0] = 0.0;\n",
    "    u[4,0] = tmp;\n",
    "\n",
    "    dout = sg.lfilter(u.ravel(),v.ravel(),dd.ravel());\n",
    "    return (gda_cvec(dout),gda_cvec(u),gda_cvec(v));\n",
    "\n",
    "# function to perform the multiplication FT (F v)\n",
    "def gdaFTFmul(v):\n",
    "    # this function is used by bicongugate gradient to\n",
    "    # solve the generalized least squares problem Fm=F.\n",
    "    # Note that \"F\" must be called \"gdaFsparse\".\n",
    "    global gdaFsparse;\n",
    "    s = np.shape(v);\n",
    "    if(len(s)==1):\n",
    "        vv = np.zeros((s[0],1));\n",
    "        vv[:,0] = v;\n",
    "    else:\n",
    "        vv=v;\n",
    "    # weirdly, \"*\" multiplies sparse matrices just fine\n",
    "    temp = gdaFsparse*vv;\n",
    "    return gdaFsparse.transpose()*temp;\n",
    "\n",
    "# function to perform the multiplication FT (G v)\n",
    "def gda_FTFmul(v):\n",
    "    # this function is used by bicongugate gradient to\n",
    "    # solve the generalized least squares problem Fm=F.\n",
    "    # Note that \"F\" must be called \"gdaFsparse\".\n",
    "    global gdaFsparse;\n",
    "    s = np.shape(v);\n",
    "    if(len(s)==1):\n",
    "        vv = np.zeros((s[0],1));\n",
    "        vv[:,0] = v;\n",
    "    else:\n",
    "        vv=v;\n",
    "    # weirdly, \"*\" multiplies sparse matrices just fine\n",
    "    temp = gdaFsparse*vv;\n",
    "    return gdaFsparse.transpose()*temp;\n",
    "\n",
    "\n",
    "# function to perform the multiplication FT f\n",
    "def gda_FTFrhs(f):\n",
    "    # this function is used by bicongugate gradient to\n",
    "    # solve the generalized least squares problem Fm=F.\n",
    "    # Note that \"F\" must be called \"gdaFsparse\".\n",
    "    global gdaFsparse;\n",
    "    return gdaFsparse.transpose()*f;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "model parameters:\n",
      "exterior bounrary 49 49 41 41\n",
      "volume            1911\n",
      "total  2091  expected  2091\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 13324 characters]",
      "text/plain": [
       "<Figure size 500x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Caption: Geometry, with each part a separate color.\n",
      "elements:  9735  expected  9735\n",
      "model parameters (recheck):\n",
      "exterior bounrary 49 49 41 41\n",
      "volume            1911\n",
      "total  2091  expected  2091\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 18952 characters]",
      "text/plain": [
       "<Figure size 500x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Caption: Temperature.\n"
     ]
    }
   ],
   "source": [
    "# Laplace's equation d2 T / dx2 + d2 T / dy2 = -f(x,y)\n",
    "# in a rectangular medium containing a rectangular hole\n",
    "# with T specified on exterior and interior surfaces\n",
    "\n",
    "# x-axis\n",
    "Lx = 41;\n",
    "xmin = 0.0;\n",
    "xmax = 1.0;\n",
    "Dx = (xmax-xmin)/(Lx-1);\n",
    "x = gda_cvec(Dx*np.linspace(0,Lx-1,Lx));\n",
    "Dxr = 1.0/Dx;\n",
    "Dxr2 = 1.0/(Dx**2);\n",
    "\n",
    "# y-axis\n",
    "Ly = 51;\n",
    "ymin=0.0;\n",
    "ymax=1.0;\n",
    "Dy = (ymax-ymin)/(Ly-1);\n",
    "y = gda_cvec(Dy*np.linspace(0,Ly-1,Ly));\n",
    "Dyr = 1.0/Dy;\n",
    "Dyr2 = 1.0/(Dy**2);\n",
    "\n",
    "# x location of the hole (the bold begins at ixA+1 and ends at ixB-1)\n",
    "ixA = 10;\n",
    "ixB = 20;\n",
    "\n",
    "# y location of the hole (the bold begins at iyA+1 and ends at iyB-1)\n",
    "iyA = 15;\n",
    "iyB = 32;\n",
    "\n",
    "# model parameters (unknowns)\n",
    "M = Lx*Ly;\n",
    "\n",
    "# lookup tables translate between (i,j) image an linear k ordering\n",
    "iofk = -np.ones((M,1),dtype=int); # x position\n",
    "jofk = -np.ones((M,1),dtype=int); # y position\n",
    "kofij = -np.ones((Lx,Ly),dtype=int); # -1 signifies in hole\n",
    "k=0;\n",
    "ks=0;\n",
    "for ix in range(Lx):\n",
    "    for iy in range(Ly):\n",
    "        iofk[k,0]=ix;\n",
    "        jofk[k,0]=iy;\n",
    "        kofij[ix,iy]=k;\n",
    "        k=k+1;\n",
    "\n",
    "# configure geometry, rectangle with hole\n",
    "geo = np.zeros((Lx,Ly));\n",
    "mel=0; mer=0; meb=0; met=0; mi=0;\n",
    "for ix in range(Lx):\n",
    "    for iy in range(Ly):\n",
    "        if( (ix==0) and (iy>0) and (iy<(Ly-1))): \n",
    "            geo[ix,iy]=5.0; # exterior left\n",
    "            mel=mel+1;\n",
    "        elif( (ix==(Lx-1)) and (iy>0) and (iy<(Ly-1)) ):\n",
    "            geo[ix,iy]=6.0; # exterior right\n",
    "            mer=mer+1;\n",
    "        elif( (iy==0)  ):\n",
    "            geo[ix,iy]=7.0;  # exterior bottom\n",
    "            meb=meb+1;\n",
    "        elif( (iy==(Ly-1)) ):\n",
    "            geo[ix,iy]=8.0;  # exterior top\n",
    "            met=met+1;\n",
    "        else:\n",
    "            geo[ix,iy]=9.0;  # interior\n",
    "            mi = mi+1;\n",
    "\n",
    "print('model parameters:')\n",
    "print('exterior bounrary', mel, mer, meb, met);\n",
    "print('volume           ', mi);\n",
    "print('total ', mi+mel+mer+meb+met, ' expected ', M);\n",
    "\n",
    "# boundary conditions (in the future, could make these variable)\n",
    "ExTop = np.zeros((met,1));\n",
    "ExBottom = np.zeros((meb,1));\n",
    "ExLeft = np.zeros((mel,1));\n",
    "ExRight = np.zeros((mer,1));\n",
    "\n",
    "# source\n",
    "f0 = 1e5;\n",
    "f = np.zeros((Lx,Ly));\n",
    "i0 = 12;\n",
    "j0 = 21;\n",
    "f[i0,j0] = f0*Dx*Dy;\n",
    "         \n",
    "fig1 = plt.figure(1,figsize=(5,4));\n",
    "IMG = geo;\n",
    "ax1 = plt.subplot(1, 1, 1);\n",
    "plt.axis([xmin,xmax,ymin,ymax]);\n",
    "mycmap=cm.jet;\n",
    "left=0; right=xmax;\n",
    "top=ymax; bottom=0;\n",
    "# for testing orientation of image\n",
    "if( 0 ):\n",
    "    IMG = np.zeros((Lx,Ly));\n",
    "    ix=2; iy=2; IMG[ix,iy]=1.0;\n",
    "    ix=Lx-2; iy=2; IMG[ix,iy]=2.0;\n",
    "    ix=2; iy=Ly-2; IMG[ix,iy]=3.0;\n",
    "dmin = np.min(IMG);\n",
    "dmax = np.max(IMG);\n",
    "plt.imshow(np.flipud(IMG.T), cmap=mycmap, vmin=dmin, vmax=dmax, extent=(left,right,bottom,top) );\n",
    "plt.colorbar(shrink=0.5);\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('y');\n",
    "plt.show()\n",
    "print(\"Caption: Geometry, with each part a separate color.\");       \n",
    "         \n",
    "Pmax = 100000;\n",
    "Dr = np.zeros((Pmax,1),dtype=int);\n",
    "Dc = np.zeros((Pmax,1),dtype=int);\n",
    "Dv = np.zeros((Pmax,1));\n",
    "s = np.zeros((M,1));\n",
    "\n",
    "nr=0;\n",
    "ne=0;\n",
    "mel=0; mer=0; meb=0; met=0; mi=0;\n",
    "for ix in range(Lx):\n",
    "    for iy in range(Ly):\n",
    "        if( (ix==0) and (iy>0) and (iy<(Ly-1))):\n",
    "            k=kofij[ix,iy]; # left exterior boundary\n",
    "            mel=mel+1;\n",
    "            Dr[ne,0] = nr;\n",
    "            Dc[ne,0] = k;\n",
    "            Dv[ne,0] = 1.0;\n",
    "            ne=ne+1;\n",
    "            s[k,0]=ExLeft[iy-1,0];\n",
    "            nr=nr+1;\n",
    "        elif( (ix==(Lx-1)) and (iy>0) and (iy<(Ly-1)) ):\n",
    "            k=kofij[ix,iy]; # right exterior boundary\n",
    "            mer=mer+1;\n",
    "            Dr[ne,0] = nr;\n",
    "            Dc[ne,0] = k;\n",
    "            Dv[ne,0] = 1.0;\n",
    "            ne=ne+1;\n",
    "            s[k,0]=ExRight[iy-1,0];\n",
    "            nr=nr+1;\n",
    "        elif( (iy==0) ):\n",
    "            k=kofij[ix,iy]; # bottom exterior boundary\n",
    "            meb=meb+1;\n",
    "            Dr[ne,0] = nr;\n",
    "            Dc[ne,0] = k;\n",
    "            Dv[ne,0] = 1.0;\n",
    "            ne=ne+1;\n",
    "            s[k,0]=ExBottom[ix,0];\n",
    "            nr=nr+1;\n",
    "        elif( (iy==(Ly-1)) ):\n",
    "            k=kofij[ix,iy]; # top exterior boundary\n",
    "            met=met+1;\n",
    "            Dr[ne,0] = nr;\n",
    "            Dc[ne,0] = k;\n",
    "            Dv[ne,0] = 1.0;\n",
    "            ne=ne+1;\n",
    "            s[k,0]=ExTop[ix,0];\n",
    "            nr=nr+1;\n",
    "        else:\n",
    "            k00=kofij[ix-1,iy-1]; # interior\n",
    "            k01=kofij[ix-1,iy];\n",
    "            k02=kofij[ix-1,iy+1];\n",
    "            k10=kofij[ix,iy-1];\n",
    "            k11=kofij[ix,iy];\n",
    "            k12=kofij[ix,iy+1];\n",
    "            k20=kofij[ix+1,iy-1];\n",
    "            k21=kofij[ix+1,iy];\n",
    "            k22=kofij[ix+1,iy+1];\n",
    "            mi=mi+1;\n",
    "            # equation d2 T / dx2 + d2 T / dy2 = -f\n",
    "            # d2 T / dx2\n",
    "            Dr[ne,0] = nr;\n",
    "            Dc[ne,0] = k01;\n",
    "            Dv[ne,0] = Dxr2;\n",
    "            ne=ne+1;\n",
    "            Dr[ne,0] = nr;\n",
    "            Dc[ne,0] = k11;\n",
    "            Dv[ne,0] = -2.0*Dxr2 -2.0*Dyr2;\n",
    "            ne=ne+1;\n",
    "            Dr[ne,0] = nr;\n",
    "            Dc[ne,0] = k21;\n",
    "            Dv[ne,0] = Dxr2;\n",
    "            ne=ne+1;\n",
    "            # d2 T / dy2\n",
    "            Dr[ne,0] = nr;\n",
    "            Dc[ne,0] = k10;\n",
    "            Dv[ne,0] = Dyr2;\n",
    "            ne=ne+1;\n",
    "            # Dr[ne,0] = nr; Absorbed into earlie element\n",
    "            # Dc[ne,0] = k11;\n",
    "            # Dv[ne,0] = -2.0*Dyr2;\n",
    "            # ne[ne,0]=ne+1;\n",
    "            Dr[ne,0] = nr;\n",
    "            Dc[ne,0] = k12;\n",
    "            Dv[ne,0] = Dyr2;\n",
    "            ne=ne+1;\n",
    "            s[k11,0]=-f[ix,iy];\n",
    "            nr=nr+1;\n",
    "\n",
    "print('elements: ',ne,' expected ',mel+mer+meb+met+5*mi);\n",
    "print('model parameters (recheck):')\n",
    "print('exterior bounrary', mel, mer, meb, met);\n",
    "print('volume           ', mi);\n",
    "print('total ', mi+mel+mer+meb+met, ' expected ', M);\n",
    "            \n",
    "Dsparse=sparse.coo_matrix((Dv[0:ne,0:1].ravel(),\n",
    "        (Dr[0:ne,0:1].ravel(),Dc[0:ne,0:1].ravel())),shape=(M,M)).tocsr();\n",
    "del Dv; del Dr; del Dc;\n",
    "mest = gda_cvec(las.spsolve(Dsparse,s.ravel()));\n",
    "\n",
    "    \n",
    "T = np.zeros((Lx,Ly));\n",
    "for k in range(M):\n",
    "    i = iofk[k,0];\n",
    "    j = jofk[k,0];\n",
    "    T[i,j] = mest[k,0];\n",
    "\n",
    "fig2 = plt.figure(2,figsize=(5,4));\n",
    "IMG = T;\n",
    "ax1 = plt.subplot(1, 1, 1);\n",
    "plt.axis([0.0,xmax,0.0,ymax]);\n",
    "mycmap=cm.jet;\n",
    "left=0; right=xmax;\n",
    "top=ymax; bottom=0;\n",
    "# for testing orientation of image\n",
    "if( 0 ):\n",
    "    IMG = np.zeros((Lx,Ly));\n",
    "    ix=2; iy=2; IMG[ix,iy]=1.0;\n",
    "    ix=Lx-2; iy=2; IMG[ix,iy]=2.0;\n",
    "    ix=2; iy=Ly-2; IMG[ix,iy]=3.0;\n",
    "dmin = np.min(IMG);\n",
    "dmax = np.max(IMG);\n",
    "plt.imshow(np.flipud(IMG.T), cmap=mycmap, vmin=dmin, vmax=dmax, extent=(left,right,bottom,top) );\n",
    "plt.colorbar(shrink=0.5);\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('y');\n",
    "plt.show()\n",
    "print(\"Caption: Temperature.\"); \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
