The Earth as Transformed by Human Action: Global and Regional Changes in the Biosphere Over the Past 300 Years,B. L. Turner, CUP Archive, 1990, 27

1. 1. Mesopotamia population Python 3 code
2. ver 00

import matplotlib import matplotlib.pyplot as plt import numpy as np import pandas as pd from scipy.interpolate import interp1d from scipy.optimize import curve_fit import math as math

def doubling_time(years1, populations1, idx1, idx2):

   t1=years1[idx1]
   t2=years1[idx2]
   p1=populations1[idx1]
   p2=populations1[idx2]
   dt1=t2-t1
   dp1=p2/p1
   r1=math.log(dp1)/dt1
   doublingtime1=math.log(2)/r1
   return(r1, doublingtime1)

1. 1. sigmoid

def sigmoid(x, L ,x0, k, b):

   y = L / (1 + np.exp(-k*(x-x0)))+b
   return (y)

1. 1. exp

def func(x, a, c, d):

   return a*np.exp(-c*x)+d

1. 1. main program

1. I made a map-thing
2. https://www.nma-fallout.com/threads/i-made-a-map-thing.218744/

1. year sigmapop biggestpop numcenters
2. 5000 8500 6500 3
3. 4700 19500 8500 4
4. 3600 81000 14000 12
5. 3300 126000 40000 12
6. 3100 131000 45000 12
7. 2900 264000 80000 12

1. years=np.array([-5000,-4700,-3600,-3300,-3100,-2900])
2. populations=np.array([8500,19500,81000,126000,131000,264000])
3. biggestplaces=np.array([6500,8500,14000,40000,45000,80000])
4. numplaces=np.array([3,4,12,12,12,12])

1. df = pd.read_csv(r'mesopotamia1.csv', delimiter=';')
2. bc;est_pop;pop_lower;pop_upper
3. years1=np.array(df['Year'])
4. populations1=np.array(df['Population'])

years1=np.array([-4100,-3500, -2500, -1900, -1600]) populations1=np.array([25000,110000, 360000,630000, 270000])

population0=populations1[0] year0=years1[0]

beginx=-4100 endx=-1900 xnum=25 xnum1=xnum+1 pxlen=25

timestep1=abs(endx-beginx)/xnum

times25=np.linspace(beginx, endx, xnum)

times1=np.linspace(beginx,endx,int(abs(endx-beginx))) timerels1=times1-times1[0]

1. 1. data turner 1990
2. The Earth as Transformed by Human Action: Global and Regional Changes in the Biosphere Over the Past 300 Years,B. L. Turner, CUP Archive, 1990, 27

origx=years1 origy=populations1

sekondx=np.linspace(beginx, endx,xnum1) xdata0=np.linspace(0, xnum,xnum1)

interp1 = interp1d(origx,origy, kind='linear')

ydata0 = interp1(sekondx)

ro1, doublingtime1=doubling_time(years1, populations1, 0, 3)

print(doublingtime1, ro1)

1. quit(-1)

1. print (sekondx)
2. rint (xdata0)
3. print (ydata0)

xdata1=xdata0 ydata1=ydata0/1000000.0

1. print(ydata1)

1. quit(-1)

1. p0 = [max(ydata1), np.median(xdata1),1,min(ydata1)] # this is an mandatory initial guess

1. popt, pcov = curve_fit(sigmoid, xdata1, ydata1,p0, method='dogbox')

1. popt, pcov = curve_fit(func, xdata1, ydata1, p0=(1, 1e-6, 1))
2. popt, pcov = curve_fit(func, xdata1, ydata1)

1. popt, pcov = curve_fit(func, xdata1, ydata1,p0=(1, 1e-6, 1), method='lm')

popt, pcov = curve_fit(func, xdata1, ydata1,p0=(1, 1e-6, 1), method='trf')

print ("Popt") print (popt) print (pcov)

1. quit(-1)

x=np.linspace(0, xnum,pxlen)

1. y = sigmoid(x, *popt)

y = func(x, *popt)*1e6

r1=0.15e-2 r2=0.175e-2 r3=0.2e-2 r4=0.25e-2

1. print("r ", r1)

1. quit(-1)

1. timerel1=x-x[0]

1. y2=origy[0]*np.exp(timerel1*r1)

yo1=origy[0]*np.exp(timerels1*ro1) y2=origy[0]*np.exp(timerels1*r1) y3=origy[0]*np.exp(timerels1*r2) y4=origy[0]*np.exp(timerels1*r3) y5=origy[0]*np.exp(timerels1*r4)

print(timerels1)

1. print(origy[0])
2. print(y2)

1. quit(-1)

annot2="Growth rate is "+str(round((ro1*100),3)) + "  %\n"+ "Doubling time is "+ str( round((doublingtime1)) ) + " years"

annot1="Initial population "+str(round((population0),2)) + " \n at year "+str(round((abs(year0)),0))+ " BC"

1. ytop1=0.63e6

ytop1=1.5e6

y3=np.where(y3>ytop1, ytop1, y3) y4=np.where(y4>ytop1, ytop1, y4) y5=np.where(y5>ytop1, ytop1, y5)

1. print(x)
2. print(y)
3. print(y2)

plt.title("Population growth of South Mesopotamia", size=22) plt.xlabel("Years (negative is BC)" , fontsize=15) plt.ylabel("Population of South Mesopotamia" , fontsize=15 ) plt.xticks(fontsize=15) plt.yticks(fontsize=15)

plt.text(-4100, 250000, "4100 BC - 25000", alpha=0.4, fontsize=17) plt.text(-2200, 630000*1.05, "1900 BC - 630000", alpha=0.4, fontsize=17) plt.text(-4000, 600000, annot2, color="green", alpha=0.4, fontsize=17)

plt.plot(origx, origy, lw=5, label="Estimation by Turner 1990")

1. plt.plot(x,y, lw=3, linestyle=":" ,label="Fitted")

plt.plot(times1,yo1, lw=3,linestyle="--" ,alpha=0.5, label="estim") plt.plot(times1,y2, lw=3,linestyle="--" ,alpha=0.5, label="r=0.15%") plt.plot(times1,y3, lw=3,linestyle="--" ,alpha=0.5, label="r=0.175%") plt.plot(times1,y4, lw=3,linestyle="--" ,alpha=0.5, label="r=0.20%") plt.plot(times1,y5, lw=3,linestyle="--" ,alpha=0.5, label="r=0.30%") plt.legend(fontsize=16) plt.grid()

plt.show()