diff --git a/params/template.ipynb b/params/template.ipynb new file mode 100644 index 0000000..cdc63b8 --- /dev/null +++ b/params/template.ipynb @@ -0,0 +1,3795 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "d007fef7-e788-45c9-a1bd-e37eff9c7e91", + "metadata": {}, + "source": [ + "# Creating paramfiles for the ameriflux PPE\n", + "- 200 latin hypercube samples of a 20-d paramspace" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "81c94223-a9cd-4683-a298-a4d188f86d03", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import xarray as xr\n", + "import os\n", + "import pandas as pd" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "ab3af885-e609-4b76-b87c-fbb1e2a8d9dd", + "metadata": {}, + "outputs": [], + "source": [ + "### import some analysis functions we wrote for this project\n", + "import sys ; sys.path.append(\"..\")\n", + "from ppe_tools import Ensemble,Member,ParamInfo\n", + "from ppe_tools.utils import get_default, parse_val" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "ab067e34-b931-443d-9b00-82c85e17afc5", + "metadata": {}, + "outputs": [], + "source": [ + "#this is the base paramfile that I will copy and edit\n", + "#important that the paramfile matches the version of CTSM\n", + "basefile = '/glade/scratch/djk2120/PPEn11/paramfiles/OAAT0000.nc'\n", + "\n", + "#where to write the new files:\n", + "paramdir = '/glade/scratch/djk2120/PPEn11/amflx_lhc/paramfiles/'\n", + "nldir = '/glade/scratch/djk2120/PPEn11/amflx_lhc/namelist_mods/'\n", + "lhckey = '/glade/scratch/djk2120/PPEn11/amflx_lhc/lhc_220506.txt'\n", + "\n", + "#make the directories, if needed \n", + "for d in [paramdir,nldir]:\n", + " if ~os.path.isdir(d):\n", + " os.system('mkdir -p '+d)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "63b3ad71-8b0f-4608-97c1-46f7f89d34e6", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + " % Total % Received % Xferd Average Speed Time Time Time Current\n", + " Dload Upload Total Spent Left Speed\n", + "100 461 0 461 0 0 1748 0 --:--:-- --:--:-- --:--:-- 1746\n", + "100 6746 0 6746 0 0 7989 0 --:--:-- --:--:-- --:--:-- 7989\n" + ] + } + ], + "source": [ + "#read in the spreadsheet\n", + "# have to publish the spreadsheet\n", + "### file > share > publish to web > csv\n", + "csv = 'may06.csv'\n", + "if not os.path.isfile(csv):\n", + " #data_url needs to be '\"foo.bar\"'\n", + " data_url = '\"https://docs.google.com/spreadsheets/d/e/2PACX-1vQbnPXdKgC4HjYddGl7qOUM6jxTqgqin2uoU2YU5wvbrpbWRuF53PU0k-Dz_Hg445v9yUsROSFEwYDq/pub?gid=1732324005&single=true&output=csv\"'\n", + " cmd = 'curl -L '+data_url+' > '+csv\n", + " os.system(cmd)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "afaaf938-d544-4737-a4a7-17d0dd857db3", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
namelocminmaxpft_minspft_maxs
0act25P40120NaNNaN
1baseflow_scalarN0.00050.1NaNNaN
2cvP0.0050.02NaNNaN
3d_maxP1060NaNNaN
4dleafPpftpft0,0.000216,0.000216,0.00072,0.0081,0.0081,0.00...0,0.00108,0.00108,0.0036,0.0567,0.0567,0.243,0...
5fffP0.025NaNNaN
6frac_sat_soil_dsl_initP0.51NaNNaN
7hksat_sfP0.110NaNNaN
8jmaxb0P0.010.05NaNNaN
9kmaxPpftpft0,3.00E-09,4.00E-09,3.40E-09,1.50E-08,1.30E-08...0,3.80E-08,3.40E-08,3.10E-08,1.70E-07,1.90E-07...
10krmaxPpftpft0,2.30E-11,2.60E-11,2.80E-11,1.10E-10,5.40E-11...0,2.50E-10,2.90E-10,1.60E-10,1.90E-08,5.10E-10...
11leafcnP50percent30percentNaNNaN
12liq_canopy_storage_scalarP0.052NaNNaN
13maximum_leaf_wetted_fractionN0.010.5NaNNaN
14medlyninterceptP15000NaNNaN
15medlynslopePpftpft9,1.53,1.53,1.53,2.91,1.7,2.05,2.05,2.05,2.805...9,4.14,4.14,4.14,9.11,9.11,7.07,7.07,7.07,5.95...
16nstemP0.030.5NaNNaN
17rootprof_betaPpftpft0.886,0.886,0.746,0.746,0.965,0.841,0.965,0.84...0.98,0.98,0.952,0.952,0.994,0.972,0.994,0.972,...
18sucsat_sfP0.0051NaNNaN
19watsat_sfP0.81.5NaNNaN
20zbedrockP14NaNNaN
21wc2wjb0P0.51.5NaNNaN
22jmaxb1N0.050.25NaNNaN
23psi50P50percent50percentNaNNaN
\n", + "
" + ], + "text/plain": [ + " name loc min max \\\n", + "0 act25 P 40 120 \n", + "1 baseflow_scalar N 0.0005 0.1 \n", + "2 cv P 0.005 0.02 \n", + "3 d_max P 10 60 \n", + "4 dleaf P pft pft \n", + "5 fff P 0.02 5 \n", + "6 frac_sat_soil_dsl_init P 0.5 1 \n", + "7 hksat_sf P 0.1 10 \n", + "8 jmaxb0 P 0.01 0.05 \n", + "9 kmax P pft pft \n", + "10 krmax P pft pft \n", + "11 leafcn P 50percent 30percent \n", + "12 liq_canopy_storage_scalar P 0.05 2 \n", + "13 maximum_leaf_wetted_fraction N 0.01 0.5 \n", + "14 medlynintercept P 1 5000 \n", + "15 medlynslope P pft pft \n", + "16 nstem P 0.03 0.5 \n", + "17 rootprof_beta P pft pft \n", + "18 sucsat_sf P 0.005 1 \n", + "19 watsat_sf P 0.8 1.5 \n", + "20 zbedrock P 1 4 \n", + "21 wc2wjb0 P 0.5 1.5 \n", + "22 jmaxb1 N 0.05 0.25 \n", + "23 psi50 P 50percent 50percent \n", + "\n", + " pft_mins \\\n", + "0 NaN \n", + "1 NaN \n", + "2 NaN \n", + "3 NaN \n", + "4 0,0.000216,0.000216,0.00072,0.0081,0.0081,0.00... \n", + "5 NaN \n", + "6 NaN \n", + "7 NaN \n", + "8 NaN \n", + "9 0,3.00E-09,4.00E-09,3.40E-09,1.50E-08,1.30E-08... \n", + "10 0,2.30E-11,2.60E-11,2.80E-11,1.10E-10,5.40E-11... \n", + "11 NaN \n", + "12 NaN \n", + "13 NaN \n", + "14 NaN \n", + "15 9,1.53,1.53,1.53,2.91,1.7,2.05,2.05,2.05,2.805... \n", + "16 NaN \n", + "17 0.886,0.886,0.746,0.746,0.965,0.841,0.965,0.84... \n", + "18 NaN \n", + "19 NaN \n", + "20 NaN \n", + "21 NaN \n", + "22 NaN \n", + "23 NaN \n", + "\n", + " pft_maxs \n", + "0 NaN \n", + "1 NaN \n", + "2 NaN \n", + "3 NaN \n", + "4 0,0.00108,0.00108,0.0036,0.0567,0.0567,0.243,0... \n", + "5 NaN \n", + "6 NaN \n", + "7 NaN \n", + "8 NaN \n", + "9 0,3.80E-08,3.40E-08,3.10E-08,1.70E-07,1.90E-07... \n", + "10 0,2.50E-10,2.90E-10,1.60E-10,1.90E-08,5.10E-10... \n", + "11 NaN \n", + "12 NaN \n", + "13 NaN \n", + "14 NaN \n", + "15 9,4.14,4.14,4.14,9.11,9.11,7.07,7.07,7.07,5.95... \n", + "16 NaN \n", + "17 0.98,0.98,0.952,0.952,0.994,0.972,0.994,0.972,... \n", + "18 NaN \n", + "19 NaN \n", + "20 NaN \n", + "21 NaN \n", + "22 NaN \n", + "23 NaN " + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#proccess csv\n", + "data = pd.read_csv(csv) # modify read_csv to account for header spanning 2 rows\n", + "included = data['final']=='1'\n", + "## make sure to match your columns with the names\n", + "cols = ['name','loc','min','max','pft_mins','pft_maxs']\n", + "params_full = data.loc[included,cols]\n", + "params = params_full.reset_index(drop=True) # reset indexing and get rid of excel row number\n", + "params" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "74fef164-1163-4f56-b49b-8e040c5d948a", + "metadata": {}, + "outputs": [], + "source": [ + "#create the lhc dictionary\n", + "lhcs ={}\n", + "paramlist = params['name'].values\n", + "for p in paramlist:\n", + " ix = params['name']==p\n", + " pinfo={}\n", + " pinfo['loc']=params['loc'][ix].values[0]\n", + " pinfo['flagged']=[] #not using flags here\n", + " for m in ['min','max']:\n", + " mval = params[m][ix].values[0]\n", + " if mval=='pft':\n", + " pftval = params['pft_'+m+'s'][ix].values[0]\n", + " formatted = np.fromstring(pftval, dtype='float', sep=',')\n", + " elif 'percent' in mval:\n", + " formatted = mval\n", + " else:\n", + " formatted = np.array(float(mval)) \n", + "\n", + " pinfo[m]=formatted\n", + "\n", + " lhcs[p] = pinfo\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "29dec9d3-c4bc-4f3f-b0e4-01ab4226de8b", + "metadata": {}, + "outputs": [], + "source": [ + "#instantiate the Ensemble object\n", + "x = Ensemble(basefile,paramdir,nldir)\n", + "\n", + "#add the new OAATS\n", + "prefix = 'LHC'\n", + "nextnum=1\n", + "n_samples=250\n", + "x.add_lhcs(lhcs,prefix,nextnum,n_samples)\n", + "\n", + "#write the ensemble\n", + "x.write(lhcfile=lhckey,default_key='LHC0000')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4d1bdfab-66e7-491d-b608-5f2356922cea", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b2a705d2-929a-4406-b8a5-5fba5f62f9a3", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 44, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAaRUlEQVR4nO3df5Rc5X3f8fdnRzI2NQbFEMzRDySIAKm0BHktmjo4isEEiI2K0rQY45pfR4fWpHaTnBbjY5+0pz4oqUmic8oplYliE1sI20DMH3bAjt1SWoK0+mUEiCALGy0iIBw3copB2tlv/5g7cDWa2ZndvXPvnTuf1zl72Ln3zs6zd4aPvvs8z72PIgIzM6uukaIbYGZm/eWgNzOrOAe9mVnFOejNzCrOQW9mVnFzim5AOyeffHIsXry46GaYmQ2Mbdu2vRIRp7TbV8qgX7x4MWNjY0U3w8xsYEj6Uad97roxM6s4B72ZWcU56M3MKs5Bb2ZWcQ56M7OK6ynoJV0q6RlJeyXd0mb/PEkPSPq+pC2Szk3t+3eSnpS0W9I9kt6a5S9gZmZT6xr0kmrAHcBlwHLgw5KWtxx2K7AzIv4x8K+A9clz5wP/FhiNiHOBGnBVds03M6uGPVu/w2NfupU9W7+T+c/uZR79SmBvROwDkLQZWA08lTpmOXAbQETskbRY0qmp13ibpCPA8cCBrBpvZlYFj3/tdlbs/hwjTHJ43xfYwz2c856LM/v5vXTdzAf2px6PJ9vSdgFrACStBE4HFkTEC8DngeeBF4G/i4iH272IpLWSxiSNHTx4cHq/hZnZoNm/hYP3/Gt23/4h3r37PzOHOjUFcznCT576bqYv1UtFrzbbWlcrWQesl7QTeALYAUxImkej+l8C/F/ga5KuiYgvH/MDIzYAGwBGR0e9GoqZVc/+LbBrE/z9QerP/AUnT05wcrJLgggIRpi3/P2ZvmwvQT8OLEw9XkBL90tEHAKuA5Ak4Lnk69eA5yLiYLLvfuCfAscEvZlZZTUDfvuXickjAIxEI9yhEfD1gElqbD/301yQYbcN9Bb0W4GlkpYAL9AYTL06fYCkk4BXI+IwcCPwSEQckvQ88E8kHQ/8DLgI8E1szKz6UtU7z34b6ocJ4o0ukqAR8ABHqLHjnR/ixF/6aOYhDz0EfURMSLoZeIjGrJmNEfGkpJuS/XcCy4C7JdVpDNLekOx7XNLXge3ABI0unQ2Z/xZmZmWRqt5Jqvcg6QOPN/u9j1Dju/Xz+bFO5JT3Xssll17RtyapjIuDj46Ohu9eaWYDY/8W+OH/gtcOwWP/FSbrpIcyI2ASMUGN79V/kVc4kT+fvJCz33Mxa1Ys4N2nz5t1EyRti4jRdvtKeZtiM7OBMfZF+ObvHBvuvPnwCDW+Wl/F/fUL2R5nMWdE/KfV53L1BYtyaaKD3sxsutL97898C6J+1O4A6jHCd+oreIUTub9+ITviLObWxEdGF2ZWxffKQW9m1qs2/e9pAUwms2c+c+RaNk9eBEBNcPUFi3IP+CYHvZnZVLr0v5M8qlPjCxOX89M4nr+aXMb2OAsg926adhz0ZmaddOh/T6trDpsnfoX7Ji58I9yhUcVftbK4Kj7NQW9mltal/x0A1fjbhRfx139/PJ9/6XzG6kuP2l2GKj7NQW9mBl373xsEIzUeX3YrH9mxjPpkHFXnl6mKT3PQm9nw6qH/HYCROfBLN/PCa3N58O/O5PM7TqQ+efRxZavi0xz0ZjZ8mtX7jk1QPwJMtj9uZC6s+Cic92E2HXgXn/3e7qSKfzPky1rFpznozWw49Fq9qwZnXwZv/3k478Nsm1zKfWPj3Lt190BV8WkOejOrvh5mzzT737n8dhi9FoBNjz/PZ7/x2MD0xXfioDezaupx9gzv/hi86zz42Y9h8YWwcCXbfvQT7ts+zr1b9w9sFZ/moDezapnG7Jl09d7UqOJ3D3wVn+agN7PBN83ZM7z1HW9U701Vq+LTHPRmNth66X9PzZ5Jh3tTFav4NAe9mQ2mZhfNtrs797+nZs+0C/gqV/FpDnozGxw9ddF07n9vagb817eNc2RispJVfJqD3swGQ7cumin639M6ddNAtar4NAe9mZXXdKZIduieaerUTSNgbk38ZgELguTFQW9m5TPLKZKtqj7Y2o2D3szK46h70BxmNl00MDyDrd046M2sWD3PgZ96imSrYa/i0xz0ZlacKQdYBbW5sPSSKadItnIVfywHvZnlK8MB1lau4ttz0JtZPjIeYE1zFT81B72Z9V9Gc+DbcRXfnYPezPqjly6aaQ6wprmK752D3syy1a2Lpod70HTjKn56HPRmlp1eumim2f+e5ip+Zhz0ZjY7fe6iaXIVP3MOejObmRy6aMBVfBYc9GY2fX3uomlyFZ8NB72Zdde8TcHb3gl/s7PzYh8ZdNGAq/is9RT0ki4F1gM14K6IWNeyfx6wETgTeA24PiJ2J/tOAu4CzqXxT//1EfFYVr+AmfXZMdW7OKqKz6iLBoZvQZC8dA16STXgDuADwDiwVdKDEfFU6rBbgZ0RcaWkc5LjL0r2rQf+IiL+uaS3AMdn+huYWX90XKovFb8ZddHAcC4IkpdeKvqVwN6I2AcgaTOwGkgH/XLgNoCI2CNpsaRTgZ8B7wOuTfYdBg5n1noz64+p+uA10gj486+ZdQUPw70gSF56Cfr5wP7U43HggpZjdgFrgEclrQROBxYAdeAg8KeSzgO2AZ+IiP/X+iKS1gJrARYt8r/aZrnrNk1yFrcp6MSDrfnoJejVZlvrX1brgPWSdgJPADuACWAusAL4rYh4XNJ64BbgM8f8wIgNwAaA0dHRNsP4ZtYXvUyTnMGdJKfiwdZ89RL048DC1OMFwIH0ARFxCLgOQJKA55Kv44HxiHg8OfTrNILezMogp2mSaa7i89dL0G8FlkpaArwAXAVcnT4gmVnzatIHfyPwSBL+hyTtl3R2RDxDY4D2KcysWB0HWhMZTZNMcxVfnK5BHxETkm4GHqIxvXJjRDwp6aZk/53AMuBuSXUaQX5D6kf8FvCVZMbNPpLK38xy1m3JvgynSbZyFV8sRZSvO3x0dDTGxsaKboZZdRTQRQOu4vMkaVtEjLbb5ytjzaqsWxdNHwZam1zFl4eD3qyqui28PYMl+3rhKr58HPRmVVLAXPg0V/Hl5KA3q4IC5sKnuYovNwe92SBrBvyOTVA/TJ4DrU2u4svPQW82qLrOpMl+Lnyaq/jB4aA3GzQdZ9IIanNh6SV9mQuf5ip+sDjozQZBLxc79bEPvslV/GBy0JuVXUEXO6V5QZDB5qA3K6sCL3ZK84Igg89Bb1ZGBV3slOYFQarDQW9WJlNV8X2+2CnNg63V4qA3K1pJBlrBg61V5aA3K1IJBlqbXMVXl4PerAglGWgFV/HDwEFvlrcSDLQ2uYofDg56s7yUZKAVXMUPGwe9WR46VfE5dtE0uYofPg56s37qVsXn1EUDruKHmYPeLGvpxT+e/faxtw92FW85c9CbZaXb4h/gKt4K4aA3y0K3mTS1uXD+Na7irRAOerPZ6DYfvs+Lf7TjKt5aOejNZmqqmTRnX9b3xT/acRVv7TjozaarRDNpmlzF21Qc9GbTUaL58OAFQaw3DnqzXpSwiveCINYrB71ZNyWt4r0giPXKQW/WyQBV8e6msak46M3aGZAqHtxNY9056M2a0rcueOZbruKtMhz0Zt1uXeAq3gacg96GW4mW8ktzFW9Z6inoJV0KrAdqwF0Rsa5l/zxgI3Am8BpwfUTsTu2vAWPACxHxwYzabjZzJbx1AbiKt/7oGvRJSN8BfAAYB7ZKejAinkoddiuwMyKulHROcvxFqf2fAJ4G3pFZy81mohnwOza1v31wQbcuAFfx1j+9VPQrgb0RsQ9A0mZgNZAO+uXAbQARsUfSYkmnRsRLkhYAvw58DvjtTFtvNh1TddMU1EUDruKt/3oJ+vnA/tTjceCClmN2AWuARyWtBE4HFgAvAX8M/HvghKleRNJaYC3AokX+UFuGOnbTFHP74DRX8ZaHXoJebba1jlqtA9ZL2gk8AewAJiR9EHg5IrZJWjXVi0TEBmADwOjoaJtRMbMZKNl8+CZX8ZanXoJ+HFiYerwAOJA+ICIOAdcBSBLwXPJ1FXCFpMuBtwLvkPTliLgmg7abdVbCq1qbXMVb3noJ+q3AUklLgBdohPfV6QMknQS8GhGHgRuBR5Lw/1TyRVLR/65D3vrOVbzZUboGfURMSLoZeIjG9MqNEfGkpJuS/XcCy4C7JdVpDNLe0Mc2m7XnKt6sLUWUrzt8dHQ0xsbGim6GDRJX8TbkJG2LiNF2+3xlrA22klbxXhDEysRBb4OrpFW8FwSxsnHQ2+ApeRXvBUGsbBz0NlgGrIp3N42VgYPeBsOAVfHgbhorDwe9lZ+reLNZcdBbebmKN8uEg97KyVW8WWYc9FYuruLNMuegt3LotiCIq3izGXPQW/G8IIhZXznorTheEMQsFw56K0ZJB1tdxVsVOegtf/u3JCE/cfR230rYrC8c9Jav/Vvgf9yWVPIJV/FmfeWgt/wc010jGKm5ijfrMwe99V/bQVfBmb8Kqz7lKt6szxz01l+dBl1HaoWEvBcEsWHkoLf+6OUK15xD3guC2LBy0Fv2SjR1slnBv/LT1/nLPS97QRAbSg56y06J7lPTDPivju1not5av7ubxoaLg96yUaIqfqouGnA3jQ0fB73NTgmr+HYzaQDm1MS/dDeNDSEHvc3cAFTxNcFFy07llBOOc8Db0HLQ28yU5DYGng9v1p2D3qavJLcx8FWtZr1x0Nv0lOA2Bq7izabHQW+9KcltDFzFm02fg96mNtUSfznexsBVvNnMOeits16W+Msh5F3Fm82Og96OVZIl/lzFm2XDQW9HK8nceFfxZtlx0NubSjA33lW8WfZ6CnpJlwLrgRpwV0Ssa9k/D9gInAm8BlwfEbslLQTuBt4FTAIbImJ9hu23rBQ8N973iTfrn65BL6kG3AF8ABgHtkp6MCKeSh12K7AzIq6UdE5y/EXABPA7EbFd0gnANknfbnmuFa3gufG+T7xZf/VS0a8E9kbEPgBJm4HVQDqslwO3AUTEHkmLJZ0aES8CLybbfyrpaWB+y3OtKAXPje/UTeP7xJtlq5egnw/sTz0eBy5oOWYXsAZ4VNJK4HRgAfBS8wBJi4HzgcfbvYiktcBagEWLXL31XcFL/Hmw1Sw/vQS92mxr/Qt7HbBe0k7gCWAHjW6bxg+Q3g7cB3wyIg61e5GI2ABsABgdHW13G3HLQsFL/Hmw1Sx/vQT9OLAw9XgBcCB9QBLe1wFIEvBc8oWkuTRC/isRcX8GbbaZ2r8FvnQFTLxGEVMnXcWbFaOXoN8KLJW0BHgBuAq4On2ApJOAVyPiMHAj8EhEHEpC/0+ApyPiDzNtuU3frnuODfkcpk66ijcrVtegj4gJSTcDD9GYXrkxIp6UdFOy/05gGXC3pDqNgdYbkqe/F/go8ETSrQNwa0R8M9tfw6bU7K7Z/mXeCPmRubDio67izYZAT/Pok2D+Zsu2O1PfPwYsbfO8R2nfx295aTvoKlhxDXzwj/r2sq7izcrDV8ZW1VT3q5nzVjjv6o5PnS1X8Wbl4qCvooIGXV3Fm5WTg76KChh0dRVvVl4O+iopYNDVVbxZ+Tnoq6KAQVdX8WaDwUE/6AoYdHUVbzZYHPSDrIBBV1fxZoPHQT+omvePr79OHoOuruLNBpeDfhC19sdrpBHwfVjL1QuCmA0+B/0g6XT/+DNW9eXWwl4QxKwaHPSDIsf7x3tBELNqcdAPgm6LdmcY8h5sNaseB33Z5bRotwdbzarLQV9mOS3a7SrerNoc9GWU06LdruLNhoODvmw6XQSV8aCrq3iz4eGgL5up7jyZQci7ijcbPg76ssjhzpOu4s2Gk4O+DPp850lX8WbDzUFftLZz5LO786SreDNz0Bepj3PkXcWbWZODvih9nCPvKt7M0hz0eevjHHlX8WbWjoM+T32cI+8q3sw6cdDnqQ9z5F3Fm1k3Dvo89GGOvBcEMbNeOej7rW13zezmyHtBEDObDgd9vx3TXTPzOfJeEMTMZsJB3y8Zd9d4sNXMZspB3w8Zdtd4sNXMZstBn7Xm1a7115ltd42reDPLgoM+S61Xu2qkMX3y/Gum1V3jKt7MsuSgz0Knq13PWDXtC6FcxZtZ1noKekmXAuuBGnBXRKxr2T8P2AicCbwGXB8Ru3t57sDL6GpXV/Fm1i9dg15SDbgD+AAwDmyV9GBEPJU67FZgZ0RcKemc5PiLenzuYMvgaldX8WbWT71U9CuBvRGxD0DSZmA1kA7r5cBtABGxR9JiSacCZ/Tw3MGUwfRJV/Fmlodegn4+sD/1eBy4oOWYXcAa4FFJK4HTgQU9PhcASWuBtQCLFpU83DKYPukq3szy0kvQq8221ivv1wHrJe0EngB2ABM9PrexMWIDsAFgdHS07TGlMYurXV3Fm1neegn6cWBh6vEC4ED6gIg4BFwHIEnAc8nX8d2eO1Bm2V3jKt7MitBL0G8FlkpaArwAXAUcVbpKOgl4NSIOAzcCj0TEIUldnzswZtFd4yrezIrUNegjYkLSzcBDNKZIboyIJyXdlOy/E1gG3C2pTmOg9YapntufX6XPZthd4yrezIqmiPJ1h4+OjsbY2FjRzWhId9dMHmls66G7xlW8meVJ0raIGG23z1fGTmUG3TVeEMTMysZBP5Vpdtd4QRAzKyMHfTvTnF3jBUHMrMwc9K2m2V3jwVYzKzsHfaseu2s82Gpmg8JBn7Z/C+zo3l3jKt7MBomDvqm5MtTkRLLh2O4aV/FmNogc9NB+ZajacUd117iKN7NB5aDfvyUJ+VQln1oZylW8mQ264Q76N7pr6m9uS60M5SrezKpgOIO+OU9+xyaoH6Ex+KpGyF9+O9sml3LfA0+4ijezShi+oG87T34EzlwFqz7FpgPv4rP//TFX8WZWGcMX9G3nyR/HnnM+zt1jb+PerbtdxZtZpQxX0HeYJ//wnF/l3zxwmPrk867izaxyhifo28yTP7j0N/mjI9dz7/9xX7yZVddwBH3LPPlghCOay8efPIetE67izazaqh/0LfPkJxGP1v8hfzzxG2yPXzjqUFfxZlZF1Q76lnnyAdRjJAn5swDfStjMqq+6Qf/GNMrXCYLJgElqfObItW+EvLtpzGwYVDfod91DTLyGCOoh/vfkuaxPVfLupjGzYVG9oE+ueq1v+zNGonHRU53aGyHvKt7Mhk21gn7/Fia/+CGov85IBBLUQ3yt/itsj7NcxZvZUKpU0D/78AbOnHiNETUHXsVh5vLnkxfykQtcxZvZcKpM0O/Z+h2WPP8AAiLgCDW+Wl/FN+J9rFm9xlW8mQ2tygT9T576Lr/AZKq7ZhVPv/s/coureDMbciNFNyAr85a/nyPMYSJGOMxcTvnla/nclf/IIW9mQ68yFf0577mYPdzDT576LvOWv59L3nNx0U0yMyuFygQ9NMIeB7yZ2VEq03VjZmbtOejNzCrOQW9mVnEOejOzinPQm5lVnIPezKziFBHdj8qZpIPAj2b49JOBVzJsTlbcrukra9vcrulxu6ZvJm07PSJOabejlEE/G5LGImK06Ha0crumr6xtc7umx+2avqzb5q4bM7OKc9CbmVVcFYN+Q9EN6MDtmr6yts3tmh63a/oybVvl+ujNzOxoVazozcwsxUFvZlZxlQl6SZdKekbSXkm3FNiOhZK+J+lpSU9K+kSy/fckvSBpZ/J1eUHt+6GkJ5I2jCXbfk7StyU9m/w319VaJJ2dOi87JR2S9MkizpmkjZJelrQ7ta3j+ZH0qeQz94ykXyugbf9F0h5J35f0gKSTku2LJf0sde7uzLldHd+7vM5Zh3bdm2rTDyXtTLbneb46ZUT/PmcRMfBfQA34AXAG8BZgF7C8oLacBqxIvj8B+GtgOfB7wO+W4Fz9EDi5ZdsfALck398C/H7B7+XfAKcXcc6A9wErgN3dzk/yvu4CjgOWJJ/BWs5tuwSYk3z/+6m2LU4fV8A5a/ve5XnO2rWrZf/twGcLOF+dMqJvn7OqVPQrgb0RsS8iDgObgdVFNCQiXoyI7cn3PwWeBuYX0ZZpWA18Kfn+S8A/K64pXAT8ICJmemX0rETEI8DftmzudH5WA5sj4vWIeA7YS+OzmFvbIuLhiJhIHv4VsKBfrz+ddk0ht3M2VbskCfgXwD39eO2pTJERffucVSXo5wP7U4/HKUG4SloMnA88nmy6OfkTe2Pe3SMpATwsaZuktcm2UyPiRWh8CIGfL6htAFdx9P98ZThnnc5P2T531wPfSj1eImmHpP8p6cIC2tPuvSvLObsQeCkink1ty/18tWRE3z5nVQl6tdlW6LxRSW8H7gM+GRGHgP8GnAn8IvAijT8bi/DeiFgBXAZ8XNL7CmrHMSS9BbgC+FqyqSznrJPSfO4kfRqYAL6SbHoRWBQR5wO/DWyS9I4cm9TpvSvLOfswRxcUuZ+vNhnR8dA226Z1zqoS9OPAwtTjBcCBgtqCpLk03sCvRMT9ABHxUkTUI2IS+AJ9/BN/KhFxIPnvy8ADSTteknRa0vbTgJeLaBuNf3y2R8RLSRtLcc7ofH5K8bmT9DHgg8BHIunUTf7M/3Hy/TYa/bpn5dWmKd67ws+ZpDnAGuDe5ra8z1e7jKCPn7OqBP1WYKmkJUlVeBXwYBENSfr+/gR4OiL+MLX9tNRhVwK7W5+bQ9v+gaQTmt/TGMjbTeNcfSw57GPAN/JuW+KoKqsM5yzR6fw8CFwl6ThJS4ClwJY8GybpUuA/AFdExKup7adIqiXfn5G0bV+O7er03hV+zoCLgT0RMd7ckOf56pQR9PNzlscoc04j2ZfTGL3+AfDpAtvxyzT+rPo+sDP5uhz4M+CJZPuDwGkFtO0MGqP3u4Anm+cJeCfwl8CzyX9/roC2HQ/8GDgxtS33c0bjH5oXgSM0Kqkbpjo/wKeTz9wzwGUFtG0vjf7b5mftzuTY30je413AduBDOber43uX1zlr165k+xeBm1qOzfN8dcqIvn3OfAsEM7OKq0rXjZmZdeCgNzOrOAe9mVnFOejNzCrOQW9mVnEOejOzinPQm5lV3P8H0cYE7UhJib0AAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "rpb0 = 0.886\n", + "rpb1 = 0.98\n", + "d0 = np.log(0.05)/np.log(rpb0)\n", + "d1 = np.log(0.05)/np.log(rpb1)\n", + "drpb = rpb1-rpb0\n", + "rpbs2 = rpb0+np.log(np.linspace(1,np.exp(1),200))*drpb\n", + "d95s2 = np.log(0.05)/np.log(rpbs)\n", + "plt.plot(rpbs,'.')\n", + "plt.plot(rpbs2,'.')\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 82, + "id": "f011346b-d9c4-42e6-9746-89e8d531eb6e", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 82, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAY50lEQVR4nO3dfYxddZ3H8c/3zu00sHZhpFUe2k6pIpGC0c6E1oBmVXSB7dIFV3lwNyYuQRNMJK7xidiYbjRxjS5/rFm2i2bNLkhRRAzRLOLDEgitnakg7WKhjJ12KMK0jEBSwjzc7/5x7xnvDPfM0z33nPM75/1KJp05c3vuj3NvP/zu9/dwzN0FAAhXJesGAADaQ5ADQOAIcgAIHEEOAIEjyAEgcNUsnnTlypW+bt26LJ4aAII1ODh4zN1XzT6eSZCvW7dOAwMDWTw1AATLzIZbHae0AgCBI8gBIHAEOQAEjiAHgMAR5AAQOIIcAAKXyfTDkAwOj2nX0HH1nNytfUdflEnacOYpGjsxrs3rT1Nfb0/WTQRQcgR5C1F4v/zKhG576Peaqrlmb/ZrkroqpusvPlsrTlpGqAPIDEE+yx27D2vbvftahnczlzRZc9364NB0qG/fer6u27Q2pZYCQB1B3jA4PKa7945o554jmqot7mYbUah/6UePa//RF3XVxtX0zgGkhiBXPcQ/ctsuvTpRm7OE8tKrkzJJK5ZXW5Zcply6ffdh7dxzhN45gNQQ5JJ+uHfkNSFenaf+/f4Np8fW0emdA0iTZXHPzv7+fs/DpllROeWugSOanKpfh2qX6er+NYsK4LiyjEla1mX60CLPBwCtmNmgu/e/5nhZg7xVOcUkXbtprb565QVLOudcA6VVBkMBtCkuyEtbWpldTjFJy5dV9MGNq5d8zus2rdW5p6/Q3XtH9IPBEU1M1ii3AOi4UvbIB4fHdO2ORzTeRjllIc8RV25Zvqyi26/fTJgDWBR65E12DR3XZCNcTdKH+9foK0ssp8Tp6+1RX2+Pzj/zlBnlFpf06kRNd+8dIcgBJKJ0QT44PKZn/viKql0VTU3VtKzaXjllPs3llmhQ1SXdNXBEJlFmAdC2UgV5NMA5PllTtWK65sK1qQRp1Ds31QdEXdLklOuO3Yd1994RyiwA2lKq3Q93DR3X+GRNNZemaq4zTz0p1QC9auNqLV9WkTV+bi6zAMBSlSbIm0sqXSYtq1a0ef1pqbahr7dHt1+/WdduWqtqVz3OozLLzfc8rsHhsVTbA6AYSlFayaqk0gplFgBJK0WPPOuSSiuUWQAkpRRB3nNytypmqmRUUmklrszyg8ERSiwAFqXwQT44PKbt9+3XVM1VMdO2LRsy741H+np79NUrL9DV/Wume+aTkzXd8sCThDmABSt8kEdlFZfk7ho7MZ51k14jKrNUJNUkPfTUMV3974/ojt2Hs24agAAUPsg3rz9N3dXsZqosRFRmueiclTL96UYV2+7dR88cwLwKH+RRSH76A+fmekZIX2+PbrrkLeqq2PSxqZpTZgEwr0IH+eDwmL71y4OSpBvf8+bchnikr7dH27eer2rFpnvmDx88po/ctoswBxCrsPPIm+eOd1fD2W0w2pvllgee1MMHj6nmbLIFYG6J9cjNrMvMfmNm9yV1znY0zx2fmKxp19DxrJu0YFGZpVphWiKA+SVZWvmUpCcSPF9bQhjknEtfb48+xLREAAuQSJCb2WpJfyXptiTOl4RQBjnnMntaIvVyAK0k1SO/RdJnVc+blszsBjMbMLOB0dHRhJ423uDwmHYNHdfm9acFGeLSzGmJFdOMejkARNoOcjPbIul5dx+c63HuvsPd+929f9WqVe0+7Zyigc5v3H8g+B4s9XIA80miR36RpCvM7JCkOyW918z+O4HzLlnIA52tUC8HMJe2g9zdv+Duq919naRrJP3C3f+u7Za1IfSBzlaolwOIU8gFQUUY6JyNejmAOIkGubv/yt23JHnOxQptNediUC8H0EqheuRFGuSMQ70cwGyFCvKiDXLGoV4OoFmhgryIg5yttKqXj0/QMwfKytw99Sft7+/3gYGBjpy7CAuBFmp6Y7CJmmqSKqagNggDsDhmNuju/bOPF6pHLtV7q0Ub5IzDTBYAUgGDvGyYyQKgUEEeTT0sW4gxkwUot8IEeRmmHs6FmSxAeRUmyMsy9TAO9XKgvAoT5GWZejgX6uVAORUmyIu4v8pSUC8HyqcwQS6Va+rhXKiXA+VSqCBHXat6eRnHDYCyIMgLKqqXd1frPXMzU8/J3Vk3C0AHFCLIyzp/fD59vT3atmWDKhVTzV3b79vPNQIKKPggL/v88fmMnRhXzZ2NtYACCz7Iyz5/fD7RtEwGPoHiCj7ImT8+NxYKAcUXfJAzf3x+LBQCii34IJeYP74QLBQCiqsQQY6FYaEQUEwEeYmwUAgoJoK8ZFgoBBQPQV5CLBQCiiX4IGdV59KwUAgojmrWDWjH9F3kJ2vcPX6Rovn34xO16YHPPYde4BoCAQq6R86qzqVjoRBQHEEHOas628NCIaAYgg5yVnW2b/ZCoakpPtkAoQk6yCVWdSYhWijUZVJXxXT0j6/QKwcC0naQm9kaM/ulmT1hZvvN7FNJNAzpiT7ZXH3hWslM3/v1YVZ8AgFJokc+Kekf3f2tkjZLutHMzkvgvEhRX2+Pzjr1JE1O1Rj4BALTdpC7+7Puvrfx/cuSnpB0VrvnRfo2rz+NgU8gQInWyM1snaR3SNrd4nc3mNmAmQ2Mjo4m+bRICDskAmFKLMjN7HWS7pZ0k7u/NPv37r7D3fvdvX/VqlVJPS0Sxg6JQHgSCXIzW6Z6iN/u7j9M4pzIBguFgPAkMWvFJH1b0hPu/s32m7Qw7LHSOSwUAsKSRI/8Ikl/L+m9ZvZo4+vyBM4bK9pj5Rv3H+Bjf4ewUAgIRxKzVh5yd3P3t7n72xtfP0micXHYYyUdLBQCwhDkyk72WEkHC4WAMAQZ5Oyxkp7ZC4XYuxzIn2D3I+/r7SHAU8Le5UC+BdkjR7qYkgjkG0GOBWFKIpBfBDkWjCmJQD4R5FgUpiQC+UOQY1GYkgjkD0GORWPvciBfCHIsCXuXA/lBkGNJ2LscyA+CHEvG3uVAPgQZ5Gxhmw+tFgqxiRmQvuCW6Edb2I5P1tRdrbBMPGPRQqE9h17Q+ERNZqaek7uzbhZQKsH1yNnCNn/6enu0bcsGVSqmmru237efT0tAioILcrawzaexE+OqubNDIpCB4EorUV1219BxbV5/GmWVnGCHRCA7wfXIpXqY3/ieNxMQOcIOiUB2ggxy5BM7JALZIMiRKHZIBNJHkCNx7JAIpIsgR+LYIRFIF0GOjmCHRCA9BDk6hh0SgXQQ5OgYdkgE0kGQo6PYIRHoPIIcHcUOiUDnEeTouGihUHe13jNnh0QgWQQ5UsEOiUDnEORIDTskAp2RSJCb2aVmdsDMDprZ55M4J4on2iGRgU8gWW0HuZl1SfqWpMsknSfpWjM7r93ztsIt3sLGDokou05lWBL7kV8o6aC7D0mSmd0paauk/0vg3NO4xVsxRAOfu4eOa3zKpxcKfXDjal5PFFonMyyJ0spZko40/TzSODaDmd1gZgNmNjA6OrroJ+EWb8XBQiGUUSczLIkgtxbH/DUH3He4e7+7969atWrRT8It3oqFhUIom05mWBKllRFJa5p+Xi3paALnnYFbvBVL9Hre8sCTevjgsRm9FF5bFM3g8Jh2DR3Xti0bNHZiPPEMSyLI90g6x8zOlvSMpGskXZfAeV+jr7eHf+QFEtXL9xx6QeMTNRYKoZDSGN9ru7Ti7pOSPinpfyQ9Iekud9/f7nlRDiwUQtGlMb6XyDxyd/+Ju7/F3d/k7l9J4pwoDxYKocjSGN9LorQCtCV6o49P1KYHPvcceoEppghep2vjEZboI3OtFgrRM0footr4N+4/oO337e/oJA2CHLkwe4dEpiQidGmufSHIkRss4UeR9JzcrYqZKimsfSHIkStRz5x7fSJkg8Nj2n7ffk3VXBUzbduyoaPjPQQ5cocl/AhdVFZxSe6usRPjHX0+ghy5xBJ+hCztLUWYfohcYgk/QpXWlMNmBDlyiyX8CE1W221TWkGusYQfIclqu22CHLnHEn6EIqvttimtIPdYwo+QXLVxtazxZ1rvT3rkyD0WCiEEUX38zl8fTv29SZAjCCwUQt5leTtKghzBYKEQ8mpweEzP/PEVVbuyuR0lQY6gsFAIedNcUpG7rrlwberjNwQ5gkK9HHnTXFKZqrnOPPWk1AfhCXIEh3o58iTNXQ7jEOQIEvVy5EHauxzGIcgRLOrlyFrauxzGIcgRLOrlyFLWM1WaEeQIGvVyZCEPM1WaEeQIHvVypC0PM1WaEeQoBOrlSFNWm2PFIchRCNTLkbarNq7OvKQSIchRGNTLkYYsN8eKQ5CjUKiXo9Oy3BwrDkGOwqFejk7J05TDZgQ5Cod6OTohb1MOmxHkKCTq5Uha3qYcNmsryM3s62b2OzP7rZndY2anJtQuoG3Uy5GkPGyOFafdHvnPJJ3v7m+T9KSkL7TfJCA51MuRhLxsjhWnrSB39/vdfbLx4y5Jq9tvEpAc6uVo1+DwmG554MlcbI4VJ8ka+cck/TTul2Z2g5kNmNnA6Ohogk8LzI16OZYqGuB86KljqrlyWVaRFhDkZvaAme1r8bW16TE3S5qUdHvcedx9h7v3u3v/qlWrkmk9sEDUy7EUzdvUViRd9OaVuZmp0qw63wPc/ZK5fm9mH5W0RdL73N2TahiQtKs2rtbde0c0PlGbrpfvOfRCLv9hIh+iAU7J1V2t6KZL3pLL90q7s1YulfQ5SVe4+4lkmgR0Rqt6+fgEPXO0lvcBzmbt1sj/VdIKST8zs0fN7NYE2gR0TFQv764ykwVzy8vdfxai3Vkrb3b3Ne7+9sbXJ5JqGNApzGTBfPK6FD8OKztRSsxkQZw8L8WPQ5CjtJjJglbyvBQ/DkGOUmPlJ2bL81L8OAQ5So16OZqFNFOlGUGO0mtVL79r4IhuvudxeuYlEsJS/DgEOaAW9fIp1x27D1NmKYlQluLHIciBhqheHoW5izJLWYSyFD8OQQ40RPXyazetVbWLaYll0jzA2b0sv0vx4xDkQJO+3h599coLdDXTEksj1AHOZgQ50ALTEssh5AHOZgQ50ALTEosv9AHOZvNuYwuUVTQtcffQcY1P+fS0RFO9xx7ax2/8yeyeeDTAGVptPEKPHJgD0xKLp1VPPMQBzmYEOTAPpiUWyw/3jujViXCnGrZCkAPzYFpicQwOj+n7A0cU3cqsmuO7/iwGQQ4sANMSi2HX0HFN1uoxbpL+tq8YYx0EObAIs6clPvTUMV3974/ojt2Hs24a5jH7ZhHLl1X0wY2rs25WIghyYBGapyXWb8krTdZc2+7dR888x0K8WcRiEOTAIkXTErsqNn1squaUWXIsGuAM6WYRi0GQA0vQ19uj7VvPV7Vi0z1zyiz5NHuAs6srzEU/cyHIgSW6btNa7fz4O3UxZZbcihb+FHGAsxlBDrSBMkt+tVr4U6QBzmYEOdAmyiz5VMSFP3EIciABlFnyY3B4TF+853HtLODCnzgEOZAQyizZi8op39t9WJNTxa6LNyPIgQS1KrM8fPCYrt3xCDdzTkFzOUWqh3hR6+LNCHIgYc1llmgv8/Ep1+27D1M376DX7KPSZbpuU7EW/sQhyIEOiMos3dU/7ZooUTfvlFbTDD/cv0ZfufKCwoe4RJADHdO8ayJ18865o/FJpwzTDOMQ5EAHRbsm/hPTEzticHhM2+7dp8la/Q5OpmJPM4yTSJCb2WfMzM1sZRLnA4ombnril370OIOgSxSVU6ZqPn2sq2KFnmYYp+0gN7M1kt4via4FMIeW0xNdDIIuQfOqzagnXq2Ytm89v3QhLiXTI/8XSZ+V5PM9ECi72dMTIwyCLs7sVZsXn7NSOz/+Tl23aW3WTctEW0FuZldIesbdH1vAY28wswEzGxgdHW3naYGgRWUWBkEXr4yrNheiOt8DzOwBSae3+NXNkr4o6QMLeSJ33yFphyT19/fTe0ep9fX2qK+3R+efeYq23btPU43BuoeeOqZHnj6u7VvPL23vMs4duw/PuFZSOVZtLsS8Qe7ul7Q6bmYXSDpb0mNmJkmrJe01swvd/Q+JthIoqOs2rdW5p6/QLQ88OV3vjQZB9x99UVdtJKSkmbNTImVZtbkQSy6tuPvj7v4Gd1/n7uskjUjaSIgDi8Mg6Nxazk4xlWbV5kLM2yMH0HnRIOjs0kHZe+ezyymm+hRDSk8zmXv65er+/n4fGBhI/XmBvBscHtPde0e0c8+RGT1Q6U/T68oQYK2ug6k+O6XMA5tmNuju/bOP0yMHciRuEFQqR+88CvAfDI5oYrI2Y05zWRf7LARBDuRQNAg6u1ca1c537jlSuN55q1kpkTIv9lkIghzIqYX0zn914HmtWrE86B56XDnJJC3rMn2of03Q/31poEYOBGCu2rlU33v76gADL64X3mXSNReuDe6/p9PiauQEORCQucoPUQ/2L859Q6576dH/lI69/Kp+/rvnSz2ou1gEOVAQURDeNXBk+r6UreStlz5fu+mFz48gBwqmuWf7qydHXzPLI1KtmK6/+GytOGmZNq8/LdWQHBwe066h43r5lQnd9tDvW36SiNpIL3x+TD8ECiYaDJXm7u1O1ly3Pjg0vZim06G+0PCW8vepIVT0yIECma/+HGkO9ZdenZRJ2nDmKRo7Mb7ogI+e0yStWF6dN7y7THrfW9+Y6zp+XtEjB0qguZc+18BotDnXrQ8OzTgeF/D7jr44I+x7Tu7WvqMv6tjLr+oXB56fs1Y/+9yUUJJHkAMFFS0qWmiZQ4oP+GbRreoWIq1yTtkR5ECBNffQ37/h9EWFepz5/g7hnT6CHCiJVqEelUgWWt+WXtsjr3aZ3tuYu77UOjvaQ5ADJdQc6s1aBfxcNXKTGLDMAYIcwLS4gEe+tXXzZQBA9ghyAAgcQQ4AgSPIASBwBDkABI4gB4DAZbJplpmNShpe4l9fKelYgs1JSl7bJeW3bbRrcfLaLim/bStau3rdfdXsg5kEeTvMbKDV7l9Zy2u7pPy2jXYtTl7bJeW3bWVpF6UVAAgcQQ4AgQsxyHdk3YAYeW2XlN+20a7FyWu7pPy2rRTtCq5GDgCYKcQeOQCgCUEOAIELKsjN7FIzO2BmB83s8xm2Y42Z/dLMnjCz/Wb2qcbxL5vZM2b2aOPr8gzadsjMHm88/0Dj2OvN7Gdm9lTjz1T3KTWzc5uuyaNm9pKZ3ZTV9TKz75jZ82a2r+lY7DUysy803nMHzOwvU27X183sd2b2WzO7x8xObRxfZ2avNF27W1NuV+xrl/H12tnUpkNm9mjjeJrXKy4fOvcec/cgviR1SXpa0npJ3ZIek3ReRm05Q9LGxvcrJD0p6TxJX5b0mYyv0yFJK2cd+2dJn298/3lJX8v4dfyDpN6srpekd0vaKGnffNeo8bo+Jmm5pLMb78GuFNv1AUnVxvdfa2rXuubHZXC9Wr52WV+vWb//hqRtGVyvuHzo2HsspB75hZIOuvuQu49LulPS1iwa4u7PuvvexvcvS3pC0llZtGWBtkr6buP770r6m+yaovdJetrdl7qyt23u/qCkF2YdjrtGWyXd6e6vuvvvJR1U/b2YSrvc/X53n2z8uEvS6k4892LbNYdMr1fEzEzShyV9rxPPPZc58qFj77GQgvwsSUeafh5RDsLTzNZJeoek3Y1Dn2x8DP5O2iWMBpd0v5kNmtkNjWNvdPdnpfqbTNIbMmhX5BrN/MeV9fWKxF2jPL3vPibpp00/n21mvzGz/zWzd2XQnlavXV6u17skPefuTzUdS/16zcqHjr3HQgpya3Es07mTZvY6SXdLusndX5L0b5LeJOntkp5V/aNd2i5y942SLpN0o5m9O4M2tGRm3ZKukPT9xqE8XK/55OJ9Z2Y3S5qUdHvj0LOS1rr7OyR9WtIdZvbnKTYp7rXLxfWSdK1mdhhSv14t8iH2oS2OLeqahRTkI5LWNP28WtLRjNoiM1um+ot0u7v/UJLc/Tl3n3L3mqT/UIc+Us7F3Y82/nxe0j2NNjxnZmc02n2GpOfTblfDZZL2uvtzjTZmfr2axF2jzN93ZvZRSVskfcQbRdXGx/Djje8HVa+rviWtNs3x2uXhelUlXSVpZ3Qs7evVKh/UwfdYSEG+R9I5ZnZ2o2d3jaQfZ9GQRv3t25KecPdvNh0/o+lhV0raN/vvdrhdf2ZmK6LvVR8o26f6dfpo42EflXRvmu1qMqOXlPX1miXuGv1Y0jVmttzMzpZ0jqRfp9UoM7tU0uckXeHuJ5qOrzKzrsb36xvtGkqxXXGvXabXq+ESSb9z95HoQJrXKy4f1Mn3WBqjuAmOBl+u+gjw05JuzrAdF6v+0ee3kh5tfF0u6b8kPd44/mNJZ6TcrvWqj34/Jml/dI0knSbp55Keavz5+gyu2cmSjks6pelYJtdL9f+ZPCtpQvXe0D/MdY0k3dx4zx2QdFnK7Tqoev00ep/d2njsBxuv8WOS9kr665TbFfvaZXm9Gsf/U9InZj02zesVlw8de4+xRB8AAhdSaQUA0AJBDgCBI8gBIHAEOQAEjiAHgMAR5AAQOIIcAAL3/3cduinUKN4UAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "d95s2 = d0+(d1-d0)*np.linspace(0,1,200)\n", + "rpbs2 = 0.05**(1/d95s2)\n", + "\n", + "lhc = np.log(np.linspace(1,70,200))\n", + "lhc = lhc/lhc.max()\n", + "rpbs3 = rpb0+(rpb1-rpb0)*lhc\n", + "\n", + "d95s3 = np.log(0.05)/np.log(rpbs3)\n", + "plt.plot(d95s3-d95s2,'.')\n" + ] + }, + { + "cell_type": "code", + "execution_count": 96, + "id": "fb9e2900-a442-46ee-8623-45c43f150b89", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 96, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAATFElEQVR4nO3df6zdd13H8ee7vSs6GK50BUd/rF0sw0Hkxy5bRWKGc9AOYtVo7DYiLpBmZlP8kcgIERX/wSDKDIWmqXWgy/oHW6AuhWFgyh/mQnsnYNtRuBa63XW6rhRUZmzv+vaPc848Oz33nnPv/Z57zvf7fT6Spvf7o+d8Pmn34sP7+/l8vpGZSJLKb9mwGyBJKoaBLkkVYaBLUkUY6JJUEQa6JFXE2LC++LLLLssNGzYM6+slqZQmJyefzszV3a4NLdA3bNjAoUOHhvX1klRKEXFitmuWXCSpIgx0SaoIA12SKsJAl6SKMNAlqSJ6BnpE7I2IpyLi8CzXIyL+KiKmIuIbEfH64pspSeqlnxH6PcCWOa5vBTY1f+0APrH4Zs1u8sQZdj48xeSJM4P8GkkqnZ7z0DPzyxGxYY5btgGfysY+vBMRcWlEXJ6ZTxbVyJbJE2e4dc8EZ2fOs2JsGfe+ezPXXLGy6K+RpFIqooa+Bni87Xi6ee4CEbEjIg5FxKFTp07N+4smjp/m7Mx5ziecmznPxPHTC2uxJFVQEYEeXc51fWtGZu7OzPHMHF+9uuvK1TltvnIVK8aWsTzgorFlbL5y1bw/Q5Kqqoil/9PAurbjtcDJAj73AtdcsZJ7372ZieOn2XzlKsstktSmiEDfD9wZEfuA64AfDKJ+3nLNFSsNcknqomegR8R9wPXAZRExDfwRcBFAZu4CDgA3AVPAM8Btg2qsJGl2/cxyubnH9QTuKKxFkqQFcaWoJFWEgS5JFWGgS1JFGOiSVBEGuiRVhIEuSRVhoEtSRZQ60N1KV5L+XxFL/4fCrXQl6flKO0J3K11Jer7SBrpb6UrS85W25OJWupL0fKUNdHArXUlqV9qSiyTp+Qx0SaoIA12SKsJAl6SKMNAlqSIMdEmqiEoFunu7SKqzUs9Db+feLpLqrjIjdPd2kVR3lQl093aRVHeVKbm4t4ukuqtMoIN7u0iqt8qUXCSp7gx0SaoIA12SKsJAl6SKqGygu2pUUt1UapZLi6tGJdVRXyP0iNgSEcciYioi7upy/cci4u8j4usRcSQibiu+qf1z1aikOuoZ6BGxHNgJbAWuBm6OiKs7brsDOJqZrwGuBz4SESsKbmvfXDUqqY76KblcC0xl5nGAiNgHbAOOtt2TwCUREcCLgO8BMwW3tW+uGpVUR/0E+hrg8bbjaeC6jns+BuwHTgKXAL+Wmec7PygidgA7ANavX7+Q9vbNVaOS6qafGnp0OZcdx28Fvga8HHgt8LGIePEFfyhzd2aOZ+b46tWr59lUSdJc+gn0aWBd2/FaGiPxdrcBD2TDFPAd4JXFNFGS1I9+Av0gsCkiNjYfdG6nUV5p9xhwA0BEvAy4CjheZEMXwznpkuqgZw09M2ci4k7gIWA5sDczj0TE7c3ru4A/Be6JiH+lUaJ5b2Y+PcB298056ZLqoq+FRZl5ADjQcW5X288ngbcU27RidJuTbqBLqqLKLv1vcU66pLqo5NL/ds5Jl1QXlQ90cE66pHqofMlFkurCQJekiqhdoDsnXVJV1aKG3uKcdElVVqsRuvukS6qyWgW6c9IlVVmtSi7OSZdUZbUKdHBOuqTqqlXJpZMzXiRVSe1G6C3OeJFUNbUdoTvjRVLV1DbQnfEiqWpqW3JxxoukqqltoIMzXiRVS21LLp2c8SKp7Go9Qm9xxoukKnCEjjNeJFWDgY4zXiRVgyUXnPEiqRoM9Kb2GS+TJ84Y7pJKx0Dv4ANSSWVlDb2DD0gllZWB3sEHpJLKypJLBx+QSiorA70LH5BKKiMDfQ4+IJVUJtbQ5+ADUkll0legR8SWiDgWEVMRcdcs91wfEV+LiCMR8U/FNnM4fEAqqUx6llwiYjmwE7gRmAYORsT+zDzads+lwMeBLZn5WES8dEDtXVI+IJVUJv3U0K8FpjLzOEBE7AO2AUfb7rkFeCAzHwPIzKeKbuiw+IBUUln0E+hrgMfbjqeB6zrueQVwUUT8I3AJcHdmfqrzgyJiB7ADYP369Qtp79D4gFTSqOunhh5dzmXH8RhwDfA24K3AH0bEKy74Q5m7M3M8M8dXr14978YOkw9IJY26fkbo08C6tuO1wMku9zydmT8EfhgRXwZeA3yrkFaOgNYD0nMz531AKmkk9RPoB4FNEbEReALYTqNm3u6zwMciYgxYQaMk85dFNnTYOh+QAux8eMp6uqSR0TPQM3MmIu4EHgKWA3sz80hE3N68viszH42IzwPfAM4DezLz8CAbPgytB6TW0yWNor5WimbmAeBAx7ldHccfBj5cXNNGV7d6uoEuadhcKboALjiSNIrcy2UBrKdLGkUG+gJZT5c0aiy5LJLz0yWNCgN9kaynSxoVllwWyXq6pFFhoBfAerqkUWDJpUDW0yUNk4FeoM56+sqLV7Dz4SkmT5wZdtMk1YAllwK119NXXryCDz54xPKLpCXjCL1g11yxkjve/BOceeas5RdJS8pAH5Bu0xknT5yxBCNpYCy5DEi36YzOgJE0SAb6ALW/j3Tnw1Pu0ChpoCy5LBFnwEgaNEfoS8QZMJIGzRH6EnIGjKRBMtCHwPKLpEGw5DIEll8kDYIj9CGx/CKpaAb6kFl+kVQUSy5DZvlFUlEcoY8Ayy+SimCgjxDLL5IWw5LLCLH8ImkxHKGPmNnKL/c/Mu1oXdKcHKGPqFb55dzMeZYvCz49Oc3Ms47WJc3OEfqIapVffu8tV/Gr4+uYedaHpZLmZqCPsFb55Zdfv9aHpZJ6suRSAj4sldQPR+gl4cNSSb30FegRsSUijkXEVETcNcd9b4iIZyPiV4protq1z1VvPSz9yBeOceueCUNdqrmeJZeIWA7sBG4EpoGDEbE/M492ue/PgIcG0VA1tJdfTn7/f7jvq489b7TeeoepZRipfvqpoV8LTGXmcYCI2AdsA4523PdbwP3AGwptoS7Qelfp5Ikz3P/ItFMbJQH9lVzWAI+3HU83zz0nItYAvwTsmuuDImJHRByKiEOnTp2ab1vVYa6pjdbWpfrpZ4QeXc5lx/FHgfdm5rMR3W5v/qHM3cBugPHx8c7P0AI4WpfU0k+gTwPr2o7XAic77hkH9jXD/DLgpoiYyczPFNFI9WZtXVI/gX4Q2BQRG4EngO3ALe03ZObG1s8RcQ/woGG+9BytS/XWM9AzcyYi7qQxe2U5sDczj0TE7c3rc9bNtfQcrUv1FJnDKWWPj4/noUOHhvLddTJ54gy37pl4brROhKN1qcQiYjIzx7tdc+l/xTlal+rDQK8Ba+tSPRjoNTLXaL21Ja8jdqm8DPSa6TZab23Je+ueCXdxlErMQK+p9tH65itXMXH89AW7ODpal8rFQK+x1mi9xVfeSeVmoAtwNoxUBQa6nuNsGKncDHRdwNG6VE4GurrqZ7T+gbe/ijPPnDXcpRHh0n/1NHnizAWj9WXAsmXB+UxLMdIScum/FqXbaD2iEeaWYqTRYaCrb+219ZUXr+CDDx6xFCONEEsuWjBLMdLSs+SigbAUI40WA12LZilGGg2WXFQ4SzHS4Fhy0ZKyFCMNh4GugbEUIy0tSy5aMpZipMWz5KKRYClGGiwDXUvOUow0GJZcNHT9lGIMd6nBkotGWq9SzNlz5/nAZw9bZ5d6MNA1MmYrxXTW2SeOnwaw1i51MNA1Utrfc3rVj19yQbhfNLaMlRev4NY9E5ydsdYutTPQNbK6hfvmK1cxcfw0Z2fOdy3HGO6qMwNdpdAe7gArxpb1rLUb7qobZ7molFozY2artbtgSVXlLBdVTq9auwuWVEd9BXpEbAHuBpYDezLzQx3XbwXe2zz8b+A3M/PrRTZUmk2vcHfBkuqiZ8klIpYD3wJuBKaBg8DNmXm07Z43Ao9m5pmI2Ar8cWZeN9fnWnLRoLlgSVW02JLLtcBUZh5vftg+YBvwXKBn5j+33T8BrF14c6VizHfBkuGususn0NcAj7cdTwNzjb7fBXyu24WI2AHsAFi/fn2fTZQWp58FS4a7qqCfQI8u57rWaSLizTQC/U3drmfmbmA3NEoufbZRWrT5PEQ13FVW/QT6NLCu7XgtcLLzpoj4KWAPsDUzTxfTPKl4hruqqp+HomM0HoreADxB46HoLZl5pO2e9cCXgF/vqKfPyoeiGjXzndtuuGsYFvVQNDNnIuJO4CEa0xb3ZuaRiLi9eX0X8AFgFfDxiACYme0LpVHlyF1l50pRqYeFrEoFd4PUYLhSVFqEhaxKfeCRaXeD1JIz0KV56Gd73wB3g9RQGOjSAs22vS/gQiYNhTV0aQCcMaNBmauGbqBLA2a4q0gGujQiDHctloEujSDDXQthoEsjznBXvwx0qUQWE+7ggqaqM9ClkppPuI8tC4jwzUwVZ6BLFdAr3Fv7XCeWaKrMQJcqplu4L2+O0J991vp7lRnoUoW1wr29hu5GYtXl5lxShbVvQdA6hsVvJAYGfdkY6FJFLWYjMR+wlpOBLtXAfDcSO/dsAknipmJlYg1dqrmiHrCCJZql4ENRSX1Z6ANWSzRLx0CXtCjOgR8dBrqkwliiGS4DXdJAWKJZega6pCVVZImm/X80DHsDXdIQLbZE88EHj3Rd+FTXcHelqKShmWsOfK+VrJ87/GTXhU/W47sz0CUtmflsU3DR2DK2vvpyDn73exeEvatauzPQJQ3dbKP4a65Y2XtPmj5XtUL1R/HW0CWVxkLr8VUaxftQVFLlzGfKZL+zalqfM8pBb6BLqo1BjOJhdILeQJdUS0WM4ketXGOgS1Kb+YziR20R1KIDPSK2AHcDy4E9mfmhjuvRvH4T8AzwG5n5yFyfaaBLGgW9RvFFLIJqfW4RQb+ohUURsRzYCdwITAMHI2J/Zh5tu20rsKn56zrgE83fJWmk9Zobv9hFUEtZsulnHvq1wFRmHgeIiH3ANqA90LcBn8rGcH8iIi6NiMsz88lCWytJS6SoRVBzzZO/992bCw31fgJ9DfB42/E0F46+u92zBnheoEfEDmAHwPr16+fbVkkauvkugpqtZHNu5jwTx08veaBHl3Odhfd+7iEzdwO7oVFD7+O7JWlkdRvF91uyuWhs2XPXitJPoE8D69qO1wInF3CPJNVGP7X5YdTQDwKbImIj8ASwHbil4579wJ3N+vp1wA+sn0vShTqDvkg9Az0zZyLiTuAhGtMW92bmkYi4vXl9F3CAxpTFKRrTFm8bSGslSbPqa7fFzDxAI7Tbz+1q+zmBO4ptmiRpPpYNuwGSpGIY6JJUEQa6JFWEgS5JFTG03RYj4hRwYoF//DLg6QKbUxZ17Hcd+wz17Hcd+wzz7/cVmbm624WhBfpiRMSh2XYbq7I69ruOfYZ69ruOfYZi+23JRZIqwkCXpIooa6DvHnYDhqSO/a5jn6Ge/a5jn6HAfpeyhi5JulBZR+iSpA4GuiRVROkCPSK2RMSxiJiKiLuG3Z5BiIh1EfFwRDwaEUci4j3N8y+JiH+IiG83fx/MHpxDFBHLI+JfIuLB5nEd+nxpRHw6Ir7Z/Dv/6Zr0+3eb/74PR8R9EfEjVet3ROyNiKci4nDbuVn7GBHva2bbsYh463y/r1SB3vbC6q3A1cDNEXH1cFs1EDPA72fmTwKbgTua/bwL+GJmbgK+2DyumvcAj7Yd16HPdwOfz8xXAq+h0f9K9zsi1gC/DYxn5qtpbM29ner1+x5gS8e5rn1s/je+HXhV8898vJl5fStVoNP2wurMPAu0XlhdKZn5ZGY+0vz5v2j8B76GRl8/2bztk8AvDqWBAxIRa4G3AXvaTle9zy8Gfhb4a4DMPJuZ36fi/W4aA340IsaAi2m85axS/c7MLwPf6zg9Wx+3Afsy838z8zs03i9x7Xy+r2yBPtvLqCsrIjYArwO+Arys9Sao5u8vHWLTBuGjwB8A59vOVb3PVwKngL9plpr2RMQLqXi/M/MJ4M+Bx2i8TP4HmfkFKt7vptn6uOh8K1ug9/Uy6qqIiBcB9wO/k5n/Oez2DFJEvB14KjMnh92WJTYGvB74RGa+Dvgh5S8z9NSsG28DNgIvB14YEe8YbquGbtH5VrZAr83LqCPiIhphfm9mPtA8/R8RcXnz+uXAU8Nq3wD8DPALEfFdGqW0n4uIv6PafYbGv+npzPxK8/jTNAK+6v3+eeA7mXkqM88BDwBvpPr9htn7uOh8K1ugP/fC6ohYQeMBwv4ht6lwERE0aqqPZuZftF3aD7yz+fM7gc8uddsGJTPfl5lrM3MDjb/XL2XmO6hwnwEy89+BxyPiquapG4CjVLzfNEotmyPi4ua/9xtoPCuqer9h9j7uB7ZHxAsiYiOwCfjqvD45M0v1i8bLqL8F/Bvw/mG3Z0B9fBON/6v1DeBrzV83AatoPBX/dvP3lwy7rQPq//XAg82fK99n4LXAoebf92eAlTXp958A3wQOA38LvKBq/Qbuo/GM4ByNEfi75uoj8P5mth0Dts73+1z6L0kVUbaSiyRpFga6JFWEgS5JFWGgS1JFGOiSVBEGuiRVhIEuSRXxf/TMYzXpCtBeAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "c=30\n", + "nx=100\n", + "wts = 1-np.log(np.linspace(1,c,nx))/np.log(c)\n", + "\n", + "plt.plot(wts,'.')\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b548d038-ab23-4243-972c-2f609771ffc9", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 55, + "id": "11ff9a42-6d2e-463e-ae16-d6aaba1f1833", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.05000000000000001" + ] + }, + "execution_count": 55, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.exp(np.log(0.05))" + ] + }, + { + "cell_type": "code", + "execution_count": 53, + "id": "d353fdbc-9960-46f6-a870-30c8e48722cf", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 53, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAWzUlEQVR4nO3df5BddX3G8fezuwmaihIhMphNsolGJaa2xp2QqT/GilaglJT4K6LVImmGGaFQ67QRpozjjCP2hyN/MM2kSCsVAVEyMg4ttFbLMGMgu7CBBBKJG5csoRDS1TANJWzy6R/3bDl7997dc5N777l7zvOa2cnec87d/ezZmyfffL7fc64iAjMzK66uvAswM7PWctCbmRWcg97MrOAc9GZmBeegNzMruJ68C6jljDPOiL6+vrzLMDObNQYHB5+PiAW19nVk0Pf19TEwMJB3GWZms4akkXr73LoxMys4B72ZWcE56M3MCs5Bb2ZWcA56M7OCyxT0ks6TtEfSXkmbauyfL2mrpEclPSRpZWrfn0naJWmnpNskvaqZP4CZmU1vxqCX1A3cCJwPrAA+KWlF1WHXAEMR8Q7gM8ANyXMXAn8K9EfESqAbWN+88s3MZr/BkTGu2foY1259jMGRsaZ//Szr6FcDeyNiGEDS7cBa4PHUMSuArwFExG5JfZLOTH2PV0t6GZgHHGhW8WZms8HgyBjbhg8xf95cdh74NQLe/sbXsfPAr3n+hZf4jz3PMX6scsv4OwdHue1P1vCuJfOb9v2zBP1CYH/q8ShwTtUxO4B1wAOSVgNLgN6IGJT0t8BTwIvAfRFxX61vImkjsBFg8eLFDf0QZmZ5GRwZ4wcPj04K7+og/+nPD/Ly+HGyvPvHy+PH2TZ8qO1Brxrbquu9HrhB0hDwGPAIMC5pPpXR/1LgV8Cdkj4dEd+Z8gUjtgBbAPr7+/1uKGbWERoZjTfDnJ4u1iw7vWlfD7IF/SiwKPW4l6r2S0QcBi4FkCRgX/LxYWBfRBxM9t0F/A4wJejNzPJSb1R+6ik93PTAPo4dj0yj8RPV0y0+8NY3sODUU1i3qrepo3nIFvTbgeWSlgJPU5lMvSR9gKTTgCMRcRTYANwfEYclPQWskTSPSuvmXMA3sTGztptuZP69gf1NHZVX6+kWn+hfVLO1I2hJuE/6/jMdEBHjkq4A7qWyaubmiNgl6fJk/2bgbOAWSceoTNJelux7UNL3gYeBcSotnS0t+UnMrPROZGQupvaiG5EejecV5DNRJ745eH9/f/julWZWT61Ab0W/HCr/EHR3iQ3vWcrhl8Y7LsQnSBqMiP5a+zryNsVmZtCeQO8S9HSJ99cZlY8dOcqaZad3RJifKAe9meWqXb3zeiPzIgT5TBz0ZtYWtQK92b3z6frlZQj0ehz0ZtYS6bbLiSxTnO64eoHeKf3yTuOgN7OT0qo++nS9cwd6Yxz0ZpZZdfulGYFe5t55uzjozWyKeqP0Ru7ZUs2Bnh8HvVnJtWKUDlP76A70/DjozUrmZCdJq3litPM56M0KrNmj9Vr3bHGgdz4HvVmBNGu07lF6sTjozWax6hH7iVxJWmuS1IFeLA56s1lkuhF71itJPUlaPg56sw7VaH+91laP1g0c9GYdZSLcX3jx5Yb767WuJPVo3cBBb5a7kwn39Ij91FfPcahbTQ56szY7mZUx7q/biXDQm7XBRLg3so7d/XVrFge9WQuc6LJHt2KsFRz0Zk1Sr9c+07JHh7u1moPe7CRkmUitfuw+u7Wbg96sQY2ukqle9ug+u7Wbg94sg0bD3e0Y6yQOerM6HO5WFA56sxSHuxWRg95Kz+FuReegt1JyuFuZOOitNNJXp2Z5k2uHuxWFg94KbyLgs1yd6nC3InLQWyE10ppxuFvROeitcL774FNc98Od04Z7+k2ufWWqFZ2D3goh3X//8e7nOHZ8asR3C849+0xfnWql46C3WSvrrX8nWjNfWbuSS85Z3N4izTpApqCXdB5wA9AN3BQR11ftnw/cDLwJ+F/gcxGxM9l3GnATsJLK/Z0+FxE/a9YPYOWTdXK1x313MyBD0EvqBm4EPgSMAtsl3R0Rj6cOuwYYioiLJb0tOf7cZN8NwL9GxEclzQXmNfUnsFJppP/u9oxZRZYR/Wpgb0QMA0i6HVgLpIN+BfA1gIjYLalP0pnAi8D7gD9O9h0FjjateiuFLP339K1/HfBmk2UJ+oXA/tTjUeCcqmN2AOuAByStBpYAvcAx4CDwj5J+CxgEroqI/6n+JpI2AhsBFi92H9VmbtF4ctUsmyxBrxrbqv/WXQ/cIGkIeAx4BBgH5gCrgCsj4kFJNwCbgL+a8gUjtgBbAPr7+7O8T7IVUNb17z2eXDXLLEvQjwKLUo97gQPpAyLiMHApgCQB+5KPecBoRDyYHPp9KkFvNsnE6P37g6PT3prA/XezxmUJ+u3AcklLgaeB9cAl6QOSlTVHkh78BuD+JPwPS9ov6a0RsYfKBO3jmKXMNMHqFo3ZyZkx6CNiXNIVwL1UllfeHBG7JF2e7N8MnA3cIukYlSC/LPUlrgRuTVbcDJOM/M0mRvF3bN9fc4LV69/NmkMRndcO7+/vj4GBgbzLsBaYqQffLVi/erFvTWDWIEmDEdFfa5+vjLW2malF4wlWs9Zw0FvLzdSimRjFu/9u1hoOemup6Ubx7sGbtYeD3lpiulG870Fj1l4Oemu6eqN4t2jM8uGgt6aY6X40nmg1y4+D3k5KlvvReBRvli8HvZ0wL5c0mx0c9NawmZZL+n40Zp3FQW8NmW6i1fejMetMDnrLZKblkm7RmHUuB71Na7rbB3ui1Wx2cNBbXdNNtnoUbzZ7OOitpsGRMa774U7GU20aAXO6xcc80Wo2qzjobYrBkTG++e8/n9SLd5vGbPZy0Nsk1e0a33jMbPZz0BtQe1WNgPcsP4OrP/gWj+LNZjEHvdVfG98lh7xZATjoSyzL2niHvNns56AvKd9K2Kw8HPQlVGvpJHhtvFlROehLxksnzcrHQV8iXjppVk4O+hLw0kmzcnPQF5yXTpqZg77AZpp0dciblYODvqA86WpmExz0BeRJVzNLc9AXTHW7xpOuZtaVdwHWPDXbNZ50NSs9j+gLYrp2jUPerNwc9LOc18ib2Uwc9LPY4MgYn7ppGy+9fNxr5M2srkw9eknnSdojaa+kTTX2z5e0VdKjkh6StLJqf7ekRyT9qFmFl91EP/7o+OSQ9xp5M6s244heUjdwI/AhYBTYLunuiHg8ddg1wFBEXCzpbcnx56b2XwU8Aby2aZWXWPVIvkuVgPebdptZLVlaN6uBvRExDCDpdmAtkA76FcDXACJit6Q+SWdGxLOSeoHfB74KfKGp1ZfUXQ+PvhLywLvf7H68mdWXpXWzENifejyabEvbAawDkLQaWAL0Jvu+CfwFcHy6byJpo6QBSQMHDx7MUFb5DI6Mcc3Wx7hjYP//t2t6eroc8mY2rSxBrxrbourx9cB8SUPAlcAjwLikC4HnImJwpm8SEVsioj8i+hcsWJChrHKZaNfc9uBTjB97ZXXNR9/lVo2ZTS9L62YUWJR63AscSB8QEYeBSwEkCdiXfKwHLpJ0AfAq4LWSvhMRn25C7aWSbtdAJeRPmdPFR1b1Tvc0M7NMQb8dWC5pKfA0lfC+JH2ApNOAIxFxFNgA3J+E/5eSDyS9H/iiQ74xE+vkv5du13SLT3ji1cwymjHoI2Jc0hXAvUA3cHNE7JJ0ebJ/M3A2cIukY1QmaS9rYc2lUWudvICP9y/iqxf/Zp6lmdkskumCqYi4B7inatvm1Oc/A5bP8DV+Cvy04QpLzO0aM2sGXxnboQZHxrjT7RozawIHfQeauOo1fatht2vM7EQ56DtMrate5/a4XWNmJ85B32F81auZNZvfeKSDTOnL+6pXM2sCB32HqNWX91WvZtYMbt10APflzayVHPQdwH15M2slB32Oat7ewH15M2syB31O6t3ewH15M2s2T8bmZNvwoUlvA+jbG5hZq3hEn5P58+bSJQHhtwE0s5Zy0OdgcGSMr/xoF8eOB91d4ssXreSScxbnXZaZFZRbN202sV5+om0TEYwdOZp3WWZWYB7Rt1Gt9fJzerpYs+z0vEszswJz0LeR18ubWR7cumkT38fGzPLioG8D38fGzPLk1k2L+T42ZpY3B32LuS9vZnlz66aF3Jc3s07goG+hbcOH3Jc3s9w56FtkcGSMp3/1Ij3dXXTL97Exs/y4R98CExOwR8eP09Ml1q9e7PvYmFluPKJvsvQtDo4HHDsevPG0VzvkzSw3HtE3kW9xYGadyEHfRF5KaWadyK2bJvFSSjPrVA76JvFSSjPrVA76JvBSSjPrZO7RnyQvpTSzTucR/UmaeJNvL6U0s07lEf1JSr/Jt5dSmlknyjSil3SepD2S9kraVGP/fElbJT0q6SFJK5PtiyT9RNITknZJuqrZP0Ce0m/y3SVx3YVv92jezDrOjEEvqRu4ETgfWAF8UtKKqsOuAYYi4h3AZ4Abku3jwJ9HxNnAGuDzNZ47K/lNvs1stsgyol8N7I2I4Yg4CtwOrK06ZgXwY4CI2A30STozIp6JiIeT7S8ATwALm1Z9TiYmYB948nmOh6+ANbPOliXoFwL7U49HmRrWO4B1AJJWA0uASesLJfUB7wQerPVNJG2UNCBp4ODBg5mKz8vEBGz6CthbN6xx28bMOlKWoFeNbVH1+HpgvqQh4ErgESptm8oXkF4D/AC4OiIO1/omEbElIvojon/BggVZas/NxARsl2DuHF8Ba2adLcuqm1FgUepxL3AgfUAS3pcCSBKwL/lA0hwqIX9rRNzVhJpzlZ6A7e7yBKyZdb4sI/rtwHJJSyXNBdYDd6cPkHRasg9gA3B/RBxOQv9bwBMR8Y1mFp6XdNvGE7BmNhvMGPQRMQ5cAdxLZTL1exGxS9Llki5PDjsb2CVpN5XVORPLKN8N/BHwAUlDyccFTf8p2qT6VgeegDWz2UAR1e32/PX398fAwEDeZUxSfauDj/Uv8q0OzKxjSBqMiP5a+3wLhIx8qwMzm618C4SMfKsDM5utPKLPwLc6MLPZzEGfgVfamNls5qCfgVfamNls5x79NPymImZWBB7RT8MrbcysCBz001iz7HTm9rhlY2azm1s3M1i3qhclf3o0b2azkYO+jnR/fm5PF+tW9c78JDOzDuTWTR3p/vzL48fZNnwo75LMzE6Ig74GL6k0syJx66aKl1SaWdF4RF/FSyrNrGgc9FW8pNLMisatmxq8pNLMisRBn+IllWZWRG7dpHhJpZkVkYM+4SWVZlZUbt3gJZVmVmwe0eMllWZWbA56vKTSzIrNrZuEl1SaWVGVPui9pNLMiq70rRsvqTSzoit90Ls/b2ZFV+rWzeDIGNuGD3HdhW9n7MhR1iw73f15Myuc0gZ9dW/+1g1rHPJmVkilbd24N29mZVHaoHdv3szKopStG/fmzaxMShf07s2bWdmUrnXj3ryZlU2moJd0nqQ9kvZK2lRj/3xJWyU9KukhSSuzPrfd3Js3s7KZsXUjqRu4EfgQMApsl3R3RDyeOuwaYCgiLpb0tuT4czM+t+18XxszK5MsPfrVwN6IGAaQdDuwFkiH9QrgawARsVtSn6QzgWUZnts2vq+NmZVRltbNQmB/6vFosi1tB7AOQNJqYAnQm/G5JM/bKGlA0sDBgwezVd8g9+fNrIyyBL1qbIuqx9cD8yUNAVcCjwDjGZ9b2RixJSL6I6J/wYIFGcpqnPvzZlZGWVo3o8Ci1ONe4ED6gIg4DFwKIEnAvuRj3kzPbad3LZnPrRvWsG34kNfOm1lpZAn67cBySUuBp4H1wCXpAySdBhyJiKPABuD+iDgsacbntsvERVJrlp3O53/3zXmUYGaWixmDPiLGJV0B3At0AzdHxC5Jlyf7NwNnA7dIOkZlovWy6Z7bmh+lPl8kZWZllunK2Ii4B7inatvm1Oc/A5ZnfW671ZqEddCbWVmU4spYT8KaWZmV5l43vkjKzMqq8EHvi6TMrOwK37rxRVJmVnaFD3r3582s7ArfuvFFUmZWdoUOel8kZWZW4KD3RVJmZhWF7dF7EtbMrKKwQe9JWDOzisK2bjwJa2ZWUdgRfXoi1iFvZmVWyBG9J2LNzF5RyBG9J2LNzF5RyKD3RKyZ2SsK2brxRKyZ2SsKF/S+GtbMbLJCBb0nYc3MpipUj96TsGZmUxUq6D0Ja2Y2VaFaN56ENTObqlBBD5Wwd8Cbmb2iUK2bwZExbvzJXgZHxvIuxcysYxRmRO8VN2ZmtRVmRO8VN2ZmtRUm6L3ixsystsK0brzixsystsIEPXjFjZlZLYVp3ZiZWW0OejOzgnPQm5kVnIPezKzgHPRmZgXnoDczKzhFRN41TCHpIDBygk8/A3i+ieU0i+tqXKfW5roa47oadyK1LYmIBbV2dGTQnwxJAxHRn3cd1VxX4zq1NtfVGNfVuGbX5taNmVnBOejNzAquiEG/Je8C6nBdjevU2lxXY1xX45paW+F69GZmNlkRR/RmZpbioDczK7jCBL2k8yTtkbRX0qYc61gk6SeSnpC0S9JVyfYvS3pa0lDycUFO9f1S0mNJDQPJttdL+jdJTyZ/tvVez5LemjovQ5IOS7o6j3Mm6WZJz0namdpW9/xI+lLymtsj6cM51PY3knZLelTSVkmnJdv7JL2YOneb21xX3d9du85ZnbruSNX0S0lDyfZ2nq96GdG611lEzPoPoBv4BbAMmAvsAFbkVMtZwKrk81OBnwMrgC8DX+yAc/VL4IyqbX8NbEo+3wR8Peff5X8BS/I4Z8D7gFXAzpnOT/J73QGcAixNXoPdba7t94Ce5POvp2rrSx+Xwzmr+btr5zmrVVfV/r8DrsvhfNXLiJa9zooyol8N7I2I4Yg4CtwOrM2jkIh4JiIeTj5/AXgCWJhHLQ1YC3w7+fzbwB/mVwrnAr+IiBO9MvqkRMT9wH9Xba53ftYCt0fESxGxD9hL5bXYttoi4r6IGE8ebgN6W/X9G6lrGm07Z9PVJUnAx4HbWvG9pzNNRrTsdVaUoF8I7E89HqUDwlVSH/BO4MFk0xXJf7Fvbnd7JCWA+yQNStqYbDszIp6ByosQeENOtQGsZ/Jfvk44Z/XOT6e97j4H/Evq8VJJj0j6T0nvzaGeWr+7Tjln7wWejYgnU9vafr6qMqJlr7OiBL1qbMt13aik1wA/AK6OiMPA3wNvAn4beIbKfxvz8O6IWAWcD3xe0vtyqmMKSXOBi4A7k02dcs7q6ZjXnaRrgXHg1mTTM8DiiHgn8AXgu5Je28aS6v3uOuWcfZLJA4q2n68aGVH30BrbGjpnRQn6UWBR6nEvcCCnWpA0h8ov8NaIuAsgIp6NiGMRcRz4B1r4X/zpRMSB5M/ngK1JHc9KOiup/SzguTxqo/KPz8MR8WxSY0ecM+qfn4543Un6LHAh8KlImrrJf/MPJZ8PUunrvqVdNU3zu8v9nEnqAdYBd0xsa/f5qpURtPB1VpSg3w4sl7Q0GRWuB+7Oo5Ck9/ct4ImI+EZq+1mpwy4GdlY/tw21/YakUyc+pzKRt5PKufpscthngR+2u7bEpFFWJ5yzRL3zczewXtIpkpYCy4GH2lmYpPOAvwQuiogjqe0LJHUnny9LahtuY131fne5nzPgg8DuiBid2NDO81UvI2jl66wds8xtmsm+gMrs9S+Aa3Os4z1U/lv1KDCUfFwA/DPwWLL9buCsHGpbRmX2fgewa+I8AacDPwaeTP58fQ61zQMOAa9LbWv7OaPyD80zwMtURlKXTXd+gGuT19we4PwcattLpX878VrbnBz7keR3vAN4GPiDNtdV93fXrnNWq65k+z8Bl1cd287zVS8jWvY68y0QzMwKriitGzMzq8NBb2ZWcA56M7OCc9CbmRWcg97MrOAc9GZmBeegNzMruP8DcTME0Ykok94AAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(rpbs,'.')" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "b253ab3b-947e-49e4-aea3-ca6065866f77", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "8ec727a8-436f-4bbc-be4c-377eb87630d3", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Help on built-in function arange in module numpy:\n", + "\n", + "arange(...)\n", + " arange([start,] stop[, step,], dtype=None, *, like=None)\n", + " \n", + " Return evenly spaced values within a given interval.\n", + " \n", + " Values are generated within the half-open interval ``[start, stop)``\n", + " (in other words, the interval including `start` but excluding `stop`).\n", + " For integer arguments the function is equivalent to the Python built-in\n", + " `range` function, but returns an ndarray rather than a list.\n", + " \n", + " When using a non-integer step, such as 0.1, the results will often not\n", + " be consistent. It is better to use `numpy.linspace` for these cases.\n", + " \n", + " Parameters\n", + " ----------\n", + " start : integer or real, optional\n", + " Start of interval. The interval includes this value. The default\n", + " start value is 0.\n", + " stop : integer or real\n", + " End of interval. The interval does not include this value, except\n", + " in some cases where `step` is not an integer and floating point\n", + " round-off affects the length of `out`.\n", + " step : integer or real, optional\n", + " Spacing between values. For any output `out`, this is the distance\n", + " between two adjacent values, ``out[i+1] - out[i]``. The default\n", + " step size is 1. If `step` is specified as a position argument,\n", + " `start` must also be given.\n", + " dtype : dtype\n", + " The type of the output array. If `dtype` is not given, infer the data\n", + " type from the other input arguments.\n", + " like : array_like\n", + " Reference object to allow the creation of arrays which are not\n", + " NumPy arrays. If an array-like passed in as ``like`` supports\n", + " the ``__array_function__`` protocol, the result will be defined\n", + " by it. In this case, it ensures the creation of an array object\n", + " compatible with that passed in via this argument.\n", + " \n", + " .. note::\n", + " The ``like`` keyword is an experimental feature pending on\n", + " acceptance of :ref:`NEP 35 `.\n", + " \n", + " .. versionadded:: 1.20.0\n", + " \n", + " Returns\n", + " -------\n", + " arange : ndarray\n", + " Array of evenly spaced values.\n", + " \n", + " For floating point arguments, the length of the result is\n", + " ``ceil((stop - start)/step)``. Because of floating point overflow,\n", + " this rule may result in the last element of `out` being greater\n", + " than `stop`.\n", + " \n", + " See Also\n", + " --------\n", + " numpy.linspace : Evenly spaced numbers with careful handling of endpoints.\n", + " numpy.ogrid: Arrays of evenly spaced numbers in N-dimensions.\n", + " numpy.mgrid: Grid-shaped arrays of evenly spaced numbers in N-dimensions.\n", + " \n", + " Examples\n", + " --------\n", + " >>> np.arange(3)\n", + " array([0, 1, 2])\n", + " >>> np.arange(3.0)\n", + " array([ 0., 1., 2.])\n", + " >>> np.arange(3,7)\n", + " array([3, 4, 5, 6])\n", + " >>> np.arange(3,7,2)\n", + " array([3, 5])\n", + "\n" + ] + } + ], + "source": [ + "help(np.arange)" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "54a2639d-edcf-4dd8-b516-1a8ffa72fb01", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAPJ0lEQVR4nO3df6zddX3H8eeLXhuHgFS4Y9I2tCwE6Iz446bp5rKQYRSYsYvZH5A5l05CmoDD/cjs9E//8cfcdKFZ0yhTIht/KCZo3HBzJmbJQG6hpWsL81oEKjgv6obZP7X2vT/O1+Xsci73e3vv7bn3w/ORnPSe7/fzvffzyUme/fZ7bs83VYUkqV3njHsCkqSVZeglqXGGXpIaZ+glqXGGXpIaNzHuCYxy8cUX15YtW8Y9DUlaMw4cOPB8VU2O2rcqQ79lyxamp6fHPQ1JWjOSPDXfPi/dSFLjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNa5X6JNcn+SJJDNJ9ozYvyHJl5I8luRbSV7Xbd+c5BtJjiU5kuSO5V6AJOmlLRj6JOuAvcANwDbg5iTb5gz7IHCwql4PvAf4VLf9FPAnVXU1sAO4bcSxkqQV1OeMfjswU1XHq+okcC+wc86YbcDXAarqcWBLkkuq6rmqeqTb/hPgGLBx2WYvSVpQn9BvBJ4Zen6CF8f6EPAugCTbgcuATcMDkmwB3gg8NOqHJLk1yXSS6dnZ2V6TlyQtrE/oM2JbzXn+EWBDkoPA+4BHGVy2GXyD5Dzgi8D7q+qFUT+kqvZX1VRVTU1OTvaZuySph4keY04Am4eebwKeHR7QxXsXQJIAT3YPkryCQeTvqar7lmHOkqRF6HNG/zBwRZKtSdYDNwH3Dw9IcmG3D+AW4JtV9UIX/c8Ax6rqL5dz4pKkfhY8o6+qU0luBx4A1gF3VdWRJLu7/fuAq4G7k/wMOAq8tzv8LcDvAYe7yzoAH6yqry7vMiRJ8+lz6YYuzF+ds23f0Nf/Blwx4rh/ZfQ1fknSWeL/jJWkxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWpcr9AnuT7JE0lmkuwZsX9Dki8leSzJt5K8ru+xkqSVtWDok6wD9gI3ANuAm5NsmzPsg8DBqno98B7gU4s4VpK0gvqc0W8HZqrqeFWdBO4Fds4Zsw34OkBVPQ5sSXJJz2MlSSuoT+g3As8MPT/RbRt2CHgXQJLtwGXApp7H0h13a5LpJNOzs7P9Zi9JWlCf0GfEtprz/CPAhiQHgfcBjwKneh472Fi1v6qmqmpqcnKyx7QkSX1M9BhzAtg89HwT8OzwgKp6AdgFkCTAk93j3IWOlSStrD5n9A8DVyTZmmQ9cBNw//CAJBd2+wBuAb7ZxX/BYyVJK2vBM/qqOpXkduABYB1wV1UdSbK7278PuBq4O8nPgKPAe1/q2JVZiiRplFSNvGQ+VlNTUzU9PT3uaUjSmpHkQFVNjdrn/4yVpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeqmnA0/9mL3fmOHAUz8e91SkRZkY9wSkteDAUz/mdz/9ICdPnWb9xDncc8sO3nzZhnFPS+rFM3qphweP/5CTp05zuuCnp07z4PEfjntKUm+GXuphx+UXsX7iHNYFXjFxDjsuv2jcU5J689KN1MObL9vAPbfs4MHjP2TH5Rd52UZriqGXenrzZRsMvNYkL91IUuMMvSQ1ztBLUuMMvSQ1ztBLUuMMvSQ1rlfok1yf5IkkM0n2jNj/6iRfTnIoyZEku4b2/VG37d+T/H2SVy7nAiRJL23B0CdZB+wFbgC2ATcn2TZn2G3A0aq6BrgW+ESS9Uk2An8ITFXV64B1wE3LOH9J0gL6nNFvB2aq6nhVnQTuBXbOGVPA+UkCnAf8CDjV7ZsAfiHJBHAu8OyyzFyS1Euf0G8Enhl6fqLbNuxO4GoGET8M3FFVp6vqe8BfAE8DzwH/XVVfG/VDktyaZDrJ9Ozs7CKXIUmaT5/QZ8S2mvP87cBB4FLgDcCdSS5IsoHB2f/Wbt+rkrx71A+pqv1VNVVVU5OTkz2nL0laSJ/QnwA2Dz3fxIsvv+wC7quBGeBJ4CrgrcCTVTVbVT8F7gN+benTliT11Sf0DwNXJNmaZD2DN1PvnzPmaeA6gCSXAFcCx7vtO5Kc212/vw44tlyTlyQtbMFPr6yqU0luBx5g8Fszd1XVkSS7u/37gA8Dn01ymMGlng9U1fPA80m+ADzC4M3ZR4H9K7MUSdIoqZp7uX38pqamanp6etzTkKQ1I8mBqpoatc//GStJjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktS4XqFPcn2SJ5LMJNkzYv+rk3w5yaEkR5LsGtp3YZIvJHk8ybEkv7qcC5AkvbQFQ59kHbAXuAHYBtycZNucYbcBR6vqGuBa4BNJ1nf7PgX8Y1VdBVwDHFumuUuSeuhzRr8dmKmq41V1ErgX2DlnTAHnJwlwHvAj4FSSC4DfAD4DUFUnq+q/lmvykqSF9Qn9RuCZoecnum3D7gSuBp4FDgN3VNVp4HJgFvjbJI8m+XSSV436IUluTTKdZHp2dnax65AkzaNP6DNiW815/nbgIHAp8Abgzu5sfgJ4E/A3VfVG4H+AF13jB6iq/VU1VVVTk5OT/WYvSVpQn9CfADYPPd/E4Mx92C7gvhqYAZ4EruqOPVFVD3XjvsAg/JKks6RP6B8GrkiytXuD9Sbg/jljngauA0hyCXAlcLyqvg88k+TKbtx1wNFlmbkkqZeJhQZU1akktwMPAOuAu6rqSJLd3f59wIeBzyY5zOBSzweq6vnuW7wPuKf7S+I4g7N/SdJZkqq5l9vHb2pqqqanp8c9DUlaM5IcqKqpUfv8n7GS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mN6xX6JNcneSLJTJI9I/a/OsmXkxxKciTJrjn71yV5NMlXlmvikqR+Fgx9knXAXuAGYBtwc5Jtc4bdBhytqmuAa4FPJFk/tP8O4NiyzFiStCh9zui3AzNVdbyqTgL3AjvnjCng/CQBzgN+BJwCSLIJ+C3g08s2a0lSb31CvxF4Zuj5iW7bsDuBq4FngcPAHVV1utv3SeDPgNO8hCS3JplOMj07O9tjWpKkPvqEPiO21ZznbwcOApcCbwDuTHJBkncAP6iqAwv9kKraX1VTVTU1OTnZY1qSpD76hP4EsHno+SYGZ+7DdgH31cAM8CRwFfAW4J1Jvsvgks9vJvn8kmctSeqtT+gfBq5IsrV7g/Um4P45Y54GrgNIcglwJXC8qv68qjZV1ZbuuH+pqncv2+wlSQuaWGhAVZ1KcjvwALAOuKuqjiTZ3e3fB3wY+GySwwwu9Xygqp5fwXlLknpK1dzL7eM3NTVV09PT456GJK0ZSQ5U1dTIfasx9ElmgafGPY9Fuhh4uf0rxjW/PLjmteGyqhr5myyrMvRrUZLp+f42bZVrfnlwzWufn3UjSY0z9JLUOEO/fPaPewJj4JpfHlzzGuc1eklqnGf0ktQ4Qy9JjTP0i5DkNUn+Kcm3uz83zDNuoRu1/GmSSnLxys96aZa65iQfT/J4kseSfCnJhWdt8ovQ4zVLkr/u9j+W5E19j12tznTNSTYn+UaSY92Nhu44+7M/M0t5nbv9a/MmSlXlo+cD+Biwp/t6D/DREWPWAd8BLgfWA4eAbUP7NzP4OImngIvHvaaVXjPwNmCi+/qjo44f92Oh16wbcyPwDww+4mMH8FDfY1fjY4lrfi3wpu7r84H/aH3NQ/v/GPg74CvjXs9iHp7RL85O4HPd158DfnvEmIVu1PJXDD6ff628C76kNVfV16rqVDfuQQaffrra9Lm5zk7g7hp4ELgwyWt7HrsanfGaq+q5qnoEoKp+wuDucXPvUbEaLeV1XtM3UTL0i3NJVT0H0P35iyPGzHujliTvBL5XVYdWeqLLaElrnuMPGJwtrTZ95j/fmL5rX22Wsub/k2QL8EbgoeWf4rJb6po/SY+bKK1GC3565ctNkn8GfmnErg/1/RYjtlWSc7vv8bYzndtKWak1z/kZH2Jwe8l7Fje7s6LPzXXmG9Pn2NVoKWse7EzOA74IvL+qXljGua2UM17z8E2Ukly73BNbaYZ+jqp663z7kvznz//p2v1z7gcjhs13o5ZfBrYChwa31mUT8EiS7VX1/WVbwBlYwTX//Hv8PvAO4LrqLnSuMn1urjPfmPU9jl2NlrJmkryCQeTvqar7VnCey2kpa/4dBjdRuhF4JXBBks/XWrm/xrjfJFhLD+Dj/P83Jj82YswEcJxB1H/+hs+vjBj3XdbGm7FLWjNwPXAUmBz3Wl5ijQu+ZgyuzQ6/Sfetxbzeq+2xxDUHuBv45LjXcbbWPGfMtayxN2PHPoG19AAuAr4OfLv78zXd9kuBrw6Nu5HBbyJ8B/jQPN9rrYR+SWsGZhhc8zzYPfaNe03zrPNF8wd2A7u7rwPs7fYfBqYW83qvxseZrhn4dQaXPB4bel1vHPd6Vvp1Hvoeay70fgSCJDXO37qRpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMb9L649IpnJtRRAAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(rpbs,'.')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "499c64ac-b495-4d7f-8ef0-87f394a6111d", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0d193f6c-87e7-4bc3-94e3-bbff6ecc646c", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "37b717ba-cb80-4ace-befc-2cdbd4b1dd4f", + "metadata": { + "tags": [] + }, + "source": [ + "### examine a pft file" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "3ac988fa-de04-44d7-af83-54b312b171d7", + "metadata": {}, + "outputs": [], + "source": [ + "ds = xr.open_dataset(basefile)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "bc331c0b-0c9b-44f2-94f0-82e92aab4161", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'kmax' (segment: 4, pft: 79)>\n",
+       "array([[0.000000e+00, 1.602140e-08, 2.234574e-08, ..., 2.000000e-08,\n",
+       "        2.000000e-08, 2.000000e-08],\n",
+       "       [0.000000e+00, 1.602140e-08, 2.234574e-08, ..., 2.000000e-08,\n",
+       "        2.000000e-08, 2.000000e-08],\n",
+       "       [0.000000e+00, 1.602140e-08, 2.234574e-08, ..., 2.000000e-08,\n",
+       "        2.000000e-08, 2.000000e-08],\n",
+       "       [0.000000e+00, 1.602140e-08, 2.234574e-08, ..., 2.000000e-08,\n",
+       "        2.000000e-08, 2.000000e-08]])\n",
+       "Coordinates:\n",
+       "    pftname  (pft) |S40 b'not_vegetated                           ' ... b'irr...\n",
+       "  * segment  (segment) |S40 b'sunlit                                  ' ... b...\n",
+       "Dimensions without coordinates: pft\n",
+       "Attributes:\n",
+       "    units:      mm h2o (transpired)/mm h2o (water potential gradient)/sec\n",
+       "    long_name:  plant segment max conductance
" + ], + "text/plain": [ + "\n", + "array([[0.000000e+00, 1.602140e-08, 2.234574e-08, ..., 2.000000e-08,\n", + " 2.000000e-08, 2.000000e-08],\n", + " [0.000000e+00, 1.602140e-08, 2.234574e-08, ..., 2.000000e-08,\n", + " 2.000000e-08, 2.000000e-08],\n", + " [0.000000e+00, 1.602140e-08, 2.234574e-08, ..., 2.000000e-08,\n", + " 2.000000e-08, 2.000000e-08],\n", + " [0.000000e+00, 1.602140e-08, 2.234574e-08, ..., 2.000000e-08,\n", + " 2.000000e-08, 2.000000e-08]])\n", + "Coordinates:\n", + " pftname (pft) |S40 ...\n", + " * segment (segment) |S40 b'sunlit ' ... b...\n", + "Dimensions without coordinates: pft\n", + "Attributes:\n", + " units: mm h2o (transpired)/mm h2o (water potential gradient)/sec\n", + " long_name: plant segment max conductance" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ds.kmax" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "5c41f87c-d99b-4bf6-9ae6-04351571c2a7", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'medlynslope' (pft: 79)>\n",
+       "array([0.  , 2.35, 2.35, 2.35, 4.12, 4.12, 4.45, 4.45, 4.45, 4.7 , 4.7 , 4.7 ,\n",
+       "       2.22, 5.25, 1.62, 5.79, 5.79, 1.79, 1.79, 5.79, 5.79, 5.79, 5.79, 5.79,\n",
+       "       5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79,\n",
+       "       5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79,\n",
+       "       5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79,\n",
+       "       5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 1.79, 1.79, 5.79, 5.79, 1.79,\n",
+       "       1.79, 1.79, 1.79, 1.79, 1.79, 5.79, 5.79])\n",
+       "Coordinates:\n",
+       "    pftname  (pft) |S40 b'not_vegetated                           ' ... b'irr...\n",
+       "Dimensions without coordinates: pft\n",
+       "Attributes:\n",
+       "    long_name:  Medlyn slope of conductance-photosynthesis relationship\n",
+       "    units:      umol H2O/umol CO2\n",
+       "    comment:    Values come from the values used in the CABLE model
" + ], + "text/plain": [ + "\n", + "array([0. , 2.35, 2.35, 2.35, 4.12, 4.12, 4.45, 4.45, 4.45, 4.7 , 4.7 , 4.7 ,\n", + " 2.22, 5.25, 1.62, 5.79, 5.79, 1.79, 1.79, 5.79, 5.79, 5.79, 5.79, 5.79,\n", + " 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79,\n", + " 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79,\n", + " 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79,\n", + " 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 1.79, 1.79, 5.79, 5.79, 1.79,\n", + " 1.79, 1.79, 1.79, 1.79, 1.79, 5.79, 5.79])\n", + "Coordinates:\n", + " pftname (pft) |S40 b'not_vegetated ' ... b'irr...\n", + "Dimensions without coordinates: pft\n", + "Attributes:\n", + " long_name: Medlyn slope of conductance-photosynthesis relationship\n", + " units: umol H2O/umol CO2\n", + " comment: Values come from the values used in the CABLE model" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ds.medlynslope" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "27f831de-8631-4f30-800f-fd1c8735ee92", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'jmaxb0' ()>\n",
+       "array(0.0311)\n",
+       "Attributes:\n",
+       "    long_name:  Baseline proportion of nitrogen allocated for electron transport\n",
+       "    units:      J
" + ], + "text/plain": [ + "\n", + "array(0.0311)\n", + "Attributes:\n", + " long_name: Baseline proportion of nitrogen allocated for electron transport\n", + " units: J" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ds.jmaxb0" + ] + }, + { + "cell_type": "markdown", + "id": "580e8d55-fbde-418c-a364-9800c0ebd6e9", + "metadata": {}, + "source": [ + "## edit a paramfile\n", + " - avoiding xarray for writing the new files, for some reason paramfiles written with xarray tend to crash the model" + ] + }, + { + "cell_type": "markdown", + "id": "f3d8f8b4-b282-4509-b99d-b5c266d9d61a", + "metadata": {}, + "source": [ + "### jmaxb0, increase to 0.05" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "id": "027b148d-3f2f-48c1-b4cc-e79044f27a0b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "## CREATE the new paramfile\n", + "newfile = 'jmaxb0_5e-2.nc'\n", + "cmd = 'cp '+basefile+' '+newfile\n", + "os.system(cmd)" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "id": "7772f5b6-99ba-4f9c-aae9-aaadfcf2d77a", + "metadata": {}, + "outputs": [], + "source": [ + "param = 'jmaxb0'\n", + "dset = netCDF4.Dataset(newfile,'r+')\n", + "dset[param][:] = 0.05\n", + "dset.close()" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "id": "6d39c684-ae26-482b-be5a-1fd9d2160e84", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'jmaxb0' ()>\n",
+       "array(0.05)\n",
+       "Attributes:\n",
+       "    long_name:  Baseline proportion of nitrogen allocated for electron transport\n",
+       "    units:      J
" + ], + "text/plain": [ + "\n", + "array(0.05)\n", + "Attributes:\n", + " long_name: Baseline proportion of nitrogen allocated for electron transport\n", + " units: J" + ] + }, + "execution_count": 26, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "xr.open_dataset(newfile)['jmaxb0']" + ] + }, + { + "cell_type": "markdown", + "id": "5dba3b01-589e-4829-8a8f-448933904552", + "metadata": {}, + "source": [ + "### kmax, double all values\n", + " - note that the shapes of the various parameters vary\n", + " - some are scalar, some have a pft-dimension, and some have pft and plant-segment dimensions\n", + " - easist to read in the example from basefile, and edit accordingly " + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "id": "93ef7e76-0041-425e-8ad6-3f4396e314fe", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "## CREATE the new paramfile\n", + "newfile = 'kmax_2x.nc'\n", + "cmd = 'cp '+basefile+' '+newfile\n", + "os.system(cmd)" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "id": "bcf2c053-b2ce-4406-a1ba-41d9684a23f2", + "metadata": {}, + "outputs": [], + "source": [ + "param = 'kmax'\n", + "kmax = xr.open_dataset(basefile)[param].values\n", + "newkmax = kmax*2" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "id": "1baf2f03-4da4-40b9-887e-f8afeb686b3d", + "metadata": {}, + "outputs": [], + "source": [ + "dset = netCDF4.Dataset(newfile,'r+')\n", + "dset[param][:] = newkmax\n", + "dset.close()" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "id": "62f01983-27c6-4c32-949b-a95580638685", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'kmax' (segment: 4, pft: 79)>\n",
+       "array([[0.000000e+00, 3.204280e-08, 4.469149e-08, ..., 4.000000e-08,\n",
+       "        4.000000e-08, 4.000000e-08],\n",
+       "       [0.000000e+00, 3.204280e-08, 4.469149e-08, ..., 4.000000e-08,\n",
+       "        4.000000e-08, 4.000000e-08],\n",
+       "       [0.000000e+00, 3.204280e-08, 4.469149e-08, ..., 4.000000e-08,\n",
+       "        4.000000e-08, 4.000000e-08],\n",
+       "       [0.000000e+00, 3.204280e-08, 4.469149e-08, ..., 4.000000e-08,\n",
+       "        4.000000e-08, 4.000000e-08]])\n",
+       "Coordinates:\n",
+       "    pftname  (pft) |S40 b'not_vegetated                           ' ... b'irr...\n",
+       "  * segment  (segment) |S40 b'sunlit                                  ' ... b...\n",
+       "Dimensions without coordinates: pft\n",
+       "Attributes:\n",
+       "    units:      mm h2o (transpired)/mm h2o (water potential gradient)/sec\n",
+       "    long_name:  plant segment max conductance
" + ], + "text/plain": [ + "\n", + "array([[0.000000e+00, 3.204280e-08, 4.469149e-08, ..., 4.000000e-08,\n", + " 4.000000e-08, 4.000000e-08],\n", + " [0.000000e+00, 3.204280e-08, 4.469149e-08, ..., 4.000000e-08,\n", + " 4.000000e-08, 4.000000e-08],\n", + " [0.000000e+00, 3.204280e-08, 4.469149e-08, ..., 4.000000e-08,\n", + " 4.000000e-08, 4.000000e-08],\n", + " [0.000000e+00, 3.204280e-08, 4.469149e-08, ..., 4.000000e-08,\n", + " 4.000000e-08, 4.000000e-08]])\n", + "Coordinates:\n", + " pftname (pft) |S40 ...\n", + " * segment (segment) |S40 b'sunlit ' ... b...\n", + "Dimensions without coordinates: pft\n", + "Attributes:\n", + " units: mm h2o (transpired)/mm h2o (water potential gradient)/sec\n", + " long_name: plant segment max conductance" + ] + }, + "execution_count": 37, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "xr.open_dataset(newfile)['kmax']" + ] + }, + { + "cell_type": "markdown", + "id": "84ba4bae-1b4e-4b5f-af7f-2411971fccd2", + "metadata": {}, + "source": [ + "### medlynslope\n", + "- substituting a value for just one pft" + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "id": "6b184d47-246d-461f-9f11-09d599a2b689", + "metadata": {}, + "outputs": [], + "source": [ + "p = xr.open_dataset(basefile)\n", + "mypft = 'broadleaf_deciduous_temperate_tree'\n", + "ix = np.array([mypft in str(n) for n in p.pftname.values]) #index vector for mypft" + ] + }, + { + "cell_type": "code", + "execution_count": 53, + "id": "006f879a-085e-4bea-904e-c0f1b4eaeb15", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0. , 2.3499999 , 2.3499999 , 2.3499999 , 4.11999989,\n", + " 4.11999989, 4.44999981, 6. , 4.44999981, 4.69999981,\n", + " 4.69999981, 4.69999981, 2.22000003, 5.25 , 1.62 ,\n", + " 5.78999996, 5.78999996, 1.78999996, 1.78999996, 5.78999996,\n", + " 5.78999996, 5.78999996, 5.78999996, 5.78999996, 5.78999996,\n", + " 5.78999996, 5.78999996, 5.78999996, 5.78999996, 5.78999996,\n", + " 5.78999996, 5.78999996, 5.78999996, 5.78999996, 5.78999996,\n", + " 5.78999996, 5.78999996, 5.78999996, 5.78999996, 5.78999996,\n", + " 5.78999996, 5.78999996, 5.78999996, 5.78999996, 5.78999996,\n", + " 5.78999996, 5.78999996, 5.78999996, 5.78999996, 5.78999996,\n", + " 5.78999996, 5.78999996, 5.78999996, 5.78999996, 5.78999996,\n", + " 5.78999996, 5.78999996, 5.78999996, 5.78999996, 5.78999996,\n", + " 5.78999996, 5.78999996, 5.78999996, 5.78999996, 5.78999996,\n", + " 5.78999996, 5.78999996, 1.78999996, 1.78999996, 5.78999996,\n", + " 5.78999996, 1.78999996, 1.78999996, 1.78999996, 1.78999996,\n", + " 1.78999996, 1.78999996, 5.78999996, 5.78999996])" + ] + }, + "execution_count": 53, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "param = 'medlynslope'\n", + "m = p[param].values\n", + "m[ix] = 6\n", + "m" + ] + }, + { + "cell_type": "code", + "execution_count": 54, + "id": "7b0a1d7d-1b76-48a1-bb8d-5b7ae496bc22", + "metadata": {}, + "outputs": [], + "source": [ + "newfile = 'medlynslope_BDT_6.nc'\n", + "cmd = 'cp '+basefile+' '+newfile\n", + "os.system(cmd)\n", + "dset = netCDF4.Dataset(newfile,'r+')\n", + "dset[param][:] = m\n", + "dset.close()" + ] + }, + { + "cell_type": "code", + "execution_count": 55, + "id": "300f8235-0c9f-4cc3-92a9-26ae32f4a31c", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'medlynslope' (pft: 79)>\n",
+       "array([0.  , 2.35, 2.35, 2.35, 4.12, 4.12, 4.45, 6.  , 4.45, 4.7 , 4.7 , 4.7 ,\n",
+       "       2.22, 5.25, 1.62, 5.79, 5.79, 1.79, 1.79, 5.79, 5.79, 5.79, 5.79, 5.79,\n",
+       "       5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79,\n",
+       "       5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79,\n",
+       "       5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79,\n",
+       "       5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 1.79, 1.79, 5.79, 5.79, 1.79,\n",
+       "       1.79, 1.79, 1.79, 1.79, 1.79, 5.79, 5.79])\n",
+       "Coordinates:\n",
+       "    pftname  (pft) |S40 b'not_vegetated                           ' ... b'irr...\n",
+       "Dimensions without coordinates: pft\n",
+       "Attributes:\n",
+       "    long_name:  Medlyn slope of conductance-photosynthesis relationship\n",
+       "    units:      umol H2O/umol CO2\n",
+       "    comment:    Values come from the values used in the CABLE model
" + ], + "text/plain": [ + "\n", + "array([0. , 2.35, 2.35, 2.35, 4.12, 4.12, 4.45, 6. , 4.45, 4.7 , 4.7 , 4.7 ,\n", + " 2.22, 5.25, 1.62, 5.79, 5.79, 1.79, 1.79, 5.79, 5.79, 5.79, 5.79, 5.79,\n", + " 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79,\n", + " 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79,\n", + " 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79,\n", + " 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 5.79, 1.79, 1.79, 5.79, 5.79, 1.79,\n", + " 1.79, 1.79, 1.79, 1.79, 1.79, 5.79, 5.79])\n", + "Coordinates:\n", + " pftname (pft) |S40 ...\n", + "Dimensions without coordinates: pft\n", + "Attributes:\n", + " long_name: Medlyn slope of conductance-photosynthesis relationship\n", + " units: umol H2O/umol CO2\n", + " comment: Values come from the values used in the CABLE model" + ] + }, + "execution_count": 55, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "xr.open_dataset(newfile)[param]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db275759-d218-4a83-af53-2fedb86563ad", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python [conda env:miniconda3-ppe-py]", + "language": "python", + "name": "conda-env-miniconda3-ppe-py-py" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}