diff --git a/Lilly1966/.ipynb_checkpoints/OSE_solve-checkpoint.ipynb b/Lilly1966/.ipynb_checkpoints/OSE_solve-checkpoint.ipynb new file mode 100644 index 0000000..f8151d2 --- /dev/null +++ b/Lilly1966/.ipynb_checkpoints/OSE_solve-checkpoint.ipynb @@ -0,0 +1,704 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "d7f053a4", + "metadata": {}, + "source": [ + "# Orr-Sommerfeld 方程式固有値解析 (Chebyshev 多項式展開)\n", + "[Orszag (1971)](https://doi.org/10.1017/S0022112071002842) による Orr-Sommerfeld 方程式の Chebyshev 多項式展開による固有値解析\n", + "\n", + "ただし, 切断波数$N$が大きい場合, 4 階微分項の評価が $O(N^7)$となり, 固有値計算時の丸め誤差が無視できなくなる可能性がある ([Dongarra et al. 1996](https://doi.org/10.1016/S0168-9274(96)00049-9)). そのため, 中間変数$\\chi $を導入して 2 階微分のみで計算する. \n", + "\n", + "\\begin{equation}\n", + "\\dfrac{d^2\\upsilon}{dy^2} -\\alpha ^2\\upsilon \\equiv \\chi , \\label{eq:Orszag1971-AA5-1}\n", + "\\end{equation}\n", + "となる変数$\\chi $を導入すると, Orr-Sommerfeld 方程式［Orszag (1971) の (2) 式］は\n", + "\\begin{equation}\n", + "\\dfrac{d^2\\chi }{dy^2} -\\alpha ^2\\chi -i\\alpha R\\left[\\overline{u} \\chi -\\dfrac{d^2\\overline{u}}{dy^2} \\upsilon \\right] +i\\alpha R\\lambda \\chi =0. \\label{eq:Orszag1971-AA5-2}\n", + "\\end{equation}\n", + "$\\upsilon ,\\; \\chi ,\\; \\overline{u} $についてそれぞれ\n", + "\\begin{equation}\n", + "\\upsilon =\\sum^{\\infty}_{n=0}{a_{\\upsilon ,n}T_n}, \\quad \\chi =\\sum^{\\infty}_{n=0}{a_{\\chi ,n}T_n} ,\\quad \\overline{u} =\\sum^{\\infty}_{n=0}{b_nT_n} \\label{eq:Orszag1971-AA5-3}\n", + "\\end{equation}\n", + "という展開を行うと, 支配方程式は \n", + "$n$次チェビシェフ多項式$T_n$の$0\\leq n\\leq N-2$について, 以下の 2 式が得られる:\n", + "\\begin{align}\n", + "&\\; \\sum^{\\infty}_{\\substack{p=n+2\\\\ p\\equiv n\\; (\\mathrm{mod} \\; 2)}}{p(p^2-n^2)a_{\\upsilon ,p}} -c_n\\alpha ^2a_{\\upsilon ,n}-c_na_{\\chi ,n} =0, \\label{eq:Orszag1971-AA5-7} \\\\\n", + "&\\; \\sum^{\\infty}_{\\substack{p=n+2\\\\ p\\equiv n\\; (\\mathrm{mod} \\; 2)}}{p(p^2-n^2)a_{\\chi ,p}} -c_n\\alpha ^2a_{\\chi ,n} -\\dfrac{i\\alpha R}{2} \\left[\\sum^{\\infty}_{m=-\\infty}{\\overline{a}_{\\chi ,m}\\overline{b}_{n-m}} -\\sum^{\\infty}_{m=-\\infty}{\\overline{a}_{\\upsilon ,m}\\sum^{\\infty}_{\\substack{p=|n-m|+2\\\\ p\\equiv |n-m|\\; (\\mathrm{mod} \\; 2)}}{p(p^2-(n-m)^2)b_p}} \\right] =-i\\alpha R\\lambda c_na_{\\chi ,n} . \\label{eq:Orszag1971-AA5-8}\n", + "\\end{align}\n", + "境界条件［Orszag (1971) の (15)-(16) 式］から得られる係数の関係式:\n", + "$$\n", + "\\sum^{N_{\\mathrm{e}}}_{\\substack{p=0\\\\ p=0\\; (\\mathrm{mod}\\; 2)}}{a_{\\upsilon ,p}} =0,\\quad \n", + "\\sum^{N_{\\mathrm{e}}}_{\\substack{p=0\\\\ p=0\\; (\\mathrm{mod}\\; 2)}}{p^2a_{\\upsilon ,p}} =0,\\quad \n", + "\\sum^{N_{\\mathrm{o}}}_{\\substack{p=1\\\\ p=1\\; (\\mathrm{mod}\\; 2)}}{a_{\\upsilon ,p}} =0,\\quad \n", + "\\sum^{N_{\\mathrm{o}}}_{\\substack{p=1\\\\ p=1\\; (\\mathrm{mod}\\; 2)}}{p^2a_{\\upsilon ,p}} =0.\n", + "$$\n", + "上の 2 つの支配方程式および境界条件から, 求める固有値方程式は\n", + "$$A\\textbf{x} =\\lambda B\\textbf{x} ,$$\n", + "ただし, 区分行列$A_{1,2}$が支配方程式 1 において, 変数 2 の係数を表す$N-1$次正方行列とすると, \n", + "$$A\\equiv \\left(\\begin{array}{cc}\n", + "A_{1,\\upsilon} & 0 & 0 & A_{1,\\chi} & 0 & 0 \\\\\n", + "\\mathrm{BC1} & & & 0 & 0 & 0 \\\\\n", + "\\mathrm{BC2} & & & 0 & 0 & 0 \\\\\n", + "A_{2,\\upsilon} & 0 & 0 & A_{2,\\chi} & 0 & 0 \\\\\n", + "\\mathrm{BC3} & & & 0 & 0 & 0 \\\\\n", + "\\mathrm{BC4} & & & 0 & 0 & 0 \\\\\n", + "\\end{array} \\right) ,\\quad B\\equiv \\left(\\begin{array}{cc}\n", + "\\textbf{0} & 0 & 0 & \\textbf{0} & 0 & 0 \\\\\n", + "0 & & & 0 & 0 & 0 \\\\\n", + "0 & & & 0 & 0 & 0 \\\\\n", + "\\textbf{0} & 0 & 0 & B_{2,\\chi} & 0 & 0 \\\\\n", + "0 & & & 0 & 0 & 0 \\\\\n", + "0 & & & 0 & 0 & 0 \\\\\n", + "\\end{array} \\right) ,\\quad \\textbf{x} \\equiv \\left[a_{\\upsilon ,0},a_{\\upsilon ,1},\\cdots ,a_{\\upsilon ,N},a_{\\chi ,0},a_{\\chi ,1},\\cdots ,a_{\\chi ,N} \\right] ^T.$$" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "7e7d6392", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Nemax = 200, Nomax = 199\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): 0.02883820761178906 - 7578.242545455161im\n", + "Schur(Eig): 0.009405756684565517 - 1314.4376797941554im\n", + "Schur(Eig): 0.014244895720953492 - 537.8574465367766im\n", + "Schur(Eig): 0.016450157282747912 - 288.51348477771575im\n", + "Schur(Eig): 0.02103346120333471 - 180.39666028187585im\n", + "Schur(Eig): 0.025848848456281236 - 123.19195139204454im\n", + "Schur(Eig): 0.03205497200377114 - 89.63935804312915im\n", + "Schur(Eig): 0.038929374231420434 - 68.13201026223997im\n", + "Schur(Eig): 0.04694294410795094 - 53.61348178339329im\n", + "Schur(Eig): 0.05575204045304491 - 43.30296211449413im\n", + "Schur(Eig): 0.06559442227475659 - 35.75067193808814im\n", + "Schur(Eig): 0.07627540357947211 - 30.033554949848465im\n", + "Schur(Eig): 0.08792850011508294 - 25.61607108907338im\n", + "Schur(Eig): 0.10043070899048237 - 22.122849306497507im\n", + "Schur(Eig): 0.11386077327612609 - 19.320140971145047im\n", + "Schur(Eig): 0.1281349105002331 - 17.03250325944365im\n", + "Schur(Eig): 0.1432998698583759 - 15.145068407245297im\n", + "Schur(Eig): 0.15929548384870762 - 13.567075557511187im\n", + "Schur(Eig): 0.17614847010943072 - 12.236839362139095im\n", + "Schur(Eig): 0.19381421003017424 - 11.103631184207575im\n", + "Schur(Eig): 0.2123065069883877 - 10.13196095432081im\n", + "Schur(Eig): 0.23159223133775889 - 9.291696984526405im\n", + "Schur(Eig): 0.25167732522671404 - 8.56126472534193im\n", + "Schur(Eig): 0.27253885904938296 - 7.92186980103505im\n", + "Schur(Eig): 0.2941796917483501 - 7.359796586150784im\n", + "Schur(Eig): 0.3165886521483531 - 6.862846950645903im\n", + "Schur(Eig): 0.3397718683313401 - 6.421976647284318im\n", + "Schur(Eig): 0.36373573381017626 - 6.029010383316242im\n", + "Schur(Eig): 0.38850049383070856 - 5.6778068207963015im\n", + "Schur(Eig): 0.41410382288308695 - 5.3627488870422075im\n", + "Schur(Eig): 0.4406013358733947 - 5.079574494956819im\n", + "Schur(Eig): 0.46809249340979986 - 4.824359741766496im\n", + "Schur(Eig): 0.4967146295049181 - 4.594112202357917im\n", + "Schur(Eig): 0.5267049464373432 - 4.386087682420243im\n", + "Schur(Eig): 0.5584293662940205 - 4.1982417641691345im\n", + "Schur(Eig): 0.5926782370272254 - 4.028910498644979im\n", + "Schur(Eig): 0.6318873363430142 - 3.878163779238189im\n", + "Schur(Eig): 0.6774124487451393 - 3.763500053650113im\n", + "Schur(Eig): 0.6724188114836468 - 3.669503887191289im\n", + "Schur(Eig): 0.6666346567856843 - 3.5475978950393094im\n", + "Schur(Eig): 0.6667369112552747 - 3.429129755491025im\n", + "Schur(Eig): 0.6668393733825629 - 3.3134147800032405im\n", + "Schur(Eig): 0.666856311072032 - 3.199687877896379im\n", + "Schur(Eig): 0.6668698179540508 - 3.0879154146174312im\n", + "Schur(Eig): 0.6668850993839182 - 2.9781031043927317im\n", + "Schur(Eig): 0.6669018359263257 - 2.870250403683462im\n", + "Schur(Eig): 0.6669201648031444 - 2.7643562438989986im\n", + "Schur(Eig): 0.6669402701188141 - 2.6604194605606213im\n", + "Schur(Eig): 0.6669623585767975 - 2.5584387838874467im\n", + "Schur(Eig): 0.6669866659509853 - 2.4584128226652346im\n", + "Schur(Eig): 0.6670134605322281 - 2.3603400490888857im\n", + "Schur(Eig): 0.667043048541936 - 2.2642187830663874im\n", + "Schur(Eig): 0.6670757804235424 - 2.1700471755920914im\n", + "Schur(Eig): 0.6671120580366842 - 2.0778231866826116im\n", + "Schur(Eig): 0.6671523431785931 - 1.9875445633351365im\n", + "Schur(Eig): 0.6671971677312326 - 1.8992088120201795im\n", + "Schur(Eig): 0.6672471464171399 - 1.8128131702983779im\n", + "Schur(Eig): 0.6673029901943451 - 1.7283545703860157im\n", + "Schur(Eig): 0.667365524784893 - 1.6458296030286281im\n", + "Schur(Eig): 0.6674357106791179 - 1.5652344670759777im\n", + "Schur(Eig): 0.6675146678396625 - 1.486564926666811im\n", + "Schur(Eig): 0.23752648881857785 + 0.0037396706222145603im\n", + "Schur(Eig): 0.6676037066353783 - 1.409816245299302im\n", + "Schur(Eig): 0.6677043627182395 - 1.3349831229686595im\n", + "Schur(Eig): 0.19005924936950183 - 0.1828219254116369im\n", + "Schur(Eig): 0.667818440284551 - 1.2620596189975923im\n", + "Schur(Eig): 0.667948064841178 - 1.1910390667367665im\n", + "Schur(Eig): 0.6680957458771399 - 1.1219139721740363im\n", + "Schur(Eig): 0.9646425098876398 - 0.0351865839248333im\n", + "Schur(Eig): 0.34910682012256655 - 0.12450197754871109im\n", + "Schur(Eig): 0.6682644469478218 - 1.0546759064547755im\n", + "Schur(Eig): 0.6684576869017671 - 0.989315378367853im\n", + "Schur(Eig): 0.9363517811244114 - 0.06325156893621248im\n", + "Schur(Eig): 0.36849847834757754 - 0.23882483171782218im\n", + "Schur(Eig): 0.6686796236220521 - 0.9258216904581784im\n", + "Schur(Eig): 0.9080563344935584 - 0.09131286178046318im\n", + "Schur(Eig): 0.668935211379344 - 0.8641827952264582im\n", + "Schur(Eig): 0.4749011869127081 - 0.20873122005003544im\n", + "Schur(Eig): 0.8797556958156186 - 0.11937073101026466im\n", + "Schur(Eig): 0.6692302988188383 - 0.8043850926397907im\n", + "Schur(Eig): 0.8514493818764631 - 0.14742560074344907im\n", + "Schur(Eig): 0.5129162046302708 - 0.28662504241366654im\n", + "Schur(Eig): 0.8231369612514765 - 0.17547807357377554im\n", + "Schur(Eig): 0.7948183878980362 - 0.20352914421365847im\n", + "Schur(Eig): 0.5872129303355416 - 0.26716171254068716im\n", + "Schur(Eig): 0.7664940772247112 - 0.23158507171718073im\n", + "Schur(Eig): 0.6695718470141386 - 0.7464132968048778im\n", + "Schur(Eig): 0.6699680740685054 - 0.6902501370427058im\n", + "Schur(Eig): 0.7381150003547364 - 0.25969189546075594im\n", + "Schur(Eig): 0.6704286339966932 - 0.6358762450837643im\n", + "Schur(Eig): 0.7088747432491715 - 0.28765535594742253im\n", + "Schur(Eig): 0.6361048559307513 - 0.3251971610110621im\n", + "Schur(Eig): 0.6828633658739772 - 0.3076116006474047im\n", + "Schur(Eig): 0.6709649733030304 - 0.5832698557777029im\n", + "Schur(Eig): 0.6715904009498063 - 0.5324065147651961im\n", + "Schur(Eig): 0.6776436663762307 - 0.3437328375681226im\n", + "Schur(Eig): 0.6723226907333578 - 0.48325970758858783im\n", + "Schur(Eig): 0.6745109420659848 - 0.3898261303001523im\n", + "Schur(Eig): 0.6732077446999831 - 0.43579516755895803im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN" + ] + }, + { + "data": { + "text/plain": [ + "Figure(PyObject
)" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): 0.028851580470429898 - 7427.412242041225im\n", + "Schur(Eig): 0.009533588970052825 - 1288.2787233690856im\n", + "Schur(Eig): 0.014396958205195097 - 527.1566441913569im\n", + "Schur(Eig): 0.016628580243660488 - 282.7769998739129im\n", + "Schur(Eig): 0.02125292187361028 - 176.81237598404454im\n", + "Schur(Eig): 0.026119154949576815 - 120.74641332223527im\n", + "Schur(Eig): 0.032385535531883815 - 87.86174245066204im\n", + "Schur(Eig): 0.039329601038889495 - 66.7825557976176im\n", + "Schur(Eig): 0.04742195916392423 - 52.553069738534845im\n", + "Schur(Eig): 0.05631896904887011 - 42.44783633451269im\n", + "Schur(Eig): 0.06625806836511919 - 35.04592966338191im\n", + "Schur(Eig): 0.0770445496186914 - 29.442672481245136im\n", + "Schur(Eig): 0.08881157642882291 - 25.113181360813634im\n", + "Schur(Eig): 0.1014361042011465 - 21.689560561534627im\n", + "Schur(Eig): 0.11499649821234285 - 18.942711695841773im\n", + "Schur(Eig): 0.1294089198181688 - 16.700686600350803im\n", + "Schur(Eig): 0.14471971920455637 - 14.850901929562848im\n", + "Schur(Eig): 0.16086866736631758 - 13.304405276029732im\n", + "Schur(Eig): 0.17788206988671296 - 12.000737477058372im\n", + "Schur(Eig): 0.195715259165424 - 10.890181166870757im\n", + "Schur(Eig): 0.21438163512998343 - 9.937951376838607im\n", + "Schur(Eig): 0.23384807257382156 - 9.114517601005865im\n", + "Schur(Eig): 0.2541201752826265 - 8.39873379878657im\n", + "Schur(Eig): 0.2751751671241231 - 7.772181247029363im\n", + "Schur(Eig): 0.29701576942380037 - 7.221416953352625im\n", + "Schur(Eig): 0.3196313288656588 - 6.73448758325907im\n", + "Schur(Eig): 0.3430283427264351 - 6.302529352368175im\n", + "Schur(Eig): 0.3672145572259243 - 5.917532503421068im\n", + "Schur(Eig): 0.39221176530007096 - 5.57347953513074im\n", + "Schur(Eig): 0.4180608555520649 - 5.264869231416005im\n", + "Schur(Eig): 0.4448216679242136 - 4.9875274895461im\n", + "Schur(Eig): 0.4726011259137812 - 4.737614362022132im\n", + "Schur(Eig): 0.501547449209317 - 4.5122027838334775im\n", + "Schur(Eig): 0.531918178839164 - 4.308614522888535im\n", + "Schur(Eig): 0.5641228380998303 - 4.124868378717572im\n", + "Schur(Eig): 0.5991030757694588 - 3.9594161375064765im\n", + "Schur(Eig): 0.640132839257879 - 3.8130758046841033im\n", + "Schur(Eig): 0.6818692445843058 - 3.70965326550958im\n", + "Schur(Eig): 0.6699810259372195 - 3.609524425068888im\n", + "Schur(Eig): 0.666593083128247 - 3.4881139930364884im\n", + "Schur(Eig): 0.6668188394688124 - 3.3709146755875516im\n", + "Schur(Eig): 0.666893879320586 - 3.2562057758916922im\n", + "Schur(Eig): 0.6669098098466136 - 3.143458355138187im\n", + "Schur(Eig): 0.6669255853951719 - 3.032666881924751im\n", + "Schur(Eig): 0.6669431586289576 - 2.9238354325428593im\n", + "Schur(Eig): 0.66696231652267 - 2.81696315510865im\n", + "Schur(Eig): 0.6669832230117946 - 2.71204894108598im\n", + "Schur(Eig): 0.6670060760038876 - 2.6090915843046405im\n", + "Schur(Eig): 0.6670310976186595 - 2.508089766507401im\n", + "Schur(Eig): 0.6670585401646524 - 2.4090420414631817im\n", + "Schur(Eig): 0.6670886906235786 - 2.31194682097146im\n", + "Schur(Eig): 0.6671218762052807 - 2.216802357546517im\n", + "Schur(Eig): 0.6671584705378463 - 2.1236067257038025im\n", + "Schur(Eig): 0.6671989015460608 - 2.032357799685315im\n", + "Schur(Eig): 0.6672436603264309 - 1.9430532304452381im\n", + "Schur(Eig): 0.6672933118803402 - 1.8556904157206815im\n", + "Schur(Eig): 0.6673485079369096 - 1.7702664692964718im\n", + "Schur(Eig): 0.667410001879772 - 1.686778184264008im\n", + "Schur(Eig): 0.6674786675328377 - 1.6052219899119025im\n", + "Schur(Eig): 0.6675555199095552 - 1.5255939061026371im\n", + "Schur(Eig): 0.6676417426648599 - 1.4478894860899223im\n", + "Schur(Eig): 0.6677387180326144 - 1.3721037552372932im\n", + "Schur(Eig): 0.27720434380851733 - 0.05089872725658317im\n", + "Schur(Eig): 0.6678480672909648 - 1.2982311373576947im\n", + "Schur(Eig): 0.21272578235378325 - 0.19936069475434753im\n", + "Schur(Eig): 0.6679716919722942 - 1.2262653745565242im\n", + "Schur(Eig): 0.6681118358356443 - 1.1561994303389072im\n", + "Schur(Eig): 0.6682711424754932 - 1.0880253863300817im\n", + "Schur(Eig): 0.9646309154505646 - 0.03516727763104578im\n", + "Schur(Eig): 0.6684527468496686 - 1.0217343142581585im\n", + "Schur(Eig): 0.4163510155763265 - 0.1382265253000814im\n", + "Schur(Eig): 0.6686603562446971 - 0.9573161471892272im\n", + "Schur(Eig): 0.9363165358812776 - 0.06320149583980876im\n", + "Schur(Eig): 0.3839876109074957 - 0.26510649960903276im\n", + "Schur(Eig): 0.9079830546296593 - 0.09122273543365653im\n", + "Schur(Eig): 0.6688983789719299 - 0.8947595190105614im\n", + "Schur(Eig): 0.8796272922085633 - 0.11923285262208567im\n", + "Schur(Eig): 0.6691720559744425 - 0.834051590850208im\n", + "Schur(Eig): 0.5320452088016702 - 0.2064652191117368im\n", + "Schur(Eig): 0.8512458401120699 - 0.1472339290642039im\n", + "Schur(Eig): 0.6694875891020631 - 0.7751778696700894im\n", + "Schur(Eig): 0.8228350407533278 - 0.1752286787491913im\n", + "Schur(Eig): 0.5222281555545437 - 0.3142970585337267im\n", + "Schur(Eig): 0.7943883846246975 - 0.20322066451492837im\n", + "Schur(Eig): 0.6698523995623807 - 0.7181220051500524im\n", + "Schur(Eig): 0.7658768124412472 - 0.2311859878020367im\n", + "Schur(Eig): 0.6367193825465925 - 0.25988571005387817im\n", + "Schur(Eig): 0.7374157502565554 - 0.2587170852291249im\n", + "Schur(Eig): 0.6702751958863984 - 0.6628655095746344im\n", + "Schur(Eig): 0.7123159247889363 - 0.2855146556769195im\n", + "Schur(Eig): 0.6707663761310181 - 0.6093876651212601im\n", + "Schur(Eig): 0.6713381853349581 - 0.5576650860648068im\n", + "Schur(Eig): 0.6941568719197801 - 0.31949859450281015im\n", + "Schur(Eig): 0.6424176412961435 - 0.35254526646010687im\n", + "Schur(Eig): 0.6720075259566941 - 0.5076720092943916im\n", + "Schur(Eig): 0.6784240677887209 - 0.36421754924582495im\n", + "Schur(Eig): 0.6728223857272042 - 0.4593698713921976im\n", + "Schur(Eig): 0.6740837771202557 - 0.41247342802481496im\n", + "Schur(Eig): NaN + NaN*im\n" + ] + }, + { + "data": { + "text/plain": [ + "PyObject Text(0.5, 1.0, 'OSE spectra: Re = 10000.0, α = 1.0')" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using OffsetArrays # 配列要素番号を任意開始にするモジュール (行列計算時は 1 始まりの通常 matrix に copy する)\n", + "using LinearAlgebra\n", + "\n", + "#-- Drawing parameters \n", + "N = 200 # wavenumber in Chebyshev\n", + "Ncal = N + 1\n", + "Ntot = 2 * (N + 1)\n", + "\n", + "# initialize working arrays\n", + "# A: non-lambda coefficients, B: lambda coefficients\n", + "A_J = reshape(zeros(Complex{Float64},Ntot,Ntot),Ntot,Ntot) # Matrix for calculation\n", + "B_J = reshape(zeros(Complex{Float64},Ntot,Ntot),Ntot,Ntot) # Matrix for calculation\n", + "### --- From here, supposing matrix element as row, col as in mathematical notation\n", + "Acal = reshape(zeros(Complex{Float64},Ntot,Ntot),Ntot,Ntot) # Matrix for calculation\n", + "Bcal = reshape(zeros(Complex{Float64},Ntot,Ntot),Ntot,Ntot) # Matrix for calculation\n", + "A_vv = OffsetArray(Acal[1:Ncal,1:Ncal],0:N,0:N) # coefficient matrix for v in v Chebyshev\n", + "A_vchi = OffsetArray(Acal[1:Ncal,1:Ncal],0:N,0:N) # coefficient matrix for chi in v Chebyshev\n", + "A_chiv = OffsetArray(Acal[1:Ncal,1:Ncal],0:N,0:N) # coefficient matrix for chi in chi Chebyshev\n", + "A_chichi = OffsetArray(Acal[1:Ncal,1:Ncal],0:N,0:N) # coefficient matrix for chi in chi Chebyshev\n", + "B_chichi = OffsetArray(Acal[1:Ncal,1:Ncal],0:N,0:N) # coefficient matrix for chi in chi Chebyshev\n", + "### --- From here, supposing matrix element as col, row as in julia notation (only using these arrays for matrix calc)\n", + "b = OffsetArray(zeros(Float64,Ncal),0:N) # A vector composed of Chebyshev coefficients for the basic state flow\n", + "c = OffsetArray(zeros(Float64,Ncal),0:N) # factorized vector (c_0=2, c_n=1, (n>0))\n", + "\n", + "# setting parameters\n", + "α = 1.0\n", + "Re = 10000.0\n", + "\n", + "# setting the Chebyshev coefficients for the basic flow \n", + "b[0] = 0.5 # For the Poiseuille flow (u = 1-y^2)\n", + "b[2] = -b[0]\n", + "\n", + "# setting the factorized vector\n", + "c[0] = 2.0\n", + "c[1:N] .= 1.0\n", + "\n", + "# setting maximum even and odd numbers from N\n", + "Nemax = 0\n", + "Nomax = 1\n", + "for i in 0:2:N\n", + " Nemax = i\n", + "end\n", + "for i in 1:2:N\n", + " Nomax = i\n", + "end\n", + "println(\"Nemax = \", Nemax, \", Nomax = \", Nomax)\n", + "\n", + "# setting matrices\n", + "# 支配方程式から構築される部分 (チェビシェフ次数 0<=n<=N-2 の範囲)\n", + "# この段階では境界条件で成り立つ係数の関係は考慮しない\n", + "for n in 0:N-2\n", + " n2 = n^2\n", + " \n", + " # In v rows, a_{v,n} and a_{chi,n}\n", + " for p in n+2:2:N\n", + " p2 = p^2\n", + " A_vv[n,p] = p * (p2 - n2)\n", + " end\n", + " \n", + " A_vv[n,n] = - c[n] * α^2\n", + " A_vchi[n,n] = - c[n]\n", + " \n", + " # In chi rows, a_{v,n} and a_{chi,n}\n", + " for p in n+2:2:N\n", + " p2 = p^2\n", + " A_chichi[n,p] = p * (p2 - n2)\n", + " end\n", + " \n", + " A_chichi[n,n] = - c[n] * α^2\n", + "\n", + " for p in -N:N\n", + " abs_n_p = abs(n - p)\n", + " abs_p = abs(p)\n", + " if abs_n_p <= N\n", + " A_chichi[n,abs_p] = A_chichi[n,abs_p] + complex(0.0, - 0.5 * α * Re * c[abs_p] * c[abs_n_p] * b[abs_n_p])\n", + " end\n", + "\n", + " if abs_n_p+2 <= N\n", + " acoef = 0.0\n", + " for m in abs_n_p+2:2:N\n", + " m2 = m^2\n", + " acoef = acoef + m * (m2 - (n - p)^2) * b[m]\n", + " end\n", + " A_chiv[n,abs_p] = A_chiv[n,abs_p] + complex(0.0, 0.5 * α * Re * c[abs_p] * acoef)\n", + " end\n", + " end\n", + " \n", + " # In chi rows, lambda coefficients\n", + " B_chichi[n,n] = complex(0.0, - α * Re * c[n])\n", + "\n", + "end\n", + "\n", + "# Set boundary conditions\n", + "for p in 0:2:N # even mode\n", + " A_vv[N-1,p] = 1.0\n", + "# A_chichi[N-1,p] = 1.0\n", + " A_chiv[N-1,p] = p^2\n", + "end\n", + "for p in 1:2:N # even mode\n", + " A_vv[N,p] = 1.0\n", + "# A_chichi[N,p] = 1.0\n", + " A_chiv[N,p] = p^2\n", + "end\n", + "\n", + "# Gather A_v, A_chi, B_chi -> Acal, Bcal\n", + "Acal[1:Ncal,1:Ncal] .= A_vv[0:N,0:N]\n", + "Acal[1:Ncal,Ncal+1:Ntot] .= A_vchi[0:N,0:N]\n", + "Acal[Ncal+1:Ntot,1:Ncal] .= A_chiv[0:N,0:N]\n", + "Acal[Ncal+1:Ntot,Ncal+1:Ntot] .= A_chichi[0:N,0:N]\n", + "Bcal[Ncal+1:Ntot,Ncal+1:Ntot] .= B_chichi[0:N,0:N]\n", + "\n", + "# Convert Transpose type to Matrix (Array) type (ドット演算して要素だけコピーしている)\n", + "A_J .= transpose(Acal) # A_J = Transpose type\n", + "B_J .= transpose(Bcal) # B_J = Transpose type\n", + "\n", + "F = schur(A_J,B_J)\n", + "Dcal = F.α ./ F.β\n", + "for i in 1:Ntot\n", + " println(\"Schur(Eig): \", Dcal[i])\n", + "end\n", + "\n", + "\n", + "#-- Drawing\n", + "##########\n", + "# Plot #\n", + "##########\n", + "\n", + "using PyCall\n", + "using PyPlot\n", + "\n", + "rc(\"font\", family=\"IPAPGothic\")\n", + "fig = figure(\"pyplot_majorminor\",figsize=(7,5))\n", + "\n", + "x = real.(Dcal[1:Ntot])\n", + "y = imag.(Dcal[1:Ntot])\n", + "\n", + "#ax = gca()\n", + "\n", + "#cp = scatter(x,y)\n", + "scatter(x,y)\n", + "xlim(0.0,1.0)\n", + "ylim(-1.0,0.1)\n", + "\n", + "#ax = gca()\n", + "\n", + "xlabel(\"Real(λ)\")\n", + "ylabel(\"Imag(λ)\")\n", + "#if i == 1 #div(ne,2)-1\n", + "# plt.colorbar(cp)\n", + "#end\n", + "grid(\"on\")\n", + "\n", + "fig.canvas.draw()\n", + "gcf() # Needed for IJulia to plot inline\n", + "\n", + "ttext = \"OSE spectra: Re = $(Re), α = $(α)\"\n", + "title(ttext)\n", + "\n", + "#fname = \"test.png\"\n", + "#savefig(fname)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "026799fc", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.6.1", + "language": "julia", + "name": "julia-1.6" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.6.1" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Lilly1966/.ipynb_checkpoints/OSE_solve_Orszag1971-checkpoint.ipynb b/Lilly1966/.ipynb_checkpoints/OSE_solve_Orszag1971-checkpoint.ipynb new file mode 100644 index 0000000..1214239 --- /dev/null +++ b/Lilly1966/.ipynb_checkpoints/OSE_solve_Orszag1971-checkpoint.ipynb @@ -0,0 +1,602 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "d7f053a4", + "metadata": {}, + "source": [ + "# 基礎方程式 Orszag 1971 の (15)-(17) 式\n", + "$n$次チェビシェフ多項式$T_n$の$0\\leq n\\leq N-4$について, \n", + "\\begin{equation}\n", + "\\begin{split}\n", + "&\\dfrac{1}{24} \\sum^{N}_{\\substack{p=n+4\\\\ p\\equiv n\\; (\\mathrm{mod} \\; 2)}}{p\\left[p^2(p^2-4)^2-3n^2p^4+3n^4p^2-n^2\\left(n^2-4 \\right) ^2-48\\alpha ^2(p^2-n^2) \\right] a_p} -8\\alpha ^2(n+2)(n+1)a_{n+2} +\\alpha ^4c_na_n \\\\\n", + "&\\quad -\\dfrac{i\\alpha R}{2} \\left[\\sum^{N}_{p=2}{a_p\\sum^{}_{\\substack{m\\equiv p\\; (\\mathrm{mod} \\; 2) \\\\ |m|\\leq p-2 \\\\ |n-m|\\leq N}}{p(p^2-m^2)c_{|n-m|}b_{|n-m|} }} -\\alpha ^2\\sum^{}_{\\substack{|p|\\leq N \\\\ |n-p|\\leq N}}{c_{|p|}a_{|p|}c_{|n-p|}b_{|n-p|}} -\\sum^{}_{\\substack{|p|\\leq N \\\\ |n-p|\\leq N}}{c_{|p|}a_{|p|}\\sum^{N}_{\\substack{m=|n-p|+2\\\\ m+n\\equiv p\\; (\\mathrm{mod} \\; 2)}}{m\\left[m^2-(n-p)^2 \\right] b_m}} \\right] \\\\\n", + "&\\quad +\\lambda \\left\\{i\\alpha R\\sum^{N}_{\\substack{p=n+2\\\\ p\\equiv n\\; (\\mathrm{mod} \\; 2)}}{p(p^2-n^2)a_p} -i\\alpha ^3Rc_na_n \\right\\} =0.\n", + "\\end{split}\n", + " \\label{eq:Orszag1971-17-A44}\n", + "\\end{equation}\n", + "および境界条件から得られる係数の関係式:\n", + "\\begin{align}\n", + "a_{N_{\\mathrm{e}}}=&\\; \\sum^{N_{\\mathrm{e}}-4}_{\\substack{p=0\\\\ p=0\\; (\\mathrm{mod}\\; 2)}}{\\dfrac{(N_{\\mathrm{e}}-2)^2-p^2}{4(N_{\\mathrm{e}}-1)} a_p} ,\\; (\\mathrm{even mode}), \\label{eq:Orszag1971-17-A45} \\\\\n", + "a_{N_{\\mathrm{e}}-2}=&\\; \\sum^{N_{\\mathrm{e}}-4}_{\\substack{p=0\\\\ p=0\\; (\\mathrm{mod}\\; 2)}}{\\dfrac{p^2-N^2_{\\mathrm{e}}}{4(N_{\\mathrm{e}}-1)} a_p} ,\\; (\\mathrm{even mode}), \\label{eq:Orszag1971-17-A46} \\\\\n", + "a_{N_{\\mathrm{o}}}=&\\; \\sum^{N_{\\mathrm{o}}-4}_{\\substack{p=1\\\\ p=1\\; (\\mathrm{mod}\\; 2)}}{\\dfrac{(N_{\\mathrm{o}}-2)^2-p^2}{4(N_{\\mathrm{o}}-1)} a_p} ,\\; (\\mathrm{odd mode}), \\label{eq:Orszag1971-17-A47} \\\\\n", + "a_{N_{\\mathrm{o}}-2}=&\\; \\sum^{N_{\\mathrm{o}}-4}_{\\substack{p=1\\\\ p=1\\; (\\mathrm{mod}\\; 2)}}{\\dfrac{p^2-N^2_{\\mathrm{o}}}{4(N_{\\mathrm{o}}-1)} a_p} ,\\; (\\mathrm{odd mode}), \\label{eq:Orszag1971-17-A48} \n", + "\\end{align}\n", + "求める固有値方程式は\n", + "$$\n", + "(A-B\\lambda )\\textbf{x} =\\textbf{0} \\quad \\rightarrow \\quad \\left(B^{-1}A-\\lambda I \\right) \\textbf{x} =\\textbf{0},\n", + "$$\n", + "$$\n", + "\\textbf{x} \\equiv \\left[a_0,a_1,\\cdots ,a_N \\right] ^T.\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "7e7d6392", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Nemax = 56, Nomax = 55\n", + "Eig: 0.4999999999334932 + 903.0589774662305im\n", + "Eig: 0.4999999999988134 - 34.897806452343005im\n", + "Eig: 0.4999999999996687 - 0.06311310408499561im\n", + "Eig: 0.49999999999970246 - 6.56472786547905im\n", + "Eig: 0.49999999999977174 - 0.01192747960575369im\n", + "Eig: 0.4999999999998186 - 0.14206996446886708im\n", + "Eig: 0.4999999999998308 - 2.684892171972676im\n", + "Eig: 0.4999999999998447 - 0.49079903717010587im\n", + "Eig: 0.4999999999998891 - 0.029692874205562946im\n", + "Eig: 0.49999999999989275 - 0.9336504164573606im\n", + "Eig: 0.49999999999990624 - 0.002057057279983828im\n", + "Eig: 0.4999999999999148 - 0.019823234068627366im\n", + "Eig: 0.4999999999999202 - 0.08891054855459826im\n", + "Eig: 0.4999999999999247 - 0.2723463450930247im\n", + "Eig: 0.49999999999992484 - 0.2073607808223159im\n", + "Eig: 0.4999999999999384 - 1.4669377227922586im\n", + "Eig: 0.4999999999999414 - 0.6543836711906931im\n", + "Eig: 0.4999999999999503 - 0.055353878991148485im\n", + "Eig: 0.49999999999995354 - 0.2373475309261251im\n", + "Eig: 0.49999999999995787 - 0.1797111270289892im\n", + "Eig: 0.4999999999999678 - 0.04153641999564741im\n", + "Eig: 0.4999999999999788 - 0.1193698705370567im\n", + "Eig: 0.49999999999998934 - 0.008829999455944851im\n", + "Eig: 0.49999999999999484 - 0.10864976174862244im\n", + "Eig: 0.4999999999999986 - 0.1540499488689494im\n", + "Eig: 0.5000000000000091 - 0.3193282346934045im\n", + "Eig: 0.5000000000000147 - 0.07114525462434143im\n", + "Eig: 0.5000000000000318 - 0.38763707816857146im\n", + "Eig: 0.5000000000000565 - 0.13036289471223642im\n", + "Eig: 0.500000000000139 - 0.006005543666073486im\n", + "Eig: 0.5000000000001434 - 0.19339246974576957im\n", + "Eig: 0.5000000000004298 - 0.34141967303665405im\n", + "Eig: 0.5000000000005788 - 0.2532778737893288im\n", + "Eig: 0.500000000000636 - 0.0798914410897462im\n", + "Eig: 0.5000000000006692 - 0.2220456958587949im\n", + "Eig: 0.5000000000007411 - 0.035478176438828495im\n", + "Eig: 0.5000000000007552 - 0.2905620974815283im\n", + "Eig: 0.5000000000008235 - 0.41539491388624555im\n", + "Eig: 0.5000000000008713 - 0.024621575876883753im\n", + "Eig: 0.500000000000906 - 0.16674398210724284im\n", + "Eig: 0.5000000000009611 - 0.04830868351926353im\n", + "Eig: 0.5000000000010048 - 0.015738866124140063im\n", + "Eig: 0.5000000000011458 - 0.5268972835178625im\n", + "Eig: 0.5000000000012533 - 0.7035033703848474im\n", + "Eig: 0.5000000000013746 - 1.0048291799593458im\n", + "Eig: 0.5000000000015431 - 0.09864369626480596im\n", + "Eig: 0.5000000000016993 - 0.0009313715062656249im\n", + "Eig: 0.500000000001738 - 0.00389477864932247im\n", + "Eig: 0.5000000000025011 + 837.2951916250561im\n", + "Eig: 0.5000000000029691 - 1.5800894926728033im\n", + "Eig: 0.5000000000040603 - 2.893758443364143im\n", + "Eig: 0.5000000000086029 - 7.078330411354216im\n", + "Schur(Eig): -12.414391055721357 + 5.226922579051332im\n", + "Schur(Eig): 13.414391055729189 + 5.2269225790526095im\n", + "Schur(Eig): 2.5871709698474388 - 2.8404406841179854im\n", + "Schur(Eig): -1.5871709698367529 - 2.8404406841191885im\n", + "Schur(Eig): 0.8364100814306935 - 1.1480587902061123im\n", + "Schur(Eig): 0.16358991858028443 - 1.1480587902078874im\n", + "Schur(Eig): 0.4023476795951801 - 0.5877092334984895im\n", + "Schur(Eig): 0.5976523204155779 - 0.587709233498583im\n", + "Schur(Eig): 0.4638958507058264 - 0.3692732435668929im\n", + "Schur(Eig): 0.5361041493017054 - 0.3692732435665236im\n", + "Schur(Eig): 0.492614646296911 - 0.26754720718976166im\n", + "Schur(Eig): 0.5073853537107162 - 0.267547207187643im\n", + "Schur(Eig): 0.5000000000037167 - 0.2200965685391654im\n", + "Schur(Eig): 0.5000000000064049 - 0.1937158540961511im\n", + "Schur(Eig): 0.5000000000082445 - 0.16671638195349062im\n", + "Schur(Eig): 0.5000000000034363 - 0.14207142621055405im\n", + "Schur(Eig): 0.500000000008706 - 0.11936982537130715im\n", + "Schur(Eig): 0.5000000000115371 - 0.09864369698714791im\n", + "Schur(Eig): 0.5000000000164841 - 0.07989144103189551im\n", + "Schur(Eig): 0.5000000000118942 - 0.06311310404272963im\n", + "Schur(Eig): 0.5000000000098264 - 0.04830868354963337im\n", + "Schur(Eig): 0.5000000000484397 - 0.03547817639410573im\n", + "Schur(Eig): 0.50000000108917 - 0.0009313706017992215im\n", + "Schur(Eig): 0.5000000002487187 - 0.0038947784487597136im\n", + "Schur(Eig): 0.5000000000475523 - 0.008829999527985556im\n", + "Schur(Eig): 0.5000000000254943 - 0.02462157590437335im\n", + "Schur(Eig): 0.5000000000323227 - 0.01573886602115532im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): -11.667938813127611 + 4.846437539795649im\n", + "Schur(Eig): 12.667938813127936 + 4.846437539795365im\n", + "Schur(Eig): -1.4266918361496632 - 2.657431916607863im\n", + "Schur(Eig): 2.4266918361502308 - 2.6574319166082403im\n", + "Schur(Eig): 0.19105911427916825 - 1.0682441582516944im\n", + "Schur(Eig): 0.8089408857204733 - 1.0682441582506212im\n", + "Schur(Eig): 0.4106042543899601 - 0.5473560620594505im\n", + "Schur(Eig): 0.5893957456120842 - 0.5473560620559172im\n", + "Schur(Eig): 0.4675035804226512 - 0.3450002460940743im\n", + "Schur(Eig): 0.5324964195766715 - 0.3450002460945023im\n", + "Schur(Eig): 0.4970967603911124 - 0.25088534250638994im\n", + "Schur(Eig): 0.5029032396105566 - 0.2508853425060887im\n", + "Schur(Eig): 0.49999999999996286 - 0.205989472415028im\n", + "Schur(Eig): 0.49999999999887046 - 0.1799064553806543im\n", + "Schur(Eig): 0.4999999999988221 - 0.15403517620866608im\n", + "Schur(Eig): 0.4999999999984104 - 0.13036356973696725im\n", + "Schur(Eig): 0.4999999999993601 - 0.10864974413355909im\n", + "Schur(Eig): 0.4999999999977332 - 0.08891054875980786im\n", + "Schur(Eig): 0.49999999999560535 - 0.07114525457409482im\n", + "Schur(Eig): 0.500000000011739 - 0.05535387894334915im\n", + "Schur(Eig): 0.49999999997772127 - 0.041536419937349245im\n", + "Schur(Eig): 0.5000000000263519 - 0.0020570567147122536im\n", + "Schur(Eig): 0.500000000030057 - 0.006005543381833193im\n", + "Schur(Eig): 0.5000000000066954 - 0.011927479562594273im\n", + "Schur(Eig): 0.4999999999753506 - 0.01982323398240479im\n", + "Schur(Eig): 0.5000000000061455 - 0.029692874128058957im\n", + "Schur(Eig): NaN + NaN*im\n" + ] + } + ], + "source": [ + "# In this procedure, two methods are used to individually determine the eigenvalues:\n", + "# Method 1: [A-λB]x = 0 -> [B^{-1}A-λI]x = 0\n", + "# here, to be regular matrix B, boundary conditions are contained in 0<=n<=N-4 (so, elements used in A/B are N-4xN-4.)\n", + "# Method 2: [A-λB]x = 0 -> Ax=λBx -> Schur decomposition\n", + "# here, elements used in A2/B2 are NxN\n", + "\n", + "#using .Lilly1966_solver\n", + "using OffsetArrays # 配列要素番号を任意開始にするモジュール (行列計算時は 1 始まりの通常 matrix に copy する)\n", + "using LinearAlgebra\n", + "\n", + "#-- Drawing parameters \n", + "N = 56 # wavenumber in Chebyshev\n", + "Ncal = N + 1\n", + "\n", + "# initialize working arrays\n", + "# A/B: using to determine the eigenvalue by Method 1,\n", + "# A2/B2: using to determine the eigenvalue by Method 2.\n", + "A_J = reshape(zeros(Complex{Float64},Ncal,Ncal),Ncal,Ncal) # Matrix for calculation\n", + "B_J = reshape(zeros(Complex{Float64},Ncal,Ncal),Ncal,Ncal) # Matrix for calculation\n", + "A2_J = reshape(zeros(Complex{Float64},Ncal,Ncal),Ncal,Ncal) # Matrix for calculation\n", + "B2_J = reshape(zeros(Complex{Float64},Ncal,Ncal),Ncal,Ncal) # Matrix for calculation\n", + "### --- From here, supposing matrix element as row, col as in mathematical notation\n", + "Acal = reshape(zeros(Complex{Float64},Ncal,Ncal),Ncal,Ncal) # Matrix for calculation\n", + "Bcal = reshape(zeros(Complex{Float64},Ncal,Ncal),Ncal,Ncal) # Matrix for calculation\n", + "A2cal = reshape(zeros(Complex{Float64},Ncal,Ncal),Ncal,Ncal) # Matrix for calculation\n", + "B2cal = reshape(zeros(Complex{Float64},Ncal,Ncal),Ncal,Ncal) # Matrix for calculation\n", + "A = OffsetArray(Acal,0:N,0:N) # coefficient matrix\n", + "B = OffsetArray(Bcal,0:N,0:N) # coefficient matrix\n", + "A2 = OffsetArray(Acal,0:N,0:N) # coefficient matrix\n", + "B2 = OffsetArray(Bcal,0:N,0:N) # coefficient matrix\n", + "### --- From here, supposing matrix element as col, row as in julia notation (only using these arrays for matrix calc)\n", + "b = OffsetArray(zeros(Float64,Ncal),0:N) # A vector composed of Chebyshev coefficients for the basic state flow\n", + "c = OffsetArray(zeros(Float64,Ncal),0:N) # factorized vector (c_0=2, c_n=1, (n>0))\n", + "\n", + "# setting parameters\n", + "α = 1.0\n", + "Re = 10000.0\n", + "\n", + "# setting the Chebyshev coefficients for the basic flow\n", + "b[0] = 0.5\n", + "#b[2] = -b[0]\n", + "\n", + "# setting the factorized vector\n", + "c[0] = 2.0\n", + "c[1:N] .= 1.0\n", + "\n", + "# setting maximum even and odd numbers from N\n", + "Nemax = 0\n", + "Nomax = 1\n", + "for i in 0:2:N\n", + " Nemax = i\n", + "end\n", + "for i in 1:2:N\n", + " Nomax = i\n", + "end\n", + "println(\"Nemax = \", Nemax, \", Nomax = \", Nomax)\n", + "\n", + "# setting matrices\n", + "# 支配方程式から構築される部分 (チェビシェフ次数 0<=n<=N-4 の範囲)\n", + "# この段階では境界条件で成り立つ係数の関係は考慮しない\n", + "for n in 0:N-4\n", + " n2 = n^2\n", + " \n", + " # a_n\n", + " for p in n+4:2:N\n", + " p2 = p^2\n", + " A[n,p] = p * (p2 * (p2 - 4)^2 - 3 * n2 * p2 * (p2 - n2) - n2 * (n2 - 4)^2 - 48 * α^2 * (p2 - n2)) / 24.0\n", + " end\n", + " \n", + " A[n,n+2] = - 8 * α^2 * (n + 2) * (n + 1)\n", + " A[n,n] = α^4 * c[n]\n", + " \n", + " for p in 2:N\n", + " acoef = 0.0\n", + " p2 = p^2\n", + " for m in -(p-2):2:(p-2)\n", + " abs_n_m = abs(n - m)\n", + " m2 = m^2\n", + " if abs_n_m <= N\n", + " acoef = acoef + p * (p2 - m2) * c[abs_n_m] * b[abs_n_m]\n", + " end\n", + " end\n", + " A[n,p] = A[n,p] - complex(0.0,0.5 * α * Re * acoef)\n", + " end\n", + "\n", + " for p in -N:N\n", + " abs_n_p = abs(n - p)\n", + " abs_p = abs(p)\n", + " if abs_n_p <= N\n", + " A[n,abs_p] = A[n,abs_p] + complex(0.0,0.5 * α^3 * Re * c[abs_p] * c[abs_n_p] * b[abs_n_p])\n", + " end\n", + "\n", + " if abs_n_p+2 <= N\n", + " acoef = 0.0\n", + " for m in abs_n_p+2:2:N\n", + " m2 = m^2\n", + " acoef = acoef + m * (m2 - (n - p)^2) * b[m]\n", + " end\n", + " A[n,abs_p] = A[n,abs_p] + complex(0.0,0.5 * α * Re * c[abs_p] * acoef)\n", + " end\n", + " end\n", + "\n", + " # for lambda\n", + " for p in n+2:2:N\n", + " p2 = p^2\n", + " B[n,p] = complex(0.0, - α * Re * p * (p2 - n2))\n", + " end\n", + " B[n,n] = complex(0.0, α^3 * Re * c[n])\n", + "\n", + "end\n", + "\n", + "# Copy A/B -> A2/B2\n", + "A2 .= A\n", + "B2 .= B\n", + "\n", + "# For Method 1, \n", + "## a_n (N-3<=n<=N) from boundary conditions\n", + "# A(B)[n,N-3:N] に含まれている値 x 境界条件から得られる係数分を A(B)[n,0:N-4] に足し合わせる\n", + "for n in 0:N-4 # For each coefficient of T_n (0<=n<=N-4)\n", + " for p in 0:2:Nemax-4 # even mode\n", + " #println(n,\", \", p, \", \", A[n,p], \", \", (((Nemax - 2)^2 - p^2)/(4 * (Nemax - 1))) * A[n,Nemax], \",\", ((p^2 - Nemax^2)/(4 * (Nemax - 1))) * A[n,Nemax-2])\n", + " A[n,p] = A[n,p] + (((Nemax - 2)^2 - p^2)/(4 * (Nemax - 1))) * A[n,Nemax]\n", + " A[n,p] = A[n,p] + ((p^2 - Nemax^2)/(4 * (Nemax - 1))) * A[n,Nemax-2]\n", + " B[n,p] = B[n,p] + (((Nemax - 2)^2 - p^2)/(4 * (Nemax - 1))) * B[n,Nemax]\n", + " B[n,p] = B[n,p] + ((p^2 - Nemax^2)/(4 * (Nemax - 1))) * B[n,Nemax-2]\n", + " end\n", + " for p in 1:2:Nomax-4 # odd mode\n", + " A[n,p] = A[n,p] + (((Nomax - 2)^2 - p^2)/(4 * (Nomax - 1))) * A[n,Nomax]\n", + " A[n,p] = A[n,p] + ((p^2 - Nomax^2)/(4 * (Nomax - 1))) * A[n,Nomax-2]\n", + " B[n,p] = B[n,p] + (((Nomax - 2)^2 - p^2)/(4 * (Nomax - 1))) * B[n,Nomax]\n", + " B[n,p] = B[n,p] + ((p^2 - Nomax^2)/(4 * (Nomax - 1))) * B[n,Nomax-2]\n", + " end\n", + "end\n", + "\n", + "# copy MATRIX to MATRIXcal\n", + "Acal[1:Ncal-4,1:Ncal-4] .= A[0:N-4,0:N-4]\n", + "Bcal[1:Ncal-4,1:Ncal-4] .= B[0:N-4,0:N-4]\n", + "A_J[1:Ncal-4,1:Ncal-4] = transpose(Acal[1:Ncal-4,1:Ncal-4]) # A_J = Transpose type\n", + "B_J[1:Ncal-4,1:Ncal-4] = transpose(Bcal[1:Ncal-4,1:Ncal-4]) # B_J = Transpose type\n", + "\n", + "Binv = inv(B_J[1:Ncal-4,1:Ncal-4])\n", + "Ccal = inv(B_J[1:Ncal-4,1:Ncal-4]) * A_J[1:Ncal-4,1:Ncal-4]\n", + "Dcal = eigvals(Ccal)\n", + "F = imag.(Dcal)\n", + "for i in 1:N-4\n", + " println(\"Eig: \", Dcal[i])\n", + "end\n", + "\n", + "# For Method 2, a_n (N-3<=n<=N) from boundary conditions\n", + "## [A-λB]x = 0 -> Ax=λBx にして A, B を QZ 分解 (一般化 Schur 分解) して上三角行列から固有値を求める\n", + "for p in 0:2:N # even mode\n", + " A2[N-3,p] = A2[N-3,p] + 1.0\n", + " A2[N-1,p] = A2[N-1,p] + p^2\n", + "end\n", + "for p in 1:2:N # odd mode\n", + " A2[N-2,p] = A2[N-2,p] + 1.0\n", + " A2[N,p] = A2[N,p] + p^2\n", + "end\n", + "\n", + "# copy MATRIX to MATRIXcal\n", + "A2cal[1:Ncal,1:Ncal] .= A2[0:N,0:N]\n", + "B2cal[1:Ncal,1:Ncal] .= B2[0:N,0:N]\n", + "# Convert Transpose type to Matrix (Array) type (ドット演算して要素だけコピーしている)\n", + "A2_J .= transpose(A2cal) # A_J = Transpose type\n", + "B2_J .= transpose(B2cal) # B_J = Transpose type\n", + "\n", + "F = schur(A2_J,B2_J)\n", + "D2cal = F.α ./ F.β\n", + "for i in 1:Ncal\n", + " println(\"Schur(Eig): \", D2cal[i])\n", + "end\n", + "\n", + "#G = eigvecs(Ccal)\n", + "#H = inv(G) * Ccal * G\n", + "#ncmax = findmax(F)[1]\n", + "#Ae = OffsetArray(A[0:div(N,2),0:div(N,2)],0:div(N,2),0:div(N,2))\n", + "#Ao = OffsetArray(A[1:div(N,2),1:div(N,2)],1:div(N,2),1:div(N,2))\n", + "#Be = OffsetArray(B[0:div(N,2),0:div(N,2)],0:div(N,2),0:div(N,2))\n", + "#Bo = OffsetArray(B[1:div(N,2),1:div(N,2)],1:div(N,2),1:div(N,2))\n", + "#for n in 0:div(N,2)\n", + "# for p in 0:div(N,2)\n", + "# Ae[n,p] = A[2*n,2*p]\n", + "# Be[n,p] = B[2*n,2*p]\n", + "# end\n", + "#end\n", + "#for n in 1:div(N,2)\n", + "# for p in 1:div(N,2)\n", + "# Ao[n,p] = A[2*n-1,2*p-1]\n", + "# Bo[n,p] = B[2*n-1,2*p-1]\n", + "# end\n", + "#end\n", + "#Acal[1:div(N,2)-1,1:div(N,2)-1] .= Ae[0:div(N,2)-2,0:div(N,2)-2]\n", + "#Bcal[1:div(N,2)-1,1:div(N,2)-1] .= Be[0:div(N,2)-2,0:div(N,2)-2]\n", + "#Ccal = inv(Bcal[1:div(N,2)-1,1:div(N,2)-1]) * Acal[1:div(N,2)-1,1:div(N,2)-1]\n", + "##Ccal = inv(Bcal) * Acal\n", + "#Dcal = eigvals(Ccal)\n", + "#F = imag.(Dcal)\n", + "#for i in 1:div(N,2)-1\n", + "# println(F[i])\n", + "#end\n", + "\n", + "#Acal[1:div(N,2)-2,1:div(N,2)-2] .= Ao[1:div(N,2)-2,1:div(N,2)-2]\n", + "#Bcal[1:div(N,2)-2,1:div(N,2)-2] .= Bo[1:div(N,2)-2,1:div(N,2)-2]\n", + "#Ccal = inv(Bcal[1:div(N,2)-2,1:div(N,2)-2]) * Acal[1:div(N,2)-2,1:div(N,2)-2]\n", + "##Ccal = inv(Bcal) * Acal\n", + "#Dcal = eigvals(Ccal)\n", + "#F = imag.(Dcal)\n", + "#for i in 1:div(N,2)-2\n", + "# println(F[i])\n", + "#end\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "1255c473", + "metadata": {}, + "outputs": [ + { + "ename": "UndefVarError", + "evalue": "\u001b[91mUndefVarError: Ncal not defined\u001b[39m", + "output_type": "error", + "traceback": [ + "\u001b[91mUndefVarError: Ncal not defined\u001b[39m", + "", + "Stacktrace:", + " [1] top-level scope at In[2]:9" + ] + } + ], + "source": [ + "#-- Drawing\n", + "##########\n", + "# Plot #\n", + "##########\n", + "using PyPlot\n", + "\n", + "rc(\"font\", family=\"IPAPGothic\")\n", + "fig = figure(\"pyplot_majorminor\",figsize=(7,5))\n", + "\n", + "x = real.(Dcal[3:Ncal-4])\n", + "y = imag.(Dcal[3:Ncal-4])\n", + "\n", + "ax = gca()\n", + "\n", + "cp = scatter(x,y)\n", + "\n", + "ax = gca()\n", + "\n", + "xlabel(\"X\")\n", + "ylabel(\"Y\")\n", + "#if i == 1 #div(ne,2)-1\n", + "# plt.colorbar(cp)\n", + "#end\n", + "grid(\"on\")\n", + "\n", + "fig.canvas.draw()\n", + "gcf() # Needed for IJulia to plot inline\n", + "\n", + "#PyPlot.title(\"Re = Re[i]\")\n", + "\n", + "fname = \"test.png\"\n", + "#savefig(\"B-init_2d.pdf\")\n", + "savefig(fname)\n", + "#psi\n" + ] + }, + { + "cell_type": "code", + "execution_count": 66, + "id": "623f8c38", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Figure(PyObject
)" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "#-- Drawing\n", + "##########\n", + "# Plot #\n", + "##########\n", + "#using PyPlot\n", + "\n", + "#rc(\"font\", family=\"IPAPGothic\")\n", + "#fig = figure(\"pyplot_majorminor\",figsize=(7,5))\n", + "\n", + "nx = Ncal-4\n", + "ny = Ncal-4\n", + "\n", + "xax = reshape(zeros(nx,ny),nx,ny)\n", + "yax = reshape(zeros(nx,ny),nx,ny)\n", + "\n", + "for i in 1:nx\n", + " xax[i,1:ny] .= i # xmin + (xmax-xmin) * (i-1) / (nx-1)\n", + "end\n", + "for i in 1:ny\n", + " yax[1:nx,i] .= i # ymin + (ymax-ymin) * (i-1) / (ny-1)\n", + "# yax[1:nx,i] .= Re[i] # ymin + (ymax-ymin) * (i-1) / (ny-1)\n", + "end\n", + "\n", + "cont_val = reshape(zeros(nx,ny),nx,ny)\n", + "shade_val = reshape(zeros(nx,ny),nx,ny)\n", + "#for i in 1:nr#div(nr,2)-1:div(nr,2)+1\n", + " cont_val = imag.(Ccal[1:nx,1:ny])\n", + " #for j in 1:ne\n", + " shade_val = imag.(A_J[1:nx,1:ny])\n", + " #end\n", + "#draw_num2 = 2\n", + "#if draw_num2 == 1\n", + "# cp = contourf(xax[1:nx,1:nx], xax'[1:nx,1:nx], Pftime[1:nx,1:nx,1], levels=[-15.0, -10.0, -5.0, 0.0, 5.0, 10.0, 15.0, 20.0], origin=\"image\", cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + " #cp = contourf(xax[1:nx,1:nx], xax'[1:nx,1:nx], Pftime[1:nx,1:nx,2], levels=[-2.0, -1.0, -0.5, -0.25, 0.25, 0.5, 1.0, 1.5, 2.0], origin=\"image\", cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + " #cp = contourf(xax[1:nx,1:no], xax'[1:nx,1:no], Pftime[1:nx,1:nx,2]*Hop'*inv(Hop*Pftime[1:nx,1:nx,2]*Hop'+Ro), origin=\"image\", cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + "#elseif draw_num2 == 2\n", + " ##cp = contour(xax[1:nx,1:ny], yax[1:nx,1:ny], shade_val[1:nx,1:ny], 8, levels=[-0.01, -0.005, -0.003, -0.002, -0.001, -0.0005, 0.0, 0.0005, 0.001, 0.002, 0.003, 0.005, 0.01])\n", + " cp = contourf(xax[1:nx,1:ny], yax[1:nx,1:ny], cont_val[1:nx,1:ny], cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + " ##cp = contourf(xax[1:nx,1:ny], yax[1:nx,1:ny], cont_val[1:nx,1:ny], 8, levels=[-1.0, -0.5, -0.2, -0.1, 0.0, 0.1, 0.2, 0.5, 1.0], cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + " #cp = contourf(xax[1:nx,1:ny], yax[1:nx,1:ny], shade_val[1:nx,1:ny], 10, cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + "#elseif draw_num2 == 3\n", + "# cp = contourf(xoax[1:nx,1:nto-1], toax[1:nx,1:nto-1], evec_max[1:nx,1:nto-1], 10, cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + "#elseif draw_num2 == 4\n", + "# cp = contourf(xoax[1:nx,1:40], toax[1:nx,1:40], x_inc[1:nx,1:40], levels=[-2.0, -1.0, -0.5, -0.25, 0.25, 0.5, 1.0, 1.5, 2.0], cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + "#elseif draw_num2 == 5\n", + "# cp = contourf(xoax[1:no,1:40], toax[1:no,1:40], y_innov[1:no,1:40], levels=[-2.0, -1.0, -0.5, -0.25, 0.25, 0.5, 1.0, 1.5, 2.0], cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + "#end\n", + "#ax.label(cp, inline=1, fontsize=10)\n", + "#legend()\n", + "ax = gca()\n", + "\n", + "xlabel(\"X\")\n", + "ylabel(\"Y\")\n", + "#if i == 1 #div(ne,2)-1\n", + " plt.colorbar(cp)\n", + "#end\n", + "grid(\"on\")\n", + "\n", + "PyPlot.title(\"Re = Re[i]\")\n", + "\n", + "#########################\n", + "# Set tick dimensions #\n", + "#########################\n", + "#ax.xaxis.set_tick_params(which=\"major\",length=5,width=2,labelsize=10)\n", + "#ax.xaxis.set_tick_params(which=\"minor\",length=5,width=2)\n", + "\n", + "fig.canvas.draw() # Update the figure\n", + "gcf() # Needed for IJulia to plot inline\n", + "#fname = \"Re\" * string(1) * \".png\"\n", + "fname = \"test.png\"\n", + "#savefig(\"B-init_2d.pdf\")\n", + "savefig(fname)\n", + "#psi\n", + "#end\n", + "#println(shade_val[1:na,1])" + ] + }, + { + "cell_type": "code", + "execution_count": 47, + "id": "af6a72bd", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ComplexF64[0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im]\n" + ] + }, + { + "data": { + "text/plain": [ + "GeneralizedSchur{ComplexF64, Matrix{ComplexF64}}\n", + "S factor:\n", + "2×2 Matrix{ComplexF64}:\n", + " 2.23607-0.0im 0.0+0.0im\n", + " 0.0+0.0im 2.23607+0.0im\n", + "T factor:\n", + "2×2 Matrix{ComplexF64}:\n", + " 2.23607+0.0im 0.0+0.0im\n", + " 0.0+0.0im 2.23607+0.0im\n", + "Q factor:\n", + "2×2 Matrix{ComplexF64}:\n", + " -0.447214-0.894427im 0.0+0.0im\n", + " 0.0+0.0im 1.0+0.0im\n", + "Z factor:\n", + "2×2 Matrix{ComplexF64}:\n", + " -1.0-0.0im 0.0+0.0im\n", + " -0.0-0.0im 0.894427+0.447214im\n", + "α:\n", + "2-element Vector{ComplexF64}:\n", + " 2.23606797749979 - 0.0im\n", + " 2.23606797749979 + 0.0im\n", + "β:\n", + "2-element Vector{ComplexF64}:\n", + " 2.23606797749979 + 0.0im\n", + " 2.23606797749979 + 0.0im" + ] + }, + "execution_count": 47, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "a33 = reshape(zeros(Complex{Float64},2,2),2,2)\n", + "c33 = reshape(zeros(Complex{Float64},2,2),2,2)\n", + "println(a33)\n", + "a33[1,1] = 1.0 + 2.0im\n", + "a33[2,2] = 2.0 - 1.0im\n", + "b33 = transpose(a33)\n", + "c33 .= b33\n", + "schur(a33,c33)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b28e1119", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.6.1", + "language": "julia", + "name": "julia-1.6" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.6.1" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Lilly1966/.ipynb_checkpoints/OSE_solve_finite-checkpoint.ipynb b/Lilly1966/.ipynb_checkpoints/OSE_solve_finite-checkpoint.ipynb new file mode 100644 index 0000000..fc0b486 --- /dev/null +++ b/Lilly1966/.ipynb_checkpoints/OSE_solve_finite-checkpoint.ipynb @@ -0,0 +1,731 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "b5e749e8", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Main.Lilly1966_solver" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "module Lilly1966_solver\n", + "\n", + "using LinearAlgebra\n", + "\n", + "export full_EQ_set, OSE_set\n", + "\n", + "function full_EQ_set( α_in, ε_in, Re_in, N_in )\n", + "# Solve the complete set of the equations (2.25)-(2.28) in Lilly (1966)\n", + " #-- arguments input:\n", + " # α_in: wavenumber for y-direction [non-dim,Vector{Float64}]\n", + " # ε_in: flow angle for Vg [degree,Vector{Float64}]\n", + " # Re_in: Reinolds number [non-dim,Vector{Float64}]\n", + " # N_in: grid numbers in vertical [1,Int64]\n", + " #-- arguments output:\n", + " # c_out: wave speed [non-dim,Array{Complex{Float64},size(α),size(ε),size(Re)}]\n", + "\n", + " # set parameters\n", + " d2r = π / 180.0 # degree to radian\n", + " na = size(α_in)[1]\n", + " ne = size(ε_in)[1]\n", + " nr = size(Re_in)[1]\n", + " N = N_in\n", + "\n", + " # initialize variables and matrix\n", + " ntot = 2*N-1 # matrix dimension\n", + " c_out = reshape(zeros(Complex{Float64},na,ne,nr),na,ne,nr)\n", + "\n", + " # initialize working arrays in this function\n", + " ### --- From here, supposing matrix element as row, col as in mathematical notation\n", + " A = reshape(zeros(Complex{Float64},ntot,ntot),ntot,ntot) # linear operator in differential equations\n", + " B = reshape(zeros(Complex{Float64},ntot,ntot),ntot,ntot) # linear operator in differential equations\n", + " A_φφ = reshape(zeros(Complex{Float64},N-1,N-1),N-1,N-1) # partial matrix (phi operator on phi EQ.) in A\n", + " A_φμ = reshape(zeros(Complex{Float64},N-1,N),N-1,N) # partial matrix (mu operator on phi EQ.) in A\n", + " A_μφ = reshape(zeros(Complex{Float64},N,N-1),N,N-1) # partial matrix (phi operator on mu EQ.) in A\n", + " A_μμ = reshape(zeros(Complex{Float64},N,N),N,N) # partial matrix (mu operator on mu EQ.) in A\n", + " B_φφ = reshape(zeros(Complex{Float64},N-1,N-1),N-1,N-1) # partial matrix (phi operator on phi EQ.) in B\n", + " B_φμ = reshape(zeros(Complex{Float64},N-1,N),N-1,N) # partial matrix (mu operator on phi EQ.) in B == 0\n", + " B_μφ = reshape(zeros(Complex{Float64},N,N-1),N,N-1) # partial matrix (phi operator on mu EQ.) in B == 0\n", + " B_μμ = reshape(zeros(Complex{Float64},N,N),N,N) # partial matrix (mu operator on mu EQ.) in B\n", + " ### --- From here, supposing matrix element as col, row as in julia notation (only using these arrays for matrix calc)\n", + " A_J = reshape(zeros(Complex{Float64},ntot,ntot),ntot,ntot) # linear operator in differential equations\n", + " B_J = reshape(zeros(Complex{Float64},ntot,ntot),ntot,ntot) # linear operator in differential equations\n", + "\n", + " U_bar = reshape(zeros(N,1),N,1) # Mean flow in x-direction defined at i+1/2\n", + " V_bar = reshape(zeros(N-1,1),N-1,1) # Mean flow in y-direction defined at i\n", + " Vh_bar = reshape(zeros(N,1),N,1) # Mean flow in y-direction defined at i+1/2\n", + " dUdz = reshape(zeros(N,1),N,1) # Mean flow shear in x-direction defined at i+1/2\n", + " dVdz = reshape(zeros(N-1,1),N-1,1) # Mean flow shear in y-direction defined at i\n", + " d2Vdz2 = reshape(zeros(N-1,1),N-1,1) # Mean flow shear in y-direction defined at i\n", + "\n", + " #-- Start loops for varying parameters (α, ε, Re)\n", + " for k in 1:nr\n", + " for j in 1:ne\n", + " for i in 1:na\n", + " α = α_in[i]\n", + " ε = ε_in[j] * d2r\n", + " Re = Re_in[k]\n", + " Δ = 1.0/sqrt(α * N) # grid spacing\n", + " \n", + " # set grid spacing and mean-wind profiles in vertical\n", + " z_p = [i*Δ for i in 1:N-1] # non-dimensional vertical coordinate for φ\n", + " zh_p = [(i-0.5)*Δ for i in 1:N] # non-dimensional vertical coordinate for μ\n", + " z_p = reshape(z_p,N-1,1)\n", + " zh_p = reshape(zh_p,N,1)\n", + " \n", + " #-- From here, considering to move them as external argument\n", + " U_bar[1:N,1] .= cos(ε) .- exp.(-zh_p[1:N,1]) .* cos.(zh_p[1:N,1] .+ ε)\n", + " V_bar[1:N-1,1] .= -sin(ε) .+ exp.(-z_p[1:N-1,1]) .* sin.(z_p[1:N-1,1] .+ ε)\n", + " Vh_bar[1:N,1] .= -sin(ε) .+ exp.(-zh_p[1:N,1]) .* sin.(zh_p[1:N,1] .+ ε)\n", + " dUdz[1:N,1] .= exp.(-zh_p[1:N,1]) .* cos.(zh_p[1:N,1] .+ ε) .+ exp.(-zh_p[1:N,1]) .* sin.(zh_p[1:N,1] .+ ε)\n", + " dVdz[1:N-1,1] .= -exp.(-z_p[1:N-1,1]) .* sin.(z_p[1:N-1,1] .+ ε) .+ exp.(-z_p[1:N-1,1]) .* cos.(z_p[1:N-1,1] .+ ε)\n", + " d2Vdz2[1:N-1,1] .= .- 2.0 .* exp.(-z_p[1:N-1,1]) .* cos.(z_p[1:N-1,1] .+ ε)\n", + "\n", + " # setting matrix elements in A and B\n", + " Δ2 = Δ*Δ\n", + " Δ4 = Δ2*Δ2\n", + " αΔ = α*Δ\n", + " V_infl = -sin(ε) + exp(-(0.5*π-ε))\n", + " # μ[1] μ[2] μ[3] μ[4]\n", + " # ----- 0 --------- 1 --------- 2 --------- 3 -- i for μ = μ[i+1]\n", + " # ------|-----------|-----------|-----------|--- for μ = μ[i+1]\n", + " # 0---------- 1 --------- 2 --------- 3 -------- i for φ = φ[i]\n", + " # |-----------|-----------|-----------|--------- for φ = φ[i]\n", + "\n", + " # For φ EQ.\n", + " for i in 3:N-3\n", + " # println(\"$i\")\n", + " A_φφ[i,i+2] = complex(1.0,α*Re*V_bar[i,1]*Δ2/12.0)\n", + " A_φφ[i,i-2] = A_φφ[i,i+2]\n", + " A_φφ[i,i+1] = complex(-4.0-2.0*((αΔ)^2),-4.0*α*Re*V_bar[i,1]*Δ2/3.0)\n", + " A_φφ[i,i-1] = A_φφ[i,i+1]\n", + " A_φφ[i,i] = complex(6.0+4.0*((αΔ)^2)+(α^4-d2Vdz2[i,1])*Δ4,α*Re*V_bar[i,1]*(2.5+(αΔ)^2)*Δ2)\n", + " A_φμ[i,i] = -2.0*(Δ^3) # A^{φ,μ}_{i,i-1/2}\n", + " A_φμ[i,i+1] = +2.0*(Δ^3) # A^{φ,μ}_{i,i+1/2}\n", + " B_φφ[i,i+2] = complex(0.0,α*Re*Δ2/12.0)\n", + " B_φφ[i,i-2] = B_φφ[i,i+2]\n", + " B_φφ[i,i+1] = -complex(0.0,4.0*α*Re*Δ2/3.0)\n", + " B_φφ[i,i-1] = B_φφ[i,i+1]\n", + " B_φφ[i,i] = complex(0.0,α*Re*(2.5+(αΔ)^2)*Δ2)\n", + " end\n", + "\n", + " # For μ EQ.\n", + " for i in 3:N-2\n", + " A_μφ[i,i] = -complex(2.0*Δ,9.0*α*Re*dUdz[i,1]*Δ2/16.0) # A^{μ,φ}_{i+1/2,i+1}\n", + " A_μφ[i,i-2] = complex(0.0,α*Re*dUdz[i,1]*Δ2/16.0) # A^{μ,φ}_{i+1/2,i-1}\n", + " A_μφ[i,i+1] = A_μφ[i,i-2] # A^{μ,φ}_{i+1/2,i+2}\n", + " A_μφ[i,i-1] = complex(2.0*Δ,-9.0*α*Re*dUdz[i,1]*Δ2/16.0) # A^{μ,φ}_{i+1/2,i}\n", + " A_μμ[i,i+1] = 1.0 # A^{μ,μ}_{i+1/2,i-1/2}\n", + " A_μμ[i,i-1] = 1.0 # A^{μ,μ}_{i+1/2,i+3/2}\n", + " A_μμ[i,i] = -complex(2.0+(αΔ)^2,α*Re*Vh_bar[i,1]*Δ2)\n", + " B_μμ[i,i] = -complex(0.0,α*Re*Δ2)\n", + " end\n", + "\n", + " # Adjacent of boundaries\n", + " # For φ EQ. (at z=z_1)\n", + " A_φφ[1,3] = complex(1.0,α*Re*V_bar[1,1]*Δ2/12.0)\n", + " A_φφ[1,2] = complex(-4.0-2.0*((αΔ)^2),-4.0*α*Re*V_bar[1,1]*Δ2/3.0)\n", + " A_φφ[1,1] = complex(7.0+4.0*((αΔ)^2)+(α^4-d2Vdz2[1,1])*Δ4,α*Re*V_bar[1,1]*(31.0/12.0+(αΔ)^2)*Δ2)\n", + " A_φμ[1,1] = -2.0*(Δ^3)\n", + " A_φμ[1,2] = +2.0*(Δ^3)\n", + " B_φφ[1,3] = complex(0.0,α*Re*Δ2/12.0)\n", + " B_φφ[1,2] = -complex(0.0,4.0*α*Re*Δ2/3.0)\n", + " B_φφ[1,1] = complex(0.0,α*Re*(31.0/12.0+(αΔ)^2)*Δ2)\n", + " # For μ EQ. (at z=z_{1/2})\n", + " A_μφ[1,1] = -complex(2.0*Δ,0.5*α*Re*dUdz[1,1]*Δ2)\n", + " A_μφ[1,2] = complex(0.0,α*Re*dUdz[1,1]*Δ2/16.0)\n", + " A_μμ[1,2] = 1.0\n", + " A_μμ[1,1] = -complex(3.0+(αΔ)^2,α*Re*Vh_bar[1,1]*Δ2)\n", + " B_μμ[1,1] = -complex(0.0,α*Re*Δ2)\n", + "\n", + " # For φ EQ. (at z=z_2)\n", + " A_φφ[2,4] = complex(1.0,α*Re*V_bar[2,1]*Δ2/12.0)\n", + " A_φφ[2,3] = complex(-4.0-2.0*((αΔ)^2),-4.0*α*Re*V_bar[2,1]*Δ2/3.0)\n", + " A_φφ[2,1] = A_φφ[2,3]\n", + " A_φφ[2,2] = complex(6.0+4.0*((αΔ)^2)+(α^4-d2Vdz2[2,1])*Δ4,α*Re*V_bar[2,1]*(2.5+(αΔ)^2)*Δ2)\n", + " A_φμ[2,2] = -2.0*(Δ^3)\n", + " A_φμ[2,3] = +2.0*(Δ^3)\n", + " B_φφ[2,4] = complex(0.0,α*Re*Δ2/12.0)\n", + " B_φφ[2,3] = -complex(0.0,4.0*α*Re*Δ2/3.0)\n", + " B_φφ[2,1] = B_φφ[2,3]\n", + " B_φφ[2,2] = complex(0.0,α*Re*(2.5+(αΔ)^2)*Δ2)\n", + " # For μ EQ. (at z=z_{1+1/2})\n", + " A_μφ[2,2] = -complex(2.0*Δ,9.0*α*Re*dUdz[2,1]*Δ2/16.0)\n", + " A_μφ[2,3] = complex(0.0,α*Re*dUdz[2,1]*Δ2/16.0)\n", + " A_μφ[2,1] = complex(2.0*Δ,-9.0*α*Re*dUdz[2,1]*Δ2/16.0)\n", + " A_μμ[2,3] = 1.0\n", + " A_μμ[2,1] = 1.0\n", + " A_μμ[2,2] = -complex(2.0+(αΔ)^2,α*Re*Vh_bar[2,1]*Δ2)\n", + " B_μμ[2,2] = -complex(0.0,α*Re*Δ2)\n", + "\n", + " # For φ EQ. (at z=z_{N-2})\n", + " A_φφ[N-2,N-4] = complex(1.0,α*Re*V_bar[N-2,1]*Δ2/12.0)\n", + " A_φφ[N-2,N-1] = complex(-4.0-2.0*((αΔ)^2),-4.0*α*Re*V_bar[N-2,1]*Δ2/3.0)\n", + " A_φφ[N-2,N-3] = A_φφ[N-2,N-1]\n", + " A_φφ[N-2,N-2] = complex(6.0+4.0*((αΔ)^2)+(α^4-d2Vdz2[N-2,1])*Δ4,α*Re*V_bar[N-2,1]*(2.5+(αΔ)^2)*Δ2)\n", + " A_φμ[N-2,N-2] = -2.0*(Δ^3)\n", + " A_φμ[N-2,N-1] = +2.0*(Δ^3)\n", + " B_φφ[N-2,N-4] = complex(0.0,α*Re*Δ2/12.0)\n", + " B_φφ[N-2,N-1] = -complex(0.0,4.0*α*Re*Δ2/3.0)\n", + " B_φφ[N-2,N-3] = B_φφ[N-2,N-1]\n", + " B_φφ[N-2,N-2] = complex(0.0,α*Re*(2.5+(αΔ)^2)*Δ2)\n", + " # For μ EQ. (at z=z_{N-1-1/2})\n", + " A_μφ[N-1,N-1] = -complex(2.0*Δ,9.0*α*Re*dUdz[N-1,1]*Δ2/16.0)\n", + " A_μφ[N-1,N-3] = complex(0.0,α*Re*dUdz[N-1,1]*Δ2/16.0)\n", + " A_μφ[N-1,N-2] = complex(2.0*Δ,-9.0*α*Re*dUdz[N-1,1]*Δ2/16.0)\n", + " A_μμ[N-1,N-2] = 1.0\n", + " A_μμ[N-1,N] = 1.0\n", + " A_μμ[N-1,N-1] = -complex(2.0+(αΔ)^2,α*Re*Vh_bar[N-1,1]*Δ2)\n", + " B_μμ[N-1,N-1] = -complex(0.0,α*Re*Δ2)\n", + "\n", + " # For φ EQ. (at z=z_{N-1})\n", + " A_φφ[N-1,N-3] = complex(1.0,α*Re*V_bar[N-1,1]*Δ2/12.0)\n", + " A_φφ[N-1,N-2] = complex(-4.0-2.0*((αΔ)^2),-4.0*α*Re*V_bar[N-1,1]*Δ2/3.0)\n", + " A_φφ[N-1,N-1] = complex(5.0+4.0*((αΔ)^2)+(α^4-d2Vdz2[N-1,1])*Δ4,α*Re*V_bar[N-1,1]*(29.0/12.0+(αΔ)^2)*Δ2)\n", + " A_φμ[N-1,N-1] = -2.0*(Δ^3)\n", + " A_φμ[N-1,N] = +2.0*(Δ^3)\n", + " B_φφ[N-1,N-3] = complex(0.0,α*Re*Δ2/12.0)\n", + " B_φφ[N-1,N-2] = -complex(0.0,4.0*α*Re*Δ2/3.0)\n", + " B_φφ[N-1,N-1] = complex(0.0,α*Re*(29.0/12.0+(αΔ)^2)*Δ2)\n", + " # For μ EQ. (at z=z_{N-1/2})\n", + " A_μφ[N,N-2] = complex(0.0,α*Re*dUdz[N]*Δ2/16.0)\n", + " A_μφ[N,N-1] = complex(2.0*Δ,-5.0*α*Re*dUdz[N]*Δ2/8.0)\n", + " A_μμ[N,N-1] = 1.0\n", + " A_μμ[N,N] = -complex(1.0+(αΔ)^2,α*Re*Vh_bar[N]*Δ2)\n", + " B_μμ[N,N] = -complex(0.0,α*Re*Δ2)\n", + "\n", + " # assigned partial matrices A^ and B^ to A and B\n", + " A[1:N-1,1:N-1] = A_φφ[1:N-1,1:N-1]\n", + " A[1:N-1,N:2*N-1] = A_φμ[1:N-1,1:N]\n", + " A[N:2*N-1,1:N-1] = A_μφ[1:N,1:N-1]\n", + " A[N:2*N-1,N:2*N-1] = A_μμ[1:N,1:N]\n", + " # non-diagonal parts in B are zero\n", + " B[1:N-1,1:N-1] = B_φφ[1:N-1,1:N-1]\n", + " B[N:2*N-1,N:2*N-1] = B_μμ[1:N,1:N]\n", + "\n", + " A_J = transpose(A) # mathematical A -> Julia A_J\n", + " B_J = transpose(B) # mathematical A -> Julia A_J\n", + " C_J = inv(B_J) * A_J\n", + " D_J = eigvals(C_J)\n", + " F_J = imag.(D_J)\n", + " E_J = inv(B_J) * B_J\n", + " nc_max = findmax(F_J)[2]\n", + " c_out[i,j,k] = D_J[nc_max]\n", + "\n", + " end\n", + " end\n", + " end\n", + "\n", + "#println(D_J[findmax(F_J)[2]])\n", + "#println(F_J)\n", + "#println(V_infl)\n", + "\n", + " return c_out\n", + "\n", + "end # function full_EQ_set\n", + " \n", + "function OSE_set( α_in, ε_in, Re_in, N_in )\n", + "# Solve the approximation (OSE) set of the equations (2.25), (2.27)-(2.28) in Lilly (1966)\n", + " #-- arguments input:\n", + " # α_in: wavenumber for y-direction [non-dim,Vector{Float64}]\n", + " # ε_in: flow angle for Vg [degree,Vector{Float64}]\n", + " # Re_in: Reinolds number [non-dim,Vector{Float64}]\n", + " # N_in: grid numbers in vertical [1,Int64]\n", + " #-- arguments output:\n", + " # c_out: wave speed [non-dim,Array{Complex{Float64},size(α),size(ε),size(Re)}]\n", + "\n", + " # set parameters\n", + " d2r = π / 180.0 # degree to radian\n", + " na = size(α_in)[1]\n", + " ne = size(ε_in)[1]\n", + " nr = size(Re_in)[1]\n", + " N = N_in\n", + "\n", + " # initialize variables and matrix\n", + " ntot = N-1 # matrix dimension\n", + " c_out = reshape(zeros(Complex{Float64},(na,ne,nr)),na,ne,nr)\n", + "\n", + " # initialize working arrays in this function\n", + " ### --- From here, supposing matrix element as row, col as in mathematical notation\n", + " A = reshape(zeros(Complex{Float64},ntot,ntot),ntot,ntot) # linear operator in differential equations\n", + " B = reshape(zeros(Complex{Float64},ntot,ntot),ntot,ntot) # linear operator in differential equations\n", + " A_φφ = reshape(zeros(Complex{Float64},N-1,N-1),N-1,N-1) # partial matrix (phi operator on phi EQ.) in A\n", + " B_φφ = reshape(zeros(Complex{Float64},N-1,N-1),N-1,N-1) # partial matrix (phi operator on phi EQ.) in B\n", + " ### --- From here, supposing matrix element as col, row as in julia notation (only using these arrays for matrix calc)\n", + " A_J = reshape(zeros(Complex{Float64},ntot,ntot),ntot,ntot) # linear operator in differential equations\n", + " B_J = reshape(zeros(Complex{Float64},ntot,ntot),ntot,ntot) # linear operator in differential equations\n", + "\n", + " U_bar = reshape(zeros(N,1),N,1) # Mean flow in x-direction defined at i+1/2\n", + " V_bar = reshape(zeros(N-1,1),N-1,1) # Mean flow in y-direction defined at i\n", + " Vh_bar = reshape(zeros(N,1),N,1) # Mean flow in y-direction defined at i+1/2\n", + " dUdz = reshape(zeros(N,1),N,1) # Mean flow shear in x-direction defined at i+1/2\n", + " dVdz = reshape(zeros(N-1,1),N-1,1) # Mean flow shear in y-direction defined at i\n", + " d2Vdz2 = reshape(zeros(N-1,1),N-1,1) # Mean flow shear in y-direction defined at i\n", + "\n", + " #-- Start loops for varying parameters (α, ε, Re)\n", + " for k in 1:nr\n", + " for j in 1:ne\n", + " for i in 1:na\n", + " α = α_in[i]\n", + " ε = ε_in[j] * d2r\n", + " Re = Re_in[k]\n", + " Δ = 1.0/sqrt(α * N) # grid spacing\n", + " \n", + " # set grid spacing and mean-wind profiles in vertical\n", + " z_p = [i*Δ for i in 1:N-1] # non-dimensional vertical coordinate for φ\n", + " zh_p = [(i-0.5)*Δ for i in 1:N] # non-dimensional vertical coordinate for μ\n", + " z_p = reshape(z_p,N-1,1)\n", + " zh_p = reshape(zh_p,N,1)\n", + " \n", + " #-- From here, considering to move them as external argument\n", + " U_bar[1:N,1] .= cos(ε) .- exp.(-zh_p[1:N,1]) .* cos.(zh_p[1:N,1] .+ ε)\n", + " V_bar[1:N-1,1] .= -sin(ε) .+ exp.(-z_p[1:N-1,1]) .* sin.(z_p[1:N-1,1] .+ ε)\n", + " Vh_bar[1:N,1] .= -sin(ε) .+ exp.(-zh_p[1:N,1]) .* sin.(zh_p[1:N,1] .+ ε)\n", + " dUdz[1:N,1] .= exp.(-zh_p[1:N,1]) .* cos.(zh_p[1:N,1] .+ ε) .+ exp.(-zh_p[1:N,1]) .* sin.(zh_p[1:N,1] .+ ε)\n", + " dVdz[1:N-1,1] .= -exp.(-z_p[1:N-1,1]) .* sin.(z_p[1:N-1,1] .+ ε) .+ exp.(-z_p[1:N-1,1]) .* cos.(z_p[1:N-1,1] .+ ε)\n", + " d2Vdz2[1:N-1,1] .= .- 2.0 .* exp.(-z_p[1:N-1,1]) .* cos.(z_p[1:N-1,1] .+ ε)\n", + "\n", + " # setting matrix elements in A and B\n", + " Δ2 = Δ*Δ\n", + " Δ4 = Δ2*Δ2\n", + " αΔ = α*Δ\n", + " V_infl = -sin(ε) + exp(-(0.5*π-ε))\n", + " # μ[1] μ[2] μ[3] μ[4]\n", + " # ----- 0 --------- 1 --------- 2 --------- 3 -- i for μ = μ[i+1]\n", + " # ------|-----------|-----------|-----------|--- for μ = μ[i+1]\n", + " # 0---------- 1 --------- 2 --------- 3 -------- i for φ = φ[i]\n", + " # |-----------|-----------|-----------|--------- for φ = φ[i]\n", + "\n", + " # For φ EQ.\n", + " for i in 3:N-3\n", + " # println(\"$i\")\n", + " A_φφ[i,i+2] = complex(1.0,α*Re*V_bar[i,1]*Δ2/12.0)\n", + " A_φφ[i,i-2] = A_φφ[i,i+2]\n", + " A_φφ[i,i+1] = complex(-4.0-2.0*((αΔ)^2),-4.0*α*Re*V_bar[i,1]*Δ2/3.0)\n", + " A_φφ[i,i-1] = A_φφ[i,i+1]\n", + " A_φφ[i,i] = complex(6.0+4.0*((αΔ)^2)+(α^4-d2Vdz2[i,1])*Δ4,α*Re*V_bar[i,1]*(2.5+(αΔ)^2)*Δ2)\n", + " B_φφ[i,i+2] = complex(0.0,α*Re*Δ2/12.0)\n", + " B_φφ[i,i-2] = B_φφ[i,i+2]\n", + " B_φφ[i,i+1] = -complex(0.0,4.0*α*Re*Δ2/3.0)\n", + " B_φφ[i,i-1] = B_φφ[i,i+1]\n", + " B_φφ[i,i] = complex(0.0,α*Re*(2.5+(αΔ)^2)*Δ2)\n", + " end\n", + "\n", + " # Adjacent of boundaries\n", + " # For φ EQ. (at z=z_1)\n", + " A_φφ[1,3] = complex(1.0,α*Re*V_bar[1,1]*Δ2/12.0)\n", + " A_φφ[1,2] = complex(-4.0-2.0*((αΔ)^2),-4.0*α*Re*V_bar[1,1]*Δ2/3.0)\n", + " A_φφ[1,1] = complex(7.0+4.0*((αΔ)^2)+(α^4-d2Vdz2[1,1])*Δ4,α*Re*V_bar[1,1]*(31.0/12.0+(αΔ)^2)*Δ2)\n", + " B_φφ[1,3] = complex(0.0,α*Re*Δ2/12.0)\n", + " B_φφ[1,2] = -complex(0.0,4.0*α*Re*Δ2/3.0)\n", + " B_φφ[1,1] = complex(0.0,α*Re*(31.0/12.0+(αΔ)^2)*Δ2)\n", + "\n", + " # For φ EQ. (at z=z_2)\n", + " A_φφ[2,4] = complex(1.0,α*Re*V_bar[2,1]*Δ2/12.0)\n", + " A_φφ[2,3] = complex(-4.0-2.0*((αΔ)^2),-4.0*α*Re*V_bar[2,1]*Δ2/3.0)\n", + " A_φφ[2,1] = A_φφ[2,3]\n", + " A_φφ[2,2] = complex(6.0+4.0*((αΔ)^2)+(α^4-d2Vdz2[2,1])*Δ4,α*Re*V_bar[2,1]*(2.5+(αΔ)^2)*Δ2)\n", + " B_φφ[2,4] = complex(0.0,α*Re*Δ2/12.0)\n", + " B_φφ[2,3] = -complex(0.0,4.0*α*Re*Δ2/3.0)\n", + " B_φφ[2,1] = B_φφ[2,3]\n", + " B_φφ[2,2] = complex(0.0,α*Re*(2.5+(αΔ)^2)*Δ2)\n", + "\n", + " # For φ EQ. (at z=z_{N-2})\n", + " A_φφ[N-2,N-4] = complex(1.0,α*Re*V_bar[N-2,1]*Δ2/12.0)\n", + " A_φφ[N-2,N-1] = complex(-4.0-2.0*((αΔ)^2),-4.0*α*Re*V_bar[N-2,1]*Δ2/3.0)\n", + " A_φφ[N-2,N-3] = A_φφ[N-2,N-1]\n", + " A_φφ[N-2,N-2] = complex(6.0+4.0*((αΔ)^2)+(α^4-d2Vdz2[N-2,1])*Δ4,α*Re*V_bar[N-2,1]*(2.5+(αΔ)^2)*Δ2)\n", + " B_φφ[N-2,N-4] = complex(0.0,α*Re*Δ2/12.0)\n", + " B_φφ[N-2,N-1] = -complex(0.0,4.0*α*Re*Δ2/3.0)\n", + " B_φφ[N-2,N-3] = B_φφ[N-2,N-1]\n", + " B_φφ[N-2,N-2] = complex(0.0,α*Re*(2.5+(αΔ)^2)*Δ2)\n", + "\n", + " # For φ EQ. (at z=z_{N-1})\n", + " A_φφ[N-1,N-3] = complex(1.0,α*Re*V_bar[N-1,1]*Δ2/12.0)\n", + " A_φφ[N-1,N-2] = complex(-4.0-2.0*((αΔ)^2),-4.0*α*Re*V_bar[N-1,1]*Δ2/3.0)\n", + " A_φφ[N-1,N-1] = complex(5.0+4.0*((αΔ)^2)+(α^4-d2Vdz2[N-1,1])*Δ4,α*Re*V_bar[N-1,1]*(29.0/12.0+(αΔ)^2)*Δ2)\n", + " B_φφ[N-1,N-3] = complex(0.0,α*Re*Δ2/12.0)\n", + " B_φφ[N-1,N-2] = -complex(0.0,4.0*α*Re*Δ2/3.0)\n", + " B_φφ[N-1,N-1] = complex(0.0,α*Re*(29.0/12.0+(αΔ)^2)*Δ2)\n", + "\n", + " # assigned partial matrices A^ and B^ to A and B\n", + " A[1:N-1,1:N-1] = A_φφ[1:N-1,1:N-1]\n", + " # non-diagonal parts in B are zero\n", + " B[1:N-1,1:N-1] = B_φφ[1:N-1,1:N-1]\n", + "\n", + " A_J = transpose(A) # mathematical A -> Julia A_J\n", + " B_J = transpose(B) # mathematical A -> Julia A_J\n", + " C_J = inv(B_J) * A_J\n", + " D_J = eigvals(C_J)\n", + " F_J = imag.(D_J)\n", + " E_J = inv(B_J) * B_J\n", + " nc_max = findmax(F_J)[2]\n", + " c_out[i,j,k] = D_J[nc_max]\n", + "\n", + " end\n", + " end\n", + " end\n", + "\n", + "#println(D_J[findmax(F_J)[2]])\n", + "#println(F_J)\n", + "#println(V_infl)\n", + "\n", + " return c_out\n", + "\n", + "end # function OSE_set\n", + " \n", + "end # module Lilly1966_solver" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "ab18554c", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "20×20×20 Array{ComplexF64, 3}:\n", + "[:, :, 1] =\n", + " 0.864005+0.749363im 0.804944+0.14258im … -0.866025-0.0236016im\n", + " 0.864131-0.0178902im 0.803465-0.0170507im -0.869369-0.0216414im\n", + " 0.863543-0.0192996im 0.802856-0.0185836im -0.871508-0.0239824im\n", + " 0.863595-0.0210307im 0.802932-0.0203377im -0.873849-0.0264434im\n", + " 0.864056-0.0230905im 0.80343-0.022402im -0.876698-0.0286047im\n", + " 0.86475-0.0254921im 0.804155-0.024811im … -0.879753-0.0293886im\n", + " 0.865561-0.0282094im 0.804976-0.0275412im -0.881353-0.0290929im\n", + " 0.866423-0.031197im 0.805824-0.0305438im -0.881612-0.0292257im\n", + " 0.86731-0.0344051im 0.806668-0.0337645im -0.881312-0.0301184im\n", + " 0.868215-0.037786im 0.807505-0.0371523im -0.880885-0.0316704im\n", + " 0.869146-0.0412972im 0.808344-0.0406621im … -0.880543-0.0337225im\n", + " 0.870114-0.044902im 0.809203-0.0442561im -0.880385-0.0361257im\n", + " 0.871132-0.0485695im 0.810098-0.0479034im -0.880442-0.0387554im\n", + " 0.872215-0.0522742im 0.811049-0.051579im -0.880708-0.0415139im\n", + " 0.873373-0.0559955im 0.812071-0.0552636im -0.881154-0.0443304im\n", + " 0.874614-0.0597173im 0.813174-0.0589422im … -0.88174-0.0471579im\n", + " 0.875942-0.063427im 0.814369-0.0626036im -0.882427-0.0499689im\n", + " 0.877361-0.0671155im 0.815659-0.0662391im -0.883176-0.0527504im\n", + " 0.878868-0.0707762im 0.817046-0.0698423im -0.883957-0.0554991im\n", + " 0.880461-0.0744049im 0.818529-0.0734083im -0.884745-0.0582175im\n", + "\n", + "[:, :, 2] =\n", + " 0.864005+0.374679im 0.804944+0.0712867im … -0.866025-0.0118008im\n", + " 0.864108-0.00886654im 0.803427-0.00846178im -0.868106-0.012169im\n", + " 0.863672-0.00911107im 0.802997-0.00866241im -0.868558-0.0141932im\n", + " 0.864041-0.00974289im 0.803468-0.00925419im -0.868931-0.0162753im\n", + " 0.864661-0.0109122im 0.804182-0.0104381im -0.888474-0.0143946im\n", + " 0.865229-0.012477im 0.804803-0.0120424im … -0.884596-0.0134702im\n", + " 0.865663-0.0142773im 0.805247-0.0138821im -0.881409-0.0141596im\n", + " 0.865975-0.0162055im 0.805536-0.0158379im -0.879248-0.01568im\n", + " 0.866211-0.018201im 0.805726-0.0178463im -0.878098-0.0175758im\n", + " 0.86642-0.020234im 0.805872-0.0198787im -0.87777-0.0195709im\n", + " 0.866646-0.0222914im 0.806024-0.0219254im … -0.878043-0.0214972im\n", + " 0.866922-0.0243683im 0.806216-0.0239847im -0.87872-0.0232656im\n", + " 0.867275-0.0264624im 0.806477-0.0260575im -0.879647-0.0248446im\n", + " 0.867722-0.0285718im 0.806824-0.0281442im -0.880704-0.0262427im\n", + " 0.868275-0.0306939im 0.807272-0.0302445im -0.881808-0.0274918im\n", + " 0.868941-0.0328251im 0.807827-0.0323562im … -0.882909-0.0286316im\n", + " 0.869724-0.0349609im 0.808496-0.0344763im -0.883979-0.0296998im\n", + " 0.870626-0.0370963im 0.809279-0.0366008im -0.885004-0.0307277im\n", + " 0.871646-0.0392262im 0.810178-0.0387252im -0.885978-0.0317393im\n", + " 0.872779-0.0413453im 0.811189-0.0408451im -0.886901-0.0327518im\n", + "\n", + "[:, :, 3] =\n", + " 0.864005+0.249782im 0.804944+0.0475207im … -0.866025-0.00786722im\n", + " 0.864071-0.0058093im 0.803362-0.00554748im -0.86723-0.0086051im\n", + " 0.863893-0.00552721im 0.80326-0.00512189im -0.867312-0.00987253im\n", + " 0.864535-0.005966im 0.804068-0.00555373im -0.889371-0.0071025im\n", + " 0.865125-0.00702275im 0.804751-0.00667512im -0.884436-0.00706937im\n", + " 0.865472-0.00835174im 0.805121-0.00806954im … -0.881078-0.00944322im\n", + " 0.865619-0.0097694im 0.80525-0.00953105im -0.879438-0.0127762im\n", + " 0.865647-0.0112068im 0.805238-0.010989im -0.869346-0.0158563im\n", + " 0.865625-0.0126475im 0.805161-0.0124314im -0.869499-0.0170308im\n", + " 0.865601-0.0140943im 0.805074-0.0138658im -0.869737-0.0184915im\n", + " 0.865611-0.015555im 0.805013-0.0153043im … -0.870152-0.0201287im\n", + " 0.865679-0.0170364im 0.805005-0.0167569im -0.870746-0.0218832im\n", + " 0.865825-0.0185426im 0.805067-0.0182306im -0.884383-0.023496im\n", + " 0.866062-0.0200751im 0.805214-0.019729im -0.885082-0.0240333im\n", + " 0.866398-0.0216328im 0.805456-0.0212529im -0.88568-0.0244282im\n", + " 0.866841-0.0232127im 0.805798-0.0228007im … -0.886227-0.0247079im\n", + " 0.867394-0.0248107im 0.806246-0.0243693im -0.886772-0.0249094im\n", + " 0.868059-0.0264218im 0.806802-0.0259549im -0.887349-0.0250805im\n", + " 0.868836-0.0280409im 0.807466-0.0275526im -0.887962-0.0252651im\n", + " 0.869724-0.0296626im 0.808239-0.0291577im -0.888601-0.0254909im\n", + "\n", + "...\n", + "\n", + "[:, :, 18] =\n", + " 0.864005+0.0415953im 0.804945+0.00788083im … -0.866025-0.0013112im\n", + " 0.865363+0.000307466im 0.805041+0.000487881im -0.866193-0.00139859im\n", + " 0.865807-0.000668222im 0.805515-0.000650438im -0.866089-0.00144926im\n", + " 0.865697-0.00121719im 0.805398-0.00123774im -0.865966-0.00163391im\n", + " 0.865517-0.00163355im 0.805221-0.00167887im -0.865855-0.00186637im\n", + " 0.865279-0.00200464im 0.804978-0.00207877im … -0.865743-0.00212822im\n", + " 0.864979-0.00232907im 0.804654-0.00243169im -0.865635-0.00241775im\n", + " 0.864644-0.00259651im 0.804275-0.00270987im -0.865543-0.00273927im\n", + " 0.864308-0.00282044im 0.803892-0.00292029im -0.86548-0.00309905im\n", + " 0.863992-0.00302886im 0.803535-0.00309986im -0.865465-0.00350208im\n", + " 0.863704-0.00324625im 0.80321-0.00328123im … -0.865511-0.0039504im\n", + " 0.863448-0.00348909im 0.802918-0.00348478im -0.865635-0.00444267im\n", + " 0.863229-0.00376737im 0.802662-0.0037223im -0.865848-0.00497461im\n", + " 0.863054-0.00408644im 0.802444-0.00400023im -0.866161-0.00553977im\n", + " 0.862928-0.00444828im 0.802272-0.00432143im -0.866581-0.00613038im\n", + " 0.862861-0.00485242im 0.802151-0.00468617im … -0.867113-0.00673813im\n", + " 0.862857-0.00529656im 0.802088-0.00509282im -0.867757-0.00735477im\n", + " 0.862926-0.00577716im 0.80209-0.00553843im -0.868511-0.00797254im\n", + " 0.86307-0.00628988im 0.802162-0.0060191im -0.899011-0.00840608im\n", + " 0.863297-0.00682993im 0.802309-0.00653043im -0.899189-0.00837619im\n", + "\n", + "[:, :, 19] =\n", + " 0.864005+0.0394021im 0.804945+0.00746168im … -0.866025-0.00124219im\n", + " 0.865446+0.000309478im 0.805133+0.000466365im -0.86618-0.00131283im\n", + " 0.86581-0.000659494im 0.805517-0.000648517im -0.866073-0.00137005im\n", + " 0.865696-0.0011821im 0.805401-0.00120656im -0.865956-0.00154959im\n", + " 0.865521-0.0015857im 0.805232-0.00163509im -0.865847-0.0017729im\n", + " 0.865285-0.00195005im 0.804994-0.00203103im … -0.865734-0.00202376im\n", + " 0.864982-0.00226814im 0.804667-0.00238349im -0.865623-0.00230038im\n", + " 0.86464-0.00252382im 0.804277-0.00265398im -0.865526-0.00260718im\n", + " 0.864298-0.00273039im 0.803883-0.00284619im -0.865458-0.00295087im\n", + " 0.863978-0.00291906im 0.80352-0.00300359im -0.865434-0.0033369im\n", + " 0.863687-0.00311628im 0.803193-0.00316276im … -0.865472-0.00376773im\n", + " 0.863428-0.00333914im 0.802899-0.00334491im -0.865585-0.00424238im\n", + " 0.863205-0.00359777im 0.80264-0.0035618im -0.865785-0.00475683im\n", + " 0.863025-0.00389753im 0.802419-0.00381974im -0.866085-0.00530484im\n", + " 0.862893-0.00424041im 0.802241-0.00412146im -0.86649-0.0058788im\n", + " 0.862818-0.00462595im 0.802114-0.00446719im … -0.867006-0.0064705im\n", + " 0.862805-0.00505187im 0.802043-0.00485529im -0.867633-0.00707177im\n", + " 0.862863-0.00551464im 0.802035-0.00528277im -0.868372-0.0076749im\n", + " 0.862997-0.00600992im 0.802097-0.00574575im -0.899194-0.00814204im\n", + " 0.86321-0.00653294im 0.802233-0.00623979im -0.899368-0.00810884im\n", + "\n", + "[:, :, 20] =\n", + " 0.864005+0.0374281im 0.804945+0.00708422im … -0.866025-0.00118008im\n", + " 0.865518+0.000303147im 0.805211+0.00043872im -0.866167-0.00123657im\n", + " 0.865812-0.00065148im 0.805519-0.000646264im -0.866059-0.00129946im\n", + " 0.865697-0.00115061im 0.805405-0.0011784im -0.865948-0.00147389im\n", + " 0.865527-0.00154344im 0.805245-0.00159608im -0.86584-0.00168897im\n", + " 0.865294-0.00190303im 0.805015-0.00198983im … -0.865726-0.00192991im\n", + " 0.864988-0.0022172im 0.804686-0.0023453im -0.865612-0.00219474im\n", + " 0.864637-0.00246286im 0.804282-0.00261118im -0.865511-0.00248806im\n", + " 0.864288-0.00265298im 0.803876-0.0027857im -0.865437-0.00281696im\n", + " 0.863964-0.00282278im 0.803507-0.00292116im -0.865407-0.00318736im\n", + " 0.863671-0.00300092im 0.803177-0.00305886im … -0.865436-0.00360212im\n", + " 0.86341-0.00320506im 0.802881-0.00322069im -0.865539-0.00406055im\n", + " 0.863184-0.00344539im 0.80262-0.00341822im -0.865729-0.00455891im\n", + " 0.862999-0.00372724im 0.802396-0.00365752im -0.866016-0.00509114im\n", + " 0.862862-0.00405259im 0.802214-0.00394116im -0.866407-0.00564977im\n", + " 0.862779-0.00442095im 0.80208-0.0042693im … -0.866908-0.00622668im\n", + " 0.862759-0.00483006im 0.802002-0.00464025im -0.867521-0.00681379im\n", + " 0.862807-0.00527639im 0.801986-0.005051im -0.868243-0.00740343im\n", + " 0.86293-0.00575561im 0.802038-0.00549767im -0.899364-0.00789997im\n", + " 0.863132-0.00626296im 0.802163-0.00597579im -0.899533-0.00786392im" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using .Lilly1966_solver\n", + "\n", + "#-- Drawing parameters \n", + "na = 20 # wavenumber interval\n", + "ne = 20 # angle interval\n", + "nr = 20 # Reinolds number interval\n", + "α_min = 0.001\n", + "α_max = 1.0\n", + "ε_min = -60.0\n", + "ε_max = 60.0\n", + "Re_min = 25.0\n", + "Re_max = 500.0\n", + "\n", + "α = [α_min+(α_max-α_min)*((i-1.0)/(na-1.0)) for i in 1:na]\n", + "ε = [ε_min+(ε_max-ε_min)*((i-1.0)/(ne-1.0)) for i in 1:ne]\n", + "Re = [Re_min+(Re_max-Re_min)*((i-1.0)/(nr-1.0)) for i in 1:nr]\n", + "\n", + "c_full = full_EQ_set( α, ε, Re, 35 )\n", + "c_ose = OSE_set( α, ε, Re, 35 )\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "68e298c5", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[-1.5844293495493304e-6, 0.0015889913115860095, 0.0022927206982129647, 0.003136899295005263, 0.0033814079172513885, 0.0035205763915490346, 0.0035213784945906062, 0.0033751176553336987, 0.0031672749972534026, 0.0028984401502025153, 0.00257101712772535, 0.0022014837123090288, 0.0017979344200940435, 0.0013628946689964404, 0.000900019673950996, 0.0004135349173473162, -9.349615665233196e-5, -0.0006189920555579306, -0.001161144563903043, -0.0017182414228851005]\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmwAAAHTCAYAAACa4UQYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABmXklEQVR4nO3df3hT5d0/8HfaNOFn6Yql0KS12CJSQXDV+mNWZh0FpBUddsDDM7XTPZW1TBYpPzor7YOsHUKrDOczwYET3bT6lbKuFJA4hlimLZvwVJ9iFShJxYBiY0ESkp7vH7WRNGmbtPlxzsn7dV29vDi5c587x+Scdz73OScKQRAEEBEREZFohQV7AERERETUNwY2IiIiIpFjYCMiIiISOQY2IiIiIpFjYCMiIiISOQY2IiIiIpFjYCMiIiISOQY2IiIiIpFjYCMiIiISOQY2IpK0kpISKBQKKBQKbNmyxeXxlpYWx+M//OEPAz9AIiIfYGAjkplt27Y5AopCoYBarUZCQgIeeOABHD16NNjDw4MPPug0vhEjRuCaa67B8uXL8cUXXwyoz3vuuQcfffQRcnJynNYBAFdeeSU++ugj/OY3v/HZayAiCjQGNiKZeuutt/DRRx/h0KFDqKiowCeffIKbbrpJFKEtLi4OH330ET766CPs27cPOp0Or776Km655RZ88803Xvc3atQoXHPNNRg1ahQAoKysDB999BEAICIiAtdccw3GjRvn09dARBRIymAPgIj8IykpCYmJiQCA66+/HnfddRdGjx6NrVu3oqKiIqhj6w5R3W666SYkJyfjzjvvxN69e3H33XcPqv9x48YxoBGRrLDCRhQiwsPDAXSFpcv95S9/wXXXXYchQ4ZgwoQJeOKJJwZU5RospVLpMj6bzYaSkhJceeWVGDZsGK6//nps3769374efPBBR1glIpIDVtiIZO6bb75Bc3MzSkpKMGrUKCxevNjx2PPPPw+dTofS0lL88Ic/xMcff4wVK1bg8OHDqKmpCcj4zGYz/vnPf+LRRx/Fbbfdhh/96EeOx+6//34cOnQI5eXluPrqq7F//348/PDDOH36NJYtWxaQ8RERiQEDG5FMJScnAwDsdjtiYmLwy1/+Elu2bMEVV1wBADh//jyWL1+Op59+Gg8//DAAIDU1FVdddRVuuukmvPfee0hLS3Ppd+fOnfjxj3/c57p37tyJu+66q9fHT5486aio2e12TJ06Fb/61a/w4IMPOips77zzDv7yl7/g//7v/3D11Vc7xmez2bBmzRo8+uijLtVCIiK5YmAjkqna2lpERUXh7bffxpo1a3Dy5ElHWAOA+vp6tLe345FHHsEjjzzi8vx//etfbgPbHXfcgX//+999rvvKK6/s8/G4uDjs3r0bX3zxBf785z9j69atCA8PdwpgdXV1EAQBKSkpTs8VBAGdnZ349NNPMXHixD7XQ0QkFwxsRDJ19dVXIzExEWlpabj22muRnZ2N66+/Hr/4xS8AAJ9//jkA4LXXXnNUsC7X20n7I0eOxOTJkwc1toiICEcf06dPx4gRI7B48WJ8//vfx7Rp0xzji4yMxMGDB9320V8oJCKSEwY2ohCQlZWF7OxsPPHEE/jpT3+KkSNHYvTo0QC6KlbeBDBfTIn2VFpaipdffhnLly/Hnj17AACjR4+G2WzGuHHjHGMlIgpVvEqUKER035j22WefBQDceuutGDlyJCorK9HZ2elo19nZiYceeggtLS1u++meEu3rLz093auxDR8+HPn5+di7dy8OHToEAJg5cyYAYMOGDU5tP//8c/z0pz8NypWsRETBwsBGFCJuu+02pKam4umnn4bFYkFkZCR++9vf4uDBg8jOzsZbb72FAwcO4J577sHevXsRHR3ttp/uKdG+/kaOHOn1+B555BEMGTIE5eXlALqC4U9+8hOUl5fjscceQ319Paqrq3HHHXfg3LlzGDp06KC2BxGRlDCwEYWQpUuX4vPPP8ef/vQnAMDixYvx+uuv48yZM8jKysK9996LkSNHor6+vtfA5i9XXHEFFi1ahJ07d6K5uRkA8PLLL+PJJ5/Ezp07MX36dCxZsgT33HMP3njjjYCOjYgo2BSCIAjBHgQR0UCVlJTgxIkT2LZtW5/ttm3bhm3btuHvf/97QMZFRORLrLARERERiRwDGxFJ3p/+9CcolUr88Y9/dHnsk08+gVKpxEMPPRSEkRER+Yaop0SPHz+OpUuXoqGhAZ2dnbj99tvx9NNPO+4PZbFYsHr1arz++uu4cOECbrzxRjz33HOIi4sL8siJKFBMJhNMJhMAQKPR4Hvf+57T41arFceOHQPQdTXq+PHjAz5GIqLBEm2F7auvvsIdd9yB7OxsGAwGHD9+HGq1Gr/73e8cbfLz81FfX4+Ghga0trYiKSkJs2fPht1uD+LIiSiQxowZ47g6tWdYAwCVSuV4nGGNiKRKtBW21atX4/Dhw/jrX//qWNbZ2YmwsK6M2draivHjx+O9995DamoqgK5v0nFxcdi6dSuys7ODMm4iIiIiXxNthc3dndK7wxoA7N+/H7GxsY6wBnR9k87MzMSuXbsCNk4iIiIifxPtT1N9/PHHiIqKws9//nO89dZbGDFiBObPn4+VK1dCqVTCaDS6PVdNo9E4zldxx2KxwGKxOP7d2dmJL7/8EqNHj4ZCofDLayEiIvIXQRDw9ddfIy4uzqmwESyNjY3Y+edMlDx1lsdVHxJtYLPb7XjyySfx+9//Hs8//zyOHTuGefPm4dy5c9iwYQMiIiLcvjEVCgX6muUtKytDaWmpP4dOREQUcKdOnYJWqw3qGARBwC9+/gO8/28LUm58DfPnzw/qeOREtIEtISEB//Vf/4Xp06cDACZOnIji4mIUFBRgw4YN0Gq1aGtrc3leW1sbNBpNr/2uWrUKOp3O8e/29nYkJCRg/G+KETZkiOsThtt67Sv2CrMXryi0qIQw/Jf1ejyv+hesis7+n0Ae47b1n57b9vOzkcEeUpfzwdlVKy/4rjqiViiwcvyVKD9+EpZ+Tp1WdgSuKjPkbMBW5Tdhdgve31oyoJ+E87XXXnsNhs9seObJGKxa/p+45557oFargz0sWRBtYEtPT3eauuzW/T8+IyMDJpMJR44cwXXXXQcAsNls0Ov1eO6553rtV61Wu33z2NQqhA1RuT5haO/l5c4Rot18QdcphGPYxWHoHBKBTgWv2vUlblv/uXzbGk0jATH8XGmHEnDzXTIgbL6bXgsPU2DYsGGwqdWwdfYe2JQdCtjc7Ir9pTP4M4g+E+zpR4vFglXL/xOly0fjgZ9EYvP2djxbroVu9ZmgjksuRPtWXblyJZ555hns378fAHDy5En893//N372s58BAGJiYpCbmwudTgez2Qy73Y6ioiJER0djzpw5wRw6ERFRyHm2XIsRw8PwwE8iER6uQPnjV+DJp7/EuXPngj00WRBtYEtOTsYrr7yC5cuXY8yYMcjIyMCCBQvwxBNPONps3LgRU6ZMQUpKCrRaLZqbm1FXVwelkpUvIhq406ZRwR5Cl47g7cuUHaI9PJAInTt3Dk8+/SXKH78C4eFdlb6ZdwzD969T4zePJwV5dPIg6mQzffp0/POf/+z1cbVajcrKSlRWVgZwVERE5A+BPHcNAIadDujqZO03jyfh+1PUmHnHMMcyhUKB3z5+BdLvNqBg+UlceeWVQRyh9PErFBGRGAWxukbkjZMnT+LZP36F3xZf4XIe3fVThmBe1ggU6aYEaXTywcBGRPQt0UyFBhmnQ8kbRbopmJc1AtdPcX91zJoVo/FmbQcaGxsDPDJ54aeSiEhsQrC6FujpUPKNw4cP483aDqxZMbrXNgnaCCx5OAqFj6b3eZ9U6hsD2yB89nlUsIdARD7Cz3OXUKmu8fy1wRMEAYWP3oYlD0chQRvRZ9uVS76HIx9a+NORgxAan0wiIqkIweoaSdOuXbvwQZMFK5d8r9+2oyLDUawbjeW6e2Gz9X5DeuodAxsRhTxW14i8Y7PZsFx3Lx7/VTRGRYZ79Jy8+0fhokXAH5/p/deIqHcMbEQU0kQV1oJcXQvWdCjPX5OerRs1uGgR8MgDUR4/R6VSYG3RaJQ89QXOnz/vv8HJFAMbERGFFJ6/Njjnz5/H6nVfYG3RaKhU3oXt+7JG4EptBDaUJgxo3du2bcPkyZOh1WqRlpaGgwcP9trWaDRi/vz5SExMhEajgU6ng9VqdWpz6NAhpKenIyEhARMmTMDmzZudHj98+DAiIiKg1Wqd/t58801HG4vFgpUrVyI5ORlxcXGYO3eu2986HywGNiIKWayuEXlvQ2kCrtRG4L6sEV4/V6FQYN0TV+Cp35/D559/7tVzt2/fjlWrVqGqqgoGgwGFhYWYM2cOjh8/7tLWarVixowZ0Gq1aGlpQVNTExobG6HT6RxtmpubkZmZiaVLl6K1tRXV1dUoLi7G66+/7mhjMBjw/e9/HwaDwenv3nvvdbTJz89HfX09Ghoa0NraiqSkJMyePRt2u29/65mBjYiIOB1KHvn888/x1O/PYd0TrjfJ9dQP0obiR7cPw+rlE716XmlpKZYtW4ZJkyYBAHJycpCeno5Nmza5tK2qqoLJZEJZWRmUSiWioqJQUVGBLVu24OzZswCA9evXY/r06Zg3bx4AICUlBYWFhSgrK3P0YzQaER8f3+uYWltbsXXrVlRUVCAqKgpKpRLl5eUwGo2ora316vX1h4GNiEKSqKprRBKxevlE/Oj2YfhB2tBB9VP26yvw4mtmHDlyBGaz2enPYrG4tD916hRaWlqQnZ3ttDw7O9vtrUL0ej0yMzOhUqkcy1JTUxEdHQ29Xu9o466/w4cPw2QyAeiqsCUk9D59u3//fsTGxiI1NdWxTKVSITMz0+e3MGENnogo2DgdGjA8f21wqus68IvK8dhxfvjgOho7HFN+cA6PPfYY3nrrLaeHVq9ejZKSEqdlRqMRABAXF+e0XKPROB7r2X7y5Mkuyy9vbzQa3fbX/diYMWNgNBqhUChw77334oMPPsDo0aPxyCOP4KGHHuq1j+5+jh071tcW8Br3EkQUclhdc8bpUPJYWAzUQ33zflENDUPGrRl44403nJar1WqXthERXTfmDQtzXrdCoXD76wkREREubXu2d9eme5q3u41CoYDJZMKmTZuQmJiIhoYGzJ07FzabDXl5eR6tx1cY2IgopIgurLG6RiFMqVQiMjKy33ZarRYA0NbWhuTkZMfytrY2R1WsZ3t3V2pe3t5dm+5/d7fZunWr0+M33ngjHn30UWzduhV5eXkercdXeA4bEVEIC5WfoiJpi42NxdSpU11O5N+9ezdmzZrl0n7mzJnYu3ev068qNDU14cyZM8jIyHC0cdfftGnTEBsbCwBuq2R2u91RicvIyIDJZMKRI0ccj9tsNuj1erfjGgx+UokoZLC6Jh7BmA7l+WvStmLFCqxbt85xbtiOHTuwZ88eFBQUuLTNyspCTEwMiouLYbfb0d7ejiVLliA3NxcxMTEAgIKCAuzbtw87d+4E0HWbj7Vr12LFihWOfrKzs7Fs2TJcuHABANDQ0IBnnnkGP//5zwEAMTExyM3NhU6ng9lsht1uR1FREaKjozFnzhyfvn4GtkES3QGAiIhIhhYuXIji4mJkZWUhLi4Oa9euRU1NDZKSkmAwGKDValFVVQWga6q1rq4OH374IeLj43Httddi6tSpeOaZZxz9JScno6amBmvWrIFGo0FWVhZKSkqwYMECR5s//OEPOHPmDCZOnIjY2Fj8x3/8B5544gn87Gc/c7TZuHEjpkyZgpSUFGi1WjQ3N6Ourg5KpW+/kIXu1zsiCimi+3Ilguoap0NJavLy8pCXl+eyXKvVwmAwuCyrrq7us7/09HS8//77vT6u0Wjw4osv9tmHWq1GZWUlKisr+2w3WPy0EpHsiS6shTheHUrkPQY2IqJAY3Ut4Hj+GkldaH1iiSjksLpGRHLAwNYfEXwTJqKBEWVYC/F9CqdDiQaGgY2IKMSE2nQokRzwU0tEssTqGnXj+WskBwxsRCQ7ogxrIsHqGpE08ZNLRBQIrK4R0SAwsPkAv80TiQc/j3Q5ToeSXDCwEZFsiDasiaS6xulQIunip5eIiPyOt/MgGhwGNiKSBVbX+sbqGpG08RNMRJIn2rBGQcXz10hOGNiIiPyF1bVv18/pUKLBYmAjIkljdY2IQgEDGxFJlqjDmkiqa6GK06EkNwxsRCRJog5rIsLpUCJ5YGAjIvI1VteIyMcY2HyE3/aJAoefN88Eu7pGRL7DTzMRSYrowxqraw7Bmg7l+WskRwxsREQyxOoakbzwE01EksHqmnTwYgMi32Jg8wR3wkRBJ/qwRqLA6VCSKwY2IiJfENEXO06HEskPP9VEJHqir66JKKyJAadDiXyPgY2IRE30YU1kQrm6xulQkrPQ/WT7AQ8sRL4lic8Uq2tEFAAMbEREMiGG6hqnQ4n8I/ifbiIiN1hdIyL6DgMbEYmOJMKayIR6dY3nr5HcBf8TTkQkRayuEVEAMbARkahIorrGsEZEAcbARkSiIYmwJkKcDg3aqokCJvifciIiSCissbpGFDTbtm3D5MmTodVqkZaWhoMHD/ba1mg0Yv78+UhMTIRGo4FOp4PVanVqc+jQIaSnpyMhIQETJkzA5s2bXfp54YUXcO2110Kj0WDSpEl4/vnnnR6vqKjAiBEjoNVqnf5On/btNwkGNh+TzEGHSET4uRk4MVTXiAJh+/btWLVqFaqqqmAwGFBYWIg5c+bg+PHjLm2tVitmzJgBrVaLlpYWNDU1obGxETqdztGmubkZmZmZWLp0KVpbW1FdXY3i4mK8/vrrjjYvvfQSSktLUVVVBaPRiDfffBMlJSX485//7GhjMBiQn58Pg8Hg9Dd27Fifvn7JfNINBgOio6Px4IMPOpZZLBasXLkSycnJiIuLw9y5c9HW1ha8QRKR1yQV1lhdc4vToRQIpaWlWLZsGSZNmgQAyMnJQXp6OjZt2uTStqqqCiaTCWVlZVAqlYiKikJFRQW2bNmCs2fPAgDWr1+P6dOnY968eQCAlJQUFBYWoqyszNHPe++9h/Xr1yMlJQUAcM0112DhwoWoqqpytDEajYiPj/fb6+4micAmCAIeeOABaLVap+X5+fmor69HQ0MDWltbkZSUhNmzZ8NutwdppEQkWyIMa6yukdTZbDaYzWanP4vF4tLu1KlTaGlpQXZ2ttPy7Oxs7Nq1y6W9Xq9HZmYmVCqVY1lqaiqio6Oh1+sdbdz1d/jwYZhMJgDA7373O/zkJz9xanP06FFERkY6/m0wGJCQkODlK/ee+PZAbmzYsAERERH48Y9/jBMnTgAAWltbsXXrVrz33nuIiooCAJSXlyMuLg61tbUu/xOISHwkVV0jt/jLBqHnJcPNGDF83KD7+bhjJ4bq9SgqKnJavnr1apSUlDgtMxqNAIC4uDin5RqNxvFYz/aTJ092WX55e6PR6La/7sfGjBnj9NilS5eg0+lQX1+P+vp6p3UdPnwYTz31FAwGA5KSklBaWoof/OAHfb18r4n+69kHH3yA8vJy/P73v3davn//fsTGxiI1NdWxTKVSITMz023aHjQRfrsmkjJJhTURfv5ZXSM5yMjIQHt7u9PfqlWrXNpFREQAAMLCnN/3CoUCgiC4bd+zbc/27tooFF1fQHr2efLkSaSnp2Pfvn145513nMKgSqXCN998g+rqarS0tCA3NxczZszAkSNHPNkEHhPfXugyFy9exKJFi1BeXo6rrrrK6TF3yRjoSsfHjh3rtU+LxeJUbjWbzQAAlSIM4Yr+doCe7SBVQrhH7eRMJYQ5/Zd8Rw7b9rRpFNTBHoQb6m8/4+qen/V+9w2BpwwTR2XL03Gov22n9uG4h5og8qNYYITZxfFeGAilUuk0vdib7lOi2trakJyc7Fje1tbmqIr1bO/unPbL27tr0/3vy/tsbGzE7Nmzcf/992Pt2rVQq533Xj0zx6JFi7B9+3a88soruO666/p9bZ4S9Vt9+fLlSEpKwsMPP+zymCfp2Z2ysjKUlpa6LF8zdjyGDRs2uAF3u+ibbuSgwHJDsIcgW5Letv3vn4NqTeQE5wUiH6+U/Pq6xGAPQXYuXLiA+j8EexT+FRsbi6lTp6K2tha//OUvHct3796NWbNmubSfOXMm8vLyYLPZoFR2RZ2mpiacOXMGGRkZjja1tbVOGWP37t2YNm0aYmNjAXRV1u666y48++yzyMnJcTu2zs5Olzxit9sd1TpfEW1g27NnD1599VUcPXrU7eOepGd3Vq1a5XRZr9lsRnx8PIpPH0f40CF9D2q4zbPBAxg7pt3jtnKkEsJQYLkBm9QNsCo6gz0cWZHytj1tGhXsIfRJjTCsiZyAYvPHsKATOC/OXaTyvDgqKt6cv6YOU+DX1yVi7ZETsHT2/qXaG0NNPulG8sLsrifpy9GKFStQWFiIWbNm4eqrr8aOHTuwZ88eHD582KVtVlYWYmJiUFxcjCeffBIdHR1YsmQJcnNzERMTAwAoKChAamoqdu7cibvvvhvNzc1Yu3YtKioqHP0sXrwYv/jFL3oNa1999RVSU1Px3//931i4cCEUCgX+9Kc/4cCBA3jmmWd8+vrFuTcCUFtbC5PJ5Ei5l3vxxRfx2muvwWQy4ciRI46So81mg16vx3PPPddrv2q12qWcCQBWoRNhQn8HP88PjlYFr1QFAKuik9vCT6S2bbvOWZNGwLSgsyuw9btPCA57Z/CnaJUdCtjhffCydAo+CWzDTgPW/puFhDAfBWCxW7hwIcxmM7KystDR0QGNRoOamhokJSXBYDDg5ptvRmVlJXJycqBUKlFXV4f8/HzEx8cjLCwMOTk5KC8vd/SXnJyMmpoa6HQ6LF68GMOGDUNJSQkWLFjgaLNr1y40Nja6vaGuwWBAVFQUXn75ZTzxxBNYvnw5LBYLJkyYgNraWsftR3xFtIHt6aefxtNPP+20rKSkBCdOnMC2bdsAAG+99RZ0Oh3+3//7fxg+fDiKiooQHR2NOXPmBH7ARNQrSV1g0E2EFxoAvNiAQlteXh7y8vJclmu1WhgMBpdl1dXVffaXnp6O999/v9fH+zrFqtvNN9+MPXv29NtusCT9yd+4cSOmTJmClJQUaLVaNDc3o66uzjFfTURERCQHkko2Pe/LolarUVlZicrKyuAMiIj6Jcnq2nklxDh9K5bqWrDvvcZfN6BQJI5PPxHJkiTDGhGRCDGwEZFfMKz5FqtrRKFNHHsAGeLBikKZZN//Ir2NB32H06EUqhjYiMinJBvWRHpVKCCe6hoRBQ/3At4Q8Q6dSAwkG9bII5wOJQoeBjYi8glJhzURfxljde07nA6lUMY9ARER9YvVNaLgYmDzI0lXHIi8IOn3OqtrRCQB3BsQ0aAwrFEgcDqUQh0DGxENmKTDmsiJqbrG6VCi4BPPHkGmeEAjuZL8e5vVNclgdY2IgY2IBkDyYU3kWF0jop7Es1cgIkmQRVhjdY2IJIaBjYg8xrDmf2KqrokBp0OJunDPEACyOMhRyOP7OPRwOpRIPBjYiKhfsglrrK4RkURx7+Atke/wiXyNYS00iaG6xulQou8wsAWIbA56FFL4vg0cVteIqC/8yklELmQX1EReXWNYI6L+cC9BRE4Y1ojToUTiw8BGRA4Ma4HH6hoReYJ7igCS3cGQZIXvTwLEUV0jIlcMbEQkz7DG6ppkcTqUyBX3FkQhjmGNurG6RiReDGwBJsuDI0mWLN+PEglrrK4RkTe4xyAKUbIMayR5nA4lco+BbSAk8g2eqDeyDWsS+WyKsbrG6VAicRPfXoOI/IphjYhIehjYgkC2B0wSPdm+9yQU1lhd6x2nQ4l6J529HBENmGyDGhFRiBDfVz0i8inZhzVW1wZFLNU1Ik9s27YNkydPhlarRVpaGg4ePNhrW6PRiPnz5yMxMREajQY6nQ5Wq9WpzaFDh5Ceno6EhARMmDABmzdvdnrcbrfj4MGDWLp0KaKjo7Ft2zaX9VgsFqxcuRLJycmIi4vD3Llz0dbW5pPXeznx7T1ChOwPoiQKsn+fMazJBqdDqT/bt2/HqlWrUFVVBYPBgMLCQsyZMwfHjx93aWu1WjFjxgxotVq0tLSgqakJjY2N0Ol0jjbNzc3IzMzE0qVL0draiurqahQXF+P11193tNm2bRt+9atfYcSIEQgPD3c7rvz8fNTX16OhoQGtra1ISkrC7NmzYbfbffr6uQchkinZhzUiCimlpaVYtmwZJk2aBADIyclBeno6Nm3a5NK2qqoKJpMJZWVlUCqViIqKQkVFBbZs2YKzZ88CANavX4/p06dj3rx5AICUlBQUFhairKzM0c9DDz2E9957D08++SSGDx/usp7W1lZs3boVFRUViIqKglKpRHl5OYxGI2pra336+hnYiGQoJMIaq2uDxulQkopTp06hpaUF2dnZTsuzs7Oxa9cul/Z6vR6ZmZlQqVSOZampqYiOjoZer3e0cdff4cOHYTKZPBrX/v37ERsbi9TUVMcylUqFzMxMt+MaDHHuRaTABweLkDioUkB99nlUaLyvJBTWiKh3NpsNZrPZ6c9isbi0MxqNAIC4uDin5RqNxvFYz/Y92/Zs766NRqNxWl9/PFmPr3CvRyQTIRHUAMmFNVbX+sfz16Tl7Jcj8LUP9jffXIyAXq9HUVGR0/LVq1ejpKTEaVlERAQAICzM+fOkUCggCIJL3xERES5te7Z310ah6PpcuOvTHU/W4yvS2vMRkVsMa0QkRRkZGXjjjTeclqnVapd2Wq0WANDW1obk5GTH8ra2NkdVrGd7d1dqXt7eXZvuf7vr0x1P1uMr4vzqF0JC5kBLfhEyU6ASxeoaUd+USiUiIyOd/twFttjYWEydOtXlRP7du3dj1qxZLu1nzpyJvXv3wmazOZY1NTXhzJkzyMjIcLRx19+0adMQGxvr0fgzMjJgMplw5MgRxzKbzQa9Xu92XIMhzr0JEfXrtGlUsIcQWKyuEYW0FStWYN26dTh27BgAYMeOHdizZw8KCgpc2mZlZSEmJgbFxcWw2+1ob2/HkiVLkJubi5iYGABAQUEB9u3bh507dwLous3H2rVrsWLFCo/HFBMTg9zcXOh0OpjNZtjtdhQVFSE6Ohpz5szxwav+DgMbEYmfBMOaWKtrYsPz18hTCxcuRHFxMbKyshAXF4e1a9eipqYGSUlJMBgM0Gq1qKqqAtBVuaurq8OHH36I+Ph4XHvttZg6dSqeeeYZR3/JycmoqanBmjVroNFokJWVhZKSEixYsMCrcW3cuBFTpkxBSkoKtFotmpubUVdXB6XSt/st6e0FiULYZ59HQY0wIDLYIwkghjWf4nQoSVleXh7y8vJclmu1WhgMBpdl1dXVffaXnp6O999/36N1nzhxwu1ytVqNyspKVFZWetTPQIl3rxJCeA4SeYLvEyKi0MXANhgS/OZP0hSyYU2CnzFW1zzH6VAiz0lvb0gUQkI2qAGSDGtERP4i3q+CISakD8zkVki/JyQa1lhdIyJ/Ee/ehSiEhXRYIyIiFwxsRCIT8mGN1TWfE2N1jeevEXlHvHuYEBTyB+oQx18tAMMaEVEvuJchEoGQD2qAZMMaEVEgcA9JFEQMat+ScFgTe3WN06FE8iDuPY0U+PhAwwN46OD/629JOKwREQUK95REAcagdhmJhzVW14goUMS9tyGSGYa1yzCsERF5jHscEeJBXZ74/5UCSazVtaGmYI+ASJpEHdheeOEFXHvttdBoNJg0aRKef/55p8ctFgtWrlyJ5ORkxMXFYe7cuWhrawvSaInc4+063GB1jYjIK6Ld67z00ksoLS1FVVUVjEYj3nzzTZSUlODPf/6zo01+fj7q6+vR0NCA1tZWJCUlYfbs2bDb7UEcOdF3GNTcYFgjIvKaaPc87733HtavX4+UlBQAwDXXXIOFCxeiqqoKANDa2oqtW7eioqICUVFRUCqVKC8vh9FoRG1tbWAH64cDEA/00saqWi8kHtakQqzToUQ0cKLde/7ud79zWXb06FHExcUBAPbv34/Y2FikpqY6HlepVMjMzMSuXbuQnZ3ttl+LxQKLxeL4t9ls7nquIgzhisHkV99nX5UQ7vM+A0UlhDn9N5ScNo2C2o/9q799r6nF+33LvfNKQOQ5Qq1QXPZf1+2rPK8Q8dfc7yjDxLeh1d+OSRUuvrFJXZid2zQUiDawXe7SpUvQ6XSor69HfX09AMBoNDrC2+U0Gg2OHTvWa19lZWUoLS11Wb5m7HgMGzbMd4P2hYvBHsDgFVhuCPYQAi8yMKtZEzkhMCvylQBtF19YM/aqYA9Bth67MzHYQ5CdCxcuoP4PwR4F+ZvoA9vJkycxf/58mM1mvPPOO5g8eTIAICIiAmFhrl91FQoFBEHotb9Vq1ZBp9M5/m02mxEfH4/i08cRPnTIwAc63Dbw5/Zh7Jh2v/TrbyohDAWWG7BJ3QCrojPYw/G706ZRAVuXGmFYEzkBxeaPYYEEtu150e9mHNQKBdaMvQrFpz+Fpcd+RHleGlUMsU6HqsMU+PV1idiw7wSs9t730eS9MLul/0YkeaLekzY2NmL27Nm4//77sXbtWqjV3000abVat1eEtrW1QaPR9NqnWq126qebVehEmDCYg59/DpxWhbQvoLAqOiX/GvrTda5a4IOTBZ3iD2wdSgRj2wxc15dAiyDActn+QNkRBjukETLsIt/cVrsAq00a21Iqwjq5PUOBaM/GOHnyJO666y48++yzWL9+vUvIysjIgMlkwpEjRxzLbDYb9Ho9Zs2aFejh8mTqEMQLC/rBz0TAibW6RkSDJ9rAtnjxYvziF79ATk6O28djYmKQm5sLnU4Hs9kMu92OoqIiREdHY86cOQEerf8wEIgT/7/0Q0ZhjbfxICIxEO1eddeuXWhsbMTmzZtdHjMYDACAjRs3YuXKlUhJSYHdbkdaWhrq6uqgVIr2ZZHEMaiRWIm9usZfOCAaHNEmm74uHOimVqtRWVmJysrKAIwoeD77PArjYr8K9jBCHsOah1hdIyLyOfnsWYn8hEHNCwxrRER+wT0SUR8Y1rwgo7AmNWKfDiWiwWNg8yU/HrAYHAKLV4B6SWZhTSr3XJOKYaeDPQIi6WNgI+qBQc1LMgtrUsPqGlFo4J6W6FsMagPAsEZEFBCssEkIA4X/cNsOgAzDGqdCiUisGNgo5DGsDYAMw5oUcTqUKHRwr+trHUpghH9+CJ58i0FtgGQa1pQdYfwK6we84IDIN7h7khiGDN/gdhwgOYc1iWF1jSi0SG8vRTQIvF3HIMg0rBGRdGzbtg2TJ0+GVqtFWloaDh482Gtbo9GI+fPnIzExERqNBjqdDlar1anNoUOHkJ6ejoSEBEyYMMHtz2H2t86KigqMGDECWq3W6e/0ad+WlxnYKGQwqA2CjMOaFKtrRKFo+/btWLVqFaqqqmAwGFBYWIg5c+bg+PHjLm2tVitmzJgBrVaLlpYWNDU1obGxETqdztGmubkZmZmZWLp0KVpbW1FdXY3i4mK8/vrrXq3TYDAgPz8fBoPB6W/s2LE+ff3cU0kQg4f3uM0GQcZhTao4HUqhqLS0FMuWLcOkSZMAADk5OUhPT8emTZtc2lZVVcFkMqGsrAxKpRJRUVGoqKjAli1bcPbsWQDA+vXrMX36dMybNw8AkJKSgsLCQpSVlXm1TqPRiPj4eL+97m4MbP4QgAMcA4hnOAU6SDIPa6yu+RcvOKD+2Gw2mM1mpz+LxeLS7tSpU2hpaUF2drbT8uzsbOzatculvV6vR2ZmJlQqlWNZamoqoqOjodfrHW3c9Xf48GGYTCaP12kwGJCQkOD9i/cS91YSxiDSN26fQWJYIyJ3Lii79g+D/bOFQa/XY9SoUU5/l1e4uhmNRgBAXFyc03KNRuN4rGf7nm17tnfXRqPROB7zdJ1GoxGHDx9Geno6xo8fjx/96Ed9nls3UPLeI4eA7lAyLvaroI5DTBjUfEDmYU3KOB1KcpKRkYE33njDaZlarXZpFxERAQAIC3P+oqVQKCAIgtv2Pdv2bO+ujULR9fkSBMHjdapUKnzzzTeorq7GqFGj8Je//AUzZszAoUOHcN1117l/4QPAvbJMfPZ5FEMbGNZ8IgTCGqtrROKgVCoRGRnZbzutVgsAaGtrQ3JysmN5W1uboyrWs31bW5vL8svbu2vT/e/L++xvnceOHXPqY9GiRdi+fTteeeUVnwY27rX8JQgHvVAOKzxXzUdCIKxJmZSqazx/jXwpNjYWU6dORW1trdPy3bt3Y9asWS7tZ86cib1798Jm++5G9k1NTThz5gwyMjIcbdz1N23aNMTGxnq8zs7OTpf12+12R7XOVxjYZCYUQ0sovma/CJGwxuoakTStWLEC69atc1S0duzYgT179qCgoMClbVZWFmJiYlBcXAy73Y729nYsWbIEubm5iImJAQAUFBRg37592LlzJ4Cu23ysXbsWK1as8HidX331FSZMmICXX34ZnZ2dEAQBL774Ig4cOID777/fp68/NPbQISZUpkcZ1HwoRMIaEUnXwoULYTabkZWVhY6ODmg0GtTU1CApKQkGgwE333wzKisrkZOTA6VSibq6OuTn5yM+Ph5hYWHIyclBeXm5o7/k5GTU1NRAp9Nh8eLFGDZsGEpKSrBgwQKP1gkAUVFRePnll/HEE09g+fLlsFgsmDBhAmprax23AvEV7qVlSs6hjUHNx0IorEm5uial6VAif8nLy0NeXp7Lcq1WC4PB4LKsurq6z/7S09Px/vvvD2id3W6++Wbs2bOnzz58IXT21MEQ5B+Cl9sVpAxqfsCwRkQkCdyDhQA5BB05vAbRYVgjP+EFB0S+Fzp77BAn1SlSBjU/CKGgJhecDiUifu30NxEdHKUUfnibDj8R0fsxUFhdIyI54J4sxIg9BDGo+RHDGhGRZHFvFoLEGojEOi5ZCMGwJhdSmw7l+WtE/sHAFggiPFiKKRyxquZnInz/BQKra0QkJ6G5JycAwb/tB0NaADCsSZrUqmtE5D/y2KvRoAQ6OLGiFiAhGtYoeDgdGhxXNJ4L9hAoALhHD5Qg30S3P/667cdp0yhY4PrDuORnIRzW5FJdI/LEmPoveCQPEfzfTA79hTZvqmJqhAGRgx8TDQDDmixwOpT6M6b+i2APgQIodPfs5BanKiUuhMMaBRenQ4n8Sz5fR6WAB1Pylw5lyL+/WF2jUMLqWuiRzx6OKFSFeFAD5BXWiPrDsBaauJcjkjKGNRIBTocGDsNa6GJgCzQeYMlXzvO9BMivusbpUOoNw1pok9eejoiIAorVNaLAYGAjkhpW1hxYXaNQweoayWtvJxWcFqWB4JWgTuQW1oh6w7BGAAMbkTQwqJEIcTrU/xjWqBsDG5HYMay5kGN1jdOh1BPDGl1Ofns9qeBBmDzB94kLOYY1IqL+cM9HJFYMayFDitU1Tof6F6tr1BMDG5EYMay5xeoahQKGNXKHe79g4kGZeuKVoL1iWBMPVtf8h2GNesM9IJFYMKiFJClOh5J/MKxRXxjYiMSAYa1PrK6R3DGsUX+4Fww2HqiJ74E+MayJC6dDiYKDe0KiYGJYC2mcDiWA1TXyDAMbUbAwrPWL1TVxYXXN9xjWyFPcG4oBD9yhhVeCElhdI4Y18g4DG1EgMah5jNU1cWF1zbcY1shb3CMSBQrDGn2L1TUi8pZHgU2v1/t7HIOybds2TJ48GVqtFmlpaTh48GCwh+Q9Hszljf9/vcLqmriwuuZbrK4NnDfHe6PRiPnz5yMxMREajQY6nQ5Wq9WpzaFDh5Ceno6EhARMmDABmzdv9nqdFosFK1euRHJyMuLi4jB37ly0tbX55gVfxqO9YmZmJh577DGXFyoG27dvx6pVq1BVVQWDwYDCwkLMmTMHx48fD/bQiLowrHlF7mGN1bXQxrA2cN4c761WK2bMmAGtVouWlhY0NTWhsbEROp3O0aa5uRmZmZlYunQpWltbUV1djeLiYrz++uterTM/Px/19fVoaGhAa2srkpKSMHv2bNjtdp++fo8rbLW1tUhNTcXRo0d9OoDBKi0txbJlyzBp0iQAQE5ODtLT07Fp06Ygj2wAeGCXF15cQDLA6prvMKwNjjfH+6qqKphMJpSVlUGpVCIqKgoVFRXYsmULzp49CwBYv349pk+fjnnz5gEAUlJSUFhYiLKyMo/X2draiq1bt6KiogJRUVFQKpUoLy+H0WhEbW2tT1+/R0eT22+/HR988AHWrFmDW2+9FcXFxVi+fLlPBzIQp06dQktLC7Kzs52WZ2dn4+mnn8aGDRtcnmOxWGCxWBz/NpvNAACVIgzhCjF8sxfDGAZP/e3rUMvk9XjtvBLwUyFFrVBc9l95bV/leUVQX5LTtvXTOJRh0qqwqXz0nUMVrnD6b6iJee9LQO2nN1WYdPcDNpvNcRzuplaroVarnZZ5e7zX6/XIzMyESqVyLEtNTUV0dDT0ej1+8pOfQK/XY8WKFS79LVu2DCaTCRaLpd917t+/H7GxsUhNTXU8rlKpkJmZiV27drk8dzA8/iiqVCqsWbMGCxYswI9+9CM8/fTTGDJkiFObTz/91GcD84TRaAQAxMXFOS3XaDSOx3oqKytDaWmpy/I1Y8dj2LBhvh9kiFsTOSHYQwiOSP+vYs3Yq/y/khC1Oml8sIcgW4/dmRjsIQRHpv/eUxcuXEDNf/itexfh3yh8cupC2KWuYFVUVOS0fPXq1SgpKXFa5u3x3mg0YvLkyS7LL29vNBrd9tf9WHdxp691uuuju82xY8dclg+GV9+dTpw4AZ1Oh/DwcBQVFWH48OE+HYy3IiIiAABhPb5dKBQKCILg9jmrVq1ymsM2m82Ij49H8enjCB86xO1zAm64LdgjGDQ1wrAmcgKKzR/Dgs5gDydwzvt/ClStUGDN2KtQfPpTWHp5n0uR8nzwKy9qhQKrk8aj9JPjftm2Ujt/bajJd32pwhV47M5EbNh3Ala7fN63noh570v/riDskn/796OMjAy88cYbTst6VtcA74/3ERERLm17tnfXRvFtlV0QBI/W6cl6fMWjo4sgCKioqEBJSQnmzp2L//3f/0VUVJRPBzIQWq0WANDW1obk5GTH8ra2NkdK7sldqRUArEInwgSxBAuxjGPwLOgMncDWoURg/t917RwsggCLaN6zg6PsCIMdIjiIf7vftQgCLJ0+3tl2KMTxGj007DTgj8vMrHYBVpt0tsNgjan/An6PU0rp7geUSiUiI/ufkvD2eK/Vat1eqXl5e3dtuv99eZ99rdOT9fiKRzXNG2+8Eb/97W/xxz/+Edu3bxdFWAOA2NhYTJ061eXEvt27d2PWrFlBGpUP8ER16eH/swGT+1WhFLp4kYHveHu8nzlzJvbu3Qub7bsZq6amJpw5cwYZGRmONu76mzZtGmJjYz1aZ0ZGBkwmE44cOeJ43GazQa/X+zyHeLSnHDduHI4ePYqcnByfrtwXVqxYgXXr1jnminfs2IE9e/agoKAgyCOjkMArQQclVMKa1KZCeWXo4DGs+Z43x/usrCzExMSguLgYdrsd7e3tWLJkCXJzcxETEwMAKCgowL59+7Bz504AXbf5WLt2rdOFCP2tMyYmBrm5udDpdDCbzbDb7SgqKkJ0dDTmzJnj09fv0ZHmr3/9q09X6ksLFy6E2WxGVlYWOjo6oNFoUFNTg6SkpGAPjeSOQY2I3GBY84++jvcGgwE333wzKisrkZOTA6VSibq6OuTn5yM+Ph5hYWHIyclBeXm5o7/k5GTU1NRAp9Nh8eLFGDZsGEpKSrBgwQKP1tlt48aNWLlyJVJSUmC325GWloa6ujoolb49RsjiiJOXl4e8vLxgD8O3OpTACOlffCBbDGuDxuqaOLG6RmLW2/Feq9XCYDC4LKuuru6zv/T0dLz//vsDWmc3tVqNyspKVFZW9tnPYIXGHpPIlxjWBi1UwhqFHlbXyF+41xQzBgPx4f8TkjFW1waHYY38iUcfIk8wqPlMKFXXpDYdSgPHsEb+Fjp7TqLB4PmEPmMbId17RnnLNkI69xtjdW3gGNYoEBjYxI6VHfEYYXP9owGxjeh0+pMzqYS2C2ODPQLpMt0yOthDoBDANEA0GO5CG0O213qGNrlNm3aHNrFPkV4Yy0obkVjJa68oVwwA0sIq3KDJtfpmGyGIvuLGStvAsMpG/sYkQORvrMINyuWhTS6VN9sIQdTVNlbaiMRHHnu/UMADvLzwfLgBkVPVTezVNlbaiMSFgY1ILNyFOAY5t+Q0ZSr24Eae47Qo+RPLNkRi5xLa+D3rcnKZMhXjhQmcGiUSD+nu3UIRp0XpcsNZjetJDpU3sVXbODXqHVbZyF+YAIjkhBc4OEi58ia2ahsrbUTBJ629GIXswZcGgVU4yVbexHR+GyttnmOVjfyBR3+iUBPiVbju0CalqpvYbwNCRP4nnT0WfSeEDq4UICFYhZNaxU0MlTZW2YiCh4FNqhjayJ9CKMBJKbgxtEkHp0XJ13jUl7IOpewPpiQSITCNahvRiXCF+L/DiuGCBF6EQBR44t87Ud9kdtAkCZFxBc42PPiVrP4Eu9rGSlv/WGUjX2JgkwOGNgo2GYY3KUyVBju0EVHgMLDJBUMbiYXMwpvYg1swQxurbESBw8AmJx1KBjcSFxldvCDm4MbQJl6cFiVfYWCTI4Y2EisZhDex3oSXoY1I3hjY5IqhjcROBtU3sQW3YP4yAkNb71hlI19gYJMzhjaSkp4BTkIh7vKqmxgCXHdwC3R4uzCWwa03DG00WDyiy5mEDnihaFzsVwCAzz6PCuo4RE2i938T0w/PB+O+bd2hjfdqc9Yd2sbUfxHkkZAUiX/PRxQk3YFKbOtRCeHARf+MRRIkFuLEEt4ur7YFKrzxBrvuMbjRQIh3L0cUAIEKZf4wdkw7rAq7R21lX8WTSIgTW3gLRHDrrrapzvp9VZJjumU0Qxt5THx7NCIfk3Io85X+toEsA53IQ5wYwlsgg9s3Y/y+CklitY08JZ69F9EgMJQNTl/bT1ZhrmeIE0mAC3Z4C3RwC2/z+2okh8GN+iOOvRVRLxjEgk/WYa5ngDuvCs44LtPzKtNABrhABTee29Y7BjfqDQObP42wieYbvBgxjElfb/8PJRvkhtsu+2+nKD6/3QFObsGNV5L2jcGNegr+3kiuur+5d/9XBDv+YGI4Cy2yCXKXV+CC/BkOxrSpbYTAaluQMbhRN9441x/cnewc6HuiiegebAxr1G1c7Fdu/yRBRDf1DeQNegNxA17ecLd/pltG8+a7/bBYLFi5ciWSk5MRFxeHuXPnoq2t7xMmDx06hPT0dCQkJGDChAnYvHmzS5tt27Zh8uTJ0Gq1SEtLw8GDB71e7913343Ro0dDq9U6/tLT0716fQxsvtbXjlxEISpQxo5pD/YQSAIkGeREEODkGNyobwxuvcvPz0d9fT0aGhrQ2tqKpKQkzJ49G3a7+9sfNTc3IzMzE0uXLkVrayuqq6tRXFyM119/3dFm+/btWLVqFaqqqmAwGFBYWIg5c+bg+PHjXq3XYDBg+/btMBgMjr8DBw549foY2HzJk512CIS2cbFfMaiRT0gqxAUxvMkpuLHa5hkGN2etra3YunUrKioqEBUVBaVSifLychiNRtTW1rp9zvr16zF9+nTMmzcPAJCSkoLCwkKUlZU52pSWlmLZsmWYNGkSACAnJwfp6enYtGmTV+s1Go2Ij48f1GtkYAsGEUyp+IuoD6gkC5IIcEEKb3ILbtQ/0y2jcTb1e8EexoDZbDaYzWanP4vF4nU/+/fvR2xsLFJTUx3LVCoVMjMzsWvXLrfP0ev1yM7OdlqWnZ2Nw4cPw2Qy4dSpU2hpaXHbprtPT9ZrtVpx5swZJCQkeP26LsfA5isD2THLKLSJ+uBJsib6KlwQg1sgwhurbaEn/IICyo7B/ylsCuj1eowaNcrp7/IKl6eMRiPi4uJclms0GhiNRo+fo9FoHI91P89dm+7HPFlvW1sbhgwZgj/84Q+4/vrrcdVVV2HRokVobW316jWG9qWLvjKYnbEMbv0hugMkhTx378mgX6EapCtOA3FbEH/fBoRXkspXRkYG3njjDadlarXapV1jYyPy8vLc9jFy5EhkZWUhLMz1Pa5QKCAI7r9UREREuDxHoeh6DwuCgIiICABw26a7T3d99GzT3t6OmJgYjBs3Du+++y46OztRVFSEjIwMfPDBBxg+fLjb8fUk7aQgBr745izR0MagRlLS8/0a1AAXhPAm9eDG+7bJk1KpRGRkZL/tUlNT0dDQ0Ovjr776qtsrQtva2hxVs560Wq3Lc7r/fflz2trakJyc7LZPd330bDN16lScPHnS6fGKigq88MILOHDgAGbNmtXr67ocp0TFwpfTJn6efhHltBORl0QzhRrg890CMVXqz/PbOEVK7mRkZMBkMuHIkSOOZTabDXq9vtdANHPmTJcLEnbv3o1p06YhNjYWsbGxmDp1qts23X16ut7OTufPnCAI6OzsdFT0PMHANhj+2MGK/Ly2oB/ciPygO7gF/epmGQY3f2Boo55iYmKQm5sLnU4Hs9kMu92OoqIiREdHY86cOW6fU1BQgH379mHnzp0Aum7zsXbtWqxYscLRZsWKFVi3bh2OHTsGANixYwf27NmDgoICj9f77rvvYuLEiXj//fcBABcvXsSjjz4KrVaLH/7whx6/Rga2gfLnTlWEoU0UlQiiABk7pj247/kAVt38Hdz8VW3jBQnU08aNGzFlyhSkpKRAq9WiubkZdXV1UCq7TjuoqqqCVquFwWAAACQnJ6OmpgZr1qyBRqNBVlYWSkpKsGDBAkefCxcuRHFxMbKyshAXF4e1a9eipqYGSUlJHq/31ltvxeOPP468vDzHTXPb2tqwZ88et+fr9UYh9HY2Xogwm80YNWoU4jesQdjQIZ49KVCBaqDntvh4fAM5aKmEcOgupqFiyHuwKtzftJAGhtvWfzzZtkE79y1A57r56xw3dZgCv0m+Ck/8+zgsnb497IT6eW1hnRY0blmF9vZ2j84HG4xx48ZBNWcBhmoGd4sKAGh74yUUP7AQhYWFPhiZ/LHCJmZBvlUIq2pEroJadQuAQN3HjYi8w8Amdt5Mi/hoh86gRtS3oH1GRHi6hLf8NT1KJHcMbFLR347ah2GNiDwj19AmxSobQxvJHQOblPS2o2ZYIwqaoFTbJB7a/P0j8kRyxMAmNZfvqH14FRnDGtHgMLR52zenRom8wcAmRT4OagxrRL7B0EZE/sLAFsIY1Ih8L+BfgiR8IQKrbESeY2ALUQxrRP4lp9DGqVGi4GNg85aEv812Y1gjCgyGNiLyFdEGNqvViuXLl2P8+PHQaDS45ZZbcODAAac2RqMR8+fPR2JiIjQaDXQ6HaxWa5BGLA0Ma0SBxdDmSb+sshH1R7SBbfHixfjggw/Q0NAAo9GI5cuXY/bs2fjkk08AdAW6GTNmQKvVoqWlBU1NTWhsbIROpwvyyMWLYY0oOBjaPOmXt/og6osoA5vVasVHH32ELVu2YPTo0QCAe++9FxMnTsTf/vY3AF0/4moymVBWVgalUomoqChUVFRgy5YtOHv2bDCHL0oMa0TBFdCLEWRw6oYvsMpGchKYXxP2kkqlwrvvvuu07Ouvv8aJEyccP2yr1+uRmZkJlUrlaJOamoro6Gjo9Xr85Cc/cdu3xWKBxWJx/NtsNnetUxGGcIUn+VWUGbdPY8e0A0J4QNepEsKc/ku+w23rP4HYtleO+RqnTaP81r/DiE7gvJ928SMB5XmFV09RKxTf/be3zRsJKDu867c/9jhgqMmnXYpOmN2324zESZSBrSeTyYT77rsPY8eOxfz58wF0nb82efJkl7YajQZGo7HXvsrKylBaWuqyfM3Y8Rg2bJjvBi0mF4O36gLLDcFbucxx2/qP37dtpH+7D/h6vLA6aXywhyA7Fy5cQP0fgj0K8regBLbGxkbk5eW5fWzkyJF4++23Hf9+++23sWjRItxwww148803MXToUABAREQEwsJcv6YpFAoIQu/nQqxatcrpPDez2Yz4+HgUnz6O8KFD+h/8cGlMNYwd0x7U9auEMBRYbsAmdQOsCl5d5kvctv4T6G0bkEob4L9KGzyvtKkVCqxOGo/ST47D0sc+GvB9lQ2Qd5UtzG7pvxFJXlACW2pqKhoaGvpt98ILL6CwsBAVFRV48MEHnR7TarVoa2tzeU5bWxs0Gk2vfarVaqjVapflVqETYYInO2jxHyDHxX4FsVwra1V0wqqwB3sYssRt6z+B2rbRsV/is8+j/L4eDLcCHf7Z3VuGAcoOD6aQv21iEQRYOvsObJZhgs9DW7g0vmsPSFg/25PkQbQnwezcuRNPPPEE3nnnHZewBgAzZ87E3r17YbN99ylsamrCmTNnkJGREcCRigsvLiCSFv48nHu+vmqUFyCQ1InyHLaOjg78/Oc/x2uvvYaUlBS3bbKyshATE4Pi4mI8+eST6OjowJIlS5Cbm4uYmJgAj1gc5LrTL0zaHbB1PfXJzICti+hy42K/8m+1bYTNb1U224hOz6psQXZhLDDsdLBHQTQwogxsjY2NOHPmDBYtWuTy2C233IKqqioolUrU1dUhPz8f8fHxCAsLQ05ODsrLy/03MBFfKi/VsBbIMOaJYI+HgTG0MbT17NP3U6NEUiXKwDZ9+nR0dvZ/rphWq0V1dXUARjQw3SHK3+eoSCGsBTsISYVH28muBj5Kw6Pj9wHhvZ9szPAnTQxtPfv0bWhjlY2kSpSBTW7cBSpf7ZAZ1qg3Pbc7A5x0SDm0EZF/8BMbJL0FLW920lIIa4+O3xfsIdC3ugMcgxv5E6tsRP7BwOaNDqXfz2PztBon9rD26Ph9wEdpwR4GuXF55Y3hTbxYZSOiy/HTKgFiD2eXc4QBu+u97kh8GN7Eze+hzU+kctUokZQwsPnRZ59HSSpsDRbPVZM2hrcQJKEqG6dFKdRJ45NKosagJj883008pFplIyLfYmCjQWFYkzdW3UKAn6psUrj4gEhKGNhowBjWQgvDW/CwyuYfnBYlKeFZoeS1wqTdDGshjv//ZcZPV7/bRvR/A3Qi8gwDm5/J7VsxD9TUjcE9sELpAqa+8EfhKVQxsJFHeHCm3vB9IRMi/q1kImJgIw/wgEz94XskMKRYZeO0KJFvMLB5SyL3LPIFVtXIG3yvBIZfQ5tEqmycFiV3LBYLVq5cieTkZMTFxWHu3Lloa2vr8zmHDh1Ceno6EhISMGHCBGzevNmlzbZt2zB58mRotVqkpaXh4MGDTo+fP38eO3fuxPz586FSqXDixAmXPoxGI+bPn4/ExERoNBrodDpYrVavXh8DWwBI8Tw2HnxpIBjyyR1W2SgQ8vPzUV9fj4aGBrS2tiIpKQmzZ8+G3W532765uRmZmZlYunQpWltbUV1djeLiYrz++uuONtu3b8eqVatQVVUFg8GAwsJCzJkzB8ePH3e0efTRR/H73/8e48ePx6VLl1zWY7VaMWPGDGi1WrS0tKCpqQmNjY3Q6XRevT4GNnLCAy75At9D/iXFqVFf83WVjaSttbUVW7duRUVFBaKioqBUKlFeXg6j0Yja2lq3z1m/fj2mT5+OefPmAQBSUlJQWFiIsrIyR5vS0lIsW7YMkyZNAgDk5OQgPT0dmzZtcrTZsmUL6urq8Mgjj7hdT1VVFUwmE8rKyqBUKhEVFYWKigps2bIFZ8+e9fg1MrARkV8wtEmURKZFfY3TotK2f/9+xMbGIjU11bFMpVIhMzMTu3btcvscvV6P7Oxsp2XZ2dk4fPgwTCYTTp06hZaWFrdteuuzt/VkZmZCpVI5lqWmpiI6Ohp6vd7jfhjYyIEHWPI1vqf8R2pVNttwVsTIlc1mg9lsdvqzWCxe92M0GhEXF+eyXKPRwGg0evwcjUbjeKz7ee7a9Nanr8bmTuicQU994oGV/IW/SypBEvlReP4gfHAM+QIYFj74fsIvdlWfioqKnJavXr0aJSUlTssaGxuRl5fntp+RI0ciKysLYWGuNSiFQgFBcP9lISIiwuU5CkXX+0kQBERERACA2za99enpegbSj/g/kTLx2edRov1GzLBGgVCYtJuhzcf4k1UkdRkZGXjjjTeclqnVapd2qampaGho6LWfV1991e0VoW1tbY6qWU9ardblOd3/vvw5bW1tSE5O9qhPT9czkH44JUpEAcMvBxISoueyUWAplUpERkY6/bkLbP3JyMiAyWTCkSNHHMtsNhv0ej1mzZrl9jkzZ850uSBh9+7dmDZtGmJjYxEbG4upU6e6bdNbn72tZ+/evbDZvvtMNTU14cyZM8jIyPC4Hwa2EMcDKAUa33O+JdbKfaDwnmwEADExMcjNzYVOp4PZbIbdbkdRURGio6MxZ84ct88pKCjAvn37sHPnTgBdt/lYu3YtVqxY4WizYsUKrFu3DseOHQMA7NixA3v27EFBQYHHY8vKykJMTAyKi4tht9vR3t6OJUuWIDc3FzExMR73w8AWwnjgpGDhey908eID8peNGzdiypQpSElJgVarRXNzM+rq6qBUdp39VVVVBa1WC4PBAABITk5GTU0N1qxZA41Gg6ysLJSUlGDBggWOPhcuXIji4mJkZWUhLi4Oa9euRU1NDZKSkjwel1KpRF1dHT788EPEx8fj2muvxdSpU/HMM8949fp4DttAdCg5XUBE8ieRiw+IgK5z3yorK1FZWen28ZycHOTk5DgtS09Px/vvv99nv3l5eb1e8HC5xMTEXi8i0Gq1qK6u7rePvrDCFqJY4aBg43uQfIU30aVQwMBGRCRxoX4em6/xPDYSIwa2EMTKBokF34uhib8tSuQ9BjYiCiqGNpHj+bpEosDAFmJ4cCSSp1CfFuXtPUjuGNgGildOEfkMv0iEHk6LEnmHgS2E8KBIYsb3p4hxWpQo6BjYiIhIFjgtSnLGwBYiWL0gKeD7dHBC/Tw2IjljYCMioqDgeWxEnmNgCwGsWpCU8P0qUjyPjSioGNgGQwJXivLgR0ShhOexkVwxsBGR6PCLxsDxPDYieWJgkzEe9EjK+P4NDTyPjcgzDGxEROQZiZzHxmlRkiMGNplidYLkgO/jgeG0KJH8MLAREVFQcVqUqH8MbIMlwitFWZUgIiKSFwY2IhI1fgERGZ7HRhQUDGyDJbKdFw9uRERE8sPARkQkQ7zwgEheGNhkhNU1IpIqXnhA1DcGNplgWCM54/ubBoLnsZGcMLAREZF3RHbuLlEoYGCTAVYfiIiI5I2BjYhIpnjhAZF8MLBJHKtrRCQXvPCAqHcMbAHij2+6DGsUSvh+p4HghQckFwxsRETkPV54QBRQDGwSxWoDERFR6JBEYPvb3/4GhUKBbdu2OS03Go2YP38+EhMTodFooNPpYLVaAzcwfsMkIpHjhQdE8iD6wGYymbBkyRIkJSU5LbdarZgxYwa0Wi1aWlrQ1NSExsZG6HS6II2UiIgGSwoXHvA8NgoG0Qe2hx56CP/1X/8FrVbrtLyqqgomkwllZWVQKpWIiopCRUUFtmzZgrNnzwZptIHB6VAKVXzv00D4+sIDomAQdWB77rnnYDAY8Nhjj7k8ptfrkZmZCZVK5ViWmpqK6Oho6PX6QA6TiCg08bQQEhGLxYKVK1ciOTkZcXFxmDt3Ltra2vp8zqFDh5Ceno6EhARMmDABmzdvdmmzbds2TJ48GVqtFmlpaTh48KDT41arFcuXL8f48eOh0Whwyy234MCBA05t7r77bowePRpardbxl56e7tXrU3rVOoCam5vx+OOP4x//+AciIiJcHjcajZg8ebLLco1GA6PR2Gu/FosFFovF8W+z2QwAUCnCEK7wNr963l4lhHvZdx/sat/15S92lfN/yXdCfNv69LPk0neY03/lQu3P7+Ye7jfVCsVl/+37OcowxWBH5fc+VSI6eobZfb+9pCg/Px8ff/wxGhoaMGLECCxfvhyzZ8/G4cOHER7uut9obm5GZmYmtm7dinnz5uHDDz9ERkYGvve97+G+++4DAGzfvh2rVq2CXq/HpEmTUFVVhTlz5uBf//oXxo8fDwBYvHgxDAYDGhoaMHr0aLz55puYPXs2PvjgA8fpXAaDAdu3b8fs2bMH/PqC8pZrbGxEXl6e28dGjhyJPXv2YNGiRfj1r3+Na6+91m27iIgIhIW5fugVCgUEoffyd1lZGUpLS12Wrxk7HsOGDfPwFQzARR/29VGaDzvzr/BjrtVR8o1Q3baBOEu1wHJDANYSQJHi6XvN2Kv8M44QduHCBdT/IdijCK7W1lZs3boV7733HqKiogAA5eXliIuLQ21tLbKzs12es379ekyfPh3z5s0DAKSkpKCwsBBlZWWOwFZaWoply5Zh0qRJAICcnBz86U9/wqZNm7BhwwZYrVZ89NFHePXVVzF69GgAwL333ouJEyfib3/7G375y18C6CoyxcfHD+o1BiWwpaamoqGhodfHi4qKEBkZiV/96le9ttFqtW5LnW1tbdBoNL0+b9WqVU4XJpjNZsTHx6P49HGEDx3i4Sv41nDPpwPGjmn3ru9ePDp+n0/68Tu7CuHHHoP96g1AeACv3A0FIb5tnzl+p9/6VglhKLDcgE3qBlgV4j/53RunTaP80/F5zw4jaoUCa8ZeheLTn8LSx5fqbsrzvq0aKTt8X4UaavJ5lwMSZrf030ikbDabY6arm1qthlrt3UzS/v37ERsbi9TUVMcylUqFzMxM7Nq1y21g0+v1WLFihdOy7OxsLFu2DCaTCRaLBS0tLS7Pzc7OxtNPP40NGzZApVLh3XffdXr866+/xokTJxAZ2fVtxmq14syZM0hISPDqNfUkoqLud2pra/HBBx+4VND279+P3NxcXLp0CTNnzkReXh5sNhuUyq6X0dTUhDNnziAjI6PXvnt7I1iFToQJ3u6gPW9vVdi97LsX4RL7YIZbpTdmqQjRbftoci2e+mSmX9dhVXT67jMrEhYv9lde8Xi/2bU/twgCLB48x97p22lcyzDB56EtXCSn8IV1Bvaiiuij7Rg17ItB9/PZWQv0ej2Kioqclq9evRolJSVe9WU0GhEXF+eyXKPR4NixYx4/p7vgYzQaHadPuWvT26lXJpMJ9913H8aOHYv58+cD6CokDRkyBH/4wx/wyiuvoL29HbfccgvKysq8CnGiPFHj3//+NwRBcPqbPn06tm7dCkEQoFQqkZWVhZiYGBQXF8Nut6O9vR1LlixBbm4uYmJigv0SnPA+SERERK4yMjLQ3t7u9Ldq1SqXdo2Njbjhhhvc/t1xxx0DOk3K3XMU355rKQiC4/x5d23c9fn2229j2rRpiIqKwj/+8Q8MHToUANDe3o6YmBiMGzcO7777Lo4ePYorrrgCGRkZOH/+vAdbqYsoK2yeUCqVqKurQ35+PuLj4xEWFoacnByUl5cHZgC8OoqIqGtf2CHZQwkFmVKpdEwd9qW/U6leffVVr0+TcndqVfe/L39OW1sbkpOT++zzhRdeQGFhISoqKvDggw86PTZ16lScPHnSaVlFRQVeeOEFHDhwALNmzer1dV1OMp+yv//97y7LtFotqqurAz8YIiIiEo2MjAyYTCYcOXIE1113HYCu8+P0ej2ee+45t8+ZOXMmamtr8fDDDzuW7d69G9OmTUNsbCyArrBVW1vruHigu83lIWvnzp144okn8M477yAlJcXtujo7O50qdYIgoLOz01HR84Qop0SJiIiIPBUTE4Pc3FzodDqYzWbY7XYUFRUhOjoac+bMcfucgoIC7Nu3Dzt37gTQdZuPtWvXOl2IsGLFCqxbt85xHtyOHTuwZ88eFBQUAAA6Ojrw85//HK+88kqvYe3dd9/FxIkT8f777wMALl68iEcffRRarRY//OEPPX6NkqmwERERicmFscCw08EeBXXbuHEjVq5ciZSUFNjtdqSlpaGurs5xYWJVVRV+9atf4dChQ9BqtUhOTkZNTQ10Oh0WL16MYcOGoaSkBAsWLHD0uXDhQpjNZmRlZaGjowMajQY1NTWO+6s1NjbizJkzWLRokct4brnlFlRVVeHWW2/F448/jry8PJhMJly8eBHp6enYs2ePV1fDMrAREZHo2EZ0QtnBSSDynFqtRmVlJSorK90+npOTg5ycHKdl6enpjspXb/Ly8nq9d+z06dPR2dn/Vc8PPPAAHnjggX7b9YWfhoHw4oIDXiFKRCQO/E1RkjIGNiKiEMAvj0TSxsBGREREJHIMbEREREQix8DmLd4wl4iIiAKMgc2PeM4IEYWEEP4ie2FssEdAoYKBjYiIiEjkGNiIiEiUbCP6v7+V933y1h4kTQxsRERERCLHwOYN3jCXiIiIgoCBjYiIiEjkGNiIiIgGgVeKUiAwsPkBp0OJSIy4byKSLgY2IiIiIpFjYPNUCN8YkoiIiIKLgY2IiAZPQl9qeS82kiIGNh/jOSJERETkawxsREQkWv74tQMiKWJgIyIiIhI5BjYiIiIikWNgk5inPpkZ7CEQERFRgDGw+RAvOCAiCk38tQPyNwY2IiIiIpFjYPOEhO4vRERERPLDwEZEREQkcgxsPsLz14iIiMhfGNiIiIiIRI6BTYJ4aw8iIqLQwsBGREREJHIMbD7A89eIAocVZvIF2wgh2EMg8goDm0TxoEVERBQ6GNiIiIiIRI6BrT+8aS4REXmAP09F/sTAJmGcFiUiIupisViwcuVKJCcnIy4uDnPnzkVbW1ufzzl06BDS09ORkJCACRMmYPPmzS5ttm3bhsmTJ0Or1SItLQ0HDx50eryiogIjRoyAVqt1+jt9+rSjjdFoxPz585GYmAiNRgOdTger1erV62NgGyRecEBERBR8+fn5qK+vR0NDA1pbW5GUlITZs2fDbre7bd/c3IzMzEwsXboUra2tqK6uRnFxMV5//XVHm+3bt2PVqlWoqqqCwWBAYWEh5syZg+PHjzvaGAwG5Ofnw2AwOP2NHdtVcrVarZgxYwa0Wi1aWlrQ1NSExsZG6HQ6r14fAxsRERFJWmtrK7Zu3YqKigpERUVBqVSivLwcRqMRtbW1bp+zfv16TJ8+HfPmzQMApKSkoLCwEGVlZY42paWlWLZsGSZNmgQAyMnJQXp6OjZt2uRoYzQaER8f3+vYqqqqYDKZUFZWBqVSiaioKFRUVGDLli04e/asx6+RgY2IiIiCwmazwWw2O/1ZLBav+9m/fz9iY2ORmprqWKZSqZCZmYldu3a5fY5er0d2drbTsuzsbBw+fBgmkwmnTp1CS0uL2zaX92kwGJCQkNDr2PR6PTIzM6FSqRzLUlNTER0dDb1e7/FrVHrckkTpqU9mojBpd7CHQUREIeDCaQu+QCNGKL4eVD+dgh3tOImGhgaMGjXK6bHVq1ejpKTEq/6MRiPi4uJclms0Ghw7dszj52g0Gsdj3cHRXRuj0ejUz+HDh/HUU0/BYDAgKSkJpaWl+MEPfuB4fPLkyW7Hdnk//WFgIyIiIo+8uusV3D17LuKEREQoVP0/oRcGfIpwKPHSSy9h69atTo+p1WqX9o2NjcjLy3Pb18iRI5GVlYWwMNdJQ4VCAUFwf5PkiIgIl+coFAoAgCAIiIiIAAC3bS7vU6VS4ZtvvkF1dTVGjRqFv/zlL5gxYwYOHTqE6667zu16+hubOwxsRCQZvDKaKLhmzpyJkYjCCfwfJuC6AfVxSbDiOD7Czl3VGDZsmEfPSU1NRUNDQ6+Pv/rqq26vCG1ra3NUzXrSarUuz+n+9+XPaWtrQ3Jycq999qzgLVq0CNu3b8crr7yC6667zu16+hubOzyHjYiIiDyiUCjw1uE6nMIn+EY4P6A+TuD/MBJRmDnTd1/AMjIyYDKZcOTIEccym80GvV6PWbNmuX3OzJkzXS5I2L17N6ZNm4bY2FjExsZi6tSpbttc3mdnZ6dL33a73VGtmzlzJvbu3Qub7bv7ujY1NeHMmTPIyMjw+DUysBEREZHHrr/+eoyBBp+gyevnfiOcxyl8grcO1zkCjS/ExMQgNzcXOp0OZrMZdrsdRUVFiI6Oxpw5c9w+p6CgAPv27cPOnTsBdN3mY+3atVixYoWjzYoVK7Bu3TpHFW3Hjh3Ys2cPCgoKAABfffUVJkyYgJdffhmdnZ0QBAEvvvgiDhw4gPvvvx8AkJWVhZiYGBQXF8Nut6O9vR1LlixBbm4uYmJiPH6NDGxERETklX+ePAATjDAL57x63idowhhocP311/t8TBs3bsSUKVOQkpICrVaL5uZm1NXVQansOvurqqoKWq0WBoMBAJCcnIyamhqsWbMGGo0GWVlZKCkpwYIFCxx9Lly4EMXFxcjKykJcXBzWrl2LmpoaJCUlAQCioqLw8ssv48UXX0R8fDxiYmLwP//zP6itrXXcCkSpVKKurg4ffvgh4uPjce2112Lq1Kl45plnvHp9CsGbM95kyGw2Y9SoUYjfsAZhQ4e4Nujjp6nEctNcUV4lalcj/KNVsE8qA8K9v0Sb+hDC29bf57CphHDoLqahYsh7sCrc32xT6j77PMp/nXf0flq0WhGGdeOSsPyzT2ARXKeQ+qLs8E9tQdnhuwpPt2Gn+2/ja2GdFjRuWYX29nZERkYGbL2Jimtgxpf4Pm73qFpmFs6hAX/H8ZOf9nkbDHKPFTYZ4InYROQJv4Y1CjkffPVPdKAdX6D/lCoIAj7GEcQjmWFtgBjYiEgS+MWESFxGjRqF8ZiEj3EUnf1UTb/AaXSgHf8+dyhAo5MfBjYiIiIakKOW99AJOz7DyV7bdFXXjmI8JiEqKipwg5MZBjYiIiIaEJVKhRdf24pP0AS74P6c7zacQCfsOGp5L8CjkxcGNiIiIhqw++67D0MxDCfh+hNQdsGGT9CEF1/b6vRbmuQ9Bra+9HGFKBEREXXdTLfmQDVO4hgswkWnx07iGIZgGO67774gjU4+RB3Y9Ho90tLSoNVqkZCQgJUrV8JqtToeNxqNmD9/PhITE6HRaKDT6ZweDyU8IZuIiILltttuQzTG4FN86FhmES7iJI7hbweqfXqT3FAl2sB26NAhLFq0CL/73e9gMBjwz3/+E2+99Rb27t0LALBarZgxYwa0Wi1aWlrQ1NSExsZG6HS6II+ciHyNX0iIxO/A/+nxGU7gvGAGABzHh4jGGNx2221BHpk8iDawrVixAsuWLcNNN90EABg3bhzee+89x09MVFVVwWQyoaysDEqlElFRUaioqMCWLVtw9uzZYA6diIgo5EycOBFxGI8W/C/OC2a04QQO/J8+2MOSDVEGtnPnzuHAgQO46667nJaHhX03XL1ej8zMTKeTGFNTUxEdHQ293v9vELH8ygEREZFYNH7+Lr6ECUdwCHEYj4kTJwZ7SLLR+++IBNEnn3wCQRBw4cIFzJ49Gx9++CHi4uKwbNkyzJs3D0DX+WuTJ092ea5Go4HRaOy1b4vFAovlu5/zaW9vBwAoLVaEKXrk1/DeLzoI6xDfBQkX2sODPYTv2MMQfuEC7O1hQLiIxiUHIbhtA/V5CxM6ccF6AWG2SwhTePfzSVKg/MbP5/he7H2bKRUKXLhwAcqLFti9/EVEpdU/5z/5o9+wILxtwuxdxzQx/NLkmDFjkIiJOIFmNJw+GOzhyEpQAltjYyPy8vLcPjZy5EiUl5cDAEpKSvDss88iMTER+/btw7333guVSoXs7GxEREQ4Vdy6KRSKPt+0ZWVlKC0tdVl+vGiNV6/hE69aB8a7wR6Ai/pgD0DGQm3bNgZsTX8P2JpCz38EewAy98UXX2DUqFHBHgZahP8N9hBkSZQ//v7ZZ58hLi4O//73vzF16lTH8sWLF+Ps2bOoqqrC4sWL8fXXX2P79u1Oz9VqtdiwYQPmz5/vtu+eFbbOzk58+eWXGD16NK9i8SGz2Yz4+HicOnUqoD9GHAq4bf2H29Z/uG39p729HQkJCTh37hx/SUDGRDklOnbsWCQnJzsFq25qtRoAMHPmTOTl5cFms0Gp7HoZTU1NOHPmDDIyMnrtW61WO/roxje4/0RGRnLn7Cfctv7Dbes/3Lb+427WieRDlP93FQoFVq9ejUceeQQnT3b9Ptnf//53vPzyy8jNzQUAZGVlISYmBsXFxbDb7Whvb8eSJUuQm5uLmJiYYA6fiIiIyKdEWWEDgP/8z//E119/jTvvvBNff/01xo4di23btuHOO+8EACiVStTV1SE/Px/x8fEICwtDTk6O4/w3IiIiIrkQbWADus5ZW7x4ca+Pa7VaVFdXB3BE5Cm1Wo3Vq1e7TD/T4HHb+g+3rf9w2/oPt21oEOVFB0RERET0HVGew0ZERERE32FgIyIiIhI5BjYiIiIikWNgowHbtm0bJk+eDK1Wi7S0NBw82PvPkJw6dQrz589HQkIC4uPjce+996K1tTWAo5UWb7bt5ZYvXw6FQoETJ074d4AS5u223bRpEyZOnAiNRoOUlBRs27YtMAOVIG+27Y4dO3DrrbciPj4eiYmJuO+++/Dxxx8HcLTSYLfbcfDgQSxduhTR0dH9vv+MRiPmz5+PxMREaDQa6HQ6WK1+/kkyCgyBaABeeuklYezYscKHH34oCIIgvPbaa8KoUaOETz/91KWt1WoVJk6cKKxcuVKwWq2CzWYTli1bJlx77bXCpUuXAj100fNm215Or9cLU6dOFQAIx48fD8BIpcfbbbthwwbhhhtuEIxGoyAIglBfXy9cddVVgsFgCNiYpcKbbfv2228LQ4cOFXbv3i0IQtc+Ys2aNYJGoxEuXLgQ0HGL3ZYtW4Qbb7xR+PWvfy1cccUVwtatW3tta7FYhEmTJgk6nU64dOmScO7cOeH2228X8vPzAzdg8hsGNhqQ5ORkYf369U7LsrKyBJ1O59L2yJEjwp133um0zGw2CwCEDz74wK/jlCJvtm23L7/8UkhISBAOHjzIwNYHb7at2WwWhg8fLjQ0NDgtt9vtfh2jVHmzbdetWydMnTrVadm5c+cEAEJjY6M/hylpV155ZZ+Bbfv27cLo0aMFi8XiWNbQ0CCo1WrhzJkzARgh+ROnRMlrp06dQktLC7Kzs52WZ2dnY9euXS7tp0yZgrfeestp2dGjRwEAI0eO9N9AJcjbbdtt8eLFyMrKwq233urvIUqWt9tWr9dj+PDhSE1NdVrOn/9x5e22/fGPf4y2tjasXbsWX3/9NVpbW7F48WLExsbi6quvDtSwZUev1yMzMxMqlcqxLDU1FdHR0dDr9UEcGfkC9zzkNaPRCACIi4tzWq7RaByP9aWxsRE5OTl48MEHMX78eL+MUaoGsm1feukl/Otf/8JTTz3l9/FJmbfb9uOPP0ZiYiJ27tyJtLQ0JCYm4q677sKRI0cCMl4p8XbbJiUl4f3338dLL72EmJgYXHnllbDZbDh48CBGjBgRkDHLkdFodPl/AHi+byZxY2Ajr0VERABwrTQoFAoI/dyHeePGjUhPT8eDDz6ILVu2+G2MUuXttj1x4gSWLl2Kl156CcOGDQvIGKXK221rt9vx8ccfo7a2Fm+99RaOHTuGO+64A+np6TAYDAEZs1R4u22NRiNmzZqFe++9F21tbTh69Cg6OztZBRqkiIgItxVgT/bNJH4MbOQ1rVYLAGhra3Na3tbWBo1G4/Y5nZ2deOihh7Bp0ya8/fbbWLt2LcLDw/0+VqnxZtt2dnbipz/9KZYsWYK0tLSAjVGqvH3fJiQkIDw8HM8++ywiIyOhUqlQWFiIuLg4/iReD95u2zVr1mDSpEkoKytDdHQ0Jk+ejO3bt6O4uBj79+8PyJjlSKvVuvw/APreN5N0MLCR12JjYzF16lTU1tY6Ld+9ezdmzZrl9jmFhYU4duwYGhoacNNNNwVimJLkzbY1m8145513UFpaCoVC4fgDgPHjx+O2224L2LilwNv37S233AKgq9LWE3+z0Zm327a1tdVl6m7o0KEYM2YM/vnPf/p1rHI2c+ZM7N27FzabzbGsqakJZ86cQUZGRhBHRj4R1EseSLJeeeUVQaPRCM3NzYIgCMKbb74pREZGCi0tLS5tDx06JFxxxRXC2bNnAz1MSfJm27oDXiXaK2+37UMPPST89Kc/FTo6OgSbzSZUVFQIV1xxhfD5558HctiS4M22ff7554Xo6GjhwIEDgiAIwqVLl4QNGzYISqWSV4n2ob+rRC9duiRce+21wsqVKwWbzSZ89dVXwh133CHk5eUFbpDkNwxsNGD/8z//I0yYMEEYN26ccMMNNwj/+Mc/BEEQhFOnTgkajUZ47bXXBEEQhJKSEmHIkCGCRqNx+duwYUMwX4Joebpt3WFg65s32/abb74RHn30UWHcuHFCbGyscOedd/JWNH3wZttu3rxZuP7664W4uDhhzJgxwo9+9CNh3759wRq6JPQMbO6266lTp4S7775bGDdunKDRaISlS5cKFy9eDMJoydcUgsAzEYmIiIjEjOewEREREYkcAxsRERGRyDGwEREREYkcAxsRERGRyDGwEREREYkcAxsRERGRyDGwEREREYkcAxsRERGRyDGwEZFP1dfXQ6lUoqqqymm51WrFlClTcN999wVpZERE0sXARkQ+dcstt+Dxxx9Hfn4+zp4961heWlqK9vZ2bN68OYijIyKSJv40FRH5nN1uR3p6OuLj4/Hqq6+ioaEBP/jBD7B3717cfvvtwR4eEZHkMLARkV98+umnmDZtGp577jmUlZXhnnvuwZNPPhnsYRERSRIDGxH5zYsvvoif/exnuPHGG/HOO+9AqVQGe0hERJLEc9iIyG/OnTuH8PBwtLe3w2q1Bns4RESSxcBGRH5x5MgRFBUVoaamBuHh4Vi6dGmwh0REJFmcEiUin/vmm29www034M4778TGjRvR2NiIm2++Ga+99hruvffeYA+PiEhyWGEjIp9btmwZLBYLysvLAQCpqakoKirCww8/DKPRGOTRERFJDytsRORTf/3rX3HPPffg7bffdrqFx6VLl3DzzTcjKioKe/fuRVgYvy8SEXmKgY2IiIhI5PgVl4iIiEjkGNiIiIiIRI6BjYiIiEjkGNiIiIiIRI6BjYiIiEjkGNiIiIiIRI6BjYiIiEjkGNiIiIiIRI6BjYiIiEjkGNiIiIiIRI6BjYiIiEjk/j9UBQD3ddd1PgAAAABJRU5ErkJggg==", + "text/plain": [ + "Figure(PyObject
)" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "#-- Drawing\n", + "##########\n", + "# Plot #\n", + "##########\n", + "using PyPlot\n", + "\n", + "rc(\"font\", family=\"IPAPGothic\")\n", + "fig = figure(\"pyplot_majorminor\",figsize=(7,5))\n", + "\n", + "nx = na\n", + "ny = ne\n", + "\n", + "xax = reshape(zeros(nx,ny),nx,ny)\n", + "yax = reshape(zeros(nx,ny),nx,ny)\n", + "\n", + "for i in 1:nx\n", + " xax[i,1:ny] .= α[i] # xmin + (xmax-xmin) * (i-1) / (nx-1)\n", + "end\n", + "for i in 1:ny\n", + " yax[1:nx,i] .= ε[i] # ymin + (ymax-ymin) * (i-1) / (ny-1)\n", + "# yax[1:nx,i] .= Re[i] # ymin + (ymax-ymin) * (i-1) / (ny-1)\n", + "end\n", + "\n", + "cont_val = reshape(zeros(nx,ny),nx,ny)\n", + "shade_val = reshape(zeros(nx,ny),nx,ny)\n", + "\n", + "for i in 1:nr#div(nr,2)-1:div(nr,2)+1\n", + " cont_val = real.(c_full[1:na,1:ne,i])\n", + " for j in 1:ne\n", + " shade_val[1:na,j] = α[1:na] .* imag.(c_full[1:na,j,i])\n", + " end\n", + "#draw_num2 = 2\n", + "#if draw_num2 == 1\n", + "# cp = contourf(xax[1:nx,1:nx], xax'[1:nx,1:nx], Pftime[1:nx,1:nx,1], levels=[-15.0, -10.0, -5.0, 0.0, 5.0, 10.0, 15.0, 20.0], origin=\"image\", cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + " #cp = contourf(xax[1:nx,1:nx], xax'[1:nx,1:nx], Pftime[1:nx,1:nx,2], levels=[-2.0, -1.0, -0.5, -0.25, 0.25, 0.5, 1.0, 1.5, 2.0], origin=\"image\", cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + " #cp = contourf(xax[1:nx,1:no], xax'[1:nx,1:no], Pftime[1:nx,1:nx,2]*Hop'*inv(Hop*Pftime[1:nx,1:nx,2]*Hop'+Ro), origin=\"image\", cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + "#elseif draw_num2 == 2\n", + " ##cp = contour(xax[1:nx,1:ny], yax[1:nx,1:ny], shade_val[1:nx,1:ny], 8, levels=[-0.01, -0.005, -0.003, -0.002, -0.001, -0.0005, 0.0, 0.0005, 0.001, 0.002, 0.003, 0.005, 0.01])\n", + " cp = contourf(xax[1:nx,1:ny], yax[1:nx,1:ny], shade_val[1:nx,1:ny], 9, levels=[-0.005, -0.0025, -0.001, -0.0005, 0.0, 0.0005, 0.001, 0.0025, 0.005], cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + " ##cp = contourf(xax[1:nx,1:ny], yax[1:nx,1:ny], cont_val[1:nx,1:ny], 8, levels=[-1.0, -0.5, -0.2, -0.1, 0.0, 0.1, 0.2, 0.5, 1.0], cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + " #cp = contourf(xax[1:nx,1:ny], yax[1:nx,1:ny], shade_val[1:nx,1:ny], 10, cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + "#elseif draw_num2 == 3\n", + "# cp = contourf(xoax[1:nx,1:nto-1], toax[1:nx,1:nto-1], evec_max[1:nx,1:nto-1], 10, cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + "#elseif draw_num2 == 4\n", + "# cp = contourf(xoax[1:nx,1:40], toax[1:nx,1:40], x_inc[1:nx,1:40], levels=[-2.0, -1.0, -0.5, -0.25, 0.25, 0.5, 1.0, 1.5, 2.0], cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + "#elseif draw_num2 == 5\n", + "# cp = contourf(xoax[1:no,1:40], toax[1:no,1:40], y_innov[1:no,1:40], levels=[-2.0, -1.0, -0.5, -0.25, 0.25, 0.5, 1.0, 1.5, 2.0], cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + "#end\n", + "#ax.label(cp, inline=1, fontsize=10)\n", + "#legend()\n", + "ax = gca()\n", + "\n", + "xlabel(\"X\")\n", + "ylabel(\"Y\")\n", + "if i == 1 #div(ne,2)-1\n", + " plt.colorbar(cp)\n", + "end\n", + "grid(\"on\")\n", + "\n", + "PyPlot.title(\"Re = Re[i]\")\n", + "\n", + "#########################\n", + "# Set tick dimensions #\n", + "#########################\n", + "#ax.xaxis.set_tick_params(which=\"major\",length=5,width=2,labelsize=10)\n", + "#ax.xaxis.set_tick_params(which=\"minor\",length=5,width=2)\n", + "\n", + "fig.canvas.draw() # Update the figure\n", + "gcf() # Needed for IJulia to plot inline\n", + "fname = \"Re\" * string(i) * \".png\"\n", + "#savefig(\"B-init_2d.pdf\")\n", + "savefig(fname)\n", + "#psi\n", + "end\n", + "println(shade_val[1:na,1])" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "8aa57851", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "20-element Vector{Float64}:\n", + " 0.001\n", + " 0.05357894736842105\n", + " 0.1061578947368421\n", + " 0.15873684210526315\n", + " 0.2113157894736842\n", + " 0.26389473684210524\n", + " 0.3164736842105263\n", + " 0.36905263157894735\n", + " 0.4216315789473684\n", + " 0.47421052631578947\n", + " 0.5267894736842105\n", + " 0.5793684210526316\n", + " 0.6319473684210526\n", + " 0.6845263157894736\n", + " 0.7371052631578947\n", + " 0.7896842105263158\n", + " 0.8422631578947368\n", + " 0.8948421052631579\n", + " 0.9474210526315789\n", + " 1.0" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "imag.(c_full[10,1:ne,1])\n", + "α" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "39f4181e", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.6.1", + "language": "julia", + "name": "julia-1.6" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.6.1" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Lilly1966/OSE_solve.ipynb b/Lilly1966/OSE_solve.ipynb index 653c929..36b858d 100644 --- a/Lilly1966/OSE_solve.ipynb +++ b/Lilly1966/OSE_solve.ipynb @@ -74,7 +74,7 @@ }, { "cell_type": "code", - "execution_count": 33, + "execution_count": 3, "id": "7e7d6392", "metadata": {}, "outputs": [ @@ -614,58 +614,33 @@ "Dcal = F.α ./ F.β\n", "for i in 1:Ntot\n", " println(\"Schur(Eig): \", Dcal[i])\n", - "end\n", - "\n" + "end\n" ] }, { "cell_type": "code", - "execution_count": 35, - "id": "1255c473", + "execution_count": 4, + "id": "d6476cfa", "metadata": {}, "outputs": [ { "data": { - "text/latex": [ - "The analogue of IPython's \\texttt{\\%matplotlib} in Julia is to use the \\href{https://github.com/stevengj/PyPlot.jl}{PyPlot package}, which gives a Julia interface to Matplotlib including inline plots in IJulia notebooks. (The equivalent of \\texttt{numpy} is already loaded by default in Julia.)\n", - "\n", - "Given PyPlot, the analogue of \\texttt{\\%matplotlib inline} is \\texttt{using PyPlot}, since PyPlot defaults to inline plots in IJulia.\n", - "\n", - "To enable separate GUI windows in PyPlot, analogous to \\texttt{\\%matplotlib}, do \\texttt{using PyPlot; pygui(true)}. To specify a particular gui backend, analogous to \\texttt{\\%matplotlib gui}, you can either do \\texttt{using PyPlot; pygui(:gui); using PyPlot; pygui(true)} (where \\texttt{gui} is \\texttt{wx}, \\texttt{qt}, \\texttt{tk}, or \\texttt{gtk}), or you can do \\texttt{ENV[\"MPLBACKEND\"]=backend; using PyPlot; pygui(true)} (where \\texttt{backend} is the name of a Matplotlib backend, like \\texttt{tkagg}).\n", - "\n", - "For more options, see the PyPlot documentation.\n", - "\n" - ], - "text/markdown": [ - "The analogue of IPython's `%matplotlib` in Julia is to use the [PyPlot package](https://github.com/stevengj/PyPlot.jl), which gives a Julia interface to Matplotlib including inline plots in IJulia notebooks. (The equivalent of `numpy` is already loaded by default in Julia.)\n", - "\n", - "Given PyPlot, the analogue of `%matplotlib inline` is `using PyPlot`, since PyPlot defaults to inline plots in IJulia.\n", - "\n", - "To enable separate GUI windows in PyPlot, analogous to `%matplotlib`, do `using PyPlot; pygui(true)`. To specify a particular gui backend, analogous to `%matplotlib gui`, you can either do `using PyPlot; pygui(:gui); using PyPlot; pygui(true)` (where `gui` is `wx`, `qt`, `tk`, or `gtk`), or you can do `ENV[\"MPLBACKEND\"]=backend; using PyPlot; pygui(true)` (where `backend` is the name of a Matplotlib backend, like `tkagg`).\n", - "\n", - "For more options, see the PyPlot documentation.\n" - ], "text/plain": [ - " The analogue of IPython's \u001b[36m%matplotlib\u001b[39m in Julia is to use the PyPlot package\n", - " (https://github.com/stevengj/PyPlot.jl), which gives a Julia interface to\n", - " Matplotlib including inline plots in IJulia notebooks. (The equivalent of\n", - " \u001b[36mnumpy\u001b[39m is already loaded by default in Julia.)\n", - "\n", - " Given PyPlot, the analogue of \u001b[36m%matplotlib inline\u001b[39m is \u001b[36musing PyPlot\u001b[39m, since\n", - " PyPlot defaults to inline plots in IJulia.\n", - "\n", - " To enable separate GUI windows in PyPlot, analogous to \u001b[36m%matplotlib\u001b[39m, do \u001b[36musing\n", - " PyPlot; pygui(true)\u001b[39m. To specify a particular gui backend, analogous to\n", - " \u001b[36m%matplotlib gui\u001b[39m, you can either do \u001b[36musing PyPlot; pygui(:gui); using PyPlot;\n", - " pygui(true)\u001b[39m (where \u001b[36mgui\u001b[39m is \u001b[36mwx\u001b[39m, \u001b[36mqt\u001b[39m, \u001b[36mtk\u001b[39m, or \u001b[36mgtk\u001b[39m), or you can do\n", - " \u001b[36mENV[\"MPLBACKEND\"]=backend; using PyPlot; pygui(true)\u001b[39m (where \u001b[36mbackend\u001b[39m is the\n", - " name of a Matplotlib backend, like \u001b[36mtkagg\u001b[39m).\n", - "\n", - " For more options, see the PyPlot documentation." + "Figure(PyObject
)" ] }, "metadata": {}, "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "Figure(PyObject
)" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ @@ -701,17 +676,17 @@ "fig.canvas.draw()\n", "gcf() # Needed for IJulia to plot inline\n", "\n", - "ttext = \"OSE spectra: Re = $(Re), α = $(α)\"\n", - "title(ttext)\n", + "#ttext = \"OSE spectra: Re = $(Re), α = $(α)\"\n", + "#title(ttext)\n", "\n", - "fname = \"test.png\"\n", - "savefig(fname)" + "#fname = \"test.png\"\n", + "#savefig(fname)" ] }, { "cell_type": "code", "execution_count": null, - "id": "026799fc", + "id": "529487ab", "metadata": {}, "outputs": [], "source": [] diff --git a/Lilly1966/OSE_solve_Orszag1971.ipynb b/Lilly1966/OSE_solve_Orszag1971.ipynb index 203db0e..1214239 100644 --- a/Lilly1966/OSE_solve_Orszag1971.ipynb +++ b/Lilly1966/OSE_solve_Orszag1971.ipynb @@ -586,15 +586,15 @@ ], "metadata": { "kernelspec": { - "display_name": "Julia 1.4.1", + "display_name": "Julia 1.6.1", "language": "julia", - "name": "julia-1.4" + "name": "julia-1.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.4.1" + "version": "1.6.1" } }, "nbformat": 4, diff --git a/Lilly1966/OSE_solve_finite.ipynb b/Lilly1966/OSE_solve_finite.ipynb index f58b196..fc0b486 100644 --- a/Lilly1966/OSE_solve_finite.ipynb +++ b/Lilly1966/OSE_solve_finite.ipynb @@ -396,7 +396,7 @@ { "data": { "text/plain": [ - "20×20×20 Array{Complex{Float64},3}:\n", + "20×20×20 Array{ComplexF64, 3}:\n", "[:, :, 1] =\n", " 0.864005+0.749363im 0.804944+0.14258im … -0.866025-0.0236016im\n", " 0.864131-0.0178902im 0.803465-0.0170507im -0.869369-0.0216414im\n",